Support for multinomial draws

[annecori]

At the moment we are using nested binomials to model multinomial draws, e.g.:

# Compute the new infections with multiple strains using nested binomials
n_inf[1] <- rbinom(S, p_inf[1])
n_inf[2:n_strains] <- rbinom(S - sum(n_inf[1:(i - 1)]), p_inf[i])

It would be great if we could write something simpler e.g.:

# Compute the new infections with multiple strains using nested binomials
n_inf[1:n_strains] <- rmultinom(S, p_inf[i])

[richfitz]

## Vector form
n_inf[1] <- rbinom(S, p_inf[1])
n_inf[2:n_strains] <- rbinom(S - sum(n_inf[1:(i - 1)]), p_inf[i])

## Matrix form, assigning into columns
n_inf[1, ] <- rbinom(S[j], p_inf[1, j])
n_inf[2:n_strains, ] <- rbinom(S[j] - sum(n_inf[1:(i - 1)], j), p_inf[i, j])

## Matrix form, assigning into rows
n_inf[, 1] <- rbinom(S[i], p_inf[i, 1]) # could be p_inf[1, i]
n_inf[, 2:n_strains] <- rbinom(S[i] - sum(n_inf[i, 1:(i - 1)]), p_inf[i, j])

[cwhittaker1000]

Whilst playing around with multinomial draws using something based on #213 (comment), I think I came across an issue, centred around the need to update the probability used in each binomial draw once you've assigned to an element of n_inf[]. I think this is what the correct version of the code in #213 (comment) should look like!

total <- rbinom(S, sum(p_inf))
n_inf[1] <- rbinom(total, p_inf[1]/sum(p_inf)) 
n_inf[2:n_strains] <- rbinom(total - sum(n_inf[1:(i - 1)]), p_inf[i]/sum(p_inf[i:n_strains]))

When doing the nested binomials as an approximation for the multinomial, each time you do one of the binomial draws (and assign to an element of n_inf[]) I think you need to update the probability used in the next binomial draw, so that it's relative to the sum of all the remaining probabilities you have left (i.e. p_inf[i:n_strains]).

Note that there's some nuance around whether you're trying to assign the entirety of S into n_inf[]'s elements (which is not what #213 (comment) is doing; or whether you're aiming to assig all of S into n_inf[]'s elements (which more closely resembles a classical multinomial draw). But the above code should be able to handle both instances.