R:mle()错误的问题:非有限的有限差分值[2]

Question

我有一个简单的x, ydata.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)

这里出了什么问题？

Answer 1

你想让我做什么？据我所知，您有一个包含 96 个“风险”值的数据集，并且您希望使用威布尔来拟合它的分布。请注意，如果是这种情况，“天”根本不相关。您有一个无序的值向量。

\n\n

上图具有误导性。您计算dweibull()风险值。该图表明，这dweibull(risk)大致等于风险。这是一个与威布尔相当不同的主张，给定的参数非常适合。

\n\n

例如，以下是数据的分布：\n hist(mydata$risk, breaks=15)\n\n而参数在相关范围内的威布尔密度如下所示：curve((function(x) dweibull(x, shape=3, scale=0.155))(x), 0, 0.0014)\n

\n\n

因此这些分布非常不同。我想说你的经验分布是均匀的加上质量为零的分布，而不是威布尔分布。

\n\n

现在到最后一个问题：由于分布不太适合，优化器会遇到数值奇点。我不太mle()清楚，但稍微调整一下maxLik::maxLik()就会发现问题：

\n\n

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

\n\n

给你

\n\n

--------------------------------------------\nMaximum Likelihood estimation\nBHHH maximisation, 43 iterations\nReturn code 2: successive function values within tolerance limit\nLog-Likelihood: 682.743 \n2  free parameters\nEstimates:\n        Estimate Std. error t value Pr(> t)    \nKappa  0.4849129  0.0473720  10.236 < 2e-16 ***\nLambda 0.0002953  0.0001028   2.873 0.00407 ** \n---\nSignif. codes:  0 \xe2\x80\x98***\xe2\x80\x99 0.001 \xe2\x80\x98**\xe2\x80\x99 0.01 \xe2\x80\x98*\xe2\x80\x99 0.05 \xe2\x80\x98.\xe2\x80\x99 0.1 \xe2\x80\x98 \xe2\x80\x99 1\n--------------------------------------------\n

\n\n

请注意，我做了一项重大更改：sum从对数似然中删除，并使用 BHHH 优化器。这通常比基于单个总似然度的优化更稳定。您还应该认真考虑为估计编写解析导数。

\n\n

您可以检查现在的分布看起来更加相似。

\n