Максимальное правдоподобие 'mle' обеспечивает стандартные ошибки NaN

Question

Я использую функцию mle из пакета stats4 для оценки моей модели с максимальной вероятностью.Я получаю разумные оценки параметров, но некоторые стандартные ошибки: NaN:

require(stats4)
df = read.csv("https://pastebin.com/raw/uDhQj4LW", stringsAsFactors = F, sep="\t")

LL <- function(b0, b1, b2, b3) {

  xmu = b0 + b1 * df$brentcpl
  xsi = b2 + b3 * df$brentcpl
  lls = log(1 / df$price / xsi / sqrt (2*pi)) - (log(df$price) - xmu)^2/ 2 / xsi^2
  lls[is.nan(lls)] = -100 - (b0 + b1 + b2 + b3)^2
  #print(paste0("b0 ", b0, "b1 ", b1, " b2 ", b2, " b3 ", b3, " mu ", min(xmu), " xsi ", min(xsi)))
  #print(sum(lls))
  -sum(lls) 
}

estres = mle(LL, start = list(b0 = 4.5171634,
                              b1 = .00949832,
                              b2 = .03892378,
                              b3 = .00069329     
))

summary(estres)
Maximum likelihood estimation

Call:
mle(minuslogl = LL, start = list(b0 = 4.5171634, b1 = 0.00949832, 
    b2 = 0.03892378, b3 = 0.00069329))

Coefficients:
       Estimate   Std. Error
b0  4.495186481 0.0193639611
b1  0.007735690 0.0006672378
b2 -0.012258204 0.0008411522
b3  0.001435819          NaN

-2 log L: 1205.609 
Warning message:
In sqrt(diag(object@vcov)) : NaNs produced

Что здесь может пойти не так?Как я могу это исправить?

Answer 1

Дисперсионно-ковариационная матрица модели может быть получена с использованием:

estres@vcov
#               b0            b1            b2            b3
# b0  3.749630e-04 -1.286236e-05  5.564893e-06 -9.836086e-11
# b1 -1.286236e-05  4.452062e-07 -1.409643e-07  3.231134e-12
# b2  5.564893e-06 -1.409643e-07  7.075369e-07  1.908407e-11
# b3 -9.836086e-11  3.231134e-12  1.908407e-11 -6.241188e-17

Дисперсия b3 (в положении [4,4]) очень мала, но отрицательна.
Вmle дисперсионно-ковариационная матрица была оценена путем вычисления обратного гессиана:

solve(oout$hessian)

, где oout$hessian равно

#               b0            b1            b2            b3
# b0  2.752547e+05  8.091507e+06 -1.634045e+04  -19890328389
# b1  8.091507e+06  2.398502e+08 -7.378089e+05 -560469634334
# b2 -1.634045e+04 -7.378089e+05  1.871352e+05   44776772461
# b3 -1.989033e+10 -5.604696e+11  4.477677e+10    7443899452

Вы можете попытаться решить проблему, изменившкала df$brentcpl:

df$brentcpl <- df$brentcpl/100
estres = mle(LL, start = list(b0 = 5,
                              b1 = .25,
                              b2 = .25,
                              b3 = .25     
))
summary(estres)

# Maximum likelihood estimation
# Call:
# mle(minuslogl = LL, start = list(b0 = 5, b1 = 0.25, b2 = 0.25, 
#     b3 = 0.25)    
# Coefficients:
#       Estimate  Std. Error
# b0  4.54064769 0.013846337
# b1  0.61318128 0.044016385
# b2  0.06951554 0.008852508
# b3 -0.16738634 0.028168022
# -2 log L: 1078.542