Déviation standard et méthode Monte Carlo

[invite156cfd77]

Bonjour,

Pour "m'amuser" un peu à prendre en main la librairie numpy de python, j'ai fais un petit programme classique qui calcule Pi via une méthode Monte Carlo.
C'est un cas d'école où l'on tire aléatoirement des coordonnées x et y d'un point dans un intervalle [0;1] et tous les points dans le quart de cercle de rayon sont ceux qui vérifient x^2+y^2 < 1. A l'infini, le rapport du nombre de point dans le cercle sur le nombre de point total correspond au rapport des surfaces et vaut donc pi/4. Cette méthode est particulièrement inefficace mais ce n'était pas le but.

J'ai codé ça en faisant b paquets de n points. Chaque paquet produit une approximation de pi. Si j'affiche la moyenne de mes b valeurs de pi, j'obtiens mon approximation de pi. Mais je peux aussi calculé la déviation standard de mes paquets. Elle m'a semblé élevée, et en y réfléchissant, ce que je veux c'est que chaque paquet renforce ma confiance dans mon résultat, hors l'écart type tel que calculé ne me donne que l'écart entre mes b valeurs de pi.

J'ai donc, pour vérifié, regardé l'évolution de la déviation standard en fonction du nombre de paquet, elle ne diminue pas. J'ai donc affiché la déviation standard divisée par la racine du nombre de paquet utilisés pour chaque paquet et là, elle me semble avoir une allure plus "juste". Je pense avoir eu raison de faire cette division, mais je n'arrive pas à trouver de justification ou de doc (gratuite du moins) la dessus. Les formules généralement indiquées ne parlent que de la somme des écarts à la moyenne, sans préciser de racine de N. Est-ce que la démarche est la bonne ?

Merci

PS : Question bonus, finalement quel intérêt de faire b paquets de n point par rapport à faire n*b points tout de suite ?

Le code :

import numpy as np
from matplotlib import pyplot as plt

n = 50000
b = 1000

T = np.square(np.random.rand(n,b))+np.square(+np.random.rand(n,b))
T = np.sum(T < 1,axis=0)/n*4

print(np.mean(T),np.std(T)/np.sqrt(b))
plt.hist(T, bins = np.arange(3.1,3.18,0.001))
plt.show()

L = [np.std(T[:i])/np.sqrt(i) for i in range(1,len(T))][1:]
plt.plot(L)
plt.show()

M = [np.mean(T[:i]) for i in range(1,len(T))][1:]
plt.plot(M)
plt.show()

-----

[gg0]

Bonjour.

Je réponds d'abord à "quel intérêt de faire b paquets de n point par rapport à faire n*b points tout de suite ?" Aucun, c'est même contre-productif. Mais c'est ce que tu as fait !
Ensuite : " l'évolution de la déviation standard en fonction du nombre de paquet, elle ne diminue pas." Normal, tu répètes une expérience aléatoire, chaque expérience te donne une estimation d'un certain paramètre, après quelques expériences, il n'y a pas de raison que ça change, sauf cas très particulier.
Puis "J'ai donc affiché la déviation standard divisée par la racine du nombre de paquet utilisés pour chaque paquet et là, elle me semble avoir une allure plus "juste". Je pense avoir eu raison de faire cette division, mais je n'arrive pas à trouver de justification" ?? Pourquoi "plus juste" ? Quel a-priori as-tu dans la tête ? En tout cas, comme tu divise un nombre à peu près constant par un nombre de plus en plus grand, ça doit diminuer. Mais quel intérêt ?
Rappel : on ne fait pas des calculs "pour voir", mais parce qu'ils correspondent à quelque chose.

Cordialement.

[invite156cfd77]

Bonjour,

Merci pour votre réponse, même si elle n'est pas forcément très aimable.

Je vais me contenter de dire que j'ai trouvé les réponses ailleurs, que apparemment faire b paquets de n jets a bien un intérêt en Monte Carlo (pas pour calculer Pi).

Ensuite : " l'évolution de la déviation standard en fonction du nombre de paquet, elle ne diminue pas." Normal, tu répètes une expérience aléatoire, chaque expérience te donne une estimation d'un certain paramètre, après quelques expériences, il n'y a pas de raison que ça change, sauf cas très particulier.

Vous n'avez pas compris ma question, à chaque itération, ma confiance dans mon résultat approché devrait s'améliorer. Il se trouve, en fouillant sur internet que ca se modélise par une 1/racine(n) (dis donc je suis inspiré), voir "sample distribution of the mean, monte carlo".

"Rappel : on ne fait pas des calculs "pour voir", mais parce qu'ils correspondent à quelque chose."
Sans être en désaccord, c'est vraiment une remarque tue curiosité.

Cordialement,

[invite9dc7b526]

si on a un bon générateur de nombres aléatoires on ne gagne rien à grouper les tirages par paquets. Dans le cas contraire, regarder ce qui se passe entre les paquets peut permettre de découvrir un problème avec le générateur de nombres aléatoires.

[A voir en vidéo sur Futura]

[gg0]

Désolé de te paraître peu aimable, mais je me suis contenté de parler mathématiques (probas, en l'affaire).
Comme il ne s'agit pas d'une méthode Monte-Carlo par tirage d'échantillons dans une petite population, mais de simulation d'un ensemble de valeurs par tirage aléatoire dans une grande population (enfin, on peut l'espérer, si le générateur pseudo-aléatoire n'est pas trop faible), faire n échantillons de b valeurs n'apporte pas de gain, et on a directement une évaluation saine en prenant le n.b valeurs.
Je n'ai pas compris ta phrase "... à chaque itération, ma confiance dans mon résultat approché devrait s'améliorer. Il se trouve, en fouillant sur internet que ça se modélise par une 1/racine(n)"; en effet, tu ne dis pas d'où sort cette confiance. Si c'est seulement la taille de l'intervalle de confiance, alors comme les écarts types n'augmentent pas ni ne diminuent (en général) d'un échantillon à l'autre, c'est normal que la moyenne de ces écarts reste stable. Après, on peut trafiquer le résultat en imitant des
calculs correspondant à autre chose ... par exemple diviser par racine(n).
Alors qu'en utilisant directement les n.b valeurs, on a un intervalle de confiance bien plus petit.

Cordialement.