﻿import numpy as np  # モジュールNumpyをnpという名前でインポート
import matplotlib.pyplot as plt  # グラフ描画のためのモジュール
a = 0.5
# 整数 200, 400, 600, ..., 200000 を格納した配列
n_values = np.arange(200, 200200, 200)

# 大数の強法則：m1
def expXvars(n):
  np.random.seed(1)  # 乱数のシード値を固定
  # U(-a,a)に従う乱数をn個生成し，Xに代入
  X = np.random.uniform(-a, a, n)
  return sum(np.exp(X)) / n

v_expXvars = np.vectorize(expXvars)  # スカラー関数をベクトル化
m1 = (np.exp(a) - np.exp(-a)) / (2 * a)
plt.xlabel("n / 200")
plt.ylabel("expXvars(n)")
plt.plot(v_expXvars(n_values), label = 'expXvars', linewidth = 0.5)
plt.axhline(y = m1, label = 'm1', linestyle = "dashed")
plt.legend()
plt.show()

# 大数の強法則：m2
def XSQvars(n):
  np.random.seed(1)  # 乱数のシード値を固定
  # U(-a,a)に従う乱数をn個生成し，Xに代入
  X = np.random.uniform(-a, a, n)
  return sum(X * X) / n

v_XSQvars = np.vectorize(XSQvars)  # スカラー関数をベクトル化
m2 = a**2 / 3
plt.xlabel("n / 200")
plt.ylabel("XSQvars(n)")
plt.plot(v_XSQvars(n_values), label = 'XSQvars', linewidth = 0.5)
plt.axhline(y = m2, label = 'm2', linestyle = "dashed")
plt.legend()
plt.show()

# 大数の強法則：m3
def XSQvars_over_expXvars(n):
  np.random.seed(1)  # 乱数のシード値を固定
  # U(-a,a)に従う乱数をn個生成し，Xに代入
  X = np.random.uniform(-a, a, n)
  return sum(X * X) / sum(np.exp(X))

# スカラー関数をベクトル化
v_XSQvars_over_expXvars = np.vectorize(XSQvars_over_expXvars)
m3 = 2 * a**3 / (3 * (np.exp(a) - np.exp(-a)))
plt.xlabel("n / 200")
plt.ylabel("XSQvars_over_expXvars(n)")
plt.plot(v_XSQvars_over_expXvars(n_values), 
         label = 'XSQvars_over_expXvars', linewidth = 0.5)
plt.axhline(y = m3, label= 'm3', linestyle = "dashed")
plt.legend()
plt.show()