本記事は、 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|$ を振幅と呼びます1。
位相
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 変換したものです。
作り方は次の通りです。
-
$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 変換する。
```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 サロゲートのパワースペクトルは、元のデータのパワースペクトルと一致します。
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 サロゲートは、データの分布を保存しながら、パワースペクトルもあまり変化しないように工夫したサロゲートです。
作り方は次の通りです。
-
標準正規分布 $\mathrm{N}(0, 1)$ に従う乱数を $N$ 個発生させ、$x$ と同じ大小関係になるように並べ替えたものを $g$ とする。すなわち、$g_0 < g_1 < \dots < g_{N-1}$.
-
$g$ のFT サロゲート $\hat{g}$ を作る。
-
元の時系列 $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 サロゲートです。
サロゲートデータのパワースペクトルが元のデータのパワースペクトルに近づくように並べ替える操作を繰り返して作ります。
-
$x$ の RS サロゲート又は AAFT サロゲートを作り、これを初期値 $s^{(0)}$ とする。
-
以下を $s^{(i)}$ が変化しなくなる（つまり $s^{(i)} = s^{(i-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)}$ とする。
```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$
- 反復後の $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.
-
正確には、振幅の定数倍です。 ↩