既然估歐國盃阿 Paul 係六估四中，假定現在是世界盃之前，我們再要阿 Paul 估八場世盃，我們的 Prior 應該是 p=4/6 信心最高，但是 p=4/6 周圍的數值信心又怎樣呢？換句話說，我們想要的是 Beta Prior 的 mean 在 4/6 ，但是 Variance 又怎樣？
查查維基．會見到 Mean 及 Variance 換算 alpha 及 beta 的公式。將此公式轉為 R 的 function 就會是：
# chew that < -
alpha <- function (mu,v) {mu*(((mu*(1-mu))/v)-1)}
beta <- function (mu,v) {(1-mu)*(((mu*(1-mu))/v)-1)}
v 的大細，其實代表我們有幾企硬 p=4/6 。
# chew that < -
p <- seq(0.00,1,by=0.01)
prior001 <- dbeta(p,alpha(4/6,0.001),beta(4/6,0.001))
prior002 <- dbeta(p,alpha(4/6,0.002),beta(4/6,0.002))
prior003 <- dbeta(p,alpha(4/6,0.003),beta(4/6,0.003))
prior03 <- dbeta(p,alpha(4/6,0.03),beta(4/6,0.03))
plot(p,prior001, type="l", ylab="Density", lty=2)
lines(p,prior002, lty=1)
lines(p,prior003, lty=3)
lines(p,prior03, lty=4)
legend(0,12,c("v=0.001","v=0.002","v=0.003","v=0.03"),lty=c(2,1,3,4))
Variance 愈細企得愈硬。
當阿 Paul 估完8場世界盃中了八場，我們可以計算不同 Variance 的 prior 所計出之 Posterior 。得出來的數字都相當細。
# chew that < -
variance.beta <- c(0.001,0.002,0.003,0.03)
round( 1- pbeta(0.5, alpha(4/6,variance.beta)+8, beta(4/6,variance.beta)+0, lower.tail=FALSE),3)
假設我自己事前不算太企硬 p=4/6 ，我用 v 值 0.03 。事前和事後的 distribution 就是這樣求出：
# Chew that < -
p <- seq(0.00,1,by=0.01)
prior <- dbeta(p,alpha(4/6,0.03),beta(4/6,0.03))
post.d <- dbeta(p,alpha(4/6,0.03)+8,beta(4/6,0.03)+0)
plot(p,post.d, type="l", ylab="Density", lty=2)
lines(p,prior, lty=1)
還可以計算 95% interval estimation 。
# Chew that < -
qbeta(c(0.05,0.95),alpha(4/6,0.03)+8,beta(4/6,0.03)+0)
得出的 p 的 95% Credible Interval 是 0.68 至 0.97 。
另一個方法是用 simulation 的方法：
# Chew that < -
ps <- rbeta (1000,alpha(4/6,0.03)+8,beta(4/6,0.03)+0)
mean(ps)
quantile(ps,c(0.05,0.95))
可求出 95% Credible interval 為 0.69 至 0.97 。更可求出 mean 值是 0.85 。