gamma分布的density的奇怪特性,如下:
Poisson的Gamma先验
h(x) 的置信区间 的 获取
> n =
> sumx=
>
> alpha=
> beta=0.01
> pmean=(alpha+sumx)/(beta+n)
> L=qgamma(0.025, alpha+sumx, beta+n) // 获得cdf的边界
> U=qgamma(0.975, alpha+sumx, beta+n)
>
> cat("Posterior mean: ", pmean, " (", L, ",",U,")")
Posterior mean: 382.8796 ( 378.1376 , 387.6506 )
// 95% 置信区间的边界值 还有期望。
Monte Carlo sampling:
N= # or 500 or 5000
L=U=rep(NA,length=)
for (i in :) {
dat=sort(rgamma(N,alpha+sumx,beta+n))
L[i]=dat[0.025*N]
U[i]=dat[0.975*N]
} // 获得某分位点的大量样本
widthL=max(L)-min(L)
widthU=max(U)-min(U)
par(mfrow=c(,))
hist(L,probability=T,xlab="Lower 95% CI bound")
points(qgamma(0.025,alpha+sumx,beta+n),,pch=,col=)
hist(U,probability=T,xlab="Upper 95% CI bound")
points(qgamma(0.975,alpha+sumx,beta+n),,pch=,col=)
cat("L interval variability (range):",widthL,"\n")
cat("U interval variability (range):",widthU,"\n")
Sampling估计的分位点,看来与True value差不多呢。
> mean(L)
[] 378.0747
> mean(U)
[] 387.5853
问题来了,N要多大才能保证要求的分位点估计精度:Sol 要 according to Central Limit Thearem.
p(y) 的预测
alpha=; beta=0.01
sumx=; n=65 theta=rgamma(,alpha+sumx,beta+n) // 后验sita
y=rpois(,theta) // poisson分布sampling,带入后验sita的f(y|sita),5000个相应的预测值
hist(y,probability=T,ylab="Density",main="Posterior predictive distribution")
// sampling法得到直方图。
// 相当吻合!
// 精确函数得到散点图。
xx=:
pr=dnbinom(xx,sumx+alpha,-/(beta+n+))
#lines(xx,pr,col=) # The (incorrect) continuous version - ok as an approx.
#The (correct) discrete version:
for (i in :length(xx)) {
lines(c(xx[i],xx[i+]),rep(pr[i],),col=,lwd=)
}