プログラミング：モンテカルロ法で遊んでみた(2)

2021年2月6日 07:08

ここまでの話の流れは次の通り。

(1) 文科省から「高等学校情報科「情報Ⅰ」教員研修用教材」が出て，その中にモンテカルロ法を用いて円周率の近似値を求める例題があった。そこでは，scatter を使っているが，「正方形内にランダムに点を打つ」のであれば，それは統計の散布図ではないのだから scatter ではなく plot で打つべきだろう，という話を書いた。

散布図は「2つの量の関係を見るために描く図」だが，いまモンテカルロ法で使っているxとyの間には何の関係もないからだ。もちろん「散布状況を見る図」だから別にいいじゃないか，という意見もあるだろうし，それはそれでよい。ただ，ちょっと気持ち悪いだけだ。
　それで，scatter ではなく plot で点を打ったプログラムを載せた。

(2) 大学入試共通テストに「情報Ⅰ」を，との流れを受けて，これからプログラミングを始める高校教員向けの入門書を書き始めた。基本的なところ（変数，関数，リスト，制御構造）はだいたい書いて，情報の授業での利用例も書いて（基数変換など），モデル化とシミュレーションの例題を書いているところ。いままでのところは note にも書いてきた。

(3) モンテカルロ法による円周率の近似について書いているところで，(1) で示したプログラムだと，50000個も点を打つとかなり時間がかかることがわかり（scatter では処理できずエラーになる），面白そうなので MacbookPro (2013 Late) と Macmini(2020 M1) とで比較したら，さすがに M1は速かった，しかもCindyscriptなら1秒で終わる，という話を書いた。すると，奥村さんから助言があって，修正したら，あっというまに速くなった，という追記をした。

これで一件落着。・・・・　といかないのが面白いところ。

　改善点は何かというと，for による 50000 回の繰り返しごとに点を打つのでは処理が遅いので，50000個の点データを作ってまとめて打つ方がよい，ということであった。
　このとき，奥村さんのページにある方法だと，点はscatter で打つのだが，そのときの色指定を True/False でやっているので，赤/青とならない。そこで，各点に対する色のブール値を 'r' と 'b' に変えるルーチンを書いたら，速度はほとんど変わらずにうまくいった，という話である。ただし，この色指定の方法は plot では使えない。scatter と plot では色指定の方法が異なるからだ。

　やはり，気に入らない。scatter でなく plot でやれないのか。問題は，「繰り返しの中で plot すると時間がかかる」ということだ。先にデータを作ってしまって，まとめて plot すればいい，というところまではわかっている。しかし，色指定の方法が scatter と plot で異なるので・・・・・
　昼食後，買い物をしながら考えた。円の内部の点と外部の点データを別々に作って plot したらどうだ。

　やってみたらうまくいった。しかも遅くなっていないようだ。となると，こんどは速さを比べたくなる。内部時計を使って計測すればいい（という知識はCindyscriptで持っている）。そこでちょっと調べて，標準ライブラリの time を使って perf_counter() で計測することにした。

　まず，色指定をブール値で行う方法。奥村さんのページからもらってきたものだ。（同一ではない）
　コードとしてはこれは一番短く，エレガントかもしれない。ただし，色は赤/青にはならない。

import numpy as np
import matplotlib.pyplot as plt
import time
plt.figure(figsize=(4, 4))
plt.axis([0, 1, 0, 1])
rng = np.random.default_rng()
totalcount = 50000
start = time.perf_counter()      # 計測開始
x = rng.random(totalcount)
y = rng.random(totalcount)
incount = np.count_nonzero(x**2 + y**2 < 1)
c = x*x + y*y < 1
plt.scatter(x, y, c=c, s=0.1)
t = np.linspace(0, np.pi/2, 50)    # 四分円を描く
x = np.cos(t)
y = np.sin(t)
plt.plot(x, y, lw=0.5)
plt.title("Monte Carlo method")
print("円周率:", format(incount * 4.0 / totalcount, ".4f"))
plt.show()
print(time.perf_counter()-start)  # 計測終了

乱数の生成の関係で毎回同じにはならないが，だいたい 0.43 s 程度。

　次は，前回，追記で書いた，ベクトルc を色コードに書き換えたもの。当然その分遅くなる・・・・　だろう。
　変更点は次の通り（同じところも含む）

def bool2col(bo):       # bo は真理値 True/False。bool が予約語なので bo とした
   if bo:
       col = 'r'
   else:
       col = 'b'
   return(col)
totalcount = 50000
start = time.perf_counter()      # 計測開始
x = rng.random(totalcount)
y = rng.random(totalcount)
incount = np.count_nonzero(x**2 + y**2 < 1)
c = []
for i in range(totalcount):
   c.append(bool2col(x[i]*x[i] + y[i]*y[i] < 1))
plt.scatter(x, y, c=c, s=0.1)

たしかにちょっと遅くなっている。

ここまでは，前回の追記までにやったもので，それに計測のための2行を追加したものだ。

　次が，plot にこだわったもの。

totalcount = 50000
start = time.perf_counter()    # 計測開始
inx = []
iny = []
outx = []
outy = []
for i in range(totalcount):
   x = rng.random()
   y = rng.random()
   if x*x + y*y < 1:
       inx.append(x)
       iny.append(y)
   else:
       outx.append(x)
       outy.append(y)
plt.plot(inx, iny, 'or', markersize=0.3)
plt.plot(outx, outy, 'ob', markersize=0.3)

へへ，やったね(^^)

なお，Cindyscriptでは，繰り返すたびにプロットする方法で（まとめてプロットする関数はない）0.8 s くらいであった。

いやー，プログラミングって楽しいですね。