R: проблема с ошибкой mle(): неконечное конечно-разностное значение [2]

У меня есть простой x, y data.frame.

mydata <- data.frame(days = 1:96, risk = c(5e-09, 5e-09, 5e-09, 1e-08, 4e-08, 6e-08, 9e-08, 1.5e-07, 4.2e-07, 
                                           7.2e-07, 1.02e-06, 1.32e-06, 1.66e-06, 2.19e-06, 2.76e-06, 3.32e-06, 
                                           3.89e-06, 4.55e-06, 5.8e-06, 7.16e-06, 8.51e-06, 9.85e-06, 1.138e-05, 
                                           1.396e-05, 1.672e-05, 1.947e-05, 2.222e-05, 2.521e-05, 2.968e-05, 
                                           3.439e-05, 3.909e-05, 4.378e-05, 4.894e-05, 5.697e-05, 6.546e-05, 
                                           7.392e-05, 8.236e-05, 9.16e-05, 0.00010573, 0.00012063, 0.00013547, 
                                           0.00015025, 0.00016642, 0.00019127, 0.00021743, 0.00024343, 0.00026924, 
                                           0.00029818, 0.00034681, 0.00039832, 0.00044932, 0.00049976, 0.0005451, 
                                           0.00056293, 0.00057586, 0.00058838, 0.0006005, 0.00061562, 0.00065079, 
                                           0.00068845, 0.00072508, 0.00076062, 0.00079763, 0.00084886, 0.00090081, 
                                           0.0009507, 0.00099844, 0.00104427, 0.00108948, 0.00113175, 0.00117056, 
                                           0.00120576, 0.00123701, 0.00126253, 0.00128269, 0.00129757, 0.00130716, 
                                           0.00131291, 0.00132079, 0.0013216, 0.00131392, 0.00129806, 0.00127247, 
                                           0.00122689, 0.00117065, 0.00110696, 0.00103735, 0.00095951, 0.00085668, 
                                           0.0007517, 0.00065083, 0.000556, 0.0004669, 0.00037675, 0.00029625, 
                                           0.00093289))

Я думаю, что Weibull(3, 0.155) довольно хорошо подходит для моих данных, судя по графику ниже.

plot(1:96, dweibull(mydata$risk, shape = 3, scale = 0.155), type = "l", xlab = "days", ylab = "risk")
lines(mydata, type = "l", col = "grey")
legend("topleft", c("Data", "Estimate"), col = c("black", "grey"), lty = c(1, 1))

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

estimate <- function(kappa, lambda){
  -sum(dweibull(mydata$y, shape = kappa, scale = lambda, log = TRUE))
}

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

> mle(estimate, start = list(kappa = 3, lambda = 0.155))
Error in optim(start, f, method = method, hessian = TRUE, ...) : 
  non-finite finite-difference value [2]
In addition: There were 50 or more warnings (use warnings() to see the first 50)

Что здесь пошло не так?

r estimation statistics mle

Adrian 20.09.2017 источник

Ответы (1)

arrow_upward
2
arrow_downward

Что ты хочешь делать? Насколько я могу судить, у вас есть набор данных из 96 значений «риска», и вы хотите сопоставить его распределение с Вейбуллом. Обратите внимание, что «дни» вообще не имеют значения, если это так. У вас есть неупорядоченный вектор значений.

Цифра выше вводит в заблуждение. Вы рассчитываете dweibull() для значений риска. Цифра показывает, что dweibull(risk) примерно равно риску. Это несколько иное утверждение, чем вейбулл, с заданными параметрами, которые хорошо подходят.

например, вот распределение ваших данных: hist(mydata$risk, breaks=15) , а плотность вейбулла с вашими параметрами в соответствующем диапазоне выглядит так: curve((function(x) dweibull(x, shape=3, scale=0.155))(x), 0, 0.0014)

Следовательно, эти дистрибутивы очень разные. Я бы сказал, что ваши эмпирические распределения равномерны плюс масса в нуле, а не вейбулл.

Теперь к вашей последней проблеме: поскольку распределения плохо подходят, оптимизатор сталкивается с числовыми особенностями. Я не слишком хорошо знаю mle(), но с небольшими изменениями maxLik::maxLik() покажет проблему:

estimate <- function(par){
   Kappa <- par[1]
   Lambda <- par[2]
   dweibull(mydata$risk, shape = Kappa, scale = Lambda, log = TRUE)
}
summary(maxLik::maxLik(estimate, start=c(Kappa=3, Lambda=0.155), method="BHHH"))

дает тебе

--------------------------------------------
Maximum Likelihood estimation
BHHH maximisation, 43 iterations
Return code 2: successive function values within tolerance limit
Log-Likelihood: 682.743 
2  free parameters
Estimates:
        Estimate Std. error t value Pr(> t)    
Kappa  0.4849129  0.0473720  10.236 < 2e-16 ***
Lambda 0.0002953  0.0001028   2.873 0.00407 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
--------------------------------------------

Обратите внимание, что я сделал одно существенное изменение: удалил sum из вашего логарифмического правдоподобия и использовал оптимизатор BHHH. Обычно это более стабильно, чем оптимизация на основе одной суммированной вероятности. Вам также следует серьезно подумать о написании аналитических производных для оценки.

Вы можете убедиться, что теперь дистрибутивы выглядят намного более похожими.

Ott Toomet 22.09.2017

comment

Спасибо. Что касается вашего последнего комментария, как я могу проверить, что теперь дистрибутивы выглядят намного более похожими? Я попробовал plot(dweibull(1:95, shape = 0.4838894, scale = 0.0002961)), но это выглядит иначе, чем распределение моих данных? - Adrian; 02.10.2017

comment

Вы можете начать с чего-то простого, например, с qqplot. Вы также можете увидеть, насколько похожи моменты и другие характеристики этих двух распределений. И, наконец, вы можете рассчитать некоторое расстояние между этими двумя распределениями, например, расстояние Кульбака-Лейблера (которое совпадает с логарифмическим правдоподобием). Последний также будет работать для определенных тестов. - Ott Toomet; 02.10.2017

comment

Есть ли быстрый графический способ увидеть, похожи ли два дистрибутива? По линиям построения плотностей. - Adrian; 02.10.2017

comment

fitdistrplus::plotdist(risk, "weibull", list(shape=0.4838894, scale=0.0002961)), на мой взгляд, неплохо справляется со своей задачей. (risk — ваша переменная риска). - Ott Toomet; 02.10.2017

comment

Спасибо. Я хотел бы уточнить, что mydata на самом деле содержит плотности (у меня нет необработанных данных), поэтому я думаю, что мой первый график (где я сравнил плотности Вейбулла с mydata$risk) хорошо подходит. В вашем ответе мы технически рассматриваем плотности mydata$risk, то есть плотности плотностей... правильно? - Adrian; 02.10.2017