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

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

【長時間相関解析】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でホワイトノイズを分析するとどうなるか

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の方法の都合で生じているのです.ですので,短いスケールの傾きから,相関の性質は議論できません.

【確率変数】変数変換したら確率密度関数はどうなるか

確率変数を変数変換する話です.

変数変換した後の確率密度関数

 ここでは,確率変数Xが,確率密度関数 f_{X}(x) に従うとします.このとき,単調関数 (単調増加か,単調減少する関数) g(x) を使って,確率変数 Y

 Y = g(X)

で与えられるとします.Yが従う確率密度関数 f_{Y}(y) を求めたいな,というのがここでの問題です.

 変数変換しても,確率は保存しないといけません.つまり,確率密度関数の対応する部分の面積が一致して,

  \begin{eqnarray} P\left(a < X < b \right) &=& P\left( g(a) < Y < g(b) \right) \\
\int_{x = a}^{x=b}\!f_{X}(x)\, dx &=& \int_{y=g(a)}^{y=g(b)}\!f_{Y}(y)\, dy \end{eqnarray}

となる必要があります.

変数変換したときの確率の保存

 左辺を,y = g(x)と置換した積分を考えると,

 \begin{eqnarray} \displaystyle \int_{x = a}^{x=b}\!f_{X}(x)\, dx = \int_{y=g(a)}^{y=g(b)}\!f_{X}\left(g^{-1}(y) \right) \frac{d \left(g^{-1}(y)\right)}{dy} \, dy 
 \end{eqnarray}

とできそうです.\frac{d \left(g^{-1}(y) \right)}{dy}は,対応する微小面積の倍率を表しています.確率は常に0以上で,マイナスになってはダメなので,プラスマイナスは無視することにすれば,

 \begin{eqnarray} \displaystyle \int_{x = a}^{x=b}\!f_{X}(x)\, dx = \int_{y=g(a)}^{y=g(b)}\!f_{X}\left(g(y) \right) \left| \frac{d \left(g^{-1}(y)\right)}{dy} \right| \, dy 
 \end{eqnarray}

となります (一般の多重積分も,向きの違いのプラスマイナスを無視するので,ヤコビアンには絶対値がついてます).

 ということで,確率 (面積)の保存を表す

  \begin{eqnarray} \int_{x = a}^{x=b}\!f_{X}(x)\, dx &=& \int_{y=g(a)}^{y=g(b)}\!f_{Y}(y)\, dy \end{eqnarray}

と対応する部分を比べると,

 \begin{eqnarray} f_{Y}(y) &=&  f_{X}(g^{-1}(y))\, \left| \frac{d \left(g^{-1}(y) \right)}{dy} \right|  \end{eqnarray}

となることがわかります.

 あらためて書くと,

確率変数X確率密度関数 f_{X}(x) に従うとき,単調関数を使ってY = g(X)と変数変換すると,Y確率密度関数は,

 \begin{eqnarray} f_{Y}(y) &=&  f_{X}(g^{-1}(y))\, \left| \frac{d \left(g^{-1}(y) \right)}{dy} \right| 
 \end{eqnarray}

ということです.

線形変換の例

 Xが標準正規分布 (平均0,分散1の正規分布)に従うとします.このとき,\sigma\muを定数として,

 \begin{eqnarray} f_{X}(x) &=&  \frac{1}{\sqrt{2 \pi}} e^{-\frac{x^2}{2}}  \end{eqnarray}

です.

 \sigma\muを定数として,確率変数 Y

 Y = \sigma X+\mu

で与えられる場合を考えます.

 \displaystyle X = \frac{Y-\mu}{\sigma} なので,\displaystyle g^{-1}(y) = \frac{y-\mu}{\sigma} です.

 したがって,確率変数 Y が従う確率密度関数は,

 \begin{eqnarray} f_{Y}(y) &=&  f_{X}(g^{-1}(y))\, \left| \frac{d \left(g^{-1}(y) \right)}{dy} \right| \\
&=&  f_{X}(g^{-1}(y))\, \frac{1}{\sigma} \\
&=&  \frac{1}{\sqrt{2 \pi}} e^{-\frac{\left( \frac{y-\mu}{\sigma} \right)^2}{2}}  \, \frac{1}{\sigma} \\
&=&  \frac{1}{\sqrt{2 \pi} \sigma} e^{-\frac{\left(y-\mu \right)^{\,2}}{2 \sigma^2}}
 \end{eqnarray}

となり,正規分布の式になります.

指数変換の例

 X正規分布に従うとします.このとき,\sigma\muを定数として,

 \begin{eqnarray} f_{X}(x) &=&  \frac{1}{\sqrt{2 \pi} \sigma} e^{-\frac{\left(y-\mu \right)^{\,2}}{2 \sigma^2}}  \end{eqnarray}

です.

 確率変数 Y

 Y = \exp X

で与えられる場合を考えます.

 \displaystyle X = \log Y なので,\displaystyle g^{-1}(y) = \log y です.

 したがって,確率変数 Y が従う確率密度関数は,

 \begin{eqnarray} f_{Y}(y) &=&  f_{X}(g^{-1}(y))\, \left| \frac{d \left(g^{-1}(y) \right)}{dy} \right| \\
&=&  f_{X}(g^{-1}(y))\, \frac{1}{y} \\
&=&  \frac{1}{\sqrt{2 \pi} \sigma y} e^{-\frac{\left(\log y-\mu \right)^{\,2}}{2 \sigma^2}}
 \end{eqnarray}

となり,対数正規分布の式になります.

【データ解析】相関係数を式で考えてみる

相関係数の話の続きです.
chaos-kiyono.hatenablog.com

 今回は確率変数を使って考えてみます.

Y=aX+bWの係数abと,XYの散布図の関係.

相関って線形な関係

 平均0,分散\sigma^2の確率変数Xと,これと独立な平均0,分散\sigma^2の確率変数Wを仮定します.XWの分布については,分散が有限ということだけで形は指定しません.

 Xと線形な依存性がある確率変数Yを,定数abを使って次の形で表します.

 \begin{eqnarray} Y = a X + b W\end{eqnarray}

この場合,a \neq 0のとき,YにはXと線形な関係があり,Wの係数の|b|が大きいほどその関係性が薄くなります.つまり,XYの相関は,aの絶対値が,bの絶対値よりも大きいほど強いと言えそうです.今回は実際に相関係数を計算してこのことを足しためてみます.

以下の計算で使う基本事項として,Yの期待値と分散を計算しておきます.Yの期待値は,

 \begin{eqnarray} \mathbb{E} [Y] &=& \mathbb{E} [a X + b W ] \\
 &=& a \mathbb{E} [X] + b \mathbb{E} [W ] \\ 
&=& a \cdot 0 + b  \cdot 0 \\ 
&=& 0 \end{eqnarray}

Yの分散は,

 \begin{eqnarray} {\rm Var} [Y] &=& {\rm Var} [a X + b W ] \\
&=& {\rm Var} [a X]  + {\rm Var} [b W ] \\
&=& a^2 {\rm Var} [X]  + b^2 {\rm Var} [ W ] \\
 &=& a^2 \sigma^2 + b^2 \sigma^2 \\ 
&=& \left( a^2 + b^2 \right) \sigma^2 \end{eqnarray}

です.

 上の例では,a=0ならXYは独立 (WXとは独立なので)です.また,b=0ならXYは単に比例しています.

 相関係数rの定義は,

 r = \frac{\displaystyle \mathbb{E} \big[ \left(X- \mathbb{E} \left[ X \right] \right)\left(Y- \mathbb{E} \left[ Y \right] \right)\big]}{\displaystyle \sqrt{\mathbb{E}  \big[ \left(X- \mathbb{E} \left[ X \right] \right)^2\big]} \sqrt{\mathbb{E}  \big[ \left(Y- \mathbb{E} \left[ Y \right] \right)^2\big]}}

です.

 上で定義した,XY相関係数を計算すると,

 \begin{eqnarray} r &=& \frac{ \mathbb{E} \big[ X \left(a X + bW \right) \big]}{\sqrt{\mathbb{E}  \big[ X^2 \big]} \sqrt{\mathbb{E}  \big[ \left(a X + b W \right)^2\big]}} \\
&=& \frac{ \mathbb{E} \big[ a X^2 + b X W \big]}{\sqrt{{\rm Var} \big[ X \big]} \sqrt{{\rm Var} \big[ a X + b W \big]}} \\
&=& \frac{ a\, \sigma^2}{\sigma \sqrt{\left( a^2 + b^2 \right) \sigma^2}} \\
&=& \frac{a}{ \sqrt{a^2+b^2}} 
\end{eqnarray}

になります.最後の式は,直角三角形をイメージしてもらって,直角をはさむ2辺の長さをabとしたときに,斜辺の長さ\sqrt{a^2+b^2}aの比になっています.

直角三角形の辺の長さ.

 上の結果を使うと,a=0のとき,つまり,YXに相関がないケースでは,r=0になります.

 また,b=0のとき,つまり,YXの定数倍のケースでは,r=1か,r=-1になります.

 そして,相関係数の範囲は,

 a > 0 のとき,0 < r \le 1

 a < 0 のとき,-1 \le r < 0

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

 ということで,相関係数がとる値の範囲は,-1 \le r \le 1です.

 また,「相関係数の推定量

 \hat{r} = \frac{\displaystyle \frac{1}{N}\sum_{i=1}^N \left(x_i-\bar{x} \right)\left(y_i-\bar{y} \right)}{\displaystyle \sqrt{\frac{1}{N}\sum_{i=1}^N \left(x_i-\bar{x} \right)^2} \sqrt{\frac{1}{N}\sum_{i=1}^N \left(y_i-\bar{y} \right)^2}}

についても,-1 \le \hat{r}  \le 1になります.上の計算はこの証明にはなっていないので,この範囲については自分で考えるか,インターネットで検索して調べてみてください.当然,いろんな教科書に載っています.

相関する確率変数の生成

 上の計算を簡単にすることを考えます.簡単化のため,互いに独立な確率変数XWは,平均0,分散1とします.そして,Yの分散も1にします.このとき,

  \begin{eqnarray} \mathbb{E} [Y^2] &=& \mathbb{E} [\left(a X + b W \right)^2 ] \\
 &=& a^2 + b^2 = 1 \end{eqnarray}

が成り立つ必要があるので,b^2=1-a^2となります.

 以上の条件でもう一度,相関係数rを計算してみると,

 \begin{eqnarray} r &=& \frac{ a}{\sqrt{a^2 +b^2}} \\
&=& \frac{a}{\sqrt{a^2 + \left( 1 - a^2 \right)}} \\
&=& a
\end{eqnarray}

になります.結局, a = rになりました.

 ということで,相関係数rの確率変数XYを作りたいとき,まず,平均0,分散1の互いに独立な確率変数XWを用意して,

 Y = r X + \sqrt{1-r^2} W

とすれば良いと言うことです.

Rで互いに相関する正規乱数の生成

 真の相関係数rになる正規乱数をRで生成してみます.使ったスクリプトはこれです.

# データの長さ
N <- 1000
# 真の相関係数
r <- 0.5
# 2変量乱数の生成
x <- rnorm(N) # 正規乱数
w <- rnorm(N) # 正規乱数
y <- r*x+sqrt(1-r^2)*w
# 相関係数の推定 
r.est <- cor(x,y)
###
# 散布図を描画
y.max <- max(abs(c(x,y)))
par(pty="s",cex.main=1.5)
plot(x,y,xlim=c(-y.max*1.3,y.max*1.3),ylim=c(-y.max*1.3,y.max*1.3),pch=16,col="blue",main=sprintf("r (est) = %.2f",r.est))
if(r >= 0){
 text(y.max*1.2,-y.max*1.2,sprintf("r (true) = %.2f",r),pos=2,cex=1.4,col=2)
 text(-y.max*1.2,y.max*1.2,sprintf("N = %d",N),pos=4,cex=1.4,col=1)
}else{
 text(-y.max*1.2,-y.max*1.2,sprintf("r (true) = %.2f",r),pos=4,cex=1.4,col=2)
 text(y.max*1.2,y.max*1.2,sprintf("N = %d",N),pos=2,cex=1.4,col=1)
}
Rスクリプトの実行結果

【データ解析】相関係数の直感的理解

相関係数の話の続きです.

chaos-kiyono.hatenablog.com

 今回は,相関係数が正になったり,負になったりすることを直感的に理解することが目的です.

相関係数と散布図の関係.100点の場合.

相関係数の推定量

 対応する2変量の実現値\{(x_i, y_i)\}N点あったとき,相関係数rの推定量は,

 r = \frac{\displaystyle \frac{1}{N}\sum_{i=1}^N \left(x_i-\bar{x} \right)\left(y_i-\bar{y} \right)}{\displaystyle \sqrt{\frac{1}{N}\sum_{i=1}^N \left(x_i-\bar{x} \right)^2} \sqrt{\frac{1}{N}\sum_{i=1}^N \left(y_i-\bar{y} \right)^2}}

が一般的です.ここで,\bar{x}\bar{y}は,それぞれ,\{x_i\}\{y_i\}の標本平均です.

散布図の構造と相関係数

 上の式で計算した相関係数が,1に近かったり,0に近かったり,-1に近かったりすることと,散布図の関係はどうなってるのでしょうか.

 標本平均\bar{x}\bar{y}を引くのは,値のばらつきを0まわりにするためなので,以下では簡単化のため,\{x_i\}\{y_i\}の標本平均が0の場合を考えます.

 標本平均が0のとき,相関係数の推定量

 r = \frac{\displaystyle \frac{1}{N}\sum_{i=1}^N x_i \, y_i}{\displaystyle \sqrt{\frac{1}{N}\sum_{i=1}^N x_i^2} \sqrt{\frac{1}{N}\sum_{i=1}^N y_i^2}}

です.この分子

 \displaystyle \frac{1}{N}\sum_{i=1}^N x_i \, y_i

の値を,散布図の象限毎に分けて考えると,相関係数の正になったり,負になったりすることの直感的な意味が見えてきます.

 下の図で,第1象限と第3象限にある点は,符号が同じなので,

 \displaystyle x_i \, y_i > 0

です.それに対し,第2象限と第4象限にある点は,符号が逆なので,

 \displaystyle x_i \, y_i < 0

です.

 下の図では,一番右の図が,象限毎の\displaystyle x_i \, y_i の和をNで割った値を示しています.

相関係数の符号の意味.

 散布図が,右上がりだと,\displaystyle x_i \, y_i > 0の割合が多いので,相関係数は正になります.特に右上がり45°の直線に近いほど,\displaystyle x_i \, y_iの値が大きくなります.

 逆に散布図が,右下がりだと,\displaystyle x_i \, y_i < 0の割合が多いので,相関係数は負になります.

 相関係数が0に近いときは,\displaystyle x_i \, y_i > 0になる領域と,\displaystyle x_i \, y_i < 0になる領域で,ほぼ対称に点が分布しているということです.

今日のまとめ

 今日,Breaking Down 6を見ました.私は昔,番組を間違えてBreaking Down 3のペーパービューを買ってしまったことがあります.そのときは,せっかくなので何試合か見ましたが,試合の面白さがまったく分かりませんでした.その後,Breaking Downへの出場のためのオーディションをYoutubeで見るようになり,Breaking Downに出場して有名になりたいという人たちの振舞に,面白みを感じるようになりました.殴り合いの喧嘩って良いことではありませんが,外から見ているとドキドキする感じはします.

 私は昔,歌舞伎町の飲み屋で見ず知らずの方々と喧嘩になり,顔を激しく殴られました.今でも顔の左に痺れが残っています.あのとき脳にダメージを受けていなければ,もっと科学者として活躍できたかもしれませんが,今では命があっただけ良かったと思います.若いときのことを思い出すと,なんであんなに愚かだったのかと,恥ずかしくなります.

【時系列解析】相関係数

今回は,相関係数のお話です.相関というのは2つの変量間の直線的な関係性の程度を表しています.直線性とは,2つの変量の値を縦軸と横軸にとった散布図に現れる構造です.この構造のイメージを持ってほしいので,いくつか図を作りました.今日は眠いので,説明は省いて,図メインで載せていきます.

2つの変量x1とx2の間の相関係数rを変化させたとき

相関係数と直線性

 相関係数の定義は,ググってください.相関係数rは,-1から1の値をとります.データ数が1000点ぐらいあるときは,相関係数と散布図の関係は以下のようになります.

相関係数と散布図.1000点の場合.右の図は点の密度で色を変えたもの.

 データ数が少なく20点ぐらいのときは,相関係数と散布図の関係は以下のようになります.

相関係数と散布図.20点の場合.右の図は点の密度で色を変えたもの.

データ点数とp値

 人によっては,検定のp値は小さいほど良いとか,値が0.05以下でないと意味がないと思っているかもしれません.私は「p値」信者ではないので,そんなこと思いませんが.

 相関係数を議論するときは,相関の強さと,p値の小ささは分けて解釈する必要があります.私が論文を書くときは,「相関係数が0.4以上のときに,相関があると判断する」とか,相関係数に基準を設けます.基準は,0.4でも,0.3でもいいと思います.

データ点数とp値の関係.相関係数が0.4の場合.

 以下の例では,相関係数を0.01に設定しました.つまり,ほぼ相関はありません.それでも,データ数を増やしていけば,ときどき,pは0.05よりも小さくなります.これで,「有意な相関がある」って言える?ってことです.

データ点数とp値の関係.相関係数が0.01の場合.

【Rで時系列解析】計測時刻のずれを推定して修正

複数のデバイスで同時計測をしたとき,計測時刻のズレが発生することがあります.今回は,そのズレを相互相関を使って修正する方法について説明します.

 下の図では,左上のグラフのように,赤のx_1と青のx_2の時系列がズレています.本当は赤と青の線がほぼ重なるはずなのに,パソコンや計測機器の時刻設定の違いでズレてしまった.ということが,結構あります.

2本の時系列の時間のズレを推定

 この時間のズレは,相互相関係数の絶対値が最大になる時間差を推定することで評価できます.上の例では,右上の相互相関のグラフが最大になる点 (lag=38)が,推定された時間差になります.相関係数というのは,2つの変量の直線的な関係性を評価していますので,右下の散布図のように,相互相関が最大になる点 (lag=38)で,2つの変数の値が最も45度の直線的に集まってきます.

 時間差が推定できれば,x_2の時間を推定されたラグだけずらせば,上の左下のグラフのように,赤と青の線がほぼ重なります.

数値実験のRスクリプト

 数値的に時系列を生成して,時間差を補正するスクリプトの例が以下です.

# 真のラグを設定
lag.true <- 38
# サンプル時系列を生成
N <- 1000
time <- 1:(N+lag.true)
x.base <- dnorm(time,N/2,N/5)*sin(2*pi*time/(N/3.5))^2
#
x1 <- x.base[1:N]+rnorm(N)*0.0001
x2 <- x.base[(1+lag.true):(N+lag.true)]+rnorm(N)*0.0001
time <- 1:N
######################################
par(mfrow=c(2,2),cex=1,cex.lab=1.6)
par(pty="m",mgp=c(2.2,0.5,0))
plot(time,x1,"l",col=2,xlab="i",ylab="",yaxt="n",ylim=c(-0.0005,0.003),lwd=2)
lines(time,x2,col=4,lwd=2)
legend("topright",legend=c(expression(paste(x[1],"[i]")),expression(paste(x[2],"[i]"))),col=c(2,4),lwd=c(2,2),cex=1.2)
#################
# 相互相関を使って時間を推定
xcor <- ccf(x1,x2,lag.max=round(N/4),main="",ylab="Cross-Correlation",col="skyblue")
# 相互相関係数の絶対値が最大になるラグを検出
lag.max <- xcor$lag[which.max(abs(xcor$acf))]
# 推定されたラグを破線としてプロット
abline(v=lag.max,col=6,lty=2,lwd=2)
###
plot(time,x1,"l",col=2,xlab="i",ylab="",main=paste("lag =",lag.max),yaxt="n",ylim=c(-0.0005,0.003),lwd=2)
# x2の時刻をlag.maxだけずらす
lines(time+lag.max,x2,col=4,lwd=2)
legend("topright",legend=c(expression(paste(x[1],"[i]")),expression(paste(x[2],"[i-lag]"))),col=c(2,4),lwd=c(2,2),cex=1.2)
########
# 散布図
par(pty="s",mgp=c(0.9,0,0))
if(lag.max >=0){
 plot(x1[(1+lag.max):N],x2[1:(N-lag.max)],pch=16,xaxt="n",yaxt="n",xlab=expression(paste(x[1],"[i]")),ylab=expression(paste(x[2],"[i-lag]")),main=paste("lag =",lag.max),col=6)
}else{
 plot(x1[1:(N+lag.max)],x2[(1-lag.max):N],pch=16,xaxt="n",yaxt="n",xlab=expression(paste(x[1],"[i]")),ylab=expression(paste(x[2],"[i-lag]")),main=paste("lag =",lag.max),col=6)
}

時刻を補正するRスクリプト

 時刻データを補正するサンプルスクリプトが以下です.このスクリプトではt1t2に時刻ベクトル,x1x2に数値データを代入しないと動きません.

# サンプリング間隔
t.samp <- 1
#####################
### 基準時系列
# 時刻
t1 <- (POSIXctの時刻)
# データ
x1 <- (時系列のベクトル)
#####################
### 時刻調整を行う時系列
# 時刻
t2 <- (POSIXctの時刻)
# データ
x2 <- (時系列のベクトル)
#####################
t.start <- max(min(t1[is.finite(x1)]),min(t2[is.finite(x2)]))
t.end <- min(max(t1[is.finite(x1)]),max(t2[is.finite(x2)]))
t.resamp <- seq(t.start,t.end,t.samp) 
###
x1.tmp <- approx(t1,x1,xout=t.resamp,method="linear")$y
x2.tmp <- approx(t2,x2,xout=t.resamp,method="linear")$y
###
par(mfrow=c(1,3),cex=1,cex.lab=1.6)
plot(t.resamp,x1.tmp,"l",col=2,xlab="Original Time",ylab="")
legend("topleft",legend=c("Reference Time Series","Target Time Series"),col=c(2,4),lwd=c(1,1),cex=0.8)
lines(t.resamp,x2.tmp,col=4)
xcor <- ccf(x1.tmp,x2.tmp,lag.max=round(length(t.resamp)/4),main="",ylab="Cross-Correlation",col="skyblue")
t.lag <- xcor$lag[which.max(abs(xcor$acf))]
abline(v=t.lag,col="blue",lwd=2,lty=2)
###
plot(t.resamp,x1.tmp,"l",col=2,xlab="Corrected Time",ylab="",main=paste("Time Lag =",Lag.sec,"sec"))
Lag.sec <- t.lag*t.samp
lines(t.resamp+Lag.sec,x2.tmp,col=4)
legend("topleft",legend=c("Reference Time Series","Target Time Series"),col=c(2,4),lwd=c(1,1),cex=0.8)

 このスクリプトを適用した例が下の図です.

スクリプトの適用例

まとめ

 私はウクライナの研究者との共同研究を,4年くらい前から続けています.ロシア侵攻が発生した後も,ほぼ毎週,オンラインで打ち合わせをしています.しかし,2週間前から,電力不足などの問題で,ウクライナからミーティングに接続できなくなりました.

 私は東日本大震災のときに,福島県に住んでおり,震災後はしばらく近所の小学校に避難していました.そのとき,燃料不足で暖房が使えず,とても寒かった経験があります.寒すぎて私は一睡もできませんでした.そして,そんな寒さから逃れるために家族を,福島県から避難させました.原発事故よりも,寒さの方がきつかったです.

 ウクライナでは多くの人が寒さに苦しんでいると思います.状況が少しでも改善することを祈っています.