如何找到下面函数中参数的 MLE 估计值?

How can I find MLE estimates for the parameters in the function below?

我有一个新函数,我想估计它的参数; a,b,alpha,vartheta 使用 MLE。我如何进行估算?

EMHL<-function(a,b,alpha,vartheta) {
  (2*a*b*alpha*vartheta*
    (x^(vartheta-1))* exp(-x^vartheta) *
    ((1-exp(-x^vartheta))^(a-1)) * 
    ((1 - (1 -( 1 - exp(-x^vartheta))^a)^b)^(alpha-1)) * 
    (1- (1 - exp(-x^vartheta))^a )^(-b*(alpha+1))
}

对于给定的数据集

x<- c(1.1, 1.4, 1.3, 1.7,1.9, 1.8, 1.6, 2.2, 
      1.7, 2.7, 4.1, 1.8, 1.5, 1.2, 1.4, 3, 
      1.7, 2.3, 1.6, 2.0)

问题中显示的函数有语法错误,所以我们使用了最后注释中的那个。

如果该函数是一个密度函数并且您想最小化相应的负对数似然,那么尝试几个不同的起始值似乎会导致收敛。

(对于问题下方评论中给出的 x,list(a = 1, b = 1, alpha = 575, vartheta = 0.01) 似乎可以作为起始值。)

NLL <- function(par) -sum(log(EMHL(par[1], par[2], par[3], par[4])))
st <- list(a = 1, b = 1, alpha = 525, vartheta = .2)
res <- optim(st, NLL); res

给予:

$par
           a            b        alpha     vartheta 
8.845296e-01 1.211526e+00 5.315759e+02 1.326975e-03 

$value
[1] -10904.36

$counts
function gradient 
     327       NA 

$convergence
[1] 0

$message
NULL

L-BFGS-B

将L-BFGS-B与上述函数一起使用会出现问题,但如果我们限制参数就可以获得答案。例如,这会在所有约束都处于活动状态时收敛,即解在可行域的边界上。请注意,从边界估计得出的检验和 p 值可能无效。

optim(c(1, 1, 1, 1), NLL, method = "L-BFGS-B", lower = 0.01, upper = 2)

其他密度

另一种可能是使用不同的分布。 descdist 制作的 Cullen & Frey 图表明伽马分布可能接近

library(fitdistrplus)
descdist(x)

(图后续)

或者我们可以试试 flexsurv 包中的广义伽玛 (dgengamma)。

library(bbmle)
library(flexsurv)
NLLgeng <- function(mu, sigma, Q) -sum(dgengamma(x, mu, sigma, Q, log = TRUE))
m <- mle2(NLLgeng, list(mu = 1, sigma = 1, Q = 1), optimizer = "nlminb")
summary(m)

library(fitdistrplus)
fit <- fitdist(x, "gengamma", start = list(mu = 1, sigma = 1, Q = 1))
summary(fit)
plot(fit)

备注

EMHL<-function(a,b,alpha,vartheta) {
  2 * a * b * alpha * vartheta * 
  (x^(vartheta-1))* exp(-x^vartheta) *((1-exp(-x^vartheta))^(a-1)) * 
  ((1 - (1 -( 1 - exp(-x^vartheta))^a)^b)^(alpha-1)) * 
  (1- (1 - exp(-x^vartheta))^a )^(-b*(alpha+1))
}
x<- c(1.1, 1.4, 1.3, 1.7,1.9, 1.8, 1.6, 2.2, 1.7, 2.7, 4.1, 1.8,
   1.5, 1.2, 1.4, 3, 1.7, 2.3, 1.6, 2.0)