# Monte Carlo Integration

#### 2018-09-23

The SI package provides several methods of MC Integrating including

• Stochastic Point Method
• Mean Value Method
• Important Sampling Method
• Stratified Sampling Method

Note that the Stochastic Point Method is only a stochastic point method to estimate integration. However, it is provided to be easily compared with other MC methods. These four methods are all belong to stochastic methods.

## Stochastic Integration

Suppose that $$h(x)\in \mathbb{C}^1$$ is well-defined and bounded on finite interval $$C=[a,b]$$. W.L.O.G., $$0\leq h(x)\leq M$$. We are tried to calculate the integral $I=\int_a^bh(x)\mathrm{d}x$ That is to say, we need to calculate the area under $$h(x)$$: $$D=\{(x,y):0\leq y\leq h(x),\ x\in C\}$$.

## Stochastic Point Method

Uniformly sampling from $$G=C\times(0,M)$$ for $$N$$ times and get the stochastic points $$Z_1,\cdots,Z_N$$ where $$Z_i=(X_i,Y_i),\ i=1,2,\cdots,N$$.

For $$i=1,2,\cdots,N$$, let $\xi_i=\begin{cases} 1&,Z_i\in D\\ 0&,\text{otherwise} \end{cases}$

Then $$\{\xi_i\}$$ are outcomes of independent duplicate trials and $$\xi_1,\cdots,\xi_N\overset{i.i.d.}{\sim}Bernoulli(1,p)$$ where $p=\mathbb{P}(Z_i\in D)=\dfrac{I}{M(b-a)}$

Thus, the stochastic point method is given by $\hat{I}_1=\hat{p}(b-a)$ where $$\hat{p}$$ is the estimator of $$p$$ and given by the proportion of points $$\{Z_i\}$$ that lie under the curve $$h(x)$$.

From Strong Law of Large Number, \begin{align*} \hat{p}=\dfrac{\sum\limits_{i=1}^N\xi_i}{N}\xrightarrow{a.s.}p\\ \hat{I}=\hat{p}M(b-a)\xrightarrow{a.s.}I \end{align*} which means that we can use $$\hat{I}$$ to approximate $$I$$ better if we try more times (with large $$N$$).

To estimate the variance (or precision), we use the Central Limit Theorem and get \begin{align*} \dfrac{\hat{p}-p}{Var(p)}&=\dfrac{\hat{p}-p}{\sqrt{\dfrac{p(1-p)}{N}}}\xrightarrow{D}N(0,1)\\ \dfrac{\hat{I}_1-I}{\sqrt{\frac{1}{N}}}&\xrightarrow{D}N(0,[M(b-a)]^2p(1-p)) \end{align*}

Therefore, the variance can be estimated by $Var(\hat{I}_1)\approx \dfrac{1}{N}[M(b-a)]^2\hat{p}(1-\hat{p})$

For $$m$$-dimension function on $$[a_1,b_1]\times \cdots \times[a_m,b_m]$$, we also have \begin{align*} \hat{I}_1&=\hat{p}\prod\limits_{i=1}^m(b_i-a_i)\\ Var(\hat{I}_1)&\approx \dfrac{1}{N}\left[M\prod\limits_{i=1}^m(b_i-a_i)\right]^2\hat{p}(1-\hat{p}) \end{align*}

In SI package, use the following code to carry out stochastic point method.

## To integrate exp(x) from -1 to 1
set.seed(0)
h <- function(x){
exp(x)
}
N <- 100000
SPMresult <- SI.SPM(h,-1,1,exp(1),N)
I1 <- SPMresult[]
VarI1 <- SPMresult[]

## Mean Value Method

Let $$U\sim Uniform(a,b)$$, then \begin{align*} \mathbb{E}[h(U)]&=\int_a^bh(x)\dfrac{1}{b-a}\mathrm{d}x\\ &=\dfrac{I}{b-a}\\ I&=(b-a)\mathbb{E}[h(U)] \end{align*}

Suppose that we first sample $$U_1,\cdots,U_N\overset{i.i.d.}{\sim}Uniform(a,b)$$, then $$Y_i=h(U_i),\ i=1,2,\cdots,N$$ are i.i.d random variables. From Strong Law of Large Number, $\overline{Y}=\dfrac{1}{N}\sum\limits_{i=1}^Nh(U_i)\xrightarrow{a.s.}\mathbb{E}[h(U)]=\dfrac{I}{b-a}$

Thus we get the mean value estimator $\hat{I}_2=\dfrac{b-a}{N}\sum\limits_{i=1}^Nh(U_i)$

To estimate the variance (or precision), we use the Central Limit Theorem and get \begin{align*} \dfrac{\hat{I}_2-I}{\sqrt{\frac{1}{N}}}&\xrightarrow{D}N(0,(b-a)^2Var[h(U)])\\ \hat{I}_2&\xrightarrow{D}N(I,\dfrac{(b-a)^2}{N}Var[h(U)]) \end{align*} where $Var[h(U)]=\int_a^b\{h(u)-\mathbb{E}[h(U)]\}^2\dfrac{1}{b-a}\mathrm{d}u$

Since it can be estiamted by sample variance $Var[h(U)]\approx \dfrac{1}{N}\sum\limits_{i=1}^N(Y_i-\overline{Y})^2$ we can estimate the variance of $$\hat{I}_2$$ by $Var(\hat{I}_2)\approx \dfrac{(b-a)^2}{N^2}\sum\limits_{i=1}^N(Y_i-\overline{Y})^2$

For $$m$$-dimension function on $$[a_1,b_1]\times \cdots \times[a_m,b_m]$$, we also have \begin{align*} \hat{I}_2&=\dfrac{1}{N}\prod\limits_{i=1}^m(b_i-a_i)\sum\limits_{i=1}^Nh(U_i)\\ Var(\hat{I}_2)&\approx \dfrac{1}{N}\left[\prod\limits_{i=1}^m(b_i-a_i)\right]^2\sum\limits_{i=1}^N(Y_i-\overline{Y})^2 \end{align*}

In SI package, use the following code to carry out mean value method.

## To integrate exp(x) from -1 to 1
set.seed(0)
h <- function(x){
exp(x)
}
N <- 100000
MVMresult <- SI.MVM(h,-1,1,N)
I2 <- MVMresult[]
VarI2 <- MVMresult[]

## Important Sampling Method

Suppose $$g(x)$$ is a density having similar shape with $$|h(x)|$$, $$h(x)=0$$ when $$g(x)=0$$ and $$h(x)=o(g(x))$$ as $$x\rightarrow\infty$$.

Suppose that $$X_i\overset{i.i.d.}{\sim}g(x),\ i=1,2,\cdots,N$$ and let $\eta_i=\dfrac{h(X_i)}{g(X_i)}$ then $\mathrm{E}\eta_i=\int_C\dfrac{h(x)}{g(x)}g(x)\mathrm{d}x=I$

Therefore, the important sampling estimator is given by $\hat{I}_3=\overline{\eta}=\dfrac{1}{N}\sum\limits_{i=1}^N\dfrac{h(X_i)}{g(X_i)}$

The variance can be estimated by \begin{align*} Var(\hat{I}_3)&=Var(\overline{\eta})\\ &=\dfrac{1}{N}Var\left(\dfrac{h(X)}{g(X)}\right)\\ &\approx \dfrac{1}{N^2}\sum\limits_{i=1}^N \left[\dfrac{h(X_i)}{g(X_i)}-\hat{I}_3\right]^2 \end{align*}

In SI package, use the following code to carry out mean value method.

## To integrate exp(x) from -1 to 1
## Use the sampling density (3/2+x)/3
set.seed(0)
h <- function(x){
exp(x)
}
N <- 100000
g <- function(x){return((3/2+x)/3)}
G_inv <- function(y){return(sqrt(6*y+1/4)-3/2)}
ISMresult <- SI.ISM(h,g,G_inv,N)
I3 <- ISMresult[]
VarI3 <- ISMresult[]

## Stratified Sampling Method

Suppose that $$C=\bigcup\limits_{j=1}^m C_j$$ and $$C_i\cap C_j=\varnothing$$. For every $$C_j$$, use mean value method to approximate $\hat{I}_{C_j}=\int_{C_j}h(x)\mathrm{d}x$

And get $\hat{I}_4=\sum\limits_{j=1}^m\hat{I}_{C_j}$

It can be shown that the varaince of $$\hat{I}_4$$ is smaller than $$\hat{I}_2$$.

In SI package, use the following code to carry out mean value method.

## To integrate exp(x) from -1 to 1
set.seed(0)
h <- function(x){
exp(x)
}
N <- 100000
SSMresult <- SI.SSM(h,-1,1,10,N)
I4 <- SSMresult[]
VarI4 <- SSMresult[]