}}
Implémenter la méthode de Monte-Carlo pour le calcul approché de \pi via les deux intégrales suivantes :
import numpy as np
from scipy import stats
n = 10000
x = np.random.uniform(0, 1, n)
y = 4 * np.sqrt(1 - x**2)
in_ = np.mean(y)
q = stats.norm.ppf(0.975)
delta = q * np.sqrt(np.var(y) / n)
ic = (in_ - delta, in_ + delta)
print("IN:", in_)
print("Confidence Interval (IC):", ic)IN: 3.140139698393205
Confidence Interval (IC): (3.12252451384522, 3.15775488294119)import numpy as np
# Define the number of points
n = 10000
# Generate random numbers for X and Y in the range [-1, 1]
x = np.random.uniform(-1, 1, n)
y = np.random.uniform(-1, 1, n)
# Calculate Z based on the condition (X^2 + Y^2 <= 1)
z = 4 * (x**2 + y**2 <= 1)
# Calculate the mean of Z
in_ = np.mean(z)
# Calculate the confidence interval (IC)
q = stats.norm.ppf(0.975)
delta = q * np.sqrt(np.var(z) / n)
ic = (in_ - delta, in_ + delta)
print("IN:", in_)
print("Confidence Interval (IC):", ic)IN: 3.1552
Confidence Interval (IC): (3.123200846461404, 3.1871991535385957)La première méthode donne une meilleure précision. Mais est-ce si concluant. Il faut regarder plus en détail les variance pour conclure.
On souhaite estimer la probabilité qu’une variable aléatoire X \sim \text{Cauchy}(0, 1) soit plus grande que 2, i.e., I = \mathbb P(X\geq 2) = \int_2^\infty \left\{ \pi (1 + x^2) \right\}^{-1} \text{d$x$} = - \frac{\arctan 2}{\pi} + \frac{1}{2}
Python la valeur exacte de l’intégrale.print(-np.arctan(2) / np.pi + 1/2)0.14758361765043326import numpy as np
# Define the number of samples
N = 10000
# Generate samples from the Cauchy distribution
X = np.random.standard_cauchy(N)
# Calculate Y based on the condition (X >= 2)
Y = (X >= 2)
# Calculate the mean of Y
in_ = np.mean(Y)
# Calculate the confidence interval (IC)
q = stats.norm.ppf(0.975)
delta = q * np.sqrt(np.var(Y) / n)
ic = (in_ - delta, in_ + delta)
print("IN:", in_)
print("Confidence Interval (IC):", ic)IN: 0.1502
Confidence Interval (IC): (0.1431976885943637, 0.1572023114056363)import numpy as np
# Define the number of samples
N = 10000
# Generate samples from the Cauchy distribution
X = np.random.standard_cauchy(N)
# Calculate Z based on the condition ((X >= 2) + (-X >= 2)) / 2
Z = ((X >= 2) + (-X >= 2)) / 2
# Calculate the mean of Z
in_ = np.mean(Z)
# Calculate the confidence interval (IC)
delta = q * np.sqrt(np.var(Z) / n)
ic = (in_ - delta, in_ + delta)
print("IN:", in_)
print("Confidence Interval (IC):", ic)IN: 0.1496
Confidence Interval (IC): (0.14511258808724478, 0.15408741191275524)L’intervalle de confiance est un peu plus resserré. La variance a été divisée par
np.var(Y) / np.var(Z)2.4349551620149925import numpy as np
# Define the number of samples
N = 10000
# Generate samples from the uniform distribution in the range [0, 2]
U = np.random.uniform(0, 2, N)
# Calculate A based on the given formula
A = 1/2 - 2 / (np.pi * (1 + U**2))
# Calculate the mean of A
in_ = np.mean(A)
# Calculate the confidence interval (IC)
q = stats.norm.ppf(0.975)
delta = q * np.sqrt(np.var(A) / n)
ic = (in_ - delta, in_ + delta)
print("IN:", in_)
print("Confidence Interval (IC):", ic)IN: 0.14915260765508548
Confidence Interval (IC): (0.14584427089376595, 0.152460944416405)L’intervalle de confiance est encore plus resserré. Par rapport à la méthode initiale, la variance a été divisée par
np.var(Y) / np.var(A)4.479849457246015import numpy as np
# Define the number of samples
N = 10000
# Generate samples from the uniform distribution in the range [0, 1/2]
U = np.random.uniform(0, 0.5, N)
# Calculate B based on the given formula
B = 1 / (2 * np.pi * (1 + U**2))
# Calculate the mean of B
in_ = np.mean(B)
q = stats.norm.ppf(0.975)
delta = q * np.sqrt(np.var(B) / n)
ic = (in_ - delta, in_ + delta)
# Calculate the confidence interval (IC)
print("IN:", in_)
print("Confidence Interval (IC):", ic)IN: 0.14756022988503692
Confidence Interval (IC): (0.14736811492855006, 0.14775234484152378)L’intervalle de confiance est encore plus resserré, on arrive presque à une précision de 10^{-3}. Par rapport à la méthode initiale, la variance a été divisée par
np.var(Y) / np.var(B)1328.496724392179Celle qui a la plus petite variance, donc la dernière.
import numpy as np
def calculate_statistics(n):
# Generate samples for x and th
x = np.random.uniform(0, 0.5, n)
th = np.random.uniform(-np.pi/2, np.pi/2, n)
# Calculate y based on the given condition
y = (np.cos(th) >= 2 * x) / 2
# Calculate the mean of y
in_ = np.mean(y)
q = stats.norm.ppf(0.975)
# Calculate the confidence interval (ic) for in_
delta = q * np.sqrt(np.var(y) / n)
ic = (in_ - delta, in_ + delta)
# Calculate pichap
pichap = 1 / in_
# Calculate icchap
icchap = (1 / ic[1], 1 / ic[0])
return in_, ic, pichap, icchap
N = 10000
IN, IC, pichap, ICchap = calculate_statistics(N)
print("IN:", IN)
print("Confidence Interval (IC) for IN:", IC)
print("pichap:", pichap)
print("Confidence Interval (IC) for pichap:", ICchap)IN: 0.32175
Confidence Interval (IC) for IN: (0.31705622639879943, 0.3264437736012005)
pichap: 3.1080031080031083
Confidence Interval (IC) for pichap: (3.0633146681536902, 3.1540147038215887)# Generate n values logarithmically spaced from 1e5 to 1e8
n_values = np.logspace(5, 8, num=10, dtype=int)
# Evaluate the function for each value of n
results = [calculate_statistics(n) for n in n_values]
# Print the results
for n, (IN, IC, pichap, ICchap) in zip(n_values, results):
print(f"n: {n:09d}, IN: {IN:.3f}, IC: ({IC[0]:.3f}, {IC[1]:.3f}), pichap: {pichap:.3f}, ICchap: ({ICchap[0]:.3f}, {ICchap[1]:.3f})")n: 000100000, IN: 0.319, IC: (0.318, 0.321), pichap: 3.133, ICchap: (3.119, 3.148)
n: 000215443, IN: 0.318, IC: (0.317, 0.319), pichap: 3.148, ICchap: (3.138, 3.158)
n: 000464158, IN: 0.317, IC: (0.317, 0.318), pichap: 3.150, ICchap: (3.144, 3.157)
n: 001000000, IN: 0.318, IC: (0.318, 0.319), pichap: 3.140, ICchap: (3.135, 3.145)
n: 002154434, IN: 0.318, IC: (0.318, 0.319), pichap: 3.141, ICchap: (3.138, 3.144)
n: 004641588, IN: 0.318, IC: (0.318, 0.319), pichap: 3.141, ICchap: (3.139, 3.143)
n: 010000000, IN: 0.318, IC: (0.318, 0.319), pichap: 3.141, ICchap: (3.139, 3.142)
n: 021544346, IN: 0.318, IC: (0.318, 0.318), pichap: 3.142, ICchap: (3.141, 3.143)
n: 046415888, IN: 0.318, IC: (0.318, 0.318), pichap: 3.142, ICchap: (3.141, 3.143)
n: 100000000, IN: 0.318, IC: (0.318, 0.318), pichap: 3.141, ICchap: (3.141, 3.142)Il faut de l’ordre d’un peu moins de 10^8 tirages.
Soit Z \sim N(0,1). Nous souhaitons évaluer la probabilité I = \mathbb P(Z >4.5).
Commençons avec un petit nombre de simulations
import numpy as np
# Define the number of samples
N = 10000
# Generate samples from the normal distribution with mean 0 and standard deviation 1
X = np.random.normal(0, 1, N)
# Calculate Y based on the condition (X > 4.5)
Y = (X > 4.5)
# Calculate the mean of Y
in_ = np.mean(Y)
# Calculate the confidence interval (ic) for in_
q = stats.norm.ppf(0.975)
delta = q * np.sqrt(np.var(Y) / n)
ic = (in_ - delta, in_ + delta)
print("IC:", in_)
print("Confidence Interval (IC) for IN:", ic)IC: 0.0
Confidence Interval (IC) for IN: (0.0, 0.0)On obtient une estimation à 0, ce qui n’est pas surprenant car on cherche la probabilité d’un événement rare
from scipy.stats import norm
# Compute 1 - pnorm(4.5)
result = 1 - norm.cdf(4.5)
print("1 - pnorm(4.5):", result)1 - pnorm(4.5): 3.3976731247387093e-06L’objectif de cet exercice est d’illustrer la difficulté de l’estimation de volumes en grande dimension, notamment par la méthode de Riemann.
On note B_d=\{x\in \mathbb R^d, \sum_{i=1}^d x_i^2 \leq 1\} la boule unité en dimension d et V_d son volume. Pour rappel un calcul exact donne V_d = \frac{\pi^{d/2}}{\Gamma(d/2+1)}. En notant f(\cdot) = {1\hspace{-3.8pt} 1}_{B_d}(\cdot) on peut donc écrire V_d = \int_{\mathbb{R}^d} f(x)dx.
On va dans la suite estimer V_d pour des valeurs de d allant de 2 à 10.
import numpy as np
from scipy import stats
from scipy.special import gamma
def f(x):
return np.sum(x**2) <= 1
def true_volume(d):
return np.pi**(d/2) / gamma(d/2 + 1)
def riemann(n, d):
n_discr = int(n**(1/d))
n_eff = n_discr**d
x = np.linspace(-1, 1, n_discr)
# Grille regulière sur [-1, 1]^d:
X = np.array(np.meshgrid(*[x]*d)).reshape(d, -1).T
Y = np.array([f(x) for x in X]) * 1.0
d_vol = (2 / n_discr)**d
return np.sum(Y) * d_vol, n_eff, Ximport matplotlib.pyplot as plt
n = 10000
d_max = 7
dimensions = range(2, d_max)
rel_err_mcs = []
rel_err_rs = []
for d in dimensions:
in_r, n_eff, _ = riemann(n, d)
in_mc, _, _ = montecarlo(n_eff, d)
true_vol = true_volume(d)
rel_err_mcs.append(np.abs(in_mc - true_vol) / true_vol)
rel_err_rs.append(np.abs(in_r - true_vol) / true_vol)
plt.plot(dimensions, rel_err_mcs, marker="s", label="Monte-Carlo", color='peru')
plt.plot(dimensions, rel_err_rs, marker="o", label="Riemann", color='dodgerblue')
plt.xlabel("Dimension")
plt.ylabel("Erreur relative")
plt.title(f"Erreur relative en fonction de la dimension\npour un nombre identique d'evaluations (> {n})")
plt.legend()
plt.show()
d_max=7
d_min=2
n=10000
fig, ax = plt.subplots(d_max - d_min + 1, 2, figsize=(8, 12), sharey='row')
for i, d in enumerate(range(d_min, d_max + 1)):
_, n_eff, X = riemann(n, d)
norms_r = np.sum(X**2, axis=1)
_, _, norms_mc = montecarlo(n_eff, d)
for j, norms, title, color in zip(range(2), [norms_r, norms_mc], ["Tirages réguliers", "Tirages MC"], ['dodgerblue', 'peru']):
ax[i, j].scatter(range(n_eff), norms, label=title, color=color, alpha=1, s=0.8)
ax[i, j].plot(range(n_eff), [1]*n_eff, label="1", color='red')
ax[i, j].set_ylabel("Carré de la norme ")
ax[i, 1].text(1.1, 0.5, f"d={d}", transform=ax[i, 1].transAxes, fontsize=12, verticalalignment='center')
ax[0, 0].set_title("Échantillon Riemann")
ax[0, 1].set_title("Échantillon Monte-Carlo")
ax[-1, 0].set_xlabel("Index")
ax[-1, 1].set_xlabel("Index")
plt.tight_layout()
plt.show()
Une limite des deux méthodes est la dimension. Regardons ce que donne la méthode de Riemann quand d varie.
Notons a=(a_1,\dots,a_d)^\top et b=(b_1,\dots,b_d)^\top deux points de \mathbb{R}^d. Tout d’abord supposons que donnons nous une fonction f:\mathbb{R}^d \to \mathbb{R}, et un hyper-rectangle \Omega=\prod_{j=1}^d [a_j, b_j] (avec b_j>a_j pour tout j\in\llbracket 1, d\rrbracket)sur lequel on veut calculer l’intégrale de f. Ainsi l’objectif est de calculer \int_{\Omega} f(x)dx.
Nous rajoutons l’hypothèse suivante: f est Lipschitzienne sur \Omega de constante L_{f} pour la norme \|\cdot\|_{\infty} (le max des amplitudes des coefficients), ainsi il existe L_f tel que pour
\forall (x, y) \in \Omega^2, \quad |f(x) - f(y)| \leq L_f \|x-y\|_{\infty} \enspace.
Pour cela on discrétise l’hyper-rectangle \Omega en m sous-intervalles de taille (b_i-a_i)/m pour i=1,\dots,d, ce qui forme une partition de n=m^d hyper-rectangles de \Omega: \Omega_{i_1,\dots,i_d} = \prod_{j=1}^d \left[a_j + \tfrac{i_j-1}{m}(b_j-a_j), a_j + \tfrac{i_j}{m}(b_j-a_j)\right]. et on a alors la partition: \Omega = \bigcup_{i_1=1}^m \cdots \bigcup_{i_d=1}^m \Omega_{i_1,\dots,i_d}. Pour simplifier on note
\mathrm{Vol}(\Omega_{i_1,\dots,i_d}) = \left(\tfrac{b_1-a_1}{m}\right) \cdots \left(\tfrac{b_d-a_d}{m}\right)= \frac{1}{m^d}\prod_{j=1}^d (b_{i_j}-a_{i_j}) \enspace. Pour tout i_1,\dots,i_d, on note \omega_{i_1,\dots,i_d} un point de \Omega_{i_1,\dots,i_d} (par exemple son centre, ou par simplicité le point “en bas à gauche de chaque hyper-rectangle”).
L’estimateur de l’intégrale de f sur \Omega par la méthode de Riemann est alors:
\widehat{I}^{\mathrm{R}} = \sum_{i_1=1}^m \cdots \sum_{i_d=1}^m f(\omega_{i_1,\dots,i_d}) \mathrm{Vol}(\Omega_{i_1,\dots,i_d}) \enspace, qui requiert n=m^d évaluations de f.
Notons E l’erreur de la méthode de Riemann, est donnée par E^R = |\int_{\Omega} f(x)dx - \widehat{I}^{\mathrm{R}}|, et on alors:
\begin{align*} E^R & = \left|\int_{\Omega} f(x)dx - \sum_{i_1=1}^m \cdots \sum_{i_d=1}^m f(\omega_{i_1,\dots,i_d}) \mathrm{Vol}(\Omega_{i_1,\dots,i_d}) \right| \\ & = \left|\sum_{i_1=1}^m \cdots \sum_{i_d=1}^m \int_{{\Omega}_{i_1,\dots,i_d}} \left( f(x) - f(\omega_{i_1,\dots,i_d}) \right) dx \right| \\ & \leq \sum_{i_1=1}^m \cdots \sum_{i_d=1}^m \int_{{\Omega}_{i_1,\dots,i_d}} \left| f(x) - f(\omega_{i_1,\dots,i_d}) \right| dx \\ & \leq \sum_{i_1=1}^m \cdots \sum_{i_d=1}^m \int_{{\Omega}_{i_1,\dots,i_d}} L_f \|x - \omega_{i_1,\dots,i_d}\|_{\infty} dx \\ & = L_f \sum_{i_1=1}^m \cdots \sum_{i_d=1}^m \int_{{\Omega}_{i_1,\dots,i_d}} \|x - \omega_{i_1,\dots,i_d}\|_{\infty} dx \enspace. \end{align*}
Pour x\in \Omega_{i_1,\dots,i_d}, alors \|x-\omega_{i_1,\dots,i_d}\|_{\infty} \leq \max_{j = 1,\dots, d}\tfrac{b_j-a_j}{m}=\tfrac{1}{m}\|b-a\|_{\infty}, et donc:
\begin{align*} E^R & = L_f \sum_{i_1=1}^m \cdots \sum_{i_d=1}^m \int_{{\Omega}_{i_1,\dots,i_d}} \tfrac{1}{m}\|b-a\|_{\infty} dx \\ & \leq L_f\mathrm{Vol}(\Omega) \tfrac{1}{m}\|b-a\|_{\infty} \\ & \leq \frac{L_f\mathrm{Vol}(\Omega)\|b-a\|_{\infty}}{n^{\tfrac{1}{d}}} \enspace. \end{align*}
Ainsi pour atteindre un niveau de précision \epsilon, il faut n tel que:
n \geq \left(\frac{L_f\mathrm{Vol}(\Omega)\|b-a\|_{\infty}}{\epsilon}\right)^d \enspace.
Concernant la méthode de Monte-Carlo l’estimateur est donnée par: \widehat{I}^{\mathrm{MC}} = \frac{1}{n}\sum_{i=1}^n f(X_i) \enspace, et le théorème centrale limite assure que \sqrt{n}\left(\widehat{I}^{\mathrm{MC}} - \int_{\Omega} f(x)dx\right) \overset{d}{\longrightarrow} \mathcal{N}(0, \mathrm{Var}(f(X))) \enspace. Cela donne donc une erreur de l’ordre de \frac{1}{\sqrt{n}} pour la méthode de Monte-Carlo.
On peut confirmer ces vitesses d’approximation en regardant les erreurs relatives en fonction du nombre d’évaluations pour d=4, pour la fonction indicatrice de la boule unité:
d = 4
fig, ax = plt.subplots(figsize=(7, 5))
n_list = np.logspace(2, 6, 15).astype(int)
true_vol = true_volume(d)
n_list_r = n_list**(1/d)
n_list_r = (n_list_r**d).astype(int)
rel_err_mc = []
rel_err_r = []
for n, n_r in zip(n_list, n_list_r):
in_mc, _, _ = montecarlo(n, d)
in_r, _, _ = riemann(n_r, d)
rel_err_mc.append(np.abs(in_mc - true_vol) / true_vol)
rel_err_r.append(np.abs(in_r - true_vol) / true_vol)
ax.plot(n_list, rel_err_mc, marker="s", color='peru', linestyle='')
ax.plot(n_list_r, rel_err_r, marker="s", color='dodgerblue', linestyle='')
ax.set_yscale('log')
ax.set_xscale('log')
ax.plot(n_list, rel_err_mc[0] * np.sqrt(n_list[0] / n_list), label=r"$\propto n^{-\frac{1}{2}}}$", color='peru', linestyle='--')
ax.plot(n_list_r, rel_err_r[0] * (n_list_r[0] / n_list_r)**(1/d), label=r"$\propto n^{-\frac{1}{d}}}$", color='dodgerblue', linestyle='--')
plt.title(f"Erreur relative en fonction du nombre d'évaluations pour $d={d}$")
plt.legend(["Monte-Carlo", "Riemann", r"$\propto n^{-\frac{1}{2}}$", r"$\propto n^{-\frac{1}{d}}$"], loc='lower left')
plt.show()
Notons que la vitesse d’approximation prévue par la théorie est conforme à ce que l’on observe même si la fonction indicatrice de la boule unité n’est pas Lipschitzienne.