2.12. Exemple numérique: Importance Sampling

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[)\).

Le changement de loi proposé est celui de passer de la loi uniforme sur \(]0,1[\) à une loi de densité \(g(u) = 2 (1-u) \mathbf{1}_{]0,1[}(u)\). En effet on sait qu’il est optimal de changer de densité avec une densité proche de \(\frac{1}{I} \cos(\frac{\pi}{2})(u) \mathbf{1}_{]0,1[}(u)\).

Si on considère \(V \sim g\) alors

\[\begin{equation*} I = \mathbf{E} \Big[ \cos\big(\frac{\pi}{2} V \big) \frac{1}{g(V)} \Big] \end{equation*}\]

et on considère l’estimateur de Monte Carlo

\[\begin{equation*} \tilde I_n = \frac{1}{n} \sum_{k=1}^n \cos\big(\frac{\pi}{2} V_k \big) \frac{1}{2(1-V_k)}, \end{equation*}\]

où \((V_k)_{k \ge 1}\) est une suite i.i.d. de copies de \(V \sim g\).

Show code cell source 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)

def g(u): 
    return 2*(1-u)

result = pd.DataFrame()

n = int(1e5)
sample_U = rng.random(size=n)
sample_MC = f(sample_U)

sample_V = 1-np.sqrt(1-sample_U)
sample_IS = f(sample_V) / g(sample_V)

result = pd.DataFrame([
    monte_carlo(sample_MC), 
    monte_carlo(sample_IS)
], index = ["MC", "MC_IS"])
result

	mean	var	lower	upper
MC	0.637571	0.094342	0.635667	0.639475
MC_IS	0.636422	0.006759	0.635913	0.636932

ratio_var = result["var"]["MC"] / result["var"]["MC_IS"]
print("Le ratio de variance est:", ratio_var)

Le ratio de variance est: 13.958903901596605

Un deuxième choix est de considérer la densité auxiliaire

\[\begin{equation*} \tilde g(u) = \frac{3}{2} (1 - u^2) \mathbf{1}_{]0,1[}(u), \end{equation*}\]

encore plus proche de \(f(u) = \cos(\frac{\pi}{2} u)\) sur \(]0,1[\). L’estimateur MC avec cette densité auxiliaire s’écrit alors

\[\begin{equation*} \tilde I_n = \frac{1}{n} \sum_{k=1}^n \cos\big(\frac{\pi}{2} \tilde V_k \big) \frac{1}{\frac{3}{2}(1- \tilde V_k^2)}, \end{equation*}\]

où \((\tilde V_k)_{k \ge 1}\) est une suite i.i.d. de copies de \(V \sim \tilde g\).

Pour simuler une variable aléatoire de densité \(\tilde g\) on utilise la méthode du rejet par rapport à la loi uniforme. En effet, sur \([0,1]\) la fonction \(\tilde g\) est majorée par \(C = \frac{3}{2}\). On considère alors une suite de v.a. auxiliaire \((Z_k)_{k \ge 1}\) uniforme sur \(]0,1[\) pour faire le test de rejet. On accepte les réalisations \(U_k\) qui vérifient \(Z_k \le 1 - U_k^2\), \(k \ge 1\).

Comme \(C = \frac{3}{2}\) on augmente le nombre \(n\) de réalisations d’un facteur \(\frac{3}{2}\) pour garder la même taille d’échantillon et la même complexité (nombre d’évaluations de la fonction \(f\)).

def g2(u): 
    return (3/2)*(1-u**2)

C = 3/2
sample_U = rng.random(size=int(C*n))
sample_Zaux = rng.random(size=int(C*n))
accepted = sample_Zaux <= 1 - sample_U**2
sample_V = sample_U[accepted]
len(sample_V)

99862

result = pd.concat([result, 
                    pd.DataFrame([monte_carlo(f(sample_V) / g2(sample_V))], 
                                 index=["MC_IS2"])])
result

	mean	var	lower	upper
MC	0.637571	0.094342	0.635667	0.639475
MC_IS	0.636422	0.006759	0.635913	0.636932
MC_IS2	0.636592	0.000994	0.636397	0.636788

ratio_var = result["var"]["MC"] / result["var"]["MC_IS2"]
print("Le ratio de variance est:", ratio_var)

Le ratio de variance est: 94.87659887164943