Sarhan and Apaloo (2013) によって導入された Exponentiated Modified Weibull Extension (EMWE) 分布の 4 つのパラメーターを最尤推定 (MLE) で推定したいと思います。この分布は、信頼性および生存分析で使用され、浴槽ハザード関数を使用してデータセットを分析します。
4 つのパラメーターの対数尤度関数の 1 次導関数が陰解を与えるため、ニュートン-ラフソン反復法を続けようとしました。私の主な焦点は、浴槽ハザード関数を取得するための正しい初期推定を選択する方法であるため、この方法の適切な初期推定を決定することに慣れていないだけです。私は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 分布に従う必要があります。
累積分布関数と経験分布関数を以下に示します。 累積分布関数と経験分布関数
したがって、ここでの問題は、私が使用する最初の推測に関連しています。初期値が悪いのかもしれません。このデータセットの適切な初期推測を選択する際に解決策を持っている人はいますか?