有没有一种简单的方法来计算 R 中参数的最大似然估计？答案

【问题标题】：Is there a simple way to calculate the maximum likelihood estimate of a parameter in R?有没有一种简单的方法来计算 R 中参数的最大似然估计？
【发布时间】：2021-12-18 19:51:13
【问题描述】：

我正在尝试计算 R 中泊松分布的 MLE。R 中是否有一个函数可以让我们这样做（例如，我知道 Stata 有一个 mlexp 函数，可以让我们很容易地进行这个计算）。我看到在 univariateML 包中有一个 mlexp 函数用于指数分布。话虽这么说，是否有一个命令允许这不仅仅是指数分布？

【问题讨论】：

泊松分布的 MLE 就是 mean。你可以在这里查看证明：statlect.com/fundamentals-of-statistics/…
是的，这是真的。话虽如此，我正在寻找可以找到 MLE 的 R 命令。我同样可以要求一个命令，该命令将从 Gamma 分布中找到参数的 MLE。

标签： r log-likelihood

【解决方案1】：

fitdistrplus 包可能对您想要的有所帮助，但它无法处理更复杂的发行版。 stats4 的 mle() 函数包也可能有帮助。我相信bbmle 包比其他选项更通用，但我自己并没有真正使用过。

【讨论】：

【解决方案2】：

这是我用于自定义分发的函数。类似于stats4::mle。它确实假定 PDF 参数名为 x（与所有经典发行版的 R 实现一样）。

mle <- function(data, fun, init, logarg = TRUE, lower = -Inf, upper = Inf) {
  if (logarg) {
    fnll <- function(...) {
      -sum(do.call(fun, l <- c(x = list(data), log = TRUE, as.list(...))))
    }
  } else {
    fnll <- function(...) {
      return(-sum(log(do.call(fun, c(x = list(data), as.list(...))))))
    }
  }
  
  params <- optim(init, fnll, method = "L-BFGS-B", lower = lower, upper = upper)$par
  names(params) <- formalArgs(fun)[-1][seq_along(init)]
  return(params)
}

# optim will check the boundaries
# set lim0 and lim1 for help in setting bounds
lim0 <- .Machine$double.eps
lim1 <- 1 + lim0

# Poisson distribution
data <- rpois(1e5, 14)
rbind(trueML = c(lambda = mean(data)),
      mle = mle(data, dpois, 1, lower = lim0))
#>          lambda
#> trueML 13.99387
#> mle    13.99387

# normal distribution
data <- rnorm(1e5, -2, 3)
rbind(trueML = c(mean = mean(data), sd = sd(data)*sqrt((length(data) - 1)/length(data))),
      mle = mle(data, dnorm, 0:1, lower = c(-Inf, lim0)))
#>             mean       sd
#> trueML -2.002658 2.993441
#> mle    -2.002657 2.993442

# gamma distribution
data <- rgamma(1e5, 0.5, 0.1)
c(mle = mle(data, dgamma,
            init = c(mean(data)^2/var(data), mean(data)/var(data)),
            lower = rep(lim0, 2)))
#> mle.shape  mle.rate 
#> 0.5007139 0.1003400

# triangular distribution
dtri <- function(x, a, b, c) {
  if (a > b) {a <- (b - a) + (b <- a)}
  if (b > c) {c <- (b - c) + (b <- c)}
  blna <- x < b
  p <- numeric(length(x))
  p[blna] <- 2*(x[blna] - a)/(c - a)/(b - a)
  p[!blna] <- 2*(c - x[!blna])/(c - a)/(c - b)
  return(p)
}

rtri <- function(n, a, b, c) {
  if (a > b) {a <- (b - a) + (b <- a)}
  if (b > c) {c <- (b - c) + (b <- c)}
  fb <- (b - a)/(c - a)
  U <- runif(n)
  blna <- U < fb
  r <-numeric(n)
  r[blna] <- a + sqrt(U[blna]*(c - a)*(b - a))
  r[!blna] <- c - sqrt((1 - U[!blna])*(c - a)*(c - b))
  return(r)
}

data <- rtri(1e5, -6, -3, 3)
mind <- min(data); maxd <- max(data)
# set logarg to FALSE because dtri doesn't have a log argument
c(mle = mle(data, dtri, logarg = FALSE,
            init = c(mind - abs(mind)*lim0, median(data), maxd + abs(maxd)*lim0),
            lower = c(-Inf, min(data), maxd + abs(maxd)*lim0),
            upper = c(mind - abs(mind)*lim0, max(data), Inf)))
#>     mle.a     mle.b     mle.c 
#> -6.003666 -3.000262  2.994473

Created on 2021-11-04 by the reprex package (v2.0.1)

【讨论】：