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

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

【インパルス応答】t = 0 に閉じ込められた物語

前の記事（以下のリンク）で，デルタ関数，あるいは，単位インパルス関数がどこに立っているのかにこだわりました．今回は，線形システムのインパルス応答を考えて，私がなぜインパルスが立っている位置（時刻）にこだわるのかを説明します．
chaos-kiyono.hatenablog.com


　インパルス応答では，下の図の上段のように，t=0の瞬間にいろんなことが起こりすぎてややこしいです．私としては，まず下の図の下段のように考えて，それをt=0に閉じ込めたとみなしたいです．

インパルス応答

　以下では具体例として，

　\begin{eqnarray} \frac{dy(t)}{dt} + a \, y(t) &=& \delta_{\varepsilon}(t) \end{eqnarray}

を考えます．ここで，aは正の定数，\delta_{\varepsilon}(t)は下の図のように定義します．

\delta_{\varepsilon}(t)の定義．右の図は，t>0に置いてt=0に近づける場合．

　この微分方程式でa=3として，\delta_{\varepsilon}(t)の\varepsilonを変えて，数値的に解いた結果が下の図です．ピンクの部分が\delta_{\varepsilon}(t)，青の線が応答（出力）y(t)で，\varepsilonを0に近づけています．

入力\delta_{\varepsilon}(t) (ピンク)に対するシステムの応答（数値計算結果）．\varepsilonを0に近づけた結果．

　通常は，\varepsilon \to 0とした結果を考えています．以下では，\varepsilon \to 0としたものを，\delta(t)として，
　
　\begin{eqnarray} \frac{dy(t)}{dt} + a \, y(t) &=& \delta(t), \quad y(0) = 0 \end{eqnarray}

を解いてみたいと思います．解き方はいくつかあります．私が思いつくのは3通りの方法です： (1) 同次方程式の解を求めてから定数変化法で特解を求める; (2) 積分因子を掛けて解く; (3) ラプラス変換で解く．

定数変化法で解く

　\begin{eqnarray} \frac{dy(t)}{dt} + a \, y(t) &=& \delta(t) \end{eqnarray}

は，非同次（非斉次）方程式（外部からの入力を表すtの関数を含む式）です．この方程式の一般解は，同次（斉次）方程式（外部からの入力を表すtの関数がない式）

　\begin{eqnarray} \frac{dy(t)}{dt} + a \, y(t) &=& 0 \end{eqnarray}

の一般解y_0 (t)と，右辺に\delta(t)を含む非同次方程式の特解y_{1} (t)の和で与えられます．なぜなら，y_0 (t)とy_{1} (t)は，それぞれ，

　\left\{\begin{array}{l} \displaystyle \frac{dy_0(t)}{dt} + a \, y_0(t) &=& 0 \\ \displaystyle \frac{dy_1(t)}{dt} + a \, y_1(t) &=& \delta(t) \end{array}\right.

を満たすので，辺々足して，

　\displaystyle \frac{d}{dt} \left(y_0(t) + y_1(t) \right) + a \left(y_0(t) + y_1(t) \right) = \delta(t)

となるからです．結局，

　 y(t) = y_0(t) + y_1(t)

も元の微分方程式の解になることが分かります．

　ここから，実際に解いていきます．まず，同次方程式の一般解を求めると，

　\begin{eqnarray} \frac{dy(t)}{dt} &=& - a \, y(t) \\ \int \frac{1}{y(t)} \, dy(t) &=& - \int a \, dt \\ \log \left| y(t) \right| &=& - a t + C_0 \\ y(t) &=& \pm e^{C_0}\, e^{- a t} \\ y(t) &=& C e^{-at} \end{eqnarray}
　
となります．ここで，C_0とCは定数で，C = \pm e^{C_0}です．

　定数変化法では，上の一般解の定数部分Cを，tの関数C(t)に変えちゃいます．つまり，y(t) = C(t) e^{-at}を非同次方程式に代入し，C(t)を求めます．

　\begin{eqnarray} \frac{dy(t)}{dt} + a \, y(t) &=& \delta (t) \\ e^{-at} \frac{dC(t)}{dt} - a C(t) e^{-at} + a C(t) e^{-at} &=& \delta(t) \\ \frac{dC(t)}{dt} &=& \delta(t) e^{at} \\ \int_{0}^{t} dC(\tau) &=& \int_{0}^{t} \delta(\tau) \, e^{a \tau} \, d\tau \\ C(t) - C(0) &=& \left\{\begin{array}{lc} 0 & (t \le 0) \\ 1 & (t > 0) \end{array}\right. \end{eqnarray}

私の流儀では，\delta(t)の面積を，t>0の領域に置いているので，

　\begin{eqnarray}\int_{0}^{\infty} \! \delta(t) \, dt &=& 1 \\ \int_{0}^{\infty} \! f(t) \, \delta(t) \, dt &=& f(0) \end{eqnarray}

になります．さらに，任意の\varepsilon > 0に対して，

　\begin{eqnarray}\int_{0}^{\varepsilon} \! \delta(t) \, dt &=& 1 \\ \int_{0}^{\varepsilon} \! f(t) \, \delta(t) \, dt &=& f(0) \end{eqnarray}

です．

　上のC(0)は，一般解に対応するので無視すれば，特解が

　y(t) = C(t) \, e^{-at} = e^{-at} \quad (t > 0)

と求まります．したがって，一般解は，Cを定数として，

　\begin{eqnarray} y(t) =\left\{\begin{array}{lc} C e^{at} & (t \le 0) \\ C e^{at} + e^{-at} & (t > 0) \end{array}\right. \end{eqnarray}

です．初期条件 y(0)=0を満たすようにCを設定すれば (C=0)，

　\begin{eqnarray} y(t) =\left\{\begin{array}{lc} 0 & (t \le 0) \\ e^{-at} & (t > 0) \end{array}\right. \end{eqnarray}

となります．インパルス \delta(t)に対する応答は，t > 0で生じています．t = 0にインパルスを与えたと考えると，変な感じがします．ですので，私は t > 0で，tが0に近い時刻でインパルスを入れたとみなしたいのです．

積分因子を使って解く

　\begin{eqnarray} \frac{dy(t)}{dt} + p(t) \, y(t) &=& q(t) \end{eqnarray}

の形の微分方程式は，両辺に

　\displaystyle \exp \left( \int p(t)\, dt \right)

をかけることで，解くことができます．ここで，上の不定積分の定数項は無視してください．これは，左辺を積の微分の形にするアイデアです．

　具体的に，

　\begin{eqnarray} \frac{dy(t)}{dt} + a \, y(t) &=& \delta(t) \end{eqnarray}

を解いてみます．この場合， p(t) = aです．

　\displaystyle \exp \left( \int a \, dt \right) = e^{at}

　\displaystyle e^{at}を微分方程式の両辺にかけると，

　\begin{eqnarray} e^{at} \frac{dy(t)}{dt} + a e^{at} \, y(t) &=& \delta(t) \, e^{at} \\ \frac{d}{dt} \left( e^{at} \, y(t) \right) &=& \delta(t) \, e^{at} \end{eqnarray}

となります．左辺が，積の微分の形になっていることに気づいてください．これを解くと，

　\begin{eqnarray} \frac{d}{dt} \left( e^{at} \, y(t) \right) &=& \delta(t) \, e^{at} \\ \int_0^{t} d \left( e^{a \tau} \, y(\tau) \right) &=& \int_0^{t} \delta(\tau) \, e^{a \tau} \, dt \\ e^{a t} \, y(t) - y(0) &=& \left\{\begin{array}{lc} 0 & (t \le 0) \\ 1 & (t > 0) \end{array}\right. \\ y(t) &=& \left\{\begin{array}{lc} 0 & (t \le 0) \\ e^{-at} & (t > 0) \end{array}\right. \end{eqnarray}

となります．

ラプラス変換で解く

　ここでは，ラプラス変換については詳しく説明しません．ラプラス変換と，その逆変換の公式の導出は，自分で調べてください．

　今回の微分方程式

　\begin{eqnarray} \frac{dy(t)}{dt} + a \, y(t) &=& \delta(t), \quad y(0) = 0 \end{eqnarray}

をラプラス変換で解いてみます．

　以下では，y(t) のラプラス変換を Y(s) とします．微分方程式の両辺をラプラス変換すると，

　\begin{eqnarray} s Y(s) - y(0) + a \, Y(s) &=& 1 \end{eqnarray}

になります．初期値 y(0)=0を使うと，

　\begin{eqnarray} s Y(s) + a \, Y(s) &=& 1 \\ Y(s) &=& \frac{1}{s+a} \end{eqnarray}

となります．両辺を逆変換すると，

　\begin{eqnarray} y(t) &=& e^{-at}\, u(t) \end{eqnarray}

がえられます．ここで，u(t)は単位ステップ関数

　\begin{eqnarray} u(t)= \begin{cases} 0 & (t \le 0) \\ 1 & (t > 0)\end{cases} \end{eqnarray}

です．

　単位ステップ関数の別の定義として，

　\begin{eqnarray} u(t)= \begin{cases} 0 & (t < 0) \\ 1 & (t \ge 0)\end{cases} \end{eqnarray}

が使われることがあります．この定義だと， y(t) = e^{-at}\, u(t)とすると，y(0) = 1になって，初期条件 y(0) = 0を満たさなくなってしまします．明らかにダメです．

まとめ

　\begin{eqnarray} \frac{dy(t)}{dt} + a \, y(t) &=& \delta(t), \quad y(0) = 0 \end{eqnarray}

を例として解いてみました．どの解き方でも解は一致します．ここで考えたかったことは，t = 0の一瞬をどうとらえるかです．私としては，t = 0の瞬間にすべてが起きたと考えるのではなく，t=0の直後の無限に短い時間の間に応答が生じたと考えたいということです．

デルタ関数，単位インパルス関数は，どこに立っている？

数学 物理

ディラック (Dirac)のデルタ関数や，単位インパルス関数（単に，インパルス関数と呼ぶことも）は，分野によって立っている位置が違うので，一つの定義にこだわっていると，計算結果に納得できないことがあります．「デルタ関数と単位インパルス関数は同じモノだよ」と明確に書かれている教科書がありますが，分野によって定義が微妙に違い，その違いを理解することが重要というのが今回のお話です．
　
　以下では，デルタ関数，あるいは，単位インパルス関数を\delta(t)とします．\delta(t)の基本性質

　\begin{eqnarray}\int_{-\infty}^{\infty} \! \delta(t) \, dt &=& 1 \\ \int_{-\infty}^{\infty} \! f(t) \, \delta(t) \, dt &=& f(0) \\ \int_{-\infty}^{\infty} \! f(t) \, \delta(t-a) \, dt &=& f(a) \end{eqnarray}

を満たすことは，以下のどの場合も共通です．

私が最初に出会ったDiracのデルタ関数

　私は物理学科出身なので，デルタ関数を以下の図を使って理解していました．

デルタ関数の定義

　つまり，上の図のように，時刻t=0について軸対称な長方形を考えます．幅\varepsilon ，高さ\displaystyle \frac{1}{\varepsilon }とするので，面積は1です．この長方形について，\varepsilon \to 0の極限をとったものがデルタ関数\delta(t)です．

　この考えた方だと，

　\begin{eqnarray} \int_{0}^{\infty} \! \delta(t)\, dt &=& \frac{1}{2} \\ \int_{0}^{\infty} \! f(t)\, \delta(t)\, dt &=& \frac{1}{2}\, f(0) \end{eqnarray}

になりそうです．実際，いくつかの教科書では，このように書かれています．

　たまに，

　\begin{eqnarray} \delta(t)=\left\{\begin{array}{cc} 0 &(t \neq 0) \\ \infty &(t=0) \end{array}\right. \end{eqnarray}

などと，強調して書いてあることがあります．しかし，\delta(t)は，積分の中に登場して初めて意味をもつのであって，こう書いても何の意味もないし，何の役にも立たないと思います．このことは，どこかの本で読んだのか，どこかの講義で聴いたのか忘れましたが，私は経験からもそう感じるようになりました．\delta(t)は，線ではなく，幅と有限な面積をもつものです．

工学系の教科書で出会った単位インパルス関数

　私が，\delta(t)の定義に違和感を覚えたのは，線形システムのインパルス応答を計算したときです．つまり，\delta(t)のラプラス変換（片側ラプラス変換）の計算です．

　\delta(t)のラプラス変換では，

　\displaystyle \int_{0}^{\infty} \! \delta(t) \, e^{-st}\, dt

を計算します．上の図をイメージすると，t \ge 0では，デルタ関数が半分しか含まれないので，結果に\frac{1}{2}が付くよね？と考えてしまいます．

　しかし，一般的な結果は，

　\displaystyle \int_{0}^{\infty} \! \delta(t) \, e^{-st}\, dt = 1

です．

　何で？と，私はしばらく悩んだ経験があります．私は\frac{1}{2}が正しいと思い込み，さんざん世間に不平不満を言ったあげくに，\delta(t)は，上とは違う図 (t=0で軸対称ではない)で定義されているということを知りました．つまり，下の図です．

単位インパルスの定義

　この図では，t=0の右側 (t \ge 0)に長方形を置き，左端を固定して右から押しつぶしていいます．「デルタ関数と単位インパルス関数は同じモノ」とか，言わないでほしかったです．個人ごとの都合に合わせて，\delta(t)の定義は変わるよ，と最初から注意しておいてほしかったです．

　結局，この定義だと，

　\begin{eqnarray}\int_{0}^{\infty} \! \delta(t) \, dt &=& 1 \\ \int_{0}^{\infty} \! f(t) \, \delta(t) \, dt &=& f(0) \end{eqnarray}

ということです．つまり，物理よく使う\delta(t)とは違うものです．

単位インパルス関数は t > 0に立っててほしい (私の願望)

　上の単位インパルス関数\delta(t)の定義を使っても，線形システムのインパルス応答を考えるときに，私には気持ち悪さが残っています．この気持ち悪さを消すために，私としては，下の図の右のように\delta(t)を，t > 0に置いてほしいです (このように\delta(t)を定義している教科書も存在します)．

\delta(t)を，t = 0 (左)じゃなくて，t > 0 (右)に立てほしい．

　
　\delta(t)を，t > 0に置くというのは，下の図で言えば，(b)か(c)のパターンです．どの場合も，\varepsilon \to 0としたものが，\delta(t)です．

\delta(t)はどこに立っているのか？

　(b)か(c)のパターンで\delta(t)を定義すると，\delta(0)の値はどうなるの？と言われそうです．上の図の(b)，(c)では，\delta(0) = 0です．おいおい，\delta(0)=\inftyだろうと，口撃されそうです．上の図の(b)，(c)を使った定義では，\delta(t) = \inftyとなるtの値を明示できません．

　最初の方でも言いましたが，\delta(0)や\delta(t)の値なんて考えても，何の意味もないし，何の役にも立たないと思います．私は数学者ではないので詳しいことは知りませんが，\delta(t)は数学的には関数じゃないそうです (面積をもつ不思議な線)．関数っぽく書くから，\delta(0)の値を覚えて満足する人が出てくるのだと思います．

　\delta(t)が，t > 0に立っていて，面積をもっているのであれば，明確に

　\begin{eqnarray}\int_{0}^{\infty} \! \delta(t) \, dt &=& 1 \\ \int_{0}^{\infty} \! f(t) \, \delta(t) \, dt &=& f(0) \end{eqnarray}

です．単純バカの私としては，この形で話をしてほしいです．

まとめ

　\delta(t)は積分の中でこそ輝くし，\delta(t)の面積が，どこに立っているのかを意識して区別しないと直感的に理解できないというのが，私の考えです．最後に注意しておきますが，今回の説明はあくまで私の個人的感想ですので，数学的な正しさは保証できません．

【Rで最小2乗フィット】QR分解 (Householder変換)を使った最小2乗法

北川先生の本「Rによる 時系列モデリング入門」で，Householder変換を使った最小2乗フィットの話が登場したので，Rでやってみます．

　ここでは，以下の式を使ってデータを生成します．

　\displaystyle y_n = \sum_{i=1}^{2} a_i x_{ni} + \varepsilon_n

ここで，\displaystyle \varepsilon_nは，平均0，標準偏差sd.epsの正規乱数とします．

　以下のスクリプトでは，データの長さをN，パラメタa_1をa1，a_2をa2に設定します．

　\vec{y} = \begin{bmatrix} y_1 \\ y_{2} \\ \vdots\\ y_{N} \\ \end{bmatrix}, \quad Z = \begin{bmatrix} x_{11} & x_{12} \\ x_{21} & x_{22} \\ \vdots & \vdots \\ x_{N1} & x_{N2} \\ \end{bmatrix}, \quad \vec{a} = \begin{bmatrix} a_1 \\ a_2 \\ \end{bmatrix} , \quad \vec{\varepsilon} = \begin{bmatrix} \varepsilon_{1} \\ \varepsilon_{2} \\ \vdots \\ \varepsilon_{N} \\ \end{bmatrix}

と表せば，

　\vec{y} = Z \vec{a} + \vec{\varepsilon}

です．ベクトルはすべて縦書きの列ベクトルとします．

逆行列を使った最小2乗法

　残差としての\{\varepsilon_{n}\} の2乗和が最小になるように \vec{a } を計算すると，

　\vec{a } = \left({}^t Z Z \right)^{-1} {}^t Z \, \vec{y}

となります．{}^t Zは，Zの転置行列を表します．説明は，過去の記事を参考にしてください．
chaos-kiyono.hatenablog.com

　以下は，N = 10，a1 = 2.4，a2 = -3.0とした例です．Rのコマンドを使って逆行列を計算して，パラメタを推定しています．

# データの長さ
N <- 10
# 誤差の標準偏差
sd.eps <- 1
eps <- rnorm(N,0,sd.eps)
##############################
a1 <- 2.4
a2 <- -3.0
##############################
# 真のモデルは
# y[i] = a1 * x1[i] + a2 * x2[i] + eps[i]
##############################
x1 <- rnorm(N,0,1)
x2 <- rnorm(N,0,1)
y <- a1*x1+a2*x2+eps 
##############################
Y <- matrix(y[1:N],　ncol=1) # ncolで列数のみ指定
Z <- matrix(c(x1[1:N],x2[1:N]),　ncol=2)
##############################
# a.hat = (t(Z) Z)^{-1} t(Z) Y 
##############################
a.hat <- solve(t(Z) %*% Z) %*% (t(Z) %*% Y)
# 残差の２乗平均
sig2.hat <- sum((Y - Z %*% a.hat)^2)/nrow(Y)
# a1, a2の推定結果
a.hat[,1]
# 残差の2乗平均 sig2
sig2.hat

　このスクリプトを実行すると，

> # a1, a2の推定結果
> a.hat[,1]
[1]  2.425303 -3.634048
> # 残差の2乗平均 sig2
> sig2.hat
[1] 1.607804

のように，パラメタと残差が計算されます．

QR分解 (Householder変換)を使った最小2乗法

　次は，QR分解 (Householder変換)でパラメタを推定してみます．pracmaパッケージを使いました．最初に，このパッケージをインストールしてください．

　使ったRスクリプトは以下です．

#install.packages("pracma")
require(pracma)
# データの長さ
N <- 10
# 誤差の標準偏差
sd.eps <- 1
eps <- rnorm(N,0,sd.eps)
##############################
a1 <- 2.4
a2 <- -3.0
##############################
# 真のモデルは
# y[i] = a1 * x1[i] + a2 * x2[i] + eps[i]
##############################
x1 <- rnorm(N,0,1)
x2 <- rnorm(N,0,1)
y <- a1*x1+a2*x2+eps 
##############################
Y <- matrix(y[1:N],　ncol=1) # ncolで列数のみ指定
Z <- matrix(c(x1[1:N],x2[1:N]),　ncol=2)
###
# X = [Z|Y] と結合
X <- cbind(Z,Y)
##############################
# Householder変換
UX <- householder(X)$R
###
m <- ncol(UX)-1
a.hat <- numeric(m)
a.hat[m] <- UX[m,m+1]/UX[m,m]
for(i in 2:m){
 a.hat[m+1-i] <- (UX[m+1-i,m+1]-a.hat[(m+2-i):m] %*% UX[m+1-i,((m+2-i)):m])/UX[m+1-i,m+1-i]
}
# 残差の２乗平均
sig2.hat <- UX[3,3]^2/nrow(Y)
# a1, a2の推定結果
a.hat
# 残差の2乗平均 sig2
sig2.hat

　このスクリプトを実行すると，

> # a1, a2の推定結果
> a.hat
[1]  2.611575 -2.950162
> # 残差の2乗平均 sig2
> sig2.hat
[1] 0.6768868

のように，パラメタと残差が計算されます．

　pracmaパッケージを使わなくても，Rのコマンドqrを使えば，QR分解できます．中身のアルゴリズムについては，何を使っているのか分かりません．Rスクリプトは以下です．

# データの長さ
N <- 10
# 誤差の標準偏差
sd.eps <- 1
eps <- rnorm(N,0,sd.eps)
##############################
a1 <- 2.4
a2 <- -3.0
##############################
# 真のモデルは
# y[i] = a1 * x1[i] + a2 * x2[i] + eps[i]
##############################
x1 <- rnorm(N,0,1)
x2 <- rnorm(N,0,1)
y <- a1*x1+a2*x2+eps 
##############################
Y <- matrix(y[1:N],　ncol=1) # ncolで列数のみ指定
Z <- matrix(c(x1[1:N],x2[1:N]),　ncol=2)
###
# X = [Z|Y] と結合
X <- cbind(Z,Y)
##############################
# QR分解でRを求める
UX <- qr.R(qr(X))
###
m <- ncol(UX)-1
a.hat <- numeric(m)
a.hat[m] <- UX[m,m+1]/UX[m,m]
for(i in 2:m){
 a.hat[m+1-i] <- (UX[m+1-i,m+1]-a.hat[(m+2-i):m] %*% UX[m+1-i,((m+2-i)):m])/UX[m+1-i,m+1-i]
}
# 残差の２乗平均
sig2.hat <- UX[3,3]^2/nrow(Y)
# a1, a2の推定結果
a.hat
# 残差の2乗平均 sig2
sig2.hat

まとめ

　QR分解について，勉強してみてください．QR分解の方法には，グラム・シュミットの正規直交化法とか，Householder変換を使う方法とか，ギブンス回転を使う方法などがあるそうです．

【付録】3変数の場合

　逆行列を使う場合．

N <- 10
sig2 <- 1
eps <- rnorm(N,0,sqrt(sig2))
##############################
a1 <- 2.4
a2 <- -3.0
a3 <- -1.8
##############################
# 真のモデルは
# y[i] = a1 * x1[i] + a2 * x2[i] + a3 * x3[i] + eps[i]
##############################
x1 <- rnorm(N,0,1)
x2 <- rnorm(N,0,1)
x3 <- rnorm(N,0,1)
y <- a1*x1+a2*x2+a3*x3+eps 
##############################
Y <- matrix(y[1:N],　ncol=1) # ncolで列数のみ指定
Z <- matrix(c(x1[1:N],x2[1:N],x3[1:N]),　ncol=3)
##############################
# a.hat = (t(Z) Z)^{-1} t(Z) Y 
##############################
a.hat <- solve(t(Z) %*% Z) %*% (t(Z) %*% Y)
# 残差の２乗平均
sig2.hat <- sum((Y - Z %*% a.hat)^2)/nrow(Y)
# a1, a2の推定結果
a.hat[,1]
# 残差の2乗平均 sig2
sig2.hat

　Householder変換を使う場合．

#install.packages("pracma")
require(pracma)
###############
N <- 10
sig2 <- 1
eps <- rnorm(N,0,sqrt(sig2))
##############################
a1 <- 2.4
a2 <- -3.0
a3 <- -1.8
##############################
# 真のモデルは
# y[i] = a1 * x1[i] + a2 * x2[i] + a3 * x3[i] + eps[i]
##############################
x1 <- rnorm(N,0,1)
x2 <- rnorm(N,0,1)
x3 <- rnorm(N,0,1)
y <- a1*x1+a2*x2+a3*x3+eps 
##############################
Y <- matrix(y[1:N],　ncol=1) # ncolで列数のみ指定
Z <- matrix(c(x1[1:N],x2[1:N],x3[1:N]),　ncol=3)
##############################
# X = [Z|Y] と結合
X <- cbind(Z,Y)
##############################
# Householder変換
UX <- householder(X)$R
###
m <- ncol(UX)-1
a.hat <- numeric(m)
a.hat[m] <- UX[m,m+1]/UX[m,m]
for(i in 2:m){
 a.hat[m+1-i] <- (UX[m+1-i,m+1]-a.hat[(m+2-i):m] %*% UX[m+1-i,((m+2-i)):m])/UX[m+1-i,m+1-i]
}
# 残差の２乗平均
sig2.hat <- UX[3,3]^2/nrow(Y)
# a1, a2の推定結果
a.hat
# 残差の2乗平均 sig2
sig2.hat

ロジスティック写像のカオス

今回は，決定論的カオス (deterministic chaos)の説明用の動画を作りました．決定論的カオス (以下では，カオス)というのは，時間発展を決定論的に決めるルールがあるのに，未来が予測できなくなる現象のことです．

カオスとは

　カオスを示す決定論的システムの具体例として，ロジスティック写像

　x_{n+1} = C x_n (1-x_n)

を考えます．ここで， 0 < x_n < 1で，Cは制御変数です (0 < C \le 4)．Cの値を変えると， x_nの漸近的な振舞が変化します．

　C = 4の場合が，カオスの例になっています．カオスの主な特徴は，

  • 短期予測可能
  • 長期予測不可能
  • 非周期パターン

です．下の図では，

　下の図は，x_{n+1} = 4 x_n (1-x_n)で，初期値x_1 = 0.3 (青)の場合と x_1 = 0.300000000000001 (赤)の場合を比較した結果です．初期値の差は，10^{-15}です．

ロジスティック写像のカオス．初期値x_1 = 0.3 (青)と x_1 = 0.300000000000001 (赤)の比較．

　図の左は，クモの巣図法と呼ばれる時間発展をグラフを使って調べる方法のアニメです．図の右は， x_nの時間発展のグラフです．

　 n= 40くらいまでは，青と赤の2つの点がほとんど同じ振舞を示します (2点は重なっています)．これが，短期予測可能ということです．初期値の値が近ければ，短期的には似た振舞になります．

　しかし， n= 40をこえると，青と赤の2つの点が全く違った振舞を示すようになります．これが，長期予測不可能ということです．カオスでは，初期値鋭敏性といって，ほんのわずかな差が，指数関数的に大きくなるため，どんなに小さい差でも，いつかはシステムサイズくらいの差に成長してしまいます．

周期倍分岐

　ロジスティック写像で制御変数Cの値を変えると，カオス以外にも，いろいろと面白い現象が見られます．例えば，同じ値が繰り返すパターンです．Cを0から次第に増加させると，1つの値に収束したり (1周期)，2つの値を繰り返したり (2周期)，4つの値を繰り返したり (4周期)と，周期が倍々になる周期倍分岐というのが現れます．下の図は，いくつかのCでの振舞です．

C=2.8，初期値x_1 = 0.3の場合，1周期に漸近．
C=3.2，初期値x_1 = 0.3の場合，2周期に漸近．
C=3.54，初期値x_1 = 0.3の場合，4周期に漸近．

無相関と独立性の違い

　x_{n+1} = 4 x_n (1-x_n)の振舞は，無相関と独立性の違いを考える上で良い例になっています．つまり，この写像で生成される時系列\{x_n\}は，無相関な時系列になっています．でも，\{x_n\}は，短期予測可能なので，異なる時刻の値同士は独立ではありません．相関と独立の違いを理解しておいてください．

まとめ

　私は，博士 (理学)をとるまでは，カオス力学系の研究をしていました．その後，生理学とか，医学への応用を指向しました．カオスについての読み物としては，ジェイムズ・グリック著の「カオス―新しい科学をつくる」 (新潮文庫，1991)をお勧めします．ただし，カオスは，一時，ブームになり私もその影響を受けましたが，現在は，グリックの本で書かれているような勢いはないです．私は，博士過程を修了した後，運良く日本学術振興会の特別研究員（ＰＤ）に採択されました．そのころから，カオスとか，自分の狭い専門性に頼った研究をするのではなく，世の中の役に立つことを実現するために，新しいことに積極的に挑戦するようになりました．

【Rで非整数ブラウン運動】拡大アニメの作成

R言語 非整数ブラウン運動 1/fゆらぎ フラクタル

非整数ブラウン運動の説明用のアニメを作ってみました．

ブラウン運動 (H = 0.5)のサンプルパスの拡大
非整数ブラウン運動 (H = 0.8)のサンプルパス
非整数ブラウン運動 (H = 0.3)のサンプルパス

高解像度な，Youtube版は以下です．

www.youtube.com

www.youtube.com

www.youtube.com

Rスクリプト

　使ったＲスクリプトはこれです．

###Rでgifアニメーションの作成
# magickパッケージを使います
# install.packages("magick")
# ライブラリの読み込み
require("magick")
# サイズを大きくすると作成に時間がかかります．
img = image_graph(width = 480,height = 240)
################
# ハースト指数の指定
H <- 0.5
################
seed <- 9
set.seed(seed)
par(mar=c(2,1,3,1))
################
# 中点法でサンプルパスの生成
n.step <- 13
i <- 2
x <- c(0,1)
y <- c(0,0)
for(j in 1:n.step){
 h <- (1/2)^j
 x.tmp <- x[1:i] #
 y.tmp <- y[1:i]
 x[(1:(i-1))*2] <- (x.tmp[1:(i-1)]+x.tmp[2:i])/2
 y[(1:(i-1))*2] <- (y.tmp[1:(i-1)]+y.tmp[2:i])/2+h^(H)*rnorm(i-1)
 x[(1:i)*2-1] <- x.tmp
 y[(1:i)*2-1] <- y.tmp
 i <- i + (i-1)
}
N <- length(x)
#######################
r <- 1/4
x1.min <- 0
x1.max <- 1
y1.min <- min(y)
y1.max <- max(y)
h1 <- (y1.max-y1.min)/20
yp1.min <- y1.min-h1
yp1.max <- y1.max+h1
h1 <- y1.max - y1.min
plot(x,y,"l",xlim=c(x1.min,x1.max),ylim=c(yp1.min,yp1.max),col=2,bty="n",xaxt="n",yaxt="n",xlab="",ylab="",main=paste("H = ",H,sep=""))
rect(x1.min,yp1.min,x1.max,yp1.max,col=rgb(1,0,0, alpha=0.01))
#######################
for(k in 1:20){
i.s <- round(N*r)
s <- x[i.s]-x[1]
h.s <- h1*(s/(x1.max-x1.min))^H
i0 <- 1
h.tmp <- 0
for(i in 1:(N-i.s)){
 tmp <- max(y[(i+1):(i+i.s)])-min(y[(i+1):(i+i.s)])
 if(tmp < h.s){
  if(tmp > h.tmp){
    i0 <- i+1
    h.tmp <- tmp
  }
 }
}
###########################
x2.min <- x[i0]
x2.max <- x[i0+i.s-1]
y2.min <- min(y[(i0+1):(i0+i.s)])
y2.max <- max(y[(i0+1):(i0+i.s)])
h2 <- (y2.max-y2.min)/20
yp2.min <- y2.min-h2
yp2.max <- y2.max+h2
h2 <- y2.max - y2.min
rect(x2.min,yp2.min,x2.max,yp2.max,col=rgb(1,0,0, alpha=0.01))
###########################
n.zoom <- 10
r.zoom <- seq(1,n.zoom,length=n.zoom)
for(j in 1:n.zoom){
 x1 <- (r.zoom[j]*x2.min+(n.zoom-r.zoom[j])*x1.min)/n.zoom
 y1 <- (r.zoom[j]*yp2.min+(n.zoom-r.zoom[j])*yp1.min)/n.zoom
 x2 <- (r.zoom[j]*x2.max+(n.zoom-r.zoom[j])*x1.max)/n.zoom
 y2 <- (r.zoom[j]*yp2.max+(n.zoom-r.zoom[j])*yp1.max)/n.zoom
 plot(x,y,"l",xlim=c(x1,x2),ylim=c(y1,y2),col=2,bty="n",xaxt="n",yaxt="n",xlab="",ylab="",main=paste("H = ",H,sep=""))
 rect(x2.min,yp2.min,x2.max,yp2.max,col=rgb(1,0,0, alpha=0.01))
}
##########################
x1.min <- x2.min
x1.max <- x2.max
y1.min <- y2.min
y1.max <- y2.max
h1 <- (y1.max-y1.min)/20
yp1.min <- yp2.min
yp1.max <- yp2.max
h1 <- h2
###
for(ii in 1:2){
 x.tmp <- x[x >= x2.min & x <= x2.max]
 y.tmp <- y[x >= x2.min & x <= x2.max]
 h <- (x.tmp[2]-x.tmp[1])/2
 x <- c()
 y <- c()
 n.x <- length(x.tmp)
 x[(1:(n.x-1))*2] <- (x.tmp[1:(n.x-1)]+x.tmp[2:n.x])/2
 y[(1:(n.x-1))*2] <- (y.tmp[1:(n.x-1)]+y.tmp[2:n.x])/2+h^(H)*rnorm(n.x-1)
 x[(1:n.x)*2-1] <- x.tmp
 y[(1:n.x)*2-1] <- y.tmp
}
plot(x,y,"l",xlim=c(x1.min,x1.max),ylim=c(yp1.min,yp1.max),col=2,bty="n",xaxt="n",yaxt="n",xlab="",ylab="",main=paste("H = ",H,sep=""))
rect(x1.min,yp1.min,x1.max,yp1.max,col=rgb(1,0,0, alpha=0.01))
###
N <- length(x)
}
#描画領域を閉じる
dev.off()
#gifの速度を設定
animation = image_animate(img, fps=5, loop=0)
#フォルダに保存
image_write(animation, path = "fBm_zoom.gif", format = "gif")

非整数ブラウン運動について

　非整数ブラウン運動やフラクタルについては，過去の記事を参照してください．
chaos-kiyono.hatenablog.com
chaos-kiyono.hatenablog.com
chaos-kiyono.hatenablog.com
chaos-kiyono.hatenablog.com
chaos-kiyono.hatenablog.com
chaos-kiyono.hatenablog.com
chaos-kiyono.hatenablog.com

【長時間相関解析】DFAでホワイトノイズを分析するとどうなるか（つづき）：数値実験

R言語 DFA

Detrended Fluctuation Analysis (DFA)で，ホワイトノイズを解析したらどうなるか，数値実験で確かめてみます．理論的な結果は，前の記事を参照してください．
chaos-kiyono.hatenablog.com

数値時系列の生成

　Rでホワイトノイズを生成するのに，以下のスクリプトを使います．スクリプトを実行すれば，xに時系列ベクトル，sigに推定した標準偏差が入ります．

# 時系列の長さ
N <- 10000
# 正規乱数の生成
x <- rnorm(N)
# 標準偏差を推定
sig <- sd(x)

1次DFA

　計算は速くないですが，Rで1次のDetrended Fluctuation Analysis (DFA1)を実装してみます．使ったスクリプトは以下です．xにホワイトノイズのサンプル時系列が入っている必要があります．

#
# xを時系列ベクトルとする
#
N <- length(x)
###############################
#　平均0化
x <- x - mean(x)
# 積分
y <- cumsum(x)
###############################
# scalesの設定
scale.min <- 3
scale.max <- N/4
scale <- unique(round(10^seq(log10(scale.min),log10(scale.max),0.05)))
n.scale <- length(scale)
###############################
F2 <- c()
for(i in 1:n.scale){
 i.sub <- seq(1,N,scale[i])
 n.sub <- length(i.sub)-1
 F2[i] <- 0
 ii <- 1:scale[i]
 for(j in 1:n.sub){
   y.sub <- y[i.sub[j]:(i.sub[j+1]-1)]
   fit <- lm(y.sub~ii)
   a0 <- fit$coefficients[[1]]
   a1 <- fit$coefficients[[2]]
   F2[i] <- F2[i] + mean((a0+a1*ii-y.sub)^2)
 }
 F2[i] <- F2[i]/n.sub
}
###############################
# 結果のプロット
# 実用の場面では，スケーリング領域を各自設定して，傾きを推定してください
par(mfrow=c(1,1),pty="s")
fit <- lm(log10(F2)/2~log10(scale))
plot(log10(scale),log10(F2)/2,col=4,,xlab="log10 s",ylab="log10 F(s)",main="DFA1")
# 理論的なスケール依存性
curve(log10(sig*sqrt(((10^x)^2-4)/(15*(10^x)))),xlim=c(log10(scale.min),4),log="xy",col=2,lty=2,lwd=2,las=1,add=TRUE)
legend("bottomright",legend=c("Numerical estimation","Theoretical prediction"),pch=c(1,NA),lty=c(NA,2),col=c(4,2),lwd=c(1,2))

　下の図が実行例です．小さいスケールで，理論的結果と数値実験結果が良く一致しています．

Rスクリプトの実行例．

2次DFA

　次は，Rで2次のDetrended Fluctuation Analysis (DFA2)を実装してみます．使ったスクリプトは以下です．今回も，xにホワイトノイズのサンプル時系列が入っている必要があります．

#
# xを時系列ベクトルとする
#
N <- length(x)
###############################
#　平均0化
x <- x - mean(x)
# 積分
y <- cumsum(x)
###############################
# scalesの設定
scale.min <- 4
scale.max <- N/4
scale <- unique(round(10^seq(log10(scale.min),log10(scale.max),0.05)))
n.scale <- length(scale)
###############################
F2 <- c()
for(i in 1:n.scale){
 i.sub <- seq(1,N,scale[i])
 n.sub <- length(i.sub)-1
 F2[i] <- 0
 ii <- 1:scale[i]
 for(j in 1:n.sub){
   y.sub <- y[i.sub[j]:(i.sub[j+1]-1)]
   fit <- lm(y.sub~ii+I(ii^2))
   a0 <- fit$coefficients[[1]]
   a1 <- fit$coefficients[[2]]
   a2 <- fit$coefficients[[3]]
   F2[i] <- F2[i] + mean((a0+a1*ii+a2*(ii^2)-y.sub)^2)
 }
 F2[i] <- F2[i]/n.sub
}
###############################
# 結果のプロット
# 実用の場面では，スケーリング領域を各自設定して，傾きを推定してください
par(mfrow=c(1,1),pty="s")
fit <- lm(log10(F2)/2~log10(scale))
plot(log10(scale),log10(F2)/2,col=4,,xlab="log10 s",ylab="log10 F(s)",main="DFA2")
# 理論的なスケール依存性
curve(log10(sig*sqrt((3*(10^x)^2-27)/(70*(10^x)))),xlim=c(log10(scale.min),4),log="xy",col=2,lty=2,lwd=2,las=1,add=TRUE)
legend("bottomright",legend=c("Numerical estimation","Theoretical prediction"),pch=c(1,NA),lty=c(NA,2),col=c(4,2),lwd=c(1,2))

　下の図が実行例です．直線からのズレは，DFA2の方が，DFA1よりも大きく見えます．

Rスクリプトの実行例．

【長時間相関解析】DFAでホワイトノイズを分析するとどうなるか

R言語 DFA

Detrended Fluctuation Analysis (DFA)で，無相関な信号，つまり，ホワイトノイズを解析したらどうなるかのお話です．DFAに詳しい人の多くは，そんなのスケールをsとして，ゆらぎ関数 (fluctuation function)は，

　F(s) \propto s^{1/2}

だろというかもしれません．sが大きい，漸近的な振舞を考えれば，これは正しいです．でも，全スケールにわたって正しいわけでもないし，具体的な関数形は分かりません．ということで，今回は，F(s)を解析に求めると，という話です．

　先に結論をまとめておくと，解析する観測時系列\{x_i\}の標本分散を\hat{\sigma}^2とすれば，fluctuation function F(s)は，以下のようになります．

  • 1次DFAでは，

　 \begin{eqnarray} F(s) &=& \sqrt{\frac{s^2-4}{15 s}} \hat{\sigma} \end{eqnarray}

  • 2次DFAでは，

　 \begin{eqnarray} F(s) &=& \sqrt{\frac{3 s^2-27}{70 s}} \hat{\sigma} \end{eqnarray}

  • 3次DFAでは，

　 \begin{eqnarray} F(s) &=& \sqrt{\frac{2 s^2- 32}{63 s}} \hat{\sigma} \end{eqnarray}

  • 4次DFAでは，

　 \begin{eqnarray} F(s) &=& \sqrt{\frac{5 s^2-125}{198 s}}\hat{\sigma} \end{eqnarray}

  • 5次DFAでは，

　 \begin{eqnarray} F(s) &=& \sqrt{\frac{3 s^2 - 108}{143 s}} \hat{\sigma} \end{eqnarray}

　上のどの場合でも，sが大きくなると，

　F(s) \propto s^{1/2}

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

　どうやって導出するかは，MarcとHolgerの論文
link.springer.com
に説明してあります．この論文中の式(12)とか、式(17)のW(s)の部分にタイプミスがありますが，計算間違いやタイプミスは誰にでもあるので，読む側が気をつければいいだけです．

解析的結果の導出

　詳しく説明するのは面倒なので，ポイントだけ言います．要は，

  • 時系列の積分は線形演算 (畳み込みで表現可能)
  • 最小二乗法でえられる回帰曲線の計算も線形演算 (畳み込みで表現可能)
  • 定義通りF^2(s)を計算すれば，自己共分散 (自己相関)との関係がでてくる
  • 自己共分散 (自己相関)を，単位インパルス関数 (ラグ0以外は全部ゼロ)として計算

ということです．つまり，面倒くさがらずに，真面目に計算するだけです．真面目に計算といっても，解いただけでは記憶に残らないので，そこにある構造や解釈のイメージをつかむように意識してください．

　n次DFAの場合について，計算の流れだけ具体的に書いておきます．

　まずは，以下の行列B(n)の逆行列B^{-1}(n)を計算します．これは，最小二乗法フィットの畳み込み係数を計算しています．

　\begin{eqnarray} B(n) &=& \left( \begin{array}{cccc} \sum _{j=1}^s 1 & \sum _{j=1}^s j & \cdots & \sum _{j=1}^s j^n \\ \sum _{j=1}^s j & \sum _{j=1}^s j^2 & \cdots & \sum _{j=1}^s j^{n+1} \\ \vdots & \vdots & \ddots & \vdots \\ \sum _{j=1}^s j^{n} & \sum _{j=1}^s j^{n+1} & \cdots & \sum _{j=1}^s j^{2 n} \\ \end{array} \right) \end{eqnarray}

　逆行列B^{-1}(n)の計算を手ですると大変ですが，数式処理ソフトのMathematicaとか，Maximaとかを使えば，屁くらいの手間です．以下では，B^{-1}(n)の要素を\{b^{-1}_{ij}\}とします．

　ここで，以下の関数を定義します．sは部分区間の長さ (スケール)です．

　 \begin{eqnarray}　P(t, k) &=& \sum_{i=1}^{n+1} t^{i-1} \sum_{j=1}^{n+1} b^{-1}_{ij} \sum_{r=k}^{s} r^{j-1} \\ &=& \left( \begin{array}{cccc} t^0 & t^1 & \cdots & t^{n} \\ \end{array} \right) \left( \begin{array}{cccc} b^{-1}_{11} & b^{-1}_{12} & \cdots & b^{-1}_{1n} \\ b^{-1}_{21} & b^{-1}_{22} & \cdots & b^{-1}_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ b^{-1}_{n1} & b^{-1}_{n2} & \cdots & b^{-1}_{nn} \\ \end{array} \right) \left( \begin{array}{c} \sum_{r=k}^{s} r^{0} \\ \sum_{r=k}^{s} r^{1} \\ \vdots \\ \sum_{r=k}^{s} r^{n} \\ \end{array} \right) \end{eqnarray}

　ということで，標本分散\hat{\sigma}^2のホワイトノイズのサンプル時系列のfluctuation functionは，

　\begin{eqnarray} F^2 (s) = \frac{1}{s} \sum_{p=1}^{s} \left\{ \hat{\sigma}^2 \sum_{k=p}^{s} \left(1 - 2 P(k,p) \right) + \sum_{k=1}^{s} P(k,p)^2 \right\} \end{eqnarray}

です．こんな知識，日本で知りたい人は5人もいないと思います．1人でも役に立てばうれしいので書いておきました．

Fluctuation functionのグラフ

　何でF(s)の全体像を考えるのかといえば，ときどき，勘違いで短時間のスケーリング指数や相関特性を議論しちゃう人がいるからです．

　上のF(s)のグラフは下の図のようになります．

ホワイトノイズをDFAで解析したときのゆらぎ関数 (fluctuation function)

　上の図では，sが小さいスケールで直線的になってません．この非直線的な構造は，時系列の性質ではなく，DFAの方法の都合で生じているのです．ですので，短いスケールの傾きから，相関の性質は議論できません．