根据累积风险函数的值计算形状参数和尺度参数

Question

威布尔概率密度函数给出如下：

累积风险函数表示为：

现在，我有一个数据集，其中 t 和 H 是这样的：

## Import libraries
from scipy.optimize import curve_fit
from scipy import stats
from scipy.optimize import minimize
import numpy as np
import matplotlib.pyplot as plt

## Load the dataset
t = np.array([4,        6,         8,       10])        ## Time data         
H = np.array([0.266919, 0.518181,  0.671102, 1.808351]) ## Cumulative Hazard

我希望找到上述数据集的形状参数（ρ）和尺度参数（λ）。

为了计算简单，我使用了以下对数公式：

因此时间向量和累积风险的对数可以这样写：

ln_t = np.log(t)   # log of time vector
ln_H = np.log(H)   # log of Cumulative Hazard

我应用了以下两种方法

Method-1：通过使用线性函数的曲线拟合：

## Define a linear function
def ln_cum_hazard(LN_T, rho, myu):
    LN_H = rho * LN_T - myu
    return LN_H

## Model parameters
popt, pcov = curve_fit(ln_cum_hazard, ln_t, ln_H)

## Shape parameter
ρ = popt[0]
print("\n ρ (shape parameter) = \n", ρ)

## Scale parameter
λ = np.exp(popt[1]/popt[0])
print("\n λ (scale parameter) = \n", λ)

## Predicted values of H
H_pred = (t/λ)**ρ

## Plot H
plt.plot(t,H,'r',label = 'Actual Data')
plt.plot(t,H_pred,'b',label = 'Curve fit: scipy.optimize')
plt.legend()
plt.title('ρ = {} and λ = {}'.format(ρ,λ))
plt.xlabel('t')
plt.ylabel('H')
plt.show()

通过 scipy.optimize 进行曲线拟合，我得到以下拟合曲线：

方法2：通过stats.norm.logpdf()进行最大似然估计

## Define the linear function
def ln_cumulative_haz(params):
    rho = params[0]               # ρ (rho):    Shape parameter
    myu = params[1]               # λ (lambda): Scale parameter
    sd  = params[2]
      
    ## Linear function
    LN_H = rho * ln_t - myu

    ## Define the negative log likelihood function
    LL = -np.sum( stats.norm.logpdf(ln_H, loc = LN_H, scale=sd ) )

    return(LL)

## Inital parameters
initParams = [1, 1, 1]

## Minimize the negative log likelihood function
results = minimize(ln_cumulative_haz, initParams, method='Nelder-Mead')
  
## Model parameters
estParms = results.x

## Shape parameter
ρ = estParms[0]
print("\n ρ (shape parameter) = \n", ρ)

## Scale parameter
λ = np.exp(estParms[1]/estParms[0])
print("\n λ (scale parameter) = \n", λ)

## Predicted values of H
H_pred = (t/λ)**ρ

## Plot H
plt.plot(t,H,'r',label = 'Actual Data')
plt.plot(t,H_pred,'b',label = 'Curve fit: stats.norm.logpdf()')
plt.legend()
plt.title('ρ = {} and λ = {}'.format(ρ,λ))
plt.xlabel('t')
plt.ylabel('H')
plt.show()

通过 stats.norm.logpdf() 通过 MLE，我得到以下曲线拟合：

通过这两种方法，我对 ρ（形状参数）和 λ（尺度参数）的精度几乎相同。

现在我有以下3个疑惑：

1. 计算 ρ 和 λ 参数的程序是否适用于这两种方法？

2. 为了定义负对数似然函数，我使用了“stats.norm.logpdf”。然而，数据的底层函数是一个二参数威布尔分布。这是正确的吗？

3. 有没有其他方法可以在 Python 中计算 ρ 和 λ？

有人可以帮我解决这些疑问吗？

Edit-1：考虑 scipy.stats.weibull_min

我已应用以下代码来最小化 weibull 函数的负对数似然。

## Define the linear function
def ln_cumulative_haz_weibull(params):
    rho    = params[0]               # ρ (rho):    Shape parameter
    myu    = params[1]               # λ (lambda): Scale parameter
    shape  = params[2]
    sd     = params[3]
    
   

    ## Linear function
    LN_H = rho * ln_t - myu

    ## Calculate negative log likelihood
    LL = -np.sum( stats.weibull_min.logpdf(x = ln_H,
                                           c = shape,
                                           loc = LN_H,
                                           scale = sd ) )

    return(LL)


## Inital parameters
initParams = [1, 1, 1, 1]   ## params[i]

## Results of MLE
results = minimize(ln_cumulative_haz_weibull, initParams, method='Nelder-Mead')


## Model parameters
estParms = results.x

ρ = estParms[0]
print("\n ρ (shape parameter) = \n", ρ)

λ = np.exp(estParms[1]/estParms[0])
print("\n λ (scale parameter) = \n", λ)


##ln_H_pred = estParms[0] * ln_t - estParms[1]
##print("\n ln_H_pred = \n",ln_H_pred)


H_pred = (t/λ)**ρ
print("\n H_pred = ",H_pred)
## Plot ln(H)
plt.plot(t,H,'r',label = 'Actual Data')
plt.plot(t,H_pred,'b',label = 'Curve fit: stats.weibull_min.logpdf')
plt.legend()
plt.title('ρ = {} and λ = {}'.format(ρ,λ))
plt.xlabel('t')
plt.ylabel('H')
plt.show()

但是用stats.weibull_min.logpdf()拟合出来的曲线拟合效果不好，出现这样的：

我还收到一条错误消息：

警告（来自警告模块）：文件 "C:\Users\user\AppData\Local\Programs\Python\Python37\lib\site-packages\scipy\optimize\optimize.py", 第597行 numpy.max(numpy.abs(fsim[0] - fsim[1:])) <= fatol): RuntimeWarning: invalid value encountered in subtract

有人可以帮我解决我哪里错了吗？

Answer 1

我同意你的结果：

ρ（形状参数）= 1.914

λ（尺度参数）= 8.3568

装修相当不错：

我不能说你画了什么。检查它。