ケィオスの時系列解析メモランダム

時系列解析,生体情報学,数学・物理などの解説です.

【Rでフラクタル解析】ランダムウォーク解析,Fluctuation analysis,2次構造関数解析

今回は,ランダムウォーク解析の導入です.detrended fluctuation analysis (DFA)とか,detrending moving average analysis (DMA)とかを使ってる研究者の方が,そこそこいると思います.ランダムウォーク解析は,DFA,DMAの原型となる,素朴な方法です.ランダムウォーク解析は,Fluctuation analysis (FA),2次構造関数解析と同じものです.

chaos-kiyono.hatenablog.com
chaos-kiyono.hatenablog.com

サンプル時系列(左)とその積分時系列(右)

ランダムウォーク解析

 解析対象は,弱定常過程の観測時系列\left\{x[i] \right\}です.事前に,平均0にしておいてください. 

 ランダムウォーク解析では,まず,観測時系列\left\{x[i] \right\}積分時系列を計算します.

 \begin{equation}
\displaystyle y[i] = \sum_{j=1}^{i} x[j]
\end{equation}

このように積分し,観測時系列をランダムウォーク軌道に変換します.

 【疑問1】何で積分する必要があるの?(自分で考えてください.近いうちに説明します.)


 次に,積分時系列の増分

 \begin{eqnarray}
 \Delta_{s} y[i] &=& y[i+s] - y[i] = \sum_{j=1}^{s} x[i+j]
 \end{eqnarray}

の2乗平均の期待値 (ゆらぎ関数: fluctuation function)

 \begin{equation}
 F^2(s) = \left\langle \Delta_{s} y^2[i] \right\rangle
 \end{equation}

を推定します.ここで, \left\langle \ \cdot \ \right\rangleは期待値を表します.今回は,確率変数も小文字で表しました.推定量は,標本の平均として計算してください.

 \left\{x[i] \right\}パワースペクトルが,低周波数側で,1/f^{\beta}に比例するとき,
 
 \begin{equation}
     F^2(s) \propto s^{2\alpha}
 \end{equation}

となります.平方根をとると,

 \begin{equation}
     F(s) \propto s^{\alpha}
 \end{equation}

です.

スケーリング指数\alpha\betaには

\begin{eqnarray}
\alpha = \frac{\beta + 1}{2} \quad (0 < \alpha < 1)
\end{eqnarray}

の関係が成り立ちます.

 【疑問2】\alpha\betaの関係をどうやって導くの?(自分で考えてください.近いうちに説明します.)

今回は,疑問の答えは与えません.これ以外にも,疑問に思う部分について考えてみてください.

数値実験

 下のRスクリプトを実行してみてください.白色ノイズのサンプル時系列xを,ランダムウォーク解析します.この例では,\betaはいくつで,推定結果の\alphaはいくつになるでしょうか?考えてみてください.

# 時系列の長さ
N <- 1000
#####################
# サンプルの時系列 (白色ノイズ)
# S(f) ~ f^0 (beta=0)の例
x <- rnorm(N)
x <- x - mean(x) # 平均0に
#####################

#####################
# ランダムウォーク解析
###
# 積分
y <- cumsum(x)
###
# 解析するスケールs
s <- unique(round(exp(seq(log(1),log(N/4),length.out=20))))
n.s <- length(s)
###
F.s <- c() # ゆらぎ関数
###
for(i in 1:n.s){
 D.y <- y[1:(N-s[i])]-y[(1+s[i]):N]
 F.s[i] <- sqrt(mean(D.y^2))
}
# log-logをとって直線をあてはめ
fit <- lm(log10(F.s)~log10(s))
###
# log-logプロット
plot(log10(s),log10(F.s),col=4,xlab="log10 s",ylab="log10 F(s)",
 main=paste("slope =",format(fit$coefficients[[2]],digits=3)))
# 直線フィット
abline(fit,col=2,lty=2)
実行例