Python手动实现GLM泊松回归发散及多框架结果对齐求助

1. 统一输入基准 #

先确认所有实现使用完全相同的输入数据，包括是否添加截距项、特征是否预处理。基准结果为：截距1.8893，特征系数0.6698。

2. 修正sklearn的输出差异 #

sklearn的PoissonRegressor默认会对特征做中心化处理，这是导致结果不符的核心原因。修正方法：

• 关闭特征标准化：设置normalize=False
• 手动添加截距项（或依赖fit_intercept，但需确保输入特征未被预处理）

示例代码：

from sklearn.linear_model import PoissonRegressor
import numpy as np

# 假设X是无截距的特征矩阵，y是目标变量
X_with_intercept = np.hstack([np.ones((X.shape[0], 1)), X])
# 关闭标准化，强制使用原始特征拟合
poisson_sklearn = PoissonRegressor(fit_intercept=False, normalize=False, max_iter=1000)
poisson_sklearn.fit(X_with_intercept, y)
# 此时coef_输出应与statsmodels一致：[1.8893, 0.6698]

3. 修复手动Newton-Raphson的发散问题 #

你怀疑的W矩阵计算错误是核心问题，泊松回归的Newton-Raphson迭代需严格遵循以下规则：

• 预测均值mu = exp(X @ beta)，需避免exp下溢导致mu为0
• 权重矩阵W是对角矩阵，对角元素为mu（泊松回归方差等于均值，权重取mu）
• 迭代更新公式需正确计算梯度与海森矩阵

修正后的手动实现代码：

def newton_raphson_poisson(X, y, max_iter=100, tol=1e-8):
    # X需包含截距项
    beta = np.zeros(X.shape[1])
    for _ in range(max_iter):
        mu = np.exp(X @ beta)
        # 防止mu过小引发数值计算错误
        mu = np.maximum(mu, 1e-10)
        # 构建权重对角矩阵
        W = np.diag(mu)
        # 计算梯度与海森矩阵
        grad = X.T @ (y - mu)
        hessian = -X.T @ W @ X
        # 更新beta：beta = beta - 海森逆矩阵 @ 梯度
        delta = np.linalg.inv(hessian) @ grad
        beta -= delta
        # 检查收敛条件
        if np.linalg.norm(delta) < tol:
            break
    return beta

# 使用带截距的特征矩阵运行
beta_manual = newton_raphson_poisson(X_with_intercept, y)
# 输出应接近[1.8893, 0.6698]

关键修正点： #

• 修正W矩阵的对角元素为mu（而非其他错误值）
• 添加mu = np.maximum(mu, 1e-10)避免数值下溢
• 正确处理海森矩阵的符号，确保beta更新方向正确

4. 验证三者结果一致 #

确保：

• 所有实现使用带截距项的特征矩阵
• sklearn关闭标准化选项
• 手动实现的W矩阵、迭代公式正确

运行后三者的系数误差应控制在数值精度范围内（如1e-4以内），与基准结果对齐。

内容的提问来源于stack exchange，提问作者JeeyCi