# 阿 Paul #2

Posted on Jul 13, 2010 by Chung-hong Chan

# 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 愈細企得愈硬。

# 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)


# 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)


# Chew that < -
qbeta(c(0.05,0.95),alpha(4/6,0.03)+8,beta(4/6,0.03)+0)


# 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))