3

我想用最大似然估计 (MLE) 估计 Sarhan 和 Apaloo (2013) 引入的指数修正威布尔扩展 (EMWE) 分布的四个参数。此分布用于可靠性和生存分析，以分析具有浴缸危险函数的数据集。

因为四个参数的对数似然函数的一阶导数给出了隐式解决方案，所以我尝试继续使用 Newton-Raphson 迭代方法。只是我对确定这种方法的良好初始猜测并不熟悉，因为我主要关注的是如何选择正确的初始猜测来获得浴缸危害函数。我正在使用 R，这是代码： 

#Parameter Estimation
xi<-c(3,4,4,4,4,4,4,4,4,4,4,4,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,8,8,8,8,8,8,8,8,9,9,9,9,9,10,10,11,12,12,14)
n<-length(xi)
log_likelihood<-function(theta){
  n*(log(theta[1])+log(theta[2])+log(theta[3])+(1-theta[2])*log(theta[4])+theta[1]*theta[4])+(theta[2]-1)*sum(log(xi))+(1/(theta[4]^theta[2]))*sum(xi^theta[2])-(theta[1]*theta[4])*sum(exp((xi/theta[4])^theta[2]))+(theta[3]-1)*sum(1-(exp((theta[1]*theta[4])*(1-(exp((xi/theta[4])^theta[2]))))))
}

result<-nlm(log_likelihood,theta<-c(0.000030763,2.2011,0.11878,6),hessian=TRUE)

我意识到在任何优化问题中，都没有系统地选择初始猜测。然而，在最大似然问题中，一个常见的选择是选择与观测数据的经验矩（均值、方差等）相匹配的初始参数，因此我选择：

theta[1]=0.0000022011
theta[2]=1.1878
theta[3]=0.21669
theta[4]=6.4941

其中theta[1]和theta[4]是比例参数，theta[2]作为theta[3]形状参数，结果如下：

result
____________________    
$minimum
[1] -1167.797

$estimate
[1] 1.540561e-06 1.187800e+00 2.166900e-01 6.494100e+00

$gradient
[1]  4.251830e+07  7.312058e+01  3.922662e+02 -1.815774e+01

$hessian
              [,1]          [,2]          [,3]          [,4]
[1,] -2.977361e+10 -7.455449e+02  1.093871e+03  3.396184e+02
[2,] -7.455449e+02 -5.215852e+01  6.316999e-04 -2.725559e+01
[3,]  1.093871e+03  6.316999e-04 -1.808591e+03 -2.906025e-04
[4,]  3.396184e+02 -2.725559e+01 -2.906025e-04  5.665718e+00

$code
[1] 2

$iterations
[1] 1

估计的参数是：

theta[1]=1.540561e-06
theta[2]=1.1878
theta[3]=2.166900e-01
theta[4]=6.4941

为了获得风险函数图，我使用以下代码：

#The estimated parameters
p1=theta[4]
p2=theta[2]
p3=theta[3]
p4=theta[1]

#Hazard Function
hazard_EMWE=(p4*p2*p3*((xi/p1)^(p2-1))*(exp(((xi/p1)^(p2))+(p4*p1*(1-(exp((xi/p1)^(p2)))))))*((1-(exp((p4*p1*(1-(exp((xi/p1)^(p2))))))))^(p3-1)))/(1-(((1-(exp((p4*p1*(1-(exp((xi/p1)^(p2))))))))^(p3))))

#Cumulative Distribution Function
cum_EMWE=((1-(exp((p4*p1*(1-(exp((xi/p1)^(p2))))))))^(p3))

#Hazard Function Plot
plot(xi,hazard_EMWE,type="l" ,main="Hazard Function Plot",xlab="x",ylab="Hazard Function",col = "black",lty=1,lwd=3,ylim=c(0.005,0.008))

这是情节： 浴缸危险情节

实际上，该图展示了浴缸风险曲线，但我认为我选择了错误的初始猜测，因此累积分布函数与经验分布函数大不相同。这意味着数据不遵循 EMWE 分布，而我的期望是数据必须遵循 EMWE 分布。

累积分布函数和经验分布函数如下图所示： 累积分布函数和经验分布函数

所以我这里的问题与我使用的初始猜测有关。这可能是较差的初始值。这里有没有人有解决方案来为这个数据集选择好的初始猜测？

4

0 回答 0