基于累积Hazard函数求解Weibull分布形状与尺度参数的疑问

背景与数据集 #

已知双参数Weibull分布的概率密度函数与累积风险函数，现有时间t和累积风险H的数据集：

## 导入库
from scipy.optimize import curve_fit
from scipy import stats
from scipy.optimize import minimize
import numpy as np
import matplotlib.pyplot as plt

## 加载数据集
t = np.array([4,        6,         8,       10])        ## 时间数据         
H = np.array([0.266919, 0.518181,  0.671102, 1.808351]) ## 累积风险数据

Weibull累积风险函数为( H(t) = (t/\lambda)^\rho )，通过对数转换可转化为线性形式：

ln_t = np.log(t)   # 时间的对数向量
ln_H = np.log(H)   # 累积风险的对数向量

参数求解方法 #

方法一：线性函数曲线拟合 #

利用对数转换后的线性关系，直接用线性拟合求解参数：

## 定义线性函数
def ln_cum_hazard(LN_T, rho, myu):
    LN_H = rho * LN_T - myu
    return LN_H

## 拟合模型参数
popt, pcov = curve_fit(ln_cum_hazard, ln_t, ln_H)

## 形状参数
ρ = popt[0]
print("\n ρ (形状参数) = \n", ρ)

## 尺度参数
λ = np.exp(popt[1]/popt[0])
print("\n λ (尺度参数) = \n", λ)

## 预测H值
H_pred = (t/λ)**ρ

## 绘图
plt.plot(t,H,'r',label = '实际数据')
plt.plot(t,H_pred,'b',label = '曲线拟合: scipy.optimize')
plt.legend()
plt.title('ρ = {} and λ = {}'.format(ρ,λ))
plt.xlabel('t')
plt.ylabel('H')
plt.show()

方法二：基于正态分布的极大似然估计（MLE） #

假设对数转换后的残差服从正态分布，构建负对数似然函数求解参数：

## 定义负对数似然函数
def ln_cumulative_haz(params):
    rho = params[0]               # ρ (形状参数)
    myu = params[1]               # 用于推导λ的中间参数
    sd  = params[2]               # 正态分布的标准差
      
    ## 线性模型预测值
    LN_H_pred = rho * ln_t - myu

    ## 计算负对数似然
    LL = -np.sum( stats.norm.logpdf(ln_H, loc = LN_H_pred, scale=sd ) )

    return(LL)

## 初始参数
initParams = [1, 1, 1]

## 最小化负对数似然
results = minimize(ln_cumulative_haz, initParams, method='Nelder-Mead')
  
## 提取估计参数
estParms = results.x

## 形状参数
ρ = estParms[0]
print("\n ρ (形状参数) = \n", ρ)

## 尺度参数
λ = np.exp(estParms[1]/estParms[0])
print("\n λ (尺度参数) = \n", λ)

## 预测H值
H_pred = (t/λ)**ρ

## 绘图
plt.plot(t,H,'r',label = '实际数据')
plt.plot(t,H_pred,'b',label = '曲线拟合: stats.norm.logpdf()')
plt.legend()
plt.title('ρ = {} and λ = {}'.format(ρ,λ))
plt.xlabel('t')
plt.ylabel('H')
plt.show()

疑问 #

1. 两种方法的参数求解流程是否正确？
2. 数据基于双参数Weibull分布，却用正态分布构建似然函数是否合理？
3. Python中还有其他求解ρ与λ的方法吗？

尝试Weibull似然拟合的错误问题 #

尝试用stats.weibull_min.logpdf构建负对数似然函数拟合，但拟合效果差，且出现警告：

Warning (from warnings module): File "C:\Users\user\AppData\Local\Programs\Python\Python37\lib\site-packages\scipy\optimize\optimize.py", line 597 numpy.max(numpy.abs(fsim[0] - fsim[1:])) <= fatol): RuntimeWarning: invalid value encountered in subtract

对应的错误代码：

## 定义负对数似然函数
def ln_cumulative_haz_weibull(params):
    rho    = params[0]               # ρ (形状参数)
    myu    = params[1]               # 用于推导λ的中间参数
    shape  = params[2]
    sd     = params[3]
    
    ## 线性模型预测值
    LN_H_pred = rho * ln_t - myu

    ## 计算负对数似然
    LL = -np.sum( stats.weibull_min.logpdf(x = ln_H,
                                           c = shape,
                                           loc = LN_H_pred,
                                           scale = sd ) )

    return(LL)


## 初始参数
initParams = [1, 1, 1, 1]

## MLE结果
results = minimize(ln_cumulative_haz_weibull, initParams, method='Nelder-Mead')


## 提取估计参数
estParms = results.x

ρ = estParms[0]
print("\n ρ (形状参数) = \n", ρ)

λ = np.exp(estParms[1]/estParms[0])
print("\n λ (尺度参数) = \n", λ)

## 预测H值
H_pred = (t/λ)**ρ
print("\n H_pred = ",H_pred)
## 绘图
plt.plot(t,H,'r',label = '实际数据')
plt.plot(t,H_pred,'b',label = '曲线拟合: stats.weibull_min.logpdf')
plt.legend()
plt.title('ρ = {} and λ = {}'.format(ρ,λ))
plt.xlabel('t')
plt.ylabel('H')
plt.show()

问题解答与错误排查 #

对三个疑问的解答 #

1. 方法流程正确性

   • 方法一：流程完全正确。Weibull累积风险函数对数转换后确实是线性关系，通过线性拟合得到参数后反推λ，逻辑符合分布的数学性质。
   • 方法二：流程逻辑成立，但本质是线性回归+正态误差假设，不是Weibull分布原生的极大似然估计，只是一种近似拟合手段。

2. 正态分布似然的合理性
   不合理。双参数Weibull分布的MLE应该基于自身的概率密度/累积分布函数构建似然，用正态分布假设相当于强行给对数残差套正态分布，和Weibull分布的统计性质不匹配，仅能作为近似拟合，不能得到严格意义上的Weibull参数MLE。

3. 其他求解方法

   • Scipy内置拟合函数：scipy.stats.weibull_min.fit()可直接拟合双参数Weibull，无需手动构建似然，返回的c对应ρ，scale对应λ。
   • 生存分析库lifelines：WeibullFitter类专门用于Weibull分布的生存分析拟合，支持处理删失数据，还能输出参数置信区间。
   • 手动构建原生Weibull似然：基于原始时间数据，用Weibull的概率密度函数构建负对数似然，直接拟合ρ和λ。

基于stats.weibull_min.logpdf拟合的错误排查 #

你的代码存在三个核心问题：

1. 似然逻辑错误：将模型预测的LN_H_pred作为Weibull分布的loc参数完全不符合分布应用逻辑。正确思路应该是假设对数残差（ln_H - LN_H_pred）服从Weibull分布，需将残差作为x传入，且loc设为0（残差均值假设为0）。
2. 多余参数导致过拟合：额外引入的shape和sd参数属于冗余参数，不仅增加优化难度，还容易导致数值不稳定（如参数过小引发对数似然为无穷大，出现NaN）。
3. 优化方法选择不当：Nelder-Mead方法对初始参数敏感，且无参数边界限制，容易出现无效参数值。建议使用L-BFGS-B方法，并给参数设置合理边界（如ρ>0，λ>0）。

修正后的示例代码：

# 基于Weibull残差的负对数似然
def ln_cumulative_haz_weibull(params):
    rho = params[0]
    myu = params[1]
    # 模型预测的对数累积风险
    LN_H_pred = rho * ln_t - myu
    # 计算残差（取绝对值，因为Weibull适用于非负数据）
    residuals = np.abs(ln_H - LN_H_pred)
    # 假设残差服从形状参数为1的Weibull分布（可根据需求调整）
    LL = -np.sum(stats.weibull_min.logpdf(x=residuals, c=1, loc=0, scale=np.std(residuals)))
    return LL

# 初始参数
initParams = [1, 1]
# 带参数边界的优化
results = minimize(ln_cumulative_haz_weibull, initParams, method='L-BFGS-B', bounds=((1e-5, None), (-10, 10)))

rho = results.x[0]
lambda_ = np.exp(results.x[1]/rho)
print("ρ (形状参数) =", rho)
print("λ (尺度参数) =", lambda_)

# 预测与绘图
H_pred = (t/lambda_)**rho
plt.plot(t,H,'r',label='实际数据')
plt.plot(t,H_pred,'b',label='曲线拟合: Weibull残差MLE')
plt.legend()
plt.title(f'ρ = {rho:.4f} and λ = {lambda_:.4f}')
plt.xlabel('t')
plt.ylabel('H')
plt.show()

更简单的原生Weibull拟合方式：

# 用scipy内置函数直接拟合双参数Weibull（固定位置参数为0）
params = stats.weibull_min.fit(t, floc=0)
rho_fit = params[0]
lambda_fit = params[2]
print("ρ (形状参数) =", rho_fit)
print("λ (尺度参数) =", lambda_fit)

# 预测与绘图
H_pred_fit = (t/lambda_fit)**rho_fit
plt.plot(t,H,'r',label='实际数据')
plt.plot(t,H_pred_fit,'g',label='曲线拟合: scipy.weibull_min.fit')
plt.legend()
plt.title(f'ρ = {rho_fit:.4f} and λ = {lambda_fit:.4f}')
plt.xlabel('t')
plt.ylabel('H')
plt.show()

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