#include <boost/math/distributions/hypergeometric.hpp>

namespace boost{ namespace math{
template <class RealType = double,
          class Policy   = policies::policy<> >
class hypergeometric_distribution;
template <class RealType, class Policy>
class hypergeometric_distribution
{
public:
   typedef RealType value_type;
   typedef Policy   policy_type;
   // Construct:
   hypergeometric_distribution(unsigned r, unsigned n, unsigned N);
   // Accessors:
   unsigned total()const;
   unsigned defective()const;
   unsigned sample_count()const;
};
typedef hypergeometric_distribution<> hypergeometric;
}} // namespaces

Гипергеометрическое распределение описывает количество «событий» k из выборки n, взятой из общей популяции N без замены .

Представьте, что у нас есть выборка N объектов, из которых r являются «дефектными», а N-r — «не дефектными» (также используются термины «успех/провал» или «красный/синий»). Если мы отобрали n элементы без замены , то какова вероятность того, что именно k элементы в образце являются дефектными? Ответ дается pdf гипергеометрического распределения f(k; n, N), в то время как вероятность k, где F(k) является CDF гипергеометрического распределения.

	Note
	В отличие от почти всех других распределений в этой библиотеке, гипергеометрическое распределение строго дискретно: оно не может быть расширено на реальные ценные аргументы его параметров или случайной величины.

Следующий график показывает, как изменяется распределение по мере изменения доли «дефектных» элементов, сохраняя при этом популяцию и размеры выборки постоянными:

Обратите внимание, что поскольку распределение симметрично по параметрам n и r, если мы изменим размер выборки и сохраним популяцию и пропорцию «дефектными», то получим в основном те же графы:

Member Functions

hypergeometric_distribution(unsigned r, unsigned n, unsigned N);

Конструирует гипергеометрическое распределение с популяцией объектов N, из которых r являются дефектными, и из которых n отобраны.

unsigned total()const;

Возвращает общее количество объектов N.

unsigned defective()const;

Возвращает число дефектных объектов r в популяции N.

unsigned sample_count()const;

Возвращает количество объектов n, которые отобраны из популяции N.

Non-member Accessors

Поддерживаются все функции доступа , которые являются общими для всех дистрибутивов: Кумулятивная функция распределения , Кумулятивная функция плотности , Квантиль , Опасная функция , Кумулятивная функция опасности , Средний , Средний , Мод , Вариантность , стандартное отклонение , Куртоз , Куртоз_избыток , Дальность и Поддержка .

Доменом случайной переменной являются неподписанные целые числа в диапазоне [max(0, n + r - N), min(n, r)]. Ошибка domain_error повышается, если случайная величина находится вне этого диапазона или не является интегральным значением.

Caution

Функция квантиля по умолчанию возвращает целочисленный результат, который был округлен наружу. То есть нижние квантили (где вероятность меньше 0,5) округлены вниз, а верхние квантили (где вероятность больше 0,5) округлены вверх. Такое поведение гарантирует, что если запрашивается квантиль X%, то по меньшей мере запрашиваемое покрытие будет присутствовать в центральном регионе, а не более запрашиваемое покрытие будет присутствовать в хвостах. Такое поведение может быть изменено таким образом, что функции квантиля округляются по-разному с помощью Политики. Настоятельно рекомендуется прочитать учебник Понимание квантилей дискретных распределений перед использованием функции квантиля на гипергеометрическом распределении. Ссылочные документы описывают, как изменить политику округления для этих распределений. Однако следует отметить, что метод реализации функции квантиля всегда возвращает интегральное значение, поэтому попытка использовать Политика, которая требует (или производит) реальный результат, приведет к ошибке компиляции времени.

Accuracy

Для малого N такой, что N < boost::math::max_factorial<RealType>::value тогда поиск по таблице результатов дает точность нескольким эпсилонам. boost::math::max_factorial<RealType>::значение составляет 170 при двойной или длинной двойной точности.

Для большего N так, что N ::мат::буст::мат::мафия тогда для расчета требуется только базовая арифметика и точность обычно составляет <20 эпсилон. Позаботится о N до 104729.

Для N > boost::prime(boost::math::max_prime), после чего точность быстро ухудшается, при этом 5 или 6 десятичных цифр теряются для N = 110000.

В общем случае для очень большого N пользователь должен ожидать потери log₁₀ N десятичных цифр точности при вычислении, при этом результаты становятся бессмысленными для N>= 10¹⁵.

Testing

Существует три набора тестов: наша реализация тестируется на основе таблицы значений, созданных Mathematica. Мы также проверяем нашу реализацию по некоторым значениям пятен, вычисленным с помощью онлайн-калькулятора здесь http://stattrek.com/Tables/Hypergeometric.aspx. Наконец, мы проверяем точность некоторых высокоточных тестовых данных с использованием этой реализации и NTL::RR.

Implementation

PDF можно рассчитать непосредственно по формуле:

Однако это может быть использовано непосредственно только тогда, когда самый большой из факториалов гарантированно не переполнит используемое представление с плавающей точкой. Эта формула используется непосредственно, когда N < max_factorial<RealType>::value в этом случае поиск таблиц факториалов дает быстрый и точный метод реализации.

Для более крупного N используется метод, описанный в «Точное вычисление функции гипергеометрического распределения», Trong Wu, ACM Transactions on Mathematical Software, Vol. 19, No. 1, March 1993, Pages 33-43. Метод опирается на то, что существует простой способ факторизации факториала в произведение простых чисел:

Где p_i - простое число i, а e_i - небольшое положительное целое число или ноль, которое можно вычислить посредством:

Далее мы можем объединить факториалы в выражении для PDF, чтобы получить PDF непосредственно как продукт простых чисел:

При этом экспоненты e_i либо положительные, либо отрицательные, либо нулевые. Действительно, такая степень отмены происходит при расчете e_i, что многие из них равны нулю, и, как правило, большинство имеют величину или не более 1 или 2.

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

Этот метод можно использовать до тех пор, пока N меньше самого большого простого числа, которое мы храним в нашей таблице простых чисел (в настоящее время 104729). Метод относительно медленный (вычисление экспонентов требует наибольшего времени), но для вычисления результата требуется только небольшое количество арифметических операций (действительно, нет более короткого метода, включающего только базовую арифметику, как только экспоненты были найдены), поэтому метод гораздо более точен, чем альтернативы.

Для гораздо большего N мы можем вычислить PDF из факториалов, используя либо lgamma, либо путем прямого объединения приближений lanczos, чтобы избежать вычисления с помощью логарифмов. Мы используем последний метод, так как обычно это 1 или 2 десятичных цифры более точные, чем вычисления с помощью логарифмов с lgamma. Однако в этой области, где N > 104729, пользователь должен ожидать потери около log₁₀ N десятичных цифр при расчете в наихудшем случае.

CDF и его дополнение рассчитываются путем непосредственного суммирования PDF. Мы начинаем с определения того, является ли CDF или его дополнение меньшим из двух, а затем вычисляем PDF при k (или k+1, если мы вычисляем дополнение) и вычисляем последовательные значения PDF через отношения повторения:

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

Квантиль рассчитывается аналогично CDF: сначала мы угадываем, к какому концу распределения мы ближе, а затем суммируем PDF-файлы, начиная с конца распределения на этот раз, пока у нас не будет некоторого значения k, которое дает требуемый CDF.

Медиана - это просто квантиль на 0,5, а остальные свойства рассчитываются через: