モンテカルロ法

0.はじめに

本記事は学校の課題でレポートを作成する際に、「scrapboxで作れば楽なんだよなぁ〜」という怠惰の下に作られたものです。

なので、構成や言葉遣いなどが記事っぽくない箇所があります。

モンテカルロ法の説明記事なんて星の数ほどあるとは思うんですが、まぁ、星がひとつやふたつ増えるくらい良いと思うんで本記事もスペースデブリとしてモンテカルロの地中海に埋もれさせて下さい。

1.概説

2.説明

3.数値解析

数値解析の分野においてのモンテカルロ法は「確率を近似的に求める手法」としてしばしば用いられる。乱択的に結果を出力し、その総計の出力回数と解析対象の出力回数を用いて確率を示す。

また、この確率を用いて積分での近似的な解を求めることも可能である。

https://upload.wikimedia.org/wikipedia/commons/thumb/8/84/Pi_30K.gif/220px-Pi_30K.gif

4.精度

4.1. 試行回数

モンテカルロ法の精度はその試行回数にも依存する。ある程度の精度でいいのであれば試行回数に問題はないと思われるが、一定の精度を求める場合には試行回数が膨大となる。なお、この試行回数の増加は指数的であり(重複対数の法則)、これは大数の弱法則からわかる。なお、この試行回数と精度が比例的になるように乱数を一様ではないようにする手法もあり、そのようなモンテカルロ法を別に準モンテカルロ法と呼ぶ。

4.2. 乱数生成

5.円周率における数値解析

モンテカルロ法における円周率の近似はある程度の精度しかでない。これは上記でも述べたようにある程度の精度を望む場合にはその試行回数が膨大になることに起因するが、それ以前に円周率の精度に対してモンテカルロ法が膨大な試行回数を必要となるまでの精度に至るのが早いことも原因となっている。これに対して、本項ではモンテカルロ法以外の円周率の近似を求められる手法を紹介する。

ライプニッツの公式

$ 1 - \frac{1}{3} + \frac{1}{5} - \frac{1}{7} + \frac{1}{9} ... = + \frac{π}{4}

マチンの公式

$ π = 16 \sum_{n=0}^{3m+2}{\frac{(-1)^n}{2n+1}(\frac{1}{5})^{2n+1}}-4 \sum_{n=0}^{m}{\frac{(-1)^n}{2n+1}(\frac{1}{239})^{2n+1}}

ラマヌジャンによる公式

$ \frac{1}{π} = \frac{2\sqrt{2}}{99^2} \sum_{n=0}^{∞}{\frac{(4n)!(1103+26390n)}{(4^n 99^n n!)^4}}

6.プログラムコード

モンテカルロ法による「事前に選択した関数」の積分による面積を求めるプログラムコードを示す。

code:MonteCarlo.py

import random

import math

def method(x):

ret = math.cos(x)

return ret

def random_method(_part,_min,_length,lcgs):

if _part == 0:

return random.uniform(_min,_length+_min)

elif _part == 1:

return (((48271*lcgs)%((2**31)-1))/((2**31)-1))*_length

_width = int(input("Please input a range of X."))

_height = int(input("Please input a range of Y."))

_small_width = int(input("Please input the minimum value of X."))

_small_height = int(input("Please input the minimum value of Y."))

_trials = int(input("How many trials do you want?"))

_part_random = int(input("Please input.\n 0.python method\n1.Linear congruential generators"))

_lcgs_x = (random_method(_part_random,_small_width,_width,_width)/2)*((2**31)-1)

_lcgs_y = (random_method(_part_random,_small_height,_height,_height)/2)*((2**31)-1)

_trials_check = 0

for num in range(_trials):

_random_x = random_method(_part_random,_small_width,_width,_lcgs_x)

_lcgs_x = (_random_x/_width)*((2**31)-1)

_output_y = method(_random_x)

_random_y = random_method(_part_random,_small_height,_height,_lcgs_y)

_lcgs_y = (_random_y/_height)*((2**31)-1)

if _output_y >= _random_y:

_trials_check = _trials_check +1

print(_width*_height*_trials_check/(num+1))

print(str("lcg = ") + str(_random_y) + " " + str(_output_y))

次に上記で説明したライプニッツの公式、マチンの公式、ラマヌジャンによる公式をそれぞれ用いて円周率の近似を求めるプログラムコードを示す。

code:pi.py

import math

_trials = int(input("How many trials do you want?"))

_leibniz = 0

_machin_a = 0

_machin_b = 0

_ramanujan = 0

for num in range(_trials):

_leibniz = _leibniz + ((-1)**(num+2))*(1/((2*num)+1))

_machin_b = _machin_b + (((-1)**num)/((2*num)+1))*(1/239)**((2*num)+1)

_ramanujan = _ramanujan + (math.factorial((4*num))*(1103+(26390*num)))/(((4**num)*(99**num)*math.factorial(num))**4)

for num in range((_trials*3)+2):

_machin_a = _machin_a + (((-1)**num)/((2*num)+1))*(1/5)**((2*num)+1)

print("Leibniz = " + str(_leibniz*4))

print("Machin = " + str((16*_machin_a)-(4*_machin_b)))

print("Ramanujan = "+ str(1/(((2*(2**0.5))/(99**2))*_ramanujan)))