Информационный портал по безопасности » Программирование » Генераторы дискретно распределенных случайных величин

 

Генераторы дискретно распределенных случайных величин

Автор: admin от 16-01-2016, 01:54, посмотрело: 1 298

Данная статья является продолжением поста Генераторы непрерывно распределенных случайных величин. В этой главе учитывается, что все теоремы из предыдущей статьи уже доказаны и генераторы, указанные в ней, уже написаны. Как и ранее, у нас имеется некий базовый генератор натуральных чисел от 0 до RAND_MAX:

unsigned long long BasicRandGenerator() {
    unsigned long long randomVariable;
    // some magic here
    ...
    return randomVariable;
}

С дискретными величинами все интуитивно понятнее. Функция распределения дискретной случайной величины:

Генераторы дискретно распределенных случайных величин

Несмотря на простоту распределений дискретных случайных величин, генерировать их подчас сложнее, нежели чем непрерывные. Начнем, как и в прошлый раз, с тривиального примера.

Распределение Бернулли


Генераторы дискретно распределенных случайных величин

Генераторы дискретно распределенных случайных величин

Пожалуй, самый быстрый из очевидных способов сгенерировать случайную величину с распределением Бернулли — это сделать подобным образом (все генераторы возвращают 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;
}



Равномерное распределение


Генераторы дискретно распределенных случайных величин

Генераторы дискретно распределенных случайных величин

Уверен, что любой из вас вспомнит, что его первый генератор случайного целого числа от 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 нацело. Однако, разница будет не столь значима, если эта длина достаточно мала.

Геометрическое распределение


Генераторы дискретно распределенных случайных величин

Генераторы дискретно распределенных случайных величин

Случайная величина с геометрическим распределением с параметром p — это случайная величина с экспоненциальным распределением с параметром -ln(1 — p), округленная вниз до ближайшего целого.


void setupRate(double p) {
    rate = -ln(1 - p);
}

double Geometric(double) {
    return floor(Exponential(rate));
}

Можно ли сделать быстрее? Иногда. Посмотрим на функцию распределения:

Генераторы дискретно распределенных случайных величин

и воспользуемся обыкновенным методом инверсии: генерируем стандартную равномерно распределенную случайную величину 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,… можно посчитать заранее и запомнить. Вопрос в том, где остановиться. И тут на сцену выходит свойство геометрического распределения, которое оно унаследовало от экспоненциального — отсутствие памяти:

Генераторы дискретно распределенных случайных величин

Благодаря такому свойству, можно запомнить лишь несколько первых значений (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;
}


Биномиальное распределение


Генераторы дискретно распределенных случайных величин

Генераторы дискретно распределенных случайных величин

По определению случайная величина с биномиальным распределением — это сумма 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 — наименьшее целое, такое что:

Генераторы дискретно распределенных случайных величин

то X имеет биномиальное распределение.

Время работы следующего кода растет только с увеличением величины 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. Идея алгоритма состоит в выборке с отклонением, используя нормальное и экспоненциальное распределения. К сожалению, рассказ о работе этого алгоритма будет слишком большим для данной статьи, но в указанной книге он исчерпывающе описан.

Распределение Пуассона


Генераторы дискретно распределенных случайных величин

Генераторы дискретно распределенных случайных величин

Если W_i — случайная величина со стандартным экспоненциальным распределением с плотностью lambda, то:

Генераторы дискретно распределенных случайных величин


Используя это свойство, можно написать генератор через сумму экспоненциально распределенных величин с плотностью rate:

double PoissonExponential(double rate) {
    int k = -1;
    double s = 0;
    do {
        s += Exponential(1);
        ++k;
    } while (s < rate);
    return k;
}

На этом же свойстве основан популярный алгоритм Кнута. Вместо суммы экспоненциальных величин, каждую из которых можно получить методом инверсии через -ln(U), используется произведение равномерных случайных величин:

Генераторы дискретно распределенных случайных величин

и тогда:

Генераторы дискретно распределенных случайных величин

Запомнив заранее значение 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;
}

Можно же использовать генерацию только лишь одной случайной величины и метод инверсии:

Генераторы дискретно распределенных случайных величин

Поиск 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», ну, или же, попросту аппроксимировать, пренебрегая точностью во имя скорости (смотря какая стоит задача).

Отрицательное биномиальное распределение


Генераторы дискретно распределенных случайных величин

Генераторы дискретно распределенных случайных величин

Отрицательное биномиальное распределение еще именуется распределением Паскаля (если 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 имеет распределение Пуассона:

Генераторы дискретно распределенных случайных величин

где плотность случайна и имеет гамма-распределение:

Генераторы дискретно распределенных случайных величин

то X имеет отрицательное биномиальное распределение. Поэтому, оно иногда называется гамма-распределением Пуассона.


Теперь, мы можем быстро написать генератор случайной величины с распределением Паскаля для больших значений параметра r.

double Pascal(double p, int r) {
    return Poisson(Gamma(r, p / (1 - p)));
}


Гипергеометрическое распределение


Генераторы дискретно распределенных случайных величин

Генераторы дискретно распределенных случайных величин

Представьте, что в урне находится N шаров и K из них белые. Вы вытаскиваете n шаров. Количество белых среди них будет иметь гипергеометрическое распределение. В общем случае лучше брать алгоритм, использующий это определение:

double HyperGeometric(int N, int n, int K) {
    int sum = 0;
    double realK = static_cast<double>(K);
    double p = realK / N;
    for (int i = 1; i <= n; ++i)
    {
        if (Bernoulli(p) && ++sum == K)
            return sum;
        p = (realK - sum) / (N - i);
    }
    return sum;
}

Или же можно воспользоваться выборкой с отклонением через биномиальное распределение с параметрами:

Генераторы дискретно распределенных случайных величин

и константой М:

Генераторы дискретно распределенных случайных величин

Алгоритм с биномиальным распределением хорошо работает в экстремальных случаях при больших n и еще больших N (таких, что n N). В остальном — лучше использовать предыдущий, несмотря на его растущую сложность с увеличением входящих параметров.

Логарифмическое распределение


Генераторы дискретно распределенных случайных величин

Генераторы дискретно распределенных случайных величин

Для данного распределения существует полезная теорема. Если U и V — равномерно распределенные случайные величины от 0 до 1, то величина

Генераторы дискретно распределенных случайных величин

будет иметь логарифмическое распределение.


Для ускорения алгоритма можно воспользоваться двумя уловками. Первая: если 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;
}


Зета распределение


Генераторы дискретно распределенных случайных величин

В знаменателе — зета-функция Римана:
Генераторы дискретно распределенных случайных величин

Генераторы дискретно распределенных случайных величин

Для зета распределения существует алгоритм, позволяющий не вычислять римановскую зета-функцию. Нужно лишь уметь генерировать распределение Парето. Доказательство могу приложить по просьбе читателей.

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


Источник: Хабрахабр

Категория: Программирование

Уважаемый посетитель, Вы зашли на сайт как незарегистрированный пользователь.
Мы рекомендуем Вам зарегистрироваться либо войти на сайт под своим именем.

Добавление комментария

Имя:*
E-Mail:
Комментарий:
  • bowtiesmilelaughingblushsmileyrelaxedsmirk
    heart_eyeskissing_heartkissing_closed_eyesflushedrelievedsatisfiedgrin
    winkstuck_out_tongue_winking_eyestuck_out_tongue_closed_eyesgrinningkissingstuck_out_tonguesleeping
    worriedfrowninganguishedopen_mouthgrimacingconfusedhushed
    expressionlessunamusedsweat_smilesweatdisappointed_relievedwearypensive
    disappointedconfoundedfearfulcold_sweatperseverecrysob
    joyastonishedscreamtired_faceangryragetriumph
    sleepyyummasksunglassesdizzy_faceimpsmiling_imp
    neutral_faceno_mouthinnocent