てさぐれ!!RStanもの

MikuHatsune2014-07-13

なんかよくわからんけどStanやってみた。導入はこちら
分からないことはだいたいここにある。

答えのあるモデルとして、y=sin(\beta x)があるとして、周期\betaと誤差N(0, \sigma ^2)を推定する感じでやろう。
よく分からんけどdata, parameters, model 部分をそれなりに書く。
文末に ; を忘れたらコンパイル時にエラーが出てイラッ☆となる。

# 元のモデル
x <- seq(0, 4*pi, length=50)
beta <- 0.6
y <- sin(beta * x) + rnorm(length(x), 0, 0.2)
plot(x, y)
lines(x, sin(beta * x))

# Stan モデルの構築
library(rstan)
library(coda)
code1 <- '
	data {
		int<lower = 1> N;
		real y[N];
		real x[N];
	}
	parameters {
		real < lower = 0 > b0;
		real < lower = 0 > noise;
	}
	transformed parameters {
		real theta[N];
		for (j in 1:N){
			theta[j] <- sin(b0 * x[j]);
		}
	}
	model {
		y ~ normal(theta, noise);
	}
'

# 入力データ
p1 <- list(N=length(x), x=x, y=y)
# 推定。chain を莫大に増やした。
fit <- stan(model_code = code1, data = p1, iter = 1000, chains = 100)

 
Stanによる推定したパラメータを取り出すにはextract関数を使う。結構いい感じに推定できるように見えて実は裏があって…

d.ext <- extract(fit, permuted=TRUE)
plot(x, y)
lines(x, sin(beta * x))
lines(x, sin(median(d.ext$b0)*x), col=2)
legend("bottomleft", legend=c("True", "Estimate"), col=1:2, bty="n", lty=1, cex=1.5, lwd=3)


 
\beta\sigmaについて各chainでの推定結果を示す。たいていのブログは4 chainでほぼ完璧に推定されているいい感じの結果しか載せてくれていないことが多いが、各chainで収束がばらつくとこんな感じになる。

fit.coda <- mcmc.list(lapply(1:ncol(fit),function(x) mcmc(as.array(fit)[,x,])))
# plot(fit)
i <- 1
par_name <- fit@sim$fnames_oi[i]
tmp <- mapply(function(x) x[, i], fit.coda)
par(mfrow=c(1, 2))
matplot(tmp, type="l")
plot(density(tmp), main=par_name)



 
というのも初期値に依存するわけで、初期値は推定して得られたStanfitオブジェクト(S4クラスに注意)のinitsに入っている。

fit@inits
$b0
[1] 1.39958

$noise
[1] 0.666027

$theta
 [1]  0.000000000  0.351273875  0.657776314  0.880442204  0.990891776
 [6]  0.975047734  0.834929474  0.588395697  0.266868243 -0.088672770
[11] -0.432912039 -0.721974721 -0.919018479 -0.998929198 -0.951521899
[16] -0.782838854 -0.514379477 -0.180360122  0.176646942  0.511139558
[21]  0.780485139  0.950354379  0.999096680  0.920499615  0.724580735
[26]  0.436310781  0.092431056 -0.263229423 -0.585340127 -0.832846601
[31] -0.974203029 -0.991392900 -0.882225288 -0.660614095 -0.354804664
[36] -0.003773784  0.347738083  0.654929166  0.878646582  0.990376540
[41]  0.975878554  0.837000457  0.591442887  0.270503262 -0.084913220
[46] -0.429507131 -0.719358426 -0.917524255 -0.998747491 -0.952675868

 
ぶっちゃけて言うと、1つのパラメータの正弦関数を推定するくらいだったらoptimを使うのが楽だろう。青木先生のサイトからパクると、これも初期値によって最終的な推定値が変わる。真の値に近いところの初期値でやれば推定結果もかなりそれっぽくなるし、Stanもそのようである。結局のところ局所解に落ち込むとこういうことになるので、推定結果は眺めて確かめたほうがいいっぽい。

resid <- function(par){
	yhat <- sin(par[1]*x)
	sum((y-yhat)^2)
}
# やれる推定法をとりあえずやっておく
o.m <- c("Nelder-Mead", "BFGS", "CG", "L-BFGS-B", "SANN")
s0 <- seq(0, 1, length=100)
s1 <- mapply(function(y) mapply(function(x) optim(x, resid, method=y)$par, s0), o.m)
init_b <- mapply(function(x) fit@inits[[x]][[1]], seq(fit@inits)) # Stanの初期値
init_b_median <- apply(tmp, 2, median) # 各chain初期値から最終的に最頻値を取って最終的な推定値とする

matplot(s0, s1, xlab="initial parameter", ylab="estimate", pch=16, ylim=c(-1, 1)*2, col=seq(o.m)+1)
points(init_b, init_b_median, col=1, pch=1)
legend("bottomright", legend=c("Stan", o.m), bty="n", pch=c(1, rep(16, length(o.m))), col=seq(length(o.m)+1), cex=1) 


 
Stanは柔軟なモデリングが〜とよく言われるけど、自分勝手に柔軟にモデルを構築するよりかは、世間一般的に認知されているものを使うほうが簡単だし受け入れられやすいっぽい。こんな論文もあったけど、はっきり言って病院の循環器内科でこれを理解できる人は0.001%もいないと思う。