Создание случайных величин, распределенных гипергеометрически - PullRequest
Новая
       30

Создание случайных величин, распределенных гипергеометрически

0 голосов
/ 07 ноября 2018

Как мы можем генерировать случайные переменные для гипергеометрического распределения в C ++? К сожалению, стандартная библиотека не предоставляет std::hypergeometric_distribution. Есть boost::math::hypergeometric_distribution; однако, похоже, что его нельзя использовать для генерации случайных переменных.

1 Ответ

0 голосов
/ 07 ноября 2018

Ну, вот простой код (возможно, с ошибками, хотя он компилируется, запускается и создает нечто, похожее на правду) - выборка Гипергеометрическое распределение . Обозначения такие же, как в статье Wiki. Внутренне, он заполняет массив значениями [0 ... N-1], перетасовывает первые n из них и считает, если перетасованные значения меньше K (успехов).

Shuffle классический Фишер-Йейтс-Кнут один. 

#include <cstdint>
#include <iostream>
#include <random>
#include <vector>
#include <numeric>
#include <cmath>

#define func auto


func nCk(uint64_t n, uint64_t k) -> double { // computes binomial coefficient
    if (k > n)
        return 0.0;

    if (k * 2ULL > n)
        k = n - k;

    if (k == 0.0)
        return 1.0;

    return exp(lgamma(n+1) - (lgamma(k+1) + lgamma(n-k+1)));
}


func PMF(uint64_t N, uint64_t K, uint64_t n, uint64_t k) -> double { // Hypergeometric distribution PMF
    return nCk(K, k) * nCk(N - K, n - k) / nCk(N, n);
}


func sample_hyper(int N, int K, int n,            // sampling from Hypergeometric distribution
                  std::mt19937_64&  rng,
                  std::vector<int>& nums) -> int {

    int rc{ 0 };
    for (int k = 0; k != n; ++k) {
        std::uniform_int_distribution<int> uni{ k, N-1 };
        auto s = uni(rng);
        std::swap(nums[k], nums[s]);

        rc += (nums[k] < K);
    }
    return rc;
}


func main() -> int {
    auto rng = std::mt19937_64{1234567891ULL};

    int N = 500; // compare with Wiki
    int K = 50;
    int n = 100;

    auto nums = std::vector<int>( N );
    std::iota(nums.begin(), nums.end(), 0); // array to shuffle, filled with 0, 1, 2, 3 ... sequence

    auto h = std::vector<int>( K, 0 ); // histogram

    int NT = 100000; // number of trials
    for(int k = 0; k != NT; ++k) {
        int q = sample_hyper(N, K, n, rng, nums);
        h[q] += 1; // count it
    }

    for (int k = 0; k != 20; ++k) { // and print it
        double e = double(h[k]) / double(NT);
        std::cout << k << "   " << e << "     " << PMF(N, K, n, k) << '\n';
    }

    return 0;
}

UPDATE

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

...