第2章：確率分布

イントロダクション

第1章では、データの特徴を要約する記述統計と、不確実性を扱う確率論の基礎（確率の公理、条件付き確率、ベイズの定理、期待値と分散）を学びました。この章では、それらの知識を土台に、確率変数がとりうる値とその確率の対応関係を表す確率分布（Probability Distribution）について詳しく見ていきます。

機械学習では、観測データがどのような確率分布から生成されたと仮定するかによって、使用するモデルや評価方法が変わります。例えば、コインの表裏のような二値データにはベルヌーイ分布、ある期間内に発生するイベントの回数にはポアソン分布、多くの自然現象や測定誤差には正規分布が適しています。適切な確率分布を選べるようになることは、データ分析やモデル設計における重要な第一歩です。

1. 離散型確率分布

離散型確率分布（Discrete Probability Distribution）は、確率変数がとりうる値が飛び飛びの値（0, 1, 2, ...など数えられる値）である場合の分布です。各値に対応する確率を表す関数を確率質量関数（Probability Mass Function, PMF）と呼びます。

1.1 ベルヌーイ分布

ベルヌーイ分布（Bernoulli Distribution）は、成功（1）か失敗（0）の2つの結果しかとらない、最も単純な確率分布です。コイン投げやクリック有無の判定など、二値の結果をモデル化する際に使われます。

成功確率を$p$とすると、PMFは次のように表せます：

$$P(X=k) = p^k (1-p)^{1-k}, \quad k \in \{0, 1\}$$

期待値と分散：

$$E[X] = p, \qquad \text{Var}(X) = p(1-p)$$

1.2 二項分布

二項分布（Binomial Distribution）は、成功確率$p$のベルヌーイ試行を$n$回独立に繰り返したときの、成功回数の分布です。「10回コインを投げて表が出る回数」のような場面で使います。

$$P(X=k) = \binom{n}{k} p^k (1-p)^{n-k}, \quad k = 0, 1, \ldots, n$$

ここで$\binom{n}{k} = \frac{n!}{k!(n-k)!}$は二項係数（$n$個から$k$個を選ぶ組み合わせの数）です。

期待値と分散：

$$E[X] = np, \qquad \text{Var}(X) = np(1-p)$$

1.3 ポアソン分布

ポアソン分布（Poisson Distribution）は、単位時間・単位面積あたりに発生する稀な事象の回数をモデル化する分布です。コールセンターへの着信数、Webサイトへのアクセス数、工場での欠陥品数など、「ある期間内に何回起こるか」を数える場面で広く使われます。

平均発生率を$\lambda$（ラムダ）とすると：

$$P(X=k) = \frac{\lambda^k e^{-\lambda}}{k!}, \quad k = 0, 1, 2, \ldots$$

期待値と分散：

$$E[X] = \lambda, \qquad \text{Var}(X) = \lambda$$

ポアソン分布の特徴的な性質は、期待値と分散が同じ値$\lambda$になることです。また、二項分布において試行回数$n$を大きく、成功確率$p$を小さくしながら$np=\lambda$を一定に保つと、二項分布はポアソン分布に近づくことが知られています。

1.4 Pythonでの実装

SciPyのscipy.statsモジュールを使うと、各分布のPMF、期待値、分散を簡単に計算できます。

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

# --- ベルヌーイ分布 Bernoulli(p) ---
p = 0.3
bernoulli = stats.bernoulli(p)
print("=== ベルヌーイ分布 Bernoulli(p=0.3) ===")
print(f"P(X=0) = {bernoulli.pmf(0):.4f}")
print(f"P(X=1) = {bernoulli.pmf(1):.4f}")
print(f"期待値 E[X] = {bernoulli.mean():.4f}")
print(f"分散 Var(X) = {bernoulli.var():.4f}")

# --- 二項分布 Binomial(n, p) ---
n, p = 10, 0.3
binomial = stats.binom(n, p)
print("\n=== 二項分布 Binomial(n=10, p=0.3) ===")
for k in [0, 3, 5, 10]:
    print(f"P(X={k}) = {binomial.pmf(k):.4f}")
print(f"期待値 E[X] = {binomial.mean():.4f}")
print(f"分散 Var(X) = {binomial.var():.4f}")

# --- ポアソン分布 Poisson(lambda) ---
lam = 3
poisson = stats.poisson(lam)
print("\n=== ポアソン分布 Poisson(λ=3) ===")
for k in [0, 2, 3, 6]:
    print(f"P(X={k}) = {poisson.pmf(k):.4f}")
print(f"期待値 E[X] = {poisson.mean():.4f}")
print(f"分散 Var(X) = {poisson.var():.4f}")

# --- 可視化 ---
fig, axes = plt.subplots(1, 3, figsize=(15, 4))

k_bern = [0, 1]
axes[0].bar(k_bern, bernoulli.pmf(k_bern), color='skyblue', edgecolor='black')
axes[0].set_title('ベルヌーイ分布 (p=0.3)')
axes[0].set_xlabel('k')
axes[0].set_ylabel('確率')
axes[0].set_xticks(k_bern)

k_binom = np.arange(0, n + 1)
axes[1].bar(k_binom, binomial.pmf(k_binom), color='lightgreen', edgecolor='black')
axes[1].set_title('二項分布 (n=10, p=0.3)')
axes[1].set_xlabel('k')
axes[1].set_ylabel('確率')

k_poisson = np.arange(0, 12)
axes[2].bar(k_poisson, poisson.pmf(k_poisson), color='salmon', edgecolor='black')
axes[2].set_title('ポアソン分布 (λ=3)')
axes[2].set_xlabel('k')
axes[2].set_ylabel('確率')

plt.tight_layout()
plt.show()

実行結果:

=== ベルヌーイ分布 Bernoulli(p=0.3) ===
P(X=0) = 0.7000
P(X=1) = 0.3000
期待値 E[X] = 0.3000
分散 Var(X) = 0.2100

=== 二項分布 Binomial(n=10, p=0.3) ===
P(X=0) = 0.0282
P(X=3) = 0.2668
P(X=5) = 0.1029
P(X=10) = 0.0000
期待値 E[X] = 3.0000
分散 Var(X) = 2.1000

=== ポアソン分布 Poisson(λ=3) ===
P(X=0) = 0.0498
P(X=2) = 0.2240
P(X=3) = 0.2240
P(X=6) = 0.0504
期待値 E[X] = 3.0000
分散 Var(X) = 3.0000

2. 連続型確率分布

連続型確率分布（Continuous Probability Distribution）は、確率変数が連続的な値（実数）をとる場合の分布です。個々の値に対する確率は常に0になるため、代わりに確率密度関数（Probability Density Function, PDF）$f(x)$を用いて、ある区間に値が入る確率を積分で表します。

$$P(a \leq X \leq b) = \int_{a}^{b} f(x)\, dx$$

2.1 正規分布

正規分布（Normal Distribution）はガウス分布（Gaussian Distribution）とも呼ばれ、身長・測定誤差・試験の得点など、多くの自然現象に現れる釣鐘型（ベル型）の分布です。平均$\mu$と分散$\sigma^2$の2つのパラメータで特徴づけられます。

$$f(x) = \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left(-\frac{(x-\mu)^2}{2\sigma^2}\right)$$

期待値と分散：

$$E[X] = \mu, \qquad \text{Var}(X) = \sigma^2$$

2.2 指数分布

指数分布（Exponential Distribution）は、ある事象が発生してから次に発生するまでの「待ち時間」をモデル化する分布です。機械の故障間隔や、顧客が到着するまでの時間などに使われます。ポアソン分布で表される事象間の待ち時間は、指数分布に従うという関係があります。

発生率$\lambda \gt 0$に対して：

$$f(x) = \lambda e^{-\lambda x}, \quad x \geq 0$$

期待値と分散：

$$E[X] = \frac{1}{\lambda}, \qquad \text{Var}(X) = \frac{1}{\lambda^2}$$

指数分布は無記憶性（Memorylessness）という特徴的な性質を持ちます。「すでに$s$時間待った」という条件は、その後さらに$t$時間待つ確率に影響しません：

$$P(X \gt s+t \mid X \gt s) = P(X \gt t)$$

2.3 ガンマ分布

ガンマ分布（Gamma Distribution）は指数分布を一般化した分布で、「独立な指数分布に従うイベントが$k$回起こるまでの待ち時間」を表します。形状パラメータ$k$（shape）と尺度パラメータ$\theta$（scale）を持ちます。

$$f(x) = \frac{1}{\Gamma(k)\theta^{k}} x^{k-1} e^{-x/\theta}, \quad x \geq 0$$

ここで$\Gamma(k)$はガンマ関数です。期待値と分散：

$$E[X] = k\theta, \qquad \text{Var}(X) = k\theta^2$$

2.4 Pythonでの実装

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

# --- 正規分布 Normal(mu, sigma) ---
mu, sigma = 0, 1
normal = stats.norm(mu, sigma)

# --- 指数分布 Exponential(lambda) ---
rate = 1.0  # レート λ（scipyでは scale = 1/λ で指定する）
exponential = stats.expon(scale=1 / rate)

# --- ガンマ分布 Gamma(k, theta) ---
shape, scale = 2.0, 2.0  # 形状パラメータk、尺度パラメータθ
gamma_dist = stats.gamma(a=shape, scale=scale)

print("=== 各分布のPDF値と要約統計量 ===")
print(f"正規分布 N(0,1):           f(0) = {normal.pdf(0):.4f}, "
      f"平均 = {normal.mean():.4f}, 分散 = {normal.var():.4f}")
print(f"指数分布 Exp(λ=1):         f(1) = {exponential.pdf(1):.4f}, "
      f"平均 = {exponential.mean():.4f}, 分散 = {exponential.var():.4f}")
print(f"ガンマ分布 Gamma(k=2,θ=2): f(2) = {gamma_dist.pdf(2):.4f}, "
      f"平均 = {gamma_dist.mean():.4f}, 分散 = {gamma_dist.var():.4f}")

fig, axes = plt.subplots(1, 3, figsize=(15, 4))

x_norm = np.linspace(-4, 4, 300)
axes[0].plot(x_norm, normal.pdf(x_norm), color='steelblue', linewidth=2)
axes[0].fill_between(x_norm, normal.pdf(x_norm), alpha=0.2, color='steelblue')
axes[0].set_title('正規分布 N(0, 1)')
axes[0].set_xlabel('x')
axes[0].set_ylabel('確率密度 f(x)')

x_exp = np.linspace(0, 6, 300)
axes[1].plot(x_exp, exponential.pdf(x_exp), color='darkorange', linewidth=2)
axes[1].fill_between(x_exp, exponential.pdf(x_exp), alpha=0.2, color='darkorange')
axes[1].set_title('指数分布 Exp(λ=1)')
axes[1].set_xlabel('x')
axes[1].set_ylabel('確率密度 f(x)')

x_gamma = np.linspace(0, 16, 300)
axes[2].plot(x_gamma, gamma_dist.pdf(x_gamma), color='seagreen', linewidth=2)
axes[2].fill_between(x_gamma, gamma_dist.pdf(x_gamma), alpha=0.2, color='seagreen')
axes[2].set_title('ガンマ分布 Gamma(k=2, θ=2)')
axes[2].set_xlabel('x')
axes[2].set_ylabel('確率密度 f(x)')

plt.tight_layout()
plt.show()

実行結果:

=== 各分布のPDF値と要約統計量 ===
正規分布 N(0,1):           f(0) = 0.3989, 平均 = 0.0000, 分散 = 1.0000
指数分布 Exp(λ=1):         f(1) = 0.3679, 平均 = 1.0000, 分散 = 1.0000
ガンマ分布 Gamma(k=2,θ=2): f(2) = 0.1839, 平均 = 4.0000, 分散 = 8.0000

3. 正規分布の性質と中心極限定理

3.1 標準正規分布と68-95-99.7ルール

平均$\mu=0$、分散$\sigma^2=1$の正規分布を標準正規分布（Standard Normal Distribution）と呼びます。任意の正規分布に従う確率変数$X$は、次の標準化（Standardization）によって標準正規分布に変換できます。

$$z = \frac{x - \mu}{\sigma}$$

この$z$の値はzスコア（z-score）と呼ばれ、「平均から標準偏差何個分離れているか」を表します。

正規分布には、平均から標準偏差の何倍以内にデータが収まるかを示す経験則があります。

import numpy as np
from scipy import stats

mu, sigma = 170, 8  # 例: 成人男性の身長（平均170cm、標準偏差8cm）
normal = stats.norm(mu, sigma)

print("=== 68-95-99.7ルールの検証 ===")
for k in [1, 2, 3]:
    lower, upper = mu - k * sigma, mu + k * sigma
    prob = normal.cdf(upper) - normal.cdf(lower)
    print(f"μ ± {k}σ の範囲 [{lower}, {upper}] に入る確率: {prob:.4f}")

# 標準化（zスコア）の計算
x = 186  # ある人の身長
z = (x - mu) / sigma
print(f"\n身長{x}cmの標準化スコア z = {z:.4f}")
print(f"標準正規分布での累積確率（この身長以下である確率）: {stats.norm.cdf(z):.4f}")

実行結果:

=== 68-95-99.7ルールの検証 ===
μ ± 1σ の範囲 [162, 178] に入る確率: 0.6827
μ ± 2σ の範囲 [154, 186] に入る確率: 0.9545
μ ± 3σ の範囲 [146, 194] に入る確率: 0.9973

身長186cmの標準化スコア z = 2.0000
標準正規分布での累積確率（この身長以下である確率）: 0.9772

3.2 中心極限定理

中心極限定理（Central Limit Theorem, CLT）は、統計学において最も重要な定理の一つです。次のように主張します。

$$\bar{X}_n = \frac{1}{n}\sum_{i=1}^{n} X_i \quad \xrightarrow{\ n \to \infty\ } \quad N\left(\mu, \frac{\sigma^2}{n}\right)$$

標準化した形では次のようになります：

$$\frac{\bar{X}_n - \mu}{\sigma/\sqrt{n}} \xrightarrow{d} N(0, 1)$$

中心極限定理が重要な理由は、元の母集団分布がサイコロの出目のような一様分布であっても、二項分布のように偏った分布であっても、標本平均を十分な数集めれば正規分布に近づくという点です。これにより、多くの統計的推測手法（信頼区間の構築や仮説検定など）で正規分布を前提とした計算が可能になります。

3.3 Pythonによるシミュレーション

サイコロを振る実験で中心極限定理を確認してみましょう。1回のサイコロの出目は1から6の一様分布に従いますが、複数回振った平均は正規分布に近づいていくはずです。

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(42)

def simulate_sample_means(sample_size, n_experiments=5000):
    """sample_size個のサイコロを振る実験をn_experiments回繰り返し、
    それぞれの標本平均を返す"""
    rolls = np.random.randint(1, 7, size=(n_experiments, sample_size))
    return rolls.mean(axis=1)

# サイコロ1個の理論的な期待値と分散
dice_mean = np.mean(np.arange(1, 7))          # 3.5
dice_var = np.var(np.arange(1, 7))            # 35/12 ≈ 2.9167

sample_sizes = [1, 2, 5, 30]
fig, axes = plt.subplots(1, 4, figsize=(18, 4))

for ax, n in zip(axes, sample_sizes):
    means = simulate_sample_means(n)
    ax.hist(means, bins=30, density=True, color='mediumpurple',
            edgecolor='black', alpha=0.7)
    ax.set_title(f'n = {n}')
    ax.set_xlabel('標本平均')
    ax.set_ylabel('密度')
    print(f"n={n}: 理論分散(σ²/n)={dice_var/n:.4f}, "
          f"シミュレーション分散={means.var():.4f}")

plt.suptitle('中心極限定理：サイコロの標本平均の分布', fontsize=14)
plt.tight_layout()
plt.show()

このシミュレーションでは、標本サイズ$n=1$のときはサイコロの出目そのものの一様な分布（1〜6の各値がほぼ等確率）になりますが、$n$が大きくなるにつれてヒストグラムは徐々に釣鐘型に近づき、$n=30$程度になるとほぼ正規分布と見分けがつかなくなります。また、標本平均の分散は理論値である$\sigma^2/n$（サイコロ1個の分散$35/12$を$n$で割った値）に近づいていくことも確認できます。

4. 確率分布のパラメータ推定

4.1 最尤推定の考え方

実際のデータ分析では、観測データがどの確率分布から生成されたかは分かっていても、その分布のパラメータ（正規分布の$\mu, \sigma$やポアソン分布の$\lambda$など）は未知であることがほとんどです。観測データから最も「もっともらしい」パラメータを求める代表的な手法が最尤推定（Maximum Likelihood Estimation, MLE）です。

観測データ$x_1, x_2, \ldots, x_n$が独立に同じ分布から得られたと仮定すると、パラメータ$\theta$のもとでこのデータが観測される同時確率（尤度関数（Likelihood Function））は次のように書けます。

$$L(\theta) = \prod_{i=1}^{n} f(x_i; \theta)$$

計算を簡単にするため、対数をとった対数尤度（Log-Likelihood）を最大化することが一般的です。

$$\ell(\theta) = \sum_{i=1}^{n} \ln f(x_i; \theta)$$

最尤推定量$\hat{\theta}$は、この対数尤度を最大にするパラメータとして定義されます。

$$\hat{\theta} = \underset{\theta}{\operatorname{argmax}}\ \ell(\theta)$$

4.2 正規分布のパラメータ推定

正規分布の場合、対数尤度を最大化すると、直感的にも自然な次の推定量が得られることが数学的に導けます。

$$\hat{\mu} = \bar{x} = \frac{1}{n}\sum_{i=1}^{n} x_i, \qquad \hat{\sigma}^2 = \frac{1}{n}\sum_{i=1}^{n} (x_i - \bar{x})^2$$

SciPyではscipy.stats.<分布名>.fit()メソッドを使うと、データから自動的に最尤推定を実行できます。

import numpy as np
from scipy import stats

np.random.seed(0)

# 真のパラメータ（通常は未知だが、検証のためあえて設定する）
true_mu, true_sigma = 50, 5
data = np.random.normal(true_mu, true_sigma, size=200)

# 手動で最尤推定量を計算
mle_mu = np.mean(data)
mle_sigma2 = np.mean((data - mle_mu) ** 2)  # MLEではnで割る
mle_sigma = np.sqrt(mle_sigma2)

# scipy.stats.norm.fit()による推定（内部で最尤推定を実行）
fit_mu, fit_sigma = stats.norm.fit(data)

print("=== 正規分布のパラメータ推定（最尤推定） ===")
print(f"真のパラメータ:        μ = {true_mu}, σ = {true_sigma}")
print(f"手動MLE推定値:         μ_hat = {mle_mu:.4f}, σ_hat = {mle_sigma:.4f}")
print(f"scipy.stats.norm.fit: μ_hat = {fit_mu:.4f}, σ_hat = {fit_sigma:.4f}")

手動計算による推定値とscipy.stats.norm.fit()の結果は一致します（乱数のシード次第で、真のパラメータとの間には標本サイズに応じたばらつきが生じます）。標本サイズ$n$を増やすほど、推定値は真のパラメータに近づいていきます。

4.3 二項分布のパラメータ推定

ベルヌーイ分布・二項分布の成功確率$p$についても、最尤推定量は観測された成功割合と一致することが導けます。

$$\hat{p} = \frac{\text{成功回数}}{\text{試行回数}}$$

import numpy as np
from scipy import stats

# 実験: 200枚のコインを投げて140枚が表だった
n_trials = 200
n_success = 140

# 最尤推定: p_hat = 成功回数 / 試行回数
p_mle = n_success / n_trials

print("=== 二項分布のパラメータ推定（最尤推定） ===")
print(f"試行回数: {n_trials}, 成功回数: {n_success}")
print(f"最尤推定値 p_hat = {p_mle:.4f}")

# 推定したpで二項分布を構築し、分布の要約統計量を確認
binom_est = stats.binom(n_trials, p_mle)
print(f"推定分布での期待値: {binom_est.mean():.4f}")
print(f"推定分布での標準偏差: {binom_est.std():.4f}")

実行結果:

=== 二項分布のパラメータ推定（最尤推定） ===
試行回数: 200, 成功回数: 140
最尤推定値 p_hat = 0.7000
推定分布での期待値: 140.0000
推定分布での標準偏差: 6.4807

5. 確率分布の可視化とシミュレーション

5.1 ヒストグラムと理論分布の重ね合わせ

観測データが特定の確率分布に従っているかを確認する基本的な方法は、データのヒストグラムに、推定したパラメータを持つ理論分布のPDFを重ねて描画することです。

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

np.random.seed(1)

# 例: 反応時間のデータ（ガンマ分布に従うと仮定）
data = np.random.gamma(shape=3.0, scale=1.5, size=500)

# ガンマ分布のパラメータをデータから推定（最尤推定）
# floc=0 で位置パラメータを0に固定し、形状・尺度パラメータのみ推定する
fit_shape, fit_loc, fit_scale = stats.gamma.fit(data, floc=0)

print("=== ガンマ分布のパラメータ推定 ===")
print(f"推定された形状パラメータ k = {fit_shape:.4f}")
print(f"推定された尺度パラメータ θ = {fit_scale:.4f}")

# ヒストグラムと推定分布の重ね合わせ
x = np.linspace(0, data.max(), 300)
fitted_pdf = stats.gamma.pdf(x, a=fit_shape, scale=fit_scale)

plt.figure(figsize=(10, 6))
plt.hist(data, bins=30, density=True, alpha=0.6, color='lightsteelblue',
         edgecolor='black', label='観測データ')
plt.plot(x, fitted_pdf, color='crimson', linewidth=2, label='推定されたガンマ分布')
plt.xlabel('値')
plt.ylabel('密度')
plt.title('ヒストグラムと推定分布の重ね合わせ')
plt.legend()
plt.grid(alpha=0.3)
plt.show()

推定された形状パラメータ$k$と尺度パラメータ$\theta$は、データ生成に使った真の値（$k=3.0$、$\theta=1.5$）に近い値となり、フィットした曲線がヒストグラムの形状によく一致することを確認できます。

5.2 Q-Qプロットによる分布の適合確認

Q-Qプロット（Quantile-Quantile Plot）は、観測データの分位点と、仮定した理論分布の分位点を散布図として比較する手法です。データが理論分布に従っていれば、点はほぼ直線上に並びます。

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

np.random.seed(2)

# 正規分布データと非正規分布（指数分布）データを比較
normal_data = np.random.normal(0, 1, 300)
exp_data = np.random.exponential(scale=1.0, size=300)

fig, axes = plt.subplots(1, 2, figsize=(12, 5))

stats.probplot(normal_data, dist="norm", plot=axes[0])
axes[0].set_title('Q-Qプロット：正規分布データ')

stats.probplot(exp_data, dist="norm", plot=axes[1])
axes[1].set_title('Q-Qプロット：指数分布データ（正規分布との比較）')

plt.tight_layout()
plt.show()

6. まとめと次のステップ

この章では、機械学習で頻繁に登場する確率分布と、それらをデータから特徴づける方法について学びました。

次のステップ

次章では、この章で学んだ確率分布の知識をもとに、標本から母集団の性質を推測する統計的推定・仮説検定について学びます。信頼区間の構築やp値の解釈など、データに基づいた意思決定に欠かせない手法を扱う予定です。

練習問題

from scipy import stats

lam = 4
poisson = stats.poisson(lam)

# 2件以下（0件、1件、2件）の確率の合計。累積分布関数(CDF)を使うと効率的
prob = poisson.cdf(2)
print(f"1時間に電話が2件以下である確率: {prob:.4f}")

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(42)

population_mean = 5.0       # Uniform(0,10)の理論平均 (a+b)/2
population_var = 100 / 12   # Uniform(0,10)の理論分散 (b-a)^2/12

n_experiments = 1000
sample_size = 50

sample_means = np.array([
    np.random.uniform(0, 10, sample_size).mean()
    for _ in range(n_experiments)
])

theoretical_mean = population_mean
theoretical_var = population_var / sample_size

print(f"理論的な標本平均の平均: {theoretical_mean:.4f}")
print(f"理論的な標本平均の分散: {theoretical_var:.4f}")
print(f"シミュレーションによる標本平均の平均: {sample_means.mean():.4f}")
print(f"シミュレーションによる標本平均の分散: {sample_means.var():.4f}")

plt.hist(sample_means, bins=30, density=True, color='teal', alpha=0.7, edgecolor='black')
plt.xlabel('標本平均')
plt.ylabel('密度')
plt.title('一様分布からの標本平均の分布（中心極限定理）')
plt.show()

import numpy as np

data = np.array([2, 4, 3, 5, 2, 6, 3, 4, 1, 5])

# ポアソン分布のパラメータλの最尤推定量は標本平均と一致する
lambda_mle = np.mean(data)

print(f"最尤推定量 λ_hat = {lambda_mle:.4f}")