本記事は、 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 変換と呼びます。
変換後の係数は大文字のアルファベットで表します：

X_k = \mathrm{DFT} \left[x_n\right] = \frac{1}{\sqrt{N}} \sum_{n=0}^{N-1} x_n \exp\left(-\frac{\mathrm{i} 2 \pi k n}{N} \right),\quad k = 1, \dots, N-1.

逆離散 Fourier 変換

x_n = \mathrm{DFT}^{-1} \left[X_k\right] = \frac{1}{\sqrt{N}} \sum_{n=0}^{N-1} X_k \exp\left(\frac{\mathrm{i} 2 \pi k n}{N} \right),\quad n = 1, \dots, N-1.

振幅
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}\}$ の順列（置換）として、

s_n = x_{p_n},\quad n = 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\left(-\frac{\mathrm{i} 2 \pi k n}{N}\right)
   

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\left(-\mathrm{i} 2 \pi u_{N-k} \right)$.
   • $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\left(-\mathrm{i} 2 \pi u_{N-k} \right)$.

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 k n}{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 < \dots < g_{N-1}$.

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

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

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^{(i)}_k = \frac{1}{\sqrt{N}}\sum_{n=0}^{N-1} s^{(i)}_n \exp\left(-\frac{\mathrm{i} 2 \pi k n}{N}\right)
      

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

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

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

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

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

      s^{(i+1)}_n = x'_{\mathrm{rank}\,\hat s^{(i)}_n}
      

   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.