C++ Multinomial distribution

Question

I'm trying to code a multinomial algorithm that will basically apply a binomial distribution to each value of an input vector, knowing the values of all previous ones. It's aimed here to generate a new population for multiple alleles knowing the initial population.

To achieve this, I'm using this recursive algorithm :

This is what my code looks like right now :

void RandomNumbers::multinomial(std::vector<unsigned int>& alleleNumbers) {

/*  In this function we need two different records of the size.
 *  We need the size from the old populations, ( N - n1 - ... - nA )
 *  and we also need the size from the newly created population,
 *  ( N - k1 - ... - kA ).
 *  In order to achieve such a task, we'll use the integer "temp" to store
 *  the value n1 before modifying it to k1 and so on.
 *
 *
 */
double totalSize = 0;

for(auto n : alleleNumbers) totalSize+=n;

double newTotalSize(totalSize);

std::cout<< newTotalSize;

for(size_t i = 0; i < alleleNumbers.size(); ++i){
    size_t temp = alleleNumbers[i];
    alleleNumbers[i] = binomial(newTotalSize,
                (alleleNumbers[i])/(totalSize));
    newTotalSize-= alleleNumbers[i];
    totalSize = temp;  
  }
}

But I'm not sure at all about this, and I was wondering if there was an already existing multinomial algorithm of that kind...

Thank you very much.

There is no direct implementation of this in the standard library and since questions suggesting a library are off-topic here, have a look e.g. here: https://scicomp.stackexchange.com/questions/382/looking-for-c-c-implementations-of-sampling-from-multinomial-and-dirichlet-dis — , Nov 24 '18 at 16:15
There is [std::discrete_distribution](https://en.cppreference.com/w/cpp/numeric/random/discrete_distribution). — eozd, Nov 24 '18 at 16:16

score 1 · Answer 1 · answered Nov 24 '18 at 16:18

You could try using the GNU Scientific Library's gsl_ran_multinomial command.

The function is called as:

gsl_ran_multinomial (const gsl_rng * r, size_t K, unsigned int N, const double p[], unsigned int n[])

where (n_1, n_2, ..., n_K) are nonnegative integers with sum_{k=1}^K n_k = N, and (p_1, p_2, ..., p_K) is a probability distribution with sum(p_i) = 1. If the array p[K] is not normalized then its entries will be treated as weights and normalized appropriately. The arrays n[] and p[] must both be of length K.

The function implements the conditional binomial method from C.S. Davis's "The computer generation of multinomial random variates" (Comp. Stat. Data Anal. 16, 1993. link), so you could implement using that approach. Let me know if you need a copy of the paper.

Hey,Thank you for this amazing response, since then I installed GSL, however, I cannot seem to manage to compile with it, and I only get a "gsl_rng does not name a type"... I will try to fix my problem, but thanks a lot ! — Camil Hamdane, Nov 25 '18 at 15:17

C++ Multinomial distribution

1 Answers1