プロビット回帰モデル

def probit_regression_irls(
    X: np.ndarray,
    y: np.ndarray,
    max_iter: int = 25,
    tol: float = 1e-8,
) -> dict:
    """
    IRLSによるプロビット回帰の最尤推定

    Parameters
    ----------
    X : np.ndarray, shape (n, p)
        計画行列（切片列を含む）
    y : np.ndarray, shape (n,)
        目的変数（0/1）
    max_iter : int
        最大反復回数
    tol : float
        収束判定の閾値

    Returns
    -------
    dict
        beta, se, log_likelihood, n_iter
    """
    n, p = X.shape
    beta = np.zeros(p)

    for iteration in range(max_iter):
        eta = X @ beta
        pi = stats.norm.cdf(eta)
        phi = stats.norm.pdf(eta)

        # 数値安定性のためクリップ
        pi = np.clip(pi, 1e-10, 1 - 1e-10)
        phi = np.clip(phi, 1e-10, None)

        # 重み: w_i = phi(eta_i)^2 / [pi_i * (1 - pi_i)]
        w = phi**2 / (pi * (1 - pi))

        # 作業従属変数: z_i = eta_i + (y_i - pi_i) / phi(eta_i)
        z = eta + (y - pi) / phi

        # 重み付き最小二乗法
        W = np.diag(w)
        XtWX = X.T @ W @ X
        XtWz = X.T @ W @ z

        beta_new = np.linalg.solve(XtWX, XtWz)

        # 収束判定
        if np.max(np.abs(beta_new - beta)) < tol:
            beta = beta_new
            break
        beta = beta_new

    # 最終結果
    eta = X @ beta
    pi = stats.norm.cdf(eta)
    pi = np.clip(pi, 1e-15, 1 - 1e-15)
    log_lik = np.sum(y * np.log(pi) + (1 - y) * np.log(1 - pi))

    # 標準誤差
    phi = stats.norm.pdf(eta)
    w = phi**2 / (pi * (1 - pi))
    XtWX = X.T @ np.diag(w) @ X
    cov = np.linalg.inv(XtWX)
    se = np.sqrt(np.diag(cov))

    return {
        "beta": beta,
        "se": se,
        "log_likelihood": log_lik,
        "n_iter": iteration + 1,
        "cov": cov,
    }

import pandas as pd
import statsmodels.api as sm

# データ生成（プロビットモデルから）
rng = np.random.default_rng(42)
n = 500
x1 = rng.normal(0, 1, n)
x2 = rng.normal(0, 1, n)
beta_true = np.array([-0.3, 0.7, -0.5])  # beta_0, beta_1, beta_2
X = np.column_stack([np.ones(n), x1, x2])
eta = X @ beta_true
pi_true = stats.norm.cdf(eta)
y = rng.binomial(1, pi_true)

print(f"サンプルサイズ: {n}")
print(f"真の係数: \u03b2\u2080={beta_true[0]}, \u03b2\u2081={beta_true[1]}, \u03b2\u2082={beta_true[2]}")
print(f"y=1 の割合: {y.mean():.3f}")

# スクラッチ実装
result_irls = probit_regression_irls(X, y)

# statsmodels
model_sm = sm.GLM(y, X, family=sm.families.Binomial(link=sm.families.links.Probit()))
result_sm = model_sm.fit()

# 比較
comparison = pd.DataFrame({
    "真の値": beta_true,
    "IRLS (自前)": result_irls["beta"],
    "statsmodels": result_sm.params,
    "SE (自前)": result_irls["se"],
    "SE (statsmodels)": result_sm.bse,
}, index=["\u03b2\u2080", "\u03b2\u2081", "\u03b2\u2082"])

print(f"IRLS 反復回数: {result_irls['n_iter']}")
print(f"対数尤度 (自前): {result_irls['log_likelihood']:.4f}")
print(f"対数尤度 (statsmodels): {result_sm.llf:.4f}")
print()
comparison

beta_hat = result_sm.params

# 平均限界効果（AME）
eta_hat = X @ beta_hat
phi_hat = stats.norm.pdf(eta_hat)
ame = phi_hat.mean() * beta_hat[1:]  # 切片を除く

# 平均での限界効果（MEM）
x_mean = X.mean(axis=0)
eta_mean = x_mean @ beta_hat
phi_mean = stats.norm.pdf(eta_mean)
mem = phi_mean * beta_hat[1:]

me_df = pd.DataFrame({
    "\u03b2": beta_hat[1:],
    "AME": ame,
    "MEM": mem,
}, index=["x\u2081", "x\u2082"])

print("限界効果の比較")
print(f"  \u03c6(X\u03b2 の平均): {phi_hat.mean():.4f}")
print(f"  \u03c6(X\u0304\u03b2):       {phi_mean:.4f}")
print()
me_df

from statsmodels.formula.api import probit, logit

# statsmodels 組み込みの Titanic データ
titanic = sm.datasets.get_rdataset("Titanic", "datasets").data

# 集計データを個票データに展開
rows = []
for _, row in titanic.iterrows():
    for _ in range(int(row["Freq"])):
        rows.append({
            "Class": row["Class"],
            "Sex": row["Sex"],
            "Age": row["Age"],
            "Survived": 1 if row["Survived"] == "Yes" else 0,
        })
df = pd.DataFrame(rows)
print(f"サンプルサイズ: {len(df)}")
print(f"生存率: {df['Survived'].mean():.3f}")

# プロビット回帰
model_probit = probit("Survived ~ C(Class) + C(Sex) + C(Age)", data=df).fit(disp=0)
print(model_probit.summary())

# ロジスティック回帰との係数比較
model_logit = logit("Survived ~ C(Class) + C(Sex) + C(Age)", data=df).fit(disp=0)

coef_comparison = pd.DataFrame({
    "Probit \u03b2": model_probit.params,
    "Logit \u03b2": model_logit.params,
    "Probit/Logit": model_probit.params / model_logit.params,
})

print(f"Probit/Logit 比の平均: {(model_probit.params / model_logit.params).mean():.3f}")
print(f"理論的近似値:          0.600")
print()
coef_comparison

# 予測確率の比較
pred_probit = model_probit.predict(df)
pred_logit = model_logit.predict(df)

fig, axes = plt.subplots(1, 2, figsize=(8, 3), constrained_layout=True)

axes[0].scatter(pred_logit, pred_probit, alpha=0.3, s=10)
axes[0].plot([0, 1], [0, 1], "r--")
axes[0].set(xlabel="Logit predicted probability",
            ylabel="Probit predicted probability",
            title="Predicted probabilities")

axes[1].hist(pred_probit - pred_logit, bins=30, edgecolor="white")
axes[1].set(xlabel="Probit \u2212 Logit", ylabel="Frequency",
            title="Difference in predicted probabilities")

plt.show()

# 限界効果（AME）
# formula APIで推定したモデルから計画行列を取得
from statsmodels.discrete.discrete_model import Probit

margeff = model_probit.get_margeff(at="overall")
print(margeff.summary())