モンテカルロ法による円周率計算

　pythonを使ったモンテカルロ法による円周率計算って初学者の演習によく用いられるので、もうみなさんにやり尽くされています。でも、自分はまだやってないのでやりましたという記事です。ただ、円周率を計算するだけでは面白くないので、計算のパラメタを弄った結果に統計的な解析を加えてみます。

計算方法

　一辺が$r$の正方形の中にランダムに点を打つと、ある頂点からの距離が$r$以下の点の確率は$\frac{\pi}{4}$のになるということを用います。これは一辺$r$の正方形と半径$r$の円の面積比が$\pi$になるということから導出できます。ランダムに打った点の数を$N$、原点からの距離が$r$以下の点の数を$p$とすると次の式が成立します。

$$\pi = 4\frac{p}{N}$$

このように、ランダム変数を用いて行うシミュレーション方法をモンテカルロ法と呼びます。

実装

ランダムに$N$点打って$p$を数えるプログラムを書きましょう。今回は計算値の試行回数依存性も見たいのでfor文で繰り返し処理も行います。計算時間の計測や計算結果の図示も含めます。

##円周率の近似値を求めよう！
import random
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import time

N = 1000 #打点数
S = 100 #試行回数

x_all_list = []
y_all_list = []
x_all_outlist = []
y_all_outlist = []
pi_all_list =[]
pi_list = []
pi_std_list = []

start = time.time()

for i in range(S):
    x_list = []
    y_list = []
    x_outlist = []
    y_outlist = []
    
    p = 0
    
    for n in range(N):
        
        xn = random.randint(0, 10e20)/10e20 #[0,1]の乱数
        yn = random.randint(0, 10e20)/10e20
        
        r = np.sqrt(xn**2 + yn**2)
        
        if r <= 1:
            p += 1
            
            x_outlist.append(1.1)
            y_outlist.append(1.1)
        else:
            x_outlist.append(xn)
            y_outlist.append(yn)
        
        
        x_list.append(xn) # xの値をx_listに格納していく
        y_list.append(yn) # yの値をx_listに格納していく
    
    pi = 4 * p / N
    
    pi_all_list.append(pi)
    
    pi_list.append(np.sum(pi_all_list)/(i+1))
    
    pi_std_list.append(np.std(pi_all_list))
    
    x_all_list.extend(x_list)
    y_all_list.extend(y_list)
    x_all_outlist.extend(x_outlist)
    y_all_outlist.extend(y_outlist)

elapsed_time = time.time() - start
print ("elapsed_time:{0}".format(elapsed_time) + "[sec]")

def reshepe_to_plot(x_all_list, S, N):
    x_all_array = np.array(x_all_list)
    x_mat = np.reshape(x_all_array, [N, S])
    x_i_mat = x_mat.T
    return x_i_mat

x_i_mat = reshepe_to_plot(x_all_list, S, N)
y_i_mat = reshepe_to_plot(y_all_list, S, N)
x_i_outmat = reshepe_to_plot(x_all_outlist, S, N)
y_i_outmat = reshepe_to_plot(y_all_outlist, S, N)

fig = plt.figure(figsize=(9,5))

ax1 = fig.add_subplot(1, 2, 1)
ax1.set_xlim(0, 1)
ax1.set_ylim(0, 1)
ax1.set_aspect('equal', adjustable='box')
ax1.set_xlabel('x')
ax1.set_ylabel('y')
titleN = ax1.text(0, 1.05, 'N = {}'.format(N), fontsize=15)

pi_renge = (max(pi_all_list) - min(pi_all_list))*0.1

ax2 = fig.add_subplot(1, 2, 2)
ax2.set_position([0.58,0.18,0.35,0.65])
ax2.set_xlim(0, S)
ax2.set_ylim(min(pi_all_list) - pi_renge, max(pi_all_list) + pi_renge)
ax2.set_xlabel('Number of trials')
ax2.set_ylabel('pi')

ims = []
xline = np.linspace(0, 1, 100)
yline = np.sqrt(1-xline**2)

for s in range(S):
    ax1 = fig.add_subplot(1, 2, 1)
    ax1.set_aspect('equal', adjustable='box')
    ax1.set_xlim(0, 1)
    ax1.set_ylim(0, 1)
    im1 = plt.plot(x_i_mat[s], y_i_mat[s], marker='o', c = 'blue' ,linestyle='None')
    title = ax1.text(0.5, 1.05, 'S = {}'.format(s+1), fontsize=15)
    im2 = plt.plot(x_i_outmat[s], y_i_outmat[s], marker='o', c = 'black' ,linestyle='None')
    plt.plot(xline, yline, color = 'red')
    
    ax2 = fig.add_subplot(1, 2, 2)
    ax2.set_position([0.58,0.18,0.35,0.65])
    ax2.set_xlim(0, S)
    ax2.set_ylim(min(pi_all_list) - pi_renge, max(pi_all_list) + pi_renge)
    im3 = plt.plot(pi_all_list[0:s], marker='o', c = 'blue' ,linestyle='None')
    im4 = plt.plot(pi_list[0:s], c = 'red')
    pia = ax1.text(1.2, 1.2, 'pi = {}'.format(pi_list[s]), fontsize=15)
    pistda = ax1.text(1.2, 1.1, 'pi_std = {}'.format(pi_std_list[s]), fontsize=15)
    
    im = im1 + im2 + im3 + im4 + [title] + [pia] + [pistda]
    
    ims.append(im)
    
    ax1.axis('off')
    ax2.axis('off')

ani = animation.ArtistAnimation(fig, ims, interval=100, repeat=False)
ani.save("caliculate_pi.gif", writer="imagemagick")

plt.show()

実行結果

総打点数$N=1000$、繰り返し回数$S=100$とした結果を上に示しました。左の図がそれぞれ試行結果で青い点数が$p$となります。右図は試行を繰り返した際の計算値$\pi$の値を示しています。青点が$S$回目の計算値でその平均値が赤線で示されています。各パラメーターを変化させてその計算結果が次の表です。

　これは打点数$N$を変化させた時の$\pi$と$\pi$の標準偏差と計算時間を示したものです。打点数を増加させると標準偏差が減少して計算時間が増加していく事が分かります。ここで、打点数が小さい時でも$\pi$真の値(3.14159...)に近くなることもありますが、打点数が小さい時には標準偏差が大きいのでそれは偶然真の値に近づいたに過ぎず、精度が低いと言えます。

　この表は試行回数$S$を変化させた時の$\pi$と$\pi$の標準偏差と計算時間です。試行回数が変化しても$\pi$の計算値、$\pi$の標準偏差が共に大きく変化しないという結果になりました。 $\pi$の近似値の平均は、アニメーションの右図の赤線の様に、試行回数が5回目くらいまではが大きく変化し、その後一定の値に収束していきます。したがって、標準偏差も数回目までは大きくみえます。実際は、各試行で偏差が変化しないため、試行数を多くしても標準偏差は変化せず、試行単体の標準偏差に収束していくはずです。

以上です。