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

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

【Rでシミュレーション】ブラウン運動のランダムウォークモデル

時系列解析 フラクタル ブラウン運動

講義用の参考資料です．今回は，ランダムウォークのサンプルパスを描きます．簡単化のためにx軸上の動きのみを考えます．

ランダムウォークのサンプルパス (p = 0.5, q=0.5)

　ランダムウォークは，ブラウン運動の確率モデルです．物理学におけるブラウン運動の重要性を知っていますか．今では奇跡の年と呼ばれる1905年に，26歳のアインシュタインは3本のノーベル賞級の論文を発表しました．その中の1本がブラウン運動（分子の大きさの決め方）についてのものです．その他の2本は，光電効果（ノーベル賞受賞）と特殊相対論（ E=mc^2の前段）です．特殊相対論の論文を博士論文として大学に提出したところ，理解されず受理されなかったので，ブラウン運動のネタの方を提出したそうです．

　アインシュタインの論文では，ブラウン運動を観測することで，分子の実在を示す理論が導かれています．当時はまだ分子の実在に否定的な意見がありました．なぜなら，「分子が存在し，その運動が古典力学で記述できる」と仮定してみても，マクロな熱現象を説明できないからです．背理法に基づけば，最初の仮定が間違っています．つまり，「分子」です．しかし，本当に間違っていたのは「分子」の部分ではなく，「古典力学」でした．当時はまだ量子力学が確立していなかったのです．

　後に，ペランが実験的にこの理論を検証して，分子の実在を示したことで1926年にノーベル賞をもらっています．興味があれば，米沢 富美子先生の本『ブラウン運動』（共立出版）を読んでみてください．ということで，ランダムウォークは物理学では重要なモデルです．さらに，時系列のフラクタル解析でも重要な基礎となりますので，ランダムウォークのイメージをもつために，各自で数値実験をしてみてください．

ランダムウォークの100本のサンプルパス（時間ステップ100まで）

　以下のRスクリプトで数値実験を実行できます．確率pで+1，確率q=1-pで-1になります．式で書くときは，確率を {\rm Pr}で表すことにして，

　 {\rm Pr}( \Delta X_i = +1)=p
　 {\rm Pr}( \Delta X_i = -1)=q = 1-p

であり，ランダムウォークの軌道は，

　 \displaystyle X_i = \sum_{j=1}^{i} \Delta X_j

です．今回は確率変数を大文字Xで書きました．実現値を表すときは小文字xを使います．

グラフの上段は， \Delta X_iの実現値 \Delta x_iで，下段はその累積（積分） x_iです．pの値を0.5から変えてみたり，Nを大きい値にしたりして，ランダムウォークの振る舞いを観察してみてください．

# +1が出る確率
p <- 0.5
# 総ステップ数
N <- 100
###############
i <- 1:N
q <- 1-p
# -1,+1の乱数の生成
tmp <- runif(N)
x <- (tmp > p)*2 - 1
###############
par(mfrow=c(2,1),mar=c(5,5,2,2))
plot(i,x,"p",pch=16,col=2,xaxs="i",las=1,xlab=expression(italic(i)),ylab=expression(paste(Delta, italic(x)[italic(i)])),
  lwd=2,xlim=c(0,N+1),main=paste("p = ",p,"; q = ",q,sep=""))
arrows(i,numeric(N),i,x,col=2,lwd=2,length=0)
abline(h=0,lty=2,col=8)
y <- cumsum(x)
plot(c(0,i),c(0,y),"o",pch=16,col=2,xaxs="i",las=1,xlab=expression(italic(i)),ylab=expression(italic(x)[italic(i)]),
  lwd=2,xlim=c(0,N+1),main=expression(paste(italic(x)[italic(i)], "=",Sigma, Delta, italic(x)[italic(i)])))
abline(h=0,lty=2,col=8)

このスクリプトを実行すると，下のような図が描かれます．

Rスクリプトの実行例

デモ用の動画も作っておきました．

www.youtube.com

【時系列解析】ラグオペレータ (lag operator)，後退オペレータ (backshift operator)の登場

時系列解析 数学

　自己回帰過程や移動平均過程を表現する際に便利な，ラグオペレータ (lag operator)，あるいは，後退オペレータ (backshift operator)と呼ばれるものを紹介します．どちらも名前が違うだけで同じものです．ラグオペレータはL，後退オペレータはBとすることが多いです．ここでは，Lを使うことにします．

ラグオペレータLのルール

　ラグオペレータLのルールは単純です．Lは時系列に対して使うので，時系列を\{x[i]\}とします．

  • Lと x[i]がかけ算でくっつくと，時間が１だけ戻る．つまり，

L x[i] = x[i-1]

  • L^j (Lのj乗)と x[i]がかけ算でくっつくと，時間がjだけ戻る．つまり，

　L^{j} x[i] = x[i-j]

  • L^{-1}と x[i]がかけ算でくっつくと，時間が1だけ進む．つまり，

　L^{-1} x[i] = x[i+1]

  • L^{-j}と x[i]がかけ算でくっつくと，時間がjだけ進む．つまり，

　L^{-j} x[i] = x[i+j]

  • Lは，変数だと思って計算して構わない．例えば，

　 (1+L)^2 = 1+2L+L^2

自己回帰過程をLを使って考えてみる

　一般形を書くのは面倒なので，私が大好きな1次自己回帰過程を考えます．

　 x[i] = a x[i-1] + w [i]　（w [i]は白色ノイズです）

Lを使うと，

　 x[i] = a L x[i] + w [i]

です．簡単です．

　これでは特にLの便利さを感じないと思いますので，ここでは，Lの素晴らしさを２つ紹介します．

Lを変数のように扱う

　上の式を変形します．Lを文字式同様に扱うと以下のようになります．

　 \begin{eqnarray} x[i] &=& a L x[i] + w [i] \\ x[i] - a L x[i] &=& w [i] \\ \left(1- a L \right) x[i] &=& w [i] \\ x[i] &=& \frac{1}{1- a L}w [i] \\ \end{eqnarray}

右辺は分数になって，最初に説明したルールでは解釈不能です．しかし，マクローリン展開の式

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

をまねてやると (収束半径なんて気にしません)，

　 \displaystyle \frac{1}{1- a L}=1+a L + a^2 L^2+a^3 L^3 + \cdots

となります．ということで，これを上の計算に使えば，

　 \begin{eqnarray} x[i] &=& \frac{1}{1- a L}w [i] \\ &=& \left( 1+a L + a^2 L^2+a^3 L^3 + \cdots \right) w [i] \\ &=& w [i] +a L w [i] + a^2 L^2 w [i]+a^3 L^3 w [i] + \cdots \\ &=& w [i] +a w [i-1] + a^2 w [i-2]+a^3 w [i-3] + \cdots \\ &=& \sum_{j=0}^{\infty} a^j w [i-j] \end{eqnarray}

となり，解釈可能な形になります．この記事では，移動平均過程を紹介していませんでしたが，自己回帰過程は，無限次の移動平均過程であることが分かります．Lを使わなくても，

　 \begin{eqnarray} x[i] &=& a x[i-1] + w [i] \\ &=& a \left(a x[i-2] + w [i-1] \right) + w [i] \\ &=& a^2 \left(a x[i-3] + w [i-2] \right) + a w [i-1] + w [i] \\ \end{eqnarray}

と，時間を戻していけば同じ結果になります．つまり，Lを変数のように考えて，マクローリン展開しても間違っていないようです．

Lは時間と周波数の世界の橋渡し役

　時系列を観測する時間iに支配される世界では，Lは時間を巻き戻す記号です．一方で，時系列をフーリエ変換すると周波数fが支配する世界になります．この時間の世界と周波数の世界の対応を考えると，

時間の世界のLは，周波数の世界ではe^{-{\bf i} 2 \pi f}になる．

ということです．{\bf i}は虚数単位です．　

時間の世界のx[i]は，周波数の世界ではX(f)になる．

ということで，

　 \begin{eqnarray} x[i] &=& \frac{1}{1- a L}w [i] \\ \end{eqnarray}

の両辺をフーリエ変換すると (Lの場所にe^{-{\bf i} 2 \pi f}を書いた)，

　 \displaystyle X(f) = \frac{1}{1- a e^{-{\bf i} 2 \pi f}} (w [i] のフーリエ変換)

これ以上書きません．ここからは，以前の記事を参考にしてもらうことにします．つまり，さらに計算を続ければ，正しいパワースペクトルが計算できます．

chaos-kiyono.hatenablog.com

【音】1/fノイズ（ピンクノイズ）は世界で人気

1/fゆらぎ ノイズ

　私の記事で扱いたいテーマの一つに，1/fノイズ（1/fゆらぎ，ピンクノイズ）があります．私は，1/fゆらぎの科学的意義を伝えたいのですが，世の中では「1/fゆらぎには，リラックス作用がある」と認識されていて，根強い人気があるようです．

　私のYoutubeチャンネルで，10日前に1/fノイズの音の動画をアップしました．１分半ほどの動画で，「ザ～～～ー～ー」と音が流れて，対応するパワースペクトルが描かれているだけです (以下の動画です)．


www.youtube.com


この動画が，私がこれまでアップしたものの中で，回数，時間とも最も再生されています．そもそも，私のYoutubeチャンネルは，講義の限定配信が中心で，一般公開している動画には，たいして人気のあるものはありません．何ヶ月たっても，視聴回数は数十程度で止まっています．

　そんな中，1/fノイズ（ピンクノイズ）の動画は，掲載後１０日経過した時点で，視聴回数3,887回，総再生時間45.4時間になってます (下の画像参照)．

Youtube動画の分析情報（動画アップから１０日後）

　過去48時間の視聴回数は，2,196回で，1時間あたり50回くらい視聴されています．下の図の青色の棒グラフは１時間ごとの視聴回数です．日本の昼夜関係なく，世界中で視聴されていることが想像できます．

過去48時間の視聴回数の推移

　何十万回，何百万回見られている動画とは比較になりませんが，私としては1/fノイズ（ピンクノイズ）が世界的人気を持つことを認識しました．ですので，今はノイズ関連の教材動画の作成に力をいれています．それと，今年度から私は，聴覚情報処理障害（Auditory Processing Disorder:APD)の研究にも取り組んでいるので，ノイズが我々の聴覚を通じた言語理解に与える影響を明らかにしたいと思います．

　1/fゆらぎのリラックス作用については，私は知りません．しかし，心臓の拍動リズム，呼吸リズム，連続血圧のゆらぎなど，体の中で観測されるゆらぎが，1/fゆらぎになっているのは事実で，1/fゆらぎであることが健康な状態であることも正しいと考えています．

【時系列解析】時系列の書き方などの慣例

時系列解析 数学

　時系列を表すときに，x(t)とか，\{x(t)\}とか，x_nとか，x[i]とか，いろいろな書き方があります．今回は，この書き方の慣例について説明します．絶対のルールがあるわけではありません．

連続時間は x(t)，整数の離散時間は下付 x_nか x[n]

　まずは，連続時間と離散時間の使い分けです．時系列（あるいは確率変数）は，時間順序をあわらす，tとか，nとかで，時刻が指定されます．

　高校の教科書では，

　t が実数（連続時間）をとるとき，x(t)
　t が整数をとるとき，x_t

になっています．高校で登場する漸化式（差分方程式）では，a_nのように下付になっています．大学の先生はそんなルールを気にしないので，大学によっては漸化式を扱った入試問題で a(n)となっていることもあります．

　高校を卒業すると，下付き，上付きの文字飾りをたくさん使うようになったり，下付文字だとそこに複雑な式を書きにくいので，時系列解析の分野では，

　n が整数をとるとき，x[n]

とします．

ですので，私の記事では，時間が整数で指定される場合に， x_nか， x[n]と書くことにしています．

時間が離散でもサンプリング間隔が \Delta tの場合は，背後に連続時間があることを意識して，

　x(0), x(\Delta t), x(2\Delta t), \cdots

とします．

中括弧は集合を表す．

　中括弧で囲んで，\{x(t)\}にするのは，たくさんの点の集合である気持ちを表しています．確率を表すときに， P(\{x_n\})としたりするのも，確率は事象の集合に対して決まるので，集合の場合は中括弧を付けます．

確率変数と実現値

　確率変数を大文字，その確率変数の実現値を小文字で表すルールがよく使われます．特に数学の場合はそのようにすることが普通だと思います．しかし，物理や工学では，大文字と小文字の使い分けはしないこともあります．

例えば，\{x(t)\}のフーリエ変換をX(f)で表すと，確率変数を大文字で表せなくなるので，時系列解析では，確率変数でも小文字のままにすることが多いです．

整数は，i, j, k

　\Sigmaで和をとるときに，i，j，kを使うことが多いです．これは，FORTRANというプログラム言語で，iからが整数型になっていたからではないかと私は思いますが，本当のことは分かりません．整数時間を表す変数についても，i, j, k, nあたりが，しっくりきます．

虚数単位はｊもある

　工学では，虚数単位をjにすることがあります．私はiを使います．整数を表すiやjとかぶるので，いやになることがあります．そんなときは，虚数単位を太字のiにしたり，他のフォントを使ったりします．

　ということで，基本は自由ですが，共通の表現ルールを使うと，理解しやすいと言うことがありますので，慣例に従うことも大切です．

【時系列解析】２次自己回帰過程のパワースペクトル（解析的導出）

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

今回から，本格的に２次自己回帰過程のパワースペクトルに入ります．自己回帰過程の理解において，最も重要なポイントの一つは

「n次自己回帰過程のパワースペクトルは，1次自己回帰過程と2次自己回帰過程のパワースペクトルの和になっている」

ということです (ほとんどの教科書には書いてありませんが，というか，書いてある教科書を私は知りません)．ですので，1次と2次自己回帰過程について理解すれば，自己回帰過程を完全にマスターした気分に近づけると思います．

今回は，前回紹介した，以下のテクニックを使います．
chaos-kiyono.hatenablog.com

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

それでは，２次自己回帰過程

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

を考えます．ここで， a_1， a_2は実数のパラメタ， \{x[n]\}は平均0，分散\sigma^2の白色ノイズです（以前はwとしたのを，xにしました）．

y [n]のフーリエ変換を， Y(f_k)，パワースペクトルの推定量をS_y(f_k)とします． x [n]についても，それぞれ， X(f_k)，S_x(f_k)とします．そして，S_x(f_k)=\sigma^2です (これは真のパワースペクトルですが)．Nを時系列の長さとすれば，f_k=k/Nです．上の式の両辺をフーリエ変換して，両辺の絶対値をとって２乗して，Nで割ると以下のようになります．
　 \begin{eqnarray} Y(f_k) &=& a_1 Y(f_k) \left( e^{{\rm \bf i} 2 \pi f_k} \right)^{-1} + a_2 Y(f_k) \left( e^{{\rm \bf i} 2 \pi f_k} \right)^{-2} + X(f_k) \\ \left(1 - a_1 e^{-{\rm \bf i} 2 \pi f_k} - a_2 e^{-2 {\rm \bf i} 2 \pi f_k} \right)Y(f_k) &=& X(f_k) \\ \left|1 - a_1 e^{-{\rm \bf i} 2 \pi f_k} - a_2 e^{-2 {\rm \bf i} 2 \pi f_k} \right|^2 \frac{ \left|Y(f_k)\right|^2}{N} &=& \frac{\left|X(f_k)\right|^2}{N} \\ \left|1 - a_1 e^{-{\rm \bf i} 2 \pi f_k} - a_2 e^{-2 {\rm \bf i} 2 \pi f_k} \right|^2_y(f_k) &=& S_x(f_k) \\ S_y(f_k) &=& \frac{S_x(f_k)}{\left|1 - a_1 e^{-{\rm \bf i} 2 \pi f_k} - a_2 e^{-2 {\rm \bf i} 2 \pi f_k} \right|^2} \\ S_y(f_k) &=& \frac{\sigma^2}{\left|1 - a_1 e^{-{\rm \bf i} 2 \pi f_k} - a_2 e^{-2 {\rm \bf i} 2 \pi f_k} \right|^2} \end{eqnarray}

複素数の絶対値の２乗が， |z|^2 = z \bar{z} であることに注意して，オイラーの公式

　 e^{{\rm \bf i} \theta} = \cos \theta + {\rm \bf i} \sin \theta

を使えば，２次自己回帰過程のパワースペクトル

　\begin{eqnarray} S_y(f_k) &=& \frac{\sigma^2}{1 + a_1^2 + a_2^2 - 2 a_1 (1 - a_2) \cos(2 \pi f_k) - 2 a_2 \cos( 4 \pi f_k)} \end{eqnarray}

がえられます．ここから，このパワースペクトルの性質を考えていきます．

パワースペクトルの例．a_1 = 0.9 \sqrt{3}，a_2 = -0.81．

パワースペクトルが1次自己回帰過程のパワースペクトルの和になる場合

２次自己回帰過程のパワースペクトルには，振動に対応するピークがあると思い込んでいませんか．違うときもあります．

　 a_1^2+4 a_2 \ge 0のときは，振動しません．
（注意：f=1/2のプラスマイナスを繰り返す振動はありえます）

ここでは，
　 \begin{eqnarray} S_y(f_k) &=& \frac{\sigma^2}{\left|1 - a_1 e^{-{\rm \bf i} 2 \pi f_k} - a_2 e^{-2 {\rm \bf i} 2 \pi f_k} \right|^2} \\ &=& \frac{\sigma^2}{\left(1 - a_1 e^{-{\rm \bf i} 2 \pi f_k} - a_2 \left(e^{-{\rm \bf i} 2 \pi f_k}\right)^2 \right)\left(1 - a_1 e^{{\rm \bf i} 2 \pi f_k} - a_2 \left(e^{{\rm \bf i} 2 \pi f_k}\right)^2 \right) } \end{eqnarray}

の分母に注目します．e^{{\rm \bf i} 2 \pi f_k} の部分を zで置き換えた式

1 - a_1 z - a_2 z^2

が実数の範囲で因数分解できて，

1 - a_1 z - a_2 z^2 = \left(1- \frac{z}{\lambda_1} \right)\left(1- \frac{z}{\lambda_2} \right)

とできるとします． a_1^2+4 a_2 > 0が仮定されています (重解はやりません)．\lambda_1と\lambda_2は， 1 - a_1 z - a_2 z^2 = 0 の解で

\lambda_1 = \displaystyle \frac{-a_1 + \sqrt{a_1^2+4 a_2} }{2 a_2}， \lambda_2 = \displaystyle \frac{-a_1 - \sqrt{a_1^2+4 a_2} }{2 a_2}

です．このとき，上の S_y(f_k)の式は， \lambda_1と\lambda_2を使って，
　 \begin{eqnarray} S_y(f_k) &=& \frac{\sigma^2}{\left(1 - \frac{ e^{-{\rm \bf i} 2 \pi f_k}}{\lambda_1} \right)\left(1 - \frac{ e^{-{\rm \bf i} 2 \pi f_k}}{\lambda_2} \right) \left(1 - \frac{ e^{{\rm \bf i} 2 \pi f_k}}{\lambda_1} \right)\left(1 - \frac{ e^{{\rm \bf i} 2 \pi f_k}}{\lambda_2} \right)} \\ &=& \frac{\sigma^2}{\left\{\left(1 - \frac{ e^{-{\rm \bf i} 2 \pi f_k}}{\lambda_1} \right)\left(1 - \frac{ e^{{\rm \bf i} 2 \pi f_k}}{\lambda_1} \right)\right\} \left\{ \left(1 - \frac{ e^{{\rm \bf i} 2 \pi f_k}}{\lambda_2} \right)\left(1 - \frac{ e^{-{\rm \bf i} 2 \pi f_k}}{\lambda_2} \right)\right\} } \\ &=& \frac{\sigma^2 \lambda_1^2 \lambda_2^2}{\left(1 + \lambda_1^2 - 2 \lambda_1 \cos \left(2 \pi f_k \right) \right) \left(1 + \lambda_2^2 - 2 \lambda_2 \cos \left(2 \pi f_k \right) \right)} \\ &=& \frac{\lambda_1^2 \lambda_2^3 \sigma^2}{\left(\lambda_1^2 \lambda_2-\lambda_1 \lambda_2^2-\lambda_1+\lambda_2\right) \left(1+\lambda_2^2-2 \lambda_2 \cos \left(\frac{2 \pi k}{N}\right)\right)}- \frac{\lambda_1^3 \lambda_2^2 \sigma^2}{\left(\lambda_1^2 \lambda_2-\lambda_1 \lambda_2^2-\lambda_1+\lambda_2\right) \left(1+\lambda_1^2-2 \lambda_1 \cos \left(\frac{2 \pi k}{N}\right)\right)} \end{eqnarray}

となります (計算間違いや，打ち間違いがあればすいません)．つまり，以下のような形になり，２つの1次自己回帰過程のパワースペクトルの足し算になってます (下図参照)．

\begin{eqnarray} S_y(f_k) &=& \frac{c_1}{1+\lambda_2^2-2 \lambda_2 \cos \left(\frac{2 \pi k}{N}\right)} + \frac{c_2}{1+\lambda_1^2-2 \lambda_1 \cos \left(\frac{2 \pi k}{N}\right)} \end{eqnarray}
(そのうち，一般的な分解法について詳しく説明します)．

パワースペクトルの例．a_1 = 0.9，a_2 = 0.7．破線が１次自己回帰過程のパワースペクトル．青と緑の破線を足すと，赤実線に一致します．

パワースペクトルが振動のピークを持つ場合

上の条件と違うとき，つまり，

　 a_1^2+4 a_2 < 0のとき，振動しちゃいます．

このとき， 4 a_2 < - a_1^2 \le 0なので， a_2 はマイナスの値です．

ここでは，最初に導いた２次自己回帰過程のパワースペクトル，

　\begin{eqnarray} S_y(f_k) &=& \frac{\sigma^2}{1 + a_1^2 + a_2^2 - 2 a_1 (1 - a_2) \cos(2 \pi f_k) - 2 a_2 \cos( 4 \pi f_k)} \end{eqnarray}

を使います． 0 < f < 1/2の領域で，パワースペクトルが振動に対応するピークを持つということは， S_y(f_k)の分母が下に凸なグラフになり，最小値を一つ持ちます (小さい値で割ると，値が大きいということだぞ)．ということで， a_2 がマイナスであることに注意して，分母のcosを含む部分の最小値を評価すると (\cos( 2 \pi f) の２次関数なので平方完成)，
　 \begin{eqnarray} - 2 a_1 (1 - a_2) \cos(2 \pi f_k) - 2 a_2 \cos( 4 \pi f_k) &=& - 2 a_2 \left\{2 \cos^2( 2 \pi f_k) - 1\right\}- 2 a_1 (1 - a_2) \cos(2 \pi f_k) \\ &=& - 4 a_2 \cos^2( 2 \pi f_k) - 2 a_1 (1 - a_2) \cos(2 \pi f_k) + 2 a_2 \\ &=& - 4 a_2 \left( \cos( 2 \pi f_k) - \frac{a_1 (a_2 - 1)}{4 a_2} \right)^2 + \frac{a_1^2 (a_2-1)^2}{4 a_2} + 2 a_2 \end{eqnarray}

となり，最小値（頂点）になる f_kは，

　 \begin{eqnarray} \cos( 2 \pi f_k) &=& \frac{a_1 (a_2 - 1)}{4 a_2} \\ 2 \pi f_k &=& \cos^{-1} \left(\frac{a_1 (a_2 - 1)}{4 a_2}\right) \\ f_k &=& \frac{1}{2 \pi}\cos^{-1} \left(\frac{a_1 (a_2 - 1)}{4 a_2}\right) \\ \end{eqnarray}

ということで，ピーク周波数は，

　\begin{eqnarray} f_k &=& \frac{1}{2 \pi}\cos^{-1} \left(\frac{a_1 (a_2 - 1)}{4 a_2}\right) \end{eqnarray}

です．英語のWikipediaには，結果のみ示してあります．
en.wikipedia.org

学生が，記事の更新が早すぎると不平を言いますが，理解してほしい内容にたどり着くまでには，このペースでも，１ヶ月くらいはかかるかもしれません．

おまけ

グラフを書くのに使用したRスクリプト．

ピークがある場合．

# ２次自己回帰仮定のパワースペクトル
f.ar2 <- function(x,a,b){
 return(1/(1+a^2+b^2-2*a*(1-b)*cos(2*pi*x)-2*b*cos(4*pi*x)))
}
f.ar11 <- function(x,a,b){
 return(a^2*b^3/(a^2*b-a*b^2-a+b)/(1+b^2-2*b*cos(2*pi*x)))
}
f.ar12 <- function(x,a,b){
 return(-a^3*b^2/(a^2*b-a*b^2-a+b)/(1+a^2-2*a*cos(2*pi*x)))
}
a1 <- 0.9*sqrt(3)
a2 <- -0.81
lambda1 <- (-a1+sqrt(a1^2+4*a2))/(2*a2)
lambda2 <-  (-a1-sqrt(a1^2+4*a2))/(2*a2)
# グラフを描いてみる
par(mfrow=c(1,2),cex.axis=1.2,cex.lab=1.5)
curve(f.ar2(x,a1,a2),xlim=c(0,1/2),n=1000,ylab="S(f)",xlab="f",col=2,lwd=2,las=1,xaxs="i")
if(a1^2+4*a2 < 0) abline(v=acos(a1*(a2-1)/(4*a2))/(2*pi),col=8,lty=2)
# 両対数目盛
curve(f.ar2(x,a1,a2),xlim=c(0.01,1/2),log="xy",n=1000,ylab="S(f)",xlab="f",col=2,lwd=2,xaxt="n",yaxt="n",xaxs="i")
if(a1^2+4*a2 < 0) abline(v=acos(a1*(a2-1)/(4*a2))/(2*pi),col=8,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))
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))
####################

１次自己回帰過程の重ね合わせになる場合．

# ２次自己回帰仮定のパワースペクトル
f.ar2 <- function(x,a,b){
 return(1/(1+a^2+b^2-2*a*(1-b)*cos(2*pi*x)-2*b*cos(4*pi*x)))
}
f.ar11 <- function(x,a,b){
 return(a^2*b^3/(a^2*b-a*b^2-a+b)/(1+b^2-2*b*cos(2*pi*x)))
}
f.ar12 <- function(x,a,b){
 return(-a^3*b^2/(a^2*b-a*b^2-a+b)/(1+a^2-2*a*cos(2*pi*x)))
}
a1 <- 0.9
a2 <- 0.7
lambda1 <- (-a1+sqrt(a1^2+4*a2))/(2*a2)
lambda2 <-  (-a1-sqrt(a1^2+4*a2))/(2*a2)
# グラフを描いてみる
par(mfrow=c(1,2),cex.axis=1.2,cex.lab=1.5)
curve(f.ar2(x,a1,a2),xlim=c(0,1/2),n=1000,ylab="S(f)",xlab="f",col=2,lwd=2,las=1,xaxs="i",ylim=c(0,3),yaxs="i")
curve(f.ar11(x,lambda1,lambda2),xlim=c(0,1/2),add=TRUE,n=1000,col=4,lwd=2,lty=2)
curve(f.ar12(x,lambda1,lambda2),xlim=c(0,1/2),add=TRUE,n=1000,col=3,lwd=2,lty=2)
if(a1^2+4*a2 < 0) abline(v=acos(a1*(a2-1)/(4*a2))/(2*pi),col=8,lty=2)
# 両対数目盛
curve(f.ar2(x,a1,a2),xlim=c(0.01,1/2),log="xy",n=1000,ylab="S(f)",xlab="f",col=2,lwd=2,xaxt="n",yaxt="n",xaxs="i",ylim=c(0.01,3))
if(a1^2+4*a2 < 0) abline(v=acos(a1*(a2-1)/(4*a2))/(2*pi),col=8,lty=2)
curve(f.ar11(x,lambda1,lambda2),xlim=c(0.01,1/2),add=TRUE,n=1000,col=4,lwd=2,lty=2)
curve(f.ar12(x,lambda1,lambda2),xlim=c(0.01,1/2),add=TRUE,n=1000,col=3,lwd=2,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))
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))
####################

【ノイズ】音で感じる確率過程

1/fゆらぎ ノイズ

1/fノイズやホワイトノイズを音として体験する動画資料です．視聴する際は音量を低めに設定してください．

1/fノイズ（ピンクノイズ）の音です．

www.youtube.com

ホワイトノイズの音です．

www.youtube.com

ホワイトノイズのハイパス信号です．カットオフ周波数が連続的に高周波数へ移動します．

www.youtube.com

1/fノイズ（ピンクノイズ）のサブバンドの音です．パワーが等しい区間に分割されています．最初は，ものすごい低音なので，スピーカーでは音が鳴らないかもしれません．逆に，高い音は年をとると聞こえないかもしれません．

www.youtube.com

クラシック音楽とパワースペクトルが同じノイズを作りました．

www.youtube.com

周波数が連続的に変化する純音です．

www.youtube.com

正弦波を重ねていくと，1/fノイズ（ピンクノイズ）になる動画です．

www.youtube.com

【時系列解析】自己回帰過程のパワースペクトルを解析的に求める

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

自己回帰過程のパワースペクトルを求めるテクニックの説明です．自己回帰過程の係数が与えられた状況から，パワースペクトルの式を導くことはとても簡単です．でも，大切なことは理屈を感じることです．

必要な基本知識

長さNの時系列\{x[i]\}_{i=0}^{N-1}のフーリエ変換とその逆変換は，

　\begin{eqnarray} x[i] &=& \frac{1}{N} \sum_{k=0}^{N-1} X(f_k) \, e^{{\rm \bf i} 2 \pi f_k i} \nonumber \\ X(f_k) &=& \sum_{i=0}^{N-1} x[i] \, e^{- {\rm \bf i} 2 \pi f_k i} \nonumber \end{eqnarray}

です． f_k = k/Nで，{\rm \bf i} は虚数単位です．

　\begin{eqnarray} x[i] &=& \frac{1}{N} \sum_{k=0}^{N-1} X(f_k) \, e^{{\rm \bf i} 2 \pi f_k i} \nonumber \end{eqnarray}

は恒等式なので， i に， i -jを代入すると，

　\begin{eqnarray} x[i-j] &=& \frac{1}{N} \sum_{k=0}^{N-1} X(f_k) \, e^{{\rm \bf i} 2 \pi f_k (i-j)} \nonumber \\ &=& \frac{1}{N} \sum_{k=0}^{N-1} X(f_k) \, e^{{\rm \bf i} 2 \pi f_k i} \, e^{-{\rm \bf i} 2 \pi f_k j} \nonumber \\ &=& \frac{1}{N} \sum_{k=0}^{N-1} X(f_k) \, e^{{\rm \bf i} 2 \pi f_k i} \, \left( e^{-{\rm \bf i} 2 \pi f_k} \right)^j \nonumber \\ &=& \frac{1}{N} \sum_{k=0}^{N-1} \left\{X(f_k) \left( e^{{\rm \bf i} 2 \pi f_k} \right)^{-j} \right\} e^{{\rm \bf i} 2 \pi f_k i} \nonumber \end{eqnarray}

フーリエ変換の式と見比べてみると，

となってます．つまり，x[i-j]の，-jの部分に注目して， \left( e^{{\rm \bf i} 2 \pi f_k} \right)の-j乗を， X(f_k)にかけてあるだけです．この部分は，各周波数成分の波の位相が周波数f_kに依存してシフトすることを表しています．iとi-jでは，位置がほんの少しずれているので，振幅は同じで，波の位置（位相）だけずれます．

練習すると，

最後に，長さNの時系列\{x[i] \}_{i=0}^{N-1}のフーリエ変換を， X(f_k) とすれば，パワースペクトルの推定量\hat{S}_x (k_f)は，

\displaystyle \hat{S}_x (f_k) = \frac{|X(f_k)|^2}{N}

以上で準備は終わりです．

差分方程式の両辺をフーリエ変換してパワースペクトルを計算

ということで，上の知識を使えば差分方程式の両辺をフーリエ変換して，フーリエ変換の関係式を導けます．パワースペクトルも簡単に求められます．

例えば，x[i]のフーリエ変換が X(f_k)で，y[i]のフーリエ変換が Y(f_k)として，1次自己回帰過程のような差分方程式，

　 y[i] = a y[i-1] + x [i]

を考えます．この両辺をフーリエ変換すると，

　 \begin{eqnarray} Y(f_k) &=& a Y(f_k) \left( e^{{\rm \bf i} 2 \pi f_k} \right)^{-1} + X (f_k) \\ \left(1- a e^{-{\rm \bf i} 2 \pi f_k} \right) Y(f_k) &=& X(f_k) \end{eqnarray}

両辺をパワースペクトルの形にするために，両辺の絶対値の２乗をとって，Nでわると，

　 \begin{eqnarray} \frac{\left|\left(1- a e^{-{\rm \bf i} 2 \pi f_k} \right) Y(f_k)\right|^2}{N} &=& \frac{\left|X(f_k)\right|^2}{N} \\ \left| 1- a e^{-{\rm \bf i} 2 \pi f_k} \right|^2 \frac{ \left|Y(f_k)\right|^2}{N} &=& \frac{\left|X(f_k)\right|^2}{N} \end{eqnarray}

x[i]のパワースペクトルの推定量は， S_x(f_k) = \frac{ \left|X(f_k)\right|^2}{N}で，
y[i]のパワースペクトルの推定量は， S_y(f_k) = \frac{ \left|Y(f_k)\right|^2}{N}なので，
パワースペクトルの関係式

　 \begin{eqnarray} \left| 1- a e^{-{\rm \bf i} 2 \pi f_k} \right|^2 S_y(f_k) &=& S_x(f_k) \\ S_y(f_k) &=& \frac{1}{\left| 1- a e^{-{\rm \bf i} 2 \pi f_k} \right|^2}S_x(f_k) \end{eqnarray}

がえられます． x [i]が白色ノイズ w [i]（分散\sigma^2）のときが，１次自己回帰過程

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

だったので， w [i]のパワースペクトル S_x [i] = S_w [i] = \sigma^2を上の式に代入すれば，

　 \begin{eqnarray} S_y(f_k) &=& \frac{\sigma^2}{\left| 1- a e^{-{\rm \bf i} 2 \pi f_k} \right|^2} \end{eqnarray}

となります．複素数の絶対値の２乗が， |z|^2 = z \bar{z} であることに注意して，オイラーの公式

　 e^{{\rm \bf i} \theta} = \cos \theta + {\rm \bf i} \sin \theta

を使えば，今まで何度もお世話になってきた，１次自己回帰過程のパワースペクトル

　\displaystyle S(f)=\frac{\sigma^2}{\left( 1- a e^{-{\rm \bf i} 2 \pi f_k} \right) \left( 1- a e^{ {\rm \bf i} 2 \pi f_k} \right)} = \frac{\sigma^2}{1+a^2-2 a \cos (2 \pi f)}

になっていることが分かります．

白色ノイズの積分過程も計算してみる

数学的にパワースペクトルが存在するとか，しないとか，私はそんなこと気にしないので，白色ノイズの積分過程

　 y[i] = y[i-1] + w [i]

を考えます．両辺をフーリエ変換して，絶対値をとって，２乗して，Ｎで割れば，

　 \begin{eqnarray} S_y(f_k) &=& \frac{1}{\left| 1- e^{-{\rm \bf i} 2 \pi f_k} \right|^2}S_w(f_k) \\ S_y(f_k) &=& \frac{\sigma^2}{2 - 2 \cos (2 \pi f)} \\ \end{eqnarray}

になります．fが小さいとして，cosの部分を展開すれば，

　 \begin{eqnarray} S_y(f_k) &=& \frac{\sigma^2}{2 - 2 \cos (2 \pi f)} \approx \frac{\sigma^2}{4 \pi^2 f^2} \\ \end{eqnarray}

となり，前にも見た結果がえられます．

これで，何次の自己回帰過程でもパワースペクトルを計算できるようになりました．今回の話は，ラグオペレータとか，後退オペレータとか呼ばれる方法が，何でうまく機能するのかを理解するのに役立ちます．この話は，近いうちにやります．