[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)
}