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

[math123]

Дана вероятность события $p$. Нужно найти точную (пусть до 8го знака после запятой) вероятность события возникнуть от $k_1$ до $k_2$ раз в $n$ испытаниях.

Математически задача решена. Есть формула Бернулли (https://ru.wikipedia.org/wiki/%D0%A4%D0 ... 0%BB%D0%B8). Можно просуммировать ее результат от $k_1$ до $k_2$.

Но как этим воспользоваться на практике, т.е. посчитать на компьютере? Речь идет о таких порядках, как $p=0.1, n=10000, k_1=1000, k_2=3000$. Вычисление $C(n,k)$ порождает огромные числа, происходит переполнение. Затем суммирование маленьких вероятностей наращивает погрешность, сумма от $0$ до $n$ не получаеся равной единице и т.д. Существует ли более оптимальное решение, чем "в лоб"?

[Pphantom]

Тема перемещена из форума «Математика (общие вопросы)» в форум «Карантин»
по следующим причинам:

- неправильно набраны формулы (краткие инструкции: «Краткий FAQ по тегу [math]» и видеоролик Как записывать формулы)
- заодно доформулируйте задачу: с какой точностью требуется получить ответ?

Исправьте все Ваши ошибки и сообщите об этом в теме Сообщение в карантине исправлено.
Настоятельно рекомендуется ознакомиться с темами Что такое карантин и что нужно делать, чтобы там оказаться и Правила научного форума.

[Pphantom]

Тема перемещена из форума «Карантин» в форум «Помогите решить / разобраться (М)»

-- 30.08.2021, 11:17 --

Например, представить факториалы в виде рядов Стирлинга, вынести (и посокращать) общие множители, ну и суммировать значения правильно - от меньших к большим.

[GAA]

Если есть быстрая и достаточно точная реализация функции бета-распределения, то проще всего, на мой взгляд, выражать искомые суммы через разность [значений] этой функции.

Если нужно несколько раз вычислить и компьютер не очень древний, т.е. не нужно писать программу, то можно немного подождать результатов в одном из пакетов, поддерживающих вычисления с большой точностью.
Например, в Matlab можно воспользоваться символьными вычислениями для получения каждого слагаемого, а затем преобразовывать в числа с плавающей точкой и далее добавлять в аккумулятор. Matlab 2013b:

Используется синтаксис Matlab M

n = sym(10000); k1 = 1000; k2 = 3000; S = 0; p= sym(0.1);
for i = k1:k2
 S = S + double(nchoosek(n, i)*p^i*(1-p)^(n-i));   
end    
format long
disp(S);

Результат: 0.504875916309321

Немного дольше будет выполняться расчёт, если не только слагаемые вычислять символьно, но и сумму накапливать в символьном виде. Matlab 2013b:

Используется синтаксис Matlab M

n = sym(10000); k1 = 1000; k2 = 3000; S = sym(0); p= sym(0.1);
for i = k1:k2
 S = S + nchoosek(n, i)*p^i*(1-p)^(n-i);   
end    
format long
disp(vpa(S, 24));

Результат: 0.504875916309322066921641

Можно воспользоваться не Matlab, а бесплатной системой Octave. Если есть Maple, то можно в нём попробовать (Maple 15):

> n:= 10000: k1:= 1000: k2:= 3000:
> P:= (k1, k2, p)->add(binomial(n, k)*p^k*(1-p)^(n-k), k=k1..k2):
> evalf[24](P(k1, k2, 0.1));
                      .504875916309322066921645

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

С практической точностью может хватить нормальной аппроксимации. По крайней мере для описанного случая

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

$p=0.1, n=10000, k_1=1000, k_2=3000$.

Есть "нормальная аппроксимация с поправками". Много в Univariate Discrete Distributions. Norman L. Johnson, Adrienne W. Kemp, Samuel Kotz. Кажется, на русский перевели.
Требование непременно считать "по формуле" не особо обосновано, но отчего нет? Барон Нэпер как раз для этого алгоритмы изобрёл. Логарифмы фигурирующих чисел будут не особо велики, а потенцировать уже после сокращения.
Ну и в порядке мемуарного анекдота - в студенческие годы работал в хоздоговорной теме, руководимой неким доцентом по прозвищу "полный доц". Вот его расчёт через нормальное распределение, через Стирлинга и т.п. не устроил, только считать точно по формуле, "мне не нужно приближённо, мы будем делать точно!". Сделали так - есть "склад числителей" и "склад знаменателей". Если накапливаемое произведение велико - берём со "склада знаменателей", деля число, если мало - со "склада числителей". $C^k_n=\frac {n!} {k!(n-k)!}$ можно сразу сократить на k!, так что "склад числителей" это числа от $n-k+1$ до n, "склад знаменателей" это числа от 1 до $n-k$. Таким образом переполнения, возникавшего при попытке в лоб всё перемножить, удалось избежать. Но при сравнении результата с расчётом по нормальной аппроксимации разницы не выявилось.

[GAA]

В конкретном случае аппроксимации Муавра — Лапласа как бы заведомо не хватает

> Digits:= 24;
> n:= 10000; k1:= 1000; k2:= 3000;
> PL:= proc(k1, k2, p, n)
             local x1, x2, q;
             q:= 1-p;
             x2:= (k2 - n*p)/sqrt(n*p*q);
             x1:= (k1 - n*p)/sqrt(n*p*q);
             evalf(1/sqrt(2*Pi)*int(exp(-z^2/2), z=x1..x2));
           end:
> PL(k1, k2, 0.1, n);
                 .499999999999999999999999

Получаем одну верную цифру две верных цифры, а по условию надо 8 верных цифр.
Через бета-распределение вычисляется (одна сумма) в Maple 15 мгновенно.

> with(Statistics):
> F:= (m, n, p)-> CDF(RandomVariable('Beta'(n-m, m+1)), 1-p):
> Digits:= 24:
> F(3000, 10000, 0.1)-F(999, 10000, 0.1);
                                  .504875916309322066932496

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

А p у нас где задаётся?

[GAA]

В смысле? В [моём] предыдущем сообщении в процедурах PL и F — это третий параметр (p).

-- Mon 30.08.2021 15:32:36 --

Банально, но на всякий случай. При нахождении разности значений функции распределения теряем верные знаки. Если вычисляем одну сумму, то можно посмотреть на значения функции распределения в двух точках и прикинуть потерю знаков. Если нужен расчёт в программе, то придётся делать оценки для получения необходимого числа цифр при вычислении значений функции распределения. Расчёт выполняется достаточно быстро и при 128 знаках, и при 256. Если вычислять много раз (при моделировании), то видимо нужно подбирать удовлетворительную аппроксимацию для функции бета-распределения (неполной бета-функции).

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

Тут одна неточность в применении аппроксимации Муавра-Лапласа. Вычитаемое надо брать не от 1000, а от 999, тогда получаем 0.504875916
Если вычитать значение для 1000 - то теряется значение вероятности для k=1000, а надо убирать меньшие 1000, но не само 1000.
При этом матожидание случайной величины 1000, то есть медиана близка к моде (а в данном случае вообще совпадает), и выбрасывание значения для k=1000 приводит к потере 0.013296956.

[GAA]

Вроде нет: Предельная теорема Муавра — Лапласа.
Если вычислить по приведенной в сообщении выше функции PL , беря $k_1=999$, то получим 0.513... (т.е. ещё большее отклонение от истинного значения).

-- Tue 31.08.2021 09:50:02 --

(В любом случае, 8 верных цифр после запятой, используя М — Л, заведомо не получить.)

[Sender]

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

"мне не нужно приближённо, мы будем делать точно!"

Можно ещё про треугольник Паскаля вспомнить, правда, тут о скорости говорить не приходится.

[Someone]

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

Тут одна неточность в применении аппроксимации Муавра-Лапласа. Вычитаемое надо брать не от 1000, а от 999, тогда получаем 0.504875916

У меня получилось $0{,}51329561381709210506961858464535$ (Wolfram Mathematica). Вы правы в том, что промежуток нужно увеличить на $1$, но конкретный способ неудачен.

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

Вроде нет: Предельная теорема Муавра — Лапласа.
Если вычислить по приведенной в сообщении выше функции PL , беря $k_1=999$, то получим 0.513... (т.е. ещё большее отклонение от истинного значения).

-- Tue 31.08.2021 09:50:02 --

(В любом случае, 8 верных цифр после запятой, используя М — Л, заведомо не получить.)

Более точная аппроксимация, чем обычно приводимая в учебниках:

$${\mathbf P}(k_1\leqslant S_n\leqslant k_2)\approx\Phi\left(\frac{k_2-np+\frac 12}{\sqrt{np(1-p)}}\right)-\Phi\left(\frac{k_1-np-\frac 12}{\sqrt{np(1-p)}}\right),$$

$${\mathbf P}(S_n\leqslant k)\approx\frac 12+\Phi\left(\frac{k-np+\frac 12}{\sqrt{np(1-p)}}\right),$$

$${\mathbf P}(S_n\geqslant k)\approx\frac 12-\Phi\left(\frac{k-np-\frac 12}{\sqrt{np(1-p)}}\right),$$

где

$$\Phi(x)=\frac 1{\sqrt{2\pi}}\int_0^xe^{-\frac{t^2}2}dt,$$

$p$ — вероятность успеха в каждом опыте, $n$ — число опытов, $S_n$ — число успехов в $n$ независимых опытах, $0\leqslant k_1\leqslant k_2\leqslant n$, $0\leqslant k\leqslant n$.

Эта формула даёт ${\mathbf P}(1000\leqslant S_{10000}\leqslant 3000)\approx 0{,}50664873019368255382207350819953$ при "точном" значении $0{,}50487591630932206692164052593040$. Так что $8$ правильных цифр тоже не получается.

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

Someone в сообщении #1530141 писал(а):

У меня получилось $0{,}51329561381709210506961858464535$ (Wolfram Mathematica). Вы правы в том, что промежуток нужно увеличить на $1$, но конкретный способ неудачен.

Ну, я не упомянул про "поправку на непрерывность", прибавку 0.5. Но да, она даёт три знака, это я ошибся, посчитал в лоб через сумму биномиальных, получил не то, понял, что значение 1000 не надо вычитать, исправил, получил искомое точное и обрадовался, поскольку понял причину расхождения в расчёте через ЦПТ. А нельзя ли попросить ТС пояснить, зачем ему 8 знаков?
Ещё упоминается более хитрая поправка, у Джонсона, Кемпа и Котца, на стр. 117, но то ли там опечатка, то ли я не понимаю, как её применить.

[GAA]

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

но то ли там опечатка, то ли я не понимаю, как её применить

В справочнике ссылаются рядом на David Peizer и John Pratt “A normal approximation for binomial, F, beta, and other common, related tail probabilities, I”, 1968, но там я этого приближения не нашёл. Во второй части статьи (John Pratt “A normal approximation for binomial, F, beta, and other common, related tail probabilities, II”, 1968) также не нашёл. Да и в целом в Univariate Discrete Distributions приводятся аппроксимации без указания погрешности, но со ссылками на источники. Как всегда нужно смотреть оригинальные статьи.

Вообще, всевозможные «продвинутые» аппроксимации в основном разрабатывались для создания и использования таблиц. David Peizer и John Pratt в начале статьи “A normal approximation for binomial, F, beta, and other common, related tail probabilities, I”, 1968 как раз явно об этом говорят.

Функция биномиального распределения [cumulative distribution function] точно выражается через неполную бета-функцию. Об этом уже давно написано в Википедии:

wikipedia in Binomial distribution писал(а):

It [cumulative distribution function, — GAA] can also be represented in terms of the regularized incomplete beta function, as follows:
$F(k;n,p) = \Pr(X \le k) = ... = (n-k) {n \choose k} \int_0^{1-p} t^{n-k-1} (1-t)^k \, dt.$

В книгах об этом написано было гораздо раньше. Например, в Большев и Смирнов «Таблицы математической статистики» (1983) указывается, что вычисления функции биномиального распределения могут быть выполнены также по таблицам неполной B-функции. У меня сейчас нет ранних изданий ТМС, но Б.Л. ван дер Варден в «Математической статистике» (1960) уже отсылает к неполной B-функции.

Для вычисления неполной B-функции разрабатывались различные алгоритмы/аппроксимации. Если нужно многократно вычислять, то можно подобрать побыстрее, но с меньшей точностью или более медленный, но с большей точностью.
Неполная B-функции или функция бета-распределения реализованы в большинстве пакетов. Мне всегда казалось, что нужны какие-то исключительные условия, чтобы понадобились «продвинутые» аппроксимации нормальным распределением.
Не «продвинутые» аппроксимации нормальным распределением страдают тем, что при близком к нулю $p$ они имеют низкую точность даже при больших $n$. Для вычисления пределов (в книге Феллера есть поучительные) «продвинутые» аппроксимации нормальным вроде не нужны, а для вычисления вероятностей для конечных $n$ нужно знать погрешность таких «продвинутых» аппроксимаций.

-- Wed 01.09.2021 06:47:57 --

Sender в сообщении #1530140 писал(а):

Можно ещё про треугольник Паскаля вспомнить, правда, тут о скорости говорить не приходится.

binomial(10000, 5000) = 159179026353243894833759727364152118865300583745761455042831910351777263712009579866326285394422221774335
859829932262055804632908708020739850879872195958489620417578664585801840995875120689143315978135317405145
347319967052139450253847727733600831205378448823951274321755502883180927364644281795459349368900235462880
547366282927213220919726803062157839769855248683450847868894994611262023360235298989458928488427591110374
321646235202929095545845304023492927787143123978410362592908300075421733055365492425368306281530729653340
889255650690875150647615944622376204326852230062678211259375951657115342848245333181068684095284004284699
504359257817996430741389422649447586626281862183757541280362546881388544759125956185871468454381861463662
350728468211441655465743993284005794170022128691686189379747227886202239788372897602049671018976190617859
305826168808117556117796960379809282174855477301204105813490546271598511886613777441541105636943056820725
244819431050256487494579628837604295079872914178005301024149340722579759834860211640098545723183096418633
688898312145597072469454456651789081935386062566029368316522506271595824234037562793787332887113614352737
971292965638066368798136853809235306441396478979814279989804419587974310478889401271971015441216840096344
652939528524306710003806696307699257220104426311836533049067512198270012436774453339363870022811792535618
814009571973175044979333952276086203573893932977683234377126461503016956149960119508206705891127875644018
328002477885570580594271739655617247279703665698618080801965541235756564655565433970795513642117996823482
940891493286717047038936158996297545140449708716896119990505242038078767450450863985246304067167020269491
256064620583001761300622284757510662566106193771435587218537809620026913816305961756296827876710659465040
754767228071475821687019166324258201685893281454941849633219010250326331594361831605955344426680189751351
988451293306946591872301020473208721181284611163964165765568933940740976656925878728168406935207314430178
725136178015792747114729015831170907171194578298482944164359840658473384707719418659651955333974514346503
817616197612616157040354559466774548777412765471478663541418800111962602957335265945686584369721309686983
612640564990207924247805354140963069566603071195931569172626802351518208786515546937379638760504643715547
953097876650816797000176926659286918757094175117347665748132703540903393455409827319346571309202004128279
611588827972847323501797969972562671972826347017756606331304016075551520523318404592756797612244679324194
846919392918520452394577675953326869067443192793756095658856432124228522403516658454319704009054696329636
363817791559641205005685702690372838060388519713403611629040056633420468941761593824568608770545269390456
038883755973215629222766342326791030991205489279359135464145696802130792488795541350742383065293811197486
421347908348956557941526999776837834147059039199747891501916363639677591945387535180152498052210450701705
508838093544209022455222930021060372371375638589078163387440553649120

Для оперирования с такими числами требуется арифметика переменной / неограниченной точности. Скорость тут дело второе.

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

Они там вообще на Лапласа ссылаются. Но проработать "Аналитическую теорию..." мне не по силам. Тем более, что не уверен, что там это действительно есть.
Но опять же у меня вопрос - зачем 8 знаков (на фига, как спрашивал поэт Андрей Вознесенский, рифмуя, однако с "бытия").
Что до конкретного расчёта - да, если нужна точность 8 знаков, бета самое лучшее. Аппроксимировать нормальным - точность ниже, хотя при p=0.1 и 10000 испытаний для большинства практических приложений удовлетворительна. Можно считать "в лоб" по формуле, если только не ставить задачу найти сперва один факториал перемножением, затем второй и третий, затем умножать на маленькие числа в больших степенях, тут без "длинной арифметики" не обойтись, но если надо 8 знаков - логарифмы такую точность обеспечат, а после потенцирования результат в массив, сортировать по возрастанию и суммировать с малых (хотя, если честно, это уже перестраховка).
В общем, настал момент, когда стоит спросить ТС "на фига?". И после ответа конкретизировать рекомендации.