【Ｒで時系列解析】Detrended Fluctuation Analysis (DFA)を使った長時間相関スケーリング解析

今回はDetrended Fluctuation Analysis (DFA)をＲでやってみます．パッケージもあると思いますが，パッケージなしで実現します．

まずは，解析する時系列の準備です．これは，longmemoパッケージを使います．分析したい時系列がある場合は，このRスクリプトは無視して，xに時系列をベクトルとして格納してください．

# パッケージのインストールが必要な場合は以下を実行
# install.packages("longmemo")
# パッケージのロード
require(longmemo)
###################################
# パラメタの設定：ハースト指数(0 < H < 1)
H <- 0.6
# 時系列の長さ
n <- 2^10
####################################
# 非整数ガウスノイズの生成 
x <- simFGN0(n,H)

以下は，DFAの説明用にグラフの作成も一緒にやるRスクリプトです．

#
# xを時系列ベクトルとする
#
N <- length(x)
###############################
#　平均0化
x <- x - mean(x)
# 積分
y <- cumsum(x)
###############################
# 作図用のパラメタ
par(mfrow=c(3,2))
par(pty="m",cex=1.1,las=1,cex.lab=1.2)
par(mar=c(5,4,2,2))
# 解析対象の時系列
plot(1:N,x,"l",col=4,xaxs="i",xlab="i",ylab="x[i]",lwd=2,main="Original time series")
plot(1:N,y,"l",col=4,xaxs="i",xlab="i",ylab="y[i]",lwd=2,main="Integrated series")
###############################
# scalesの設定
scale.min <- 8
scale.max <- N/4
scale <- unique(round(10^seq(log10(scale.min),log10(scale.max),0.05)))
n.scale <- length(scale)
# プロットする3つのスケールの指定
p.scale <- c(5,round(5+(n.scale-5)/2),n.scale-2)
###############################
F2 <- c()
for(i in 1:n.scale){
 i.sub <- seq(1,N,scale[i])
 n.sub <- length(i.sub)-1
 F2[i] <- 0
 if(is.element(i, p.scale)){
  plot(1:N,y,"l",col=4,xaxs="i",lwd=2,xlab="i",ylab="y[i]",main=paste("scale =",scale[i]))
  abline(v=i.sub-1,lty=2,col=8,lwd=1)
 }
 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(pty="s")
fit <- lm(log10(F2)/2~log10(scale))
plot(log10(scale),log10(F2)/2,xlab="log10 n",ylab="log10 F(n)",main=paste("slope =",format(fit$coefficients[[2]],digit=2)))
points(log10(scale[p.scale]),log10(F2[p.scale])/2,pch=16,col=2)
abline(fit,lty=2,col=2,lwd=2)

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

上の図では，左上に解析する時系列x，その右が積分時系列，下の時系列は，各スケールで区分的に直線（赤実線）を当てはめた結果，右下がスケーリング解析（ $\log_{10} n$ に対する $\log_{10} F(n)$ ）の結果です．

実際の時系列解析では，こんなに図は必要ないので，スクリプトはこんな感じになります．

#
# xを時系列ベクトルとする
#
N <- length(x)
###############################
#　平均0化
x <- x - mean(x)
# 積分
y <- cumsum(x)
###############################
# scalesの設定
scale.min <- 8
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,xlab="log10 n",ylab="log10 F(n)",main=paste("slope =",format(fit$coefficients[[2]],digit=2)))
abline(fit,lty=2,col=2,lwd=2)

DFAの参考文献については，初心者はこれかな，
Frontiers | Detrended Fluctuation Analysis: A Scale-Free View on Neuronal Oscillations