[Math] Integration using Monte Carlo Method

monte carlosimulation

I'm trying to solve this integral using the Monte Carlo Method.

$$I=\int_0^\pi \frac{1}{\sqrt{2\pi}}e^\frac{-sin(x)^2}{2}dx$$

Now it seems to me that there is a normal probability density function in there, but I'm not sure because of the sine function.

If there is a normal, then it's easy to simulate $N$ random numbers from the standard normal distribution and compute $I$. It's is also easy to find the size of the sample to get a 95% confidence interval.

So is there a normal probability density function? Or from what distribution should I get the random numbers to compute $I$?.

Best Answer

Recall that if $Y$ is a random variable with density $g_Y$ and $h$ is a bounded measurable function, then $$\mathbb E[h(Y)] = \int_{\mathbb R} h(y)g_Y(y)\,\mathsf dy. $$ Moreover, if $Y\sim\mathcal U(0,1)$, then $a+(b-a)U\sim\mathcal U(a,b)$. So applying the change of variables $x=a+(b-a)u$ (with $a=0$, $b=\pi$) to the given integral, we have $$I = \int_0^1 \frac{\pi}{\sqrt{2\pi}} e^{-\frac12\sin^2 (\pi u) }\,\mathsf du=\int_0^1 h(u)\,\mathsf du, $$ with $h(u)=\sqrt{\frac\pi 2} e^{-\frac12\sin^2 (\pi u) }$. It follows then that $I=\mathbb E[h(U)]$ with $U\sim\mathcal U(0,1)$. Let $U_i$ be i.i.d. $\mathcal U(0,1)$ random variables and set $X_i=h(U_i)$, then for each positive integer $n$ we have the point estimate $$\newcommand{\overbar}[1]{\mkern 1.75mu\overline{\mkern-1.75mu#1\mkern-1.75mu}\mkern 1.75mu} \widehat{I_n} =: \overbar X_n= \frac1n \sum_{i=1}^n X_i$$ and the approximate $1-\alpha$ confidence interval $$\overbar X_n\pm t_{n-1,\alpha/2}\frac{S_n}{\sqrt n}, $$ where $$S_n = \sqrt{\frac1{n-1}\sum_{i=1}^n \left(X_i-\overbar X_n\right)^2} $$ is the sample standard deviation.

Here is some $\texttt R$ code to estimate an integral using the Monte Carlo method:

# Define "h" function
hh <-function(u) {
  return(sqrt(0.5*pi) * exp(-0.5 * sin(pi*u)^2))
}

n <- 1000 # Number of trials
alpha <- 0.05 # Confidence level
U <- runif(n) # Generate U(0,1) variates
X <- hh(U) # Compute X_i's
Xbar <- mean(X) # Compute sample mean
Sn <- sqrt(1/(n-1) * sum((X-Xbar)^2)) # Compute sample stdev
CI <- (Xbar + (c(-1,1) * (qt(1-(0.5*alpha), n-1) * Sn/sqrt(n)))) # CI bounds

# Print results
cat(sprintf("Point estimate: %f\n", Xbar))
cat(sprintf("Confidence interval: (%f, %f)\n", CI[1], CI[2]))

For reference, the value of the integral (as computed by Mathematica) is $$e^{-\frac14}\sqrt{d\frac{\pi }{2}} I_0\left(\frac{1}{4}\right) \approx 0.991393, $$ where $I_\cdot(\cdot)$ denotes the modified Bessel function of the first kind, i.e. $$I_0\left(\frac14\right) = \frac1\pi\int_0^\pi e^{\frac14\cos\theta}\,\mathsf d\theta. $$

Related Question