﻿from scipy import stats  # 高度な統計分析を行うためのサブモジュール
import numpy as np  # モジュールNumpyをnpという名前でインポート
import math  # 数学関係の機能をまとめたモジュールをインポート
import matplotlib.pyplot as plt  # 図の出力に使用するモジュールをインポート
# 負の二項分布の確率関数
def probf(r, p, k):
    val = math.factorial(r+k-1)
    val = val / (math.factorial(r-1) * math.factorial(k))
    val = val * (p**r) * ((1-p)**k)
    return val

# 負の二項分布のパラメータ
r = 10
p = 0.3
n = 100000  # 試行回数
np.random.seed(1)  # 乱数のシード値を固定
Yr = stats.nbinom.rvs(r, p, size = n)  # 負の二項分布に従う乱数
# 期待値E(Yr)を表示
print(np.mean(Yr))  # モンテカルロ法
print(r * (1 - p) / p)  # 参考：理論値
# 分散V(Yr)を表示
print(np.var(Yr))  # モンテカルロ法
print(r * (1 - p) / (p**2))  # 参考：理論値
# Yrの乱数の最大値以下の非負整数を格納した配列
x = np.arange(0, max(Yr)+1, 1)
P = [0] * (x.size)  # xの各値に対する確率関数の値
for i in range(x.size):
    P[i] = probf(r, p, x[i])

plt.title("Yr")
plt.hist(Yr, bins = 15, density = True)  # Yrのヒストグラム
plt.plot(x, P)  # Yrの確率関数
plt.show()  # ヒストグラムと確率関数の比較