双参数对数似然最大化异常：ODE模型参数优化调试

以下是针对你的对数似然最大化问题的核心修正点及完整代码：

核心问题分析 #

1. 目标函数参数不匹配：scipy.optimize.minimize 仅接受单个数组作为参数输入，但你的 loglik 函数要求两个独立参数。
2. 最大化/最小化逻辑颠倒：你需要最大化对数似然，等价于最小化负对数似然，但原代码返回的是似然值的指数结果，方向完全错误。
3. 无真实观测数据：原代码用模型输出作为观测值的均值，形成循环依赖，优化失去拟合目标。
4. 初始参数格式错误：minimize 要求初始猜测为单个数组，而非两个独立值。

修正后完整代码 #

import numpy as np
from scipy.integrate import odeint
from scipy.optimize import minimize
import math

# 模型基础参数
N = 1
I0, R0 = 0.001, 0
U0 = N - I0 - R0
J0 = I0  # 初始累计发病率
Lf0, Ls0 = 0, 0
true_beta, true_gamma = 8, 0.4  # 真实参数（用于生成模拟数据）
mu, muTB, sigma_rec, rho = 1/80, 1/6, 1/6, 0.03
u, v, w = 0.88, 0.083, 0.0006
t = np.linspace(0, 500, 501)

# SIR扩展模型微分方程
def deriv(y, t, N, beta, gamma, mu, muTB, sigma_rec, rho, u, v, w):
    U, Lf, Ls, I, R, cInc = y
    b = (mu * (U + Lf + Ls + R)) + (muTB * I)
    lamda = beta * I
    clamda = 0.2 * lamda
    dU = b - ((lamda + mu) * U)
    dLf = (lamda*U) + ((clamda)*(Ls + R)) - ((u + v + mu) * Lf)
    dLs = (u * Lf) - ((w + clamda + mu) * Ls)
    dI = w*Ls + v*Lf - ((gamma + muTB + sigma_rec) * I) + (rho * R)
    dR = ((gamma + sigma_rec) * I) - ((rho + clamda + mu) * R)
    cI = w*Ls + v*Lf + (rho * R)  # 日发病率
    return dU, dLf, dLs, dI, dR, cI

# 生成模拟观测数据（用真实参数生成）
solve_true = odeint(deriv, (U0, Lf0, Ls0, I0, R0, J0), t, args=(N, true_beta, true_gamma, mu, muTB, sigma_rec, rho, u, v, w))
_, _, _, I_true, _, cInc_true = solve_true.T

true_prev = I_true[-1] * 100000  # 真实终末患病率（每10万人口）
true_inc = (cInc_true[1:] - cInc_true[:-1])[-1] * 100000  # 真实最后一日发病率（每10万人口）

# 添加对数正态噪声生成模拟观测值
sigmaPrev = 40  # 观测患病率标准差
sigmaInc = 30   # 观测发病率标准差

log_mu_prev = np.log(true_prev**2 / np.sqrt(true_prev**2 + sigmaPrev**2))
log_sigma_prev = np.sqrt(np.log(1 + (sigmaPrev**2 / true_prev**2)))
observed_prev = np.random.lognormal(log_mu_prev, log_sigma_prev)

log_mu_inc = np.log(true_inc**2 / np.sqrt(true_inc**2 + sigmaInc**2))
log_sigma_inc = np.sqrt(np.log(1 + (sigmaInc**2 / true_inc**2)))
observed_inc = np.random.lognormal(log_mu_inc, log_sigma_inc)

print(f"模拟观测患病率：{observed_prev:.2f} / 10万")
print(f"模拟观测发病率：{observed_inc:.2f} / 10万")

# 负对数似然函数（用于最小化）
def neg_loglik(params):
    beta, gamma = params
    # 用当前参数求解ODE
    solve = odeint(deriv, (U0, Lf0, Ls0, I0, R0, J0), t, args=(N, beta, gamma, mu, muTB, sigma_rec, rho, u, v, w))
    _, _, _, I, _, cInc = solve.T
    
    # 模型预测值（转换为每10万人口）
    model_prev = I[-1] * 100000
    model_inc = (cInc[1:] - cInc[:-1])[-1] * 100000
    
    # 患病率对数似然计算
    log_mu_prev = np.log(model_prev**2 / np.sqrt(model_prev**2 + sigmaPrev**2))
    log_sigma_prev = np.sqrt(np.log(1 + (sigmaPrev**2 / model_prev**2)))
    L_prev = -0.5 * ((np.log(observed_prev) - log_mu_prev)/log_sigma_prev)**2 - np.log(observed_prev * log_sigma_prev * np.sqrt(2*math.pi))
    
    # 发病率对数似然计算
    log_mu_inc = np.log(model_inc**2 / np.sqrt(model_inc**2 + sigmaInc**2))
    log_sigma_inc = np.sqrt(np.log(1 + (sigmaInc**2 / model_inc**2)))
    L_inc = -0.5 * ((np.log(observed_inc) - log_mu_inc)/log_sigma_inc)**2 - np.log(observed_inc * log_sigma_inc * np.sqrt(2*math.pi))
    
    # 返回负对数似然
    return -(L_prev + L_inc)

# 运行优化
initial_guess = [7, 0.3]  # 初始猜测（与真实值略有差异）
result = minimize(neg_loglik, x0=initial_guess, method='L-BFGS-B', bounds=[(0.1, 20), (0.01, 1)])

print("\n优化结果：")
print(result)
print(f"\n估计beta值：{result.x[0]:.4f}")
print(f"估计gamma值：{result.x[1]:.4f}")
print(f"真实beta值：{true_beta}，真实gamma值：{true_gamma}")

关键修正说明 #

1. 参数格式统一：将 loglik 改为接受单个数组参数，内部拆分为 beta 和 gamma。
2. 目标函数调整：返回负对数似然值，让 minimize 实现最大化对数似然的效果。
3. 添加观测数据：用真实参数生成带噪声的模拟观测数据，给优化提供明确的拟合目标。
4. 参数边界约束：通过 bounds 参数确保 beta 和 gamma 始终为正（符合流行病学意义）。

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