2022-06-30

【時系列解析】ランダムウォーク解析の解析的考察（その4：スケーリング指数の関係）

「ランダムウォーク解析」の話の続きです．前回の記事は，これです．
chaos-kiyono.hatenablog.com

　ランダムウォーク解析では，(1) 観測時系列の平均を０にして，(2) 積分して，(3) その増分の2乗平均の平方根 $F(s)$ をいろんなスケール $s$ で計算して，(4) $\log s$ 対 $\log F(s)$ のプロットの傾きをスケーリング指数 $\alpha$ として求める，のでした．ここでは，

　 $F(s) \propto s^{\alpha}$

の関係がみられる場合を考えています．

　これまで説明してなかったかもしれませんが，「スケーリング」とはべき関数（ $a$ を定数として $x^a$ の形）で表される関係のことで，「スケーリング指数」とはべき指数（ $x^a$ のときは $a$ ）のことです．

　 $F(s) \propto s^{\alpha}$ の関係がみられる代表的なケースに，以下の「長期記憶」，「長時間相関」，「フラクタル性 (自己アフィン性)」があります．

時系列の長期記憶，長時間相関，フラクタル性

長期記憶過程では，自己共分散関数 $C(k)$ が，

　 $C(k) \propto |k|^{-\gamma} \quad (0 < \gamma < 1)$

となり，パワースペクトル $S(f)$ が，

　 $S(f) \propto |f|^{-\beta} \quad (0 < \beta < 1)$

となります．スケーリング指数 $\alpha$ と $\beta$ の間に

　 $\gamma+\beta=1$

の関係が成り立ちます．

長時間相関過程は長期記憶過程を含みますが，長期記憶過程に対応する正の長時間相関だけでなく，負の長時間相関もあります．長時間相関過程では，パワースペクトル $S(f)$ が，

　 $S(f) \propto |f|^{-\beta} \quad (-1 < \beta < 1)$

となります．負の長時間相関では，自己共分散関数 $C(k)$ が負の値もとるので，絶対値をとると，

　 $|C(k)| \propto |k|^{-\gamma} \quad (1 < \gamma < 2)$

となってます． $\beta=0$ と， $\beta=1$ の場合を除いて， $\gamma+\beta=1$ は成り立ってます．

フラクタル過程 (自己アフィン性をもつ時系列)では，パワースペクトル $S(f)$ が，

　 $S(f) \propto |f|^{-\beta} \quad (1 < \beta < 3)$

となります．非定常なので，自己共分散は， $C(k)$ のように，ラグ $k$ の関数として書くことはできません．

　ここまでに，３つのスケーリング指数， $\alpha$ ， $\beta$ ， $\gamma$ が登場しました．今回は，ランダムウォーク解析で， $\alpha$ と $\beta$ ，あるいは， $\alpha$ と $\gamma$ の関係を解析的に導きます．

基本事項の確認

1. $F(s)$ との関係式
　これまでの記事で，以下の関係を導きました．今回は，これらの関係式を使って計算していきます．

自己共分散関数 $C(k)$ とランダムウォーク解析のゆらぎ関数 $F(s)$ の関係
　

$\displaystyle F^2(s) = \sum_{k=-(s-1)}^{s-1} (s-|k|) C(k)$

パワースペクトル $S(f_k)$ とランダムウォーク解析のゆらぎ関数 $F(s)$ の関係 ( $f_k=k/N$ です)

$\begin{equation} \displaystyle F^2(s) = \frac{1}{N}\sum_{k=0}^{N-1} \frac{\sin^2(\pi f_k s)}{\sin^2(\pi f_k)} S(f_k) \end{equation}$

2. 和を積分で近似
　区分求積法で，幅を無限小にしない場合，和と積分は等しくないけれど，近似になってます．つまり，

$\displaystyle \sum_{k=k_1}^{k_2} f(x[k]) \Delta x[k] \approx \int_{x[k_1]}^{x[k_2]} f(x) \, dx$

3. $\sin x \approx x$
　よく使う近似式です．マクローリン展開で，１次の項まで残してます．

$\alpha$ と $\beta$ の関係

　基本事項1にある

　 $\begin{equation} \displaystyle F^2(s) = \frac{1}{N}\sum_{k=0}^{N-1} \frac{\sin^2(\pi f_k s)}{\sin^2(\pi f_k)} S(f_k) \end{equation}$

において， $\Delta f_k = f_{k+1}-f_{k}=1/N$ なので，

　 $\begin{eqnarray} \displaystyle F^2(s) &=& \sum_{k=0}^{N-1} \frac{\sin^2(\pi f_k s)}{\sin^2(\pi f_k)} S(f_k) \Delta f_k \\ \end{eqnarray}$

と書けて， $N \to \infty$ で，

　 $\begin{eqnarray} \displaystyle F^2(s) &=& \int_{0}^{1} \frac{\sin^2(\pi f s)}{\sin^2(\pi f)} S(f) df \\ &=& \int_{-1/2}^{1/2} \frac{\sin^2(\pi f s)}{\sin^2(\pi f)} S(f) df \end{eqnarray}$

となります．積分区間は，ナイキスト周波数を意識して変更しました．

　上の積分に $S(f) \propto f^{-\beta}$ を代入して計算できればいいですが，私にはできないので， $\frac{\sin^2(\pi f s)}{\sin^2(\pi f)}$ を近似することで，無理やり計算します．

　 $f$ が，ものすごく小さいときは，上の基本事項３を使って，

　 $\displaystyle \frac{\sin^2(\pi f s)}{\sin^2(\pi f)}=\frac{(\pi f s)^2}{(\pi f)^2} = s^2$

とできそうです．

　 $f$ が， $1/2$ に近づいてくると， $s \gg 1$ ( $s$ がでかい)が影響してきて，分子の $\pi f s$ が小さい値とは言えなくなります．ですので， $\sin x \approx x$ は使えません．別の方法で近似する必要があります．
　
　 $\sin^2(\pi f s)$ は， $s$ が大きい値のとき， $f$ が変わると激しく振動します．それと， $\sin^2(\pi f s)$ は，最小値0，最大値1で振動するだけです．ということで，分子を平均値で近似してもいいですが面倒なので，最大値１で近似します．つまり， $f$ が $1/2$ に近い領域では，分子はほぼ一定で１くらい，分母は $\sin^2 (\pi f) = (\pi f)^2$ と近似します．

　 $\displaystyle \frac{\sin^2(\pi f s)}{\sin^2(\pi f)}=\frac{1}{(\pi f)^2}$

私の言っていることが信じられないと思いますので，グラフで近似結果を見てみます．下の図で，近似なしが青実線，近似が赤破線です．特徴は近似できています．

　この近似を使うことにして，２つの近似曲線の交点は， $\displaystyle f = \frac{1}{\pi s}$ です．したがって，

　 $\displaystyle \frac{\sin^2(\pi f s)}{\sin^2(\pi f)} = \left\{ \begin{array}{ll} \displaystyle s^2 & \displaystyle \left(0 < |f| < \frac{1}{\pi s} \right) \\ \displaystyle\frac{1}{(\pi f)^2} & \displaystyle \left(\frac{1}{\pi s} < |f| < \frac{1}{2} \right) \end{array} \right.$

と近似します． $\displaystyle F^2(s)$ の式に $S(f) \propto f^{-\beta}$ を代入し，さらにこの近似を使って，計算してみます．

　 $\begin{eqnarray} \displaystyle F^2(s) &=& \int_{-1/2}^{1/2} \frac{\sin^2(\pi f s)}{\sin^2(\pi f)} S(f) df \\ &=& 2 \int_{0}^{1/2} \frac{\sin^2(\pi f s)}{\sin^2(\pi f)} S(f) df \\ &=& 2 \int_{0}^{1/2} \frac{\sin^2(\pi f s)}{\sin^2(\pi f)} f^{-\beta} df \\ &\approx& 2 \int_{0}^{1/(\pi s)} \! f^{-\beta} s^2 \, df + \int_{1/(\pi s)}^{1/2} \! f^{-\beta} (\pi f)^{-2} \, df \end{eqnarray}$

この計算では， $-1 < \beta < 1$ のとき，積分の公式が使えるので，積分の部分を別々に計算すると，

$\begin{eqnarray} \int_{0}^{1/(\pi s)} \! f^{-\beta} s^2 \, df &=& s^2 \left[\frac{f^{-\beta+1}}{-\beta+1} \right]_0^{1/(\pi s)}\propto s^{\beta+1} \\ \int_{1/(\pi s)}^{1/2} \! f^{-\beta} (\pi f)^{-2} \, df &=& \frac{1}{\pi^2} \left[\frac{f^{-\beta-1}}{-\beta-1} \right]_{1/(\pi s)}^{1/2} \propto s^{\beta+1} \end{eqnarray}$

です ( $s$ の何乗になるかだけ気にしてください)．
ということで， $-1 < \beta < 1$ の範囲では，

　 $\begin{eqnarray} \displaystyle F^2(s) &\propto& s^{\beta+1} \end{eqnarray}$

になってます． $F^2(s) \propto s^{2 \alpha}$ と比較して，

　 $\displaystyle \alpha = \frac{\beta + 1}{2} \quad (-1 < \beta < 1)$

　 $\beta < -1$ では，積分で $f = 1/2$ の近くの面積が大きくなるので， $s$ の依存性は消えて，一定値，つまり， $F^2(s) \propto s^{0}$ ( $\alpha = 0$ )になっちゃいます．

　 $\beta > 1$ では，第一項の積分において， $f = 0$ に近い領域の面積がどでかくなるので，積分に関係ない $s^2$ のみが効いてきます．したがって， $F^2(s) \propto s^2$ ，すなわち， $\alpha = 1$ です．

以上をまとめますと，以下のようになります．

　 $\displaystyle \alpha = \left\{ \begin{array}{ll} \displaystyle \frac{\beta + 1}{2} & \displaystyle (-1 < \beta < 1) \\ \displaystyle 0 & \displaystyle (\beta < -1) \\ 1 & (1 < \beta) \end{array} \right.$

結局，使える $\alpha$ の範囲は，0から1です．狭すぎです．

$\alpha$ と $\gamma$ の関係

　これは，前に導きました．
chaos-kiyono.hatenablog.com

長期記憶過程に対応する， $0 < \gamma < 1$ の範囲では，

　 $\displaystyle \alpha = \frac{2-\gamma}{2} \quad (0 < \gamma < 1)$

負の長時間相関の， $1 < \gamma < 2$ の範囲でも，成り立ちますが，導出について今は思いつかないので，機会があれば説明します．
　

2022-06-29

【時系列解析】ランダムウォーク解析の解析的考察（その３：パワースペクトルとの関係）

時系列解析 1/fゆらぎパワースペクトルフラクタル長時間相関

今回も，ランダムウォーク解析についてです．ランダムウォーク解析 (fluctuation analysisと呼ばれることも)は，Natureに出版された論文でも使われています (以下のリンク)．
www.nature.com
とはいえ，ランダムウォーク解析は，今では実用的な時系列解析法とはいえません．それでも，教育的な題材としては役に立ちます．

以下では，ゆらぎ関数 $F(s)$ とパワースペクトル $S(f)$ の関係を導いていきます．前回の記事は，これです．
chaos-kiyono.hatenablog.com

基本事項の確認

1. ラグオペレータ $L$ を，時系列 $\{x[i]\}$ に作用させると，

$\displaystyle L x [i] = x[i-1]$ ，　 $\displaystyle L^m x [i] = x[i-m]$

となります．つまり， $L$ が1つかかるごとに，時間が１戻ります．

2. 時系列 $\{x[i]\}$ のフーリエ変換を $X(f)$ とします．時間の世界でのラグオペレータ $L$ をフーリエ変換すると， $e^{-{\bf i } 2 \pi f}$ になります．つまり，

$\displaystyle L x [i]$ をフーリエ変換すると， $\displaystyle e^{{\bf i} 2 \pi f} X(f)$

となります．フーリエ変換すると， $L$ が $e^{{\bf i} 2 \pi f}$ に変わると覚えてください．

3. 時系列 $\{x[i]\}$ のパワースペクトルの推定量 $S_x(f)$ は，長さ $N$ の時系列 $\{x[i]\}$ のフーリエ変換を $X(f_k)$ (ここで， $f_k = k/N$ )として，

$\displaystyle S_x(f_k) = \frac{|X(f_k)|^2}{N}$

です．

4. $S_x(f_k)$ (ここで， $f_k = k/N$ )の全面積は，時系列 $\{x[i]\}$ の2乗平均 (平均0なら分散)になります．

$\displaystyle \sum_{k=0}^{N-1} S_x(f_k) \Delta f_k=\sum_{k=0}^{N-1} \frac{S_x(f_k)}{N} = \left\langle x[i]^2 \right\rangle$

これは，パーセヴァルの等式 (Parseval's identity)と呼ばれる関係から来ています．

ランダムウォーク解析を周波数の世界で考える

　ランダムウォーク解析では，平均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}$

を求めるのでした．

　増分の計算の部分の計算から考えます．計算が楽になるように和をとる区間を１シフトします．つまり，「1からs」を，「0から (s-1)」にしちゃいました．そうすると，

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

です．これについて考えます．

ラグオペレータ $L$ を使えば，

　 $\begin{eqnarray} \Delta_{s} y[i] &=& x[i]+L^{-1} x[i] + \cdots + L^{-(s-1)} x[i] \\ &=& (1+L^{-1}+L^{-2}+\cdots + L^{-(s-1)} ) x[i] \\ &=& \frac{1-L^{s}}{1-L^{-1}} x[i] \end{eqnarray}$

です．最後の行になるのは，等比数列の和の公式を使ったからです．ちなみに等比数列の和の公式は，(初項－末項×公比)/(1－公比)でしたよね．

次に，両辺をフーリエ変換します． $\Delta_{s} y[i]$ のフーリエ変換を $Y_{\Delta s} (f_k)$ ， $x[i]$ のフーリエ変換を $X(f_k)$ と表せば，

　 $\begin{eqnarray} Y_{\Delta s} (f_k) &=& \frac{1-\left(e^{-{\bf i} 2 \pi f} \right)^{-s}}{1-\left(e^{-{\bf i} 2 \pi f} \right)^{-1}} X(f_k) \end{eqnarray}$

となります．これを使えば，パワースペクトルの関係式を求めることができます．つまり，「両辺の絶対値をとり，２乗して， $N$ でわる」の３ステップです．

　 $\begin{eqnarray} \left|Y_{\Delta s} (f_k)\right| &=& \left|\frac{1-\left(-e^{{\bf i} 2 \pi f_k} \right)^{-s}}{1-\left(-e^{{\bf i} 2 \pi f_k} \right)^{-1}} X(f_k) \right| \\ \left|Y_{\Delta s} (f_k)\right|^2 &=& \left|\frac{1-e^{-{\bf i} 2 \pi f_k s}}{1-e^{-{\bf i} 2 \pi f_k}} X(f_k) \right|^2 \\ \left|Y_{\Delta s} (f_k)\right|^2 &=& \left|\frac{1-e^{-{\bf i} 2 \pi f_k s}}{1-e^{-{\bf i} 2 \pi f_k}}\right|^2 \left|X(f_k) \right|^2 \\ \frac{\left|Y_{\Delta s} (f_k)\right|^2}{N} &=& \frac{\left|1-e^{-{\bf i} 2 \pi f_k s}\right|^2}{\left|1-e^{-{\bf i} 2 \pi f_k}\right|^2} \frac{\left|X(f_k) \right|^2}{N} \end{eqnarray}$

ここで， $\Delta_{s} y[i]$ のパワースペクトルを $S_{\Delta s} (f_k)$ ， $x[i]$ のパワースペクトルを $S(f_k)$ と表せば，

　 $\begin{eqnarray} \frac{\left|Y_{\Delta s} (f_k)\right|^2}{N} &=& \frac{\left|1-e^{-{\bf_k i} 2 \pi f_k s}\right|^2}{\left|1-e^{-{\bf i} 2 \pi f_k}\right|^2} \frac{\left|X(f_k) \right|^2}{N} S_{\Delta s} (f_k) &=& \frac{\left|1-e^{-{\bf i} 2 \pi f_k s}\right|^2}{\left|1-e^{-{\bf i} 2 \pi f_k}\right|^2} S(f_k) \\ \end{eqnarray}$

さらに，

　 $\begin{eqnarray} \frac{\left|1-e^{-{\bf i} 2 \pi f_k s}\right|^2}{\left|1-e^{-{\bf i} 2 \pi f_k}\right|^2} &=& \frac{(1-e^{-{\bf i} 2 \pi f_k s})(1-e^{{\bf i} 2 \pi f_k s})}{(1-e^{-{\bf i} 2 \pi f})(1-e^{{\bf i} 2 \pi f_k})} \\ &=& \frac{2-2 \cos(2 \pi f_k s)}{2 - 2 \cos (2 \pi f_k)} \\ &=& \frac{\sin^2(\pi f_k s)}{\sin^2(\pi f_k)} \\ \end{eqnarray}$

となるので，結局，

　 $\begin{eqnarray} S_{\Delta s} (f_k) &=& \frac{\sin^2(\pi f_k s)}{\sin^2(\pi f_k)} S(f_k) \\ \end{eqnarray}$

　ここで，ゆらぎ関数

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

に登場してもらい，上の基本事項の4を使えば，

　 $\begin{equation} \displaystyle F^2(s) = \frac{1}{N}\sum_{k=0}^{N-1} \frac{\sin^2(\pi f_k s)}{\sin^2(\pi f_k)} S(f_k) \end{equation}$

という関係式が得られます．これで，時系列のパワースペクトルとゆらぎ関数を解析的に結びつけることができました．正しいか心配なので，数値実験してみました．下の図は，同じ時系列 (H = 0.8の非整数ガウスノイズ)のランダムウォーク解析 (左)とパワースペクトル解析 (中)をした結果です．右の図は，ランダムウォーク解析の結果 (青丸)と，上の式を使ってパワースペクトル解析の結果を $F(s)$ に変換したもの (赤バツ)の比較です．わずかに違いましたが，丸め誤差の影響だと思います．ポイントは，「ランダムウォーク解析とパワースペクトル解析は同じ情報を持っている」ということです．とはいえ，ランダムウォーク解析の方が，ギザギザしないので，直線の当てはめはやりやすいです．一方で，パワースペクトル解析の結果は，平滑化してやる必要があるということです．

今日は，明日締め切りの原稿があるので，ここまでにしておきます．ミスを見つけたら教えてください．

おまけ

　作図に使用したRスクリプトはこれです．

# パッケージのインストールが必要な場合は以下を実行
# install.packages("longmemo")
# パッケージのロード
require(longmemo)
###################################
# パラメタの設定：ハースト指数(0 < H < 1)
H <- 0.8
###################################
# 時系列の長さ
N <- 10000
######################
# 解析するスケールs
s <- unique(round(exp(seq(log(1),log(N/10),length.out=20))))
n.s <- length(s)
######################
# 時系列の生成
x <- simFGN0(N,H)
x <- x - mean(x) # 平均0に
#####################
# ランダムウォーク解析
###
# 積分
y <- cumsum(x)
###
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))
}
freq <- (0:(N-1))/N
S.f <- abs(fft(x))^2/N
par(mfrow=c(1,3))
###
# log-logプロット
par(pty="s")
# log-logをとって直線をあてはめ
logs <- log10(s)
logFs <- log10(F.s)
fit <- lm(logFs~logs)
#
plot(logs,logFs,col=4,xlab="log10 s",ylab="log10 F(s)",las=1,
 main="Random Walk Analysis")
abline(fit,col=1,lty=2,lwd=2)
# log-logをとって直線をあてはめ
fit <- lm(log10(S.f[freq>0 & freq<1/10])~log10(freq[freq>0 & freq<1/10]))
#
plot(log10(freq[freq>0 & freq<1/2]),log10(S.f[freq>0 & freq<1/2]),"l",col=2,xlab="log10 f",ylab="log10 S(f)",las=1,
 main="Spectral Analysis")
abline(fit,col=1,lty=2,lwd=2)

F2.Sf <- c()
for(i in 1:n.s){
 F2.Sf[i] <- 0
 for(k in 1:N){
  fk <- k/N
  F2.Sf[i] <- F2.Sf[i]+sin(pi*fk*s[i])^2/sin(pi*fk)^2*S.f[k]
 }
}
plot(logs,logFs,col=4,xlab="log10 s",ylab="log10 F(s)",las=1,main="Comparison")
points(log10(s),log10(F2.Sf/N)/2,pch=4,col=2,xlab="log10 s",ylab="log10 F(s)",las=1)

2022-06-29

【時系列解析】ランダムウォーク解析の解析的考察（その２：自己相関関数との関係）

時系列解析フラクタル 1/fゆらぎ長時間相関

ランダムウォーク解析について前回に引き続き考えていきます．今回は自己相関（自己共分散）関数との関係です．前回の記事は以下です．

chaos-kiyono.hatenablog.com

ゆらぎ関数と自己相関（自己共分散）関数との関係

　ここでは．平均0，分散 $\sigma^2$ の弱定常過程の時系列 $\left\{x[i] \right\}$ を考えます．自己相関（自己共分散）関数は，

　 $C(k) = \left\langle x[i] x[i+k] \right\rangle$

です．折れたかっこ（ブラケット）は，期待値を表します．

　前回の記事に書いた，2乗平均の期待値 (ゆらぎ関数: fluctuation function)

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

の計算結果を見直して，書き直すと，

　 $\begin{eqnarray} F^2(s) &=& \left\langle \Delta_{s} y^2[i] \right\rangle = \displaystyle \left\langle \left(\sum_{j=1}^{s} x[i+j]\right)^2 \right\rangle \\ &=& s \left\langle x[i]^2 \right\rangle + 2 \sum_{k=1}^{s-1} (s-k) \left\langle x[i]x[i+k] \right\rangle \\ &=& \sum_{k=-(s-1)}^{s-1} (s-|k|) \left\langle x[i]x[i+k] \right\rangle \\ &=& \sum_{k=-(s-1)}^{s-1} (s-|k|) C(k) \\ \end{eqnarray}$

が成り立つことがわかります．公式的に示しておくと，こうなってます．

　 $\displaystyle F^2(s) = \sum_{k=-(s-1)}^{s-1} (s-|k|) C(k)$

長期記憶過程を解析すると

　長期記憶過程の時系列をランダムウォーク解析するとどうなるでしょうか．長期記憶過程では，自己相関関数が，

　 $C(k) \propto |k|^{-\gamma} \quad (0 < \gamma < 1)$

となります．これは， $0.5 < H < 1$ の非整数ガウスノイズに対応します．

　上の関係式を使うと，

　 $\begin{eqnarray} \displaystyle F^2(s) &=& \sum_{k=-(s-1)}^{s-1} (s-|k|) C(k) \\ &\propto& \sum_{k=-(s-1)}^{s-1} (s-|k|) |k|^{-\gamma} \\ &=& s \sum_{k=-(s-1)}^{s-1} |k|^{-\gamma} - \sum_{k=-(s-1)}^{s-1} |k|^{-\gamma+1} \\ &\propto& 2 s \sum_{k=1}^{s-1} k^{-\gamma} - 2 \sum_{k=1}^{s-1} k^{-\gamma+1} \\ &\approx& 2 s \int_1^{s} k^{-\gamma}\, dk - 2 \int_1^{s} k^{-\gamma+1}\, dk \\ &=& 2 \left(\frac{s^{2-\gamma}}{(1-\gamma)(2-\gamma)} -\frac{s}{1-\gamma} + \frac{1}{2-\gamma}\right) \end{eqnarray}$

と評価できます ( $\propto$ は比例する部分のみ残した， $\approx$ は近似を表します)． $s$ が大きくなると，第1項が支配的になるので，

　 $\begin{eqnarray} F^2(s) &\propto& s^{2-\gamma} \\ F(s) &\propto& \sqrt{s^{2-\gamma}} \\ F(s) &\propto& s^{(2-\gamma)/2} \end{eqnarray}$

となります． $\displaystyle F(s) \propto s^{\alpha}$ と比較して，以下の関係が成り立つことがわかります．

　 $\displaystyle \alpha = \frac{2-\gamma}{2} \quad (0 < \gamma < 1)$

伝統的な時系列解析は，弱定常中心主義，自己相関中心主義なので，「相関」が重んじられるような印象があります．しかし，上の式の $\gamma$ の値の範囲 ( $0 < \gamma < 1$ )は狭すぎです． $\alpha$ で言えば， $0.5 < \alpha < 1$ の範囲です．そのため，今回の議論では， $\alpha$ について，限られた知識しか教えてくれません．

次回は，パワースペクトルとの関係を説明します．いよいよ，ラグオペレータを使います．

2022-06-29

【時系列解析】ランダムウォーク解析の解析的考察（その１：白色ノイズの場合）

時系列解析短時間相関ブラウン運動フラクタルノイズ DFA DMA

ランダムウォーク解析のつづきです．今回は，この解析法について解析的に考えてみます．

chaos-kiyono.hatenablog.com

白色ノイズを分析すると

　白色ノイズは，無相関な時系列です．平均0，分散 $\sigma^2$ の白色ノイズの時系列 $\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}$

は，どんな感じでしょうか．

計算してみます．今回は，確率変数について大文字，小文字の区別はしません．全部，小文字です．弱定常であること，つまり，時間シフトしても，平均と自己共分散が変らないことを使えば，

　 $\begin{eqnarray} F^2(s) &=& \left\langle \Delta_{s} y^2[i] \right\rangle = \displaystyle \left\langle \left(\sum_{j=1}^{s} x[i+j]\right)^2 \right\rangle \\ &=& \left\langle x[i+1]^2 + x[i+2]^2 + \cdots + 2 \left( x[i+1] x[i+2]+ \cdots + \right) \right\rangle \\ &=& s \left\langle x[i]^2 \right\rangle + 2 \sum_{k=1}^{s-1} (s-k) \left\langle x[i]x[i+k] \right\rangle \end{eqnarray}$

となります．白色ノイズの場合， $\left\langle x[i]^2 \right\rangle = \sigma^2$ ， $\left\langle x[i]x[i+k] \right\rangle = 0$ ( $k \neq 0$ )なので，

　 $\begin{eqnarray} F^2(s) &=& \sigma^2 s \\ F(s) &=& \sigma s^{1/2} \\ \end{eqnarray}$

になります．ですので， $\log s$ に対して $\log F(s)$ をプロットして，傾きを推定すると $1/2=0.5$ になるはずです．

今回の計算には，弱定常の性質と，和の期待値の性質を使いました．

chaos-kiyono.hatenablog.com

次回の考察では，ラグオペレータに活躍してもらいますので，復習しておいてください．

chaos-kiyono.hatenablog.com

おまけ

図の作成に使ったRスクリプトです．

# 時系列の長さ
N <- 1000
# サンプル数
N.samp <- 100
######################
# 解析するスケールs
s <- unique(round(exp(seq(log(1),log(N/10),length.out=20))))
n.s <- length(s)
###
# ゆらぎ関数の準備
F.s <- c()
logFs.table <- data.frame()

par(mfrow=c(1,3))
for(j in 1:N.samp){
 #####################
 # サンプルの時系列 （白色ノイズ）
 # S(f) ~ f^0 (beta=0)の例
 x <- rnorm(N)
 x <- x - mean(x) # 平均0に
 #####################
 # ランダムウォーク解析
 ###
 # 積分
 y <- cumsum(x)
 ###
 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))
 }
 if(j == 1){
  par(pty="m")
  plot(1:N,x,"l",xaxs="i",xlim=c(0,N),main="サンプル時系列 (解析対象)",col=4,xlab="i")
  plot(1:N,y,"l",xaxs="i",xlim=c(0,N),main="積分時系列",col=4,xlab="i")
  logFs.table <- data.frame(log10(F.s))
 }else{
  logFs.table <- cbind(logFs.table,data.frame(log10(F.s)))
 }
}
logFs <- rowMeans(logFs.table)
# log-logをとって直線をあてはめ
logs <- log10(s)
fit <- lm(logFs~logs)
###
# log-logプロット
par(pty="s")
plot(logs,logFs,col=4,xlab="log10 s",ylab="log10 F(s)",las=1,
 main=paste("解析結果： 推定スロープ =",format(fit$coefficients[[2]],digits=2)))
# 直線フィット
abline(fit,col=2,lty=2)

2022-06-28

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

1/fゆらぎ DFA DMA R言語ノイズパワースペクトルフラクタルブラウン運動時系列解析

今回は，ランダムウォーク解析の導入です．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)

2022-06-28

【確率統計】確率変数の和の期待値と分散，積の期待値

統計数学 B4

ここでは，公式として，２つの確率変数 $X$ と $Y$ の和，積，分散の性質をまとめおきます．証明は，たくさんの方が与えてくれているのでそれを参照してください．以下の関係は，今後，ランダムウォークの軌道の性質を議論するために使います．

期待値の表現

　数学系の方々は，確率変数 $X$ の期待値を

　 ${\rm E} [ X]$ ， ${\rm E} (X)$ ， ${\rm E} X$

と表すことが多い気がします． ${\rm E}$ ではなく，特別なフォント $\mathbb{E}$ (TeXで\mathbb{E})を使うこともあります．

　物理系の私は，確率変数 $X$ の期待値を，ブラケットを使って

　 $\left\langle X \right\rangle$

とします．物理の偉大な先輩たちがこの表現を使ってきたので，私も論文を書くときはこの表現を使います．以下でもブラケットで期待値を表します．

期待値の定義

　離散確率変数 $X$ の確率質量関数を

　 $p(x) = {\rm Pr} (X=x)$

とすれば， $X$ の期待値は

　 $\displaystyle \left\langle X \right\rangle = \sum_{i} x_i p(x_i)$

です． $X$ の関数 $g(X)$ の期待値は

　 $\displaystyle \left\langle g(X) \right\rangle = \sum_{i} g(x_i) p(x_i)$

です．
　連続確率変数 $X$ の確率密度関数 $f(x)$ を

　 $f(x)\, dx = {\rm Pr} (x < X < x + dx)$

とすれば， $X$ の期待値は

　 $\displaystyle \left\langle X \right\rangle = \int_{-\infty}^{\infty} x \, f(x) \, dx$

です． $X$ の関数 $g(X)$ の期待値は

　 $\displaystyle \left\langle g(X) \right\rangle = \int_{-\infty}^{\infty} g(x) \, f(x) \, dx$

です．大文字と小文字の使い分けをマスターしてください．

２つの確率変数の和の期待値

　確率変数の和 $X + Y$ の期待値は，条件なしに，

　 $\displaystyle \left\langle X + Y \right\rangle = \left\langle X \right\rangle + \left\langle Y \right\rangle$

です．どんなときも，足し算，引き算は，分けることができます．

つまり，線形なので， $a$ ， $b$ を定数として，

　 $\displaystyle \left\langle a X + b Y \right\rangle = a \left\langle X \right\rangle + b \left\langle Y \right\rangle$

です．

２つの確率変数の積の期待値

　確率変数の積 $XY$ の期待値は，互いに独立なときのみ，

　 $\displaystyle \left\langle XY \right\rangle = \left\langle X \right\rangle \left\langle Y \right\rangle$

です．

２つの確率変数の分散

　確率変数 $X$ の分散を，

　 $V(X) = \left\langle (X - \left\langle X \right\rangle )^2 \right\rangle$

とします．
　確率変数の和 $X + Y$ の分散は，互いに独立なときのみ，

　 $V(X+Y) = V(X) + V(Y)$

です． $a$ ， $b$ を定数として，

　 $V(a X+b Y) = a^2 V(X) + b^2 V(Y)$

です．

参考になるページ

　証明は面倒なのでしません．以下のホームページがとても参考になるお思います．

2022-06-28

【フラクタル解析】１次元ランダムウォークの自己アフィン性

フラクタル DFA DMA 時系列解析長時間相関

　時系列のフラクタル解析に分類されるdetrended fluctuation analysis (DFA)とか，detrending moving average analysis (DMA)で何をやっているのかを理解するためには，無相関な時系列を積分すると，その積分時系列がどのような特性を持つようになるのかを理解することが必要です．

下の図は，＋１と－１が，ともに確率1/2で出る確率過程のサンプル時系列を生成し，それを積分したものです．積分時系列を２００本重ね書きしてあります．Rを使って，線が透明なグラフにしましたので，たくさん線が重なっているところほど色が濃くなっています．

上の３つの図は，横軸の幅が，それぞれ，100，1000，10000になっています．縦軸の幅は，100，1000，10000の平方根に比例するようにとりました（具体的には，平方根の6.4倍の幅です）．上の３つの図は，横幅が全然違うのに，目をぼやかして見ると全部同じようなグラフに見えませんか．１本のサンプル軌道を見ていても分かりづらいですが，たくさん重ねて描くと，軌道の広がりの様子がとても似ています．これが，ランダムウォークの自己アフィン性です．
自己アフィン性とは，縦横の倍率比が違う拡大をして，拡大の前後で形が同じになる性質です．ランダムウォークには統計的な意味で自己アフィン性があります．無相関な時系列を積分すると，自己アフィン性を持つ軌道になる．このことを理解することが，時系列のフラクタル解析の基礎を理解する出発点になります．正規分布に従う白色ノイズを積分しても，同じ特性をもつ軌道になります．今後，少しずつ説明していきます．

ランダムウォークの記事も参考にしてください．
chaos-kiyono.hatenablog.com

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

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

【時系列解析】ランダムウォーク解析の解析的考察（その4：スケーリング指数の関係）

時系列の長期記憶，長時間相関，フラクタル性

基本事項の確認

$\alpha$ と $\beta$ の関係

$\alpha$ と $\gamma$ の関係

【時系列解析】ランダムウォーク解析の解析的考察（その３：パワースペクトルとの関係）

基本事項の確認

ランダムウォーク解析を周波数の世界で考える

おまけ

【時系列解析】ランダムウォーク解析の解析的考察（その２：自己相関関数との関係）

ゆらぎ関数と自己相関（自己共分散）関数との関係

長期記憶過程を解析すると

【時系列解析】ランダムウォーク解析の解析的考察（その１：白色ノイズの場合）

白色ノイズを分析すると

おまけ

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

ランダムウォーク解析

数値実験

【確率統計】確率変数の和の期待値と分散，積の期待値

期待値の表現

期待値の定義

２つの確率変数の和の期待値

２つの確率変数の積の期待値

２つの確率変数の分散

参考になるページ

【フラクタル解析】１次元ランダムウォークの自己アフィン性