Article épinglé

vendredi 23 décembre 2022

Méthode de Monte-Carlo, pi et python

Comme il en est question dans le livre "Lune rouge" de Robinson, je creuse et je trouve:

https://fr.wikipedia.org/wiki/M%C3%A9thode_de_Monte-Carlo 

Présentation 

Cette méthode est proche de l'expérience de l'aiguille de Buffon.

Soit un point M de coordonnées (x, y), où 0<x<1 et 0<y<1. On tire aléatoirement les valeurs de x et y entre 0 et 1 suivant une loi uniforme. Le point M appartient au disque de centre (0,0) de rayon R=1 si et seulement si x2+y2≤1. La probabilité que le point M appartienne au disque est π/4, puisque le quart de disque est de surface σ=π R2/4=π/4, et le carré qui le contient est de surface S=R2=1 : si la loi de probabilité du tirage de point est uniforme, la probabilité de tomber dans le quart de disque vaut σ/S=π/4.
En faisant le rapport du nombre de points dans le disque au nombre de tirages, on obtient une approximation du nombre π/4 si le nombre de tirages est grand.


Script python

Rq
Résultat du script du dessous.
Vous noterez avec raison que les boucles sont redondantes et que l'usage du dictionnaire n'est pas obligatoire. Donc vous optimiserez le script ;)
Notez les commentaires dans les premières lignes du script. Je vois deux problème:
  • Dans la méthode x et y sont strictement sup à 0 et inf à 1 alors que la fonction random() génère des nombres entre 0 et 1 inclus.
  • L'autre problème concerne encore cette fonction, elle génère des nombres pseudo-aléatoires et non aléatoires. 



---- script en dessous que vous devriez optimiser

# script découpé pour comprendre les articulations de la demarche
# evidemment, une boucle devrait suffire

import random as rd
tirage=5000

point={}#dictionnaire

for i in range (tirage):# création des couples x et y et enregistrement dans une bibliothèque
    # pas fait - verifier les limites de random car x et y strict compris entre 0 et 1
    # a noter pb nbre aleat avec random car pseudo aleatoire
    x=rd.random()
    y=rd.random()
    point[i]=(x,y)

# pour comprendre le fonctionnement du dictionnaire, décommenter les 3 lignes du dessous  
# print(point)
# print("Point0=",point[0])
# print("X0=",point[0][0],"Y0=",point[0][1])

# Ncercle: comptage des points qui sont DANS le cercle donc X²+Y² <=1
Ncercle=0
for i in range(tirage):
    c=point[i][0]**2+point[i][1]**2
    if c<=1:
        Ncercle=Ncercle+1

# calcul de l'approximation de  = 4*Ncercle/tirage
print ("Pi approx=",4*Ncercle/tirage)

# mise en place graphique avec matplotlib
# creation des points
import matplotlib.pyplot as plt
for i in range (tirage):
    c=point[i][0]**2+point[i][1]**2
    if c<=1:
        couleur='r'
    else:
        couleur='b'
    plt.scatter(point[i][0],point[i][1],color=couleur)

# creation du quart de cercle
import numpy as np
x = np.linspace(0,1, 100)
y = np.sqrt(1-x**2)

plt.plot(x, y)

plt.show()




Aucun commentaire:

Enregistrer un commentaire

Tout commentaire nous engage ;)