﻿import numpy as np  # モジュールNumpyをnpという名前でインポート
import matplotlib.pyplot as plt  # 図の出力に使用するモジュール
n = 400  # 試行回数
np.random.seed(1)  # 乱数のシード値を固定
# P_(i+1)=(xi[i], eta[i])
# U(0,1)に従う乱数をn個生成し，xiに代入
xi = np.random.uniform(0, 1, n)  
# U(0,1)に従う乱数をn個生成し，etaに代入
eta = np.random.uniform(0, 1, n)
def ND(n):  # 自然数nに対してNn(D)を返す関数
    s = 0
    for i in range(n):
        if xi[i]**2 + eta[i]**2 <= 1:
            s += 1
    return s

Z = np.zeros(n)  # Nn(D)/nを格納するベクトル
for i in range(n):
    Z[i] = ND(i+1) / (i+1)

plt.subplot(1, 2, 1)  # 1行2列のうち1列目
plt.xlabel('n')
plt.ylabel('Nn(D)/n')
plt.plot(Z, label = 'Nn(D)/n')
# 理論値μの定数関数を点線で描画
plt.axhline(y = np.pi / 4, label = 'μ', linestyle = "dashed")
plt.legend()
plt.subplot(1, 2, 2)  # 1行2列のうち2列目
plt.xlabel('x')
plt.ylabel('y')
plt.scatter(xi, eta)  # Pのプロット
x_values = np.linspace(0, 1, 100)
y_values = np.sqrt(1 - x_values**2)
plt.plot(x_values, y_values)  # x^2 + y^2 = 1 のグラフ
plt.tight_layout()  # グラフの配置を調節
plt.show()