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

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

【時系列解析】なぜ,自己共分散(自己相関)が急激に減衰しないとパワースペクトルは定義できないのか?

北川 源四郎先生の本「Rによる 時系列モデリング入門」(33 ページ)では,

 「ラグkが大きくなるとき自己共分散関数C_kが急激に減衰し」

という条件が,パワースペクトル

 \displaystyle p(f) = \sum_{k=-\infty}^{\infty} C_k e^{- 2 \pi i k f}

の定義を与える部分に書いてあります (以下では,\displaystyle p(f)\displaystyle S(f)と書きます).何で? この条件「自己共分散関数C_kが急激に減衰」の部分を読み飛ばしてはいけないというのが今回のお話です.

 Wikipediaのウィーナー=ヒンチンの定理 (Wiener–Khinchin theorem)の説明
ja.wikipedia.org

にはそんなこと書いてないけど,と思うかもしれません.

1次自己回帰過程 (a=1\sigma^2=1)の自己共分散とパワースペクトル

自己共分散とパワースペクトルの関係

 時系列\{y[i]\}の自己共分散関数を,

 \begin{eqnarray}C(\, j\, ) = E \left[ y [ i ]\, y[i -  j ]\right]
\end{eqnarray}

とします.ここでは,Eで期待値を表しています.このとき,次のような有限の長さNの時系列のパワースペクトルの推定量の期待値\bar{S}(f_k)を計算してみます.ここでは,時系列のフーリエ変換を使ってパワースペクトルを計算すると (パワースペクトルの定義は複数あります (関連記事1)),

 \begin{eqnarray}
\bar{S}(f_k) &=& {\rm E} \left[ \frac{\left|Y_k \right|^2}{N} \right] \nonumber \\
&=& {\rm E} \left[ \frac{1}{N} \left(\sum_{n=0}^{N-1} y[n] \, e^{- i 2 \pi f_k n} \right)\left(\sum_{n=0}^{N-1} y[n] \, e^{+ i 2 \pi f_k n} \right) \right] \nonumber \\
&=& \frac{1}{N} \sum_{j=-(N-1)}^{N-1} \left(N-|\, j \,|\right) C(\, j \,)\, e^{- i 2 \pi f_k j}  \nonumber \\
&=& \sum_{j=-(N-1)}^{N-1} \left(1- \frac{|\, j \,|}{N}\right) C(\, j \,)\, e^{- i 2 \pi f_k j}
\end{eqnarray}

となります (ちょっと自信がありませんので自分で確かめてください).途中で,

 \displaystyle \frac{1}{N}\sum_{k=0}^{N-1} e^{i 2 \pi n k/N } = \left\{
\begin{array}{ll}
\displaystyle 1 & \displaystyle  (n = 0) \\
\displaystyle 0 & \displaystyle  (n \neq 0)
\end{array}\right.

の公式を使いました.これは,離散の単位インパルス (ディラックじゃありません)

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

フーリエ変換

 \displaystyle \sum_{n=0}^{N-1}  \delta [n]  e^{- i 2 \pi n k/N } = 1

になることの逆変換です.フーリエ変換と逆フーリエ変換には1対1の関係がありますので片方が分かれば,逆も言えます.

 \begin{eqnarray}
\bar{S}(f_k) &=& \sum_{j=-(N-1)}^{N-1} \left(1- \frac{|\, j \,|}{N}\right) C(\, j \,)\, e^{- i 2 \pi f_k j}
\end{eqnarray}
となったので,さらに,N \to \inftyとすれば,

 \begin{eqnarray}
S(f_k) &=& \sum_{j=-\infty}^{\infty} C(\, j \,)\, e^{- i 2 \pi f_k j}
\end{eqnarray}

になりそうですが,本当にそうでしょうか.問題は,

 \displaystyle \left(1- \frac{|\, j \,|}{N}\right) C(\, j \,)

の部分です.|\, j \,|がすごーく大きいときは,Nがそこそこ大きくても,

 \displaystyle \left(1- \frac{|\, j \,|}{N}\right) \approx 1

と近似するには無理があります.

 このような無理をしないために,|\, j \,|がすごーく大きくなる前に,|\, j \,|の増加にともない急速に

 \displaystyle C(\, j \,) \approx 0

となってくれる必要があるのです.つまり,

 「ラグkが大きくなるとき自己共分散関数C_kが急激に減衰し」

の条件が必要なのです.とはいえ,なぜ,自己共分散関数C_kフーリエ変換パワースペクトルを定義しているのかについては,私は十分に理解できていません.弱定常過程の前提ありきなのかなと想像しています.どなたか,詳しい方に教えていただきたいです.

練習問題 (すごく重要だと思います.お気に入りの問題です)

 次の確率過程ではない例

 \begin{eqnarray}
y [i]&=& A \sin \left(2 \pi f_k i \right) \quad (f_k = k/N)
\end{eqnarray}

の自己共分散関数 (減衰せずに振動し続けとる (下図参照))
 \begin{eqnarray}
C(\, j \,) &=& \frac{A^2}{2} \cos \left(2 \pi f_k j \right) 
\end{eqnarray}
を使って,パワースペクトルを計算してみてください.この自己共分散関数は,

 \begin{eqnarray}
C(\, j \,) &=& \lim_{N \to \infty}\frac{1}{N} \cdot A \sin \left(2 \pi f_k i + \theta \right) \cdot A \sin \left(2 \pi f_k (i+j)  + \theta  \right) \\
 &=& \lim_{N \to \infty} A^2 \left( \frac{1}{2} \cos\left(2 \pi f_k j \right) +  \frac{1}{N} \frac{\sin(\cdots) - \sin (\cdots)}{\sin (2 \pi f_k)} \right) \\
 &=& \frac{A^2}{2} \cos \left(2 \pi f_k j \right)
\end{eqnarray}
として,無理矢理計算しました.\thetaは任意の値です.

正弦波の自己共分散関数とパワースペクトル (f_k=50/NN=1000)

 長さNの時系列 \{ y [i] \}フーリエ変換を使ったパワースペクトルの計算結果 (これは,N \to \inftyで収束しないぞ!)は,

 \begin{eqnarray}
\hat{S}(f_k) = \left\{
\begin{array}{ll}
\displaystyle \frac{A^2}{4} N & (k=j, k=N-j) \\ \\
\ 0 & ({\rm otherwise})
\end{array}
\right.
\end{eqnarray}

となります (上の図右側).「えー,なんでNがかかってんの?」と感じること,疑問に思うことが大切です.パワースペクトルというのは,正式にはパワースペクトル密度です.パワースペクトル密度の面積が,時系列の分散に一致するのでした.この場合,長方形で近似することにすれば (三角形でもいいですが),周波数の幅\Delta f=1/Nを使って,全面積は
 
  \displaystyle \frac{A^2}{4} N \cdot \frac{1}{N} + \frac{A^2}{4} N \cdot \frac{1}{N} = \frac{A^2}{2}

となります.これは,時系列の分散\displaystyle C(0)=\frac{A^2}{2}と一致します (正しそうです).

 ということで,上に示したパワースペクトル\hat{S}(f_k) と,自己共分散関数C(\, j \,) フーリエ変換 (ウィーナー=ヒンチンの定理)を使ったパワースペクトルの計算を比較してください (発散しちゃうぞ!).C(\, j \,) フーリエ変換が存在しないことに気づいてください.

 ちなみに,上で導いた関係式

 \begin{eqnarray}
\bar{S}(f_k) &=& \sum_{j=-(N-1)}^{N-1} \left(1- \frac{|\, j \,|}{N}\right) C(\, j \,)\, e^{- i 2 \pi f_k j}
\end{eqnarray}

で計算すると,
 \begin{eqnarray}
\sum_{j=-(N-1)}^{N-1} \left(1- \frac{|\, j \,|}{N}\right) \left\{ \frac{A^2}{2} \cos \left(2 \pi f_k j \right) \right\} e^{- i 2 \pi \frac{k}{N} j} &=& \left\{
\begin{array}{ll}
\displaystyle \frac{A^2}{4} N & (k=j, k=N-j) \\ \\
\ 0 & ({\rm otherwise}) \nonumber \\
\end{array}
\right.
\end{eqnarray}

となり,上に書いたパワースペクトル\hat{S}(f_k)と一致します.

おまけ

 作図に使ったRスクリプト

sig2 <- 1
a <- 0.9
par(mfrow=c(1,2))
########
# 解析的に計算した自己共分散関数
curve(sig2/(1-a^2)*a^abs(x),col=4,lwd=2,n=1000,xlim=c(-4*pi/log(a),4*pi/log(a)),xlab="k",ylab="C(k)",main="Autocovariance",xaxs="i")
# 解析的に計算したパワースペクトル
curve(sig2/(1-2*a*cos(2*pi*x)+a^2),col=4,lwd=2,n=1000,xlim=c(-1/2,1/2),xlab="f",ylab="S(f)",main="Power spectral density",xaxs="i")
k <- 50
N <- 1000
A <- 1
par(mfrow=c(1,2))
########
# 解析的に計算した自己共分散関数
curve(1/2*cos(2*pi*k/N*x),xlim=c(-N/k*4,N/k*4),col=4,lwd=2,xlab="k",ylab="C(k)",main="Autocovariance",xaxs="i",n=1000)
# 解析的に計算したパワースペクトル
curve(0*x,col=4,lwd=2,n=100,xlim=c(-1/2,1/2),ylim=c(0,N/3),xlab="f",ylab="S(f)",main="Power spectral density",xaxs="i")
points(c(k/N,(N-k)/N-1),c(N/4*A^2,N/4*A^2),pch=16,col=4)
points(c(k/N,(N-k)/N-1),c(0,0),pch=16,col="white")
points(c(k/N,(N-k)/N-1),c(0,0),pch=1,col=4)

【時系列解析】長期記憶,長時間相関,フラクタル

時系列の長期記憶 (long memory),長時間相関 (long-range correlation),フラクタル性 (fractality)について整理しておきます.

長期記憶過程 (long memory process)

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

\begin{eqnarray}
\sum_{k=0}^{\infty} \left |C(k) \right| = \infty
\end{eqnarray}

 となる確率過程です.何で全部足すの?という疑問がわくと思います.以前に1次自己回帰過程で,相関が0とみなせるまでの時間の話をしました.
 chaos-kiyono.hatenablog.com

1次自己回帰過程のように,相関が0とみなせるまでの時間が見積もれる (特徴的な時間スケールを持つ)場合は,上の積分は有限の値になります.その値が記憶の長さと比例します.上の積分が有限の値になる場合は,短期記憶過程です.

 長期記憶過程では,

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\quad (0 < \gamma < 1, 0 < \beta < 1)

の関係が成り立ちます.とはいえ,\gammaが0に近かったり,1に近かったりすると,実際の時系列からの推定では,C(k) \propto |k|^{-\gamma}からかなりずれています.

 長期記憶過程の代表的モデルは,

  • ハースト (Hurst)指数Hが,0.5 < H < 1の領域の非整数ガウスノイズ (fractional Gaussian noise)
  • ARFIMA (Autoregressive fractionally integrated moving average)過程のうち,パラメタd0 < d < 0.5の領域のARIFIMA(0,d,0)過程

があります.非整数ガウスノイズは,非整数ブラウン運動 (fractional Brownian motion)の増分過程です (差分をとった時系列)ので,ハースト (Hurst)指数Hは対応する,非整数ブラウン運動の値で指定されます.

長時間相関過程 (long-range correlated process)

 長時間相関過程は長期記憶過程を含みます.長期記憶過程は,正の長時間相関と同じです.長時間相関には,負の長時間相関もあります.長時間相関過程では,パワースペクトル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は成り立ってます.

 負の長時間相関過程では,

 \begin{eqnarray}
\sum_{k=0}^{\infty} \left |C(k) \right|
\end{eqnarray}

は有限の値になりますので,長期記憶過程ではありません.

 長時間相関の代表的モデルは,

  • 非整数ガウスノイズ (0 < H < 0.5, 0.5 < H < 1)
  • ARIFIMA(0,d,0)過程 (-0.5 < d  <0,  0< d< 0.5)

があります.非整数ガウスノイズは,非整数ブラウン運動 (fractional Brownian motion)の増分過程です (差分をとった時系列).

フラクタル過程

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

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

となります.非定常なので,自己共分散は,C(k)のように,ラグkの関数として書くことはできません (弱定常ではないし,上の2つの過程と比較できるような自己共分散関数は存在しません).

 ハースト (Hurst)指数H\betaには,

 \beta = 2H + 1

の関係が成り立ちます.

 フラクタル過程の代表的モデルは,

です.長時間相関過程を積分したものも,フラクタル過程になります.

 S(f) \propto |f|^{-\beta}を示す時系列を積分すると,S(f) \propto |f|^{-(\beta+2)}となることを,以前,説明しました.
chaos-kiyono.hatenablog.com
chaos-kiyono.hatenablog.com
ですので,何回も時系列を積分すれば,S(f) \propto |f|^{-\beta} \quad (\beta > 3)の時系列も作ることができます.

解析方法

 ここで説明した特徴をもつ時系列については,自己共分散関数C(k)の推定はほとんどの場合,役に立ちません.時系列の解析では,その他の方法を使ってください.

  • パワースペクトルの推定 (スペクトル解析)は,どの場合でも使うことができます.ただし,トレンド成分を除くなどの前処理が必要です.
  • ランダムウォーク解析,detrended fluctuation analyis (DFA), deterending moving average analysis (DMA)には,スケーリング指数 (多くの場合\alphaで表す)を推定できる下限と上限があります.
  • DFA,DMAの高次の方法を使えば,時系列に含まれる多項式トレンド (多項式で近似できるトレンド)を除くことができます.
  • ウェーブレット解析でも,スケーリング指数を推定できます.バニッシングモーメントを持つウェーブレットを使えば,多項式トレンドも除去できます.

【確率過程・時系列解析】自己回帰過程の特性方程式の根 (数値解)

前回は,自己回帰過程の特性方程式の話をしました.
chaos-kiyono.hatenablog.com
つまり,M次自己回帰過程

  \begin{eqnarray} y[n] &=& a_1 y[n-1] + a_2 y[n-2] + \cdots + a_m y[n-m] + w[n] \\ 
&=& \sum_{m=1}^M a_m y[n-m] + w[n] \end{eqnarray}

について,

特性方程式
 \displaystyle 1- a_1 z- a_2 z^{2} - \cdots - a_M z^{M} = 0
の解 (根)の絶対値がすべて1より大きければ,自己回帰過程は発散しない.

ということを説明しました.今回は,Rをつかって実際に特性方程式の根を求めてみます.

 使うスクリプトは,以下です.

# 自己回帰過程の係数を与える
# a <- c(a1,a2,a3, ... )
a <- c(0.9*sqrt(3),-0.81)
##############################
# 根を計算
roots <- polyroot(c(-1,a))
##############################
# 根を複素数平面にプロット
r <- abs(roots)
r.max <- max(1,r)
x <- Re(roots)
y <- Im(roots)
par(pty="s")
plot(x,y,pch=4,col=2,xlim=c(-r.max,r.max),ylim=c(-r.max,r.max),xlab="Re(z)",ylab="Im(z)")
# 単位円の描画
curve( sqrt(1-x^2),xlim=c(-1,1),col=8,add=TRUE)
curve(-sqrt(1-x^2),xlim=c(-1,1),col=8,add=TRUE)
# 根
roots

 このスクリプトでは,例として,2次自己回帰過程で,a_1=0.9 \sqrt{3}a_2=-0.81の場合を計算しています.実行すると以下の図が描かれます.円は半径1の単位円|z|=1です.赤いxが根ですので,xが円の外側にあれば,この過程は発散しません.

2次自己回帰過程の根を複素数平面に描いた図.赤いxが根の位置.

 以前の記事 
chaos-kiyono.hatenablog.com
で例として使った,a <- c(2.373596,-2.528090,1.404494)で試してみると,

3次自己回帰過程の根の位置

となります.すいません,発散してしまいます.安定になる条件を,間違えました.

# 3次自己回帰過程の係数
a <- c(2.373596,-2.528090,1.404494)
# 繰り返し回数
N <- 30
########################################
y <- numeric(N)
w <- rnorm(N)
y[1:3] <- w[1:3]
for(n in 4:N){
 y[n] <- a[1]*y[n-1]+a[2]*y[n-2]+a[3]*y[n-3]+w[n]
}
plot(1:n,y,"l",las=1)

数値実験してみると,下の図のようになります.この例では,マイナス側に発散しています (常にプラスかマイナスに発散します).

数字実験の例.発散しちゃいます.

【確率過程・時系列解析】自己回帰過程の特性方程式の解と定常性

今回は,自己回帰過程の特性方程式の解と定常性の話です.特性方程式の解が,単位円の内側とか外側とかいう話の説明をします.

M次自己回帰過程 (以下のw[n] は,平均0,分散\sigma^2の白色ノイズです)

  \begin{eqnarray} y[n] &=& a_1 y[n-1] + a_2 y[n-2] + \cdots + a_m y[n-m] + w[n] \\ 
&=& \sum_{m=1}^M a_m y[n-m] + w[n] \end{eqnarray}

をラグオペレータLを使って書けば,

  \begin{eqnarray} y[n] &=& \sum_{m=1}^M a_m L^{m} y[n] + w[n] \\
\left(1- \sum_{m=1}^M a_m L^{m} \right) y[n] &=&  + w[n] \\
y[n] &=& \frac{1}{\displaystyle 1- \sum_{m=1}^M a_m L^{m}}w[n]  \end{eqnarray}

と書けます.これは,以前にパワースペクトルを求めるときにも書きました.今回は分母のLを,なぜかzで置き換えた式,

 \displaystyle 1- \sum_{m=1}^M a_m z^{m} = 0

の解 (根)についてのお話です.(根って何?解と違うの?という疑問があるかたは,「根」では重解も,1個ずつメンバとして,大切に数えてあげるのだと覚えてください.同じものです)

 ということで,この式を特性方程式と呼びます.

\displaystyle 1- a_1 z- a_2 z^{2} - \cdots - a_M z^{M} = 0

ほとんどの教科書や,ホームページでは,この方程式の解(根)が,すべて単位円の外側にあれば (絶対値が1よりも大きければ),自己回帰過程は安定とか,弱定常とかさらっと書いてあります.何で?って疑問に思う人は多いと思います.ということで,説明します.

 特性方程式の根 (複素数もありえます)を

 \lambda_1, \lambda_2, \cdots, \lambda_M

とすれば (根なので必ずメンバはM個です),特性方程式は,

 \displaystyle \left( 1 - \frac{z}{\lambda_1} \right) \left( 1 - \frac{z}{\lambda_2} \right) \cdots \left( 1 - \frac{z}{\lambda_M} \right) = 0

因数分解できます.これを,元の自己回帰過程の式に反映させると (zをLにして),

  \begin{eqnarray} y[n] &=& \frac{1}{\displaystyle 1- \sum_{m=1}^M a_m L^{m}}w[n] \\
&=& \frac{1}{\displaystyle \left( 1 - \frac{L}{\lambda_1} \right) \left( 1 - \frac{L}{\lambda_2} \right) \cdots \left( 1 - \frac{L}{\lambda_M} \right) }w[n]  \end{eqnarray}

となります.右辺の分数の部分は必ず部分分数分解できるので (これは認めてください),

  \begin{eqnarray} y[n] &=& \frac{1}{\displaystyle \left( 1 - \frac{L}{\lambda_1} \right) \left( 1 - \frac{L}{\lambda_2} \right) \cdots \left( 1 - \frac{L}{\lambda_M} \right) }w[n] \\ 
&=& \left( \frac{c_1}{\displaystyle 1 - \frac{L}{\lambda_1}} + \frac{c_2}{\displaystyle 1 - \frac{L}{\lambda_2}}+ \cdots + \frac{c_M}{\displaystyle 1 - \frac{L}{\lambda_M}} \right)  w[n] \\ \\
&=& \frac{c_1}{\displaystyle 1 - \frac{L}{\lambda_1}} w[n]+ \frac{c_2}{\displaystyle 1 - \frac{L}{\lambda_2}} w[n]+ \cdots + \frac{c_M}{\displaystyle 1 - \frac{L}{\lambda_M}} w[n] 
\end{eqnarray}

となります. \{c_1, c_2, \cdots\}には定数が入ります.

 \displaystyle c_k = \left. \frac{ \displaystyle \left( 1 - \frac{z}{\lambda_k} \right)}{\displaystyle \left( 1 - \frac{z}{\lambda_1} \right) \left( 1 - \frac{z}{\lambda_2} \right) \cdots \left( 1 - \frac{z}{\lambda_M} \right) } \right|_{z = \lambda_k}

右辺の各項が収束すれば,この過程は安定とか,定常とか,そういうことです.要は,y[n]が,n \to \inftyで発散しなければ,めでたしめでたしです.

 ということで,

  \displaystyle \frac{c_k}{\displaystyle 1 - \frac{L}{\lambda_k}} w[n]

について考えてみます.前にもやりましたが,分数の部分について,マクローリンをまねて展開します.

  \begin{eqnarray} \frac{c_k}{\displaystyle 1 - \frac{L}{\lambda_k}} w[n] &=& \displaystyle \left(1 + \frac{L}{\lambda_k} + \left( \frac{L}{\lambda_k} \right)^2 + \cdots  \right) w[n] \\ 
&=& \displaystyle w[n] + \frac{L}{\lambda_k} w[n] + \left( \frac{L}{\lambda_k} \right)^2 w[n] + \cdots \\ 
&=& \displaystyle w[n] + \frac{1}{\lambda_k} w[n-1] + \frac{1}{\lambda_k^2} w[n-2] + \cdots \end{eqnarray}
 
これの分散を計算すれば (分散を {\rm Var}で表します)

 \begin{eqnarray} && {\rm Var} \left[ w[n] + \frac{1}{\lambda_k} w[n-1] + \frac{1}{\lambda_k^2} w[n-2] + \cdots \right] \\
&=&  {\rm Var} \left[  w[n] \right] + {\rm Var} \left[ \frac{1}{\lambda_k} w[n-1] \right] + \cdots \\ 
&=&  \sigma^2 + \frac{\sigma^2}{\lambda_k^2}  + \frac{\sigma^2}{\lambda_k^4}  + \cdots  \end{eqnarray}

となります.公比が,1/\lambda_k^2等比数列になっています (\lambda複素数の場合,1/|\lambda_k|^2にした方がいいかも).これが,収束する条件は,

 \displaystyle \left| \frac{1}{\lambda_k^2}\right| < 1

なので,

 \displaystyle \left|\lambda_k \right| > 1

であれば (複素数でも実数でも,絶対値が1より大きければ),上で計算していた分散は,

 \displaystyle \frac{ \sigma^2}{1-1/|\lambda_k|^2}

に収束します.これと似た計算がM個ありますが,特性方程式の解 (根)の絶対値が,すべて1より大きければ,全体も収束します.ということで,

特性方程式
 \displaystyle 1- a_1 z- a_2 z^{2} - \cdots - a_M z^{M} = 0
の解 (根)の絶対値がすべて1より大きければ,自己回帰過程は発散しねー

ということがわかります.途中がわからなければ自分で勉強してください.他にもっとエレガントな方法や,間違いがあるかもしれません.他人を安易に信じないでください.

【心拍変動解析】RR間隔時系列の再サンプリング (等間隔サンプリング化)

心拍変動解析では,パワースペクトルの推定を使った心拍変動指標がよく使われます.HFパワーとか,LFパワーとか呼ばれるやつです.今回は,これら指標について説明しません.心拍変動指標を計算する前に,絶対に理解しておかなければならないことがあるからです.それは,心拍変動時系列の再サンプリング (resampling: リサンプリング)です.今回は,心拍変動のパワースペクトルを計算する前処理の一つである,RR間隔時系列の再サンプリングについて説明します.

RR間隔とは

 そもそもRR間隔とは,下の図のように,心電図に見られるR波とR波の間隔のことです.心拍変動解析では,正常同調律を構成するR波のみを対象として,RR間隔のゆらぎを解析します.不整脈がある場合は,不整脈がなかったとしたらどの位置にR波が来るかを予想して,RR間隔を補正します.これも,心拍変動解析に必要な前処理です (今回は説明しません).

心電図とRR間隔の関係

RR間隔の再サンプリング

 パワースペクトルは,周波数fの関数です.通常,fの単位はHzで,これは,時間 (秒)の逆数の次元 (単位)を持っています.

 RR間隔のグラフを描くとき横軸の単位は何でしょうか?横軸は,最初のRR間隔,2番目のRR間隔,3番目のRR間隔というふうに,拍動の順番になっています.このままでは,横軸の単位が時間になっていないので,パワースペクトルを計算しても,周波数fの単位はHzになりません.

 ということで,周波数fの単位をHzにしたいので,RR間隔データを描く横軸の単位を時間にする必要があります.そのために,下の図のように,R波の発生時刻にその時刻のRR間隔を配置します.

RR間隔を時間軸上に配置

 これで,横軸は時間になりました.しかし,RR間隔は一定値ではないので,このままでは不等時間間隔データになっています.一般的なフーリエ変換では,等時間間隔にサンプリングされたデータが前提です.したがって,追加の作業として,等時間間隔でサンプリングされたデータにする必要があります.そこで,再サンプリング (resampling)の登場です.

 下の図のように,実際の点の間を補間し,一定間隔になる位置の点をとっていきます.

等間隔サンプリング

 代表的な補間方法は,下の図のような階段補間(左)と線形補間(右)です.スプライン補間は響きがかっこいいですが,不自然に振動する場合があるのでお勧めしません.

階段補間(左)と線形補間(右)

Rで再サンプリング

 Rを使って実際にRR間隔時系列を再サンプリングしてみます.データは以下の記事にあるリンクからダウンロードできます.フォルダへのパスの確認法もこの記事を参照してください.

chaos-kiyono.hatenablog.com

 線形補間を使うと下の図のようになります.元のRR間隔時系列を補間した(線でつないだ)結果が青実線で,再サンプリングされた時系列が赤丸です.

線形補間の例

 階段補間を使うと下の図のようになります.

階段補間の例

 使ったRスクリプトは以下です.setwd("...")の部分は,"chaosRRI.csv"が保存されているフォルダのパスを指定してください.

線形補間のRスクリプト

# サンプリング間隔 [s]
dt <- 0.5
#####################################
setwd("...")
FN <- "chaosRRI.csv"
TMP <- read.csv(FN,header=TRUE)
RRI <- TMP$RRI
# 時刻の設定 [s]
time <- (cumsum(RRI)-RRI[1])/1000
#####################################
# 出力時刻の設定
t.out <- seq(0,tail(time,1),dt)
#####################################
# 再サンプリング (線形補間)
RRI.resamp <- approx(time,RRI,xout=t.out,method="linear")$y
#####################################
# プロット (一部のみ)
n.sub <- 20
plot(time[1:n.sub],RRI[1:n.sub],"l",col=4,xlab="time [s]",ylab="RRI [ms]",las=1)
points(t.out[t.out <= time[n.sub]],RRI.resamp[t.out <= time[n.sub]],pch=16,col=2)
legend("topright", legend=c("original RRI","resampled RRI"), col = c(4,2), pch = c(NA,16),lty=c(1,NA))
#####################################
# t.out: 再サンプリング時刻
# RRI.resamp: 再サンプリングされたRR間隔

階段補間のRスクリプト

# サンプリング間隔 [s]
dt <- 0.5
#####################################
setwd("...")
FN <- "chaosRRI.csv"
TMP <- read.csv(FN,header=TRUE)
RRI <- TMP$RRI
# 時刻の設定 [s]
time <- (cumsum(RRI)-RRI[1])/1000
#####################################
# 出力時刻の設定
t.out <- seq(0,tail(time,1),dt)
#####################################
# 再サンプリング (線形補間)
RRI.resamp <- approx(time,RRI,xout=t.out,method="constant")$y
#####################################
# プロット (一部のみ)
n.sub <- 20
plot(time[1:n.sub],RRI[1:n.sub],"s",col=4,xlab="time [s]",ylab="RRI [ms]",las=1)
points(t.out[t.out <= time[n.sub]],RRI.resamp[t.out <= time[n.sub]],pch=16,col=2)
legend("topright", legend=c("original RRI","resampled RRI"), col = c(4,2), pch = c(NA,16),lty=c(1,NA))
#####################################
# t.out: 再サンプリング時刻
# RRI.resamp: 再サンプリングされたRR間隔

【時系列解析】自己回帰過程のパラメタとパワースペクトルの関係

今回は,一般にm次自己回帰過程

  y[n] = a_1 y[n-1] + a_2 y[n-2] + \cdots + a_m y[n-m] + w[n]

のパラメタ \{a_1, a_2, \cdots, a_m\}と, w[n] の分散\sigma^2 (平均は0とします)がわかったときに,そのパワースペクトルを導きます.というか,パワースペクトルをとりあえず式で書きます.

 今回は,パワースペクトルパラメトリック推定を理解するための準備になってます.

 これまでの2次自己回帰過程の話で,すでに方法はわかったと思います.まず,基本事項 (使う道具)の確認です.

基本事項の確認

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

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

となります.つまり,Lが1つかかるごとに,時間が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)

となります.フーリエ変換すると,Le^{{\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}

です.

パワースペクトルの導出

 以下の自己回帰過程

  y[n] = a_1 y[n-1] + a_2 y[n-2] + \cdots + a_m y[n-m] + w[n]

を,ラグオペレータLを使って表すと,

 \begin{eqnarray}
 y[n] &=& a_1 L y[n] + a_2 L^2 y[n] + \cdots + a_m L^m y[n] + w[n] \\
 y[n] - \left(a_1 L y[n] + \cdots + a_m L^m y[n] \right) &=& w[n] \\
 \left(1 - a_1 L + \cdots + a_m L^m \right) y[n] &=& w[n] \\
 y[n] &=& \frac{1}{1 - a_1 L + a_2 L^2 + \cdots + a_m L^m} w[n] \\
\end{eqnarray}

 ここから,両辺のフーリエ変換をして,両辺の絶対値をとって,両辺を2乗して,両辺をNでわると.....
左辺は,

  y[n] が, Y(f_k) になって, |Y(f_k)| になって, |Y(f_k)| ^2になって, \displaystyle \frac{|Y(f_k)|^2}{N} になるから,

最後は, y[n] パワースペクトルS_y(f_k)になります.

左辺は, w[n] パワースペクトル\sigma^2になります.L \displaystyle e^{{\bf i} 2 \pi f}になります.

 ということで,

 \displaystyle 
 S_y (f_k) = \frac{\sigma^2}{\left|1 - a_1 e^{{\bf i} 2 \pi f} + a_2 e^{{\bf i} 4 \pi f} + \cdots + a_m e^{{\bf i} 2 \pi f m} \right|^2}

です.

 こんなこと,どの時系列解析の教科書にも書いてあるので,つまらないですが,ここからは, \{a_1, a_2, \cdots, a_m\}と, w[n] の分散\sigma^2が与えられたときのパワースペクトルをRで描いてみます.ポイントは,

  \{a_1, a_2, \cdots, a_m\}\sigma^2が決まるとパワースペクトルが計算できる

ということです.

下のスクリプトで,

a <- c(2.373596,-2.528090,1.404494)

の部分の数値や長さを変えていろいろと実験してみてください (すいません,この例は不安定で,良くない係数でした).

Rでパワースペクトルを描いてみる.
# 自己回帰過程のパラメタを与える
# a <- c(a_1,a_2,...,a_m)とします
a <- c(2.373596,-2.528090,1.404494)
# 白色ノイズの分散
sig2 <- 1
####################################
m <- length(a)
# 周波数の設定 (length.outで点数を指定)
f <- seq(0,1/2,length.out=100)
n.f <- length(f)
# パワースペクトルの計算
Sf <- c()
for(k in 1:n.f){
 Sf[k] <- 1/Re((c(1,-a) %*% exp(2*pi*f[k]*(0:m)*1i))*(c(1,-a) %*% exp(-2*pi*f[k]*(0:m)*1i)))
}
# パワースペクトルのプロット
par(mfrow=c(1,2))
plot(f,Sf,"l",col=4,lwd=2,xlab="f",ylab="S(f)",main="線形目盛")
plot(f[f>0],Sf[f>0],"l",log="xy",col=4,lwd=2,xlab="f",ylab="S(f)",main="両対数目盛")

注意

 上で例として使った係数は,良くありません.良くないというのは,この過程が不安定,つまり,値が発散してしまうからです.このことについては,以下の記事を参考にしてください. chaos-kiyono.hatenablog.com chaos-kiyono.hatenablog.com

【Rで時系列解析】心拍変動データを読み込み

今回は実際の時系列データを読み込んでグラフを描いでみます.

 濃厚接触者として自宅待機中(現在)の私の心拍変動 (RRI間隔)時系列をサンプルとして使ってください.以下のリンクからダウンロードできます.
https://drive.google.com/file/d/1i0VIooI1-IKAUI9fHi0pEwS9AF6AsHYk/view?usp=sharing
(私は,コロナ関連の研究もしているので,自分が感染する前後の心拍変動の変化を記録しようと,ずっと心拍変動を計測しています.しかし,PCR陰性でしたので,今回は感染しない可能性が高いです.)

基本事項の確認

 上のリンクにあるファイルをダウンロードし,自分のパソコンに保存してください.そして,このファイルを保存したフォルダのパスを確認してください (下図).

Windowsでのアドレスバーの位置.アドレスバーの部分をクリックすればパスが表示されます.これをコピーしてください.

ここでは,フォルダのパス(仮)を

"C:\Users\chaos\Documents\時系列解析"

読み込むファイル名 (上のリンクからダウンロードしたファイル)を

"chaosRRI.csv"

と仮定して話を進めます.フォルダのパスは仮ですので,各自の環境に合わせて変える必要があります.そのままでは,絶対に動きません.

ファイルを読み込んで,とりあえずグラフ

 ファイルを読み込んで,グラフを描いてみます.以下のRスクリプトを実行すると,下の図が描かれます.

RRI時系列のプロット

使ったRスクリプトは以下のものです.注意点は,

 setwd("C:\Users\chaos\Documents\時系列解析")

"C:\Users\chaos\Documents\時系列解析"の部分を,自分の環境に合わせて書き換えることです.フォルダの区切り文字「\」を,「/」にかえる必要あります.

# "chaosRRI.csv"を保存したフォルダのパスを指定 【重要】 「\」を「/」にかえること
setwd("C:\Users\chaos\Documents\時系列解析")
# ファイルの読み込み.FNにファイル名を代入
FN <- "chaosRRI.csv"
TMP <- read.csv(FN,header=TRUE)
# 時刻データの処理
time <- as.POSIXct(TMP$time,tz="Japan")
# RRIデータの処理 (as.numericはなくても大丈夫)
RRI <- as.numeric(TMP$RRI)
# グラフの描画
plot(time,RRI,"l",col=4,xlab="time",ylab="RR interval [ms]",xaxs="i")

使用するコマンドについて,今後,すこしずつ説明していきます.

心拍変動って何?RR間隔って何?

 今回の時系列データについて簡単に説明します.このデータは,心臓の拍動時間間隔を表しています.心電図を測ると下の右の図のような波形がえられます.心電図にある上向きの鋭いピークがR波と呼ばれるものです.この波形は心室が電気を発生したときに現れます.このR波のピークの間隔がRR間隔です.

心電図のR波

 今回描いたRR間隔時系列のグラフは,意外とギザギザしているな,と感じたかもしれません.これは,私が不健康だからではありません.健康な人は,RR間隔が結構ばらつく(ゆらぐ)のです.逆に,年を取ったり,病気があったりすると,このゆらぎは小さくなります.

心臓の拍動リズムのゆらぎは,「心拍変動 (heart rate variability)」と呼ばれています

私は20年にわたり心拍変動の研究をしてきました.心拍変動の特徴を数理的に読み解くと,健康や病気についていろんなことがわかります.興味があれば,大学院から私の研究室にきてください.