モンテカルロ法でガウス関数の積分を求める

を求める。
とすると、乱数r(-a≤r≤a)に対する関数の値の値の平均値を計算し、積分区間の幅をこのにかけることで、もとの定積分の近似値が求まる

import random
import math
import numpy as np
import matplotlib.pyplot as plt
from tqdm import tqdm

a_s = [2,10,100,1000,10000]
Ks = [1000,5000,10000,20000,40000,100000]
areass = []
for a in tqdm(a_s):
    areas = []
    for K in Ks:
        vals = np.array([])
        for _ in range(K):
            r = random.random()*a
            val = 1/math.sqrt(2*math.pi)*math.exp(-0.5*r**2)
            vals = np.append(vals,val)

        average = vals.mean()
        width = 2*a
        area = average*width
        # print(area)
        areas.append(area)
    areass.append(areas)

for areas in areass:
    plt.plot(Ks,areas)

plt.xlabel('K')
plt.ylabel('area')
plt.legend(a_s)
plt.show()

Kの値を増加させることで定積分の値は収束する。
積分区間の幅aが大きいほど、収束は悪くなる。

区間幅aの値が1000と10000の場合にフォーカスしてみる。

import random
import math
import numpy as np
import matplotlib.pyplot as plt
from tqdm import tqdm

# a_s = [2,10,100,1000,10000]
a_s = [10,1000,10000]
# Ks = [1000,5000,10000,20000,40000,100000]
Ks = np.linspace(0,20000,100,dtype=int)
areass = []
for a in tqdm(a_s):
    areas = []
    for K in Ks:
        vals = np.array([])
        for _ in range(K):
            r = random.random()*a
            val = 1/math.sqrt(2*math.pi)*math.exp(-0.5*r**2)
            vals = np.append(vals,val)

        average = vals.mean()
        width = 2*a
        area = average*width
        # print(area)
        areas.append(area)
    areass.append(areas)

for areas in areass:
    plt.plot(Ks,areas)

plt.xlabel('K')
plt.ylabel('area')
plt.legend(a_s)
plt.show()

グラフがカクカクしているのは、以下の理由による。

以下が被積分関数f(x)のグラフだが、ほとんど|x|>2でゼロに漸近している。したがって、1000や10000という大きく幅をとっても、遠方では関数の値はほとんどゼロであり、積分の値に寄与しない。たまたま|x|<2の有益な値が乱数として出力されたときに値が入るので、急激な値のジャンプとして観測される（カクカクする）。