統計学 – Python で学ぶ 最尤推定法

from itertools import product

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

np.random.seed(0)


def likelihood(mu, sigma, x):
    # μ=mu, σ=sigma の正規分布の結合確率密度関数の値を計算する
    return norm.pdf(x, loc=mu, scale=sigma).prod()


# 母集団分布が μ=0, σ=4 の正規分布であるとき、大きさ10のランダム標本の観測値を得る。
x = norm.rvs(0, 4, size=10)
print(x)

fig = plt.figure(figsize=(9, 9))
fig.subplots_adjust(wspace=0.3, hspace=0.3)
for i, (mu, sigma) in enumerate(product([-5, 0, 5], [3, 4, 5])):
    ax = fig.add_subplot(3, 3, i + 1)

    # μ=mu, σ=sigma の正規分布の確率密度関数を描画する
    xs = np.linspace(-15, 15, 100)
    ys = norm.pdf(xs, loc=mu, scale=sigma)
    ax.plot(xs, ys, "g")

    # ランダム標本の観測値を描画する。
    y = norm.pdf(x, loc=mu, scale=sigma)
    markerline, stemline, baseline = ax.stem(
        x, y, linefmt="k:", markerfmt="ko", basefmt="none", use_line_collection=True
    )
    plt.setp(markerline, ms=3)

    l = likelihood(mu, sigma, x)
    ax.set_title(fr"$L({mu}, {sigma};\mathbf{{x}})$ = {l:.2e}")
    ax.set_xticks([-10, 0, 10])
    ax.set_yticks([0.0, 0.1, 0.2])
    ax.grid()
    
    if i % 3 == 0:
        ax.set_ylabel(fr"$\mu = {mu}$", fontsize=12)
    if i // 3 == 2:
        ax.set_xlabel(fr"$\sigma = {sigma}$", fontsize=12)

plt.show()

fig, ax = plt.subplots()

x = np.linspace(0.01, 1, 100)
y = np.log(x)
ax.plot(x, y)

ax.grid()
ax.set_xlabel(r"$L(\mathbf{\theta};\mathbf{x})$")
ax.set_ylabel(r"$\log L(\mathbf{\theta};\mathbf{x})$")

plt.show()

import numpy as np
from scipy.stats import bernoulli

np.random.seed(0)

# 母集団分布が p=0.3 の正ベルヌーイ分布であるとき、大きさ50のランダム標本の観測値を得る。
p = 0.3
x = bernoulli.rvs(p, size=50)
print(x)

# 最尤推定値を計算する。
p_pred = x.mean()
print(p_pred)

np.random.seed(0)

for n in np.geomspace(1, 10000000, 8, dtype=int):
    # 母集団分布が p=0.3 の正ベルヌーイ分布であるとき、大きさnのランダム標本の観測値を得る。
    x = bernoulli.rvs(0.3, size=n)

    # 最尤推定値を計算する。
    p_pred = x.mean()
    print(f"N={n:<8} p={p_pred:.6f}")

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

np.random.seed(0)


def likehood(x, mu, sigma):
    return np.prod(norm.pdf(x, loc=mu, scale=sigma))


# 集団分布が平均0、標準偏差40の正規分布であるとき、大きさ100のランダム標本の観測値を得る。
mu, sigma = 0, 40
x = norm.rvs(loc=mu, scale=sigma, size=100)  # size を増やしすぎると、float で表せなくなります。

# 各点の尤度関数の値を計算する。
X, Y = np.mgrid[-10:10:30j, 20:60:30j]
Z = []
for m, t in zip(X.ravel(), Y.ravel()):
    Z.append(likehood(x, mu=m, sigma=t))
Z = np.array(Z).reshape(*X.shape)

# 最尤推定量を計算する。
mu_pred, sigma_pred = x.mean(), x.std()

fig = plt.figure(figsize=(12, 6))

# 3D グラフを描画する。
ax1 = fig.add_subplot(121, projection="3d")

ax1.plot_wireframe(X, Y, Z)

ax1.set_xlabel(r"$\mu$", fontsize=15)
ax1.set_ylabel(r"$\sigma$", fontsize=15)

# 等高線を描画する。
ax2 = fig.add_subplot(122)

ax2.contour(X, Y, Z)
ax2.scatter(mu, sigma, c="g", label=fr"True: $\mu$={mu:.2f}, $\sigma$={sigma:.2f}")
ax2.scatter(
    mu_pred,
    sigma_pred,
    c="r",
    label=fr"Estimation: $\mu$={mu_pred:.2f}, $\sigma$={sigma_pred:.2f}",
)

ax2.set_xlabel(r"$\mu$", fontsize=15)
ax2.set_ylabel(r"$\sigma$", fontsize=15)
ax2.grid()
ax2.legend()

plt.show()