> For the complete documentation index, see [llms.txt](https://zavalab.gitbook.io/statschemebook/llms.txt). Markdown versions of documentation pages are available by appending `.md` to page URLs; this page is available as [Markdown](https://zavalab.gitbook.io/statschemebook/chapter-4-estimation-for-random-variables/06-monte-carlo-approximations.md).
# Monte Carlo Approximations
* *Monte Carlo* (MC) is a computational technique that uses *random samples* to estimate properties of an RV and of its derived quantities (e.g., statistics).
* Consider a random sample $$x\_\omega \in \mathcal{S}$$ of the multivariate $$X$$. MC uses the samples to compute a sample (empirical) approximation of an *expectation*:
$$
\hat{\mathbb{E}}^S\_{\phi(X)}=\frac{1}{S}\sum\_{\omega \in \mathcal{S}}\phi(x\_\omega)\approx \mathbb{E}\[\phi(X)]
$$
where $$\phi(X)$$ is vector-valued function of $$X$$.
* We recall that the true (theoretical) expectation for a continuous RV is given by:
$$
\mathbb{E}\[\phi(X)]=\int\_{x\in \mathcal{D}\_X}\phi(x)f\_X(x)dx,
$$
where $$f\_X(x)$$ is the joint pdf of $$X$$.
* Computing expectation thus involves computing a multi-dimensional integral (which usually cannot be carried out analytically).
* MC provides a numerical technique to approximate such an integral.
* Many quantities of interest can be expressed as expectation operations and thus can be computed via MC.
* Some examples include:
* Expectation of $$X$$: $$\displaystyle\hat{\mathbb{E}}*X^S:=\frac{1}{S}\sum*{\omega \in \mathcal{S}}x\_\omega\approx \mathbb{E}\[X]$$
* Variance of $$X$$: $$\displaystyle\hat{\mathbb{V}}*X^S:=\frac{1}{S}\sum*{\omega \in \mathcal{S}}(x\_\omega-\hat{\mathbb{E}}\_X^S)^2\approx \mathbb{V}\[X]$$
* CDF of $$X$$: $$\displaystyle\hat{F}*X^S(x):=\frac{1}{S}\sum*{\omega \in \mathcal{S}}\mathbf{1}\[x\_\omega \leq x]\approx F\_X(x)$$
* Expectation of function: $$\displaystyle\hat{\mathbb{E}}*\varphi^S:=\frac{1}{S}\sum*{\omega \in \mathcal{S}}\varphi(x\_\omega,u)\approx \mathbb{E}\[{\varphi(X,u)}]$$.
* Probability of event: $$\hat{\mathbb{P}}(X\in\mathcal{A}):=\frac{1}{S}\sum\_{\omega \in \mathcal{S}}\mathbf{1}\[x\_\omega \in \mathcal{A}] \approx \mathbb{P}(X\in\mathcal{A})$$.
* Note dependence of empirical approxs on number of samples $$S$$.
* When empirical approx use data, these are called data-driven approxs.
* Empirical approximations are also often called sample average approxs.
* Important questions related to MC approximations are:
* Do the approximations become exact as we accumulate data?
* How accurate are these approximations when we have a limited amount of data?
* How much data do we need to obtain a desired accuracy?
* How do approximation errors behave as we change our sample data?
* We will now explore tools to answer these questions; the use of these tools is enabled by the fact that samples are collected at random.
* As part of our discussion we will also see that RVs exhibit “universal" behavior as we accumulate data (i.e., as $$S\to \infty$$).
## Law of Large Numbers
* Let's assess quality of empirical approximations as $$S\to \infty$$.
* Consider an i.i.d. random sequence $$X\_1,X\_2,...,X\_S$$ for univariate RV $$X$$. Since the samples are i.i.d., they have same underlying pdf with $$\mathbb{E}\[X]$$.
* Since samples are random, they are unbiased $$\mathbb{E}\[X\_1]=\mathbb{E}\[X\_2]=\cdots=\mathbb{E}\[X\_S]=\mathbb{E}\[X]$$.
* Consider now the MC approximation of $$\mathbb{E}\[X]$$:
$$
\hat{\mathbb{E}}*X^S=\frac{1}{S}\sum*{\omega \in \mathcal{S}}X\_\omega
$$
Because the samples $$X\_\omega,,\omega \in \mathcal{S}$$ are random, approx $$\hat{\mathbb{E}}\_X^S$$ is also random (i.e., different samples give different approxs).
* There is thus inherent variability in approx $$\hat{\mathbb{E}}\_X^S$$ (variability is a measure of approximation accuracy).
* Law of large numbers (LLN) indicates that:
$$
\lim\_{S\to \infty}\hat{\mathbb{E}}\_X^S=\mathbb{E}\[X]
$$
* LLN guarantees *stable long-run* behavior of sample averages.
## Example: Monte Carlo Approximations
{% hint style="danger" icon="flask" %}
* We generate random samples $$X\_\omega,; \omega \in \mathcal{S}$$ from $$X\sim \textrm{Weibull}(1,2)$$.
* Compute MC approxs for $$\mathbb{E}\[X]$$, $$\mathbb{E}\[\log(X)]$$, $$\mathbb{V}\[\exp(X)+X^2]$$, and $$\mathbb{P}(X\geq 1)$$ for different $$S$$.
{% endhint %}
#### Python implementation
{% hint style="success" %}
Notebook [`ch4_montecarlo_approx.ipynb`](https://colab.research.google.com/github/vzavala/StatsBook/blob/master/Python/ch4_montecarlo_approx.ipynb) — open it in Colab and run the cells alongside this section.
MATLAB live script: [`ch4_montecarlo_approx.mlx`](https://github.com/vzavala/StatsBook/blob/master/Matlab-Live/ch4_montecarlo_approx.mlx)
{% endhint %}
Convergence of Monte Carlo approximations (Victor Z, UW-Madison, 2024). Samples come from a Weibull RV with scale $$\beta=2$$ and shape $$\xi=1$$ (as in the slide figure; the slide text writes "Weibull(2,1)" in Matlab's (scale, shape) order).
| Slides | Python | Meaning |
| ------------------------------------------------------------------------ | -------------------- | ----------------------------------------------------------- |
| $$S$$ | `S_values` | sample sizes |
| $$\hat{\mathbb{E}}^S\_{\phi(X)}=\frac{1}{S}\sum\_\omega\phi(x\_\omega)$$ | `E_hat`, `E_log_hat` | MC estimates of $$\mathbb{E}\[X]$$, $$\mathbb{E}\[\log X]$$ |
| $$\hat{\mathbb{V}}^S$$ | `V_hat` | MC estimate of $$\mathbb{V}\[e^{-X}+X^2]$$ |
| $$\hat{\mathbb{P}}(X\in\mathcal{A})$$ | `P_hat` | MC estimate of $$\mathbb{P}(X\leq 1)$$ |
We use $$e^{-X}$$ instead of $$e^{X}$$: with this RV, $$\mathbb{V}\[e^{X}]$$ is infinite.
**Note:** This notebook reproduces Example 4.6 (Section 4.3.2) of *Statistics for Chemical Engineers*: Monte Carlo approximations of several quantities as the number of samples S grows.
```python
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import weibull_min
xi, beta = 1, 2
X = weibull_min(xi, scale=beta)
np.random.seed(1)
S_values = np.arange(1, 5001, 10)
E_hat, E_log_hat, V_hat, P_hat = [], [], [], []
for S in S_values:
x_omega = X.rvs(size=S)
E_hat.append(np.mean(x_omega))
E_log_hat.append(np.mean(np.log(x_omega)))
V_hat.append(np.var(np.exp(-x_omega) + x_omega**2))
P_hat.append(np.mean(x_omega <= 1))
```
```python
exact = {'E': X.mean(), 'P': X.cdf(1)}
fig, axs = plt.subplots(2, 2, figsize=(10, 8))
panels = [(E_hat, r'$\hat{\mathbb{E}}[X]$', exact['E']), (E_log_hat, r'$\hat{\mathbb{E}}[\log X]$', None),
(V_hat, r'$\hat{\mathbb{V}}[e^{-X}+X^2]$', None), (P_hat, r'$\hat{\mathbb{P}}(X\leq 1)$', exact['P'])]
for ax, (vals, label, ref) in zip(axs.ravel(), panels):
ax.plot(S_values, vals, 'o', markerfacecolor='w', markeredgecolor='black', markersize=4)
if ref is not None:
ax.axhline(ref, color='k')
ax.set_xlabel('$S$')
ax.set_ylabel(label)
ax.grid(True)
plt.tight_layout()
plt.show()
print(f"exact E[X] = {exact['E']:.3f}, exact P(X <= 1) = {exact['P']:.3f}")
```
```
exact E[X] = 2.000, exact P(X <= 1) = 0.393
```
**Note:** As in the book, the approximations are unstable for small samples (S below about 100) and become increasingly accurate as S grows (law of large numbers).
