﻿import numpy as np  # モジュールNumpyをnpという名前でインポート
m = 100000  # 試行回数
n = 5  # 景品の種類
Yn = np.zeros(m)  # 景品を全種類集めるまでの購入回数を格納する配列
np.random.seed(1)  # 乱数のシード値を固定
for i in range(m):
    num_of_purchases = 0  # 購入回数
    coupon_array = np.zeros(n)  # n次元配列(初期値は(0,0,...,0))   
    while sum(coupon_array) != n:
        # 0から(n-1)までの整数から無作為に1つの整数を取得
        k = np.random.choice(n, size=1)
        # k番目の要素を1とし，(k+1)番目の景品を取得済みとみなす
        coupon_array[k] = 1
        num_of_purchases += 1
    Yn[i] = num_of_purchases

# 平均と分散の理論値を計算
t = np.sum(1 / np.arange(1, n+1))
mean_theoretical = n * t
u = np.sum(np.arange(1,n) / (n - np.arange(1, n)) ** 2)
var_theoretical = n * u
# 平均E(Yn)を表示
print(np.mean(Yn))  # モンテカルロ法
print(mean_theoretical)  # 参考：理論値
# 分散V(Yn)を表示
print(np.var(Yn))  # モンテカルロ法
print(var_theoretical)  # 参考：理論値