如何用极大似然估计反推非线性模型参数(α,ξ,β1,β2)

你拥有由以下非线性模型生成的数据集(T,H)，生成逻辑与代码如下：

import numpy as np
from matplotlib import pyplot as plt

## 定义模型曲线
def curve(t,α,ξ,β1,β2):  
    ## 模型参数：(α, ξ, β1, β2)
    ## 协变量1（指数分布）
    X1 = np.exp(t/20) 
    ## 协变量2（指数分布）
    X2 = np.exp(t/10) 
    ## 双参数Weibull基准危险率
    ho = (α/ξ) * (t/ξ)**(α-1) 
    ## 最终危险率
    h = ho * np.exp(β1*X1 + β2*X2)
    return h

## 时间轴
T = np.linspace(0, 5, 1000)

## 生成数据集
H = curve(T, α = 2.3 ,ξ = 4.8e-4, β1 = 1.1, β2 = 2.1)

## 可视化
plt.rcParams["figure.figsize"] = [7.50, 7.50]
plt.rcParams["figure.autolayout"] = True
plt.plot(T, H, color='red')
plt.show()

数据集的可视化结果如下：

核心思路：极大似然估计的实现步骤 #

你的模型是带协变量的Weibull危险率模型，形式为：
$h(t; \alpha, \xi, \beta_1, \beta_2) = \frac{\alpha}{\xi} \left( \frac{t}{\xi} \right)^{\alpha-1} \exp\left( \beta_1 e^{t/20} + \beta_2 e^{t/10} \right)$

要通过MLE反推参数，需遵循以下步骤：

1. 定义似然函数
   假设观测数据带有独立同分布的正态噪声：$H_i = h(t_i; \theta) + \epsilon_i$，其中$\epsilon_i \sim N(0, \sigma^2)$，$\theta = (\alpha, \xi, \beta_1, \beta_2)$。对应的对数似然函数为：
   $\mathcal{L}(\theta, \sigma) = -\frac{N}{2}\ln(2\pi\sigma^2) - \frac{1}{2\sigma^2}\sum_{i=1}N \left( H_i - h(t_i; \theta) \right)^2$
   最大化该似然等价于最小化残差平方和（可忽略常数项），同时需满足参数约束：$\alpha>0, \xi>0$（Weibull参数的物理意义要求）。

2. 数值优化求解
   使用支持边界约束的数值优化算法（如L-BFGS-B），最小化负对数似然函数（优化器通常默认找最小值）。

3. 验证拟合结果
   将估计出的参数代入原模型，对比拟合曲线与观测数据的匹配度。

完整实现代码 #

import numpy as np
from scipy.optimize import minimize
from matplotlib import pyplot as plt

# 原模型函数
def curve(t, α, ξ, β1, β2):
    X1 = np.exp(t/20)
    X2 = np.exp(t/10)
    ho = (α/ξ) * (t/ξ)**(α-1)
    h = ho * np.exp(β1*X1 + β2*X2)
    return h

# 生成带噪声的模拟数据（模拟真实观测场景）
np.random.seed(42)
T = np.linspace(0, 5, 1000)
true_params = [2.3, 4.8e-4, 1.1, 2.1]
H_true = curve(T, *true_params)
# 添加10%均值水平的正态噪声
H = H_true + np.random.normal(0, 0.1*np.mean(H_true), size=len(T))

# 定义负对数似然函数（含参数合法性约束）
def neg_log_likelihood(params, t, h_obs):
    α, ξ, β1, β2 = params
    # 排除无效参数（α和ξ必须大于0）
    if α <= 0 or ξ <= 0:
        return np.inf
    h_pred = curve(t, α, ξ, β1, β2)
    # 正态噪声下的负对数似然（忽略常数项）
    return np.sum((h_obs - h_pred)**2)

# 设置参数初始值（尽量接近真实值或合理范围）
init_params = [2.0, 1e-4, 1.0, 2.0]

# 运行优化：L-BFGS-B支持边界约束
result = minimize(
    neg_log_likelihood,
    init_params,
    args=(T, H),
    method='L-BFGS-B',
    # 参数边界：α>0, ξ>0，β1/β2设合理范围
    bounds=[(1e-5, None), (1e-6, None), (-5, 5), (-5, 5)]
)

# 输出估计结果
print("参数估计结果：")
print(f"α: {result.x[0]:.4f} (真实值: {true_params[0]})")
print(f"ξ: {result.x[1]:.6f} (真实值: {true_params[1]})")
print(f"β1: {result.x[2]:.4f} (真实值: {true_params[2]})")
print(f"β2: {result.x[3]:.4f} (真实值: {true_params[3]})")

# 可视化拟合效果
plt.figure(figsize=(7.5,7.5))
plt.plot(T, H, 'r.', label='观测数据', alpha=0.3)
plt.plot(T, curve(T, *result.x), 'b-', label='MLE拟合曲线')
plt.plot(T, H_true, 'g--', label='真实曲线')
plt.legend()
plt.show()

关键注意事项 #

• 参数约束：必须保证$\alpha>0$和$\xi>0$，否则Weibull基准危险率无意义，优化时通过边界约束或函数内判断实现。
• 初始值选择：合理的初始值能大幅提升优化效率和准确性，建议基于领域知识或数据特征设置。
• 噪声假设：上述代码假设噪声为正态分布，若实际噪声分布不同，需对应修改似然函数的形式。

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