[[R News 2008年2号 (PDF)|http://cran.r-project.org/doc/Rnews/Rnews_2008-2.pdf]] や [[useR!2008での面白そうな発表をいくつか上げてみた|http://d.hatena.ne.jp/syou6162/20081004/1223092624]]で紹介されていたR の [[animation library|http://animation.yihui.name/]] を使ってみたよ。
まずは、非線形モデル y = mesor + a * cos(2*pi*(t-acrophase)/P) にN(0,0.5)のガウシアンノイズを加えたデータを用意し、これを測定されたデータと考える。そのデータに対して非線形回帰するために、Nelder-Mead法で目的関数sum( (y-yhat)^2 )を最大化する。いわゆる cosinor analysis ですね。その計算過程を各ステップごとにプロットし、animation libraryを使ってアニメ化する。
アニメーションの動きが激しくなるように、わざと初期値をおおげざにはずしてある。青が測定データで、緑が予測した回帰曲線。緑のほうの線が動いてみえるはず。収束間近ではあまり動かないのでしばらくみつめていると、ループしてまた初期値から計算しなおす。
{{image 0, 'cosinor analysis, cosin fitting', nil, [480,480]}}
コードは以下の通り。
<<<
library(animation)
t <- seq(0,48,1)
mesor <- 0
a <- 3
acrophase <- 12
P <- 24
y <- mesor + a * cos(2*pi*(t-acrophase)/P)
y.noise <- mesor + a * cos(2*pi*(t-acrophase)/P) + rnorm(length(t), mean=0, sd=1)
resid <- function(par) {
mesor <- par[1]
a <- par[2]
acrophase <- par[3]
P <- par[4]
yhat <- mesor + a * cos(2*pi*(t-acrophase)/P)
matplot(t,
matrix(c(y.noise, yhat), ncol=2),
col=c("#1E5692", "#3E9A3B"),
type="l", lty=1, pch=21, lwd=5, cex=2,
ylim=c(-5,5))
sum( (y-yhat)^2 )
}
saveMovie(optim(c(5,10,2,18), resid), interval = 0.01, movietype = "gif", outdir = getwd())
>>>
ようするにユーザはプロットを何枚か出力するような関数を用意するだけ。それをsaveMovie()関数に渡せば、内部でImageMagickやらSWF Toolsががんばってくれてアニメになる、という仕組みになっている。ImageMagickを MacPortsあたりで入れておく必要がある。
追記: あやまって hatena_bookmark_nocomment.rb を有効にしていたせいで、はてブのコメント一覧が見えなくなってしまいました。このプラグインを無効にしたのですが、最初にブクマされた時点で非表示に設定されると公開することができないようです。コメント頂いたかた、ごめんなさい。見れないので直接、ここにコメントして頂けると幸いです><
まずは、非線形モデル y = mesor + a * cos(2*pi*(t-acrophase)/P) にN(0,0.5)のガウシアンノイズを加えたデータを用意し、これを測定されたデータと考える。そのデータに対して非線形回帰するために、Nelder-Mead法で目的関数sum( (y-yhat)^2 )を最大化する。いわゆる cosinor analysis ですね。その計算過程を各ステップごとにプロットし、animation libraryを使ってアニメ化する。
アニメーションの動きが激しくなるように、わざと初期値をおおげざにはずしてある。青が測定データで、緑が予測した回帰曲線。緑のほうの線が動いてみえるはず。収束間近ではあまり動かないのでしばらくみつめていると、ループしてまた初期値から計算しなおす。
{{image 0, 'cosinor analysis, cosin fitting', nil, [480,480]}}
コードは以下の通り。
<<<
library(animation)
t <- seq(0,48,1)
mesor <- 0
a <- 3
acrophase <- 12
P <- 24
y <- mesor + a * cos(2*pi*(t-acrophase)/P)
y.noise <- mesor + a * cos(2*pi*(t-acrophase)/P) + rnorm(length(t), mean=0, sd=1)
resid <- function(par) {
mesor <- par[1]
a <- par[2]
acrophase <- par[3]
P <- par[4]
yhat <- mesor + a * cos(2*pi*(t-acrophase)/P)
matplot(t,
matrix(c(y.noise, yhat), ncol=2),
col=c("#1E5692", "#3E9A3B"),
type="l", lty=1, pch=21, lwd=5, cex=2,
ylim=c(-5,5))
sum( (y-yhat)^2 )
}
saveMovie(optim(c(5,10,2,18), resid), interval = 0.01, movietype = "gif", outdir = getwd())
>>>
ようするにユーザはプロットを何枚か出力するような関数を用意するだけ。それをsaveMovie()関数に渡せば、内部でImageMagickやらSWF Toolsががんばってくれてアニメになる、という仕組みになっている。ImageMagickを MacPortsあたりで入れておく必要がある。
追記: あやまって hatena_bookmark_nocomment.rb を有効にしていたせいで、はてブのコメント一覧が見えなくなってしまいました。このプラグインを無効にしたのですが、最初にブクマされた時点で非表示に設定されると公開することができないようです。コメント頂いたかた、ごめんなさい。見れないので直接、ここにコメントして頂けると幸いです><
コメント
コメントを投稿