warnings in gamma distribution (*ml* vignette)

[thibautjombart]

See the ml vignette, warnings when using ML estimation to fit a discretised gamma distribution. Reproducible code:

shape <- 10
rate <- 1.2
y <- distcrete("gamma", interval = 1L, shape, rate)
set.seed(1)
sim2 <- y$r(300)

ll2 <- function(param) {
   d <- distcrete("gamma", interval = 1L, param[1], param[2])$d
   sum(d(sim2, log = TRUE))
}

dev2 <- function(param) -2 * ll2(param)

optim(c(1,1), dev2)

Gives:

## Warning in cdf(x, ..., log.p = log): NaNs produced

## Warning in cdf(x, ..., log.p = log): NaNs produced

## Warning in cdf(x, ..., log.p = log): NaNs produced

## Warning in cdf(x, ..., log.p = log): NaNs produced

## Warning in cdf(x, ..., log.p = log): NaNs produced

## Warning in cdf(x, ..., log.p = log): NaNs produced

[richfitz]

This warning is probably reasonable, IMO. It comes from pgamma itself:

pgamma(1, 2, -0.1, log.p = TRUE)
# [1] NaN
# Warning message:
# In pgamma(1, 2, -0.1, log.p = TRUE) : NaNs produced

Suppressing that warning seems suboptimal, partly because mimicking the behaviour of the underlying probability functions is easiest to think about. It would be possible to add an option to suppress warnings but there is a small overhead there and it's just as easy to do in user code;

ll2 <- function(param) {
  d <- distcrete("gamma", interval = 1L, param[1], param[2])$d
  sum(suppressWarnings(d(sim2, log = TRUE)))
}

or, (probably better) using things we already know about the distribution

ll2 <- function(param) {
  if (param[[2]] <= 0) { # really, the '=' part is contentious...
    -Inf
  } else {
    d <- distcrete("gamma", interval = 1L, param[1], param[2])$d
    sum(d(sim2, log = TRUE))
  }
}