きっかけ
特にこれといった理由は無いのですが、ネットサーフィンをしていたらモンテカルロ法のWikipediaにたどりつき、以下の画像を目にしました。
 |
| Wikipediaのモンテカルロ法のページより引用 |
いつものことですが間違っていたり、改善点がありましたらご教示お願いいたします。
モンテカルロ法とは
モンテカルロ法と円周率の近似計算によるとモンテカルロ法とは乱数を使って何かしらの値を推定する方法のことのようです。
詳しくは後程説明しますが、そのモンテカルロ法を用いた推定の中で最も有名なのが上記図のGifのように円の内側に入った点の数を求めて円周率を推定する方法です。
今回はモンテカルロ法を使って円周率を求める方法について説明します。
上図のようにN個の点を一辺がrの正方形(上図、青＋緑部分)の中に打っていき、その中で4分の1の円の中(上図青の部分)に入る点の数Sを数えます。
Nの数が十分に大きいと、Nは正方形の面積に対応し、Sは4分の1の円の面積に対応するようになります。
Nの数が十分に大きいと、Nは正方形の面積に対応し、Sは4分の1の円の面積に対応するようになります。
4分の1の円の面積は$\frac{\pi r^2}{4}$で表され、正方形の面積は$r^2$で表されるので、点Sと点Nの個数の比は次のようになります。
\begin{equation}
\frac{S}{N} = \frac{\pi r^2/4}{r^2} \nonumber \\
=\frac{\pi^2}{4}
\end{equation}
よって点Sと点Nの比を求めて、4をかけたものが円周率の近似値として表すことが出来ます。
Gifを作る前に静止画でモンテカルロ法で円周率を求める
今回の目的はGif形式の画像を作ることですが、その前に静止画でモンテカルロ法を使って円周率を求めてみます。
ソースコードは以下の通りです。
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
def judge_in_circle(x,y):
return x**2+y**2 <= 1#円の内側かどうかの判定
if __name__ == '__main__':
#初期化・パラメータ設定
#--------------------------------------------------------------------------------------------------------------------
Num = 10**6 #打つ点の数
Z = []#円の内側か、外側かを格納する配列(1が内側、０が外側)
X,Y = np.random.rand(2,Num)#0以上1未満の乱数を返す。2は次元数(今回はX,Y 2つあるので)
in_circle_count = 0 #円の内側に入った点の数を数える
fig = plt.figure(figsize = (9,9))#図の描画
#--------------------------------------------------------------------------------------------------------------------
#実際の処理
#--------------------------------------------------------------------------------------------------------------------
for x,y in zip(X,Y):
if judge_in_circle(x,y):
Z.append(float(1)) #色を指定する際にfloat型でないとダメなので、float型に変換し格納
in_circle_count += 1
else:
Z.append(float(0))
pi = 4*in_circle_count/Num
#--------------------------------------------------------------------------------------------------------------------
#画像出力
#--------------------------------------------------------------------------------------------------------------------
plt.scatter(X,Y, s = 0.5, c = Z, cmap = "seismic") #描画
plt.savefig("monte_carlo.png")
plt.show()
print(pi)
#--------------------------------------------------------------------------------------------------------------------
出力された画像は以下のようになりました。
いよいよGifを作成する。
事前準備(ImageMagickのインストール)
Gif形式で保存するためにはImageMagickというソフトをインストールする必要があるみたいです。ImageMagickのインストールについては以下のサイトを参考にしました。
今回Gif形式にして保存する際に、matplotlibのanimation.FuncAnimationを利用しました。
matplotlib.animation.FuncAnimation(fig, func, frames=None, init_func=None, fargs=None, save_count=None, **kwargs)参考文献によると、引数に利用するupdate関数を利用するのに以下のような注意点があるそうです。
- Funcanimation関数より前に定義しておく
- 2回目以降にこの関数が呼び出されるときに、図が初期化されるようにしておく。
- plt.cla()を用いる
- 今回は前の図に上書きしていくようなイメージでGifを作成したため、初期化処理は行っていない
- 少なくとも1つ引数を設ける
- update関数の第一引数はこの関数が呼び出される毎に0から1ずつ増加していく。
- update関数に2つ以上の引数を設けたい場合にはfargsで指定する。
ソースコード
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
def judge_in_circle(x,y):
return x**2+y**2 <= 1#円の内側かどうかの判定
def plot(i,X,Y,Z,circle_count_list):#描画用の関数
# plt.cla() #1枚ずつ作るなら必要そうだけど、今回は前回の画像にどんどん追加していくイメージ？
if i != 0:
plt.title(f"n={i*frame_point_num} π ≒ {(4*circle_count_list[i]/((i+1)*frame_point_num)):.5}")
plt.scatter(X[(i-1)*frame_point_num:i*frame_point_num],Y[(i-1)*frame_point_num:i*frame_point_num],c = Z[(i-1)*frame_point_num:i*frame_point_num],s= 0.5,cmap="seismic") #2500点毎に打っていく
if __name__ == '__main__':
#初期化・パラメータ設定
#--------------------------------------------------------------------------------------------------------------------
total_point_num = 3*10**5 #打つ点の数
frame_point_num = 2500 #1フレーム毎に打つ点の数
interval_num = 10**(2) #アニメーションで図が切り替わる時間間隔の指定
X,Y = np.random.rand(2,total_point_num)
in_circle_count = 0
Z = []
circle_count_list = [] #frame_point_num毎に円の内側に入った点の数を格納する配列
fig = plt.figure(figsize = (9,9))
#--------------------------------------------------------------------------------------------------------------------
#実際の処理
#--------------------------------------------------------------------------------------------------------------------
for i,(x,y) in enumerate(zip(X,Y)):
if judge_in_circle(x,y):
Z.append(float(1)) #色を指定する際にfloat型でないといけないので、変換し配列に格納
in_circle_count += 1
else:
Z.append(float(0))
if i%frame_point_num == 0 and i != 0: #frame_point_num毎に円の内側に入った点の数を格納。0の時は意味が無いので、例外処理
circle_count_list.append(in_circle_count)
#--------------------------------------------------------------------------------------------------------------------
#画像の出力
#--------------------------------------------------------------------------------------------------------------------
print("画像の出力を始めます")
ani = animation.FuncAnimation(fig, plot,interval = interval_num,frames = int((total_point_num/frame_point_num)-1),fargs=(X,Y,Z,circle_count_list))
ani.save('anim_tmp.gif', writer="imagemagick")
print("画像の出力が終わりました")
#--------------------------------------------------------------------------------------------------------------------
これを実行すると次のようなGifが出来ます。
当初の目的としていたGifを作ることが出来たので、とりあえず満足しました。
苦戦したところ
12行目の
plt.scatter(X[(i-1)*point_num:i*point_num],Y[(i-1)*point_num:i*point_num],c = Z[(i-1)*point_num:i*point_num],s= 0.5,cmap="seismic")
ですが、最初は以下のように1フレーム毎に画像を初期化して描画していました。
plt.cla() #場所は上記コード参照
plt.scatter(X[(:i*point_num],Y:i*point_num],c = Z[:i*point_num],s= 0.5,cmap="seismic")
このようにすると、打点数が増えるごとに処理に時間が掛かってしまうため、打点数が$10^6$でも相当時間が掛かってしまいました。
画像を初期化せずに点をどんどん追加していくことで、時間を短縮することが出来たと思います。
後は円の内側と外側で色分けする方法が分からずに色々調べていたら、実装するまでに時間が掛かってしまいました。
0 件のコメント:
コメントを投稿