【Rでフラクタル】非整数ブラウン運動Zoomingアニメの作成

非整数ブラウン運動は，横と縦方向の拡大倍率を， $\lambda$ 対 $\lambda^H$ とすれば，時系列の統計的性質が拡大前後で同じになる過程です． $H$ は，ハースト（Hurst）指数と呼ばれるパラメタで，非整数ブラウン運動は $H$ の値で特徴付けられます．数式で理解する前に，非整数ブラウン運動のイメージを持ってもらうことが，今回の目的です．

例えば， $H=0.8$ では，以下の図のような感じです．図では，サンプルの軌道をだんだん拡大してみています．横と縦方向の拡大倍率が，異なる拡大になっていることに注意してください．

この図（GIFアニメは）以下のRスクリプトで作成しました．ここでは，中点法という方法[Fournier et al., Communications of the ACM, 25, 371 (1982)]で非整数ブラウン運動の軌道を作成しています（詳細は近いうちに説明します）．

#install.packages("magick")
###Rでgifアニメーションplotの作成
#ライブラリの読み込み
require("magick")
################
# ハースト指数の設定
H <- 0.8
# 中点法の最初のステップ
n.step <- 12
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)
i.x <- round(N/2)
x0 <- x[i.x]
y0 <- y[i.x]
d00.x <- -x0
d00.y <- -y0+min(y)-1
d01.y <- 1+max(y)-y0
d10.x <- 1-x0
r.x <- 1/(2^(1/16))
r.y <- r.x^H
################
img = image_graph()
par(pty="m")
# 余白をなしに設定
par(mar=c(0,0,0,0))
par(ps = 18)
#####################################
n.d <- 0
for(k in 1:10){
for(i in 1:16){
n.d <- n.d + 1
plot(x,y,"l",xlim=c(x0+r.x^n.d*d00.x,x0+r.x^n.d*d10.x),ylim=c(y0+r.y^n.d*d00.y,y0+r.y^n.d*d01.y),col=2,bty="n",xaxt="n",yaxt="n",xlab="",ylab="")
i.x1 <- min(which(x < x0+r.x^n.d*d00.x))
i.x2 <- min(which(x > x0+r.x^n.d*d10.x))
x.TMP <- x[i.x1:i.x2]
y.TMP <- y[i.x1:i.x2]
remove(x)
remove(y)
x <- x.TMP
y <- y.TMP
}
###
 n.x <- length(x)
 n.d <- n.d + 1
 h <- (x[2]-x[1])/2
 x.tmp <- x[1:n.x]
 y.tmp <- y[1:n.x]
 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(x0+r.x^n.d*d00.x,x0+r.x^n.d*d10.x),ylim=c(y0+r.y^n.d*d00.y,y0+r.y^n.d*d01.y),col=2,bty="n",xaxt="n",yaxt="n",xlab="",ylab="")
###
i.x1 <- min(which(x < x0+r.x^n.d*d00.x))
i.x2 <- min(which(x > x0+r.x^n.d*d10.x))
x.TMP <- x[i.x1:i.x2]
y.TMP <- y[i.x1:i.x2]
remove(x)
remove(y)
x <- x.TMP
y <- y.TMP
}
#描画を閉じる
dev.off()
#gifの速度を設定
animation = image_animate(img, fps=10, loop=0)
#保存
image_write(animation, path = "fBm_zoom08c9.gif", format = "gif")

PCパワーによっては，実行時間が何分もかかるので，注意してください．

ハースト指数は以下の箇所に， $0.5 \le H < 1$ の範囲で設定してください（ $0 < H < 0.5$ の範囲も可能ですが，アニメとしてイマイチです）．