import pandas as pd
import numpy as np
import pymc as pm
import arviz as az
import matplotlib.pyplot as plt

# サンプルコード（データの読み込み部分は適宜修正してください）
# data_csv = pd.read_csv('your_data.csv')

# 仮のデータを作成（実際のデータを使用してください）
np.random.seed(42)
data_csv = pd.DataFrame({
    'well': np.repeat(['A', 'B', 'C'], 10),
    'DIV': np.tile(np.arange(1, 11), 3),
    'Number of Spikes': np.random.poisson(lam=20, size=30)
})

# サンプルコードの続き
data_csv["sample_code"] = data_csv["well"].astype("category").cat.codes
sample_idx = data_csv["sample_code"].values
n_samples = data_csv["sample_code"].nunique()
print(f"Number of samples: {n_samples}")

# GLMMを構築（ランダム切片を削除）
with pm.Model() as model:
    # 固定切片
    a = pm.Normal("a", mu=0, sigma=10)

    # ランダムスロープのハイパーパラメータ
    mu_b = pm.Normal("mu_b", mu=0, sigma=10)
    sigma_b = pm.HalfNormal("sigma_b", sigma=10)

    # ランダムスロープ
    b = pm.Normal("b", mu=mu_b, sigma=sigma_b, shape=n_samples)

    # 線形予測子
    log_mu = a + b[sample_idx] * data_csv["DIV"].values
    mu = pm.math.exp(log_mu)

    # 負の二項分布のパラメータ
    alpha = pm.Exponential("alpha", 1)
    Y_obs = pm.NegativeBinomial('Y_obs', mu=mu, alpha=alpha, observed=data_csv["Number of Spikes"].values)

    # サンプリング
    trace = pm.sample(2000, tune=2000, target_accept=0.9, chains=4, cores=4, random_seed=42)

# 予測データの作成
unique_DIV = np.unique(data_csv["DIV"])
pred_df = pd.DataFrame({
    "well": np.repeat(data_csv["well"].unique(), len(unique_DIV)),
    "DIV": np.tile(unique_DIV, n_samples)
})
pred_df["sample_code"] = pred_df["well"].astype("category").cat.codes
pred_sample_idx = pred_df["sample_code"].values

# 事後分布からパラメータを抽出
posterior = trace.posterior

# 'a' と 'b' のサンプルを取得
a_samples = posterior['a'].stack(sample=("chain", "draw")).values  # shape: (total_draws,)
b_samples = posterior['b'].stack(sample=("chain", "draw")).values  # shape: (n_samples, total_draws)

# 総サンプル数
n_posterior_samples = a_samples.shape[0]

# 予測データポイントの数
n_predictions = pred_df.shape[0]

# 予測スパイク数を格納する配列を初期化
mu_pred = np.zeros((n_posterior_samples, n_predictions))

for i in range(n_posterior_samples):
    # 固定効果の切片
    a_i = a_samples[i]  # スカラー値
    # 各サンプルのランダムスロープ
    b_i = b_samples[:, i]  # shape: (n_samples,)

    # 対数リンク関数を使用した線形予測子の計算
    log_mu_i = a_i + b_i[pred_sample_idx] * pred_df["DIV"].values

    # 期待値の計算
    mu_pred[i, :] = np.exp(log_mu_i)

# 平均値と95%信用区間の計算
mu_mean = mu_pred.mean(axis=0)
mu_hpd = az.hdi(mu_pred, hdi_prob=0.99)  # 信用区間を99%に設定

# 予測結果をデータフレームに追加
pred_df['mu_mean'] = mu_mean
pred_df['mu_lower'] = mu_hpd[:, 0]
pred_df['mu_upper'] = mu_hpd[:, 1]

# プロット
for sample in pred_df["well"].unique():
    sample_df = pred_df[pred_df["well"] == sample]
    plt.figure(figsize=(10, 6))
    plt.plot(sample_df["DIV"], sample_df["mu_mean"], color="black", label="Predicted")
    plt.fill_between(sample_df["DIV"], sample_df["mu_lower"], sample_df["mu_upper"],
                     color="gray", alpha=0.5, label="99% CI")
    plt.scatter(data_csv[data_csv["well"] == sample]["DIV"], data_csv[data_csv["well"] == sample]["Number of Spikes"],
                color="red", label="Observed")
    plt.title(f"Sample: {sample}")
    plt.xlabel("DIV")
    plt.ylabel("Number of Spikes")
    plt.legend()
    plt.show()


上記のコードで、40日とかの日にちの99%信頼区間は十分に裾幅があるけど、手前の日にちの裾幅が小さい。
それを改良したい。