モンテカルロ法で円周率πと任意の図形の面積を求める方法をExcelで実演！
円周率を求める＝円の面積を求める
モンテカルロ法を理解するための入門事例として最も有名なのは、円周率πの近似値を求める問題です。
ここではそのやり方を解説した後、変な図形の面積の近似値を求める問題に応用していきます。
円周率を求めるということは、実は円の面積を求めているため、同じやり方でどんな図形の面積でも求められるのです。
図形を表す式は分かっているという「但し付き」ではありますが。
もし図形が手書きで式が分からないけど面積を求めたいという場合は、こちらを参照下さい。
モンテカルロ法で円周率πを求める方法
円周率πは円の面積を求める時に登場します。
円の面積＝半径×半径×π
でしたね。
従って、半径が1なら円の面積は円周率πそのものになります。
円に接する正方形内にランダムに多くの点を打つ
ここに、下図のように原点に中心がある半径が1の円があるとします。
この円に接する正方形は1辺が2になりますので、面積は4です。
この正方形の中にランダムに点を沢山打っていきます。
これをExcelでやらせるには、RAND関数を使います。
RAND関数は、「=RAND()」と入力するだけで、0から1までの乱数を作ってくれます。
>> エクセルを使って正規分布の【疑似乱数】を生成する方法を具体例で実演！
この関数だけでは0から1の乱数しかできませんが、2倍すれば0から2の乱数ができます。
欲しいのは、x、y共に-1から+1の乱数ですから、-1に0から2の乱数を足せば良いでしょう。
つまり、x、y共に
=-1+RAND()*2
と入力すれば、-1から+1の乱数が作れます。
Excelシートで10,000個の乱数を作ると、次のようになります。
できた乱数を散布図にすると、次のようになります。
円の中に入っている点の数を数える
正方形の内側に10,000個の点がランダムに入りますので、ほぼ塗りつぶされます。
この10,000個の点のうち、半径1の円に入っているのはいくつあるでしょうか？
それは、それぞれの点について√(x2+y2)を計算してみれば分かります。
原点からの距離です。
√(x2+y2)<=1
を満たす点はすべて半径1の円に入っているはずです。
Excelでは次のように計算できます。
円の中に入った点の割合に正方形の面積を掛ければそれが円の面積
E列に円に入っている点には1が入力されます。
例えば、この1が7,861個あったとしましょう。
すると、10,000個の点のうち78.61%が円の中に入ったことになります。
正方形の面積は4ですから、
4×78.61%＝3.1444
が円の面積になります。
つまり、π＝3.1444です。
点を多く打つほど精度は上がる
以上が、モンテカルロ法による円周率πの求め方になりますが、当然ながら乱数が変わる度にπの値は変わります。
また、今回は10,000個の乱数で計算しましたが、より多くの乱数で計算するほどπの精度が上がることが予想されます。
これを検証してみましょう。
10,000個、1,000個、100個の乱数で計算した時の計算結果は、次のようになりました。
やはり、乱数の数が多い方が精度は高いように見えます。
しかしこれもたった1回の試行の結果に過ぎません。
100回行った場合の結果は、次のようになりました。
100回繰り返せば乱数の数に拘わらず大体同じ値になりますが、ばらつきは乱数の数を増やせば小さくなることが分かります。
モンテカルロ法で任意の関数で囲まれた図形の面積を求める方法
以上が円周率を求めるやり方ですが、結局は円の面積を求めています。
ですので、このやり方は一般的な図形の面積を求めることに応用できます。
その図形の式さえ分かっていれば、どんな図形の面積も求めることができます。
面積を求めたい図形を囲む長方形内にランダムに多くの点を打つ
例題として、次の2つの関数で囲まれた部分の面積を求めてみましょう。
ここで、
f(x)＝{sin(30x2)/x – cos(20x)/ex + √(5x)}2
g(x)＝-40x + 40
です。
【Excelで機械学習】モンテカルロ積分①難しい積分が誰でも簡単にExcelで解ける
でf(x)とx軸で囲まれた部分の面積を求めましたが、今回のやり方だと任意の関数同士で囲まれた部分の面積を求められます。
f(x)はとても複雑な関数ですが、先ほどの円周率と同じやり方で解けます。
まず、求めたい部分を囲うような長方形の領域を決めます。
例えば、下図で影掛けした長方形に決めます。
この長方形はx座標は0から1、y座標は0から40なので、
x：0～1の乱数
y：0～40の乱数
を沢山作れば、長方形の内側に満遍なく点を打つことができます。
面積を求めたい図形の中に入った点の割合に長方形の面積を掛ける
そして、y座標がf(x)より小さくてg(x)より大きい点の数を数えれば、それが求めたい部分に入る点の数です。
その点の数を全体の点の数で割れば、求めたい部分の全体面積に対する比率が求まります。
全体面積は40ですので、40×比率が求めたい部分の面積です。
これをExcelで表現すると、次のようになります。
この面積の正解は分かりませんが、乱数10,000個で計算すると0.8になりました。
多分、大体これくらいなのでしょう。
でも、実はこれ、大変効率が悪いのです。
F列の1の数を数えると、たったの200しかありません。
つまり、生成した10,000個の乱数のうち、2%しか求める部分の図形に入らなかったのです。
この比率は高ければ高いほど効率が良く、また計算精度も高くなります。
点は数多く打つほど精度は上がる
参考までに、乱数の数を10,000個、1,000個、100個それぞれの場合について100回ずつ試行してみた結果が次の通りです。
先ほどの円周率を求める場合と比べて、随分と精度が落ちていることが分かります。
円周率の場合は約80%の乱数が求めたい図形の中に入ったのに対し、今回は約2%しか入らなかったためです。
このようにモンテカルロ法によりどんな図形の面積も求めることはできますが、精度は落ちることがあります。
しかしその場合でも、乱数を沢山打つことによりカバーすることができます。