I am trying to calculate the MLE of a poisson distribution in R. Is there a function in R that allows us to do this (eg. I know Stata has a mlexp function that allows us to make this calculation quite easily). I see that there is a mlexp function in the univariateML package for the exponential distribution. That being said, is there a command that allows this for more than just exponential distributions?

CodePudding user response：

The fitdistrplus package might help with what you want, though it can't handle more complicated distributions. The mle() function of the stats4 package might also help. I believe the bbmle package is more general than the other options, but I haven't really used it myself.

CodePudding user response：

Here's a function I've used for custom distributions. Similar to stats4::mle. It does assume that the PDF argument is named x (as with the R implementaion of all the classic distributions).

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)