2.6. Exemple numérique: MC avec variables antithétiques#
On considère le calcul de \(I = \mathbf{E} \Big[\cos\big(\frac{\pi}{2} U\big) \Big]\) où \(U \sim \mathcal{U}(]0,1[)\) qu’on réécrit comme
\[\begin{equation*}
I = \mathbf{E} \bigg[\frac{1}{2} \Big( \cos\big(\frac{\pi}{2} U\big) + \cos\big(\frac{\pi}{2} (1-U)\big) \Big) \bigg]
\end{equation*}\]
L’estimateur de Monte Carlo avec variables antithétiques s’écrit alors
\[\begin{equation*}
I_n = \frac{1}{n} \sum_{k=1}^n \frac{1}{2} \Big(\cos\big(\frac{\pi}{2} U_k \big)
+ \cos\big(\frac{\pi}{2} (1-U_k) \big)\Big),
\end{equation*}\]
où \((U_k)_{k \ge 1}\) est une suite i.i.d. uniforme sur \(]0,1[\). Voici le code de cet estimateur.
Show code cell source
from scipy import stats
def monte_carlo(sample, proba = 0.95):
"""
Computes the mean, variance, and a confidence interval (CI) of a
given sample data set using the Monte Carlo method.
Parameters:
-----------
sample : array-like
The data set to be analyzed
proba : float, optional
The probability that the true mean of the population is
within the calculated interval. Default is 0.95
Returns:
--------
dict of { 'mean': float, 'var': float, 'lower': float, 'upper': float }
The mean, variance, lower bound of the CI and upper bound of the CI
"""
mean = np.mean(sample)
var = np.var(sample, ddof=1)
alpha = 1 - proba
quantile = stats.norm.ppf(1 - alpha/2) # fonction quantile
ci_size = quantile * np.sqrt(var / sample.size)
result = { 'mean': mean, 'var': var,
'lower': mean - ci_size,
'upper': mean + ci_size }
return result
Show code cell source
from math import pi
def f(u):
return np.cos(0.5*pi*u)
# on crée un DataFrame vide pour stocker nos résultats
result = pd.DataFrame()
n = int(1e5)
sample_U = rng.random(size=n)
result = pd.DataFrame([
monte_carlo(f(sample_U)),
monte_carlo(0.5 * (f(sample_U) + f(1-sample_U)))
], index= ["MC", "MC_VA"])
result
| mean | var | lower | upper | |
|---|---|---|---|---|
| MC | 0.635883 | 0.094952 | 0.633973 | 0.637793 |
| MC_VA | 0.636591 | 0.003863 | 0.636206 | 0.636976 |
ratio_var = result["var"]["MC"] / result["var"]["MC_VA"]
print("Le ratio de variance est de:" , ratio_var)
Le ratio de variance est de: 24.57896856396313