Алгоритм для быстрых и точных вычислений по формуле Бернулли

[Vince Diesel]

math123 в сообщении #1529916 писал(а):

Речь идет о таких порядках, как $p=0.1, n=10000, k_1=1000, k_2=3000$.

В матпакетах при таких параметрах без проблем считается с любой точностью. Надо только задавать вероятность без десятичной точки: p=1/10, и с нужной точностью считать каждое слагаемое. Например, в математике при вычислениях с 50-ю знаками

p = 1/10; n = 10000; k1 = 1000; k2 = 3000;
res = Sum[N[Binomial[n, k] p^k (1 - p)^(n - k), 50], {k, k1, k2}]

выдает

0.50487591630932206692164052593040040929639368687300

меньше, чем за секунду.

[GAA]

И без символьных (или с переменной точностью) вычислений можно вычислить с требуемой точностью (8 знаков или около)
Matlab 2013b:

p = 0.1; n = 10000; k1 = 1000; k2 = 3000;
format long
cdf('beta', 1-p, n-k2, k2+1) - cdf('beta', 1-p, n-(k1-1), (k1-1)+1)

Результат: 0.504875916314735 (Вычисляет мгновенно, но верных цифр после запятой только 11.)

(старый) Matlab 6.5:

>> p = 0.1; n = 10000; k1 = 1000; k2 = 3000;
>> format long
>> betacdf(1-p, n-k2, k2+1) - betacdf(1-p, n-(k1-1), (k1-1)+1)
ans =
   0.50487591631200

Аналогично, вычисляет мгновенно, 11 цифр [верных] после запятой. Это же можно выполнить в Octave, но нужно загрузить пакет statistics.

Аналогично мгновенно вычисляет Mathcad (проверено в 13 и 15 версиях. В этой части движок численных вычислений не менялся, результаты совпадают)

Вложение:

MathCAD13.PNG [ 2.54 Кб | Просмотров: 782 ]

Всё это означает, что если пишется какая-то программа, например на C/Pascal/Fortran/… и требуется порядка 8 верных цифр, то можно юзать аппаратные типы данных (числа с плавающей точкой). Т.е. имеет смысл обсуждать алгоритмы вычисления B-функции (погрешности, скорость, и т.п.)

-- Wed 01.09.2021 09:53:20 --

Евгений Машеров в сообщении #1530193 писал(а):

В общем, настал момент, когда стоит спросить ТС "на фига?".

Может и не дождёмся ответа. Но помимо вычисления констант для критериев (нужна обратная), простейшие способы генерации биномиально распределённых основываются на функции B-распределения (или обратной к ней). Тут уже арифметика переменной точности слишком затратна, СКА не подходят. (Но лучше, конечно, обойтись без функции B-распределения или обратной.)

[Sender]

GAA в сообщении #1530187 писал(а):

Скорость тут дело второе.

Намёк на то, что важен также расход памяти? Ну нам не обязательно хранить весь треугольник, мы можем последовательно пройтись по диагоналям, отбрасывая значения, которые уже не будут использоваться, т.е. расход получится $O(n^2\log n)$. А вот расход вычислительных ресурсов $O(n^3\log n)$.

[GAA]

Sender, сразу напишу, что я так глубоко не думал. Я сразу написал: «требуется арифметика переменной точности». Для меня это плохо. :)
Быстро и достаточно точно нужно при имитационном моделировании. Программа, считающая что-то за разумное время, должна быть скомпилирована. Прилинковывать модули арифметики с переменной точностью (неограниченной точностью) без особой необходимости вроде как-то не разумно. Что это даёт?
Возьмём малые значения параметров (чтобы не загромождать текст). Maple 15:

> n:= 100: p:= 1/10: k1:= 10: k2:= 30:
> P:= (k1, k2, p)->add(binomial(n, k)*p^k*(1-p)^(n-k), k=k1..k2):
> S:= P(k1, k2, p);
   S = 1371774571277394842940725790581354557154632652158897364471415676574650904343047885951752764676932069/
   2500000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000

Потом можем преобразовать в число с плавающей точкой (с какой угодно точностью). Например, 64 цифры. Maple 15:

> evalf[64](S);
         .5487098285109579371762903162325418228618530608635589457885662706

(И для значений параметров из начального сообщения одна сумма вычисляется быстро, но привести точный результат в сообщении затруднительно: он будет занимать много-много места. Конечно, точный результат можно было и не выводить, а сразу преобразовать в число с плавающей точкой.) Если нужно посчитать несколько сумм, то можно просто воспользоваться одной из СКА.

Для меня быстрый расчёт нужен, если этих сумм вычисляется большое число (сотни миллионов). Основной расчёт в имитационной программе идёт на аппаратных типах данных (переменный с плавающей точкой одинарной или двойной точности, int32, int64). Т.е. все равно к этим типам всё и преобразовывать. Тут и задумываются о числе верных знаков для получения требуемых результатов. И вообще, не получится ли уйти от бета-функции / обратной бета-функции при помощи более хитрого алгоритма. :)

А моё отступление в этой теме в сторону символьной арифметике — это, в первую очередь, на случай необходимости тестирования программы. [Ну, может и ТС задал не столь точно вопрос и ему инженерного пакета достаточно.] Из примера в моём первом сообщении ветки также видно, что сортировка при данных значений параметров как-то особо ничего не даёт. [Но с увеличением диапазона суммирования может и дать.]
И инженерные пакеты при численном (а не символьном расчёте) не выдают результат с 14-15 верными знаками, и в справке об этом не пишут. И, посчитав символьно, можно сравнить. :)

В общем, то о чем Вы пишете, видимо, нужно для СКА. Я этим не интересовался; совершенно не владею информацией. Увы.

[Евгений Машеров]

Сугубое ИМХО.
Совет, который можно дать топикстартеру, зависит от того, каковы цели расчёта.
Для большинства практических задач хватит 2-3 знаков, которые вполне обеспечивают Муавр с примкнувшим к нему Лапласом. Для запрошенных восьми знаков

(Оффтоп)

("на фига?" - спрашивает поэт Андрей Вознесенский, рифмуя с "бытия")

действительно лучше аппроксимация через бета-распределение, но можно и напрямую считать, только не перемножая до упора, чтобы потом делить длиннейшие целые, а прологарифмировав и заменив факториалы Стирлингом (с поправкой по степеням 1/n). А задачи, где надо всё перемножить, а потом потерять почти всю точность, поделив - просто не вижу.

[GAA]

При указанных в начальном сообщении значений параметров выполняется вычисление неполной бета-функции

$I_x(a, b) = \frac 1 {B(a,b)} \int_0^xt^{a-1}(1-t)^{b-1}$

для $x \approx a/(a+b)$. Действительно, $x=1-p = 0.9$, $a = n-k$, $b = k+1$, при $k=999$ получим $a/(a+b) \approx 0.90001$.

В статье Didonato A. R., Morris A. H. "Algorithm 708 Significant Digit Computation of the Incomplete Beta Function Ratios" (1992) приводится алгоритм, обеспечивающий высокую точность для случая $x \approx a/(a+b)$. (Может кому-то будет интересно, не всегда в библиотеках (или в справках к функциям) приводят ссылки на источники. Ес-но в статье приведены ссылки и на более ранние алгоритмы.)