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

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

【時系列解析】パワースペクトルの定義の確認

時系列解析 パワースペクトル

弱定常が仮定できる場合のパワースペクトルの定義について，まとめておきます．これまでの私の記事では，時系列のフーリエ変換から計算するものと，時系列の自己共分散関数（自己相関関数）から計算するものを，私のかってな都合にあわせて使い分けていました．全部同じものですので，ここでもう一度整理したいと思います（厳密さよりわかりやすさを優先します）．面倒なので，サンプリング間隔は1とします．

前提

弱定常ということは，確率過程のサンプル時系列が無限の長さになったとしても，

  • 平均値が有限の値をとり，いつでも同じ

　（無限個のサンプル時系列を使って計算した各時刻の期待値がすべて同じ）

  • 分散が有限の値をとり，いつでも同じ

　（無限個のサンプル時系列を使って計算した各時刻の分散がすべて同じ）

  • 自己共分散（自己相関）関数が，2点間の時間差（ラグ）のみの関数で表される

を満たします．

時間無限，サンプル数無限の説明のための図．左はサンプル時系列，右はそのパワースペクトルの推定量．

真のパワースペクトル

弱定常な確率過程で，無限の長さのサンプル時系列を，無限本発生させて，期待値をとったパワースペクトルが「真のパワースペクトル」です（上図参照）．実際の時系列解析では，有限の長さの時系列で，1本のサンプル時系列からパワースペクトルを計算すると思います．これは，真のパワースペクトルではなく，それを推定したものです．長さNの1本の時系列から計算したパワースペクトルの推定量を \hat{S}_N (f)とすれば，真のパワースペクトルは，

　 \displaystyle S (f) = {\rm E} \left[ \lim_{N \to \infty} \hat{S}_N (f) \right]　
　（期待値は，たくさんの \hat{S}_N (f)を使って，fごとに求めます）

です．

パワースペクトルの推定量

まずは，時系列を離散フーリエ変換するパターンです．長さNの時系列 \{y [n]\}_{n=0}^{N-1}のパワースペクトルは，

　\displaystyle \hat{S}_N (f_k)= \frac{1 }{N} \left|\sum_{n=0}^{N-1} y [n] \, e^{- i 2 \pi f_k n} \right|^2　 \displaystyle \left(f_k = \frac{k}{N}\right)

です．

つぎは，まず長さNの時系列の自己共分散（自己相関）関数の推定量 \hat{C}[k]を計算し，その後に \hat{C}[k]のフーリエ変換を計算するパターンです．このパターンは，使える条件があることに注意してください．
\hat{C}[k]が， |k|の増加にともない，そこそこ早く0に近づく場合のみ，パワースペクトルは，

　 \displaystyle \hat{S}_N (f_k)= \sum_{n=-(N-1)}^{N-1} \hat{C}[n ] \, e^{- i 2 \pi f_k n} = \sum_{n=-(N-1)}^{N-1} \hat{C}[n ] \, \cos (2 \pi f_k n) 　 \displaystyle \left(f_k = \frac{k}{N}\right)

です．ここで，

　\displaystyle \hat{\mu} = \frac{1}{N} \sum_{n=0}^{N-1} y [ n]

　\displaystyle \hat{C} [ k ] = \frac{1}{N} \sum_{n=k}^{N-1} (y [ n ] - \hat{\mu}) (y [ n - k ] -\hat{\mu})

です． \hat{C}[k]の推定量は，正定値性とが，何か小難しいことを言う場合の形です． N-kで割っても，たいした違いはありません．

パワースペクトルを推定する方法は他にもあります．今後，それを理解するための知識を説明してから，追加していきます．

【Rで時系列解析】1/fゆらぎを積分すると傾き2だけ変化（数値実験）

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

前回は，時系列の積分がパワースペクトルにもたらす変化

「1/f^{\beta}ゆらぎを積分すると，1/f^{\beta+2}ゆらぎになる」

を解析的に示しました．

chaos-kiyono.hatenablog.com

今回は，数値実験です．

非整数ガウスノイズ（左）とその積分（右），下は対応するパワースペクトルのlog-logプロット．

使ったRスクリプトは以下です．パワースペクトルは，N.sampsで指定した数のサンプル時系列の平均で求めています．時系列は，1例を示しているだけです．

# longmemoがないときは以下を実行
# install.packages("longmemo")
# パッケージのロード
require(longmemo)
####################################
# ハースト指数(0 < H < 1)
H <- 0.6
# 時系列の長さ
n <- 2^10
####################################
# サンプル数
N.samps <- 100
####
S.x.mean <- numeric(n)
S.y.mean <- numeric(n)
for(i in 1:N.samps){
 ####################################
 # 非整数ガウスノイズの生成 
 x <- simFGN0(n,H)
 ####################################
 # 積分 (xの平均をゼロにする必要がある．なーんでか？)
 x <- x-mean(x)
 y <- cumsum(x)
 ####################################
 # フーリエ変換
 X <- fft(x)
 Y <- fft(y)
 # パワースペクトルの推定
 S.x <- abs(X)^2/n
 S.y <- abs(Y)^2/n
 #
 S.x.mean <- S.x.mean+S.x
 S.y.mean <- S.y.mean+S.y
}
S.x.mean <- S.x.mean/N.samps
S.y.mean <- S.y.mean/N.samps
################################################
f <- (0:(n-1))/n
log.S.x <- log10(S.x.mean[f > 0 & f <= 1/2])
log.S.y <- log10(S.y.mean[f > 0 & f <= 1/2])
y1 <- min(log.S.x,log.S.y)
y2 <- max(log.S.y)
log.f <- log10(f[f > 0 & f <= 1/2])
# 直線フィット
fit.x <- lm(log.S.x[log.f < -1]~log.f[log.f < -1])
fit.y <- lm(log.S.y[log.f < -1]~log.f[log.f < -1])
# プロット
par(mfrow=c(2,2))
par(pty="m")
plot(0:(n-1),x,"l",col=2,xlab="i",ylab="x[i]",xaxs="i")
plot(0:(n-1),y,"l",col=4,xlab="i",ylab="y[i]",xaxs="i")
par(pty="s")
plot(log.f,log10(S.x.mean[f > 0 & f <= 1/2]),col=2,pch=16,"l",lwd=2,ylim=c(y1,y2),las=1,
  xlab="log f",ylab="log Sx(f)",xaxs="i",main=paste("slope =",format(fit.x$coefficients[[2]],digit=3)),xaxs="i")
abline(fit.x,lty=2,col=1,lwd=2)
plot(log.f,log10(S.y.mean[f > 0 & f <= 1/2]),col=4,pch=16,"l",lwd=2,ylim=c(y1,y2),las=1,
  xlab="log f",ylab="log Sy(f)",xaxs="i",main=paste("slope =",format(fit.y$coefficients[[2]],digit=3)),xaxs="i")
abline(fit.y,lty=2,col=1,lwd=2)
########################################

ちなみに，上のRスクリプトで

 ####################################
 # 積分 (xの平均をゼロにする必要がある．なーんでか？)
 x <- x-mean(x)

の部分にある
x <- x-mean(x)
を消して，H>0.5のときの数値実験をしてみてください．積分時系列の傾きが-2になってしまうと思います．なぜでしょうか？

積分時系列のパワースペクトルが，予想に反して-2になってしまう例．

上の図を作成したRスクリプトです．

# longmemoがないときは以下を実行
# install.packages("longmemo")
# パッケージのロード
require(longmemo)
####################################
# ハースト指数(0 < H < 1)
H <- 0.9
# 時系列の長さ
n <- 2^10
####################################
# サンプル数
N.samps <- 100
####
S.x.mean <- numeric(n)
S.y.mean <- numeric(n)
for(i in 1:N.samps){
 ####################################
 # 非整数ガウスノイズの生成 
 x <- simFGN0(n,H)
 ####################################
 # 積分 (xの平均をゼロにする必要がある．なーんでか？)
 y <- cumsum(x)
 ####################################
 # フーリエ変換
 X <- fft(x)
 Y <- fft(y)
 # パワースペクトルの推定
 S.x <- abs(X)^2/n
 S.y <- abs(Y)^2/n
 #
 S.x.mean <- S.x.mean+S.x
 S.y.mean <- S.y.mean+S.y
}
S.x.mean <- S.x.mean/N.samps
S.y.mean <- S.y.mean/N.samps
################################################
f <- (0:(n-1))/n
log.S.x <- log10(S.x.mean[f > 0 & f <= 1/2])
log.S.y <- log10(S.y.mean[f > 0 & f <= 1/2])
y1 <- min(log.S.x,log.S.y)
y2 <- max(log.S.y)
log.f <- log10(f[f > 0 & f <= 1/2])
# 直線フィット
fit.x <- lm(log.S.x[log.f < -1]~log.f[log.f < -1])
fit.y <- lm(log.S.y[log.f < -1]~log.f[log.f < -1])
# プロット
par(mfrow=c(2,2))
par(pty="m")
plot(0:(n-1),x,"l",col=2,xlab="i",ylab="x[i]",xaxs="i")
plot(0:(n-1),y,"l",col=4,xlab="i",ylab="y[i]",xaxs="i")
par(pty="s")
plot(log.f,log10(S.x.mean[f > 0 & f <= 1/2]),col=2,pch=16,"l",lwd=2,ylim=c(y1,y2),las=1,
  xlab="log f",ylab="log Sx(f)",xaxs="i",main=paste("slope =",format(fit.x$coefficients[[2]],digit=3)),xaxs="i")
abline(fit.x,lty=2,col=1,lwd=2)
plot(log.f,log10(S.y.mean[f > 0 & f <= 1/2]),col=4,pch=16,"l",lwd=2,ylim=c(y1,y2),las=1,
  xlab="log f",ylab="log Sy(f)",xaxs="i",main=paste("slope =",format(fit.y$coefficients[[2]],digit=3)),xaxs="i")
abline(fit.y,lty=2,col=1,lwd=2)
########################################

【時系列解析】1/fゆらぎを積分するとパワースペクトルの指数が2だけ変化（解析的導出）

1/fゆらぎ 時系列解析 長時間相関

今回は，「1/f^{\beta}ゆらぎを積分すると，1/f^{\beta+2}ゆらぎになる」ことの説明です．

つまり，時系列 \left\{ x[i] \right\}_{i=0}^{N-1}のパワースペクトルS_x(f)が，低周波数側で

　 \displaystyle S_x(f) \propto \frac{1}{f^{\beta}}

に従うとき， \left\{ x[i] \right\}_{i=0}^{N-1}を積分した時系列

　 y[i] = \displaystyle \sum_{k=0}^{i-1} x[k]　　(i=1, 2, \cdots, N)

のパワースペクトルS_y(f)が，低周波数側で

　 S_y(f) \propto \displaystyle \frac{1}{f^{\beta+2}}

になることの説明です（下図参照）．逆に，差分をとると \displaystyle \frac{1}{f^{\beta-2}}になります．

これは， \frac{1}{f^{\beta}}ゆらぎを生成するためには，積分や差分の拡張して（2ずつの傾き変化でなく），非整数積分や非整数差分を導入すればいいな（自由に傾きをかえられる），というアイデアにつながる重要な事実です．

積分前後での時系列（上）とパワースペクトルの変化（下）
積分するとlog-logプロットの傾きが2変化

既に知っている1次自己回帰過程で考察

1次自己回帰過程

　y [n ]=a y [n-1] + w [n ]

で，a=1として，x [n] = w [n+1 ]とすれば，

　 y[i] = \displaystyle \sum_{k=0}^{i-1} x[k] 　　(i=1, 2, \cdots, N)

の形になってます．このとき， \left\{ x[i] \right\}_{i=0}^{N-1}は，白色ノイズなので，

　 \displaystyle S_x(f) \propto \frac{1}{f^{0}} = (定数)

です． y[i] のパワースペクトルについては， a < 1で， aが1に近いときは，高周波数側で

　\displaystyle S_y (f) \approx \frac{\sigma^2}{ 4 \pi^2 a f^2} \propto \frac{1}{f^{2}}

と近似できます（以前の記事で説明しました）．
chaos-kiyono.hatenablog.com

1次自己回帰過程のパワースペクトル

ですので，上の図のように， a \to 1とすれば，これがどんどん低周波数側に広がるので，最終的に，低周波数側でも，

　 S_y(f) \propto \displaystyle \frac{1}{f^{2}} = \frac{1}{f^{0+2}}

となることが期待できます．ですので，f^0からf^{-2}に変化することが分かります．

数学者のように細かいことを言うと，\{ y[i] \}のパワースペクトルは存在しません．なぜなら，n \to \inftyで時系列の分散が発散するので，パワースペクトル密度の全面積も発散するからです．

ということで，有限長の時系列を考える必要があります．

ここでの問題設定は，有限長の確率過程のサンプル時系列を積分して，そのパワースペクトルの推定値を計算するとどうなるか？です．

有限長時系列を積分したときのパワースペクトルの変化（解析的にやる）

　 \left\{ x[i] \right\}_{i=0}^{N-1}のフーリエ変換をX(f_k)とします．ここで，f_k = k/Nです．Nが奇数のときのみ正しい説明をします（偶数のときはf=1/2の扱いのみ違います）．

フーリエ変換の結果X(f_k)を使えば，

　 \displaystyle x[i] = \sum_{k=-(N-1)/2}^{k=(N-1)/2} A_k \cos (2 \pi f_k i + \theta_k)

と表せます．ここで，A_k = \frac{|X(f_k)|}{N}，\theta_k = {\rm arg} X(f_k)です．

y[i] = \displaystyle \sum_{j=0}^{i-1} x[j] の計算では，

　 \displaystyle y[i] = \sum_{j=0}^{i-1} \sum_{k=-(N-1)/2}^{k=(N-1)/2} A_k \cos (2 \pi f_k j + \theta_k) = \sum_{k=-(N-1)/2}^{k=(N-1)/2} \left(\sum_{j=0}^{i-1} A_k \cos (2 \pi f_k j + \theta_k) \right)

となるので，周波数成分ごとに和をとることにします．

\cosの和の公式

　\displaystyle \sum_{k=0}^{i-1} \cos\left( a k \right) = \frac{\sin \frac{a}{2}i}{\sin \frac{a}{2}} \cos \frac{a(i-1)}{2}

を使えば，

　\displaystyle \sum_{j=0}^{i-1} A_k \cos (2 \pi f_k j + \theta_k) =\frac{A_k}{2 \sin (\pi f_k)} \left\{\sin(2 \pi f_k i - \pi f_k + \theta_k ) +\sin (\pi f_k - \theta)\right\}

となります．この中で定数項は，パワースペクトルの傾きの議論には寄与しないので，

　\displaystyle \sum_{k=0}^{i-1} A_k \cos (2 \pi f_k i + \theta_k) \propto \frac{A_k}{2 \sin (\pi f_k)} \sin(2 \pi f_k i - \pi f_k + \theta_k )

を考えれば十分です．もっといえば，各周波数成分の振幅が，積分の前後で，

　 A_k

から，

　 \displaystyle \frac{A_k}{2 \sin (\pi f_k)}

になる結果を使うだけで，結論が導けます．周波数fが小さい領域を考えるので，ローラン展開（複素関数論ででてきます）をつかって，主要項のみ残せば（分母にsin x = xの近似でも構いません），

　 \displaystyle \frac{A_k}{2 \sin (\pi f_k)} \approx \frac{A_k}{2 \pi f_k}

となります．

ということで， S_x(f_k) \propto f_k^{-\beta} \propto A_k^2なので，

　 S_y(f_k) \propto \left( \frac{A_k}{2 \pi f_k} \right) ^2 = \frac{A_k^2}{\left(2 \pi f_k\right) ^2}= \frac{f_k^{-\beta}}{\left(2 \pi f_k\right) ^2} \propto f_k^{-\beta-2}

が導けました．パワースペクトルが存在するとかの話ではなく，有限長時系列でパワースペクトルの推定量の計算をすれば，これが成り立つということです．

【時系列解析】自己共分散，自己相関，自己相関係数，何が違う？

時系列解析

自己共分散 (autocovariance)，
自己相関 (autocorrelation)，
自己相関係数 (autocorrelation coefficient)，
この3つの違いを知ってますか？

ここで言いたいことは，これらの定義は分野によって違うので，名前にこだわるな，思い込みは捨てろ，教科書や論文ごとに最初に定義を確認しろ，ということです．これらの用語は，世界で統一されていないし，今後も統一されないと思います．

自己相関係数は，一般的な相関係数とのアナロジーで，[-1, 1]の値をとるものを呼ぶことが多いです．しかし，これを自己相関と呼ぶこともよくあります．平均なんて除かず，単純な積の期待値を自己相関関数と呼ぶことも，よくあります．最近私が自己共分散と呼んでいるものは，分野によっては自己相関と呼ばれます．今年度，学生と一緒に読んでいる本では，自己共分散と呼んでいるので，私も同じ用語を使っているだけです．私が論文を書くときには，「自己相関関数」と書くと思います．

【時系列解析】1次自己回帰過程の記憶（自己相関）はいつ消える？

時系列解析 パワースペクトル 短時間相関

私の研究室の学生は，私があまりにしつこく1次自己回帰過程の話をするのでうんざりしているかもしれません．さらに，しつこくやっていきます．

1次自己回帰過程

　y_{n}=a y_{n-1} + w_n　(-1 < a < 1，w_nは平均0，分散\sigma^2の白色ノイズ)

のパワースペクトルS(f)は，

　\displaystyle S(f)=\frac{\sigma^2}{1+a^2-2 a \cos (2 \pi f)}=\frac{\sigma^2}{(1-a)^2-2 a( \cos (2 \pi f)-1)}

です．今回は，このパワースペクトルのグラフの構造と時系列の記憶（自己相関）が消える時間の長さ（ラグ）の関係について説明します．

aの値を変えてグラフを描いたもの（両対数目盛だぞ）が以下です．ここでは，aは1に近いとします．

1次自己回帰過程のパワースペクトルの両対数プロット

パワースペクトルの特徴は，低い周波数で平らな丘があり，高い周波数側で下り坂があります．aを1に近づけていく（0.9 -> 0.99 -> 0.999）と，高い周波数側の下り坂の領域が広がり，f^{-2}に比例する構造が見えてきます．

1次自己回帰過程は，短期記憶過程の典型例です．時系列解析の分野では，「記憶」とは「自己相関（自己共分散）」のことです．短期記憶過程では，特徴的な時間の長さがあり，その時間以上では記憶が消える，つまり，自己相関が0とみなせます．では，この自己相関が0とみなせる特徴的な時間の長さを，どのようにすれば見積もることができるのでしょうか．

指数関数的に減衰する関数では，時定数とか，緩和時間とか呼ばれる特徴量があります．
　 x(t) \propto e^{-t/\tau_0}
の形のとき，時定数は\tau_0です．肩に乗っている部分が-1になるとして，tを求めれば，それが時定数です．

1次自己回帰過程の自己共分散関数C(k)は，

　\displaystyle C(k)=\frac{\sigma^2}{1-a^2} a^{k} \propto e^{(\log a) k}

ですので，時定数は

　\displaystyle -\frac{1}{\log a}

です（\displaystyle \log aはマイナスの値です）．時定数は値がほぼ0とみなせる時間だと，理解している人がいると思います．では，

　k > \displaystyle -\frac{1}{\log a}

で，記憶（自己相関）が消えるのでしょうか．違います．オーダーとしては正しい（10倍も違わない）ですが，無相関とみなせる時間としてもっと正確な範囲は，

　\displaystyle k > \frac{2 \pi \sqrt{a}}{1-a}，あるいは，\displaystyle -\frac{2 \pi}{\log a}

です．

なぜかを説明するために，パワースペクトルを考えます．パワースペクトルでは，自己相関が0とみなせる低周波数の領域で，平らな丘が現れます（両対数目盛でのグラフだぞ）．平らな丘は，白色ノイズと同じ構造です．これは，1次自己回帰過程のパワースペクトルに見られる低周波数側の平地のことです．ですので，パワースペクトルの低周波数側で，平らな丘に達する周波数を見積もってやれば，記憶が消えるスケールを見積もることができます．

解析的にこの周波数を計算するために，低周波数側の丘の部分と，高周波数側の坂道の近似関数を求めます．近似のために，aは1に非常に近いことを仮定します．さらに，

　 \displaystyle \cos x = 1 - \frac{x^2}{2!}+ \cdots

を使います．fは，1/2までの値をとるので，fの2乗まで考えます．

このとき，パワースペクトルは

　\displaystyle S(f) \approx \frac{\sigma^2}{(1-a)^2+ 4 \pi^2 a f^2}

となり，分母の(1-a)^2が，4 \pi^2 a f^2よりも大きいときは，

　\displaystyle S(f) \approx \frac{\sigma^2}{(1-a)^2}

分母の4 \pi^2 a f^2が，(1-a)^2よりも大きいときは，

　\displaystyle S(f) \approx \frac{\sigma^2}{ 4 \pi^2 a f^2}

と近似できそうです．数学的な厳密さとか言われても私は数学者ではないので，数値的に検証します．検証した結果が下の図の右上です．水色の破線（近似）と赤の線（近似なし）が一致しているので，正しそうです．

1次自己回帰過程のパワースペクトルと自己共分散関数の構造

坂から丘に変化する点を，上の2つの近似式の交点を使って見積もると，

　\displaystyle f = \frac{1-a}{2 \pi \sqrt{a}}

になります．これを，特徴的な周波数\displaystyle f_cとします．この逆数が特徴的な時間になるので，これを k_cと表せば，

　\displaystyle k_c = \frac{2 \pi \sqrt{a}}{1-a}

となります．今日は疲れたので説明しませんが，

　\displaystyle k_c = -\frac{2 \pi}{\log a}

もほぼ同じ値になります（正直，どうやって導いたか今は思い出せません）．どちらも，1-aでテーラー展開すると，

　\displaystyle \frac{2 \pi \sqrt{a}}{1-a} \approx (1-a) + \frac{1}{2} (1-a)^2+\cdots
　\displaystyle -\frac{2 \pi}{\log a} \approx (1-a) + \frac{1}{2} (1-a)^2+\cdots

のように2乗まで一致するので，ほとんど違いはありません．

上の図の右下の図は自己共分散のグラフです．垂直な灰色破線が時定数の位置，青色破線が\displaystyle k_c = 1/f_cの位置です．灰色破線の時間の長さ（ラグ）だと値が0には見えないけど，青色破線の時間ではほぼ0に見えます．

どのくらいの時間で記憶が消えるか，わかりましたか．その式を覚えるよりも，パワースペクトルの読み方を理解することが大切です．

おまけ

図の作成に使用したRスクリプト

a <- 0.98
par(mfrow=c(2,2))
par(pty="m",cex.axis=1.3,cex.lab=1.5,mar=c(5,5,2,2))
####################
curve(sig2/(1-2*a*cos(2*pi*x)+a^2),log="xy",col=2,lwd=2,n=1000,xlim=c(0.00001,1/2),ylim=c(0.2,10000),xaxs="i",xlab="f",ylab="PSD S(f)",xaxt="n",yaxt="n")
# 対数目盛を描く
axis(1,at=10^(-5:5)%x%(1:9),label=FALSE,tck=-0.02)
tmp<-paste(paste(sep="","expression(10^",-5:5,")"),collapse=",")
v.label<-paste("axis(1,las=1,at=10^(-5:5),label=c(",tmp,"),tck=-0.03)",sep="")
eval(parse(text=v.label))
axis(2,at=10^(-5:5)%x%(1:9),label=FALSE,tck=-0.02)
tmp<-paste(paste(sep="","expression(10^",-5:5,")"),collapse=",")
v.label<-paste("axis(2,las=1,at=10^(-5:5),label=c(",tmp,"),tck=-0.03)",sep="")
eval(parse(text=v.label))
####################
curve(sig2/(1-2*a*cos(2*pi*x)+a^2),log="xy",col=2,lwd=2,n=1000,xlim=c(0.00001,1/2),ylim=c(0.2,10000),xaxs="i",xlab="f",ylab="PSD S(f)",xaxt="n",yaxt="n")
# 対数目盛を描く
axis(1,at=10^(-5:5)%x%(1:9),label=FALSE,tck=-0.02)
tmp<-paste(paste(sep="","expression(10^",-5:5,")"),collapse=",")
v.label<-paste("axis(1,las=1,at=10^(-5:5),label=c(",tmp,"),tck=-0.03)",sep="")
eval(parse(text=v.label))
axis(2,at=10^(-5:5)%x%(1:9),label=FALSE,tck=-0.02)
tmp<-paste(paste(sep="","expression(10^",-5:5,")"),collapse=",")
v.label<-paste("axis(2,las=1,at=10^(-5:5),label=c(",tmp,"),tck=-0.03)",sep="")
eval(parse(text=v.label))
########################
abline(h=sig2/(1-a)^2,col=4,lwd=2,lty=2)
curve(sig2/(pi*x)^2/(4*a),add=TRUE,col=4,lwd=2,lty=2,xlim=c(0.00001,1/2))
#curve(sig2/(pi*x)^2/(4*a)*(pi*x)^2/sin(pi*x)^2,add=TRUE,col=5,lwd=2,lty=2,xlim=c(0.00001,1/2))
abline(v=(1-a)/(2*pi*sqrt(a)),col=4,lwd=2,lty=2)
#abline(v=-log(a)/(2*pi),col=1,lwd=2,lty=2)
####################################################
# 解析的に計算した自己共分散関数
curve(sig2/(1-a^2)*a^x,col=2,lwd=2,lty=1,xlim=c(1,10000),log="x",xaxs="i",xlab="k",ylab="Autocovariance C(k)",xaxt="n")
axis(1,at=10^(-5:5)%x%(1:9),label=FALSE,tck=-0.02)
tmp<-paste(paste(sep="","expression(10^",-5:5,")"),collapse=",")
v.label<-paste("axis(1,las=1,at=10^(-5:5),label=c(",tmp,"),tck=-0.03)",sep="")
eval(parse(text=v.label))
####################################################
# 解析的に計算した自己共分散関数
curve(sig2/(1-a^2)*a^x,col=2,lwd=2,lty=1,xlim=c(1,10000),log="x",xaxs="i",xlab="k",ylab="Autocovariance C(k)",xaxt="n")
axis(1,at=10^(-5:5)%x%(1:9),label=FALSE,tck=-0.02)
tmp<-paste(paste(sep="","expression(10^",-5:5,")"),collapse=",")
v.label<-paste("axis(1,las=1,at=10^(-5:5),label=c(",tmp,"),tck=-0.03)",sep="")
eval(parse(text=v.label))
abline(h=0,col=8,lwd=2,lty=2)
abline(v=-1/log(a),col=8,lwd=2,lty=2)
abline(v=-2*pi/log(a),col=4,lwd=2,lty=2)

【非線形力学系の数値解析】レスラー方程式のカオス

カオス 力学系 R言語

今日の講義資料です．散逸力学系のカオスアトラクター（下図）をＲを使って描画してみます．ただし，ここで説明する例は，ただの3次元プロットです．考える微分方程式はレスラー方程式です．パラメタは，カオスになる代表的な値を使います．方程式については，以下のWikipediaを参考にしてくださ（詳しく書いていないですが）．
レスラー方程式 - Wikipedia

レスラー方程式の時間発展の様子
レスラー方程式のカオスアトラクタ．ストレンジアトラクタとも呼ばれます．

Ｒを使って数値解を3次元プロット

以下のＲスクリプトの
　plt1 <- persp(x,y,z,theta=30,phi=10,.....)
に含まれる，thetaとphiの値を変えれば，見る角度を変えられます．thetaで水平面内の角度，phiで上下方向の角度を変えられます．

# 数値計算にdeSolveパッケージを使いました
#install.packages("deSolve")
require(deSolve)
### レスラーモデルの定義 ###################
Lorenz.eq <- function(t,state,parameters) {
 with(as.list(c(state, parameters)), {
  dxdt <- -y - z
  dydt <- x + a * y
  dzdt <- b + x * z - c * z
  list(c(dxdt, dydt, dzdt))
 })
}
# よく使われている制御変数の値を使用
control.parms <- c(a=0.2,b=0.2,c=5.7)
# 初期条件はトランジェントを捨てずにすませるためアトラクタに近い点をとった．
ini.cond <- c(x=-4.2842519245, y=-2.438439520, z=0.01965899)
# 時間ステップの設定
time.step <- seq(0, 200, 0.01)
#######################################
# deSolveを使って数値計算
PROF <- ode(y=ini.cond,times=time.step,func=Lorenz.eq,parms=control.parms)
###
# グラフの描画
x.range <- 12
x <- y <- c(-x.range, x.range)
z0 <- 0
z <- matrix(rep(z0,4), nrow=2)
###
# thetaとphiの値を変えれば見る角度が変る
par(mfrow=c(2,2))
plt1 <- persp(x,y,z,theta=30,phi=10,xlim=c(-x.range,x.range),ylim=c(-x.range,x.range),zlim=c(z0,2*x.range),expand=1,cex.lab=1.4)
lines(trans3d(PROF[,"x"],PROF[,"y"],PROF[,"z"], pmat=plt1), col = 2)
#
plt1 <- persp(x,y,z,theta=80,phi=10,xlim=c(-x.range,x.range),ylim=c(-x.range,x.range),zlim=c(z0,2*x.range),expand=1,cex.lab=1.4)
lines(trans3d(PROF[,"x"],PROF[,"y"],PROF[,"z"], pmat=plt1), col = 2)
#
plt1 <- persp(x,y,z,theta=20, phi=10,xlim=c(-x.range,x.range),ylim=c(-x.range,x.range),zlim=c(z0,2*x.range),expand=1,cex.lab=1.4)
lines(trans3d(PROF[,"x"],PROF[,"y"],PROF[,"z"], pmat=plt1), col = 2)
#
plt1 <- persp(x,y,z,theta=-45,xlim=c(-x.range,x.range),ylim=c(-x.range,x.range),zlim=c(z0,2*x.range), phi=30,expand=1,cex.lab=1.4)
lines(trans3d(PROF[,"x"],PROF[,"y"],PROF[,"z"], pmat=plt1), col = 2)

このスクリプトの実行結果が以下です．

Rスクリプトの実行結果


www.youtube.com

【時系列解析】２次自己回帰過程のサンプル時系列から自己共分散関数とパワースペクトルを推定

時系列解析 パワースペクトル R言語

本日のセミナーの課題は，２次自己回帰過程

　 y [n] = a_1 y [n-1] + a_2 y [n-2] + w[n ]

のサンプル時系列から自己共分散関数とパワースペクトルを推定することです． w[n]は，平均0，分散\sigma^2の白色ノイズです．とりあえず， a_1 = 0.9 \sqrt{3}， a_2 = -0.81でやってみて，他のパラメタでも調べてみてください．

これまでの知識を使えば簡単に実行できると思います．追加の課題として，白色ノイズ w[n]の部分を，単位インパルス

　 u[n] = \left\{\begin{array}{ll} 1 & (n=0) \\ 0 & (n \neq 0)\end{array} \right.

に置き換えた場合の振舞も数値的に調べてみてください．インパルスをいれるタイミングは0でなくても構いません。例えば，初期値をy[1] = 0，y[2] = 0として，ｗ[n]を{0, 0, 1, 0, 0, 0, 0, 0, 0, 0,... ,0}にしてください．

2次自己回帰過程のサンプル時系列．時間幅をだんだん短くしています．

課題の解答例

ということで，Rを使って， a_1 = 0.9 \sqrt{3}， a_2 = -0.81，分散1の白色ノイズ w[n]で，パワースペクトルと自己共分散関数を推定し，解析的な結果と比較します．
a1とa2がモデルのパラメタで，Nが時系列の長さ，N.sampsがサンプルの数です．自分で値を変更して，振舞をみてください．

# 2次自己回帰過程 (2nd-order autoregressive process) 
# y[n] = a1 y[n-1] + a2 y[n-2] + w[n]
# w[n]: white noise with zero mean and variance sig^2
####################
# パラメタ (parameter)
a1 <- 0.9*sqrt(3)
a2 <- -0.81
# w[n]の標準偏差sig^2
sig2 <- 1
# 時系列の長さ
N <- 2^13
####
y <- c()
####################
# 平均をとるサンプル時系列数
N.samps <- 2^6
####################
# 最大ラグの決定
m <- N-1
####################
# 周波数の定義
freq <- (0:m)/(2*m)
####################
PSD0 <- c()
####################
for(ii in 1:N.samps){
 ####################
 # 空の変数を事前に準備
 # 初期値をy[1]=w[1]とする
　# 過渡状態を少しまわす
 w <- rnorm(20,0,sqrt(sig2))
 y[1] <- w[1]
 y[2] <- w[2]
 for(n in 3:20){
  y[n] <- a1*y[n-1]+a2*y[n-2]+w[n]
 }
 #　ここが本番
 ####################
 # 白色ノイズの生成 
 w <- rnorm(N,0,sqrt(sig2))
 y[1] <- y[19]
 y[2] <- y[20]
 for(n in 3:N){
  y[n] <- a1*y[n-1]+a2*y[n-2]+w[n]
 }
 ######################
 # 自己共分散の推定 (sample autocovariance)
 a.cov <- acf(y,lag.max=m,type="covariance",plot=FALSE)$acf[,,1]
 # fftでフーリエ・コサイン変換
 PSD0 <- Re(fft(c(a.cov,rev(a.cov[-1]))))
 ####################
 # 平滑化
 PSD <- numeric(m+1)
 PSD[2:m] <- PSD0[2:m]*0.5+0.25*PSD0[1:(m-1)]+0.25*PSD0[3:(m+1)]
 PSD[1] <- (PSD0[1]+PSD0[2])/2
 PSD[m+1] <- (PSD0[m]+PSD0[m+1])/2
 ######################
 if(ii == 1){
  PSD.samps <- data.frame(PSD)
  AutoCov.samps <- data.frame(a.cov)
 }else{
  PSD.samps <- cbind(PSD.samps,data.frame(PSD))
  AutoCov.samps <- cbind(AutoCov.samps,data.frame(a.cov))
 }
}
f <- freq 
Spec <- rowMeans(PSD.samps)
AutoCov <- rowMeans(AutoCov.samps)
#################
# 図を描く
par(mfrow=c(1,3))
#################
# パワースペクトル（横軸のみ対数スケール）
plot(f[-1],Spec[-1],"l",col=4,log="x",main="Power spectrum density",xlab="f",xaxt="n",xaxs="i",ylab="",las=1)
abline(h=0,lty=2)
# 対数目盛を描く
axis(1,at=10^(-5:5)%x%(1:9),label=FALSE,tck=-0.02)
tmp<-paste(paste(sep="","expression(10^",-5:5,")"),collapse=",")
v.label<-paste("axis(1,las=1,at=10^(-5:5),label=c(",tmp,"),tck=-0.03)",sep="")
eval(parse(text=v.label))
########
# 解析的に計算したパワースペクトル
curve(sig2/(1-2*a1*(1-a2)*cos(2*pi*x)-2*a2*cos(4*pi*x)+a1^2+a2^2),add=TRUE,col=2,lwd=2,lty=2,n=1000,xlim=c(1/N,1/2))
legend("topleft", legend=c("Analytical PSD"), col=c(2), lty=c(1),lwd=c(2))
#######
# パワースペクトル（両対数スケール）
plot(f[-1],Spec[-1],"l",col=4,log="xy",main="Power spectrum density",xlab="f",xaxt="n",yaxt="n",xaxs="i",ylab="")
# 対数目盛を描く
axis(1,at=10^(-5:5)%x%(1:9),label=FALSE,tck=-0.02)
tmp<-paste(paste(sep="","expression(10^",-5:5,")"),collapse=",")
v.label<-paste("axis(1,las=1,at=10^(-5:5),label=c(",tmp,"),tck=-0.03)",sep="")
eval(parse(text=v.label))
axis(2,at=10^(-5:5)%x%(1:9),label=FALSE,tck=-0.02)
tmp<-paste(paste(sep="","expression(10^",-5:5,")"),collapse=",")
v.label<-paste("axis(2,las=1,at=10^(-5:5),label=c(",tmp,"),tck=-0.03)",sep="")
eval(parse(text=v.label))
########
# 解析的に計算したパワースペクトル
curve(sig2/(1-2*a1*(1-a2)*cos(2*pi*x)-2*a2*cos(4*pi*x)+a1^2+a2^2),add=TRUE,col=2,lwd=2,lty=2,n=1000,xlim=c(1/N,1/2))
legend("bottomleft", legend=c("Analytical PSD"), col=c(2), lty=c(1),lwd=c(2))
###########################################
# 自己共分散関数（横軸のみ対数スケール）
plot(1:m,AutoCov[-1],"l",log="x",col=4,main="Autocovariance",xlab="f",xaxt="n",xaxs="i",ylab="",las=1)
abline(h=0,lty=2)
# 対数目盛を描く
axis(1,at=10^(-5:5)%x%(1:9),label=FALSE,tck=-0.02)
tmp<-paste(paste(sep="","expression(10^",-5:5,")"),collapse=",")
v.label<-paste("axis(1,las=1,at=10^(-5:5),label=c(",tmp,"),tck=-0.03)",sep="")
eval(parse(text=v.label))
########
# 解析的に計算した自己共分散関数
acov <- c()
acov[1] <- sig2/((1-a2)^2-a1^2)*(1-a2)/(1+a2)
acov[2] <- a1/(1-a2)*acov[1]
for(i in 3:(m+1)){
 acov[i] <- a1*acov[i-1]+a2*acov[i-2]
}
lines(1:m,acov[-1],col=2,lwd=2,lty=2)
legend("topright", legend=c("Analytical Autocovariance"), col=c(2), lty=c(1),lwd=c(2))

このスクリプトを実行すれば，下の図が描かれます．

Rスクリプトの実行例．ただし，平均をとるサンプル数を増してあります．青がサンプル時系列からの推定，赤破線が解析的な結果です．

教科書の式(3.6)は，情報が抜けていて，それだけでは自己共分散関数C(k)は計算できません．たぶん，正しいものは，
C(0) = \frac{(1-a_2) \sigma^2}{(1+a_2) \left[ (1-a_2)^2-a_1^2 \right]}
C(1) = \frac{a_1}{1-a_2} C(0)
C(k) = a_1 C(k-1) + a_2 C(k-2)
です．来週までに，導いておいてください．

おまけ

単位インパルスを上記の過程に入れると，何が起こるか？というのが追加の課題でした．2次の自己回帰過程（特性方程式が複素数解を持つ場合）が，2階の線形微分方程式の減衰振動に対応することを感じてください．1次自己回帰過程は，1階の線形微分方程式の振舞や，過減衰（振動せずに0に漸近）に対応しています（この数値実験もしてみてください）．

2次自己回帰過程の単位インパルスに対する応答．赤が単位インパルスが入った箇所です．
# 2次自己回帰過程 (2nd-order autoregressive process) 
# y[n] = a1 y[n-1] + a2 y[n-2] + w[n]
####################
# パラメタ (parameter)
a1 <- 0.9*sqrt(3)
a2 <- -0.81
# 
y[1] <- 0
y[2] <- 0
# wに{0,0,1,0,0,....,0}を代入
w <- numeric(N) # N個0が代入される
w[3] <- 1 #3番目に1を代入
###
for(n in 3:N){
  y[n] <- a1*y[n-1]+a2*y[n-2] + w[n]
}
par(mfrow=c(1,1))
plot(1:N-1,y,"o",pch=16,col=4,xaxs="i",las=1,xlab="i",ylab="y[i]",lwd=2)
abline(h=0,col=8,lty=2)
lines(c(2,2),c(0,1),col=2,lwd=2)
points(c(2),c(1),pch=16,col=2)