На примере хорошей задачки по теории вероятностей мы сравним что круче - аналитические решения или моделирование. Вначале я хотел написать только аналитическую часть и попросить заинтересованных читателей запрограммировать симуляцию задачи. Но нехорошо перекладывать на других свои дела, поэтому, в меру способностей, постарался сделать и вторую часть. Попутно знакомясь с Питончиком. Если что-то в коде получилось коряво, милости просим наставить на путь истинный.
Задачка такая:
В ящике 4 белых 10 черных шаров. Из него наудачу вынимают шар, фиксируют его цвет и возвращают шар назад в ящик. Назовем «белым пулом» любую максимальную цепочку подряд вынутых белых шаров. Найти математическое ожидание количества «белых пулов» при извлечении из ящика 20 шаров.
Аналитический подход
Чтобы решить такую задачку, нужно прочувствовать ее изнутри - представить себя белым шаром, черным шаром, всеми шарами сразу вместе... Картинка поможет настроиться на нужный лад и не подглядывать в решение, если захотелось порешать сначала самому.
Внимание, мой вариант решения. Возможны арифметические ошибки, но ответ похож на правду и будет проверен моделированием.
Мне кажется, задачка нетривиальная. По крайне мере, все виденные мной на просторах интернета попытки решения были весьма плачевны. Люди получали в ответе математическое ожидание числа белых пулов, большее чем их максимально возможное число...
Первый шаг, который я сделал - интуитивно почувствовал, что нужно использовать индуктивный подход, начав с небольшого числа извлеченных шаров. Хорошо бы еще формализовать задачу. Сделать это проще всего перейдя к рассмотрению двоичных последовательностей см.также предыдущий пост про слоистые среды. Ноль у меня будет обозначать черный шар, единичка - белый. Тогда, например, в последовательности 1-0-1-0-111-0-11-0-1-000 будет 5 белых пулов с длиной от 1 до 3. Изолированная с двух сторон нулями единичка - уже белый пул! Поэтому при вытаскивании 20 шаров максимальное число пулов будет равно 10 чередование нулей и единиц. Соответственно ответ, который хитрый студент найдет в гугле = 100/7, будет весьма далек от правды.
Идем дальше. При извлечении 2х шаров нам нужно рассмотреть 4 варианта: 0=00, 1=01, 2=10, 3=11. Соответственно, при извлечении 3х шаров нужно рассматривать двоичные последовательности задающие числа от 0 до 7, и т.д.
В случае последовательности из 2х шаров матожидание числа белых пулов найти просто. У нас есть ноль пулов, и три раза по одному пулу. А дальше мы делаем твист, который может сбить с толку. Обозначим это матожидание как <N(2)>, но вычислять его не будем - вдруг и не понадобится потом!
Если взять большее число испытаний, то возникает сложность с одновременным подсчетом числа пулов и соответствующих вероятностей. Понятно, что для каждой двоичной последовательности вероятности появления нулей и единиц нужно перемножить. Но их нужно как-то подсчитать. К счастью, в задаче есть повторяющаяся структура. При увеличении числа испытаний сначала все двоичные последовательности повторяются, потом слева появляется дополнительная единичка и снова все повторяется.
Логично попробовать вывести рекуррентную формулу для матожидания <N((n))>. Сделать это можно так — при n испытаниях вы должны выписать двоичные представления чисел от 0 до
Обозначим
Следующий шаг уже вполне предсказуем. Наша стратегия заключается в том, чтобы не выводить формулы для числа пулов и числа единичек! Используем аналог ленивых вычислений - ленивую аналитику. Запишем матожидание числа пулов
Тоже самое сделаем для матожидания при увеличенном на 1 числе испытаний
Теперь нужно в формуле для <N((n+1))> попытаться выделить <N((n))>. Для этого нужно понять как работает число пулов, при добавлении еще одного двоичного разряда. Если число i по прежнему может быть представлено n разрядами, то добавление нового двоичного разряда просто соответствует приписыванию 0 слева. Число белых пулов и число единичек не меняется
Дальше в двоичном разложении слева появляется изолированная единичка, что увеличивает число белых пулов и число единиц на, простите, единицу
И в оставшемся диапазоне, левая единичка сливается с предыдущей единичкой, в итоге
Теперь все это добро надо подставить в формулу для математического ожидания длины белых пулов при n+1 испытаниях. После некоторых алгебраических преобразований, у меня получилось, что
Это уже замечательно! Число белых пулов вычислять не надо. Но есть сумма, в которую входит выражение для числа единиц. Для оставшейся суммы я не придумал ничего лучше, чем снова найти рекуррентное соотношение
Теперь если мы увеличим число испытаний на единичку, число слагаемых вырастет в два раза. Первая половина суммы не поменяется, а у вторая половина будет точно такая же как первая, только с дополнительным множителем за счет добавления одной единички
Получилась геометрическая прогрессия!
Теперь можно выписать искомое рекуррентное соотношение без наших вспомогательных функций p и u:
Это арифметическая прогрессия! Значит см. UPD1.
И, окончательно, <N(20)>=200/49~4, то есть в среднем мы ожидаем 4 белых пула. Ответ похож на правду, если учесть что число возможных пулов от 0 до 10.
Моделирование
Давайте теперь решим задачу с помощью моделирования! Надо запустить генератор двоичных последовательностей длиной в 20 символов с вероятностью появления нулей и единичек из задачи. Я просто сделал ящик=последовательность с 10 нулями и 4 единичками и случайным образом вытаскиваю 20 раз по одному элементу. Появление и исчезновение белого пула я подсчитываю по изменению следующего элемента. При этом каждый пул кроме, возможно, последнего, будет учтен два раза.
import random
box = [1,1,1,1,0,0,0,0,0,0,0,0,0,0]
def poolnumber():
unitcount=0
poolcount=0
pooldelim=0
for i in range(20):
ball=random.choice(box)
#print(ball)
unitcount=unitcount+ball
if ball!= pooldelim:
poolcount+=1
pooldelim=ball
res=round((poolcount+0.1)/2)
#print(unitcount)
#print(round((poolcount+0.1)/2))
return resНу а теперь много раз вызовем функцию подсчета пулов и найдем среднее
totalpoolnumber=0
size=5000
for i in range(size):
totalpoolnumber=totalpoolnumber+poolnumber()
print(totalpoolnumber/size)
print(204/49)Удивительно, но вероятные арифметические ошибки в аналитической части не случились! И 500 и 20000 испытаний дают 4.0xxx - близкий ответ к аналитическому решению! Однако, насколько проще и быстрее реализовать модель! К тому же, можно уже полученную модель применить для решения более сложных задач - искать, например, среднюю длину белого пула, или среднюю длину пула максимальной длины. Так что, хоть я и занимаюсь "точно решаемыми" задачами, но должен признаться - использование программирования и моделирования является мощным и универсальным инструментом.
PS: Моделирование сначала давало заниженный ответ, близкий к 3.8. И только знание точного ответа, хотя я в нем и засомневался, позволило мне найти ошибку - шарик я вытаскивал не 20 раз, а всего 19. Когда это поправил, сразу добился замечательного согласия теории и эксперимента! И в защиту аналитического метода решения, скажу, что ответ у нас получился для любого числа испытаний - мы нашли, что среднее число белых пулов растет с числом вытащенных шаров линейно. В принципе, этого и следует ожидать от аналитического подхода - попадание в общие свойства решения.
UPD 1. В развертывании арифметической прогрессии я потерял начальное слагаемое. Правильная формула
приводит к ответу 204/49~4.16.
Да, такая простая формула наталкивает на мысли, что все можно сделать проще, и действительно в комментариях есть простой вывод рекуррентного соотношения. Нужно сразу думать о том, что может происходить с математическим ожиданием числа пулов, если будет вытягиваться на один шарик больше. А также есть альтернативное решение строго в парадигме теории вероятностей.
В моделировании я поправил пределы в цикле для подсчета пулов, но опять ошибся на единичку в подсчете среднего. Эта ошибка не должна была влиять для большого размера выборки, но ответ получался систематически ниже. Конечно, это из-за корявого использования округления, которое я исправил "регуляризацией" вместо корректного решения. Зато теперь точно все сходится. Спасибо за примеры корректных реализаций модели на самых разных языках программирования!
Комментарии (61)
gorilych
27.08.2021 07:43+3Кмк, неправильная (путающая) формулировка
Назовем «белым пулом» любую максимальную цепочку подряд вынутых белых шаров
Вместо "любая максимальная" лучше "граничащая с началом, концом, или вынутым чёрным шаром".
vesper-bot
27.08.2021 09:37+2Слово "максимальная" запутало и меня — для меня "максимальная" цепочка это та, которая максимальной длины из имеющихся в последовательности, т.е. если вытащили 10101000, длина максимальной цепочки 1 и всего их 3, а если 10101100, то соответственно 2 и 1. И увидев решение, я всерьез прифигел, а оказалось, считали просто количество переходов с черного/ничего на белое при вытаскивании.
Tiriet
27.08.2021 11:22имеется в виду, что в цепочке 11110 можно разглядеть и 1'1'1'1'0 и 11'11'0 и 11110, и что такую цепочку надо считать не как четыре цепочки по одному шару подряд, а как одна цепочка из четырех шаров подряд. Но в целом я с Вами согласен- определение несколько "корявое"
ProRunner
27.08.2021 10:01+2Более понятно будет, наверное, "любая цепочка белых шаров (из одного и более шара), разделенная черными шарами"
vesper-bot
27.08.2021 10:04В такой формулировке для вытаскивания одних белых шаров число цепочек 0.
ProRunner
27.08.2021 10:21-1Вероятность один из 100 миллионов, если правильно посчитал (0.4 перемножить 20 раз), я не думаю, что как-то может повлиять на результат)
Tiriet
27.08.2021 11:24+1а как Вам такой поворот: назовем "пулом" непрерывную последовательность белых шаров, не являющуюся частью более длинной такой последовательности.
annalen_der_physik Автор
27.08.2021 17:59Да, в формулировке задачи есть сложности с интерпретацией. На 99% я уверен, что интерпретировал правильно, потому что нужно считать число этих самых пулов и цепочка должна быть «любой максимальной подряд».
В любом случае, можно подойти творчески и попробовать варианты задачи!
Tiriet
27.08.2021 09:42+5как-то у Вас очень все усложнено. Рекурсии, суммы какие-то дикие.
Обратите внимание на простой факт- белый пул оканчивается сочетанием 10, еще один белый пул возможен, если последний шар- 1 (тогда после него не вытаскивается новый шар и это тоже пул). Таким образом, количество белых пулов равно числу пар вида "10" в последовательности + число окончаний на "1". В последовательности на 20 шаров есть 19 пар и один финальный беспарный шар.
Каково матожидание числа пар "10" из партии в 20 шаров? вероятность вытащить случайно белый шар- 4/14. вероятность вытащить после этого черный шар- 10/14. вероятность парного события- 4/14*10/14=40/196, то есть, вероятность того, что случайна пара шаров будет "белым пулом"- 40/196. Всего мы вытаскиваем 19 пар (последний шар- беспарный!), вероятность того, что из 20 шаров мы вытащим сколько-то полноценных белых пулов равна 19*40/196. Вероятность того, что последний шар тоже будет концом белого пула- 4/14, итого, мое мат.ожидание числа белых пулов- 19*40/196 + 4/14 = (760+56)/196 = 816/196 ~= 4.163.
vesper-bot
27.08.2021 09:49Правильно так — белый пул начинается сочетанием 01, и первый шар может быть белым пулом. И пропущена информация о том, что если некая пара — белый пул, следующая пара заведомо НЕ белый пул, т.е. количество попыток вытащить пару всё ещё 19, но условная вероятность следующей пары быть белым пулом если предыдущая им оказалась есть ноль.
Tiriet
27.08.2021 10:19че-то я промазал с ответом- ниже мои расчеты. Ясен пень, подтверждающие мою правоту. Жду теперь Ваших расчетов, подтверждающих Вашу правоту :-).
vesper-bot
27.08.2021 12:09Ну, откровенно говоря, ниже вы измеряете именно то, что требовалось, т.е. количество пар в 20-символьной строке плюс последняя. А первая фраза решается разворотом строки, что очевидно не меняет количество последовательностей единиц в ней. Другое дело, что я запилил на коленке "полный анализ" тупым методом "в лоб" и получил ME= 4.163265306113811… что отличается от 816/196 на машинный нуль. Прекрасно :)
ProRunner
27.08.2021 09:52+1Долго думал, почему моя реализация давала ответ в 4.16..., когда судя по статье должно быть 4.0....
mentin
06.09.2021 02:26К этой формуле наверное надо добавить что-то про эргодичность и возможность подсчета мат-ожидания таким образом. Потому что формула предполагает что все 19 пар независимы и даже 19 пулов могу реализоваться одновременно с вероятностью (40/196)^19, что конечно не так. Описанные здесь вероятности 19 пар не независимы - если первая пара оказалась пулом, то вероятность что вторая окажется пулом равна нулю. Формула похоже работает, но надо показать почему.
Tiriet
06.09.2021 11:43рассмотрим три шара. это две пары. вероятность вытащить первой парой пул-
40/196- только БЧ*. два таких случая на все поле событий.
второй парой пул будет только в двух случаях- ББЧ и ЧБЧ. вероятность того, что вторая пара- пул равна вероятность того, что первый шар белый или черный, помноженная на вероятность того, что второй шар белый и третий черный. то есть, (4/14+10/14)*4/14*10/14= 40/196.
добавим еще один шар в выборку, и посмотрим, какова вероятность того, что третья пара из четыре шаров- пул: в нашей выборке таких четыре события, которые реализуются в случаях, если первые два шара- ББ, БЧ, ЧБ, ЧЧ. однако, в любом случае первые два шара реализуют одну из этих комбинаций, и вероятность этого равна 1.0. то есть, вероятность того, что третья пара из четырех шаров- пул, равна 1,0*Рб*Рч = 40/196. и мне все равно, что я вытащил второй парой. ну и по индукции дальше.
Tiriet
27.08.2021 09:59+1Назовем «белым пулом» любую максимальную цепочку подряд вынутых белых шаров
где тут сказано, что перед белым пулом должен быть ноль? нету такого.
Про условные вероятности- а зачем мне рассматривать их? Я просто хочу знать среднее число "хвостов" на цепочке из 19пар+1. без всяких условностей.
в наших с Вами "аналитических значениях" разница наблюдается в третье цифре, при этом дисперсия случайной величины у нас достаточно большая, поэтому для проверки аналитики экспериментом нужно гонять не по 5к прогонов, а по 50к (а лучше- по 500к).
Tiriet
27.08.2021 10:12+3procedure TForm1.Button1Click(Sender: TObject); const L = 20; const Pwhite = 4/14; const Count = 500000; var i, j: integer; b_next,b_prev: boolean; // true= White; false = black Pool_count: integer; begin randomize; Pool_count := 0; for i := 1 to Count do begin b_next := random < Pwhite; // цвет первого шара. for j := 2 to L do begin b_prev := b_next; // шар становится предыдушим b_next := random < Pwhite; // цвет нового шара if b_prev and (not b_next) // White-Black- the end of pool. then inc(Pool_count); end; if b_next then inc(Pool_count); // Last ball is white- it is pool too; end; Memo1.Lines.Add( Format('%6.5g', [Pool_count / Count]) ); end;4,1653 4,1621 4,164 4,1609 4,162 4,1635 4,1622 4,1618 4,1613 4,1649 4,1645
vesper-bot
27.08.2021 12:10За Дельфи отдельный респект.
saltpepper
27.08.2021 18:12+1А за программирование в обработчике клика кнопки и вывод в Memo1… наверное, приз зрительских симпатий.
vesper-bot
29.08.2021 07:53Не только лишь все умеют в многопоточность, мало кто может это делать.
PS: в Дельфи это довольно неудобная вещь, когда-то пилил, там надо было и написать много чего, и всерьез заняться контролем границ, да и вообще передачей данных между потоками, так что для задач на "разок посчитать" подход "посчитать прямо в тыкалке" валиден.
pavelgein
27.08.2021 10:39+1Могу предложить другое решение.
Пусть
— случайная величина, равная числу белых пулов. Требуется найти математическое ожидание этой случайной величины. Обозначим его через 
Обозначим вероятность вытащить белый шар через
Рассмотрим всевозможные белые пулы и перенумеруем их каким-нибудь образом. Пусть всего их
тогда через 
обозначим случайную величину, которая равна 1, если этот белый пул реализовался при извлечении шаров, и 0 в противном случае. Ясно, что 
Воспользуемся линейностью математического ожидания:
Теперь осталось вычислить
Пусть 
отвечает за пул длины 
Рассмотрим случаи 
Возможно два случая: если пул стоит с краю, то чтобы он реализовался, достаточно чтобы только один из соседей был черным. Это происходит с вероятностью 
и возможно два таких случая (с одного и другого края). Если пул стоит в середине, то его местоположение можно выбрать 
способов, оба его соседа должны быть черными, это реализуется с вероятностью 
В этом случае есть два способа выбрать белый пул, а оставшийся шар должен быть черным. Это реализуется с вероятностью 
В этом случае все шары белые, это возможно в единственном случае, который реализуется с вероятностью 
Суммируя, получаем математическое ожидание,
В случае автора, как раз получается
Tiriet
27.08.2021 11:18+2случаи 1 и 2 у Вас идентичные, просто число n-k-1 для случая 2 равно нулю. но это так, замечание, которое лишь слегка облегчает Ваши выкладки. А вот не мелочь:
Давайте в Вашу формулу подставим n=3 и p=0.5. тогда получим, что среднее число пулов равно 3*0.5*0.5 = 0.75. Простой перебор вариантов 000 001 010 011 100 101 110 111 говорит, что на 8 равновероятных раскладов у меня 8 пулов, и МО равно 1.0. не сходится Ваша теория с моей практикой.
pavelgein
27.08.2021 11:47+2Спасибо, вы правы!
Я неправильно преобразовал последнюю сумму.
Правильно вот так:
Разница в слагаемом
В случае автора — 
Осталось найти несоответствие в рассуждениях.Tiriet
27.08.2021 14:54Собственно, ровно то, что и у меня- только я не расписывал.
(n-1)p(1-p)+p= (p-p^2)(n-1)+p = pn-p-np^2+p^2+p = np(p-1)+p^2
annalen_der_physik Автор
27.08.2021 19:11Со всеми задачами по теории вероятностей теперь сразу к вам за консультацией!
sashagil
28.08.2021 21:11+1Меня этот вывод от pavelgein не убеждает (наоборот, вижу путаницу в определениях случайных величин, а "догадаться, что хотел сказать автор, но оговорился" нет времени, желания), а вот разбор от Pavgran ниже по комментарием выглядит корректно, плюс, в результате - компактная формула с очевидной линейной асимптотикой. Кстати, догадка о линейности асимптотики - довольно "физическая", можно рассматривать, как ведёт себя величина при возрастании N вдвое при больших значениях N - разделив испытание с 2N шарами на комбинацию двух испытаний с N шарами. При больших N одной возможной серией белых шаров, продолжающейся с первого N-испытания на второе можно пренебречь, посчитав её два раза - линейность тогда становится наглядной (что не умаляет необходимости доказать её строго, с формулой от Pavgran).
Интуиция подсказывает мне, что прав Pavgran, попробую найти время / силы перепроверить.
GarryC
27.08.2021 10:46+2Что то я недопонял - в последней программе проводят 0.5*10**6 опытов несколько раз, хотя полное поле событий насчитывает 2**20=(2**10)**2=(10**3)**2=1*10**6, что всего в 2 раза больше одного опыта. Так почему не провести исчерпывающий анализ?
Tiriet
27.08.2021 10:48:-) мы все замерли в ожидании!
vesper-bot
27.08.2021 12:13+2WHITEPROB=4.0/14.0 me=0.0 tp=0.0 for i in range(2**20): st=bin(i)[2:] pc=st.count("1") prob=WHITEPROB**pc*(1.0-WHITEPROB)**(20-pc) prev="0" v=0 for c in st: if (c=="1") and (prev=="0"): v+=1 prev=c me+=v*prob tp+=prob print("ME=",me,"total=",tp)Результат: ME= 4.163265306113811 total= 0.9999999999959284
total здесь — прямая сумма посчитанных вероятностей с погрешностями на уровне машинного нуля, при математических подсчетах total=1.GarryC
27.08.2021 14:39+1А по времени счета сколько заняло ? Ну и в порядке улучшения точности - считать вероятности в длинных целых числах (4 и 10) а потом делить сумму на 14**20
vesper-bot
27.08.2021 14:47+2Секунд пять (много). Ну и дробь правильнее в этом случае упростить до 2/7, тогда результат даже влезет в 64 бита "чистого" целого (720<820==2**60). Но уже влом :)
speshuric
28.08.2021 02:58+1Я вот тоже согласен, что в этой задаче проще в лоб посчитать, чем колхозную монтекарлу городить и чем в формулах путаться.
Я посчитал всё в целых числах (64-битных). Там, действительно, можно всё посчитать точно (в рамках рациональных чисел, где числитель и знаменатель целые 64-битные). Писал не на пайтоне, а на котлине (привычнее), получилось так:
fun pow(n: Long, p: Int) = generateSequence { n }.take(p).fold(1L) { acc, v -> acc * v } tailrec fun countPool(v: Long, lastInPool: Int = 0, cnt: Int = 0): Int = if (v == 0L) cnt else countPool(v / 2, v.mod(2), cnt + (1 - lastInPool) * v.mod(2)) fun main() { val size = 20 val denominator = pow(7L, size) val black = 5L val white = 2L val ps = MutableList(size / 2 + 1) { 0L } //0..10 for (i in (0L until pow(2L, 20))) { val popcnt = i.countOneBits() ps[countPool(i)] += pow(black, size - popcnt) * pow(white, popcnt) } println("count pool\tnumerator\tprobability") ps.forEachIndexed { ind, v -> println("$ind\t$v\t${v.toDouble() / denominator.toDouble()}") } val m = ps.reduceIndexed { i, acc, v -> acc + v * i } // числитель матожидания println("$m/$denominator~=${m.toDouble() / denominator.toDouble()}") }Но у меня почему-то получились другие результаты:
count pool numerator probability 0 95367431640625 0.0011951964277455186 1 1229180230500756 0.015404753963449494 2 6287186216221020 0.07879443093859312 3 16493725701609600 0.20670832484053933 4 24061484648040000 0.30155158844962027 5 19906942939800000 0.24948461628537263 6 9188837373800000 0.11515949853494763 7 2251181376000000 0.028213027157337086 8 265696380000000 0.0033298512791828257 9 12514000000000 1.5683224177797937E-4 10 150000000000 1.8798814341295274E-6 332291741405372221/79792266297612001~=4.164460502550194Здесь колонки:
count pool - счетчик пулов белых
numerator - числитель вероятности
probability - вероятность
В сумме numerator дают 7 ^ 20, т.е. вероятность суммы - единица. Т.е. кажется, я посчитал вероятности правильно.
От 204/49 отличается ровно на 95367431640625/79792266297612001, т.е. на вероятность того, что пулов не будет. То ли я ошибся, то ли остальные - не пойму.
Время работы варианта - около 0,4 с, но это я не упарывался по производительности и ПК древний.
speshuric
28.08.2021 03:07+1Ну да, если закешировать степени, то получается 0,1 с.
fun pow(n: Long, p: Int) = generateSequence { n }.take(p).fold(1L) { acc, v -> acc * v } tailrec fun countPool(v: Long, lastInPool: Int = 0, cnt: Int = 0): Int = if (v == 0L) cnt else countPool(v / 2, v.mod(2), cnt + (1 - lastInPool) * v.mod(2)) fun main() { val size = 20 val denominator = pow(7L, size) val black = 5L val white = 2L val ps = MutableList(size / 2 + 1) { 0L } //0..10 val probQuants = List(size + 1) {pow(black, size - it) * pow(white, it)} for (i in (0L until pow(2L, 20))) { ps[countPool(i)] += probQuants[i.countOneBits()] } println("count pool\tnumerator\tprobability") ps.forEachIndexed { ind, v -> println("$ind\t$v\t${v.toDouble() / denominator.toDouble()}") } val m = ps.reduceIndexed { i, acc, v -> acc + v * i } // числитель матожидания println("$m/$denominator~=${m.toDouble() / denominator.toDouble()}") }speshuric
28.08.2021 23:10+1А не надо было Хабр читать и кодить в 2 часа ночи. Вместо fold был reduce и считалось неправильно.
Правильный код:
fun pow(n: Long, p: Int) = generateSequence { n }.take(p).fold(1L) { acc, v -> acc * v } tailrec fun countPool(v: Long, lastInPool: Int = 0, cnt: Int = 0): Int = if (v == 0L) cnt else countPool(v / 2, v.mod(2), cnt + (1 - lastInPool) * v.mod(2)) fun main() { val size = 20 val denominator = pow(7L, size) val black = 5L val white = 2L val ps = MutableList(size / 2 + 1) { 0L } //0..10 val probQuants = List(size + 1) {pow(black, size - it) * pow(white, it)} for (i in (0L until pow(2L, 20))) { ps[countPool(i)] += probQuants[i.countOneBits()] } println("count pool\tnumerator\tcount event\tprobability") ps.forEachIndexed { ind, v -> println("$ind\t$v\t${v.toDouble() / denominator.toDouble()}") } val numerator = ps.foldIndexed(0L) { i, acc, v -> acc + v * i } // числитель матожидания var numeratorReduced = numerator var denominatorReduced = denominator // деноминатор это 7^20, поэтому gcd считать просто while (numeratorReduced.mod(7)==0) { numeratorReduced/=7L denominatorReduced/=7L } println("$numerator/$denominator = $numeratorReduced/$denominatorReduced ~= ${numerator.toDouble() / denominator.toDouble()}") }Считает около 0,1 секунды и выводит:
count pool numerator probability 0 95367431640625 0.0011951964277455186 1 1229180230500756 0.015404753963449494 2 6287186216221020 0.07879443093859312 3 16493725701609600 0.20670832484053933 4 24061484648040000 0.30155158844962027 5 19906942939800000 0.24948461628537263 6 9188837373800000 0.11515949853494763 7 2251181376000000 0.028213027157337086 8 265696380000000 0.0033298512791828257 9 12514000000000 1.5683224177797937E-4 10 150000000000 1.8798814341295274E-6 332196373973731596/79792266297612001 = 204/49 ~= 4.163265306122449annalen_der_physik Автор
01.09.2021 03:37Надо бы все комментарии собрать в виде пособия — N подходов к решению задач по теории вероятностей!
Tiriet
27.08.2021 15:19+1тогда это уже не будет "Аналитика vs моделирование", это будет аналитика vs полный перебор. А может просто ошибся человек- ноль потерял, вместо 200 написал 20. На цепочке из 200 шаров его аналитика и мое моделирование будут работать, а вот полное поле событий мы уже не потянем.
Pavgran
27.08.2021 18:11+2Какое-то у Вас сложноватое аналитическое решение. Через рекурсию проще:
Пусть E(n|b) - условное мат. ожидание количества пулов, если последний шар чёрный, и E(n|w) - условное мат. ожидание кол-ва пулов, если последний шар белый. Пусть p - вероятность вытянуть белый шар.
E(1|b)=0, E(1|w)=1
E(n+1|b) = E(n) =E(n|b)*(1-p)+E(n|w)*p - если последний шар чёрный, то мат. ожидание такое же, как если бы шара не было.
E(n+1|w) = E(n,w)*p + (E(n,b)+1)*(1-p) = E(n)+(1-p) - если последний шар белый, то мат. ожидание не меняется, если до этого тоже был белый шар, и увеличивается на 1, если до этого был чёрный шар.
E(n+1)=E(n+1|b)*(1-p)+E(n+1|w)*p=E(n)*(1-p)+E(n)*p+(1-p)*p=E(n)+p(1-p)
Учитывая E(1)=p, получаем E(n)=p+(n-1)p(1-p)=p^2+np(1-p)
E(20)=4/49+20*2/7*5/7=204/49
sashagil
28.08.2021 21:26+1Изящное и, по-моему, правильное решение, не вижу изъянов - лучший комментарий здесь (я, правда, пролистал бегло)! Бонус - очевидная линейная асимптотика в результирующей компактной формуле (автор заметки упоминает этот аспект, хорошо понятный при "физическом" рассмотрении того, что происходит при масштабировании N, например, рассмотрении испытания с 2N шарами как двух испытаний с N шарами).
serebrennikov1973
28.08.2021 22:09+2Эту линейность легко объяснить в рамках теории вероятностей. Дело в том, что искомую случайную величину N(n) (число белых пулов в n испытаниях) можно представить как сумму n случайных величин с распределением Бернулли (принимающих значение 1, если в данном испытании создаётся новый белый пул, и 0 в противном случае), по одной для каждого испытания. В силу свойства линейности математическое ожидание N(n) равно сумме математических ожиданий этих n случайных величин. При этом оказывается, что значение p (вероятность значения 1) для всех испытаний, кроме первого, одинаково и равно 5/7×2/7 (новый белый пул возникает в конкретном испытании тогда и только тогда, когда в данном испытании выпадает белый шар, при этом в предыдущем испытании выпал чёрный шар, а цвета шаров более ранних испытаний не важны(!), что и даёт одинаковые значения p для всех испытаний, кроме первого). Для первого испытания p просто равно 2/7 (вероятность выпадения белого шара).
serebrennikov1973
27.08.2021 18:11+1По моим подсчётам, у Вас есть ошибка в формуле, должно получиться <N(n)>=2/7+(n-1)10/49, и <N(20)>=204/49~4.16327, где конкретно ошибка - не искал. Можно легко проверить правильность формулы для небольших n, посчитав матожидание через пространство элементарных событий (<N(1)>=1×2/7+0×5/7=2/7, <N(2)>=1×2/7×2/7+1×2/7×5/7+1×5/7×2/7+0×5/7×5/7=24/49, <N(3)>=1×2/7×2/7×2/7+1×2/7×2/7×5/7+2×2/7×5/7×2/7+1×2/7×5/7×5/7+1×5/7×2/7×2/7+1×5/7×2/7×5/7+1×5/7×5/7×2/7+0×5/7×5/7×5/7=238/243 и т.д.). В этой задаче проще всего считать матожидание аддитивно от вклада каждого последующего испытания, учитывая, что один новый белый пул появляется в первом испытании при белом шаре, при последующих испытаниях - тогда и только тогда, когда белый шар следует непосредственно за чёрным, причём цвет шаров до чёрного не важен(!). Первое испытание увеличивает искомое матожидание на 1×2/7, второе и все последующие - на 1×5/7×2/7, отсюда и полученная формула.
annalen_der_physik Автор
01.09.2021 03:47Спасибо! Вот в момент обдумывания условий, мне эта аддитивность совсем не была очевидна. Аддитивность — это мы независимо рассматриваем пулы разной длины, а вдруг часть из их реализаций перекрывается, не будем ли мы учитывать лишние события и т.д. Потом уже, после решения, раз я пришел к той же формуле (по модулю ошибки), аддитивность мне стала действительно понятна — так как те функции (число пулов и число единиц в данном событии), которые казались мне важными для решения, просто исчезли.
serebrennikov1973
02.09.2021 14:08+1Да, использование свойства линейности математического ожидания - очень полезный метод, который часто позволяет обойтись без громоздких вычислений при решении задач на нахождение математического ожидания. Правда, нужно аккуратно выбрать разбиение искомой случайной величины на более элементарные (обычно с распределением Бернулли), чтобы вычисление p для них получалось простым (как в данной задаче), что бывает не всегда просто сделать. А некоторые задачи без использования линейности математического ожидания вообще не решить либо решить крайне сложно, в качестве примера приходит в голову такая вполне классическая задача: в закрытой коробке лежат n перепутанных шнурков, 2n их концов выведены через дырки наружу. Все концы произвольным образом попарно связывают, при этом образуются петли разной длины (возможно, зацепленные). Определить математическое ожидание общего числа петель.
serebrennikov1973
02.09.2021 15:33+1Нет. Пока, наверное, не буду приводить здесь ответ и решение - задача красивая, если не видели её раньше - будет интересно порешать еще. Если нужно, готов прислать ответ и/или решение.
sashagil
08.09.2021 00:49Я тоже, наверно, порешаю. Пока догадка (после беглой ручной проверки n от 1 до 3, в которой я мог допустить ошибку) такая: 2 - 1/n!. Реализую и погоняю перебор проверить или заменить догадку, потом буду разбираться с результатом. Не припомню, чтобы такую задачу встречал, но вспомнились задачи с похожими формулировками - про подброшенные в воздух шляпы и про сумасшедшую старушку и прочих пассажиров авиарейса.
serebrennikov1973
08.09.2021 16:55Удачи!
sashagil
12.09.2021 20:05(вы уже второй раз помещаете новый комментарий в корень обсуждения :) )
Исправляю догадку, теперь сделав аналитические выкладки с рекурсивным вычислением величины от n. Теперь у меня получается Hn, частичная сумма первых n элементов гармонического ряда. Асимптотическая формула есть, log n + 0.5772... (постоянная Эйлера-Маскерони). Посчитаю моделированием для проверки - в следующий раз. Рекурсивные выкладки: данный шнурок с вероятностью 1/n добавляет новую петлю, а с вероятностью (n-1)/n не добавлет петлю к ситуации с n-1 шнурками. После приведения получаем X(n) = 1/n + X(n-1). Кажется, у меня есть на полке книга, в которой я мог бы подсмотреть похожие рассуждения - в одном из трёх томов "Искусства программирования" Кнута попадалось. Сейчас должен заняться другими вопросами - буду проверять тред на днях, не появилось ли что-либо от вас.serebrennikov1973
12.09.2021 23:37Да, увидел, спасибо! Теперь у Вас получилось практически все верно, только вероятность петли на первом выбранном шнурке будет 1/(2n-1): первый конец берём любой, а вот второй нам подходит только 1 из оставшихся (2n-1). Можно через рекурсию, а можно и сразу в явном виде всю сумму записать: первый завязанный узел даст петлю с вероятностью 1/(2n-1), но в любом случае после его зявязывания остаётся (n-1) свободных шнурков (один, возможно, двойной длины) и ( 2n-2) свободных концов, так что второй узел даст петлю с вероятностью 1/(2n-3) и так далее до 2 свободных концов, когда последний, n-й, узел даст петлю с вероятностью 1, в результате получим 1/(2n-1)+1/(2n-3)+...+1/1. В обоих вариантах мы используем свойство линейности математического ожидания: искомая случайная величина (число петель) представляется как сумма n случайных величин с распределением Бернулли (число петель, которые завязываются каждым из n узлов). А вот без использования линейности решать очень сложно: число элементарных исходов посчитать легко (1/(2n-1)×1/(2n-3)×...×1/1), но вот прямой подсчёт числа петель в каждом исходе крайне затруднён.
sashagil
13.09.2021 19:47Спасибо за разбор! Путаница n с 2n случилась у меня в голове с формулировкой задачи - представлял коробку с n концами шнурков слева и n концами шнурков справа, и случайное связывание левых с правыми. Исходное условие проще формулируется, а я, похоже, переключился на другое условие, которое мне было проще решать.
lxsmkv
Мне показалось необычным, что в тексте я не нашел слов "распределение", "Бернулли", "биномиальный", "случайная величина", поэтому подумал, что Вы решили не искать вероятностную модель подходящую для решения этой задачи, а сочли более интересным придумать свою? Это не упрек, просто интересно почему.
Да, программное моделирование довольно интересная штука. Я как-то делал программный эксперимент где на квадратном поле, со 100 одинаковыми ячейками, стоят люди, в каждой ячейке по одному, а сверху на них одновременно кидают 100 кирпичей. Вопрос: сколько человек выживет. Получился для меня неожиданный результат, 63~. Оказывается соотношение ячеек в которые кирпич не попадает - стремится к четкому значению в 36~ процентов. Так я открыл для себя Закон Малых Чисел, или "Закон двух третей", найденый в свое время Владиславом Борткевичем.
Также, компьютерное моделирование вероятностных экспериментов помогает с задачами, где нас подводит интуиция. Ну и вообще помогает структурировать задачу, потому что для компьютера нужно четко определить переменные, как минимум.
sergeyns
Это вы гауссиану открыли
Tiriet
С компьютерными экспериментами по статистике еще очень желательно определить необходимый размер выборки. 200/49 отличается от 4,0 в третьем знаке. Среднеквадратичное отклонение в данной задаче- ~1,3 (30%!). Для отлова третьего знака МО при таких сигмах надо сотни тысяч экспериментов, а не тысячи.
longclaps
Кирпичи настигнут 1/e от популяции, что довольно очевидно.