C ++ Полиномиальное распределение

Question

Я пытаюсь закодировать полиномиальный алгоритм, который будет в основном применять биномиальное распределение к каждому значению входного вектора, зная значения всех предыдущих.Здесь он нацелен на то, чтобы создать новую популяцию для нескольких аллелей , зная начальную популяцию.

Чтобы добиться этого, я использую этот рекурсивный алгоритм:

Вот как выглядит мой код прямо сейчас:

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;  
  }
}

Но я совсем не уверен в этом, и я былинтересно, существовал ли уже такой многочленный алгоритм такого рода ...

Большое спасибо.

Answer 1

Вы можете попробовать использовать команду gsl_ran_multinomial из научной библиотеки GNU.

Функция называется:

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

, где (n_1, n_2, ..., n_K) - неотрицательные целые числас sum_ {k = 1} ^ K n_k = N, а (p_1, p_2, ..., p_K) - это распределение вероятностей с sum(p_i) = 1.Если массив p[K] не нормализован, то его записи будут обрабатываться как веса и соответственно нормализоваться.Массивы n[] и p[] должны иметь длину K.

Функция реализует условный биномиальный метод из книги Дэвиса "Компьютерное генерирование полиномиальных случайных величин" (Comp. Stat. Data Anal16, 1993. link ), так что вы можете реализовать этот подход.Дайте мне знать, если вам нужна копия документа.