时间戳

首页 > 解决方案 > 在R中获取伽马分布参数（尤其是位置参数）的MLE

r - 在R中获取伽马分布参数（尤其是位置参数）的MLE

问题描述

伽玛分布pdf

嗨，我想手动估计伽马分布参数！我知道很多 R 函数来估计形状和比例参数，但似乎很难找到关于估计位置参数的代码。

x <- c(108,91,62,59,84,60,71,105,70,69,66,65,78,83,82,68,107,68,68,69,80,
       75,89,68,64,68,70,57,62,87,51,55,56,57,75,98,60,68,81,47,76,48,63,
       58,40,62,61,58,38,40,45,68,56,64,49,53,50,39,54,47,37,50,54,70,49,
       57,52,47,43,52,57,46,63,56,50,51,50,42,46,56,52,59,45,50,59,44,52,
       54,53,63,45,56,55,53,56,46,45,49,63,50,41,42,53,50,58,50,37,53,58,
       49,53,51,64,44,53,53,55,43,50,60,51,55,56,52,51,45,49,51,63,48,51,
       60,45,40,50,66,62,69,53,54,49,47,63,55,62,57,58,51,50,57,62,45,47,
       52,35,41,53,48,59,45,41,52,36,84,62,31,41,48,47,50,50,57,53,37,46,
       41,56,51,39,59,53,51,49,45,42,32,55,34,43,35,48,33,41,38,57,37,40,
       34,44,43,62,36,41,51,48,31,28,33,35,48,31)
# estimate shape and scale parameter
gamma_likelihood <- function(para){
  sum (  (para[2] -1)*log(x) - para[2]*log(para[1]) - log(gamma(para[2])) - x/para[1] + 1/para[1])
}

MLE = optim(c(10,10), 
            fn = gamma_likelihood, 
            method = "L-BFGS-B", 
            lower = 0.00001, 
            control = list(fnscale = -1), 
            hessian = T 
)
MLE$par


# estimate location, shape and scale parameter
gamma_likelihood <- function(para){
  x = x[x > para[1]]
  sum (  (para[3] -1)*log(x - para[1]) - para[3]*log(para[2]) - 
           log(gamma(para[3])) - x/para[2] + para[1]/para[2] )
}

MLE = optim(c(23,6,7), 
            fn = gamma_likelihood,
            method = 'L-BFGS-B',
            lower = 0.00000001,
            control = list(fnscale = -1)
)
MLE$par

这是我的代码，我可以估计形状和比例参数。

但是，在将位置参数添加到对数似然度时。结果似乎不正确。TRUE 参数是 c(21.4, 5.47, 6.0)。

标签： rmle

解决方案

如果您的任何观察值小于或等于您的位置参数，则该 lambda 值的全部可能性必须为 0（请记住，它是参数的函数，而不是观察值）。

x = x[x > para[1]]正在削减对特定位置参数没有意义的观察，使您的函数返回一个有效数字，-Inf如果任何一个x“无效”，它应该返回，因为你有0可能。

这是对数似然函数的更正版本：

# estimate location, shape and scale parameter
gamma_likelihood <- function(para){
  if(min(x) < para[1]) return(-Inf)
  sum (  (para[3] -1)*log(x - para[1]) - para[3]*log(para[2]) - 
           log(gamma(para[3])) - x/para[2] + para[1]/para[2] )
}

MLE = optim(c(23,6,7), 
            fn = gamma_likelihood,
            method = 'L-BFGS-B',
            lower = 0.00000001,
            control = list(fnscale = -1)
)
MLE$par

结果：[1] 21.161109 5.394343 6.136862

推荐阅读