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