﻿import numpy as np  # モジュールNumpyをnpという名前でインポート
import matplotlib.pyplot as plt  # 図の出力に使用するモジュール
def f(x):  # 被積分関数f(x)
    return np.sin(x)

# モンテカルロ法による積分計算
def monte_carlo(n, a):
    sigma = 0
    np.random.seed(1)  # 乱数のシード値を固定
    for i in range(n):
        Z = np.random.uniform(0,a)  # U(0,a)に従う乱数を1個生成
        sigma = sigma + f(Z)
    return sigma / n

# プロット関数
def mc_plot(n, a):
    y = np.zeros(n)
    for i in range(n):
        # ベクトルyにmonte(i+1, a) (i=0,1,...,n-1)を順に代入
        y[i] = monte_carlo(i+1, a)
    plt.title('n=%d' % n)
    plt.xlabel('n', fontname = 'MS Gothic')
    plt.ylabel('Monte Carlo integration', fontname = 'MS Gothic')
    plt.plot(y)  # (i,y[i])をプロットし，直線で結ぶ

a = np.pi  # 円周率π
mc_plot(1000, a)
integral_val = 2  # ∫[0,a]f(x)dx=∫[0,π]sin(x)dxの理論値
# モンテカルロ積分の系列が理論値に収束する様子を描画
plt.axhline(y = integral_val / a, color = 'blue', linestyle = 'dashed') 
plt.show()