Find an estimate for $\int_{\pi/2}^\pi \sin(x) dx$ using the Monte Carlo Simulation

definite integralsintegrationmonte carlosimulationtransformation

I want to find an estimate for

$$\int_{\pi/2}^\pi \sin(x) dx$$

I want to use the monte carlo simulation method. I've plotted the graph of $\sin(x)$ in the given interval. The total area

$$\begin{align*}
A &= 1\cdot \frac{\pi}{2}\\[5pt]
&=\frac{\pi}{2}\\[5pt]
\end{align*}$$

The target area $a = \big(\frac{n}{N}\big) \cdot A$ where $n$ is the number of hits and $N$ is the total number of 'throws'. I've been given the following random numbers:

$$\begin {align}
8442 \quad 5653\\
2887 \quad 7703\\
6412 \quad 3965\\
4941 \quad 7376\\
0646 \quad 6931\\
7556 \quad 9697\\
\end {align}$$

What I need to do now is transform my $x$ values so they are in the range $\big(\frac{\pi}{2}, \pi\big)$, the graph only goes up to $y = 1$ so they $y$ values don't need to be transformed. Where I'm struggling is what transformation I should use on the random $x$-values, I tried $x_{transformed} = x_{random} +\big(\frac{\pi}{2}\big)$ but I didn't get any hits, my $\sin(x_{transformed})$ values were very low.

If anyone could help out and show me what transformation I should use on the $x$ values to get them in the range I require it would be much appreciated.

Thanks in advance

(I'll add the graph of $\sin(x)$ in the interval I require so you can see the rectangle and the target area:)

Best Answer

It looks like your $x_{random}$ values are from $0$ to $10000$. So dividing by $10000$, $x_{random}/10000$ is in the $[0,1]$ interval. Now multiply by the length of the desired interval, and shift it to the beginning of that interval: $$\frac{x_{random}}{10000}\frac{\pi}2$$ is in the $[0,\pi/2]$ interval $$\frac{x_{random}}{10000}\frac{\pi}2+\frac\pi2$$ is in the $[\pi/2,\pi]$ interval

Related Solutions

[Math] Integration using Monte Carlo Method

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

The Monte Carlo Method For Simulation Example Problem

In your stated problem, "counting" is counting weeks where any tie was chosen twice.

After you do your trials, rolling the dice for X number of week sequences, you will end up with a number of weeks that a tie was chosen twice / number of week sequences you ran for your Monte Carlo series, and this value will be your estimated probability that in a future week a tie will be chosen twice.

The idea of Monte Carlo simulations is that when you have some sort of random set of events that are too large (or perhaps infinite) to exhaustively calculate all finite probabilities, so by running a number of random tests you can start to approximate the probability. The larger your number of tests, the better the approximation is.

Your particular assigned problem is pretty straightforward to calculate exhaustively, though to list out the 3125 different permutations would be a bit tedious. What the assignment is trying to show you is how many "weeks" you have to simulate to start to approach the final probability.

Related Question