モンテカルロ法は、与えられた関数を再現するように、関数の積分をサンプリングした点の和に変換する方法です。
Wikipedia の説明にもあるように、これは乱数生成の品質に依存します。
モンテカルロ法の原理
モンテカルロ法による数値計算の原理は、アンケートにおいて、多くの人員が所属する部署からは多数のサンプルを抽出し、少ない人員が所属する部署からは少数のサンプルを抽出し、これで会社全体の意見をアンケートに反映させることができるという考え方に似ています。たとえば円周率の計算については 0 から 1 の範囲の乱数の組を生成して、それを x-y 座標上の点としてみて単位円の第 1 象限に入っているかどうか、その割合から求めるといった具合です。
連続的な値を取る確率変数の場合、確率分布は確率密度関数 f(x) で表現されます。このとき平均値 u と分散 σ^2 は次のようになります。
\mu = \int{xf(x)dx} \\
\sigma^2 = \int{(x - \mu)}^2f(x)dx
実際にコードで試す
2 つのスポーツチーム A と B がありました。この 2 つのスポーツチームは成績はまったく同点でしたがその分散が異なります。これをモンテカルロシミュレーションで表現します。
from pylab import *
from scipy.stats import *
import matplotlib.pyplot as plt
runs = 10000
# 各チームの成績
teamperfA = [0.9,0.8,0,9] # チーム A のスコア
teamperfB = [0.9,0.8,0.9] # チーム B のスコア
# チームごとの成績の分散
teamvarianceA = [0.03, 0.04, 0.02]
teamvarianceB = [0.05, 0.08, 0.09] # チーム B のほうが分散が高い
weights = [1.0, 0.95, 0.8]
def result(perf,variance,weights):
res = 0.0
for i in range(len(perf)-1):
res += perf[i] * weights[i] * norm(1,variance[i]).rvs()
return res
resultsA = zeros(shape=(runs,), dtype=float)
resultsB = zeros(shape=(runs,), dtype=float)
for i in range(runs):
resultsA[i] = result (teamperfA,teamvarianceA,weights)
resultsB[i] = result (teamperfB,teamvarianceB,weights)
subplot(211)
width = 2
height=runs
title('Team A')
hist(resultsA, bins=50)
axis([1.4,1.9,0,height/10])
subplot(212)
title('Team B')
hist(resultsB, bins=50)
axis([1.4,1.9,0,height/10])
show()
plt.savefig("image.png")
参考
Python matplotlib, Monte Carlo simulation, and basic statistics
http://softwaredevelopmentperestroika.wordpress.com/2013/12/06/python-matplotlib-monte-carlo-simulation-and-basic-statistics/