【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}$

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

　【疑問１】何で積分する必要があるの？（自分で考えてください．近いうちに説明します．）

　次に，積分時系列の増分

　 $\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}$

の関係が成り立ちます．

　【疑問２】 $\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)