Implementation of dmultinom ?

[tiboloic]

For revising a paper I am porting some TMB code to RTMB. It is a random effect model with a mutltinomial response, but I cannot get it to work. Since I didn't find a implementation of dmultinom in distributions.R, I have also tried my own implementation of dmultinom without success.

Here is a snippet of code:

# RTMB multinomial test

n = 10
p = 20

# simulate data
obs <- t(rmultinom(n, 100, prob=p:1))

par <- list(
  mu = 2,
  sig = 1,
  bs = rep(2, n)
)

nLL <- function(params) {

  # own implementation of log probability mass function for multinomial
  dm <- function(x,p) {
    lgamma(sum(x) + 1) - sum(lgamma(x+1)) + sum(x*log(p))
  }

  getAll(par)

  nLL <- 0.0

  # random effects
  nLL <- nLL - sum(dnorm(bs, mu, sig, TRUE))

  for (i in 1:nrow(obs)) {
    beta <- bs[i]
    preds <- 1/(1:p)^beta
    preds <- preds/sum(preds)

    #nLL <- nLL - dmultinom(obs[i,], prob = preds, log = TRUE)
    nLL <- nLL - dm(obs[i,], preds)
  }
  nLL
}

obj <- MakeADFun(nLL, par, random = c("bs"))

fit <- nlminb(obj$par, obj$fn, obj$gr)

[kaskr]

dmultinom is available in RTMB (see ?Distributions). The problem is in your nLL where getAll(par) should actually be getAll(params) :)

Edit: I think getAll could be improved to avoid this kind of mistake, by requiring that the argument exists in the calling frame.

[tiboloic]

Apologies for the trivial error. Thank you so much for looking into it.