二項分布のMCMC

ある産地のミカンを3個食べました。3個とも「甘い」ミカンでした。改めてこの産地のミカンをひょいと１つとって食べたとき、このミカンが「甘い」確率は？ただし、ミカンは「甘い」「酸っぱい」の２種類しかないとし、その確率は、それぞれ $\theta,1-\theta$。

共役分布を使った計算

$x_i$は、$i$番目のミカンが甘い($x_i=1$)か酸っぱいか($x_i=0$)の離散変数。
3個分のデータ得られた$\theta$の事後確率は、
$$
\Pr(\theta|x_1,x_2,x_3)=\frac{\Pr(x_1,x_2,x_3|\theta)\Pr(\theta)}{\Pr(x_1,x_2,x_3)}
$$
このうち、尤度は、
$$
\Pr(x_1,x_2,x_3|\theta)={}_3C_3\theta^3(1-\theta)^{3-3}=\theta^3
$$
の二項分布となる。

事前分布$\Pr(\theta)$に一様分布$\Pr(\theta)=1$を仮定すると、事後分布は、
$$
\Pr(\theta|x_1,x_2,x_3)=k\theta^3\propto \theta^3
$$
であり、$Be(4,1)$のベータ分布となる。

curve(dunif,0,1,ylim=c(0,4))
curve(4*x^3,0,1,add=T,col=2)
legend(0,4,c("事前分布","事後分布"),lty=c(1,1),col=c(1,2),cex=0.8)

事前分布$\Pr(\theta)$にベータ分布$Be(2,2)$を仮定すると、
$$
\Pr(\theta)=k\theta^{2-1}(1-\theta)^{2-1}
$$
であり、事後分布は、
$$
\Pr(\theta|x_1,x_2,x_3)
\propto \theta^3\times k\theta^{2-1}(1-\theta)^{2-1}
\propto \theta^{5-1}(1-\theta)^{2-1}
$$
であり、$Be(5,2)$のベータ分布となる。

f01<-function(x) dbeta(x,2,2)
f02<-function(x) 30*x^4*(1-x)
curve(f01,0,1,ylim=c(0,4))
curve(f02,0,1,add=T,col=2)
legend(0,4,c("事前分布","事後分布"),lty=c(1,1),col=c(1,2),cex=0.8)

MCMCによる計算

まずは事前分布を一様分布で

RMarkdownで以下を記述。
model{}にtheta~uniform(0,1)と事前分布を記述してもよいが、Stanはデフォルトで、一様分布なので省略できる。

```{stan,output.var="mikan01"}
data{
  int N;
  int y;
}

parameters{
  real<lower=0,upper=1> theta;
}

model{
  y~binomial(N,theta);
}
```

データを作成!

d01<-list(N=3,y=3)

MCMC実行！

fit01<-sampling(mikan01,d01)

結果のサマリー

summary(fit01)$summary

            mean     se_mean        sd      2.5%        25%        50%        75%
theta  0.7945553 0.004674643 0.1689501  0.389770  0.6963226  0.8373022  0.9323246
lp__  -3.0953976 0.023468577 0.7917435 -5.441765 -3.3238426 -2.7864350 -2.5677014
           97.5%    n_eff    Rhat
theta  0.9927365 1306.231 1.00069
lp__  -2.5024634 1138.140 1.00602

事後分布

stan_dens(fit01,pars="theta")

事前分布をベータ分布$Beta(2,2)$で

```{stan,output.var="mikan02"}
 data{
   int N;
   int y;
 }

 parameters{
   real<lower=0,upper=1> theta;
 }

 model{
   y~binomial(N,theta);
   theta~beta(2,2);
 }
```

MCMC実行！

fit02<-sampling(mikan02,d01)

結果のサマリー

summary(fit02)$summary

            mean     se_mean        sd       2.5%        25%        50%        75%
theta  0.7120068 0.004321131 0.1600423  0.3420686  0.6162984  0.7319016  0.8357492
lp__  -4.7194705 0.019542178 0.7538810 -6.8582568 -4.8917791 -4.4233716 -4.2449272
           97.5%    n_eff     Rhat
theta  0.9480958 1371.749 1.001601
lp__  -4.1883397 1488.194 1.000210

事後分布

stan_dens(fit02,pars="theta")