search

mrc-ide / odin

ᚩ A DSL for describing and solving differential equations in R
https://mrc-ide.github.io/odin
Other
106 stars 13 forks source link

use rmultinom with arrays #255

Open qu-cheng opened 2 years ago

qu-cheng commented 2 years ago

Dear odin team,

I've reat that rmultinom can not work as y[1, ] <- rmultinom(10, p[i, ]), but what if I have a metapopulation model? In that case, I'll have multiple sizes. I need to do something like n_Eout[i,] <- rmultinom(n_Eout[i], p_Eout_all), where p_Eout_all is a vector with the same number of elements as the number of columns of n_Eout.

Thanks, Qu

richfitz commented 2 years ago

Please see https://github.com/mrc-ide/odin/issues/134 and https://github.com/mrc-ide/odin/issues/213 for discussion on this. One day we'll sort it out, but it's unlikely to be soon I'm afraid as it needs some new syntax to deal with the ambiguity, plus some work to avoid additional vector allocations.

qu-cheng commented 2 years ago

Thanks a lot! Using rbinom multiple times could do the trick, although a little tedious

richfitz commented 2 years ago

Agreed, quite tedious. In #213 there's a trick to re-write the multinomial draw as a set of binomial draws - that will allow expansion to work on a matrix (or higher-order structure). We use this on some large models and it will be just as efficient.

Also: if you're using odin for stochastic models, you should look at our dust package (via odin.dust) as that will compile models that can be run in parallel (see this paper, this vignette and this walkhrough)

thibautjombart commented 1 year ago

+1 to implementing this, alongside hypergeometric distributions. I am not sure an alternative such as the nested binomials outlined in https://github.com/mrc-ide/odin/issues/213 exists for the latter, but would love to hear about a similar trick if it exists.

richfitz commented 1 year ago

yes, same trick would likely work for multivariate hypergeometric (note that R does not actually support this distribution). The implementation at present is here: https://github.com/mrc-ide/odin/blob/master/inst/library.c#L358-L375 so one could adapt that same logic iterating across the axes of your array.

While I recognise that this feature would be of use, I would not expect it to be implemented anytime soon based on our current priorities (certainly not this calendar year). It's likely a lot of work and something that can at least be worked around, unlike some of the other things we hope to get done this year.

thibautjombart commented 1 year ago

Makes sense, thanks.