錯誤最大似然估計，使用R

Question

我目前正在整理用下面的pdf關於對Sarhan和Apaloo（2013）引入指數化修改威布爾擴展（EMWE）分佈參數估計我的本科畢業論文：

f(x,theta)=theta[1]*theta[2]*theta[3]*((x/theta[4])^(theta[2]-1))*(exp(((x/theta[4])^theta[2])+(theta[1]*theta[4]*(1-(exp(x/theta[4])^theta[2])))))*(1-(exp(theta[1]*theta[4]*(1-(exp(x/theta[4])^theta[2])))))^(theta[3]-1) 

這分佈具有使用最大似然估計估計的四個參數。由於隱式估計參數，我試圖繼續使用Newton-Raphson迭代方法。對於我的計算過程，我正在使用統計軟件R語言。我使用的包是「maxLik」，Newton-Raphson方法的初始值是(theta [1] = 0.00007181, theta [2] = 3,148, theta [3] = 0.145, theta [4] = 49.05)。

這是對數似然函數：

l(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])))))) 

但在有R語言的幫助參數估計的過程中，我絕境因爲我所獲得的結果的不與在參考紙估計結果相似我用。這是下述R語言的語法使用：

xi<-c(0.1,0.2,1,1,1,1,1,2,3,6,7,11,12,18,18,18,18,18,21,32,36,40,45,46,47,50,55,60,63,63,67,67,67,67,72,75,79,82,82,83,84,84,84,85,85,85,85,85,86,86); 
n <-length (xi); 
parameter <-function (theta, xi) { 
logL<-(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]))))))) 
return (-logL) 
}; 
library(maxLik); 
output <-maxLik (parameter, start = c (0.00007181,3.148,0.145,49.05), xi = xi); 

基於該語法，參數估計我得到的結果是：

theta [1] = 4.785855e-03 
theta [2] = 1.759048e-04 
theta [3] = 2.983679e + 04 
theta [4] = 9.139192e + 02 

雖然在紙上屬於Sarhan和Apaloo（2013），所述結果應該如下：

theta [1] = 2.506924e-06 
theta [2] = 3.148000e + 00 
theta [3] = 1.450000e-01 
theta [4] = 4.905000e + 01 

我很困惑我的錯誤在上面的程序中。以前，如果我困擾你們，我很抱歉。我非常感謝你幫助我完成本科畢業論文。很快我會介紹這篇論文，並發現了很多的僵局。我真的希望得到你的幫助，不管它多麼小的幫助，它將不勝感激。非常感謝你

• 說實話，我很抱歉我的英語語法不好。我不說英語

Answer 1

你的可能性函數有些問題。我無法閱讀它，但請注意，maxLik最大化目標函數，因此您必須返回loglik，而不是-loglik。我以一種更可讀的形式重新編寫了它（請參閱Sarhan & Apaloo 2013）（對不起，請爲參數命名，添加空間，將長方程拆分爲多行......），而且我也不想使用名稱「參數「爲對數似然函數...

loglik <-function(theta, xi) { 
    lambda <- theta[1] 
    beta <- theta[2] 
    gamma <- theta[3] 
    alpha <- theta[4] 
    xi.a <- xi/alpha 
    A <- log(lambda) + log(beta) + log(gamma) + (beta - 1)*log(xi.a) 
    LA1 <- lambda*alpha*(1 - exp(xi.a^beta)) 
    B <- xi.a^beta + LA1 
    C <- log(1 - exp(LA1)) 
    logL <- A + B + (gamma - 1)*C 
    return(logL) 
} 

library(maxLik); 
start <- c(2.506924e-6, 3.148,0.145,49.05) 
m <- maxLik(loglik, start=start, xi = xi); 

它有些作品。主要問題似乎是數值不穩定。不同的優化方法發揮，特別是BFGS似乎讓你很接近：

summary(maxLik(loglik, start = start, method="bfgs", xi = xi)) 
-------------------------------------------- 
Maximum Likelihood estimation 
BFGS maximization, 337 iterations 
Return code 0: successful convergence 
Log-Likelihood: -213.3168 
4 free parameters 
Estimates: 
     Estimate Std. error t value Pr(> t)  
[1,] 1.213e-05 5.976e-06 2.030 0.0423 * 
[2,] 3.133e+00 3.709e-02 84.456 < 2e-16 *** 
[3,] 1.255e-01 1.897e-02 6.615 3.72e-11 *** 
[4,] 4.496e+01   NA  NA  NA  
--- 
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

另外，如果你修復第一個參數，那麼你得到的確切值與BHHH：

summary(maxLik(loglik, start = start, method="bhhh", xi = xi, fixed=1)) 
-------------------------------------------- 
Maximum Likelihood estimation 
BHHH maximisation, 13 iterations 
Return code 2: successive function values within tolerance limit 
Log-Likelihood: -213.5116 
3 free parameters 
Estimates: 
     Estimate Std. error t value Pr(> t)  
[1,] 2.507e-06 0.000e+00  NA  NA  
[2,] 3.091e+00 9.577e-03 322.74 < 2e-16 *** 
[3,] 1.153e-01 2.234e-02 5.16 2.46e-07 *** 
[4,] 4.189e+01 1.592e+00 26.32 < 2e-16 *** 

它提示剩下的問題與源自第一組分（λ）的數值不穩定性有關。我可以建議兩種補救措施：

• 向maxLik函數提供分析梯度。我知道，這是一個地獄，但可能只是提供lambda就足夠了。
• 重新參數化問題。即使指定lambda <- theta[1]/1e6，並相應地調整起始值，似乎也會改善收斂性。

請注意，我刪除了似然函數的求和：現在你也可以使用BHHH方法，它通常比NR更強大。