サロゲートデータ生成

本記事は、 [U-TOKYO AP Advent Calendar 2018](https://qiita.com/advent-calendar/2018/ut_ap) の3日目です。

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

**目次**
- [はじめに](#はじめに)
- [Random shuffle surrogates](#random-shuffle-rs-surrogates)
- [Fourier transform surrogates](#fourier-transform-ft-surrogates)
- [Amplitude adjusted Fourier transform surrogates](#amplitude-adjusted-fourier-transform-aaft-surrogates)
- [Iterated amplitude adjusted Fourier transform surrogates](#iterated-amplitude-adjusted-fourier-transform-iaaft-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 変換と呼びます。
変換後の係数は大文字のアルファベットで表します：

```math
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 変換**

```math
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|$ を振幅と呼びます[^amp]。

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

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

[^amp]: 正確には、振幅の定数倍です。

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

```math
s_n = x_{p_n},\quad n = 0, 1, \dots, {N-1}
```

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

```matlab:rssurrogate.m
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 変換する。

    ```math
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 変換する。

    ```math
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 サロゲートのパワースペクトルは、元のデータのパワースペクトルと一致します。

```matlab:ftsurrogate.m
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 サロゲートは、元の時系列を並べ替えたものなので、分布は保存されています。
さらに、パワースペクトルも、元のデータのパワースペクトルに近いものが得られます。

```matlab:aaftsurrogate.m
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 変換する。
    
        ```math
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)}$ とする。

        ```math
\hat S^{(i)}_k = \frac{\left|X_k\right|}{\left|S^{(i)}_k\right|} S^{(i)}_k
```
    3. $\hat S^{(i)}$ を逆 Fourier 変換する。

        ```math
\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$ を並べ替える。

        ```math
s^{(i+1)}_n = x'_{\mathrm{rank}\,\hat s^{(i)}_n}
```
    5. $i \leftarrow i+1$

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

```matlab:iaaftsurrogate.m
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.