非整数ブラウン運動の説明用のアニメを作ってみました．
高解像度な，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")
