Python で二項分布・ポアソン分布をシミュレーションする

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

# 二項分布 B(n=10, p=0.5)
n, p = 10, 0.5
binom_dist = stats.binom(n=n, p=p)

# P(X = 5) を計算
print(f"P(X = 5) = {binom_dist.pmf(5):.4f}")

# P(X <= 3) を計算
print(f"P(X <= 3) = {binom_dist.cdf(3):.4f}")

# 期待値と分散
print(f"期待値: {binom_dist.mean():.2f}")
print(f"分散: {binom_dist.var():.2f}")

# 理論的な確率分布
x = np.arange(0, n + 1)
pmf = binom_dist.pmf(x)

# 乱数を 10000 回生成してシミュレーション
np.random.seed(42)
samples = binom_dist.rvs(size=10000)

# グラフの描画
plt.figure(figsize=(10, 6))
plt.bar(x - 0.2, pmf, width=0.4, label='理論値', alpha=0.7)
plt.hist(samples, bins=np.arange(-0.5, n + 1.5, 1), density=True, 
         width=0.4, align='mid', label='シミュレーション', alpha=0.7)
plt.xlabel('表の回数')
plt.ylabel('確率')
plt.title(f'二項分布 B({n}, {p})')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

# ポアソン分布 Po(λ=4)
lam = 4
poisson_dist = stats.poisson(mu=lam)

# P(X = 3) を計算
print(f"P(X = 3) = {poisson_dist.pmf(3):.4f}")

# P(X >= 6) を計算
print(f"P(X >= 6) = {1 - poisson_dist.cdf(5):.4f}")

# 期待値と分散（ポアソン分布では両方とも λ）
print(f"期待値: {poisson_dist.mean():.2f}")
print(f"分散: {poisson_dist.var():.2f}")

# 理論的な確率分布
x = np.arange(0, 15)
pmf = poisson_dist.pmf(x)

# 乱数を 10000 回生成
np.random.seed(42)
samples = poisson_dist.rvs(size=10000)

# グラフの描画
plt.figure(figsize=(10, 6))
plt.bar(x - 0.2, pmf, width=0.4, label='理論値', alpha=0.7)
plt.hist(samples, bins=np.arange(-0.5, 15.5, 1), density=True,
         width=0.4, align='mid', label='シミュレーション', alpha=0.7)
plt.xlabel('イベント数')
plt.ylabel('確率')
plt.title(f'ポアソン分布 Po({lam})')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

# B(n=1000, p=0.004) と Po(λ=4) を比較
n_large, p_small = 1000, 0.004
binom_large = stats.binom(n=n_large, p=p_small)
poisson_approx = stats.poisson(mu=n_large * p_small)

x = np.arange(0, 15)
plt.figure(figsize=(10, 6))
plt.bar(x - 0.2, binom_large.pmf(x), width=0.4, label=f'B({n_large}, {p_small})', alpha=0.7)
plt.bar(x + 0.2, poisson_approx.pmf(x), width=0.4, label=f'Po({n_large * p_small})', alpha=0.7)
plt.xlabel('x')
plt.ylabel('確率')
plt.title('二項分布のポアソン近似')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()