﻿import numpy as np  # モジュールNumpyをnpという名前でインポート
from scipy import stats  # 高度な統計分析を行うためのサブモジュール
p0 = 0.65  # pの真のパラメータ値：p0
np.random.seed(1)  # 乱数のシード値を固定
# 幾何分布Ge(p0)に従う乱数を10000個生成 (※-1を補正)
x = np.random.geometric(p = p0, size = 10000) - 1
p_hat = 1 / (1 + np.mean(x))  # 最尤推定値を計算
print("p_hat =", p_hat)  # 最尤推定値を表示