ベイズ更新とピカチュウ

【1】目的（What to Estimate）

くさむらでの「ピカチュウ出現確率」 θ を推定する。
出現を1、非出現を0とした観測データを D = {x₁, x₂, …, xₙ} とする。

xᵢ ∈ {0, 1},  θ = P(ピカチュウ出現)

---

【2】確率モデル（Bernoulli–Beta Model）

各観測はベルヌーイ分布に従う：

p(xᵢ | θ) = θ^(xᵢ) * (1 - θ)^(1 - xᵢ)

観測全体の尤度：

p(D | θ) = θ^k * (1 - θ)^(n - k)

（k = 出現回数）

---

【3】事前分布（Prior）

事前に「ピカチュウは出にくい」と思っている場合、
ベータ分布 Beta(α, β) を仮定。

p(θ) = (1 / B(α, β)) * θ^(α - 1) * (1 - θ)^(β - 1)

---

【4】ベイズ更新（Posterior）

観測 D 後の事後分布は：

p(θ | D) = Beta(α + k, β + n - k)

更新後の期待値：

E[θ | D] = (α + k) / (α + β + n)

---

【5】くさむらでの逐次観測（Online Bayesian Updating）

観測が1回ごとに入るたび：

αₜ = αₜ₋₁ + xₜ
βₜ = βₜ₋₁ + (1 - xₜ)
E[θₜ] = αₜ / (αₜ + βₜ)

→ 「ピカチュウが何回出たか」に応じて確率が連続的に修正される。

---

【6】時間連続モデル（Poisson–Gamma 拡張）

もし「出現イベント」を時間で数えるなら：

観測モデル：

N(t) ~ Poisson(λt)

事前分布：

λ ~ Gamma(α, β)

更新後の事後分布：

λ | D ~ Gamma(α + N(t), β + t)
E[λ | D] = (α + N(t)) / (β + t)

→ 「くさむら1分あたりのピカチュウ出現率」をリアルタイムで推定できる。

---

【7】要約（ピカチュウ × ベイズ推定）

概念	数式	くさむらでの意味
θ	出現確率	ピカチュウが現れる確率
α, β	信念の強さ	過去の経験（出た／出なかった）
k, n	観測データ	実際の遭遇回数と試行回数
E[θ	D]	期待出現率	次の探索で出る確率の推定
λ	出現頻度	単位時間あたりのピカチュウ発見率

---

したがって：

ピカチュウのくさむら出現 = 経験を通じて確率を更新するベイズ学習過程

観測を重ねるごとに、
「この草むらにはピカチュウがどれくらい出やすいか」を数学的に学んでいく。

# ===========================================================
# Pikachu Bayesian Inference in the Grass Field
# ===========================================================
# ピカチュウのくさむら出現確率 θ をベイズ推定で推定する
# Bernoulli–Betaモデル + 逐次更新 + Poisson–Gamma拡張
# ===========================================================

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import beta, gamma

# --- パラメータ初期化 ---
alpha, beta_param = 2, 8   # 事前: ピカチュウはやや出にくい
n_obs = 50                 # 試行回数（くさむら探索回数）
true_theta = 0.25          # 真の出現確率（仮定）

# --- 観測データ生成（出現=1, 非出現=0）---
data = np.random.binomial(1, true_theta, n_obs)

# --- ベイズ逐次更新 ---
alpha_t, beta_t = alpha, beta_param
posterior_means = []

for x_t in data:
    alpha_t += x_t
    beta_t += (1 - x_t)
    posterior_means.append(alpha_t / (alpha_t + beta_t))

# --- 結果 ---
print(f"Final posterior mean (E[θ|D]): {posterior_means[-1]:.3f}")
print(f"Total Pikachu sightings: {np.sum(data)} / {n_obs}")

# --- プロット1: 出現確率の逐次更新 ---
plt.figure(figsize=(8, 4))
plt.plot(posterior_means, color='gold', label="Posterior mean of θ")
plt.axhline(true_theta, color='red', linestyle='--', label='True θ (Ground Truth)')
plt.xlabel("Encounter count (n)")
plt.ylabel("Estimated P(Pikachu appears)")
plt.title("Online Bayesian Update of Pikachu Appearance Probability")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

# --- プロット2: 事前分布と最終事後分布の比較 ---
theta = np.linspace(0, 1, 200)
prior_pdf = beta.pdf(theta, alpha, beta_param)
posterior_pdf = beta.pdf(theta, alpha_t, beta_t)

plt.figure(figsize=(8, 4))
plt.plot(theta, prior_pdf, label=f"Prior Beta({alpha}, {beta_param})", color="gray")
plt.plot(theta, posterior_pdf, label=f"Posterior Beta({alpha_t:.0f}, {beta_t:.0f})", color="gold")
plt.title("Pikachu Appearance Probability Distribution (Prior → Posterior)")
plt.xlabel("θ = P(Pikachu appears)")
plt.ylabel("Probability Density")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

# ===========================================================
# 【Poisson–Gammaモデル】連続時間観測（くさむら出現頻度）
# ===========================================================

# 1分あたり出現率 λ の推定：N(t) ~ Poisson(λt), λ ~ Gamma(α, β)
alpha_g, beta_g = 2, 3  # 事前
t_values = np.linspace(1, 10, 10)  # 10分観測
observed_counts = np.random.poisson(lam=2.0 * t_values / 10)  # 仮想観測

posterior_alpha = alpha_g + np.sum(observed_counts)
posterior_beta = beta_g + t_values[-1]
lambda_est = posterior_alpha / posterior_beta

print(f"\nEstimated Pikachu appearance rate λ ≈ {lambda_est:.3f} per minute")

# --- プロット3: 事前・事後のGamma分布 ---
λ = np.linspace(0, 3, 300)
prior_pdf = gamma.pdf(λ, a=alpha_g, scale=1/beta_g)
posterior_pdf = gamma.pdf(λ, a=posterior_alpha, scale=1/posterior_beta)

plt.figure(figsize=(8, 4))
plt.plot(λ, prior_pdf, color='gray', label=f"Prior Gamma({alpha_g},{beta_g})")
plt.plot(λ, posterior_pdf, color='orange', label=f"Posterior Gamma({posterior_alpha:.0f},{posterior_beta:.0f})")
plt.title("Continuous Bayesian Update of Pikachu Appearance Rate λ")
plt.xlabel("λ = Expected sightings per minute")
plt.ylabel("Probability Density")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()