Estimate the value of the following integral using Monte Carlo integration:
\[ I = \int_0^1 e^{-x^2}\,dx. \]
This integral does not have a closed-form solution in terms of elementary functions, making it a good candidate for Monte Carlo approximation. We can rewrite the integral as an expectation: \[ I = \int_0^1 e^{-x^2} \cdot 1\,dx = \mathbb{E}[e^{-X^2}], \] where \(X \sim \text{Uniform}(0,1)\) and \(f(x) = 1\) for \(x \in [0,1]\). Therefore, the Monte Carlo estimator for \(I\) is \[ \hat I_n = \frac{1}{n}\sum_{i=1}^n e^{-X_i^2}. \]
We can implement this estimator as follows:
import numpy as np
from scipy.integrate import quad
np.random.seed(1234)
n = 10000
x = np.random.uniform(0, 1, n)
I_hat = np.mean(np.exp(-x**2))
I, _ = quad(lambda x: np.exp(-x**2), 0, 1)
print("Monte Carlo estimate: ", round(I_hat, 5),
"\nAnalytical value: ", round(I, 5), sep="")Monte Carlo estimate: 0.74638
Analytical value: 0.74682Estimate the value of the following probability using Monte Carlo integration:
\[ P(X > 2.5), \quad X \sim N(0,1). \]
This probability can be expressed as an integral of the normal density function: \[ P(X > 2.5) = \int_{2.5}^\infty \frac{1}{\sqrt{2\pi}} e^{-x^2/2}\,dx. \] This integral does not have a closed-form solution, but we can estimate it using Monte Carlo integration. To do this, we can generate samples from the standard normal distribution and compute the proportion of samples that exceed 2.5. \[ P(X > 2.5) = \int_{2.5}^\infty f(x)\,dx = \mathbb{E}[\mathbf{1}_{\{X > 2.5\}}], \] where \(\mathbf{1}_{\{X > 2.5\}}\) is an indicator function that equals 1 if \(X > 2.5\) and 0 otherwise. The Monte Carlo estimator for \(P(X > 2.5)\) is therefore \[ \hat P_n = \frac{1}{n} \sum_{i=1}^n \mathbf{1}_{\{X_i > 2.5\}}, \] where \(X_i\) are i.i.d. samples drawn from the standard normal distribution.
We can implement this estimator as follows:
import numpy as np
from scipy.stats import norm
np.random.seed(1234)
n = 5000
x = np.random.normal(0, 1, n)
p_hat = np.mean(x > 2.5)
p_true = 1 - norm.cdf(2.5)
print("Monte Carlo estimate: ", round(p_hat, 5),
"\nAnalytical value: ", round(p_true, 5), sep="")Monte Carlo estimate: 0.0052
Analytical value: 0.00621Estimate the area under the curve:
\[ y = \sin(x), \quad 0 \le x \le \pi. \]
The area under the curve can be computed analytically using integration: \[ A = \int_0^\pi \sin(x)\,dx = \left[-\cos(x)\right]_0^\pi = 2. \]
Since \(X\sim\text{Uniform}(0, \pi)\), we can rewrite the integral as an expectation: \[ A = \int_0^\pi \pi\sin(x) \cdot \frac{1}{\pi}\,dx = \mathbb{E}[\pi\sin(X)], \] where \(X \sim \text{Uniform}(0, \pi)\). To estimate the area under the curve using Monte Carlo integration, we can use the following estimator: \[ \hat A_n = \frac{\pi}{n} \sum_{i=1}^n \sin(X_i), \] where \(X_i\) are i.i.d. samples drawn uniformly from the interval \([0, \pi]\).
We can implement this estimator as follows:
import numpy as np
from scipy.integrate import quad
np.random.seed(1234)
n = 5000
x = np.random.uniform(0, np.pi, n)
A_hat = (np.pi/n) * np.sum(np.sin(x))
A_true, _ = quad(np.sin, 0, np.pi)
print("Monte Carlo estimate: ", round(A_hat, 5),
"\nAnalytical value: ", round(A_true, 5), sep="")Monte Carlo estimate: 1.98564
Analytical value: 2.0Estimate \(\mathbb{E}[X^2]\) where \(X \sim \text{Exp}(\lambda = 1)\).
The expectation \(\mathbb{E}[X^2]\) for an exponential distribution with rate parameter \(\lambda = 1\) can be computed using the formula for the \(n\)-th moment of an exponential distribution: \[ \mathbb{E}[X^n] = \frac{n!}{\lambda^n}. \] For \(n=2\), we have \[ \mathbb{E}[X^2] = \frac{2!}{1^2} = 2. \]
To estimate \(\mathbb{E}[X^2]\) using Monte Carlo integration, we can generate samples from the exponential distribution and compute the average of their squares. The Monte Carlo estimator is given by: \[ \hat E_n = \frac{1}{n} \sum_{i=1}^n X_i^2, \] where \(X_i\) are i.i.d. samples drawn from \(\text{Exp}(1)\).
We can implement this estimator as follows:
import numpy as np
np.random.seed(1234)
n = 5000
x = np.random.exponential(scale=1, size=n)
E_hat = np.mean(x**2)
E_true = 2
print("Monte Carlo estimate: ", round(E_hat, 5),
"\nAnalytical value: ", round(E_true, 5), sep="")Monte Carlo estimate: 1.96882
Analytical value: 2Estimate the value of the following integral using Monte Carlo integration:
\[ I = \int_{[0,1]^5} \exp\left(-\sum_{i=1}^5 x_i^2\right) dx_1 dx_2 dx_3 dx_4 dx_5. \]
This is a five-dimensional integral, which can be challenging to compute using traditional numerical integration methods. Use Monte Carlo integration to estimate the value of \(I\).
We can rewrite the integral as an expectation: \[ I = \int_{[0,1]^5} \exp\left(-\sum_{i=1}^5 x_i^2\right) dx_1 dx_2 dx_3 dx_4 dx_5 = \mathbb{E}\left[\exp\left(-\sum_{i=1}^5 X_i^2\right)\right], \] where \(X_i \sim \text{Uniform}(0,1)\) for \(i=1,\ldots,5\). The Monte Carlo estimator for \(I\) is therefore \[ \hat I_n = \frac{1}{n} \sum_{j=1}^n \exp\left(-\sum_{i=1}^5 X_{ij}^2\right), \] where \(X_{ij}\) are i.i.d. samples drawn uniformly from the interval \([0, 1]\) for each dimension.
We can implement this estimator as follows:
import numpy as np
np.random.seed(123)
n = 10000
X = np.random.uniform(0, 1, (n, 5))
g = np.exp(-np.sum(X**2, axis=1))
I_hat = np.mean(g)
print("Monte Carlo estimate: ", round(I_hat, 5))Monte Carlo estimate: 0.2307This illustrates an important idea:
Monte Carlo integration complexity depends mainly on the number of simulations, not strongly on the dimension of the integral.
This makes Monte Carlo methods especially valuable in:
Suppose a machine fails if the stress S exceeds the strength R.
Let
\[ S \sim N(50, 10^2), \quad R \sim N(70, 15^2), \]
independetly. Estimate the probability that the machine fails, i.e., \(P(S > R)\), using Monte Carlo simulation.
We can express the probability of failure as an expectation: \[ P(S > R) = \mathbb{E}[\mathbf{1}_{\{S > R\}}], \] where \(\mathbf{1}_{\{S > R\}}\) is an indicator function that equals 1 if \(S > R\) and 0 otherwise. To estimate this probability using Monte Carlo simulation, we can generate pairs of samples \((S_i, R_i)\) from their respective distributions and compute the proportion of pairs for which \(S_i > R_i\). The Monte Carlo estimator for \(P(S > R)\) is therefore \[ \hat P_n = \frac{1}{n} \sum_{i=1}^n \mathbf{1}_{\{S_i > R_i\}}, \] where \(S_i \sim N(50, 10^2)\) and \(R_i \sim N(70, 15^2)\) are i.i.d. samples.
We can implement this estimator as follows:
import numpy as np
from scipy.stats import norm
np.random.seed(1234)
n = 10000
S = np.random.normal(50, 10, n)
R = np.random.normal(70, 15, n)
p_hat = np.mean(S > R)
p_true = 1 - norm.cdf(0, loc=50-70, scale=np.sqrt(10**2 + 15**2)) # independent normals
print("Monte Carlo estimate: ", round(p_hat, 5),
"\nAnalytical value: ", round(p_true, 5), sep="")Monte Carlo estimate: 0.1349
Analytical value: 0.13363