Блог

В настоящее время я борюсь с вычислением биномиальных коэффициентов для очень больших чисел, скажем, «n выберите k» с n < 10 000 000 и n < k. Это необходимо в контексте вычисления гипергеометрических распределений вероятностей.

До этого момента я перепробовал множество подходов для обработки больших чисел, которые получаются в результате этих вычислений. Однако проблема в том, что мне не нужно вычислять эти биномиальные коэффициенты один раз, а сотни тысяч раз. Это означает, что обычные подходы к вычислению факториалов слишком дороги, а стандартные типы данных, такие как long long int , слишком ограничены и не могут содержать эти числа.

Я уже пробовал использовать многоточные типы данных из Boost библиотеки, но, как я упоминал, выполнение вычислений так много раз приводит к чрезвычайно низкой производительности. Я также пробовал использовать многопоточность OpenMP , но выигрыш в производительности все еще был слишком низким. В результате я переключился на вычисление логарифма биномиальных коэффициентов, чтобы сохранить числа небольшими. Хотя это решило проблему больших чисел, это не ускорило процесс. Вот почему я попробовал аппроксимацию Стирлинга логарифмических биномиальных коэффициентов. Мое текущее решение выглядит так:

 #include <math.h>

long double calc_hgeom(unsigned int k, unsigned int n, unsigned int K, unsigned int N)
{
    long double hprob = std::exp((log_C(K, k)   log_C(N-K, n-k)) - log_C(N, n));
    return hprob;
}

long double log_C(unsigned int u, unsigned int m)
{
    long double C = u * std::log(u) - m * std::log(m) - (u-m) * std::log(u-m))   0.5 * (std::log(u) - std::log(m) - std::log(u-m) - std::log(2*M_PI));
    return C;
}
  

Однако результаты довольно сильно отличаются от фактических значений, вплоть до 7 %. Отсюда мой вопрос: существует ли эффективный способ вычисления логарифма биномиальных коэффициентов или можно улучшить мое приближение для повышения точности?

Любая помощь была бы очень признательна, поскольку это вычисление является основой всего моего алгоритма.

Рассмотрим функцию lchoose R …

 > choose(10000, 5000) 
[1] Inf
> lchoose(10000, 5000)
[1] 6926.641
  

Репозиторий исходных текстов базового языка R является отличным источником идей для подобных задач.

См https://github.com/wch/r-source/blob/trunk/src/nmath/choose.c

Хитрость здесь заключается в том, чтобы работать с ln-преобразованными входными данными, чтобы избежать переполнений.

Пожалуйста, обратите внимание, что код находится под лицензией GNU.

Вы должны использовать формулу аппроксимации Стерлинга для n! , которая, примененная к биномиальным коэффициентам, дает вам:

для самого биномиального коэффициента и для логарифма просто поместите логарифм правой части этого равенства; большая часть этого материала станет намного проще достаточно скоро. У вас все равно будет k! хотя, что — для больших k — вам снова понадобится формула аппроксимации. В конечном итоге у вас будет что-то более работоспособное (т. Е. Более численно стабильное).

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