MC Approximation
For Bayesian
直观上理解，已知后验分布，但是该分布很难进行积分，而我们想计算均值、概率、分位数等。这时我们可以采用 MC 近似。也就是从后验分布中产生规模为
的样本，
- 用样本的均值代替总体均值
- 用样本的经验分布代替该分布
- 用样本的分位数近似总体的分位数
总结起来，就是先从后验分布中产生样本，然后用样本的量近似总体的量。
1
## ########################################
2
## prior distribution: gamma(a, b)
3
## Y_1, ..., Y_n\mid \theta: iid Poisson(\theta)
4
## posterior distribution: gamma(a + \sum y_i, b+n)
5
## ########################################
6
7
## ########################################
8
## Expectation
9
## posterior mean: (a+\sum y_i)/(b+n)=1.51
10
## ########################################
11
12
a = 2; b = 1
13
sy = 66; n = 44
14
15
theta.mc10 = rgamma(10, a+sy, b+n)
16
theta.mc100 = rgamma(100, a+sy, b+n)
17
theta.mc1000 = rgamma(1000, a+sy, b+n)
18
19
mean(theta.mc10)
20
mean(theta.mc100)
21
mean(theta.mc1000)
22
23
## ########################################
24
## Probabilities
25
##
26
## ########################################
27
28
## posterior Probabilities
29
30
pgamma(1.75, a+sy, b+n)
31
32
## MC approximations
33
mean(theta.mc10 < 1.75)
34
mean(theta.mc100 < 1.75)
35
mean(theta.mc1000 < 1.75)
36
37
## ########################################
38
## quantiles
39
##
40
## ########################################
41
42
## posterior quantiles
43
qgamma(c(.025, .975), a+sy, b+n)
44
45
## MC approximations
46
47
quantiles(theta.mc10, c(.025, .975))
48
quantiles(theta.mc100, c(.025, .975))
49
quantiles(theta.mc1000, c(.025, .975))
50
51
52
## #######################################
53
## Log-odds
54
## #######################################
55
56
a = 1; b = 1
57
theta.prior.mc = rbeta(10000, a, b)
58
gamma.prior.mc = log(theta.prior.mc/(1-theta.prior.mc))
59
60
n0 = 860-441; n1 = 441
61
theta.post.mc = rbeta(10000, a+n1, b+n0)
62
gamma.post.mc = log(theta.post.mc/(1-theta.post.mc))
63
64
## #######################################
65
## Functions of two parameters
66
## #######################################
67
68
a = 2; b = 1
69
sy1 = 217; n1 = 111
70
sy2 = 66; n2 = 44
71
theta1.mc = rgamma(10000, a+sy1, b+n1)
72
theta2.mc = rgamma(10000, a+sy2, b+n2)
73
74
mean(theta1.mc > theta2.mc)
75
76
## ######################################
77
## posterior
78
##
79
## ######################################
80
81
a = 2; b = 1
82
sy1 = 217; n1 = 111
83
sy2 = 66; n2 = 44
84
85
theta1.mc = rgamma(10000, a+sy1, b+n1)
86
theta2.mc = rgamma(10000, a+sy2, b+n2)
87
y1.mc = rpois(10000, theta1.mc)
88
y2.mc = rpois(10000, theta2.mc)
89
90
mean(y1.mc > y2.mc)
91
92
## #####################################
93
## ratio
94
##
95
## #####################################
96
a = 1; b = 2
97
t.mc = NULL
98
99
for (s in 1:10000){
100
theta1 = rgamma(1, a+sy1, b+n1)
101
y1.mc = rpois(n1, theta1)
102
t.mc = c(t.mc, sum(y1.mc==2)/sum(y1.mc==1))
103
}
Copied!
For Integration
为了计算积分 d\mathbf%20x)
我们通常采用MC模拟:
- 1.计算区域的volume：
- 2.近似：
根据大数律有
并且由中心极限定理有
其中))
1
## estimate pi
2
##
3
## I = \int H(x, y)dxdy
4
## where H(x, y) = 1 when x^2+y^2 <= 1;
5
## otherwise H(x, y) = 0
6
7
## volume
8
V = 4
9
10
n = 100000
11
x = runif(n, -1, 1)
12
y = runif(n, -1, 1)
13
H = x^2+y^2
14
H[H<=1] = 1
15
H[H>1] = 0
16
I = V* mean(H)
17
cat("I = ", I, "\n")
18
19
## n = 100, I = 2.96
20
## n = 1000, I = 3.22
21
## n = 10000, I = 3.1536
22
## n = 100000, I = 3.14504
Copied!
Last modified 3yr ago
Copy link