積分法の数値計算をプログラミングしてみよう：数学×Pythonプログラミング入門（3/5 ページ）

[羽山博，著]

目標3： モンテカルロ法により定積分の近似計算を行う

　積分などの近似計算を行う場合には、乱数を利用したモンテカルロ法と呼ばれる方法も使われます。同じ例ばかりで退屈かも知れませんが、結果が簡単に確認できるので、ここでも、

で考えてみましょう。

　モンテカルロ法の考え方は極めてシンプルです。この場合であれば、0 ≤ x < 1, 0 ≤ y < 1 の範囲で乱数のペア(x,y)を作ります。そのペアで表される点には以下のような2通りの場合があります。具体例で見てみましょう。

• (x, y) ＝ (0.7, 0.2)の場合
  　　→ x²＝0.49なので、この点のy座標はx²の曲線よりも下にある
  　　　　 → アミカケの部分に入っている
• (x, y) ＝ (0.4, 0.8)の場合
  　　→ x²＝0.16なので、この点のy座標はx²の曲線よりも上にある
  　　　　→ アミカケの部分には入っていない

　この2つの点をグラフ上に描画したのが図4です。

図4　モンテカルロ法により面積を求めるための考え方
モンテカルロ法では、乱数のペアを多数作成し、アミカケの部分に入った個数÷全体の個数で面積を近似する。

　このようにして、多数の乱数を作っていき、アミカケの部分に入った個数を数え、全体の個数で割れば、アミカケの部分の割合が分かります。つまり、「y < x²なら、カウントを増やす」という処理を繰り返すわけです*2。この例では0 ≤ x < 1, 0 ≤ y < 1 の四角い範囲の面積が1なので、アミカケの部分の割合がそのまま面積となります。

*2　ここで利用する乱数は、範囲内のどの値も同じ程度の確からしさで現れる「一様乱数」です。

　これは、0 ≤ x < 1, 0 ≤ y < 1の四角い的に向かって、ダーツの矢を数多く投げて、アミカケの部分に突き刺さった個数を数えるのと似ています（そういう状況を脳裏に思い浮かべてみてください）。ダーツの矢はランダムに投げますが、必ず四角い範囲のどこかに突き刺さるものとすれば、「アミカケの部分に突き刺さった矢の数÷全体の矢の数」が、アミカケの部分の割合に近くなるというわけです。

3. モンテカルロ法による積分をプログラミングする

　では、プログラミングに取りかかりましょう（リスト12）。乱数はrandomモジュールを使って作成することもできますが、NumPyのrandomモジュールでも作成できるので（その方が何かと便利なので）、そちらを利用します。

import numpy as np

def monte(n):
  count = 0
  for i in range(n):
    x = np.random.random() # 0以上1未満の乱数を作成する
    y = np.random.random()
    if y < x**2:
      count += 1
  return count / n

リスト12　モンテカルロ法により面積の近似値を求める関数monteの定義
引数nには、発生させる乱数のペアの個数を指定する。NumPyのrandom.random関数は0以上1未満の一様乱数を作成するための関数。x²よりもyの値が小さければ、二次曲線よりも下にあるので、個数countを増やす。最後にcountを全体の個数nで割った値を返す。

　実行例は以下の通りです。

print(monte(1000))
# 出力例：0.335
print(monte(10000))
# 出力例：0.3382
print(monte(100000))
# 出力例：0.33369

リスト13　モンテカルロ法で面積の近似値を求めた結果
乱数を使っているので、結果は毎回変わるが、10万回の繰り返しで0.33369となった。

　乱数は毎回異なるものが使われるので、結果は実行するたびに異なります。毎回同じ結果を得たいのであれば、seedを決めておいてやるといいでしょう。seed（種）とは乱数を作成する際の初期設定のために使われる値のことです（seedを指定しない場合はシステムの時刻が使われるので、毎回異なる並びの乱数が作られます）。以下のコードを、リスト12のcount=0の次に書いてみてください。

np.random.seed(0)

　このように書いておけば、毎回同じ乱数が作成されます（値はランダムですが、値の並びは毎回同じになるということです）。その場合の実行例は以下のようになります。

print(monte(1000))
# 出力例：0.327
print(monte(10000))
# 出力例：0.3336
print(monte(100000))
# 出力例：0.33347

リスト14　乱数の初期値を固定した場合の実行結果
乱数のseedを固定すると、毎回同じ結果が求められる。サンプルプログラムと同じ結果になるか動作確認する場合などに便利。

モンテカルロ法の可視化

　せっかくなので、アミカケの部分に入った点を可視化しておきましょう。 y < x²のときに1つずつ点を打つと処理に時間がかかるので、y < x²となる全ての点をPythonのリストに記録しておき、そのリストを使って散布図を作成します。リスト15の関数monte_mkdataは、アミカケの部分に入る全ての点の座標を返す関数です。

import numpy as np

def monte_mkdata(n):
  np.random.seed(0)
  px = []
  py = []
  for i in range(n):
    x = np.random.random()
    y = np.random.random()
    if y < x**2:
      px.append(x)
      py.append(y)
  return px, py

リスト15　モンテカルロ法を可視化するためのデータを作成する関数monte_mkdata
引数nには繰り返し数を指定する。y < x²となるxとyをそれぞれpx、pyに追加していく。返り値はpxとpy。この関数はデータを作成するだけで、面積の近似値は求めていない。

　関数monte_mkdataを使ってデータを作成し、散布図を描いてみましょう（リスト16）。

import matplotlib.pyplot as plt

# 散布図を描く
plt.scatter(*monte_mkdata(10000), s=1, c="red")

# y=x**2のグラフを重ねて描く
xrange = np.arange(0, 1, 0.01)
plt.plot(xrange, xrange**2)

plt.show()

リスト16　モンテカルロ法を可視化する
scatter関数の最初の引数がmonte_mkdataによって作成した点の座標。s=1は点のサイズを1とするという指定。c="red"は点の色の指定。y＝x²のグラフも合わせて描いておいた。

　ここでは、monte_mkdata(10000)の前に書かれている*に注意してください。これは、リストやタプルの要素を展開するという意味です。つまり、monte_mkdata(10000)が返した(x, y)というタプルを個別の要素であるxとyに分け、plt.scatterの引数としてxとyをそれぞれ指定することになります。実行例は以下の通りです。

図5　モンテカルロ法によって作成された点をプロットする
プロットした全ての点が、アミカケの部分の面積をほぼ覆い尽くしていることが分かる。

コラム　NumPyの機能を活用して繰り返し処理を排除する

　NumPyのrandom.random関数に引数を指定すると、指定した個数の乱数をNumPyの配列（ndarray）として作成してくれます。その機能を利用すれば、繰り返し処理を使わずにコードを書くこともできます（リスト17）。

import numpy as np

def monte_np(n):
  np.random.seed(0)
  x = np.random.random(n) # n個の乱数
  y = np.random.random(n)
  result = y < x**2 # TrueかFalseかの配列が返される
  return np.count_nonzero(result) / n # Trueの個数を全体の個数で割る

リスト17　乱数を全て作成し、x**2よりも小さいかを一気に調べる
まずn個の乱数を作成し、それをxとする。次にn個の乱数を作成し、それをyとする。resultには、y < x**2の全ての結果（TrueまたはFalse）が配列として返される。Trueの数を数えるには、count_nonzero関数が便利。この関数はゼロでない値の個数を数えるものだが、Falseはゼロと見なされるので、Trueの個数が数えられる。

　ただし、繰り返し処理を使った例（リスト13）では、乱数を、

　　xの値、yの値、次のxの値、次のyの値……

という順に作成していましたが、上のコードでは、

　　xの値をn個作成した後、yの値をn個作成している

ので、実行結果は異なります（例えば、monte_np(1000)を実行すると、結果は0.331になります）。リスト13の実行例（リスト14）と完全に同じ結果（0.327）が欲しければ、リスト18のように、乱数を2n個作成し、偶数番目をxに、奇数番目をyに代入するといいでしょう。

import numpy as np

def monte_np(n):
  np.random.seed(0)
  data = np.random.random(2*n) # 2n個の乱数を一気に作る
  x = data[0::2] # 0番から増分値2で順に取り出す
  y = data[1::2] # 1番から増分値2で順に取り出す
  result = y < x**2
  return np.count_nonzero(result) / n

リスト18　スライス機能を使って交互に値を取り出す
NumPyの配列では、[start:stop:step]の形式で、部分配列を取り出すことができる（スライス機能と呼ばれる）。startを省略すると0が指定されたものと見なされ、stopを省略すると末尾までが指定されたものと見なされる。stepを省略すると1が指定されたものと見なされる。

　とはいえ、本来の目的は全く同じ結果を得ることではなく、あくまでも近似値を得ることなので、特段の理由がなければ、リスト17のままで構いません。また、リスト18でx = data[0:n]; y = data[n:]として前半の乱数をxに代入し、後半の乱数をyに代入しても構いません。

　モンテカルロ法も簡単でしたね。では、別の関数を使った例も見てみましょう。次は積分というよりは、指定された分布に従ってデータを取り出す（サンプリングする）ことを目的とした例です。もちろん、サンプリングされたデータを使って、面積（確率）を求めることもできます。