# Integration using Monte Carlo Method

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$?.

• Ok. How did you find the confidence intervals? – Fawcett512 Aug 31 '16 at 4:04
• I know how to get the confidence intervals for a sample of uniform random variables, but I'm not sure how to find them in this case given that I don't have a parameter but rather the approximate value of $I$ – Fawcett512 Aug 31 '16 at 4:06
• I think I'm making the problem more complicate that it should. So I could generate $N$ uniform random variables $x_i$ and then compute $\mathbb{E}=\frac{\pi}{N}\sum f(x_i)$. As $N$ tends to infinity, we should get closer to the integral. Yet I don't see how to get the confidence intervals. – Fawcett512 Aug 31 '16 at 4:40
• I think that I was trying to do some importance sampling (hence the question if we could use a normal probability density function), but seems that in this case that is not possible, Furthermore, if I just compute the above formula for $\mathbb{E}$,I should get a good estimate for $N$ provided that $N$ is big enough. The only thing that's not very clear is how to compute a 95% confidence interval. – Fawcett512 Aug 31 '16 at 4:47

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.$$

• Thank you. This was very clear. – Fawcett512 Aug 31 '16 at 4:54
• You're welcome! The change of variables is not strictly necessary, but it simplifies the computations considerably. – Math1000 Aug 31 '16 at 4:58
• Just one more question because I have trouble understanding the concept of confidence intervals. How do I know the size of the sample if now I want to get a 98% confidence interval? – Fawcett512 Aug 31 '16 at 6:54
• I don't understand your question - the sample size ($n$) and confidence level ($\alpha$) are parameters which influence the width of the confidence interval, but are not directly related to each other. – Math1000 Aug 31 '16 at 7:49
• Never mind. That was a bad question. I was wondering how big the sample must be (because I was making a simulation) to get a 98% confidence interval, but as you said is a parameter which influences the width of the confidence interval. All is clear know (or at least I think so). – Fawcett512 Sep 1 '16 at 5:59