Support for Multinomial distribution

[erumpler]

Hi and thanks for a great package !

I noticed the DSL doesn’t support the Multinomial distribution yet. Are there any plans to add this soon?

monty_dsl_distributions()
#> name args
#> 1 Beta a, b
#> 2 BetaBinomial size, pr....
#> 3 BetaBinomial size, a, b
#> 4 Binomial size, prob
#> 5 Cauchy location....
#> 6 Exponential rate
#> 7 Exponential mean
#> 8 Gamma shape, rate
#> 9 Gamma shape, scale
#> 10 Hypergeometric m, n, k
#> 11 LogNormal meanlog,....
#> 12 NegativeBinomial size, prob
#> 13 NegativeBinomial size, mu
#> 14 Normal mean, sd
#> 15 Poisson lambda
#> 16 TruncatedNormal mean, sd....
#> 17 Uniform min, max
#> 18 Weibull shape, scale

Thank you !

[edknock]

Hi @erumpler. Thanks for the question - can you clarify which DSL are you interested in using the multinomial distribution in? The monty DSL or odin (odin2) DSL? They are different but related - e.g. the odin DSL uses the distribution functions built into the monty package.

[erumpler]

Hi @edknock and thank you for your message. I'd like to use a multinomial distribution to compare a model output to data and compute a likelihood in odin2::odin.

[edknock]

So multinomial distributions is something we've discussed for a while but need to settle on syntax (we currently support no multivariate distributions). With regards to using it for likelihood components, the extra issue there is that we currently only support scalar datastreams, so we'd need support for array data (and how they'd be inputted in a data frame).

So the multinomial distribution is not currently explicitly available, and it may be in future, though we can't really say when.

It's not all bad news however as a multinomial distribution can be expressed in terms of nested, conditional binomial distributions, which shouldn't be too bad to express if the dimensions of the multinomial are not too large.

This is illustrated here for a 4d multinomial with a toy example:

mod <- odin2::odin({
  update(X[]) <- Poisson(10)
  initial(X[]) <- 0
  dim(X) <- 4
  
  y1 <- data()
  y2 <- data()
  y3 <- data()
  y4 <- data()
  
  prob[] <- X[i] / sum(X)
  dim(prob) <- dim(X)
  
  y1 ~ Binomial(y1 + y2 + y3 + y4, prob[1])
  y2 ~ Binomial(y2 + y3 + y4, prob[2] / sum(prob[2:4]))
  y3 ~ Binomial(y3 + y4, prob[3] / sum(prob[3:4]))
})

sys <- dust2::dust_system_create(mod, seed = 1)

dust2::dust_system_run_to_time(sys, 1)

y <- c(8, 11, 5, 6)
d <- list(y1 = y[1], y2 = y[2], y3 = y[3], y4 = y[4])
dust2::dust_system_compare_data(sys, d)

which if you compare to the following you should find you get the same result

X <- dust_system_state(sys)
prob <- X / sum(X)
dmultinom(y, prob = prob, log = TRUE)

In creating a data frame to go dust_filter_create, it would then need columns for y1, y2, y3 and y4.

An alternative formulation would be to instead use the cumulative sum values over the multinomial vector as the data, you'll find this provides the same outcome

mod2 <- odin2::odin({
  update(X[]) <- Poisson(10)
  initial(X[]) <- 0
  dim(X) <- 4
  
  z1 <- data()
  z2 <- data()
  z3 <- data()
  z4 <- data()
  
  prob[] <- X[i] / sum(X)
  dim(prob) <- dim(X)
  
  z1 ~ Binomial(z2, prob[1] / sum(prob[1:2]))
  z2 ~ Binomial(z3, sum(prob[1:2]) / sum(prob[1:3]))
  z3 ~ Binomial(z4, sum(prob[1:3]) / sum(prob[1:4]))
})

z <- cumsum(y)

sys2 <- dust2::dust_system_create(mod2, seed = 1)

dust2::dust_system_run_to_time(sys2, 1)

d2 <- list(z1 = z[1], z2 = z[2], z3 = z[3], z4 = z[4])
dust2::dust_system_compare_data(sys2, d2)

[edknock]

If you need a multinomial distribution for random draws in the odin model, that's slightly easier as the array support there should allow it to be done in two lines with nested binomials:

y[1] <- Binomial(size, prob[1])
y[2:length(y)] <- Binomial(size - sum(y[1:(i - 1)]), prob[i] / sum(prob[i:length(y)]))

[erumpler]

Great, thank you very much @edknock !

[edknock]

@erumpler Support for array data is forthcoming mrc-ide/odin2#151

Multinomial distribution might follow from this but not necessarily soon. However, support for array data should be able to simplify things so a multinomial likelihood component could be expressed in terms of a loop of nested binomials in one line I think!

[edknock]

Hey @erumpler, you can now use vector data if you update to the latest version of odin2! For guidance for how to set this up in your odin code and your data frame, see here https://mrc-ide.github.io/odin2/articles/fitting.html#comparing-against-vectors-of-data

Revisiting the toy multinomial example, this can now be written as a nested binomial loop in one line!

mod <- odin2::odin({
  update(X[]) <- Poisson(10)
  initial(X[]) <- 0
  dim(X) <- n

  n <- parameter(4)

  y[] <- data()
  dim(y) <- n

  prob[] <- X[i] / sum(X)
  dim(prob) <- dim(X)

  y[1:(n - 1)] ~ Binomial(sum(y[i:n]), prob[i] / sum(prob[i:n]))
})

sys <- dust2::dust_system_create(mod, seed = 1)

dust2::dust_system_run_to_time(sys, 1)

y <- c(8, 11, 5, 6)
d <- list(y = y)
dust2::dust_system_compare_data(sys, d)

and you will see that it proudces the same as

X <- dust2::dust_system_state(sys)
prob <- X / sum(X)
dmultinom(y, prob = prob, log = TRUE)