In [1]:
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import *

1 - Méthode de rejet

Illustration du temps d'attente pour tirer de la mesure uniforme sur la boule $B_d$ de rayon un en dimension $d$, $B_d=\{x\in\mathbb R^d : \; \|x\|_2=1\}$, partant de la distribution uniforme sur l'hypercube $[-1,1]^d$.

In [2]:
def M(d):
    return(2**d*gamma(d/2+1)/(np.pi)**(d/2))
In [3]:
n = 30
x = [100/M(i) for i in range(n)]

plt.plot(x)
plt.xlabel('Dimension')
plt.ylabel("Probabilité d'accepter (%)")
plt.show()
In [4]:
1/M(10)
Out[4]:
0.00249039457019272
In [5]:
def RejetSphere(d,Tstop):
    tentative = 0
    proposition = 2*np.random.uniform(size=d)-1  # uniforme [-1,1]^d
    while (sum(proposition**2) > 1 and tentative < Tstop): 
        tentative += 1
        proposition = 2*np.random.uniform(size=d)-1
    if (sum(proposition**2) <= 1):
        print(proposition[1:7],"\n","nombre de tentatives =",tentative)
    else:
        print("Pas accepté","\n","nombre de tentatives =",Tstop)
In [6]:
RejetSphere(10,1e6)

[ 0.00235523  0.51764811 -0.59333443 -0.1313208   0.22281556  0.08204258] 
 nombre de tentatives = 184
In [7]:
RejetSphere(20,1e6)

Pas accepté 
 nombre de tentatives = 1000000.0