今回は、Pythonで簡単なランダム実験のモンテカルロシミュレーションを行う方法について学習します。
注意：モンテカルロ・シミュレーションは数学的に複雑な分野です。そのため、私たちはMCの詳細には踏み込んでいません。その代わり、モンテカルロシミュレーションの必要性と実装を理解するために、いくつかの直感と例を用いて、数学的な背景があまりない人でも、数学的なものをあまり使わずに確率に触れることができるようにしました。
Python でモンテカルロシミュレーション
まず、Randomライブラリを使ってコイントス実験のシミュレーションを行い、モンテカルロ実験への勘所を構築します。
1. random モジュール
まず、random モジュールをインポートします。
# Import the random moduleimport random
|
ここでは、random モジュールの uniform 関数を広範に使用します。この関数はユーザが指定した下限と上限の間の数を与える。上限と下限の間の各数値の出現確率は等しくなります。
# Generates uniform random number between 4 and 6random.uniform(4, 6)
|
を結果を出力すると、以下の様になります。
5.096077749225385 |
さて、この一様関数を用いて簡単なコイン投げをシミュレートしてみましょう。下限と上限の間の各数字の出現確率は等しいとお伝えしました。
つまり、0と1の間の一様な値をとれば、ある数字が0.5より大きくなる確率も0.5より小さくなる確率も等しくなるのです。この特徴を利用します。0から1までの乱数をとり、その値が0.5より大きければ、結果は表となり、そうでなければ裏となる。
a = random.uniform(0, 1)
if a>0.5:
print("Head")
else:
print("Tail")
|
結果は、以下の通りになります。
Head |
2. 偏りのないコイントスをシミュレートする関数を定義する
さて、前節で学んだことを活かして、偏りのないコイントスをシミュレートする関数を書いてみましょう。関数を書くことで、コードがより読みやすく、よりモジュール化されます。
def unbiased_coin_toss():
# Generate a random number from 0 to 1
x = random.uniform(0, 1)
# Probability that the number falls between 0 and 0.5 is 1/2
if x > 0.5:
# Heads for True
return True
else:
# Tails for False
return False
|
この関数の結果をテストしてみましょう。
for i in range(10):
print(unbiased_coin_toss())
|
出力してください。
FalseTrueFalseFalseFalseFalseTrueFalseFalseFalse |
3. 少ない回数でコインを投げる
さて、これで本物のコイン投げをシミュレートできました。では、一連のランダムなコイン投げで、表が出る確率を検証してみましょう。具体的には、呼び出すたびに表か裏が出るような関数を定義します。
さて、コインを何回か投げて、その結果をリストに格納します。表が出る確率はリストから計算される。
N=10
# List to store the result(True/False)results = []
# Toss the coin 10 times and store that in a listfor i in range(N):
result = unbiased_coin_toss()
results.append(result)
# Find the total number of headsn_heads = sum(results)
# Find probability of heads in the experimentp_heads = n_heads/N
print("Probability is {:.1f}".format(p_heads))
|
出力せよ。
Probability is 0.9 |
おっとっと! これはなかなかうまくいかなかったようです。このブロックは何度でも実行できますが、この実験では頭の確率が期待される0.5から大きく変化していることがわかります。
シミュレーションに問題はありませんか？
実を言うと、「はい」と「いいえ」の両方があります。先ほど定義した関数が完璧に動作しなかったために、このような不完全な結果になったと思うかもしれません。実際の問題は、シミュレーションの方法にあるのです。
大数の法則により、実験回数が多いと実験確率は実際の確率と期待される確率に近くなります。
上記の説明は少し変な感じがします。これを検証するための数学的な証明や仮説には踏み込まず、簡単な直感に基づく考え方をします。
例えば、インドで小麦が消費される確率を求める仕事を与えられたとします。理想的には、あなたは一人一人に会いに行き、小麦を消費しているかどうかを尋ねる必要があります。小麦を消費する確率は
prob = []
# Make 1000 experimentsfor i in range(1000):
# Each experiment have 10 coin tosses
N = 10
results = []
# Toss the coin 10 times and store that in a list
for i in range(N):
result = unbiased_coin_toss()
results.append(result)
n_heads = sum(results)
p_heads = n_heads/N
prob.append(p_heads)
# average the probability of heads over 1000 experimentsp_heads_MC = sum(prob)/1000
print("Probability is {:.3f}".format(p_heads_MC))
|
しかし、13億人に聞くのは面倒な作業です。そこで、その国の全人口を代表する100人を選んで、実験をするのです。確率を求める作業は、かなり楽になります。しかし、そうでしょうか？
もし、パンジャブ州のような小麦消費量の多い州から多くの人を集め、西ベンガル州のような小麦消費量の少ない州から少ない人を集めたり、その逆の場合は、実験確率がかなりずれてしまうかもしれません。
これは、実験のために無作為に選んだ100人が、母集団全体を正しく代表することができないために起こることです。ですから、結果は常に誤差を含みやすいのです。
コイン投げゲームも同じような考え方です。コイン投げの回数が足りなかったので、性急に解決に至ってしまいました。これを直そう!
Pythonでモンテカルロシミュレーションを行う。
モンテカルロシミュレーションは、この問題を回避するための最良の方法の一つです。
簡単に言うと、モンテカルロシミュレーションでは、異なる入力値から始まる値から異なる実験結果を取り出し、その結果を平均化（期待値）します。
その結果得られた平均値が、ここで求めているエラーフリー（エラーが少ない）な答えとなるのです。
Probability is 0.502 |
結果は、以下の通りになります。
このコードのブロックを実行するたびに、期待値に非常に近い確率の値が得られる。我々は実験回数を増やすことで精度を上げている(1000)。この数を増やして、自分自身で結果をテストすることができる。
結論
コイントスのモンテカルロ・シミュレーションについて、今回で終わりです。0.5ではなく、頭の確率が偏ったコインでモンテカルロシミュレーションを試してみてください。実は、他のランダムな実験でもモンテカルロ・シミュレーションを試して結果を得ることができます。