Pythonで〇×ゲームのAIを一から作成する その224 強化学習とモンテカルロ法

目次と前回の記事

Python のバージョンとこれまでに作成したモジュール

本記事のプログラムは Python のバージョン 3.13 で実行しています。また、numpy のバージョンは 2.3.5 です。

リンク	説明
marubatsu.py	Marubatsu、Marubatsu_GUI クラスの定義
ai.py	AI に関する関数
mbtest.py	テストに関する関数
util.py	ユーティリティ関数の定義
tree.py	ゲーム木に関する Node、Mbtree クラスなどの定義
gui.py	GUI に関する処理を行う基底クラスとなる GUI クラスの定義

AI の一覧とこれまでに作成したデータファイルについては、下記の記事を参照して下さい。

強化学習

これまでに作成した AI は、局面の優劣を表す評価値の計算 と、何手先までの局面を考慮するか を表す ゲーム木の探索 という 2 種類の方法を組み合わせて着手を計算 しました。

AI は 評価値の精度が高い程、探索するゲーム木 の深さが 深い程強く なります。例えば すべての局面 で 100 % 正しい評価値を計算 することができれば 1 手先の局面の評価値を調べるだけ で 強解決の AI を作成 することができます。また、ゲーム終了時まで の すべての局面を探索 することができれば、決着がついた局面だけ を 正しく評価できる評価関数 で 強解決の AI を作成 することができます。ただし、残念ながら下記の理由から ほとんどのゲーム では 全ての局面での 100 % 正しい評価値の計算 や、完全なゲーム木の探索 を行うことは 不可能または困難 です。

完全なゲーム木の探索 が 困難な理由 は多くのゲームでは ゲーム木の規模が大きすぎるから です。そのため、αβ 法 などのさまざまな ゲーム木を効率よく探索するアルゴリズム や、深さを制限したゲーム木の探索のアルゴリズム が利用されます。

これまでの記事で紹介 した局面の 評価値を計算 する 静的評価関数 は、人間がゲームの性質を考慮して考えたも のです。〇× ゲームのような単純なゲームでは 100 % 正しい評価値を計算する評価関数を作成することができましたが、オセロや将棋などの 多くのゲーム は 複雑すぎるため、100 % 正しい評価値を計算する静的評価関数 を 人間が考え出す ことは 非常に困難または不可能 です。

この問題を解決する方法の一つに 強化学習 があります。強化学習は プログラム が 周囲の環境を観察 することで どのように行動すべきを学習 するアルゴリズムで、〇× ゲーム の場合は ゲームの対戦を観察 することで、どの着手を行うべきかを学習 します。強化学習が発達する前 は、将棋や囲碁などの複雑なゲームの AI は 人間のプロよりもはるかに弱かった のですが、強化学習の研究の発達 によってに ゲームの AI が一気に強くなり、現在ではプロをはるかに超える実力を持つ ようになりました。

今回の記事から 〇× ゲームの AI を 強化学習で作成する様々な方法 について説明します。

参考までに下記に Wikipedia の強化学習のリンクを示します。

なお、強化学習 はプログラムが 経験から学習するアルゴリズムの総称 である 機械学習 の 分類の一つ で、機械学習には他にも人間が用意したお手本から学習する 教師あり学習、お手本なしで学習する 教師なし学習など があります。

参考までに下記に Wikipedia の機械学習のリンクを示します。

原始モンテカルロ法による 〇× ゲームの AI のアルゴリズム

〇× ゲーム のような 二人零和有限確定完全情報ゲーム に対する 強化学習 の中でおそらく 最も簡単な AI のアルゴリズム は 原始モンテカルロ法¹ではないかと思います。モンテカルロ法 は、乱数を利用した計算 を行うことで 近似値を求めるアルゴリズムの総称 で、〇× ゲームの原始モンテカルロ法 場合は以下のようなアルゴリズムで 最善手を計算 します。

1. 現在の局面 から 合法手を着手したすべての局面 に対して下記の計算を行う
   1. ゲームの決着がつくまで 乱数を利用して ランダムな着手を行い続け、その結果を記録 する。この作業を プレイアウト と呼ぶ
   2. あらかじめ決めておいた回数 または、あらかじめ決めておいた時間 になるまで プレイアウトを繰り返し、その 勝率を計算 する
2. 最も高い勝率が計算 された局面になる 合法手を最善手 とする

モンテカルロ法を知らなかった方は 上記のような乱数を使った方法 で 最善手を計算できるわけがないと思った人が多いのではないか と思いますが、モンテカルロ法 は ゲーム以外のざまざまな分野 でも 実際に利用されています。今回の記事では ゲーム以外の分野 で良く取り上げられる モンテカルロ法の代表例 を紹介します。

なお、モンテカルロ法 は 近似値を求める アルゴリズムであるため、一般的には 正確な値を求めることはできない 点に注意して下さい。

参考までに下記に Wikipedia のモンテカルロ法のリンクを示します。

モンテカルロ法はカジノで有名なモナコ公国のモンテカルロという地区の名前から付けられており、乱数を用いた計算のアルゴリズムがカジノで行われるギャンブルの偶然性に似ていると考えられたことが由来です。

モンテカルロ法による円周率の計算

モンテカルロ法 の具体例として良く取り上げられるのが 円周率 の 近似値を求める という例です。具体的には円周率である $\boldsymbol{π}$ = 3.1415926... は以下の 乱数を利用 した モンテカルロ法のアルゴリズム で計算することができます。

1. 下図のような 1 辺の長さが 2 で 原点を中心 とする 正方形 と、その 正方形に内接 する 原点を中心 とする 半径が 1 の円 を考えると、正方形の面積 は $2 × 2 = \boldsymbol{4}$、円の面積 は $π × 1^2 = \boldsymbol{π}$ となる
2. 乱数を利用 して図の 正方形の内部 に ランダムに $\boldsymbol{n}$ 個点を配置 する
3. ランダムに配置した点 が 円の内部に配置 される 確率 は 以下の式 で表される
   円の面積 ÷ 正方形の面積 = $\boldsymbol{\frac{π}{4}}$
4. 円の 内部に配置されたの点の数 を $\boldsymbol{i}$ とすると、大数の法則 から $\boldsymbol{n}$ が大きくなればなるほど $\boldsymbol{\frac{i}{n}}$ は $\boldsymbol{\frac{π}{4}}$ に近づく
5. 従って $\boldsymbol{π}$ は $\boldsymbol{\frac{4i}{n}}$ で近似 でき、$\boldsymbol{n}$ が大きい程 その 精度は高くなる

大数の法則 の 厳密な定義 は 次回の記事で詳しく説明する予定 なので、今回の記事では簡単に説明 します。また、今回の記事 では 大数の法則が正しいという前提 で説明を行います。

大数の法則 はサイコロの出目のような 特定の確率で発生する事象 を 数多く試行 することで、その 発生率 が 真の確率に近づいていく という現象のことです。例えば 出目の偏りのないサイコロ を 複数回振った際 に 出目が 1 となる割合 は、振った回数が少ない場合 は ばらつくことが多い ですが、多く振れば振る程 1 が出る確率 の 1/6 に近づいていきます。

近づくのは発生回数ではなく、割合を表す発生率 である点に注意して下さい。

このような現象 はおそらくほとんどの方が 実生活で経験 しているのではないかと思います。例えば勝負事は 一発勝負 では 弱い人が勝つこともある かもしれませんが、対戦回数が多くなる ほど 実力差に応じた勝敗結果 になります。野球などのスポーツが、長い時間をかけて多くの対戦を行うことで優勝者を決める ことが多いのは 実力をより正確に測ることができるから です。また、よほど大きな実力差がない限り、数回の対戦 では 優勝チームが最下位のチームに負け越す ことは それほど珍しいことではありません。

参考までに下記に Wikipedia の大数の法則のリンクを示します。

%timeit で処理時間を計算する場合も、なるべく多くの回数の処理時間の平均を計算することで大数の法則によって精度を高めています。

numpy による一様乱数の計算

上記のアルゴリズム で 円周率の近似値を計算できる ことを Python のプログラムで確認 することにします。また、numpy を利用することで プログラムを簡潔に記述でき 、高速に実行できる ので numpy を利用 することにします。

このアルゴリズムの 正方形 の x 座標と y 座標の範囲 は $\boldsymbol{-1 ≦ x ≦ 1}$、$\boldsymbol{-1 ≦ y ≦ 1}$ なので、正方形の内部 の ランダムな点 は x 座標と y 座標 の それぞれ に対して、その範囲の数値 を 均等な確率発生 させる 一様乱数（uniform random number） と呼ばれる 乱数を利用 することで計算できます。一様乱数 は numpy の 乱数を扱う random モジュール で定義された、下記の表の 仮引数 を持つ uniform という 関数で計算 することができます。なお、uniform は numpy を np という名前でインポート した場合は np.random.uniform と記述 します。

仮引数	意味
low	計算する一様乱数の 最小値
high	計算する一様乱数の 上界。ただし high は含まない
size	計算する一様乱数の 個数。整数を代入 すると その個数の要素 を持つ 1 次元の ndarray が計算 される。tuple を代入 すると その shape の ndarray が計算される

uniform は 複数の一様乱数をまとめて計算できる という点が 便利 です。また、ndarray が計算 されるので、後述するように 複数の要素をまとめて計算を行う ことができます。

下記は -1 以上 1 未満 の 一様乱数 を要素として持つ 下記の ndarray を計算 するプログラムです。実行結果からわかるように、np.random.uniform が 複数の乱数をまとめて作成 することが確認できます。

• 5 個の要素 を持つ 1 次元の ndarray
• (2, 3) の shape の 2 次元の ndarray

import numpy as np

print(np.random.uniform(low=-1.0, high=1.0, size=5))
print(np.random.uniform(low=-1.0, high=1.0, size=(2, 3)))

実行結果（実行結果はランダムなので下記と異なる場合があります。以後も同様です）

[-0.67055191 -0.72736821 -0.22868012  0.71370765  0.22886755]
[[-0.201856    0.8363786   0.33227609]
 [ 0.54750427 -0.51160777 -0.15898894]]

正方形の内部のランダムな点の計算と描画

下記は 正方形の内部 にある ランダムな 10 個の点を計算 するプログラムです。変数 X には 10 個の点の x 座標 を要素として持つ 1 次元の ndarray が計算 されます。変数 Y も同様 です。numpy を利用したプログラムでは、複数の要素を持つ ndarray が代入 された 変数の名前 を 大文字で表現することが多い ので、本記事でも 複数の x 座標が代入 された 変数の名前 を 大文字の X としました。以後はそのように変数の名前をつけることにします。

X = np.random.uniform(low=-1.0, high=1.0, size=10)
Y = np.random.uniform(low=-1.0, high=1.0, size=10)
print(X)
print(Y)

実行結果

[-0.30879553 -0.02810647 -0.24487785  0.71597128  0.69743244 -0.23973726
  0.52680899  0.57072444  0.17023537 -0.45887281]
[ 0.80941134 -0.91153367 -0.83717383 -0.32630709  0.40881786  0.08874042
  0.94306893 -0.50149512  0.41753724 -0.24005425]

uniform では仮引数 high に代入された値は計算されない ので上記では -1 以上 1 以下ではなく、-1 以上 1 未満の乱数が計算 されます。ただし、仮に uniform が -1 以上 1 以下の乱数を計算した場合 でも 1 が実際に計算される確率はほぼ 0 になる² ので、-1 以上 1 未満の乱数を計算しても 円周率の計算結果に影響はありません。

下記は上記で 計算した 10 個の点 を matplotlib で表示 するプログラムです。matplotlib について忘れた方は以前の記事を復習して下さい。実行結果から、先程計算した点のうち 8 つが円の内部に配置される ことが確認できました。なお、ランダムに配置 するので 下記のプログラムを実行するたび に 配置される点の位置 や 円の内部に配置される点の個数 は 変化する 点に注意して下さい。

• 4 行目：subplots で Figure と Axes を計算する
• 5 ~ 8 行目：Rectangle と Circle で正方形と円の図形を作成して Axes に登録する
• 9 行目：scatter で X、Y に代入された 10 個の点を描画する

1  import matplotlib.pyplot as plt
2  import matplotlib.patches as patches
3  
4  fig, ax = plt.subplots(figsize=(5, 5))
5  r = patches.Rectangle(xy=(-1, -1), width=2.0, height=2.0, fill=False, ec="k")
6  c = patches.Circle(xy=(0, 0), radius=1, fill=False,ec="r")
7  ax.add_patch(r)
8  ax.add_patch(c)
9  ax.scatter(X, Y)

行番号のないプログラム

import matplotlib.pyplot as plt
import matplotlib.patches as patches

fig, ax = plt.subplots(figsize=(5, 5))
r = patches.Rectangle(xy=(-1, -1), width=2.0, height=2.0, fill=False, ec="k")
c = patches.Circle(xy=(0, 0), radius=1, fill=False,ec="r")
ax.add_patch(r)
ax.add_patch(c)
ax.scatter(X, Y)

実行結果

円の内部にあるかどうかの判定方法

下記は (x, y) の点が 原点を中心とする半径 1 の円の内部 にある³ことを表す 不等式 です。

$\boldsymbol{x^2 + y^2 ≦ 1^2 = 1}$

下記は 先程計算した 10 個の点 が 円の内部にあるかどうか を 一つずつ上記の式で計算 するプログラムです。実行結果に numpy の True を表す データである np.True_ が 9 つ表示される ことから、円の内部に 9 つの点があることが確認 できました。

• 1 行目：各点が円の内部（inner）にあるかどうかを計算する変数 I を空の list で初期化する
• 2 行目：組み込み関数 zip を利用して X と Y の各要素を順番に x、y に代入しながら繰り返し処理を行う
• 3 行目：上記の不等式の計算結果を I の要素に追加する

1  I = []
2  for x, y in zip(X, Y):
3      I.append(x * x + y * y <= 1)
4  print(I)

行番号のないプログラム

I = []
for x, y in zip(X, Y):
    I.append(x * x + y * y <= 1)
print(I)

実行結果（はみ出るので途中で改行しました）

[np.True_, np.True_, np.True_, np.True_, np.True_, 
 np.True_, np.False_, np.True_, np.True_, np.True_]

I = [x * x + y * y <= 1 for x, y in zip(X, Y)] のように list 内包表記を利用して 1 行で記述することもできます。

上記のプログラムでは for 文 で 点を一つずつ確認 しましたが、以前の記事で説明した numpy の ブロードキャスト の機能を利用することで、下記のプログラムのように 1 行で上記と同様の計算 を行うことができます。下記のプログラムでは ndarray の各要素 に対して それぞれ個別 に $\boldsymbol{x^2 + y^2 ≦ 1}$ の計算が行われ てその結果が X や Y と同じ形状 の 1 次元の ndarray として計算 されます。実行結果から 上記と同じ意味を持つデータが計算される ことが確認できました。こちらの方が 簡潔に記述できる だけでなく、処理速度も高速 なので本記事ではこちらの方法を採用することにします。

I = X * X + Y * Y <= 1
print(I)

実行結果

[ True  True  True  True  True  True False  True  True  True]

list の要素 に numpy で計算したデータを代入 したものを print で表示 すると、先程の実行結果の np.True_ のように list の要素に numpy のデータ型が代入されていることが分かるような表示 が行われます。

一方、ndarray を print で表示 すると上記の実行結果のように 単に True のような値が表示 されますが、それらの値は numpy のデータ型 のデータです。

円の内外を区別した描画

円の内外を区別 できるように 円の内部 の点を 赤色 で、外部 の点を 黒色で描画 する方法を紹介します。その方法の一つに、下記のプログラムのように先ほど計算した 円の内部の点であるかどうか を表す I に代入された ndarray を利用 して、円の内部の点 を表す x 座標と y 座標 を要素として持つ list をそれぞれ計算 するという方法が考えられます。円の外部 の点についても 同様 です。実行結果から 正しく描画されることが確認 できます。

• 1 ~ 4 行目：円の内部（inner）の点の x, y 座標の一覧 を表す IX と IY、円の外部（outer）の点の x, y 座標の一覧 を表す OX と OY をそれぞれ 空の list で初期化 する
• 5 行目：組み込み関数 zip を利用して X、Y、I の 各要素 を順番に x、y、inner に代入 しながら 繰り返し処理 を行う
• 6 ~ 11 行目：内部の点 の場合は IX と IY、外側の点 の場合は OX と OY の要素に 点の座標のデータを追加 する
• 18、19 行目：内部の点 を ax.scatter の実引数に c="r" を記述 して 赤色で描画 し、外部の点 を ax.scatter の実引数に c="k" を記述 して 描画で描画 する。色の指定方法 について忘れた方は以前の記事を復習すること

 1  IX = []
 2  IY = []
 3  OX = []
 4  OY = []
 5  for x, y, inner in zip(X, Y, I):
 6      if inner:
 7          IX.append(x)
 8          IY.append(y)
 9      else:
10          OX.append(x)
11          OY.append(y)
12  
13  fig, ax = plt.subplots(figsize=(5, 5))
14  r = patches.Rectangle(xy=(-1, -1), width=2.0, height=2.0, fill=False, ec="k")
15  c = patches.Circle(xy=(0, 0), radius=1, fill=False,ec="r")
16  ax.add_patch(r)
17  ax.add_patch(c)
18  ax.scatter(IX, IY, c="r")
19  ax.scatter(OX, OY, c="k")

行番号のないプログラム

IX = []
IY = []
OX = []
OY = []
for x, y, inner in zip(X, Y, I):
    if inner:
        IX.append(x)
        IY.append(y)
    else:
        OX.append(x)
        OY.append(y)

fig, ax = plt.subplots(figsize=(5, 5))
r = patches.Rectangle(xy=(-1, -1), width=2.0, height=2.0, fill=False, ec="k")
c = patches.Circle(xy=(0, 0), radius=1, fill=False,ec="r")
ax.add_patch(r)
ax.add_patch(c)
ax.scatter(IX, IY, c="r")
ax.scatter(OX, OY, c="k")

実行結果

masked array を利用した描画方法

上記のプログラムは 円の内部 と 外部の点 の 座標を表す list を 繰り返し処理を利用して計算 する必要 がありましたが、ndarray の 特定の要素を隠す（mask） ことができる masked array を利用することで もっと簡単に上記と同じ処理 を行うことができます。

masked array は numpy の ma という モジュールで定義 された MaskedArray という名前の クラスのインスタンス で、ndarray とほぼ同じ性質 を持ちますが、マスク（mask）した要素 が 隠されて計算の際に利用されない という性質があります。なお、このマスクという用語は以前の記事で説明した ビットマスク のマスクと 同じ用語 です。

masked array は 様々な方法で作成 することができますが、本記事では MaskedArray の array という メソッドを利用して作成 する方法を紹介します。numpy を np という名前でインポート した場合は np.ma.array と記述 するので、以後はそのように表記します。

np.ma.array は 仮引数 として list または ndarray を代入 する data と、data の中 で マスクして隠す要素 を表す mask を持ちます。mask には一般的に data と同じ形状 の list や ndarray を代入 し、True が代入された要素に対応 する data の要素が隠されて利用されなく なります。

下記は [1, 2, 3] という list のうちの 1 番の要素 である 2 をマスクして隠した masked array を作成 するプログラムです。仮引数 mask には [1, 2, 3] と同じ 3 つの要素 を持ち、隠したい要素に対応 する 1 番の要素が True で それ以外の要素が False である [False, True, False] という list を代入 します。

masked array を print で表示 すると、実行結果のように マスクして隠された要素 は -- が表示 されます。また、np.sum で 要素の合計を計算 すると マスクされた要素が無視される ので実行結果のように 1 + 3 = 4 が計算 されます。

ma = np.ma.array(data=[1, 2, 3], mask=[False, True, False])
print(ma)
print(np.sum(ma))

実行結果

[1 -- 3]
4

先程計算した I には 円の内部 にある点に 対応する要素が True となる ndarray が代入 されているので、下記のプログラムで 円の内部の点 を マスクして隠した、円の外部（outer）の点の x 座標の一覧を計算 することができます。実行結果から、円の内部 の点に 対応する 9 つの要素が隠されて -- のように表示される ことが確認できます。

X_outer = np.ma.array(X, mask=I)
print(X)
print(X_outer)

実行結果

[-0.30879553 -0.02810647 -0.24487785  0.71597128  0.69743244 -0.23973726
  0.52680899  0.57072444  0.17023537 -0.45887281]
[-- -- -- -- -- -- 0.5268089917275507 -- -- --]

円の内部 の点の x 座標の一覧を計算 する場合は、I の要素 の True と False を反転した ndarray を mask に代入 する必要があります。そのような ndarray は下記のプログラムの 2 行目のように 以前の記事 で説明した True と False を反転 させる ~ という NOT 演算子 で計算することができます。従って、円の内部（inner）の点の x 座標の一覧 は下記の 4 行目のプログラムで計算することができます。実行結果から、円の内部 の点に 対応する 1 つの要素が隠されて -- のように表示される ことが確認できます。

1  print(I)
2  print(~I)
3  print()
4  X_inner = np.ma.array(X, mask=~I)
5  print(X)
6  print(X_inner)

行番号のないプログラム

print(I)
print(~I)
print()
X_inner = np.ma.array(X, mask=~I)
print(X)
print(X_inner)

実行結果

[ True  True  True  True  True  True False  True  True  True]
[False False False False False False  True False False False]

[-0.30879553 -0.02810647 -0.24487785  0.71597128  0.69743244 -0.23973726
  0.52680899  0.57072444  0.17023537 -0.45887281]
[-0.3087955278719061 -0.02810647135091693 -0.2448778462923149
 0.7159712825357611 0.6974324374039698 -0.23973725640451948 --
 0.5707244381212662 0.17023537223780005 -0.45887281453520856]

下記は上記で計算された masked array を利用して 円の内部 の点を 赤色 で、円の外部 の点を 黒色で描画 するプログラムです。matplotlib の x 座標と y 座標を指定して描画を行う関数 は、片方の座標がマスクされて隠されていれば、もう片方の座標がマスクされていなくても その点は描画されない ようになっているので、y 座標の一覧 は マスクしていない Y のままで構いません。実行結果から 先程と同じ画像が描画される ことが確認できました。

fig, ax = plt.subplots(figsize=(5, 5))
r = patches.Rectangle(xy=(-1, -1), width=2.0, height=2.0, fill=False, ec="k")
c = patches.Circle(xy=(0, 0), radius=1, fill=False,ec="r")
ax.add_patch(r)
ax.add_patch(c)
ax.scatter(X_inner, Y, c="r")
ax.scatter(X_outer, Y, c="k")

実行結果

masked array、masked array の作成方法、np.ma.array の詳細は下記のリンク先を参照して下さい。

円周率の計算

先程説明したように、モンテカルロ法 による 円周率の近似値の計算式 は、配置した 点の総数を $\boldsymbol{n}$、円の内部 の 点の数 を $\boldsymbol{i}$ とした場合に $\boldsymbol{\frac{4i}{n}}$ という 式で計算 できます。円の内部の点の数 は I の True の要素の数 で True は 0 ではない値 なので、下記のプログラムのように 0 でない要素の数を数える np.count_nonzero で計算することができます。

i = np.count_nonzero(I)
print(i)

実行結果

9

従って、下記のプログラムで 円周率の近似値を計算 することができます。変数の名前 は モンテカルロ法（monte carlo) で計算した値なので mc_pi としました。また、numpy モジュールでは $\boldsymbol{π}$ の値 が np.pi に代入 されている⁴ので 誤差も計算 しました。

n = 10
mc_pi = 4 * i / n
print(f"mc_pi = {mc_pi}")
print(f"np.pi = {np.pi}")
print(f"誤差  = {mc_pi - np.pi}") 

実行結果

mc_pi = 3.6
np.pi = 3.141592653589793
誤差  = 0.458407346410207

実行結果から、10 個の点をランダムに配置 して計算した 今回の場合 の 円周率の近似値は 3.6、円周率との誤差 は 約 0.45 であることが確認できました。乱数の値によって異なる結果になる ので興味がある方は 上記のプログラムを何度か実行 してみて下さい。

大数の法則の確認

次に、配置する点の数を増やす ことで 大数の法則 によって計算される 円周率の近似値の精度が高くなる ことを確認することにします。

配置する点の個数と精度の最大値の関係

先程説明したように、円周率の近似値 は配置した 点の総数を $\boldsymbol{n}$、円の内部 の 点の数 を $\boldsymbol{i}$ とした場合に $\boldsymbol{\frac{4i}{n}}$ という 式で計算 できます。 $\boldsymbol{i}$ は整数 なので、計算される 円周率の近似値 は $\boldsymbol{i}$ が 1 増える と $\boldsymbol{\frac{4}{n}}$ 増える ことになります。下記は $\boldsymbol{n}$ を 10、100、1,000、10,000、100,000、1,000,000 とした場合の $\boldsymbol{\frac{4}{n}}$ の値を表す表です。また、$\boldsymbol{\frac{4i}{n}}$ という式で 計算できる $π$ = 3.1415926... 以下と以上 である 最も精度の高い値 を示しました。

$\boldsymbol{n}$	$\boldsymbol{\frac{4}{n}}$	$π$ 以下の最も精度の高い値	$π$ 以上の最も精度の高い値
10	0.4	2.8	3.2
100	0.04	3.12	3.16
1,000	0.004	3.140	3.144
10,000	0.0004	3.1412	3.1416
100,000	0.00004	3.14156	3.14160
1,000,000	0.000004	3.141592	3.141596

上記の表から、点の個数を 10 倍 にすることで 最も精度が高い近似値が計算 された場合に 正しい値が計算される小数点以下の桁数が 1 増える ので 誤差が約 1/10 になる ことが確認できました。また、点の個数を 100 万 にした場合でも、最大 で 3.14159 という 小数点以下 5 桁までの精度 の近似値 しか計算できない ことも確認できました。なお、実際に計算される円周率の近似値の精度 は それより低くなる ので そのことを確認する ことにします。

calc_mc_pi の定義

点の数を増やした場合 に計算されれる 円周率の近似値 と 配置された点を図示 できるように、任意の数の点 でモンテカルロ法で 円周率の近似値を計算 し、配置された点を図示 する calc_mc_pi を下記のプログラムように 定義 することにします。なお、仮引数 num は 配置する点の数 を表します。

• 1 行目：仮引数 num を持つ calc_mc_pi を定義する
• 3、4 行目：正方形内の num 個のランダムな点の座標を計算する
• 5 ~ 10 行目：円の内部の点の数を計算し、円周率の近似値と誤差を表示する
• 12 ~ 20 行目：先程と同じ方法でランダムに配置された点を表す画像を描画する

 1  def calc_mc_pi(num):
 2      print(f"点の個数 = {num}")
 3      X = np.random.uniform(low=-1.0, high=1.0, size=num)
 4      Y = np.random.uniform(low=-1.0, high=1.0, size=num)    
 5      I = X * X + Y * Y <= 1
 6      i = np.count_nonzero(I)
 7      mc_pi = 4 * i / num
 8      print(f"mc_pi = {mc_pi}")
 9      print(f"np.pi = {np.pi}")
10      print(f"誤差  = {mc_pi - np.pi}")     
11  
12      X_inner = np.ma.array(data=X, mask=~I)
13      X_outer = np.ma.array(data=X, mask=I)
14      fig, ax = plt.subplots(figsize=(5, 5))
15      r = patches.Rectangle(xy=(-1, -1), width=2.0, height=2.0, fill=False, ec="k")
16      c = patches.Circle(xy=(0, 0), radius=1, fill=False,ec="r")
17      ax.add_patch(r)
18      ax.add_patch(c)
19      ax.scatter(X_inner, Y, c="r")
20      ax.scatter(X_outer, Y, c="k")

行番号のないプログラム

def calc_mc_pi(num):
    print(f"点の個数 = {num}")
    X = np.random.uniform(low=-1.0, high=1.0, size=num)
    Y = np.random.uniform(low=-1.0, high=1.0, size=num)    
    I = X * X + Y * Y <= 1
    i = np.count_nonzero(I)
    mc_pi = 4 * i / num
    print(f"mc_pi = {mc_pi}")
    print(f"np.pi = {np.pi}")
    print(f"誤差  = {mc_pi - np.pi}")     

    X_inner = np.ma.array(data=X, mask=~I)
    X_outer = np.ma.array(data=X, mask=I)
    fig, ax = plt.subplots(figsize=(5, 5))
    r = patches.Rectangle(xy=(-1, -1), width=2.0, height=2.0, fill=False, ec="k")
    c = patches.Circle(xy=(0, 0), radius=1, fill=False,ec="r")
    ax.add_patch(r)
    ax.add_patch(c)
    ax.scatter(X_inner, Y, c="r")
    ax.scatter(X_outer, Y, c="k")

様々な個数の点を配置した場合の計算

下記は 10、100、1,000、10,000、100,000、1,000,000 個 の点を ランダムに配置 した場合の 円周率の近似値の計算 と 配置された点を描画 するプログラムです。下記の実行結果には近似値だけを表記し、描画された画像は後でまとめます。

for num in [10, 100, 1000, 10000, 100000, 1000000]:
    calc_mc_pi(num)

実行結果

点の個数 = 10
mc_pi = 3.2
np.pi = 3.141592653589793
誤差  = 0.05840734641020706
点の個数 = 100
mc_pi = 3.2
np.pi = 3.141592653589793
誤差  = 0.05840734641020706
点の個数 = 1000
mc_pi = 3.192
np.pi = 3.141592653589793
誤差  = 0.050407346410207055
点の個数 = 10000
mc_pi = 3.1404
np.pi = 3.141592653589793
誤差  = -0.0011926535897930357
点の個数 = 100000
mc_pi = 3.1322
np.pi = 3.141592653589793
誤差  = -0.00939265358979302
点の個数 = 1000000
mc_pi = 3.140156
np.pi = 3.141592653589793
誤差  = -0.0014366535897929467

下記は上記の実行結果の 点の数、円周率の近似値、誤差 をまとめた表です。先程説明したように、点の数が多くなる ほど より小さな小数点以下の値が計算される ようになることが確認できました。また、点の数が多くなる ほど 誤差が小さくなる傾向がある ことが確認できます。なお、点の数が少なくても 運よく誤差が小さい近似値が計算される場合がある ので、点の数が多くなる と 必ずしも誤差が少なくなるとは限らない 点に注意して下さい。実際に点の数を 100,000 とした場合の誤差は 10,000 の場合よりも多くなっています。

点の数	円周率の近似値	誤差
10	3.2	0.05840734641020706
100	3.2	0.05840734641020706
1,000	3.192	0.050407346410207055
10,000	3.1404	-0.0011926535897930357
100,000	3.1322	-0.00939265358979302
1,000,000	3.140156	-0.0014366535897929467

また、上記の表から 小数点以下第 3 桁 で 近似値を四捨五入 した場合に 数学で良く利用される 3.14 を計算するためには、点の数が 10000 個必要 になることが確認できました。このように、モンテカルロ法 による 円周率の近似値の計算 は、精度が非常に高い近似値を求めるため には 非常に大きな数の点が必要 になります。一方で、実用的な 3.14 という円周率の近似値を計算 するためにはコンピューターでは 一瞬で計算を行える 10000 個の点で十分に計算できる ことも確認できました。

点の個数と誤差 の より正確な関係 については 次回の記事で説明 する予定です。

2 次元のヒストグラムによる配置された点のばらつきの描画

下記は上記のプログラムで描画された ランダムに配置された点の画像 をまとめたものです。

点の数	画像	点の数	画像
10		100
1,000		10,000
100,000		1,000,000

精度の高い円周率の近似値を計算 するためには、ランダムに配置された点が 正方形の中 にまんべんなく 均等に配置される必要 があります。上記の画像 から以下の事がわかるため、点の数が多くなる と計算された 近似値の精度が高くなる ことがわかります。

• 点の数が小さい 場合は 点の配置が特定の場所に片寄る傾向が高い
• 点の数が多い 程点がまんべんなく 均等に配置される ようになる

先程説明したように、大数の法則 は 特定の確率で発生する事象 が、数多く試行 することでその 発生率 が 真の確率に近づいていく という現象です。上記のようになる ことは、大数の法則 によって 以下のように説明 できます。

• 正方形 を 同じ面積になる ように小さく $\boldsymbol{m}$ 個に区切る と、それぞれの範囲 に 点が配置される確率 は $\boldsymbol{\frac{1}{m}}$ になる
• 全体に配置する点の個数 を $\boldsymbol{n}$ とすると、大数の法則 により それぞれの範囲に配置される点を $\boldsymbol{n}$ で割って計算 する 割合 は、$\boldsymbol{n}$ が大きくなればなるほど $\boldsymbol{\frac{1}{m}}$ に近くなる
• このことから、$\boldsymbol{n}$ が小さい場合 はそれぞれの範囲に 配置される点の割合がばらつく ので、点が均等に配置されない可能性が高く なる
• $\boldsymbol{n}$ が大きい程 それぞれの範囲に 配置される点の割合 が$\boldsymbol{\frac{1}{m}}$ に近くなる ので、点が均等に配置される可能性が高く なる

このことを 実際に図で確認する ことにします。先ほどの図 は、点の数が多くなる と 同じ場所に同じ色で多くの点が描画される ので どのくらいの数の点が配置されているかがわかりづらくなるという欠点 があります。そのことは 1,000,000 個の点を配置した場合 の図が 真っ赤に塗りつぶされる ことから明らかでしょう。そこで、どこにどれくらいの数の点が配置されたかがわかる ような 図を描画 することにします。

データの 特定の範囲の値ごとの個数 は ヒストグラム（度数分布表）で表現することができます。以前の記事で紹介したヒストグラムは 1 種類のデータの集合 の 特定の範囲の個数をグラフで表現 するというものでしたが、上記の点の集合は x 座標と y 座標 という 2 種類のデータの組み合わせの集合 です。そのような場合は 2 次元の座標 を 正方形の範囲で区切り、その範囲内のデータの個数を表現 する 2 次元のヒストグラム を利用できます。

2 次元のヒストグラムの代表例 としては、天気予報でよくみかける 地図を正方形の範囲で区切り、各範囲内の雨量 を 正方形を塗りつぶした色で表す というものがあります。下図は日本気象協会のホームページに掲載されていた降雨量の図の一部です。著作権を考慮して図の一部のみを抜粋しました。図からわかるように、地図が正方形で区切られ、その正方形に降水量を表す色が塗りつぶされています。

matplotlib には上記のような 2 次元のヒストグラム を Axes の hist2d というメソッドで描画することができるのでその方法を紹介します。

網目状に区切られたもののことをメッシュと呼び、上記のような正方形で網目状に区切られた地図のことをメッシュ地図と呼びます。

数学やプログラミングでは、計算する領域を正方形などで有限個に分割した領域の事をメッシュと呼びます。

下図のように正方形で区切られた部分に棒グラフを表示するという、3 次元の立体的な画像で 2 次元のヒストグラムを表すという方法もあります。なお、下図は東海テレビのホームページにあった図の一部を抜粋しました。

上記のような図はデータの量が直観的に分かりやすくなるという利点がありますが、棒グラフの後ろにあるデータが隠れてしまうという欠点があるので本記事では採用しません。matplotlib では上記のような 3 次元のグラフを描画する機能があるので興味がある方は調べてみると良いでしょう。

正方形ではなく、正六角形で区切る場合もあります。matplotlib では Axis の hexbin メソッドでそのようなヒストグラムを描画することができます。詳細は下記のリンク先を参照して下さい。

hist2d メソッドの利用方法

Axes の hist2d メソッドには下記のような 仮引数 があります。なお、以下の説明では 2 次元のヒストグラムで個数を表現する データ を、x 座標と y 座標 で表される 2 次元の座標で表されるデータとみなす ことにします。

仮引数	意味
x	2 次元の座標のうちの x 座標の一覧
y	2 次元の座標のうちの y 座標の一覧
bins	2 次元の x 座標と y 座標を区切る個数5
range	図示する x 座標と y 座標の範囲。下記の 2 次元の list で指定する [[x 座標の最小値, x 座標の最大値], [y 座標の最小値, y 座標の最大値]] 省略すると x と y で指定したデータから自動的に設定される。
weights	点の数を数える際に乗算する重み（weights) の一覧。点の個数だけ指定する必要がある。省略するとすべての重みが 1 で計算される
cmap	個数を表す色を表すカラーマップ。省略するとあらかじめ決められたカラーマップが使用される。詳細はこの後で説明する
vmin、vmax	カラーマップが表現する値の最小値と最大値 指定しない場合はデータから自動的に計算した値が設定される

Axes の hist2d メソッドの詳細は下記のリンク先を参照して下さい。

下記は先ほどの 10 個の点を配置した場合 の 2 次元のヒストグラムを描画 するプログラムです。hist2d の 最初の 2 つの仮引数 は x と y なので、ax.hist2d(x=X, y=Y, bins=20) の代わりに 下記のプログラムのように記述することが多い のではないかと思います。2 行目の最後に ; を記述 したのは hist2d の大量の返り値を VS Code のセルに表示しないようにするため です。なお、hist2d の返り値の一部については後述します。

bins=20 を記述することで x 座標と y 座標 を それぞれ 20 個に区切った 2 次元のヒストグラムが描画されます。また、実行結果のように 点が配置された範囲は黄色 で、配置されていない範囲は紫色 で塗りつぶされます。bins=20 とした理由 は -1 ~ 1 の範囲を 20 等分 することで、区切られた正方形の一辺の長さ を 0.1 という きりの良い長さにするため です。興味がある方は別の値を設定してみて下さい。

fig, ax = plt.subplots(figsize=(5, 5))
ax.hist2d(X, Y, bins=20);

実行結果

仮引数 range の利用方法

仮引数 range の値を指定しない場合 は 図示する x 座標と y 座標の範囲 は X と Y のデータから自動的に決まり ます。下記のプログラムの実行結果からわかるように X の最小値と最大値 は 約 -0.46 と 0.71 なので、図の x 軸に表示される数値の範囲 から x 座標の範囲 も そのような範囲 になっていることが確認できます。Y の範囲の確認は省略しますが y 座標に関しても同様 です。そのため、正方形の範囲の一部が描画されない という問題が発生します。

print(min(X), max(X))

実行結果

-0.45887281453520856 0.7159712825357611

正方形の範囲を全て描画する ためには、下記のプログラムのように 仮引数 range に x 座標と y 座標の範囲 としてそれぞれ -1 ~ 1 の範囲を表すデータを代入 します。実行結果の x 軸と y 軸に表示される数値の範囲 から、2 次元のヒストグラムに描画される範囲が 正方形の範囲と同じになる ことが確認できました。

fig, ax = plt.subplots(figsize=(5, 5))
ax.hist2d(X, Y, bins=20, range=[[-1, 1], [-1, 1]]);

実行結果

カラーバーの描画

下図は 点を 1000 個配置した場合 の 2 次元のヒストグラムを描画するプログラムです。実行結果のように正方形には その範囲の点の個数 を表す 様々な色で塗りつぶされます が、それぞれの色が表す個数 が 何であるかがわからない という問題があります。

X = np.random.uniform(low=-1.0, high=1.0, size=1000)
Y = np.random.uniform(low=-1.0, high=1.0, size=1000) 
fig, ax = plt.subplots(figsize=(5, 5))
ax.hist2d(X, Y, bins=20, range=[[-1, 1], [-1, 1]]);

実行結果

この問題は 色と個数の対応 を表す カラーバーを描画 することで解決することができます。カラーバー は Figure の colorbar メソッドで描画 することができます。また、hist2d で描画した場合の カラーバー は、下記のプログラムのように hist2d の返り値の 3 番の要素を実引数に記述 することで描画されます。実行結果から 2 次元のヒストグラムの横 に 色と個数の対応を表すカラーバーが描画される ことが確認できました。

fig, ax = plt.subplots(figsize=(5, 5))
h = ax.hist2d(X, Y, bins=20, range=[[-1, 1], [-1, 1]])
fig.colorbar(h[3])

実行結果

Figure の colorbar メソッドの詳細は下記のリンク先を参照して下さい。

Axes の軸のスケールを同じにする方法

上記のプログラムで カラーバーを描画 すると先程の実行結果のように 2 次元のヒストグラム がその分だけ 左右に縮んで描画されてしまう という 問題が発生 します。Axes には x 軸と y 軸 の スケール（倍率） の 縦横比（aspect）を設定 する set_aspect というメソッドがあり、下記のプログラムのように 実引数 に 軸のスケールを同じ（equal) にする "equal" を記述 して呼び出すことで、実行結果のようにこの問題を解決することができます。

fig, ax = plt.subplots(figsize=(5, 5))
h = ax.hist2d(X, Y, bins=20, range=[[-1, 1], [-1, 1]])
fig.colorbar(h[3])
ax.set_aspect("equal")

実行結果

2 次元のヒストグラムとカラーバーの高さをあわせたい場合は、Figure の横幅を大きくすると良いでしょう。

Axes の set_aspect メソッドの詳細は下記のリンク先を参照して下さい。

仮引数 weights の利用方法

先程の実行結果の図の カラーバーの横に表示される数値 からわかるように、作成される 2 次元のヒストグラム は色が 配置された点の個数を表します が、先程説明したように 大数の法則 は個数ではなく 割合を対象 とする法則です。割合 はそれぞれの範囲に 配置した点の数 を、全体に配置した点の総数 で 除算する ことで計算することができます。また、そのような計算は 配置した点の数を数える際 に、1 個として数える代わりに 1 / 全体に配置する点の総数 個として数えるという 方法でも計算 することができます．

A、B、C、A、A という 5 つの文字の中の A の割合の計算は、下記の 2 種類の方法で計算できます。

• A の数を数えてから 5 で割り算する。3 / 5 = 0.6 となる
• A の数を数える際に、1 ではなく 1 / 5 個として数える。1 / 5 * 3 = 0.6 となる

そのような計算は hist2d の 仮引数 weights を利用することで行うことができます。仮引数 weights には 仮引数 x、y で指定した それぞれの座標の点を数える際 の 個数 を list や ndarray で指定 します。例えば、3 つの点の座標 を x = [1, 2, 3]、y = [4, 5, 6] で指定した場合に、1 つ目の点を 2 個、2 つ目の点を 3 個、3 つ目の点を 0.5 個 として 数えたい場合 は weights = [2, 3, 0.5] を指定 します。

例えば 1,000 個の点を配置 する場合に、それぞれの範囲に 配置された点の割合を計算 するためには、1,000 個のそれぞれの点を数える際 に 1 / 1000 個として数える 必要があるので、1,000 個 の 1 / 1000 を要素 とする list として weights=[1 / 1000] * 1000 のように記述します。別の記述方法として、以前の記事で説明した 同じ要素を持つ ndarray を作成 する np.full を利用して weights=np.full((1000, ), 1 / 1000) のように記述しても良いでしょう。本記事では後者を採用することにします。

下記は先ほどの 点を 1000 個配置 した場合の 個数の割合 の 2 次元のヒストグラムを描画するプログラムです。実行結果の左の図から カラーバーの横 に 割合を表す 0 以上 1 未満の数値 が描画されるようになったことが確認できます。先程の実行結果の右図 と見比べてみればわかると思いますが、個数とその割合 は 同じ色に割り当てられる ので どちらも同じ範囲が同じ色で塗りつぶされます。

fig, ax = plt.subplots(figsize=(5, 5))
weights = np.full((1000, ), 1 / 1000)
h = ax.hist2d(X, Y, bins=20, range=[[-1, 1], [-1, 1]], weights=weights)
fig.colorbar(h[3])
ax.set_aspect("equal")

実行結果（右は先ほどの図です）

　

hist2d の仮引数 density に True を代入することで同様のグラフを描画することができますが、その場合に描画されるのは確率密度（probability density）とよばれるグラフなので上記の結果とは異なる値が計算されます。具体的には、確率密度の場合は、「それぞれの範囲の面積 × その範囲に配置された個数」の合計が 1 になるような値が計算される点が異なります。

下記は density=True を記述した場合のプログラムです。それぞれの範囲の面積は 0.1 × 0.1 = 0.01 なので、それぞれの範囲の個数は実行結果のカラーバーの右の数値から確認できるように上記の場合の 100 倍の値が計算されます。

fig, ax = plt.subplots(figsize=(5, 5))
h = ax.hist2d(X, Y, bins=20, range=[[-1, 1], [-1, 1]], density=True)
fig.colorbar(h[3])
ax.set_aspect("equal")

実行結果

参考までに Wikipedia の確率密度の項目のリンクを下記に示します。

仮引数 vmin、vmax の利用方法

下記は、点を 10,000 個 配置した場合と 1,000,000 個 配置した場合の 個数の割合 の 2 次元のヒストグラムを描画するプログラムです。

for size in [10000, 1000000]:
    X = np.random.uniform(low=-1.0, high=1.0, size=size)
    Y = np.random.uniform(low=-1.0, high=1.0, size=size) 
    fig, ax = plt.subplots(figsize=(5, 5))
    weights = np.full((size, ), 1 / size)
    h = ax.hist2d(X, Y, bins=20, range=[[-1, 1], [-1, 1]], weights=weights)
    fig.colorbar(h[3])
    ax.set_aspect("equal")

実行結果

　

実行結果からわかるように、どちらも 同じような配色の 2 次元のヒストグラムが描画 されます。2 次元のヒストグラム を 描画した目的 は、先程説明した下記が 実際に成り立つ ことを 視覚的に表現 するためでしたが、点の数を多くしても 描画される図が変わらない のであればその 目的を達成することはできません。同じような図が描画される理由について少し考えてみて下さい。

• 正方形 を 同じ面積になる ように小さく $\boldsymbol{m}$ 個に区切る と、それぞれの範囲 に 点が配置される確率 は $\boldsymbol{\frac{1}{m}}$ になる
• 全体に配置する点の個数 を $\boldsymbol{n}$ とすると、大数の法則 により それぞれの範囲に配置される点を $\boldsymbol{n}$ で割って計算 する 割合 は、$\boldsymbol{n}$ が大きくなればなるほど $\boldsymbol{\frac{1}{m}}$ に近くなる
• $\boldsymbol{n}$ が大きい程 それぞれの範囲に 配置される点の割合 が$\boldsymbol{\frac{1}{m}}$ に近くなる ので、点が均等に配置される可能性が高く なる

その理由は、hist2d が 個数の割合の最小値から最大値までの範囲 の値に 色を割り当てる からです。そのことは、カラーバーの右に表示される数値 が 左上図 の場合は 0.0015 ~ 0.0040 の範囲 に、右上図 の場合は 0.00240 ~ 0.00260 の範囲 になっていることから確認できます。また、前者の範囲の幅が 0.0025、後者の範囲の幅が 0.0002 のように 範囲の幅が 10 倍以上も大きく異なる ことかわかりました。これが 同じような図が描画される原因 で、個数の割合のばらつき比較 する場合は、同じような範囲の割合で比較 する必要があります。

hist2d では、仮引数 vmin と vmax を利用して 色を割り当てる数値の範囲を指定 することができます。下記は先ほどと同じ 10,000、1,000,000 個の点を配置 した場合に、割合の全体の範囲 を表す 0 ~ 1 の範囲 の数値に 色を割り当てた 2 次元のヒストグラムを描画 するプログラムです。

for size in [10000, 1000000]:
    X = np.random.uniform(low=-1.0, high=1.0, size=size)
    Y = np.random.uniform(low=-1.0, high=1.0, size=size) 
    fig, ax = plt.subplots(figsize=(5, 5))
    weights = np.full((size, ), 1 / size)
    h = ax.hist2d(X, Y, bins=20, range=[[-1, 1], [-1, 1]], weights=weights, vmin=0, vmax=1)
    fig.colorbar(h[3])
    ax.set_aspect("equal")

実行結果

　

実行結果から どちらも同じ色で塗りつぶされた図形が描画 されてしまうことがわかります。そのようなことが起きる理由について少し考えてみて下さい。

その理由は以下の通りです。

• 上記のプログラムでは 正方形 を x 軸と y 軸で 20 等分 して区切ったので、区切られた正方形の個数 は 20 × 20 = 400 個 である
• 従って、それぞれの範囲に 点が配置される確率 は $\boldsymbol{\frac{1}{400} = 0.0025}$ である
• 大数の法則 から $\boldsymbol{n}$ が大きくなる と 区切られた範囲 に 配置される点の数の割合 は 0.0025 に近い値 になる
• カラーバー は 一番下が 0、一番上が 1 を表すので 0.0025 に 割り当てられた色 は $\boldsymbol{n}$ の値に関わらず、カラーバーの下から 0.0025 = 0.25 % の位置の色 になる
• 配置される点の割合 は 0.0025 に近い値 なので、全ての色 が カラーバーの下から $\boldsymbol{\frac{1}{400}}$ 付近の位置 の ほぼ同じ色 で塗りつぶされる

実際に一つ前の vmin と vmax を指定しなかった場合の図 から、配置される点の割合 は 0.0015 ~ 0.004 と 0.0024 ~ 0.0026 のように 0.0025 に近い非常に狭い範囲の値 です。また、その 最大値 は 0.004 = 0.4 % なので、カラーバーの色の中 のたった 0.4 % の色しか使われていない ことがわかります。このように 使われる色が極端に少ないことが、ほぼ同じ色で塗りつぶされてしまう原因 です。

使われる色の範囲を増やす方法 の一つは、色を割り当てる数値の最大値 を 配置される点の割合の最大値 とすることです。具体的には下記のプログラムのように、仮引数 vmax に値を代入しない ようにします。また、先ほどと異なり 色が割り当てられる数値の範囲 は 0 以上 0.0025 付近の値 になるので、先程のように 割り当てられる範囲 が 10 倍以上のように 大きく異なるようなことはおきません。実行結果から 配置した点の数が多い右図 の方が 色のばらつきが少ない ことから、より均等に点が配置される ことが 視覚的に確認 できます。

なお、2 つの図の 色が異なる のは 割合の最大値が異なる ため、最も数が多くなる 0.0025 に対応する色が異なる からです。

for size in [10000, 1000000]:
    X = np.random.uniform(low=-1.0, high=1.0, size=size)
    Y = np.random.uniform(low=-1.0, high=1.0, size=size) 
    fig, ax = plt.subplots(figsize=(5, 5))
    weights = np.full((size, ), 1 / size)
    h = ax.hist2d(X, Y, bins=20, range=[[-1, 1], [-1, 1]], weights=weights, vmin=0)
    fig.colorbar(h[3])
    ax.set_aspect("equal")

実行結果

　

配置された点の分散の計算

数値の分布のばらつき は以前の記事で説明した 分散 という指標で表すことができます。hist2d の 返り値の最初の要素 には、下記のプログラムの実行結果のように 分割された範囲 の 個数を表すデータ が代入さるので この値を利用して分散を計算 することができます。

h = ax.hist2d(X, Y, bins=20, range=[[-1, 1], [-1, 1]], vmin=0)
print(h[0])

実行結果

[[0.002466 0.002619 0.002528 0.002483 0.002522 0.002518 0.002472 0.002528
  0.002509 0.002466 0.002531 0.002548 0.002452 0.002461 0.002505 0.002512
  0.002448 0.002547 0.002417 0.002448]
 [0.002533 0.002594 0.002495 0.002513 0.002482 0.002478 0.002467 0.002636
  0.002487 0.002515 0.002474 0.002459 0.002547 0.002449 0.002539 0.002532
  0.002431 0.002515 0.002526 0.002568]
 [0.002446 0.00257  0.002526 0.00252  0.002447 0.002512 0.002556 0.002509
  0.002566 0.002454 0.002499 0.002555 0.002451 0.002486 0.002462 0.002449
  0.002573 0.002504 0.002457 0.002526]
 [0.002463 0.002437 0.002463 0.002454 0.002468 0.002491 0.002418 0.002481
  0.002401 0.002526 0.002539 0.002444 0.00254  0.002526 0.002454 0.002527
  0.002512 0.002431 0.002463 0.002492]
 [0.002507 0.002443 0.002477 0.002506 0.002533 0.002474 0.002457 0.002549
  0.002574 0.002509 0.002425 0.002547 0.002518 0.002446 0.002475 0.002513
  0.002617 0.002518 0.002512 0.002586]
 [0.002454 0.002554 0.002446 0.002452 0.002526 0.002458 0.002526 0.002541
  0.002478 0.002479 0.0025   0.002514 0.002441 0.002469 0.002538 0.002465
  0.002389 0.002462 0.002385 0.002567]
 [0.002441 0.002582 0.002505 0.002495 0.002481 0.002474 0.00237  0.002491
  0.002491 0.0024   0.002513 0.002338 0.002466 0.002506 0.002446 0.002461
  0.00253  0.002502 0.002463 0.002561]
 [0.002568 0.002488 0.002498 0.002507 0.002478 0.002505 0.002533 0.002532
  0.002492 0.002602 0.002537 0.002502 0.002509 0.00247  0.002443 0.002508
  0.002543 0.002431 0.002545 0.00244 ]
 [0.002475 0.002462 0.002449 0.002501 0.002546 0.002452 0.002537 0.002459
  0.002465 0.002503 0.002446 0.002505 0.002531 0.002445 0.002605 0.002496
  0.002482 0.002582 0.002521 0.00248 ]
 [0.002468 0.002541 0.002513 0.002465 0.002602 0.002541 0.002549 0.002469
  0.002497 0.002474 0.002541 0.002445 0.00246  0.002515 0.002459 0.002514
  0.002491 0.002537 0.002593 0.002487]
 [0.002512 0.002469 0.002509 0.002457 0.002463 0.002554 0.002446 0.002481
  0.002537 0.002494 0.002505 0.002419 0.002479 0.002483 0.002449 0.002508
  0.002426 0.002572 0.002567 0.002481]
 [0.002521 0.002535 0.002457 0.00253  0.002575 0.002469 0.002545 0.00252
  0.002422 0.002455 0.00251  0.002525 0.002581 0.002559 0.002506 0.002531
  0.002408 0.002547 0.002496 0.002477]
 [0.002531 0.002505 0.002541 0.002436 0.002526 0.002522 0.002499 0.002522
  0.002447 0.002492 0.002533 0.002484 0.002475 0.002615 0.002527 0.002485
  0.002468 0.002579 0.002466 0.002477]
 [0.002473 0.002559 0.002433 0.002491 0.002574 0.002503 0.002523 0.002527
  0.002467 0.002549 0.002537 0.002572 0.002423 0.002533 0.002399 0.002575
  0.002466 0.002564 0.002537 0.002521]
 [0.002502 0.002536 0.002503 0.002538 0.002496 0.00251  0.002482 0.002494
  0.002516 0.002524 0.002566 0.002556 0.00252  0.002507 0.002455 0.002603
  0.002528 0.002514 0.002492 0.00251 ]
 [0.002484 0.002503 0.002435 0.002509 0.002543 0.002582 0.002518 0.002483
  0.002574 0.002418 0.002481 0.00238  0.002494 0.002468 0.002482 0.002563
  0.002402 0.002543 0.002511 0.00253 ]
 [0.002516 0.002472 0.002451 0.002478 0.00254  0.002426 0.002552 0.002529
  0.002478 0.002521 0.002515 0.002533 0.002494 0.002586 0.002575 0.002544
  0.00249  0.002496 0.002513 0.002499]
 [0.002447 0.002504 0.002555 0.002467 0.002439 0.002501 0.002503 0.002513
  0.002498 0.002522 0.002522 0.002434 0.002551 0.002426 0.002589 0.002506
  0.002515 0.002403 0.002597 0.00253 ]
 [0.002447 0.002545 0.002528 0.002517 0.002538 0.002567 0.00244  0.002493
  0.002491 0.00252  0.002537 0.002531 0.00249  0.002508 0.002484 0.002536
  0.002489 0.002549 0.002475 0.002491]
 [0.002543 0.002472 0.002484 0.002422 0.002458 0.002542 0.002527 0.002447
  0.002439 0.002468 0.002514 0.002496 0.002506 0.002505 0.002518 0.002479
  0.002425 0.002459 0.002451 0.00245 ]]

ndarray や list に記録されたデータの 分散（variance) は numpy の var で計算することができます。下記は 配置した点の数が 10,000 と 1,000,000 の場合に配置された 点の割合の分散を計算 するプログラムです。以前の記事で説明したように e-07 は $\boldsymbol{10^
{-7}}$ を表す指数表記なので、e-09 の方が 1 / 100 程小さい値 を表します。従って、実行結果から 配置した点の数が多い 1,000,000 のほうが 分散が小さくなっている ため、より均等に点が配置される ことが確認できました。

for size in [10000, 1000000]:
    X = np.random.uniform(low=-1.0, high=1.0, size=size)
    Y = np.random.uniform(low=-1.0, high=1.0, size=size) 
    fig, ax = plt.subplots(figsize=(5, 5))
    weights = np.full((size, ), 1 / size)
    h = ax.hist2d(X, Y, bins=20, range=[[-1, 1], [-1, 1]], weights=weights, vmin=0)
    print(np.var(h[0]))

実行結果（表示される画像は分散とは関係ないので省略します）

2.5469999999999923e-07
2.5100100000006664e-09

numpy の var の詳細については下記のリンク先を参照して下さい。

仮引数 cmap の利用方法

これまでの hist2d では 数値が小さい場合は青色 に、大きい場合は黄色 が 割り当てられました が、異なる配色を割り当てる こともできます。数値と色の対応表 のことを カラーマップ と呼び、自分でオリジナルのカラーマップを作成することもできますが、matploblib の cm モジュール で あらかじめ定義されているカラーマップを利用 するのが 一般的 ではないかと思います。あらかじめ定義されているカラーマップ には 名前がついており、下記のプログラムのようにその 名前を表す文字列を仮引数 cmap に代入 することで指定することができます。"jet" は 小さい数値を青色 で、大きな数値を赤色 で表すカラーマップです。

X = np.random.uniform(low=-1.0, high=1.0, size=size)
Y = np.random.uniform(low=-1.0, high=1.0, size=size) 
fig, ax = plt.subplots(figsize=(5, 5))
weights = np.full((size, ), 1 / size)
h = ax.hist2d(X, Y, bins=20, range=[[-1, 1], [-1, 1]], weights=weights, vmin=0, cmap="jet")
fig.colorbar(h[3])
ax.set_aspect("equal")

実行結果

カラーマップとしては "jet" がわかりやすく、良く使われると思いますので 本記事 では "jet" を利用 することにします。

matploblib の cm モジュール の詳細と、あらかじめ定義されている カラーマップの一覧 については下記のリンク先を参照して下さい。

calc_mc_pi の改良

上記で説明した内容で calc_mc_pi を改良 することにします。下記はそのプログラムです。最後に 2 次元のヒストグラムの描画 と 分散の計算処理を追加 しただけなので説明と修正箇所は省略します。

def calc_mc_pi(num):
    print(f"点の個数 = {num}")
    X = np.random.uniform(low=-1.0, high=1.0, size=num)
    Y = np.random.uniform(low=-1.0, high=1.0, size=num)    
    I = X * X + Y * Y <= 1
    i = np.count_nonzero(I)
    mc_pi = 4 * i / num
    print(f"mc_pi = {mc_pi}")
    print(f"np.pi = {np.pi}")
    print(f"誤差  = {mc_pi - np.pi}")     

    X_inner = np.ma.array(data=X, mask=~I)
    X_outer = np.ma.array(data=X, mask=I)
    fig, ax = plt.subplots(figsize=(5, 5))
    r = patches.Rectangle(xy=(-1, -1), width=2.0, height=2.0, fill=False, ec="k")
    c = patches.Circle(xy=(0, 0), radius=1, fill=False,ec="r")
    ax.add_patch(r)
    ax.add_patch(c)
    ax.scatter(X_inner, Y, c="r")
    ax.scatter(X_outer, Y, c="k")
    
    fig, ax = plt.subplots(figsize=(5, 5))
    weights = np.full((num, ), 1 / num)
    h = ax.hist2d(X, Y, bins=20, range=[[-1, 1], [-1, 1]], weights=weights, vmin=0, cmap="jet")
    fig.colorbar(h[3])
    ax.set_aspect("equal")
    print(f"分散  = {np.var(h[0])}")

下記は先ほどと同様に 10、100、1,000、10,000、100,000、1,000,000 個 の点を 配置して円周率の近似値を計算 してグラフなどを描画するプログラムです。描画される画像は後でまとめます。実行結果の計算結果は先ほどと同様になるので省略します。

calc_mc_pi(100)for num in [10, 100, 1000, 10000, 100000, 1000000]:
    print(f"num = {num}")
    calc_mc_pi(num)

実行結果

点の個数 = 10
mc_pi = 3.6
np.pi = 3.141592653589793
誤差  = 0.458407346410207
分散  = 0.00024374999999999993
点の個数 = 100
mc_pi = 2.96
np.pi = 3.141592653589793
誤差  = -0.18159265358979315
分散  = 2.775e-05
点の個数 = 1000
mc_pi = 3.26
np.pi = 3.141592653589793
誤差  = 0.11840734641020667
分散  = 2.3950000000000004e-06
点の個数 = 10000
mc_pi = 3.1496
np.pi = 3.141592653589793
誤差  = 0.00800734641020684
分散  = 2.4224999999999925e-07
点の個数 = 100000
mc_pi = 3.14576
np.pi = 3.141592653589793
誤差  = 0.004167346410206996
分散  = 2.387050000000012e-08
点の個数 = 1000000
mc_pi = 3.141392
np.pi = 3.141592653589793
誤差  = -0.0002006535897929318
分散  = 2.433290000000646e-09

下記は上記の実行結果をまとめた表です。誤差は小数点以下第 6 桁で、分散は誤差は小数点以下第 10 桁で四捨五入しました。表から 点の数が多くなる と 誤差と分散が小さくなる傾向がある ことが確認できました。

点の数	円周率の近似値	誤差	分散
10	3.6	0.45841	0.0002437500
100	2.96	-0.18159	0.0000277500
1,000	3.26	0.11841	0.0000023950
10,000	3.1496	0.00801	0.0000002422
100,000	3.14576	0.00417	0.0000000239
1,000,000	3.141392	-0.00020	0.0000000024

下記は描画された 2 次元のヒストグラムの画像をまとめたもの です。点の配置の画像は先ほどとほぼ同じになるので省略します。画像から、点の数が多くなるほど 配置された 点の数の割合を画像の表す色 が より均一になる 事から、配置された点の数の 割合のばらつきが減り、正方形の中に 均等に点が配置されるようになることが確認 できます。

点の数	画像	点の数	画像
10		100
1,000		10,000
100,000		1,000,000

今回の記事のまとめ

今回の記事では、乱数を利用した計算 を行うことで 近似値を求める アルゴリズムである モンテカルロ法 について説明し、円周率の近似値 を モンテカルロ法で計算する方法 を紹介しました。次回の記事では大数の法則について説明し、原始モンテカルロ法で AI を作成する方法について説明する予定です。

本記事で入力したプログラム

リンク	説明
marubatsu.ipynb	本記事で入力して実行した JupyterLab のファイル

次回の記事

1. モンテカルロ法を利用したゲームの AI のアルゴリズムにはいくつかの種類があります。今回の記事で紹介するのはその中で最も単純なものなので原始モンテカルロ法と呼ばれるようです ↩

2. 詳細は省略しますが、その場合に 1 が計算される確率は $10^{-10}$ 以下です ↩

3. この式の場合は円周上の点は円の内部にあるとみなされます。円周上の点を円の内部とみなしたくない場合は $x^2 + y^2 < 1^2 = 1$ とする必要があります。ただし、実際に円周上の点が計算される確率は極めて 0 に近いのでどちらでも円周率の計算結果に影響はありません ↩

4. 円周率は無理数で小数点以下に無限の数値が並ぶ値なので、np.pi は正確な円周率の値ではありません。np.pi は円周率の小数点以下 16 桁目を四捨五入した値です ↩

5. x 座標と y 座標を区切る個数を別々に指定することもできます。詳細は hist2d のドキュメント を参照して下さい ↩