プログラミング:モンテカルロ法で遊んでみた
モンテカルロ法で円周率の近似値を求める。図も描く。それを教材化するためのPythonのプログラムを書いていた。
1000個ではいい値が出ないので,個数を10000個にしたら,止まってしまった・・・? あれれ,何かミスったかな。
そうではなくて,単に時間がかかるだけだった。
マシンは,MacbookPro13' Late 2013。普段使いで特に不満はないマシンだ。大量の画像処理とか映像作品は扱っていない。
まてよ,12月に買った M1 Macmini ならどうだろう。
測ってみた。点を50000個打つ。このくらいなら結構いい近似値が出るのだ。
MacbookPro13' Late 2013 ・・・・ 1分16秒
Macmini (M1,2020) ・・・・ 31秒
さすがに M1 は速いねえ。
そうだ,2,3年前に授業でやったCinderellaではどうだろう。言語はCindyscript
Macmini(M1) でやってみた・・・・・ 1秒 れれれ?
間違いない。Cindyscriptはインタプリタなので遅いかと思ったがやたら速いじゃないか。そういえば,Java で作られていて,インタプリタだけど速い,と書いてあった。
まてよ,Pythonは Jupyter Notebook だから遅いのか。
そこで,IDLE でやってみた 20秒。 うん,少し速い。
直接やったら? ターミナルで,python3 monte.py として実行・・・ 18秒
それにしても,Cindyscriptが速いのは意外だった。
プログラムコードも見てわかる通り,ずっと簡単。
リストのインデックスが1スタートなのがなんとも。来るべき共通テストのプログラミングの言語仕様はインデックスが0スタートなのだ。たいしてかわらない?
いやいや,これは生徒にとってかなりのハードルだということが予想される。制御ブロックがインデントになっているのもPythonだし。Cindyscriptはインデントが推奨だが,なくても動く。いやはや・・・ (Cindyscriptの方が簡単なのに,と,またボヤキ)
----------------------------------------------------------
と,ここまで書いて下書き保存。その後,「まてよ,アルゴリズムを変えたらどうなる」と考えた。というのは,show() でも時間がかかっているからだ。その前のprint で円周率を表示した後さらに時間がかかっている。plotを変えたらどうか。
文科省の「高等学校情報科「情報Ⅰ」教員研修用教材」では,plot ではなく scatter を使っている。
plt.plot(x, y, 'ro',markersize = 0.5)
を
plt.scatter(x,y,c = 'r',s = 0.5)
に変えてやってみた。
すると,10000個の点でも40秒。50000個の点ではなかなか終わらないので停止ボタンで停止してみると,なんと plt.show でエラーになっていた。
ここでは scatter はやめましょうね ->文科省担当者殿
次に,逐次点を打つのではなく,座標を色データを配列に格納しておいて,まとめてプロットしたらどうか。とはいえ,プロットするのに,点によって赤・青が異なるのでやはり for で1点ずつ繰り返せざるを得ない。やってみると34秒(Jupyter)で改善はされなかった。
==== さらに追記 =====
三重大の奥村先生から助言をいただいた。以前にも拝見した,奥村研究室のページ。(忘れていた)
これにならってやってみたら,なんと1秒だった。速っ。
ただし,色は赤と青ではない。
ポイントは,for ループでその都度プロットするのではなく,50000回分の乱数のリスト(ベクトルデータ)を作ってしまうことだ。
そのことは考えた。Numpyで,random の引数に個数を渡すとその個数だけリストにして返すことは知っていたが(マニュアルに記載してある。忘れていたけど),四分円の内部と外部で色データも必要なのでできないと思っていた。まさか,色指定を真理値(True/False)で与えられるとはねえ。それに,count_nonzero は知らなかったし。
なお,scatter と plot では色指定の方法が異なるので,plot では色の(真理値の)ベクトル c で色分けすることはできないようだ。
これだけ多くの知識が必要となると,(しかも scatter の色指定が 真理値でもいいなどとは matplotlib のWebページを見てもわからない)今作っている,「これからプログラミングを始める高校教員向けの入門書」に入れるかどうかはなんともいえない。<参考>として書いておけばよいか。
=== 追記その2 ===
赤と青にする方法を思いついた。True/False を 'r' と 'b' に変えたベクトルを作ればいいのだ。(あたりまえ?)
def bool2col(bo):
if bo:
col = 'r'
else:
col = 'b'
return(col)
を定義しておいて(bool は予約語なのでbo)
c = []
for i in range(totalcount):
c.append(bool2col(x[i]*x[i] + y[i]*y[i] < 1))
とすればよい。速度はほとんど落ちない。