文科省教材:確率モデルのシミュレーション

高等学校情報科「情報Ⅰ」教員研修用教材(本編)の「第3章 コンピュータとプログラミング」の乱数を使ったシミュレーション(p.138〜)についてコメントする。

まず,Pythonでの乱数の作り方:

図表7 乱数発生のプログラム

この方法は互換性維持のため現在でも使えるが,2013年にリリースされた NumPy 1.17.0 からは Generator を使った新しいインターフェースが推奨されている:

import numpy as np

rng = np.random.default_rng()
ransuu = rng.random()
print(ransuu)

古い方法より1行余分に必要なので採用しなかったのかもしれないが,とりあえず指摘しておく。なお,ここでは不要であるが,実際のシミュレーションでは再現性が大切なので,乱数のタネ(seed)を例えば次のように設定する:

rng = np.random.default_rng(12345)

授業の課題では各自の学籍番号・出席番号をタネとして使わせてもよかろう。

次はさいころを100回振るプログラムである:

図表8 サイコロのプログラム

とりあえず Generator を使う方式に直して,PEP 8 に合わせて書き直す:

import numpy as np

rng = np.random.default_rng()
saikoro = rng.integers(1, 7, size=100)
print(saikoro)
deme = []
for i in range(6):
    deme.append(np.count_nonzero(saikoro == i + 1))
print("出現数:", deme)

np.count_nonzero() は,np.sum() や,単なる sum() でもよいが,この場合は np.count_nonzero() が最も速く,次に np.sum() が速い。

最後の4行は次のようにできる:

deme = [np.count_nonzero(saikoro == i + 1) for i in range(6)]
np.unique() を使えばもっとスマートにできる:

x, deme = np.unique(saikoro, return_counts=True)

これを使えば,次の度数分布図も簡単になる:

import matplotlib.pyplot as plt

plt.bar(x, deme)

さて,度数分布図である:

図表10 度数分布(棒グラフ)表現のプログラム

これはこれで問題ないが,align="center" はデフォルトなので省いて,次でもよい:

plt.bar(range(1, 7), deme)

E.V.ジュニアさんのnoteでは次の賢い方法が使われている:

plt.hist(saikoro, bins=np.arange(0.5, 7.5), rwidth=0.8)

seaborn を使えばもっと簡単である:

import seaborn as sns

sns.histplot(saikoro, discrete=True)

次はモンテカルロ法で円周率を求めるプログラムである:

図表12  円周率を求めるプログラム

Generator を用いる方式に直す:

import numpy as np

rng = np.random.default_rng()
totalcount = 1000
incount = 0
for i in range(totalcount):
    x = rng.random()
    y = rng.random()
    if x**2 + y**2 < 1:
        incount += 1
print("円周率:", incount * 4 / totalcount)

最後の7行はループをベクトル化すると次のようになる:

x = rng.random(totalcount)
y = rng.random(totalcount)
incount = np.count_nonzero(x**2 + y**2 < 1)
print("円周率:", incount * 4 / totalcount)

散布図を描くプログラム:

図表13 散布図を作成するプログラム

plt.scatter() を点ごとに呼び出すのはたいへん能率が悪い。ループ全体は次のようにできる:

x = rng.random(totalcount)
y = rng.random(totalcount)
c = x*x + y*y < 1
plt.scatter(x, y, c=c)
plt.axis('scaled')  # これをしないとアスペクト比が1にならない

色を選ぶ方法はいろいろ考えられそう。パレットをいじらずに済ませる方法:

col = ['red' if i else 'blue' for i in c]
plt.scatter(x, y, color=col)

別解

plt.scatter(x[c], y[c], color="red")
plt.scatter(x[np.logical_not(c)], y[np.logical_not(c)], color="blue")
# または plt.scatter(x[~c], y[~c], color="blue")

seaborn を使う例(デフォルトでうまい具合の色になる):

sns.scatterplot(x=x, y=y, hue=c)

ちなみに,文科省の研修用教材には他プログラミング言語版としてJavaScript版,VBA版,ドリトル版,swift版(なぜかSwiftだけ全角文字で小文字で始まっている)がある。2019年11月21日に出たVBA版では次のようになっていた:

図表13 散布図を作成するプログラム(VBA版)

実行結果は次の通り:

図表14 実行結果(VBA版)

これを見ればバグっていることは自明である。どこがおかしいかわかるだろうか。ちなみに私たちがツイッターで騒いだせいか,2019年12月11日にはリンク切れとなり,12月27日には修正版が出た。


Last modified: