MC Approximation
直观上理解，已知后验分布，但是该分布很难进行积分，而我们想计算均值、概率、分位数等。这时我们可以采用 MC 近似。也就是从后验分布中产生规模为
的样本，
- 用样本的均值代替总体均值
- 用样本的经验分布代替该分布
- 用样本的分位数近似总体的分位数
总结起来，就是先从后验分布中产生样本，然后用样本的量近似总体的量。
## ########################################
## prior distribution: gamma(a, b)
## Y_1, ..., Y_n\mid \theta: iid Poisson(\theta)
## posterior distribution: gamma(a + \sum y_i, b+n)
## ########################################
## ########################################
## Expectation
## posterior mean: (a+\sum y_i)/(b+n)=1.51
## ########################################
a = 2; b = 1
sy = 66; n = 44
theta.mc10 = rgamma(10, a+sy, b+n)
theta.mc100 = rgamma(100, a+sy, b+n)
theta.mc1000 = rgamma(1000, a+sy, b+n)
mean(theta.mc10)
mean(theta.mc100)
mean(theta.mc1000)
## ########################################
## Probabilities
##
## ########################################
## posterior Probabilities
pgamma(1.75, a+sy, b+n)
## MC approximations
mean(theta.mc10 < 1.75)
mean(theta.mc100 < 1.75)
mean(theta.mc1000 < 1.75)
## ########################################
## quantiles
##
## ########################################
## posterior quantiles
qgamma(c(.025, .975), a+sy, b+n)
## MC approximations
quantiles(theta.mc10, c(.025, .975))
quantiles(theta.mc100, c(.025, .975))
quantiles(theta.mc1000, c(.025, .975))
## #######################################
## Log-odds
## #######################################
a = 1; b = 1
theta.prior.mc = rbeta(10000, a, b)
gamma.prior.mc = log(theta.prior.mc/(1-theta.prior.mc))
n0 = 860-441; n1 = 441
theta.post.mc = rbeta(10000, a+n1, b+n0)
gamma.post.mc = log(theta.post.mc/(1-theta.post.mc))
## #######################################
## Functions of two parameters
## #######################################
a = 2; b = 1
sy1 = 217; n1 = 111
sy2 = 66; n2 = 44
theta1.mc = rgamma(10000, a+sy1, b+n1)
theta2.mc = rgamma(10000, a+sy2, b+n2)
mean(theta1.mc > theta2.mc)
## ######################################
## posterior
##
## ######################################
a = 2; b = 1
sy1 = 217; n1 = 111
sy2 = 66; n2 = 44
theta1.mc = rgamma(10000, a+sy1, b+n1)
theta2.mc = rgamma(10000, a+sy2, b+n2)
y1.mc = rpois(10000, theta1.mc)
y2.mc = rpois(10000, theta2.mc)
mean(y1.mc > y2.mc)
## #####################################
## ratio
##
## #####################################
a = 1; b = 2
t.mc = NULL
for (s in 1:10000){
theta1 = rgamma(1, a+sy1, b+n1)
y1.mc = rpois(n1, theta1)
t.mc = c(t.mc, sum(y1.mc==2)/sum(y1.mc==1))
}
为了计算积分 我们通常采用MC模拟:
- 1.计算区域的volume：
- 2.近似：
根据大数律有
并且由中心极限定理有
其中
## estimate pi
##
## I = \int H(x, y)dxdy
## where H(x, y) = 1 when x^2+y^2 <= 1;
## otherwise H(x, y) = 0
## volume
V = 4
n = 100000
x = runif(n, -1, 1)
y = runif(n, -1, 1)
H = x^2+y^2
H[H<=1] = 1
H[H>1] = 0
I = V* mean(H)
cat("I = ", I, "\n")
## n = 100, I = 2.96
## n = 1000, I = 3.22
## n = 10000, I = 3.1536
## n = 100000, I = 3.14504
Last modified 4yr ago