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

高等学校情報科「情報Ⅰ」教員研修用教材（本編）の「第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だけ全角文字で小文字で始まっている）がある。2019年11月21日に出たVBA版では次のようになっていた：

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

実行結果は次の通り：

図表14 実行結果（VBA版）

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

Last modified: