モンテカルロ法とは乱数を用いて行うシミュレーションの手法のことです。本記事の内容は、ググったら出てくるような内容をPythonを用いて実行できる雛形コードを載せました。本コードの特徴は次の2つです。
1. 乱数を一様乱数の場合と、正規乱数の場合で実行できる仕様にしました
2. 計算過程をmatplotlibのグラフ描画し、その動画をmp4ファイルに保存する仕様にしました。特に、matplotlibのグラフテクニックとして、軸目盛の間隔や表記角度などを盛り込んでいます。
▼本コードの実施例は下記リンク先です （動画編集には別ソフトを使用しています)
■本プログラム
#!/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")
以上
<広告>
リンク
