Python モンテカルロシミュレーションにより円周率π=3.14を導く

Python 数値計算 グラフ化 散布図

　モンテカルロ法とは乱数を用いて行うシミュレーションの手法のことです。本記事の内容は、ググったら出てくるような内容をPythonを用いて実行できる雛形コードを載せました。本コードの特徴は次の2つです。
1. 乱数を一様乱数の場合と、正規乱数の場合で実行できる仕様にしました
2. 計算過程をmatplotlibのグラフ描画し、その動画をmp4ファイルに保存する仕様にしました。特に、matplotlibのグラフテクニックとして、軸目盛の間隔や表記角度などを盛り込んでいます。

▼本コードの実施例は下記リンク先です （動画編集には別ソフトを使用しています)

www.youtube.com

■本プログラム

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# https://unit.aist.go.jp/diversity/ja/jst/teens/montecarlo.htm
import math
import random
import matplotlib.pyplot as plt
import matplotlib.patches as patches
import matplotlib.ticker as ticker
import matplotlib.animation as animation

# 一様乱数を生成する関数
def random_uniform_func(n):
    min = -1 * my_radius
    max = my_radius
    a = random.uniform(min, max)
    b = random.uniform(min, max)

    return (a, b)

# 正規乱数を生成する関数
def random_normalvariate_func(n):
    mu = 0 # 分布の中心をゼロ
    sigma = my_radius/3 # 3σが半径とした場合の乱数
    a = random.normalvariate(mu, sigma)
    b = random.normalvariate(mu, sigma)

    return (a, b)

# アニメーションのためのアップデート関数
def update_ani(i, X, Y, pai_list):
    # plt.cla()
    if (i % 20) == 0:
        print(i, '/', N)

    ax.set_title(str(i) + '/' + str(N) + '   ' + 'pai=' + str(pai_list[i]), fontsize = my_fontsize)
    if math.sqrt(X[i]*X[i] + Y[i]*Y[i]) <= my_radius:
        ax.scatter(X[i], Y[i], s=80, color='b', marker='.')
    else:
        ax.scatter(X[i], Y[i], s=80, color='r', marker='.')

    #return img,

# メイン関数
def main():
    global ax

    # 乱数の自作関数をリストで設定。ひとつだけでなくて、複数を順番に実行するため。
    my_random_func_list = [random_uniform_func, random_normalvariate_func]
    for k, my_random_func in enumerate(my_random_func_list):
        ############## グラフの設定、32行
        fig, ax = plt.subplots()
        ax.scatter([], [])
        myCircle = patches.Circle(xy=(0, 0), radius=my_radius, ec='k', fill=False) # fc='b',
        ax.add_patch(myCircle)

        tick_spacing = my_radius/2 # 目盛り表示する間隔
        
        ### X軸表記の設定
        ax.set_xlabel("x", fontsize = my_fontsize)
        # レンジの範囲
        ax.set_xlim(-1*my_radius, my_radius)
        # デフォルトの目盛り表記をゲットする
        x_ticklabels = ax.get_xticklabels()
        # 目盛り表記を90度回転。#フォントサイズの指定する場合
        plt.setp(x_ticklabels, rotation=45, fontsize=my_fontsize)
        # X軸目盛の表示間隔を間引く
        ax.xaxis.set_major_locator(ticker.MultipleLocator(tick_spacing))

        ### Y軸表記の設定
        ax.set_ylabel("y", fontsize = my_fontsize)
        ax.set_ylim(-1*my_radius, my_radius)
        y_ticklabels = ax.get_yticklabels()
        plt.setp(y_ticklabels, rotation=0, fontsize=my_fontsize)
        ax.yaxis.set_major_locator(ticker.MultipleLocator(tick_spacing))

        # 格子
        ax.grid()
        # 縦横比
        ax.set_aspect('equal')
        # ラベルレイアウトを調整する。文字がはみ出ないようにする。
        plt.tight_layout()
        ##############

        cnt = 0 # 円の中にプロットされた数を数えるための変数の初期化
        X = [] # プロット点のX座標を乱数でリストで設定するための変数の初期化
        Y = [] # プロット点のY座標を乱数でリストで設定するための変数の初期化
        pai_list = [] # 各時点での計算結果πを格納するリストの初期化
        # 乱数を生成して、X, Yリストへ格納
        for n in range(N):
            x, y = my_random_func(n)
            if math.sqrt(x*x + y*y) <= my_radius:
                cnt = cnt + 1 # cntの数は円の中にプロットされた数
            X.append(x)
            Y.append(y)
            # モンテカルロ法（乱数を使用したシミュレーション）によるπの計算結果
            P = cnt / (n+1) # cntの数は円の中にプロットされた数
            calc_pai = 4 * P
            pai_list.append(round(calc_pai, 2)) # 小数点2桁でリストへ格納
        print("pai:", calc_pai)

        # アニメーションの実行と、それをオブジェクト化
        ani = animation.FuncAnimation(fig,
                                      update_ani,
                                      frames = N, #アップデート関数を呼ぶ回数
                                      interval = 200,
                                      blit = False,
                                      repeat = False,
                                      fargs = (X,Y, pai_list)) #アップデート引数を二つ以上使用する場合

        # 動画ファイルに保存
        ani.save('anim_' + str(k) + '_N' + str(N) + '_pai' + str(calc_pai) + '.mp4', writer='ffmpeg', fps=30)
        #ani.save('anim_' + str(k) + '_N' + str(N) + '_pai' + str(calc_pai) + 'gif', writer="imagemagick", fps=30)
        #plt.show()
        plt.close("all")

if __name__ == '__main__':
    my_radius = 10
    my_diameter = my_radius * 2
    my_fontsize = 15

    N = 100
    random.seed(200727)
    
    main()
    print("finished")

以上

<広告>