# MC Approximation

## For Bayesian

直观上理解，已知后验分布，但是该分布很难进行积分，而我们想计算均值、概率、分位数等。这时我们可以采用 MC 近似。也就是从后验分布中产生规模为 $$n$$ 的样本，

* 用样本的均值代替总体均值
* 用样本的经验分布代替该分布
* 用样本的分位数近似总体的分位数

总结起来，就是先从后验分布中产生样本，然后用样本的量近似总体的量。

```r
## ########################################
## 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))
    }
```

## For Integration

为了计算积分 ![](https://latex.codecogs.com/gif.latex?I%20%3D%20\int%20_D%20g\(\mathbf%20x\)d\mathbf%20x) 我们通常采用MC模拟:

1. 计算区域的volume：![](https://latex.codecogs.com/gif.latex?V%20%3D%20\int_D%20d\mathbf%20x)
2. 近似：![](https://latex.codecogs.com/gif.latex?\hat%20I_m%3DV\frac{1}{m}\sum\limits_{i%3D1}^mg\(\mathbf%20x^{\(m\)}\))

根据大数律有

![](https://latex.codecogs.com/gif.latex?\lim_{m\rightarrow%20\infty}%20\hat%20I_m%3DI)

并且由中心极限定理有

![](https://latex.codecogs.com/gif.latex?\frac{1}{V}{}\sqrt{m}\(\hat%20I_m-I\)\rightarrow%20N\(0%2C%20\sigma^2\))

其中![](https://latex.codecogs.com/gif.latex?\sigma^2%3Dvar\(g\(\mathbf%20x\)\))

```r
## 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
```


---

# Agent Instructions: Querying This Documentation

If you need additional information that is not directly available in this page, you can query the documentation dynamically by asking a question.

Perform an HTTP GET request on the current page URL with the `ask` query parameter:

```
GET https://mc.hohoweiya.xyz/mcapproximation.md?ask=<question>
```

The question should be specific, self-contained, and written in natural language.
The response will contain a direct answer to the question and relevant excerpts and sources from the documentation.

Use this mechanism when the answer is not explicitly present in the current page, you need clarification or additional context, or you want to retrieve related documentation sections.