Help us understand the problem. What is going on with this article?

サロゲートデータ生成

More than 1 year has passed since last update.

本記事は、 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 サロゲートは時系列の分布を保存します。

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

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

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

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

      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 サロゲートである。

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.

  1. 正確には、振幅の定数倍です。 

Why do not you register as a user and use Qiita more conveniently?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
Comments
Sign up for free and join this conversation.
If you already have a Qiita account
Why do not you register as a user and use Qiita more conveniently?
You need to log in to use this function. Qiita can be used more conveniently after logging in.
You seem to be reading articles frequently this month. Qiita can be used more conveniently after logging in.
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away