1
$\begingroup$

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

$\endgroup$
  • $\begingroup$ Ok. How did you find the confidence intervals? $\endgroup$ – Fawcett512 Aug 31 '16 at 4:04
  • $\begingroup$ 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$ $\endgroup$ – Fawcett512 Aug 31 '16 at 4:06
  • $\begingroup$ 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. $\endgroup$ – Fawcett512 Aug 31 '16 at 4:40
  • $\begingroup$ 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. $\endgroup$ – Fawcett512 Aug 31 '16 at 4:47
1
$\begingroup$

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

$\endgroup$
  • 1
    $\begingroup$ Thank you. This was very clear. $\endgroup$ – Fawcett512 Aug 31 '16 at 4:54
  • $\begingroup$ You're welcome! The change of variables is not strictly necessary, but it simplifies the computations considerably. $\endgroup$ – Math1000 Aug 31 '16 at 4:58
  • $\begingroup$ 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? $\endgroup$ – Fawcett512 Aug 31 '16 at 6:54
  • $\begingroup$ 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. $\endgroup$ – Math1000 Aug 31 '16 at 7:49
  • 1
    $\begingroup$ 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). $\endgroup$ – Fawcett512 Sep 1 '16 at 5:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.