Ctrl+K
LPSM

Vincent Lemaire

Probabilités numériques
vers les tps

Exemple d’estimateur de Monte Carlo

2.3. Exemple d’estimateur de Monte Carlo#

On considère le calcul de \(I = \int_0^1 \cos\big(\frac{\pi}{2} u\big) \operatorname{d}\!u\). On suppose qu’on ne connait pas \(I\), et on veut approcher cette quantité par une méthode de Monte Carlo. Pour cela, on utilise la représentation probabiliste \(I = \mathbf{E} \Big[\cos\big(\frac{\pi}{2} U\big) \Big]\) où \(U \sim \mathcal{U}(]0,1[)\).

L’estimateur de Monte Carlo s’écrit alors

\[\begin{equation*} I_n = \frac{1}{n} \sum_{k=1}^n \cos\big(\frac{\pi}{2} U_k \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.

Hide 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

from math import pi
def f(u):
    return np.cos(0.5*pi*u)

n = int(1e3)
sample_U = rng.random(size=n)
sample = f(sample_U)
result = pd.DataFrame([monte_carlo(sample)], index=[n])
result
mean var lower upper
1000 0.62861 0.092616 0.609748 0.647472 
result["upper"][n] - result["lower"][n]

0.037724304864996894

size_IC = result["upper"][n] - result["lower"][n]
print("La taille de l'IC est:", size_IC)

La taille de l'IC est: 0.037724304864996894

2.3.1. On veut réduire l’erreur#

D’après le TCL la vitesse est en \(\sqrt{n}\) ce qui signifie que pour diviser l’erreur par \(p=10\) il faut multiplier \(n\) par \(p^2=100\). Vérifions empiriquement ce résultat.

sample_U = rng.random(size=100*n)
result = pd.concat([result, 
                    pd.DataFrame([monte_carlo(f(sample_U))], index=[100*n])])
result
mean var lower upper
1000 0.628610 0.092616 0.609748 0.647472
100000 0.637107 0.095160 0.635195 0.639019 
size_IC = result["upper"][100*n] - result["lower"][100*n]
print("La taille de l'IC est:", size_IC)

La taille de l'IC est: 0.0038238882417132025

2.3.2. On veut atteindre une erreur#

Si on veut que la taille de l’intervalle de confiance soit d’un ordre \(\varepsilon\) fixé, il faut choisir le bon nombre de réalisations \(n\). Si on a une bonne approximation de la variance \(\sigma^2\) il suffit de choisir

\[\begin{equation*} n \simeq (2 \times 1.96)^2 \sigma^2 \varepsilon^{-2}. \end{equation*}\]

Le terme 1.96 vient du quantile de la loi normale pour construire l’IC à 95% et le facteur 2 car on veut la taille totale de l’IC.

Voici un exemple pour atteindre une erreur de \(10^{-4}\).

epsilon = 1e-4
n_eps = int((2*1.96/epsilon)**2 * result["var"][100*n])
print(n_eps)
sample_U = rng.random(size=n_eps)
result = pd.concat([result, 
                    pd.DataFrame([monte_carlo(f(sample_U))], index=[n_eps])])
result

146226586
mean var lower upper
1000 0.628610 0.092616 0.609748 0.647472
100000 0.637107 0.095160 0.635195 0.639019
146226586 0.636618 0.094724 0.636568 0.636668 
size_IC = result["upper"][n_eps] - result["lower"][n_eps]
print("La taille de l'IC est:", size_IC)

La taille de l'IC est: 9.97687322641383e-05
Contents