In [1]:
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as sps
import scipy.integrate as integrate
from time import time
In [2]:
def f(x):
return sps.norm.pdf(x)*(1+np.sin(3*x)) # comme sin(3x) est impair, l'intégrale de f vaut un.
In [3]:
# Graphe de la densité
plt.figure()
x = np.linspace(-5,5,int(1e4))
y = f(x)
plt.plot(x,y)
plt.title("Densité f")
plt.show()
In [4]:
def Rejet():
Y = np.random.randn() # Y ~ N(0,1)
gY = sps.norm.pdf(Y)
U = 2*gY*np.random.rand() # U ~ Unif([0,2*g(Y)])
while U > f(Y):
Y = np.random.randn()
U = 2*sps.norm.pdf(Y)*np.random.rand()
return Y
In [5]:
# Histogramme Rejet
Echantillon_Rejet = [Rejet() for i in range(1000)]
plt.figure()
x = np.linspace(-5,5,int(1e4))
y = f(x)
plt.plot(x,y,"r")
plt.hist(Echantillon_Rejet, bins=70, density=True, histtype="step")
plt.title("Echantillon Rejet")
plt.show()
In [6]:
# Algorithme d'échantillonnage Metropolis-Hastings
# on fait des pas ~ N(0,sigma^2) autour de l'ancienne valeur x,
# c'est à dire Q(x,y)=g((x-y)/sigma)/sigma où g est la densité d'une N(0,1)
sigma = 3 # on choisit sigma de maniere à avoir environ 25% des pas acceptés
def Metropolis_Hastings(Nsteps):
x = np.random.randn()
count = 0
for i in range(Nsteps):
y = x + sigma*np.random.randn()
u = np.random.rand()
if u < f(y)/f(x):
x = y
count += 1
return x, count
In [7]:
Metropolis_Hastings(100)
Out[7]:
In [8]:
# Histogramme MCMC
Nstep = 50
Echantillon_MH = [Metropolis_Hastings(Nstep)[0] for i in range(1000)]
plt.figure()
x = np.linspace(-5,5,int(1e4))
y = f(x)
plt.plot(x,y,'r')
plt.hist(Echantillon_MH, bins=70, density=True, histtype="step")
plt.title("Echantillon Metropolis-Hastings")
plt.show()
In [9]:
# NB : Pour approcher la densité f(x) par un histogramme,
# il vaut mieux utiliser le théorème ergodique (approche MCMC) :
def MCMCHist(Nsteps):
Ech = []
x = np.random.randn()
for i in range(Nsteps):
y = x + sigma*np.random.randn()
u = np.random.rand()
if u < f(y)/f(x):
x = y
Ech.append(x)
return Ech
plt.figure()
x = np.linspace(-5,5,int(1e4))
y = f(x)
plt.plot(x,y,'r')
plt.hist(MCMCHist(50000), bins=70, density=True, histtype="step")
plt.title("histogramme de la chaîne")
plt.show()
Question 3 : estimation du moment d'ordre 2¶
Concernant la valeur théorique, on a par parité : \begin{align*} \int x^2 (1+\sin(3 x))\,e^{-x^2/2}\frac{d x}{\sqrt{2\pi}} & =\int x^2\,e^{-x^2/2}\frac{d x}{\sqrt{2\pi}} =1 \end{align*}
In [10]:
# Approche "Rejet + Monte Carlo"
def EstRejet(n):
L = [Rejet()**2 for i in range(n)]
return(np.sum(L)/n)
In [11]:
EstRejet(10000)
Out[11]:
In [12]:
# Précision/temps : Rejet + MC
def Precision_RejetMC(Nsteps):
t = time()
Est = EstRejet(Nsteps)
Time = time()-t
Error = abs(Est-1)
print('Steps =', Nsteps,'\n','Error =',Error,'\n','Time =',Time)
In [13]:
Precision_RejetMC(10000)
In [14]:
# Approche MCMC
def EstMCMC(Nsteps):
S = 0
x = np.random.randn()
for i in range(Nsteps):
y = x + sigma*np.random.randn()
u = np.random.rand()
if u < f(y)/f(x):
x = y
S += x**2
return S/Nsteps
In [15]:
EstMCMC(1000)
Out[15]:
In [16]:
# Précision/temps : MCMC
def Precision_MCMC(Nsteps):
t = time()
Est = EstMCMC(Nsteps)
Time = time()-t
Error = abs(Est-1)
print('Steps =', Nsteps,'\n','Error =',Error,'\n','Time =',Time)
In [17]:
Precision_MCMC(20000)
In [18]:
# BONUS : Version analyse numérique "classique" (quadratures)
# [en dimension d=1, imbattable ! mais dans les choux quand d > 10]
integrate.quad(lambda x: x**2*f(x), -np.inf, np.inf)
Out[18]:
Question 4 : Même chose en plus grande dimension¶
On condidère maintenant la densité sur $\mathbb R^d$ donnée par $$ f(x_1,\ldots,x_d)=\prod_{i=1}^d (1+\sin(3x_i))\frac{e^{-x_i^2/2}}{\sqrt{2\pi}} $$ et on cherche à résoudre l'analogue de la question 3, à savoir évaluer l'intégrale $$ \int_{\mathbb R^d} \|x\|^2 f(x) dx. $$ Sa valeur théorique est donnée par $$ \int_{\mathbb R^d} \|x\|^2 f(x) dx = \sum_{i=1}^d \int_{\mathbb R^d}x_i^2 f(x) dx = \sum_{i=1}^d \int_{\mathbb R}x_i^2(1+\sin(3 x_i) dx_i=d. $$
In [19]:
d = 10
In [20]:
def f_D(x): return(np.prod(f(x)))
In [21]:
# Methode Rejet : dimension quelconque
def Rejet_D():
Y = np.random.multivariate_normal(np.zeros(d),np.identity(d)) # Y ~ g = N(0,I_d)
gY = np.prod(sps.norm.pdf(Y))
U = (2**d)*gY*np.random.rand() # U ~ Unif([0,2^d*g(Y)])
while U > f_D(Y):
Y = np.random.multivariate_normal(np.zeros(d),np.identity(d))
U = (2**d)*gY*np.random.rand()
return Y
def EstRejet_D(n):
L = [sum(Rejet_D()**2) for i in range(n)]
return(sum(L)/n)
def Precision_RejetMC_D(Nsteps):
t = time()
Est = EstRejet_D(Nsteps)
Time = time()-t
Error = abs(Est-d)
print('Steps =', Nsteps,'\n','Error =',Error,'\n','Time =',Time)
In [22]:
Precision_RejetMC_D(100) # C'est long (et mauvais), c'est normal : Proba de rejet = 1 - 1/2^d !
In [23]:
# Methode MCMC : dimension quelconque
sigma = 0.3 # on a encore choisi sigma de manière à avoir environ 25% de pas acceptés (d=10)
def MCMC_D(Nsteps):
x = np.random.multivariate_normal(np.zeros(d),np.identity(d))
count = 0
for i in range(Nsteps):
y = x + sigma*np.random.multivariate_normal(np.zeros(d),np.identity(d))
u = np.random.rand()
if u < f_D(y)/f_D(x):
x = y
count += 1
return x, count
In [24]:
MCMC_D(100)[1]
Out[24]:
In [25]:
def EstMCMC_D(n):
S = 0
x = np.random.multivariate_normal(np.zeros(d),np.identity(d))
for i in range(n):
y = x + sigma*np.random.multivariate_normal(np.zeros(d),np.identity(d))
u = np.random.rand()
if u < f_D(y)/f_D(x):
x = y
S += sum(x**2)
return S/n
def Precision_MCMC_D(Nsteps):
t = time()
Est = EstMCMC_D(Nsteps)
Time = time()-t
Error = abs(Est-d)
print('Steps =', Nsteps,'\n','Error =',Error,'\n','Time =',Time)
In [26]:
# MCMC >> Rejet en dimension >10 !
Precision_MCMC_D(10000)