二次元ヒートカーネル(熱核)の時間発展の可視化
二次元ヒートカーネルなんていうと難しく聞こえるけど、要するにラプラス方程式の解である、時間発展するガウシアンの可視化です。物理系だと微分方程式の基本解という意味で、熱核じゃなくてグリーン関数って呼んでるものと等価（なはず）。
具体的には
1. 二次元のドリフト付ブラウン運動を複数生成する
2. 1の結果に対し、各時点において二次元カーネル密度推定を実行
3. 2.の結果をアニメーションとして可視化
という流れになる。まず必要となるライブラリを読み込んでおく。なかったらinstall.packages関数しましょう。
library(dplyr) library(MASS) library(pipeR) library(animation)
次に二次元ドリフト付ブラウン運動を生成させる。これは大体以下のコードを使っている。
#二次元ドリフト付ブラウン運動 #タイムステップ数 size <- 35 #適当な運動を作る randomwalk2D <- function(){cbind(x=cumsum(1+rnorm(size)), y=cumsum(1+rnorm(size)))} #本シミュレーションを走らせて結果をdata.frameに統合 df <- lapply(1:100, function(i)data.frame(randomwalk2D(), group=LETTERS[i], time_step=1:size)) %>>% rbind_all
次に、後の描画範囲指定で使う上下限の範囲を設定しつつ、
でやったように、二次元カーネル回帰させる。
#x, yの上下限計算 xlim <- df %>>% (x) %>>% range ylim <- df %>>% (y) %>>% range #カーネル回帰 kde2ds <- lapply(1:size, function(i){ df %>>% filter(time_step==i) %>>% with({ #バンド幅の計算 band_width <- c(bandwidth.nrd(x),bandwidth.nrd(y)) #二次元カーネル密度推定 kde2d(x, y, band_width, n=50, lims=c(xlim, ylim)) }) })
んで、animationパッケージの中ではImageMagickを使うので、それはそれで別途いれておく。入れておくとして、Rとの連携の際には一応こんなかんじでパスを設定しておくとうまくいった。昔使った時にはいらなかった記憶があるんだけどなぁ…
#animationパッケージ用の設定 ani.options(convert = 'C:/Program Files/ImageMagick-6.9.0-Q16/convert.exe')
以下のコードで二次元ヒートカーネル(熱核)の時間発展の可視化を実行している。z方向の値(確率密度)がすべての時間で共通になるようにうまくスケーリングしている。パラメーターの設定は黒魔術的に適当にやった。
#スケーリングあり ani.record(reset=TRUE) zlim_max <- c(kde2ds[[1]]$z, kde2ds[[size]]$z) %>>% max %>>% `*`(0.1) for(i in 1:size){ kde2ds[[i]] %>>% image(xlim=xlim, ylim=ylim, breaks=seq(from=0, to=zlim_max, length.out=13)) ani.record() } ani.options(interval = 0.1) ani.replay() saveGIF(ani.replay(), movie.name="2dheatkernel_scale.gif")
同様に今度はスケーリングなしのケース。従って、、「(異なる時点(アニメーションの1コマ))における図での明るさが等しい = 確率密度が等しい」とはならない。
#スケーリングなし ani.record(reset = TRUE) for(i in 1:size){ kde2ds[[i]] %>>% image(xlim=xlim, ylim=ylim) ani.record() } ani.options(interval = 0.1) ani.replay() saveGIF(ani.replay(), movie.name="2dheatkernel_noscale.gif")
頑張って3Dにしてみた。色がなんかおかしいぞ・・・？
pal <- colorRampPalette(c("red", "yellow"))(50) animation3d <- function(){ zmax <- max(kde2ds[[1]]$z) for(i in 1:size){ kde2ds[[i]] %>>% with({ col.ind <- cut(z,50) persp(x,y,z, zlim=c(0, zmax), col=pal[col.ind], theta=45)}) } } saveGIF(animation3d(), movie.name="2dheatkernel_3d.gif", interval=0.1)
