Данная статья является продолжением поста Генераторы непрерывно распределенных случайных величин. В этой главе учитывается, что все теоремы из предыдущей статьи уже доказаны и генераторы, указанные в ней, уже реализованы. Как и ранее, у нас имеется некий базовый генератор натуральных чисел от 0 до RAND_MAX:
С дискретными величинами все интуитивно понятнее. Функция распределения дискретной случайной величины:
![](https://habrastorage.org/r/w1560/files/99e/65a/4a5/99e65a4a51ed4235aa4e234029198627.png)
Несмотря на простоту распределений дискретных случайных величин, генерировать их подчас сложнее, нежели чем непрерывные. Начнем, как и в прошлый раз, с тривиального примера.
![](https://habrastorage.org/r/w1560/files/594/995/aef/594995aef57c42a597c5c2770afcf282.png)
![](https://habrastorage.org/r/w1560/files/792/3a5/cb7/7923a5cb7e36429cacaad61739956285.png)
Пожалуй, самый быстрый из очевидных способов сгенерировать случайную величину с распределением Бернулли — это сделать подобным образом (все генераторы возвращают double лишь для единства интерфейса):
Иногда можно сделать быстрее. «Иногда» означает «в случае, при котором параметр p является степенью 1/2». Дело в том, что если целое число, возвращаемое функцией BasicRandGenerator() является равномерно распределенной случайной величиной, то равномерно распределен и каждый его бит. А это значит, что в двоичном представлении число состоит из битов, распределенных по Бернулли. Так как в данных статьях функция базового генератора возвращает unsigned long long, то мы имеем 64 бита. Вот какой трюк можно провернуть для p = 1/2:
Если же время работы функции BasicRandGenerator() недостаточно мало, а криптоустойчивостью и размером периода генератора можно пренебречь, то для таких случаев существует алгоритм, использующий лишь одну равномерно распределенную случайную величину для любого размера выборки из распределения Бернулли:
![](https://habrastorage.org/r/w1560/files/bfb/a14/4df/bfba144dfaee41ffb3548a43afbf87c6.png)
![](https://habrastorage.org/r/w1560/files/424/7d1/6f8/4247d16f823c4a12bd6b16409cf41c30.png)
Уверен, что любой из вас вспомнит, что его первый генератор случайного целого числа от a до b выглядел похожим на это:
Что ж, и это вполне правильное решение. С единственным и не всегда верным предположением — у вас хороший базовый генератор. Если же он дефективный как старый С-шный rand(), то вы будете получать четное число за нечетным каждый раз. Если вы не доверяете своему генератору, то лучше пишите так:
Еще следует заметить, что распределение не будет совсем равномерным, если RAND_MAX не делится на длину интервала b — a + 1 нацело. Однако, разница будет не столь значима, если эта длина достаточно мала.
![](https://habrastorage.org/r/w1560/files/e49/155/542/e4915554219f4d52aa6eb7601f321dda.png)
![](https://habrastorage.org/r/w1560/files/132/0d9/9a3/1320d99a3b0048b2bb89f3ea13ec3dd8.png)
Случайная величина с геометрическим распределением с параметром p — это случайная величина с экспоненциальным распределением с параметром -ln(1 — p), округленная вниз до ближайшего целого.
Можно ли сделать быстрее? Иногда. Посмотрим на функцию распределения:
![](https://habrastorage.org/r/w1560/files/b23/6ac/c0b/b236acc0b05a4583aeb7ff9e4e0843fb.png)
и воспользуемся обыкновенным методом инверсии: генерируем стандартную равномерно распределенную случайную величину U и возвращаем минимальное значение k, для которого эта сумма больше U:
Такой последовательный поиск довольно эффективен, учитывая, что значения случайной величины с геометрическим распределением сконцентрированы возле нуля. Алгоритм, однако, быстро теряет свою скорость при слишком маленьких значениях p, и в таком случае лучше все же воспользоваться экспоненциальным генератором.
Есть секрет. Значения p, (1-p) * p, (1-p)^2 * p,… можно посчитать заранее и запомнить. Вопрос в том, где остановиться. И тут на сцену выходит свойство геометрического распределения, которое оно унаследовало от экспоненциального — отсутствие памяти:
![](https://habrastorage.org/r/w1560/files/744/3a9/b76/7443a9b762ba48439317d5436683a37b.png)
Благодаря такому свойству, можно запомнить лишь несколько первых значений (1-p)^i * p и затем воспользоваться рекурсией:
![](https://habrastorage.org/r/w1560/files/261/815/4b6/2618154b61574a3097c1e904b04004a8.png)
![](https://habrastorage.org/r/w1560/files/362/a20/9be/362a209be1a24fa384af9b5c4c833759.png)
По определению случайная величина с биномиальным распределением — это сумма n случайных величин с распределением Бернулли:
Однако, есть линейная зависимость от n, которую нужно обойти. Можно воспользоваться теоремой: если Y_1, Y_2,… — случайные величины с геометрическим распределением с параметром p и X — наименьшее целое, такое что:
![](https://habrastorage.org/r/w1560/files/d26/c9d/6d3/d26c9d6d356a499d8d2a7d3a2aee389a.png)
то X имеет биномиальное распределение.
Время работы следующего кода растет только с увеличением величины n * p.
И все же, временная сложность растет.
У биномиального распределения есть одна особенность. При росте n оно стремится к нормальному распределению или же, если p ~ 1/n, то к распределению Пуассона. Имея генераторы для этих распределений, можно ими заменить генератор для биномиального в подобных случаях. Но что, если этого недостаточно? В книге Luc Devroye «Non-Uniform Random Variate Generation» есть пример алгоритма, работающего одинаково быстро для любых больших значений n * p. Идея алгоритма состоит в выборке с отклонением, используя нормальное и экспоненциальное распределения. К сожалению, рассказ о работе этого алгоритма будет слишком большим для данной статьи, но в указанной книге он исчерпывающе описан.
![](https://habrastorage.org/r/w1560/files/e17/cba/d3e/e17cbad3e4314e10983c3907ce8d15bc.png)
![](https://habrastorage.org/r/w1560/files/7b5/657/993/7b56579932ae4300bef692441cde30fd.png)
Если W_i — случайная величина со стандартным экспоненциальным распределением с плотностью lambda, то:
![](https://habrastorage.org/r/w1560/files/1a9/96b/c0e/1a996bc0e7224352a5b5ffe99b684347.png)
Используя это свойство, можно написать генератор через сумму экспоненциально распределенных величин с плотностью rate:
На этом же свойстве основан популярный алгоритм Кнута. Вместо суммы экспоненциальных величин, каждую из которых можно получить методом инверсии через -ln(U), используется произведение равномерных случайных величин:
![](https://habrastorage.org/r/w1560/files/1e9/692/637/1e96926372c143c0aeaf39b155fce5c0.png)
и тогда:
![](https://habrastorage.org/r/w1560/files/9a6/b83/d5a/9a6b83d5a9994388a5595b7cd1ea4c39.png)
Запомнив заранее значение exp(-rate), последовательно перемножаем U_i, пока произведение его не превысит.
Можно же использовать генерацию только лишь одной случайной величины и метод инверсии:
![](https://habrastorage.org/r/w1560/files/23a/6a6/2e5/23a6a62e57f64899ac89aa894f01abe6.png)
Поиск k, удовлетворяющего такому условию лучше начинать с наиболее вероятного значения, то есть с floor(rate). Сравниваем U с вероятностью, что X < floor(rate) и идем вверх, если U больше, или вниз в другом случае:
Проблема всех трех генераторов в одном — их сложность растет вместе с параметром плотности. Но есть спасение — при больших значениях плотности распределение Пуассона стремится к нормальному распределению. Можно также использовать непростой алгоритм из вышеуказанной книги «Non-Uniform Random Variate Generation», ну, или же, попросту аппроксимировать, пренебрегая точностью во имя скорости (смотря какая стоит задача).
![](https://habrastorage.org/r/w1560/files/0f6/4cf/fcd/0f64cffcd84941a8aa0f04b62b054a5f.png)
![](https://habrastorage.org/r/w1560/files/e4e/efb/eb0/e4eefbeb0bf9469fa4d4e4337ad8c1af.png)
Отрицательное биномиальное распределение еще именуется распределением Паскаля (если r целое) и распределением Поля (если r может быть вещественным). Используя характеристическую функцию, легко доказать, что распределение Паскаля — это сумма r геометрически распределенных величин с параметром 1 — p:
Проблема такого алгоритма налицо — линейная зависимость от r. Нам нужно что-то такое, что будет работать одинаково хорошо при любом параметре. И в этом нам поможет отличное свойство. Если случайная величина X имеет распределение Пуассона:
![](https://habrastorage.org/r/w1560/files/034/dc5/ad9/034dc5ad968c4e5da9f47de05bb0733d.png)
где плотность случайна и имеет гамма-распределение:
![](https://habrastorage.org/r/w1560/files/701/434/8af/7014348af9754038a57c4803e33aced8.png)
то X имеет отрицательное биномиальное распределение. Поэтому, оно иногда называется гамма-распределением Пуассона.
Теперь, мы можем быстро написать генератор случайной величины с распределением Паскаля для больших значений параметра r.
![](https://habrastorage.org/r/w1560/files/83a/154/d0f/83a154d0f6c8490e8cab2e895e28645d.png)
![](https://habrastorage.org/r/w1560/files/559/598/597/559598597ca64dbaa54e538478634c7b.png)
Представьте, что в урне находится N шаров и K из них белые. Вы вытаскиваете n шаров. Количество белых среди них будет иметь гипергеометрическое распределение. В общем случае лучше брать алгоритм, использующий это определение:
Или же можно воспользоваться выборкой с отклонением через биномиальное распределение с параметрами:
![](https://habrastorage.org/r/w1560/files/abb/f41/a17/abbf41a175ac4dfba28858c1a0d77408.png)
и константой М:
![](https://habrastorage.org/r/w1560/files/ebf/6bc/472/ebf6bc472d3941e993a2932d31280c9c.png)
Алгоритм с биномиальным распределением хорошо работает в экстремальных случаях при больших n и еще больших N (таких, что n << N). В остальном — лучше использовать предыдущий, несмотря на его растущую сложность с увеличением входящих параметров.
![](https://habrastorage.org/r/w1560/files/24c/8b3/ee3/24c8b3ee3215449e8d480b619e8709b8.png)
![](https://habrastorage.org/r/w1560/files/723/d1f/72b/723d1f72bc93492fa87fc03928a4dc04.png)
Для данного распределения существует полезная теорема. Если U и V — равномерно распределенные случайные величины от 0 до 1, то величина
![](https://habrastorage.org/r/w1560/files/d86/562/ae1/d86562ae14504011a577f64b6f286a3e.png)
будет иметь логарифмическое распределение.
Для ускорения алгоритма можно воспользоваться двумя уловками. Первая: если V > p, то X = 1, т.к. p >= 1-(1-p)^U. Вторая: пускай q = 1-(1-p)^U, тогда если V > q, то X = 1, если V > q^2, то X = 2 и т.д. Таким образом, можно возвращать наиболее вероятные значения 1 и 2 без частых подсчетов логарифмов.
![](https://habrastorage.org/r/w1560/files/be4/534/640/be4534640928488294915f54a3460bac.png)
В знаменателе — зета-функция Римана:
![](https://habrastorage.org/r/w1560/files/89f/ac3/841/89fac38416a343dca149dc83388cc7eb.png)
![](https://habrastorage.org/r/w1560/files/39a/17d/4bb/39a17d4bb3894bb8b27c18f4c8e548b0.png)
Для зета распределения существует алгоритм, позволяющий не вычислять римановскую зета-функцию. Нужно лишь уметь генерировать распределение Парето. Доказательство могу приложить по просьбе читателей.
Напоследок маленькие алгоритмы для прочих сложных распределений:
Реализация этих и других генераторов на C++17
unsigned long long BasicRandGenerator() {
unsigned long long randomVariable;
// some magic here
...
return randomVariable;
}
С дискретными величинами все интуитивно понятнее. Функция распределения дискретной случайной величины:
![](https://habrastorage.org/files/99e/65a/4a5/99e65a4a51ed4235aa4e234029198627.png)
Несмотря на простоту распределений дискретных случайных величин, генерировать их подчас сложнее, нежели чем непрерывные. Начнем, как и в прошлый раз, с тривиального примера.
Распределение Бернулли
![](https://habrastorage.org/files/594/995/aef/594995aef57c42a597c5c2770afcf282.png)
![](https://habrastorage.org/files/792/3a5/cb7/7923a5cb7e36429cacaad61739956285.png)
Пожалуй, самый быстрый из очевидных способов сгенерировать случайную величину с распределением Бернулли — это сделать подобным образом (все генераторы возвращают double лишь для единства интерфейса):
void setup(double p) {
unsigned long long boundary = (1 - p) * RAND_MAX; // we storage this value for large sampling
}
double Bernoulli(double) {
return BasicRandGenerator() > boundary;
}
Иногда можно сделать быстрее. «Иногда» означает «в случае, при котором параметр p является степенью 1/2». Дело в том, что если целое число, возвращаемое функцией BasicRandGenerator() является равномерно распределенной случайной величиной, то равномерно распределен и каждый его бит. А это значит, что в двоичном представлении число состоит из битов, распределенных по Бернулли. Так как в данных статьях функция базового генератора возвращает unsigned long long, то мы имеем 64 бита. Вот какой трюк можно провернуть для p = 1/2:
double Bernoulli(double) {
static const char maxDecimals = 64;
static char decimals = 1;
static unsigned long long X = 0;
if (decimals == 1)
{
/// refresh
decimals = maxDecimals;
X = BasicRandGenerator();
}
else
{
--decimals;
X >>= 1;
}
return X & 1;
}
Если же время работы функции BasicRandGenerator() недостаточно мало, а криптоустойчивостью и размером периода генератора можно пренебречь, то для таких случаев существует алгоритм, использующий лишь одну равномерно распределенную случайную величину для любого размера выборки из распределения Бернулли:
void setup(double p) {
q = 1.0 - p;
U = Uniform(0, 1);
}
double Bernoulli(double p) {
if (U > q)
{
U -= q;
U /= p;
return 1.0;
}
U /= q;
return 0.0;
}
Почему этот алгоритм работает?
![](https://habrastorage.org/r/w1560/files/b58/624/591/b5862459180b4fd197b65e5b78d7f2b3.png)
![](https://habrastorage.org/r/w1560/files/f96/acd/832/f96acd832ac2415db05edab4eef37014.png)
![](https://habrastorage.org/files/b58/624/591/b5862459180b4fd197b65e5b78d7f2b3.png)
![](https://habrastorage.org/files/f96/acd/832/f96acd832ac2415db05edab4eef37014.png)
Равномерное распределение
![](https://habrastorage.org/files/bfb/a14/4df/bfba144dfaee41ffb3548a43afbf87c6.png)
![](https://habrastorage.org/files/424/7d1/6f8/4247d16f823c4a12bd6b16409cf41c30.png)
Уверен, что любой из вас вспомнит, что его первый генератор случайного целого числа от a до b выглядел похожим на это:
double UniformDiscrete(int a, int b) {
return a + rand() % (b - a + 1);
}
Что ж, и это вполне правильное решение. С единственным и не всегда верным предположением — у вас хороший базовый генератор. Если же он дефективный как старый С-шный rand(), то вы будете получать четное число за нечетным каждый раз. Если вы не доверяете своему генератору, то лучше пишите так:
double UniformDiscrete(int a, int b) {
return a + round(BasicRandGenerator() * (b - a) / RAND_MAX);
}
Еще следует заметить, что распределение не будет совсем равномерным, если RAND_MAX не делится на длину интервала b — a + 1 нацело. Однако, разница будет не столь значима, если эта длина достаточно мала.
Геометрическое распределение
![](https://habrastorage.org/files/e49/155/542/e4915554219f4d52aa6eb7601f321dda.png)
![](https://habrastorage.org/files/132/0d9/9a3/1320d99a3b0048b2bb89f3ea13ec3dd8.png)
Случайная величина с геометрическим распределением с параметром p — это случайная величина с экспоненциальным распределением с параметром -ln(1 — p), округленная вниз до ближайшего целого.
Доказательство
Пускай W — случайная величина, распределенная экспоненциально с параметром -ln(1 — p). Тогда:
![](https://habrastorage.org/r/w1560/files/4dd/ce6/96e/4ddce696e19d485d8301a89310b28a82.png)
![](https://habrastorage.org/files/4dd/ce6/96e/4ddce696e19d485d8301a89310b28a82.png)
void setupRate(double p) {
rate = -ln(1 - p);
}
double Geometric(double) {
return floor(Exponential(rate));
}
Можно ли сделать быстрее? Иногда. Посмотрим на функцию распределения:
![](https://habrastorage.org/files/b23/6ac/c0b/b236acc0b05a4583aeb7ff9e4e0843fb.png)
и воспользуемся обыкновенным методом инверсии: генерируем стандартную равномерно распределенную случайную величину U и возвращаем минимальное значение k, для которого эта сумма больше U:
double GeometricExponential(double p) {
int k = 0;
double sum = p, prod = p, q = 1 - p;
double U = Uniform(0, 1);
while (U < sum) {
prod *= q;
sum += prod;
++k;
}
return k;
}
Такой последовательный поиск довольно эффективен, учитывая, что значения случайной величины с геометрическим распределением сконцентрированы возле нуля. Алгоритм, однако, быстро теряет свою скорость при слишком маленьких значениях p, и в таком случае лучше все же воспользоваться экспоненциальным генератором.
Есть секрет. Значения p, (1-p) * p, (1-p)^2 * p,… можно посчитать заранее и запомнить. Вопрос в том, где остановиться. И тут на сцену выходит свойство геометрического распределения, которое оно унаследовало от экспоненциального — отсутствие памяти:
![](https://habrastorage.org/files/744/3a9/b76/7443a9b762ba48439317d5436683a37b.png)
Благодаря такому свойству, можно запомнить лишь несколько первых значений (1-p)^i * p и затем воспользоваться рекурсией:
// works nice for p > 0.2 and tableSize = 16
void setupTable(double p) {
table[0] = p;
double prod = p, q = 1 - p;
for (int i = 1; i < tableSize; ++i)
{
prod *= q;
table[i] = table[i - 1] + prod;
}
}
double GeometricTable(double p) {
double U = Uniform(0, 1);
/// handle tail by recursion
if (U > table[tableSize - 1])
return tableSize + GeometricTable(p);
/// handle the main body
int x = 0;
while (U > table[x])
++x;
return x;
}
Биномиальное распределение
![](https://habrastorage.org/files/261/815/4b6/2618154b61574a3097c1e904b04004a8.png)
![](https://habrastorage.org/files/362/a20/9be/362a209be1a24fa384af9b5c4c833759.png)
По определению случайная величина с биномиальным распределением — это сумма n случайных величин с распределением Бернулли:
double BinomialBernoulli(double p, int n) {
double sum = 0;
for (int i = 0; i != n; ++i)
sum += Bernoulli(p);
return sum;
}
Однако, есть линейная зависимость от n, которую нужно обойти. Можно воспользоваться теоремой: если Y_1, Y_2,… — случайные величины с геометрическим распределением с параметром p и X — наименьшее целое, такое что:
![](https://habrastorage.org/files/d26/c9d/6d3/d26c9d6d356a499d8d2a7d3a2aee389a.png)
то X имеет биномиальное распределение.
Доказательство
По определению, случайная величина с геометрическим распределением Y — это количество опытов Бернулли до первого успеха. Есть альтернативное определение: Z = Y + 1 — количество опытов Бернулли до первого успеха включительно. Значит, сумма независимых Z_i — количество опытов Бернулли до X + 1 успеха включительно. И эта сумма больше n тогда и только тогда, когда среди первых n испытаний не более чем X успешных. Соответственно:
![](https://habrastorage.org/r/w1560/files/842/a22/0f0/842a220f063049bcba3cc3b3071e6c86.png)
q.e.d.
![](https://habrastorage.org/files/842/a22/0f0/842a220f063049bcba3cc3b3071e6c86.png)
q.e.d.
Время работы следующего кода растет только с увеличением величины n * p.
double BinomialGeometric(double p, int n) {
double X = -1, sum = 0;
do {
sum += Geometric(p) + 1.0;
++X;
} while (sum <= n);
return X;
}
И все же, временная сложность растет.
У биномиального распределения есть одна особенность. При росте n оно стремится к нормальному распределению или же, если p ~ 1/n, то к распределению Пуассона. Имея генераторы для этих распределений, можно ими заменить генератор для биномиального в подобных случаях. Но что, если этого недостаточно? В книге Luc Devroye «Non-Uniform Random Variate Generation» есть пример алгоритма, работающего одинаково быстро для любых больших значений n * p. Идея алгоритма состоит в выборке с отклонением, используя нормальное и экспоненциальное распределения. К сожалению, рассказ о работе этого алгоритма будет слишком большим для данной статьи, но в указанной книге он исчерпывающе описан.
Распределение Пуассона
![](https://habrastorage.org/files/e17/cba/d3e/e17cbad3e4314e10983c3907ce8d15bc.png)
![](https://habrastorage.org/files/7b5/657/993/7b56579932ae4300bef692441cde30fd.png)
Если W_i — случайная величина со стандартным экспоненциальным распределением с плотностью lambda, то:
![](https://habrastorage.org/files/1a9/96b/c0e/1a996bc0e7224352a5b5ffe99b684347.png)
Доказательство
Пускай f_k(y) — плотность стандартного распределения Эрланга (суммы k стандартных экспоненциально распределенных случайных величин):
![](https://habrastorage.org/r/w1560/files/7ae/6fd/f29/7ae6fdf29edf4ac69bbf41f0f10f8690.png)
тогда
![](https://habrastorage.org/r/w1560/files/907/a99/8b9/907a998b9d874cea9a3d91bff8b66065.png)
и вероятность получить k:
![](https://habrastorage.org/r/w1560/files/2a7/4e9/c73/2a74e9c73f894fc0b34290921151efeb.png)
совпадает с распределением Пуассона. q.e.d.
![](https://habrastorage.org/files/7ae/6fd/f29/7ae6fdf29edf4ac69bbf41f0f10f8690.png)
тогда
![](https://habrastorage.org/files/907/a99/8b9/907a998b9d874cea9a3d91bff8b66065.png)
и вероятность получить k:
![](https://habrastorage.org/files/2a7/4e9/c73/2a74e9c73f894fc0b34290921151efeb.png)
совпадает с распределением Пуассона. q.e.d.
Используя это свойство, можно написать генератор через сумму экспоненциально распределенных величин с плотностью rate:
double PoissonExponential(double rate) {
int k = -1;
double s = 0;
do {
s += Exponential(1);
++k;
} while (s < rate);
return k;
}
На этом же свойстве основан популярный алгоритм Кнута. Вместо суммы экспоненциальных величин, каждую из которых можно получить методом инверсии через -ln(U), используется произведение равномерных случайных величин:
![](https://habrastorage.org/files/1e9/692/637/1e96926372c143c0aeaf39b155fce5c0.png)
и тогда:
![](https://habrastorage.org/files/9a6/b83/d5a/9a6b83d5a9994388a5595b7cd1ea4c39.png)
Запомнив заранее значение exp(-rate), последовательно перемножаем U_i, пока произведение его не превысит.
void setup(double rate)
expRateInv = exp(-rate);
}
double PoissonKnuth(double) {
double k = 0, prod = Uniform(0, 1);
while (prod > expRateInv) {
prod *= Uniform(0, 1);
++k;
}
return k;
}
Можно же использовать генерацию только лишь одной случайной величины и метод инверсии:
![](https://habrastorage.org/files/23a/6a6/2e5/23a6a62e57f64899ac89aa894f01abe6.png)
Поиск k, удовлетворяющего такому условию лучше начинать с наиболее вероятного значения, то есть с floor(rate). Сравниваем U с вероятностью, что X < floor(rate) и идем вверх, если U больше, или вниз в другом случае:
void setup(double rate)
floorLambda = floor(rate);
FFloorLambda = F(floorLamda); // P(X < floorLambda)
PFloorLambda = P(floorLambda); // P(X = floorLambda)
}
double PoissonInversion(double) {
double U = Uniform(0,1);
int k = floorLambda;
double s = FFloorLambda, p = PFloorLambda;
if (s < U) {
do {
++k;
p *= lambda / k;
s += p;
} while (s < U);
}
else {
s -= p;
while (s > U) {
p /= lambda / k;
--k;
s -= p;
}
}
return k;
}
Проблема всех трех генераторов в одном — их сложность растет вместе с параметром плотности. Но есть спасение — при больших значениях плотности распределение Пуассона стремится к нормальному распределению. Можно также использовать непростой алгоритм из вышеуказанной книги «Non-Uniform Random Variate Generation», ну, или же, попросту аппроксимировать, пренебрегая точностью во имя скорости (смотря какая стоит задача).
Отрицательное биномиальное распределение
![](https://habrastorage.org/files/0f6/4cf/fcd/0f64cffcd84941a8aa0f04b62b054a5f.png)
![](https://habrastorage.org/files/e4e/efb/eb0/e4eefbeb0bf9469fa4d4e4337ad8c1af.png)
Отрицательное биномиальное распределение еще именуется распределением Паскаля (если r целое) и распределением Поля (если r может быть вещественным). Используя характеристическую функцию, легко доказать, что распределение Паскаля — это сумма r геометрически распределенных величин с параметром 1 — p:
double Pascal(double p, int r) {
double x = 0;
for (int i = 0; i != r; ++i)
x += Geometric(1 - p);
return x;
}
Проблема такого алгоритма налицо — линейная зависимость от r. Нам нужно что-то такое, что будет работать одинаково хорошо при любом параметре. И в этом нам поможет отличное свойство. Если случайная величина X имеет распределение Пуассона:
![](https://habrastorage.org/files/034/dc5/ad9/034dc5ad968c4e5da9f47de05bb0733d.png)
где плотность случайна и имеет гамма-распределение:
![](https://habrastorage.org/files/701/434/8af/7014348af9754038a57c4803e33aced8.png)
то X имеет отрицательное биномиальное распределение. Поэтому, оно иногда называется гамма-распределением Пуассона.
Доказательство
![](https://habrastorage.org/r/w1560/files/758/e85/9a9/758e859a9b3d48a2a12c88689b4825c5.png)
![](https://habrastorage.org/files/758/e85/9a9/758e859a9b3d48a2a12c88689b4825c5.png)
Теперь, мы можем быстро написать генератор случайной величины с распределением Паскаля для больших значений параметра r.
double Pascal(double p, int r) {
return Poisson(Gamma(r, p / (1 - p)));
}
Гипергеометрическое распределение
![](https://habrastorage.org/files/83a/154/d0f/83a154d0f6c8490e8cab2e895e28645d.png)
![](https://habrastorage.org/files/559/598/597/559598597ca64dbaa54e538478634c7b.png)
Представьте, что в урне находится N шаров и K из них белые. Вы вытаскиваете n шаров. Количество белых среди них будет иметь гипергеометрическое распределение. В общем случае лучше брать алгоритм, использующий это определение:
double HyperGeometric(int N, int n, int K) {
double sum = 0, p = static_cast<double>(K) / N;
for (int i = 1; i <= n; ++i)
{
if (Bernoulli(p) && ++sum == K)
return sum;
p = (K - sum) / (N - i);
}
return sum;
}
Или же можно воспользоваться выборкой с отклонением через биномиальное распределение с параметрами:
![](https://habrastorage.org/files/abb/f41/a17/abbf41a175ac4dfba28858c1a0d77408.png)
и константой М:
![](https://habrastorage.org/files/ebf/6bc/472/ebf6bc472d3941e993a2932d31280c9c.png)
Алгоритм с биномиальным распределением хорошо работает в экстремальных случаях при больших n и еще больших N (таких, что n << N). В остальном — лучше использовать предыдущий, несмотря на его растущую сложность с увеличением входящих параметров.
Логарифмическое распределение
![](https://habrastorage.org/files/24c/8b3/ee3/24c8b3ee3215449e8d480b619e8709b8.png)
![](https://habrastorage.org/files/723/d1f/72b/723d1f72bc93492fa87fc03928a4dc04.png)
Для данного распределения существует полезная теорема. Если U и V — равномерно распределенные случайные величины от 0 до 1, то величина
![](https://habrastorage.org/files/d86/562/ae1/d86562ae14504011a577f64b6f286a3e.png)
будет иметь логарифмическое распределение.
Доказательство
Вспомнив генераторы для экспоненциального и геометрического распределений, несложно увидеть, что X, заданный формулой выше, удовлетворяет следующему распределению:
![](https://habrastorage.org/r/w1560/files/e8b/9b6/e99/e8b9b6e99b53427897627aed86ff540c.png)
Докажем, что X имеет логарифмическое распределение:
![](https://habrastorage.org/r/w1560/files/482/8aa/712/4828aa7122df4de695d7c8e37e162652.png)
q.e.d.
![](https://habrastorage.org/files/e8b/9b6/e99/e8b9b6e99b53427897627aed86ff540c.png)
Докажем, что X имеет логарифмическое распределение:
![](https://habrastorage.org/files/482/8aa/712/4828aa7122df4de695d7c8e37e162652.png)
q.e.d.
Для ускорения алгоритма можно воспользоваться двумя уловками. Первая: если V > p, то X = 1, т.к. p >= 1-(1-p)^U. Вторая: пускай q = 1-(1-p)^U, тогда если V > q, то X = 1, если V > q^2, то X = 2 и т.д. Таким образом, можно возвращать наиболее вероятные значения 1 и 2 без частых подсчетов логарифмов.
void setup(double p) {
logQInv = 1.0 / log(1.0 - p);
}
double Logarithmic(double p) {
double V = Uniform(0, 1);
if (V >= p)
return 1.0;
double U = Uniform(0, 1);
double y = 1.0 - exp(U / logQInv);
if (V > y)
return 1.0;
if (V <= y * y)
return floor(1.0 + log(V) / log(y));
return 2.0;
}
Зета распределение
![](https://habrastorage.org/files/be4/534/640/be4534640928488294915f54a3460bac.png)
В знаменателе — зета-функция Римана:
![](https://habrastorage.org/files/89f/ac3/841/89fac38416a343dca149dc83388cc7eb.png)
![](https://habrastorage.org/files/39a/17d/4bb/39a17d4bb3894bb8b27c18f4c8e548b0.png)
Для зета распределения существует алгоритм, позволяющий не вычислять римановскую зета-функцию. Нужно лишь уметь генерировать распределение Парето. Доказательство могу приложить по просьбе читателей.
void setup(double s) {
double sm1 = s - 1.0;
b = 1 - pow(2.0, -sm1);
}
double Zeta(double) {
/// rejection sampling from rounded down Pareto distribution
int iter = 0;
do {
double X = floor(Pareto(sm1, 1.0));
double T = pow(1.0 + 1.0 / X, sm1);
double V = Uniform(0, 1);
if (V * X * (T - 1) <= T * b)
return X;
} while (++iter < 1e9);
return NAN; /// doesn't work
}
Напоследок маленькие алгоритмы для прочих сложных распределений:
double Skellam(double m1, double m2) {
return Poisson(m1) - Poisson(m2);
}
double Planck(double a, double b) {
double G = Gamma(a + 1, 1);
double Z = Zeta(a + 1)
return G / (b * Z);
}
double Yule(double ro) {
double prob = 1.0 / Pareto(ro, 1.0);
return Geometric(prob) + 1;
}
Реализация этих и других генераторов на C++17