読者です 読者をやめる 読者になる 読者になる

廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

RStan でベイズ逆問題もどき

y_i = p_i A +\varepsilon_i,  0\le p_i \le 1 という系列を考える。A は定数。

(この系列になんらかの解釈があればおもしろかったんだけど思いつかなかった。)

いま、観測されるのは y_i のみである。

f:id:abrahamcow:20151014032900p:plain

p_i が既知ならば A を推定するのはかんたん。A が既知ならば p_i を推定するのはかんたん。

しかし、いま観測されるのは y_i のみ。

A \times p_iy_i になる組み合わせは無数にあり、p_iA を同時に推定するのは困難に思える。

しかし、以下の仮定を置くことで、p_iA を推定することができる。

推定に用いた Stan のコードはこう。

#inverse.stan
data{
	int<lower=0> N;
	real <lower=0>y[N];
}
parameters{
	real <lower=0,upper=1> p[N];
	real <lower=0>sigma;
	real <lower=0>A;
}
transformed parameters  {
    real mu1[N];
   for (i in 1:N) {
    	mu1[i] <- p[i]*A;
   	}
}
model{
	for (i in 1:N) {
 	y[i] ~ normal(mu1[i], sigma);
	}
}

R のコードはこう。

##データ生成##
set.seed(123)
p <- runif(20)
A <- 100
y <- rnorm(20,p*A)
plot(y,type="b")
dat4stan <- list(N=20,y=y)
############
model1 <- stan_model("inverse.stan")
fit <- sampling(model1, data=dat4stan)
traceplot(fit,pars=c("A","sigma"))
ex<-extract(fit)
##図示##
plot(y,type="b")
lines(apply(ex$mu,2,mean),col="orange")
plot(p,apply(ex$p,2,mean),xlim=c(0,1),ylim=c(0,1),ylab="estimates")
abline(0,1,lty=2)

Aσ について、MCMC サンプリングのトレースラインを描いてみる。収束してくれたようだ。

f:id:abrahamcow:20151014032924p:plain

次に、原系列(黒)と推定値(オレンジ)を重ねてみる。

f:id:abrahamcow:20151014033109p:plain

最後に、横軸に乱数で生成した真の p_i、縦軸に p_i の推定値をとりプロットした。

f:id:abrahamcow:20151014033755p:plain

ほぼ一直線上に並んでいることから、p_i の推定値の妥当さが分かる。

しかしなんでこれが求まってしまうんだろう……?