[PYTHON] Saupoudrer de grains de riz pour trouver le rapport de circonférence (méthode de Monte Carlo)

Postscript

Il s'agit d'un article qui a tenté d'implémenter la méthode Monte Carlo. J'ai complètement oublié.

1.Tout d'abord

Quand je lisais un livre de statistiques pour enfants [1], j'ai vu un article sur l'aspersion de grains de riz et le calcul du rapport de circonférence. Cela semblait intéressant, j'ai donc calculé le rapport de circonférence en utilisant les données de coordonnées x et y générées aléatoirement.

[procédure] C'est presque comme ça

  1. Préparez un carré et un cercle inscrits dans ce carré.
  2. Tracez-y au hasard (au hasard). Le fait est que cette ** merde **!
  3. Obtenez le rapport de circonférence à partir du nombre de coordonnées à l'intérieur et à l'extérieur du cercle. image.png [Objectif] Voyons combien de coordonnées (échantillon, c'est-à-dire, grains de riz) sont nécessaires pour calculer une valeur proche du rapport de circonférence réel ou une valeur proche du rapport de circonférence.

2. Mettre en œuvre

Il était difficile de semer des grains de riz et d'expérimenter, alors je l'ai essayé sur un ordinateur.

2-1. Modules d'importation.

import numpy as np
import random
import matplotlib.pyplot as plt
import matplotlib.patches as patches

Ensuite, définissez des listes et des constantes. Cette fois, nous utiliserons un carré avec un côté et un cercle avec un rayon de 0,5.

Le but de ce temps est d'étudier combien de grains de riz devraient être nécessaires pour calculer une valeur proche du vrai rapport de circonférence, donc un exemple du nombre d'échantillons a été provisoirement défini dans n_list.

#Nombre de coordonnées générées
n_list = [10, 20, 50, 70, 100, 500, 1000, 5000, 10000, 100000,1000000]

#Rayon du cercle
radius = 0.5

#Initialisation de la liste des taux de circonférence
pi_list = []

2-2. Générer des coordonnées aléatoires

Définit une fonction qui génère des données de coordonnées aléatoires (correspondant à un grain de riz enroulé autour d'un désordre).

def gen_random_coordinate(n):
    #Initialiser la liste de coordonnées générée
    x=[]
    y=[]

    #Générer des coordonnées avec des nombres aléatoires
    x = np.random.rand(n)
    y = np.random.rand(n)
        
    return x,y

2-3. Visualisez

Visualisez comment vous avez réellement tracé des coordonnées aléatoires.

#Fonction de dessin graphique
def make_graph(x,y):
    
    #Créer un objet figure et un objet Axes lui appartenant en même temps
    fig, ax = plt.subplots()
    #Créer un cercle
    circle = patches.Circle(xy=(0.5,0.5), radius=0.5, ec='r', fill=False)

    #Ajuster au bon rapport hauteur / largeur
    ax.set_aspect('equal')
    ax.add_patch(circle)
    ax.scatter(x,y)
    
    plt.axis('scaled')
    plt.grid()
    plt.xlim(0, 1.0)
    plt.ylim(0, 1.0)

    plt.show()

2-4. Calculer le rapport de circonférence à partir des données de coordonnées

Calculez le rapport de circonférence en utilisant le nombre de points à l'intérieur et le nombre de points à l'extérieur du cercle.

def calc_circular_constant(x,y):
    #Initialisation du nombre de points à l'intérieur du cercle
    inner=0
    
    #Initialisation du nombre de points en dehors du cercle
    outer=0
    
    #Comptez le nombre de points, que chaque point soit à l'intérieur ou à l'extérieur du cercle
    for xx , yy in zip(x,y):
        
        #Calculez la distance du centre du cercle
        dist=np.sqrt(np.square(xx-radius)+np.square(yy-radius))
        
        if dist > radius:
            outer +=1
        else:
            inner +=1
            
    print('inner:outer = ', str(inner), ':', str(outer))
    
    #Calculer le rapport de circonférence
    pi = inner /(outer+inner) *4
    pi_list.append(pi)
     
    return pi

2-5. Exécuter

Exécutez chacune des fonctions définies ci-dessus. Vous avez vérifié le nombre d'échantillons dont vous devriez disposer pour vous rapprocher du rapport de circonférence réel. Je vais donc les essayer dans l'ordre en utilisant une liste d'échantillons n_list.

for n in n_list:
    x,y = gen_random_coordinate(n)
    make_graph(x,y)
    print(calc_circular_constant(x,y))

3. Résultats expérimentaux

Voici le résultat de la sortie. Le nombre d'échantillons commence à 10 et le nombre de points augmente à mesure que vous descendez.

image.png Point intérieur du cercle: Point extérieur du cercle = 8: 2 Rapport de circonférence calculé: 3,2

image.png Point intérieur du cercle: Point extérieur du cercle = 16: 4 Rapport de circonférence calculé: 3,2

image.png Point intérieur du cercle: Point extérieur du cercle = 39: 11 Rapport de circonférence calculé: 3,12

image.png Point intérieur du cercle: Point extérieur du cercle = 53:17 Rapport de circonférence calculé: 3,0285714285714285

image.png Point intérieur du cercle: Point extérieur du cercle = 79:21 Rapport de circonférence calculé: 3,16

image.png Point intérieur du cercle: Point extérieur du cercle = 399: 101 Rapport de circonférence calculé: 3,192

image.png Point intérieur du cercle: Point extérieur du cercle = 793: 207 Rapport de circonférence calculé: 3,172

image.png Points à l'intérieur du cercle: Points à l'extérieur du cercle = 3894: 1106 Rapport de circonférence calculé: 3.1152

image.png Points à l'intérieur du cercle: Points à l'extérieur du cercle = 7818: 2182 Rapport de circonférence calculé: 3,1272

image.png Points à l'intérieur du cercle: Points à l'extérieur du cercle = 78611: 21389 Rapport de circonférence calculé: 3,14444

image.png Points à l'intérieur du cercle: Points à l'extérieur du cercle = 785814: 214186 Rapport de circonférence calculé: 3,143256

4. Résumé

Dans cette expérience, il a été confirmé qu'à mesure que le nombre d'échantillons augmentait, il s'approchait du rapport de circonférence réel. Le graphique devient bleu foncé lorsque le nombre de données dépasse 5 000. Au final, il était de 3,143256, ce qui était proche du vrai ratio de circonférence (3,141592653589 ..).

Si le nombre de données est d'environ 100, il est encore loin de 3,16. (C'est vrai, mais j'ai écrit la raison en 5.)

Les données à générer doivent également être des données uniformément aléatoires. Avec la méthode ci-dessus, si vous générez des données qui suivent une distribution normale, vous ne pourrez pas calculer un rapport de circonférence décent. (Je me demande si je peux corriger l'extérieur en utilisant le nombre de données à l'intérieur du cercle, je vais l'essayer à un autre moment.)

5. Enfin

Même si je l'écris maintenant, peu importe à quel point cela va avec 100, il sera calculé en 0,04 pas, donc après 3,12, ce sera 3,16. (Je n'ai pas vraiment abordé ce point car je voulais faire l'expérience ci-dessus.)

Si vous voulez obtenir des résultats comme le ratio de circonférence que vous apprenez à l'école primaire, vous avez besoin de 400 données. (Parce que le pas du rapport de circonférence à calculer est de 0,01). Si vous le faites encore et encore, le nombre 3.14 est susceptible de sortir.

Même si vous ne comptez les grains de riz qu'à l'extérieur du cercle, il semble être une expérience difficile de saupoudrer 400 grains de riz et d'essayer plusieurs fois. Étant donné que la taille du grain de riz lui-même est également grande, il semble que vous deviez préparer un papier de la taille appropriée.

Si vous avez des enfants, que diriez-vous d'une étude gratuite des «vacances de printemps» soudainement longues?

Les références

[1] À quoi servent les statistiques? (Auteur) Yoshiyuki Wakui (Edit) Département de rédaction scientifique pour enfants / Seibundo Shinkosha https://www.seibundo-shinkosha.net/book/kids/20689/

Recommended Posts

Saupoudrer de grains de riz pour trouver le rapport de circonférence (méthode de Monte Carlo)
Méthode #Monte Carlo pour trouver le rapport de circonférence en utilisant Python
Introduction à la méthode Monte Carlo
Trouvez le ratio de la superficie du lac Biwa par la méthode de Monte Carlo
Méthode de Monte Carlo
Essayez d'implémenter la méthode Monte Carlo en Python
La première méthode de Monte Carlo en chaîne de Markov par PyStan
Calcul de l'itinéraire le plus court selon la méthode de Monte Carlo
Augmentez la vitesse de la méthode Monte Carlo de l'implémentation de découpage Cython.
[Statistiques] Visualisez et comprenez la méthode Hamiltonian Monte Carlo avec animation.
Simuler la méthode Monte Carlo en Python
Estimation de π par la méthode de Monte Carlo