Incompréhension de l'estimation de Monte Carlo Pi

9

Je suis assez sûr de comprendre comment fonctionne l'intégration de Monte-Carlo, mais je ne comprends pas la formulation de la façon dont elle est utilisée pour estimer Pi. Je vais suivre la procédure décrite dans la 5ème diapositive de cette présentation http://homepages.inf.ed.ac.uk/imurray2/teaching/09mlss/slides.pdf

Je comprends les étapes préliminaires. Pi est égal à 4 fois l'aire d'un quart du cercle unitaire. Et l'aire du quart supérieur droit du cercle unité centrée sur (0,0) est équivalente à l'intégrale de la courbe qui est le quart supérieur droit du cercle unité dans 0<x<1 et 0<y<1 .

Ce que je ne comprends pas, c'est comment cette intégrale est

I((x2+y2)<1)P(x,y)dxdy

où est uniformément distribué dans le carré unitaire autour du quart de cercle (c'est-à-dire qu'il est toujours égal à 1 si et et 0 sinon). Cela signifie donc que est la fonction qui est le quadrant supérieur droit du cercle unitaire à et mais je ne comprends pas comment cela est vrai car la fonction d'indicateur ne peut être que 1 ou 0. Je comprends qu'il est probablement écrit de cette façon pour faciliter l'échantillonnage de Monte Carlo (c'est-à-dire qu'il s'agit d'une attente, il suffit donc d'échantillonner à partir de et obtenir la moyenne des échantillons appliqués àP(x,y)0<x<10<y<1I((x2+y2)<1)P(x,y)
0<x<10<y<1P(x,y)I((x2+y2)<1)), mais cela n'a tout simplement pas de sens intuitif pour moi pourquoi cette intégrale représente l'aire sous cette courbe.

Quelqu'un pourrait-il fournir une explication intuitive de cela. Peut-être montrer comment cette intégrale a été dérivée étape par étape?

ÉDITER:

J'ai pu acquérir une meilleure compréhension en reliant l'attente à un domaine. Je vais l'expliquer ici au cas où cela aiderait quelqu'un. Commencez d'abord par lier Pi à l'aire du quadrant supérieur droit du cercle unitaire

π=4×Atr

Ensuite, nous plaçons le quadrant supérieur droit dans le carré de l'unité. Et sous une distribution uniforme sur le carré unitaire, l'aire du quadrant circulaire est proportionnelle à la probabilité d'en obtenir un échantillon. Il s'ensuit que l'égalité suivante s'applique

P(x2+y2<1)=AtrAsquare

et doncAsquare=1

P(x2+y2<1)=Atr

Et en remplaçant l'équation originale

π=4×P(x2+y2<1)

et il est également vrai que qui est égale à la double intégrale d'origine.P(x2+y2<1)=E[I(x2+y2<1)]

Je l'ai donc compris en reliant l'aire à une probabilité puis en reliant cette probabilité à une attente équivalente à l'intégrale. Faites-moi savoir si j'ai fait des erreurs.

user1893354
la source

Réponses:

8

L'aire d'un cercle cercle de rayon est égale à . Cela signifie qu'un quart de cercle a l'aire . Cela signifie que le carré avec le côté du rayon du cercle comme .lπl2l2π/4area=l2

Cela signifie que le rapport entre l'aire d'un quart de cercle et l'aire du carré est . π/4

Un point est dans le carré si . et il est dans le quart de cercle si . (x,y)0<x<1,0<y<10<x<1,0<y<1,x2+y2<1

Votre intégrale est donc C'est exactement la zone décrite par un quart de cercleI((x2+y2)<1)P(x,y)=I((x2+y2)<1)I(0<x<1)I(0<y<1)

entrez la description de l'image ici

Donbeo
la source
Je suppose que j'ai juste du mal à établir une connexion entre les termes à l'intérieur de l'intégrale et la courbe elle-même. Si vous représentiez I (x ^ 2 + y ^ 2 <1) I (0 <x <1) (0 <y <1) pour différentes valeurs de x et y, vous n'obtiendrez pas la courbe. Pourquoi donc?
user1893354
1
{(x,y):(x2+y2<1),(0<x<1),(0<y<1)} sont les points sur le quart de cercle. Je vous suggère d'essayer de tracer ces points
Donbeo
Je suis d'accord avec ça. Mais lorsque vous appliquez la fonction d'indicateur I (.), Ils sont tous poussés à 1 ou à 0.
user1893354
Que voulez-vous dire?
Donbeo
1
La fonction d'indicateur dans une intégrale est juste une autre façon de définir les courbes où calculer l'intégrale. quarter of circle=1(x2+y2<1)1(0<x<1)1(0<y<1)
Donbeo
4

L'explication intuitive la plus simple repose sur la compréhension que . Ainsi, . Une fois que vous réalisez que le double entier est simplement une probabilité, il devrait être intuitif que vous puissiez échantillonner et partir du carré unitaire et calculer la proportion de tirages pour lesquels . E(I(A))=P(A)I(x2+y2<1)dxdy=P(x2+y2<1)xyx2+y2<1

Peut-être que l'autre intuition qui manque à votre compréhension est le lien entre la zone et la probabilité. Etant donné que la zone de l'ensemble du carré unitaire est de 1 et des points sont uniformément répartis à l'intérieur de la place, la surface de toute la région l' intérieur du carré unité correspondrait à la probabilité qu'un point choisi au hasard serait à l'intérieur de .(x,y)AA

jsk
la source
C'est aussi ainsi que je le comprends. Mais j'ai du mal à le connecter à la formulation Pi = 4x (zone du quart de cercle). Il n'est pas vraiment intuitif de comparer des zones à des échantillons. Je suppose que le lien est que sous une distribution uniforme, le nombre d'échantillons est proportionnel à la zone.
user1893354
1
@ user1893354 Réponse révisée. Faites-moi savoir si cela aide votre intuition.
jsk
0

J'ai atterri sur ce CV de surf, et je vois que le code du Monte Carlo est en Octave. Il se trouve que j'ai une simulation dans R qui rend l'idée de dériver le nombre comme une distribution uniforme bivariée dans le plan sous les contraintes des intégrales dans l'OP très intuitive:π[0,1]

Étant donné que le quart de cercle est enfermé dans un carré d'une unité, l'aire est . Donc, générer des points uniformément répartis dans le carré finira par recouvrir le carré entier, et calculer la fraction remplissant équivaudra à intégrer car nous ne faisons que sélectionner la fraction de points dans le cercle par rapport au carré unitaire:π/4(x,y)1<(x2+y2)1((x2+y2)<1)1(0<x<1)1(0<y<1)

x <- runif(1e4); y <- runif(1e4)
radius <- sqrt(x^2 + y^2)
# Selecting those values within the circle is obtained with radius[radius < 1]:
(pi = length(radius[radius < 1]) / length(radius)) * 4     =    3.1272

Nous pouvons tracer les valeurs comprises dans le rayon parmi 10 000 tirages:

entrez la description de l'image ici

Et nous pouvons, naturellement, nous rapprocher de plus en plus en sélectionnant plus de points. Avec 1 million de points, nous obtenons:

(pi = length(radius[radius < 1]) / length(radius)) * 4 [1] 3.141644

un résultat très approximatif. Voici l'intrigue:

entrez la description de l'image ici

Antoni Parellada
la source