本記事は、 U-TOKYO AP Advent Calendar 2018 の3日目です。

サロゲートデータの作り方を解説します。
解説の後に Matlab のソースコードを付けています。

目次
- はじめに
- Random shuffle surrogates
- Fourier transform surrogates
- Amplitude adjusted Fourier transform surrogates
- Iterated amplitude adjusted Fourier transform surrogates

はじめに

サロゲートデータとは

時系列データ $x=(x_0,x_1,\dots ,x_{N-1})$ が与えられたとします。
サロゲートデータとは、 $x$ の「特徴」の一部だけを保存した時系列 $s=(s_0,s_1,\dots ,s_{N-1})$ のことです。

記号・用語の定義

時系列データ
時系列データを小文字のアルファベットで表します。
実数値を取る時系列 $x=(x_0,x_1,\dots ,x_{N-1})\in {\mathbf{R}}^N$ を考えます。
以下、元のデータを $x$, サロゲートデータを $s$ と書きます。

虚数単位
$\mathrm{i}=\sqrt{-1}$ とアップライト体で書きます。（イタリック体の $i$ は反復回数のインデックスとして使います。）

離散 Fourier 変換
本記事では離散 Fourier 変換しか登場しないので、離散 Fourier 変換のことを単に Fourier 変換と呼びます。
変換後の係数は大文字のアルファベットで表します：

逆離散 Fourier 変換

振幅
Fourier 変換の絶対値 $\left|X_k\right|$ を振幅と呼びます¹。

位相
Fourier 変換の偏角 $\arg X_k$ を位相 (phase) と呼びます。

パワースペクトル
Fourier 変換の絶対値の2乗 ${\left|X_k\right|}^2$ をパワースペクトルと呼びます。
パワースペクトルはその周波数成分のエネルギーを表しています。

Random shuffle (RS) surrogates

Random shuffle surrogates は、 $x_0,x_1,\dots ,x_{N-1}$ をランダムに入れ替えたサロゲートデータです。
つまり、 $p=(p_0,p_1,\dots ,p_{N-1})$ を $\{0,1,\dots ,N-1\}$ の順列（置換）として、

で表されます。
RS サロゲートは時系列の分布を保存します。

function y = rssurrogate(x)
    perm = randperm(length(x));
    y = x(perm);
end

Fourier transform (FT) surrogates

RS サロゲートはパワースペクトルを保存しません。
つまり、 $x$ の Fourier 変換の振幅 $|X|$ と $s$ の Fourier 変換の振幅 $|S|$ は全く別のものとなります。
FT サロゲート（phase-randomised surrogates とも言います）は、元の時系列を Fourier 変換し、振幅は保存しつつ位相をランダムに変えて、逆 Fourier 変換したものです。
作り方は次の通りです。

1. $x$ を Fourier 変換する。

   $$X_k=\frac{1}{\sqrt{N}}\sum _{n=0}^{N-1}{x}_n\exp (-\frac{\mathrm{i}2\pi kn}{N})$$

2. 各 $X_k$ の位相をランダムに変換する。
   次のステップ3で逆 Fourier 変換した際にきちんと実数値の時系列になるように、前半はランダムに位相を変換し、後半は逆位相となるように変換する。
   実装上は、時系列の長さ $N$ が偶数か奇数かによって分けて考える。
   $u_k$ を独立な $[0,1]$ 上の一様乱数とする。

   • $N$ が偶数のとき
     • $k=0,N/2$ はそのまま ${\hat{X}}_0=X_0$, ${\hat{X}}_{N/2}=X_{N/2}$,
     • $k=1,\dots ,N/2-1$ に対して ${\hat{X}}_k=X_k\exp \left(\mathrm{i}2\pi u_k\right)$,
     • $k=N/2+1,\dots ,N-1$ に対して ${\hat{X}}_k=X_k\exp (-\mathrm{i}2\pi u_{N-k})$.
   • $N$ が奇数のとき
     • $k=0$ はそのまま ${\hat{X}}_0=X_0$,
     • $k=1,\dots ,(N-1)/2$ に対して ${\hat{X}}_k=X_k\exp \left(\mathrm{i}2\pi u_k\right)$,
     • $k=(N+1)/2,\dots ,N-1$ に対して ${\hat{X}}_k=X_k\exp (-\mathrm{i}2\pi u_{N-k})$.

3. ${\hat{X}}_k$ を逆 Fourier 変換する。

   $$s_n=\frac{1}{\sqrt{N}}\sum _{n=0}^{N-1}{\hat{X}}_k\exp \left(\frac{\mathrm{i}2\pi kn}{N}\right)$$

ステップ2で、位相は変わりますが、振幅は $\left|X_k\right|=\left|{\hat{X}}_k\right|$ と保存されます。
従って、 FT サロゲートのパワースペクトルは、元のデータのパワースペクトルと一致します。

function z = ftsurrogate(x)
    y = fft(x);
    n = length(x);
    if mod(n, 2) == 0
        l = n/2 - 1;
        r = exp(2i*pi*rand(l,1));
        v = [1; r; 1; flip(conj(r))];
    else
        l = (n-1)/2;
        r = exp(2i*pi*rand(l,1));
        v = [1; r; flip(conj(r))];
    end
    z = ifft(y.*v);
end

Amplitude adjusted Fourier transform (AAFT) surrogates

FT サロゲートは、パワースペクトルを保存する一方で、データの分布を変えてしまいます。
AAFT サロゲートは、データの分布を保存しながら、パワースペクトルもあまり変化しないように工夫したサロゲートです。
作り方は次の通りです。

1. 標準正規分布 $\mathrm{N}(0,1)$ に従う乱数を $N$ 個発生させ、$x$ と同じ大小関係になるように並べ替えたものを $g$ とする。すなわち、$g_0<g_1<\cdots <g_{N-1}$.

2. $g$ のFT サロゲート $\hat{g}$ を作る。

3. 元の時系列 $x$ を、大小関係が $\hat{g}$ と同じになるように並べ替えたものが AAFT サロゲートである。
   すなわち、 $x$ を小さい順に並べたものを $x^{\prime }$ とし、 ${\hat{g}}_n$ が $\hat{g}$ の $i$ 番目に小さい要素であるとき $\mathrm{rank}{\hat{g}}_n=i-1$ と書くことにすると、 $s_n=x_{\mathrm{rank}{\hat{g}}_n}^{\prime }$.

AAFT サロゲートは、元の時系列を並べ替えたものなので、分布は保存されています。
さらに、パワースペクトルも、元のデータのパワースペクトルに近いものが得られます。

function z = aaftsurrogate(x)
    n = length(x);
    r = sort(randn(n,1));
    [y, I] = sort(x);
    [I2, J] = sort(I);
    g = r(J);
    h = ftsurrogate(g);
    [h2, K] = sort(h);
    [K2, L] = sort(K);
    z = y(L);
end

Iterated amplitude adjusted Fourier transform (IAAFT) surrogates

通常の AAFT サロゲートよりも、さらにパワースペクトルを元の時系列に近づけたのが IAAFT サロゲートです。
サロゲートデータのパワースペクトルが元のデータのパワースペクトルに近づくように並べ替える操作を繰り返して作ります。

1. $x$ の RS サロゲート又は AAFT サロゲートを作り、これを初期値 $s^{(0)}$ とする。

2. 以下を $s^{(i)}$ が変化しなくなる（つまり $s^{(i)}=s^{(i-1)}$ になる）まで繰り返す。

   1. $s^{(i)}$ を Fourier 変換する。

      $$S_k^{(i)}=\frac{1}{\sqrt{N}}\sum _{n=0}^{N-1}{s}_n^{(i)}\exp (-\frac{\mathrm{i}2\pi kn}{N})$$

   2. 振幅を元のデータの振幅で置き換えたものを ${\hat{S}}^{(i)}$ とする。

      $${\hat{S}}_k^{(i)}=\frac{\left|X_k\right|}{\left|S_k^{(i)}\right|}{S}_k^{(i)}$$

   3. ${\hat{S}}^{(i)}$ を逆 Fourier 変換する。

      $${\hat{s}}_n^{(i)}=\frac{1}{\sqrt{N}}\sum _{n=0}^{N-1}{\hat{S}}_k^{(i)}\exp \left(\frac{\mathrm{i}2\pi kn}{N}\right)$$

   4. ${\hat{s}}^{(i)}$ の大小関係を使って $x$ を並べ替える。

      $$s_n^{(i+1)}=x_{\mathrm{rank}{\hat{s}}_n^{(i)}}^{\prime }$$

   5. $i\leftarrow i+1$

3. 反復後の $s^{(i)}$ が IAAFT サロゲートである。

function z = iaaftsurrogate(x)
    [sx, I] = sort(x);
    c = abs(fft(x));
%    y = rssurrogate(x); % 初期値RSサロゲート
    y = aaftsurrogate(x); % 初期値AAFTサロゲート
    count = 0;
    while(count < 100) % 繰り返しは100回まで
        count = count + 1
        f = fft(y);
        g = c.*f./abs(f);
        h = ifft(g);
        [h2, K] = sort(h);
        [K2, L] = sort(K);
        z = sx(L);
        if(z == y)
            break
        end
        y = z;
    end
end

参考文献

• Theiler, J., Eubank, S., Longtin, A., Galdrikian, B., & Farmer, J. D. (1992). Testing for nonlinearity in time series: the method of surrogate data. Physica D: Nonlinear Phenomena, 58(1-4), 77-94.
• Schreiber, T., & Schmitz, A. (1996). Improved Surrogate Data for Nonlinearity Tests. Physical Review Letters, 77(4), 635.
• Schreiber, T., & Schmitz, A. (2000). Surrogate time series. Physica D: Nonlinear Phenomena, 142(3-4), 346-382.