[Bayes] qgamma & rgamma: Central Credible Interval

gamma分布的density的奇怪特性，如下：

---

 

 

Poisson的Gamma先验

 

 h(x) 的置信区间 的 获取

> n = 65
> sumx=24890
> 
> alpha=1
> 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=500 # or 500 or 5000
L=U=rep(NA,length=250)

for (i in 1:250) {
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(1,2))
hist(L,probability=T,xlab="Lower 95% CI bound")
points(qgamma(0.025,alpha+sumx,beta+n),0,pch=16,col=2)

hist(U,probability=T,xlab="Upper 95% CI bound")
points(qgamma(0.975,alpha+sumx,beta+n),0,pch=16,col=2)

cat("L interval variability (range):",widthL,"\n")
cat("U interval variability (range):",widthU,"\n")

 

Sampling估计的分位点，看来与True value差不多呢。

> mean(L)
[1] 378.0747
> mean(U)
[1] 387.5853

 

问题来了，N要多大才能保证要求的分位点估计精度：Sol 要 according to Central Limit Thearem.

 

 p(y) 的预测

alpha=1; beta=0.01
sumx=24890; n=65

theta=rgamma(5000,alpha+sumx,beta+n)　　// 后验sita
y=rpois(5000,theta)　　// poisson分布sampling，带入后验sita的f(y|sita),5000个相应的预测值
hist(y,probability=T,ylab="Density",main="Posterior predictive distribution")
// sampling法得到直方图。
// 相当吻合！
// 精确函数得到散点图。
xx=300:500
pr=dnbinom(xx,sumx+alpha,1-1/(beta+n+1))
#lines(xx,pr,col=2) # The (incorrect) continuous version - ok as an approx.
#The (correct) discrete version:
for (i in 1:length(xx)) {
　　lines(c(xx[i],xx[i+1]),rep(pr[i],2),col=2,lwd=2)
}