Кто такие Dракоши смотрите в первой части.

1. Зачем нам вероятности?
2. Постановка задачи
3. Подсчет способов расстановки агентов
4. Моделирование вероятностей методом Монте-Карло
5. Что получилось?
1. Зачем нам вероятности?
Как уследить за Dракошами
Прежде чем начать запускать моделирование надо бы понять, а как мы будем смотреть на мир Dракош? Ведь их будут тысячи, каждый будет заниматься своими делами, куда-то идти, что-то кушать и т.д.
Поэтому в самом начале пришлось задумался о параметрах, которые будут описывать текущее состояние моделируемой системы и показывать динамику происходящих в ней изменений. Получилось всего 24 величины, за которыми интересно проследить, помимо самого пространства с распределением агентов и пищи. Среди этих параметров есть вполне очевидные: количество агентов, количество пищи в пространстве, рождаемость и смертность Но о них более подробно будет рассказано в третей части. Сейчас же я больше остановлюсь на специальной группе параметров, которые позволят оценить степень осознанности действий агентов. Другими словами, отличить случайное поведение от осознанного.
Напомню, агенты могут совершать шесть разных действий: Pass, Step 1, Step 2, Eat, Attack, Sex. Для каждого из них можно рассчитывать долю агентов, которые его совершает в течении каждого такта времени:
Все действия кроме Pass могут быть успешными или неуспешными, в зависимости от сложившейся ситуации в текущей и в соседних ячейках, а также от желаний соседа по занимаемой ячейке. Для отслеживания успешности используются четыре параметра:
Индексы осознанности или зачем нам вероятности
На начальном этапе, когда нейронная сеть агентов еще не обучена, действия агентов будут носить случайный характер. А по мере эволюционного развития, под действием отбора, можно ожидать увеличения степени осознанности и следовательно повышения успешности действий. Но необходимо учесть, что успешность действий зависит так же и от окружающей обстановки.
Так успешность действий Step 1/2 определяется наличием в ячейке, которую агент выбрал для перемещения, свободного места (всего в ячейках есть по два места для агентов). И будет зависеть как от правильности выбора агента, так и от количества агентов в пространстве. Текущее количество агентов удобно отслеживать с помощью параметра
Здесь
Таким же образом успешность действий Attack и Sex будет зависеть от вероятности (
Осознанность совершения атаки можно определить как превышение успешности
А вот успешность
2. Постановка задачи
Для вычисления введенных выше индексов осознанности необходимо научится рассчитывать две вероятности:
Сформулируем задачу следующим образом.
Имеется пространство изячеек. В одной ячейке может находится не более
агентов. В этих ячейках располагаются
агентов (
). Необходимо определить вероятность того, что в выбранной ячейке с одним агентом будет находится еще один агент
. И вероятность того, что в выбранной ячейке будет
агентов
.
Для решения данной задачи я применил два подхода. Первый это расчет вероятности в классическом ее понимании, как отношение количества различных расстановок агентов, которые благоприятствуют искомому событию, к полному количеству расстановок агентов по ячейкам. Для вычислений воспользовался методом производящих функций.
Во втором подходе использовалось понятие частотной вероятности. Расчеты выполнялись с применением метода Монте-Карло. Вероятно, есть и другие способы вычислить эти вероятности, и возможно даже более простые. Буду благодарен за наводку на таковые.
3. Подсчет способов расстановки агентов
Вывод рекуррентной формулы
Предложенный ниже подход к расчету вероятностей позволяет решить поставленную задачу даже в более общей постановке, когда вместимости ячеек различны. Зададим вместимости ячеек пространства в виде последовательности
Поставленную задачу можно свести к подсчёту количества способов расставить агентов по ячейкам пространства. Тогда вероятность
Вероятность
Обозначим через
Рассмотрим процесс, при котором мы последовательно перебираем все ячейки пространства (от
Количество способов заполнить ячейку агентами это сочетание без повторений из
Поместив
Для того что бы воспользоваться полученной рекуррентной формулой необходимо задать граничные условия:
- Если
, т.е. у нас не осталось агентов которых можно поместить в ячейку, тогда у нас есть всего один способ выполнить расстановку агентов никого не помещать в ячейку:
- Для последней ячейки
так же есть только один способ расставить агентов оставить все
агентов в ней, но только если им хватит места:
- Если места не хватает для размещения
агентов в ячейках с
по
, то способов заполнить ячейку нету:
Теперь можем записать выражения для вероятности
где
Обратите внимание, здесь начинаем суммировать с единицы, т.к. нас интересуют те расстановки, при которых в ячейку
Выражения для вероятности
где
Результаты вычислений сведены в таблицу:


Упрощение и переход приведенному виду
Количество расстановок очень быстро растет и уже для
Первое что можно сделать это вынести из-под знака суммы
Но прежде введем обозначение приведенного количества расстановок
Полное количество способов расстановки теперь будет:
Аналогично для
И можем сократить выражения для вероятностей на факториал:
Где
После введения новых обозначения требуется переписать выражения для граничных условий. Для дальнейших выкладок удобнее представить граничные условия в следующем виде:
-
, при
,
;
-
, при
,
;
-
, при
,
;
-
, при
.
Полученные выражения для приведенных количеств расстановок смотрятся поинтереснее т.к. коэффициенты в рекуррентной формуле постоянны и не растут со скоростью факториала как в формулах выше. Но все равно быстро и надежно вычислить вероятности при размерах пространства порядка 104 не очень получится.

Следующим шагом предлагаю избавится аж от трех видов выражений для приведенного количества расстановок:
Готовые выражения для вероятностей теперь зависят только от приведенного количества способов расстановки агентов одного вида
Построение производящей функции
Осталось научится рассчитывать элементы последовательности
Хорошим пособием для начинающих комбинаторов по производящим функциям можно назвать учебник Ландо С.К. [3]. Еще один ресурс типа справочник по производящим функциям. И еще вот этот пост на хабре. Но наиболее полезными для меня были видео лекции Игоря Клейна по производящим функциям, выложенные тут.
Но в перечисленных ресурсах в основном речь ведется о производящих функциях для одномерных последовательностей, с одним индексом. А наша последовательность двумерная, с двумя индексами. Производящие функции от двух переменных лишь упоминаются. Поэтому пришлось рискнуть и попробовать вывести производящую функцию для нашей двумерной последовательности на основе одномерных примеров.
Идея метода в том, что получив каким-то из способов выражение для производящей функции
Суммируя то, что имеется на данный момент по нашей последовательности приведенных количеств расстановок запишем рекуррентную формулу для
$$display$$\left\{ \begin{array}{} f_{0,0}=1 \\ f_{i,0}=f_{i-1,0} \\ f_{i,1}=f_{i-1,1}+f_{i-1,0} \\ f_{i,n}=f_{i-1,n}+f_{i-1,n-1}+\frac{1}{2}f_{i-1,n-2}, \text{ при } n\geqslant 2 \\ f_{i,n}=0, \text{ при } n>2i \end{array} \right.$$display$$
Наша последовательность
Используя специальный технический примем можно из рекуррентного выражения для элементов последовательности получить выражение для производящей функции. У меня получилось вот это:
$$display$$\begin{array}{cccccccc} i=0: \qquad & x^0y^0 \times \qquad & f_{0,0}= & \color{Blue}{1} \\ i=1: \qquad & x^1y^0 \times \qquad & f_{1,0}= & \color{Blue}{f_{0,0}} \\ \qquad & x^1y^1 \times \qquad & f_{1,1}= & \color{Blue}{0} & + & \color{Green}{f_{0,0}} \\ \qquad & x^1y^2 \times \qquad & f_{1,2}= & \color{Blue}{0} & + & \color{Green}{0} & + & \color{Red}{\frac{1}{2}f_{0,0}}\\ i=2: \qquad & x^2y^0 \times \qquad & f_{2,0}= & \color{Blue}{f_{1,0}} \\ \qquad & x^2y^1 \times \qquad & f_{2,1}= & \color{Blue}{f_{1,1}} & + & \color{Green}{f_{1,0}} \\ \qquad & x^2y^2 \times \qquad & f_{2,2}= & \color{Blue}{f_{1,2}} & + & \color{Green}{f_{1,1}} & + & \color{Red}{\frac{1}{2}f_{1,0}}\\ \qquad & x^2y^3 \times \qquad & f_{2,3}= & \color{Blue}{0} & + & \color{Green}{f_{1,2}} & + & \color{Red}{\frac{1}{2}f_{1,1}}\\ \qquad & x^2y^4 \times \qquad & f_{2,4}= & \color{Blue}{0} & + & \color{Green}{0} & + & \color{Red}{\frac{1}{2}f_{1,2}}\\ i=3: \qquad & x^3y^0 \times \qquad & f_{3,0}= & \color{Blue}{f_{2,0}} \\ \qquad & x^3y^1 \times \qquad & f_{3,1}= & \color{Blue}{f_{2,1}} & + & \color{Green}{f_{2,0}} \\ \qquad & x^3y^2 \times \qquad & f_{3,2}= & \color{Blue}{f_{2,2}} & + & \color{Green}{f_{2,1}} & + & \color{Red}{\frac{1}{2}f_{2,0}}\\ \qquad & x^3y^3 \times \qquad & f_{3,3}= & \color{Blue}{f_{2,3}} & + & \color{Green}{f_{2,2}} & + & \color{Red}{\frac{1}{2}f_{2,1}}\\ \qquad & x^3y^4 \times \qquad & f_{3,4}= & \color{Blue}{f_{2,4}} & + & \color{Green}{f_{2,3}} & + & \color{Red}{\frac{1}{2}f_{2,2}}\\ \qquad & x^3y^5 \times \qquad & f_{3,5}= & \color{Blue}{0} & + & \color{Green}{f_{2,4}} & + & \color{Red}{\frac{1}{2}f_{2,3}}\\ \qquad & x^3y^6 \times \qquad & f_{3,6}= & \color{Blue}{0} & + & \color{Green}{0} & + & \color{Red}{\frac{1}{2}f_{2,4}}\\ \ldots \qquad & \ldots \qquad & \ldots \end{array}$$display$$
Просуммируем полученные выражения в столбик следующим образом. Отдельно суммируем левый от равно столбец и получим собственно производящую функцию
Суммируя синие слагаемые справа от равно, получим:
Замете, что начиная со второго слагаемого можно вынести за скобки множитель
Суммируя зеленые слагаемые, получим:
Суммируя красные слагаемые, получим:
После выполненных суммирований получим уравнение относительно искомой функции:
$$display$$\mathscr{F}(x,y)=\color{Blue}{1+x\mathscr{F}(x,y)}+\color{Green}{xy\mathscr{F}(x,y)}+\color{Red}{\frac{1}{2}xy^2\mathscr{F}(x,y)}$$display$$
Решая которое относительно
Теперь дело за малым, разложить полученную производящую функцию в ряд по степеням
С учетом этих обозначений производящая функция запишется:
Вот тут лежит аж семь объяснений того почему
Теперь разберемся с
А это уже бином Ньютона:
Аналогично для
Осталось подставить полученное выражение для
Теперь все что нам нужно это выделить коэффициент при
Для этого необходимо сгруппировать все слагаемые, у которых степень
В результате получим два выражения, одно для чётного
и одно для нечетного
Полученные выражения уже можно попробовать вычислить, но так как биноминальные коэффициенты, которые входят в них, всё же растут довольно быстро, для расчетов необходимо воспользовался инструментами символьных вычислений MatLab.
function out = f(i,n)%% Символьное вычисление f(i,n)nu = ceil(n/2); % нижний предел суммированияsyms kif mod(n,2)==0 % если n чётное syms S(t,h) S(t,h) = symsum(nchoosek(t,k)*nchoosek(k,2*h-k)/2^(2*h-k), k, h, 2*h); out = S(i,nu);else % если n нечётное syms S(t,h) S(t,h) = symsum(nchoosek(t,k)*nchoosek(k,2*h-k-1)/2^(2*h-k-1), k, h, 2*h-1); out = S(i,nu);end
Вычисление
function out = W1sym(ii,N)%% Символьное вычисление вероятности обнаружить второго агента в ячейке пространстваout = 1/(1 + f(ii-1,N-1)/f(ii-1,N-2));
Вычисление
function out = W2sym(ii,N)%% Символьное вычисление вероятности обнаружить второго агента в ячейке пространстваout = 1/(1 + 2*f(ii-1,N)/f(ii-1,N-2) + 2*f(ii-1,N-1)/f(ii-1,N-2));
4. Моделирование вероятностей методом Монте-Карло
Метод Монте-Карло или метод статистических испытаний это большая группа численных методов решения математических задач (причем не только вероятностной природы) при помощи моделирования случайных величин. Один из вариантов метода Монте-Карло предполагает что смоделированная случайная величина с известным распределением (зачастую равномерно распределённая на отрезке 0..1) используется для установления статистических характеристик других случайных величин или для приближённой оценки некоторого значения
Для установления статистических характеристик случайных величин, например, для определения вероятности
На сегодня метод имитационного моделирования развился в самостоятельное направление методов исследований. Но для решения поставленной в разделе 2 задачи использовался подход, который с одной стороны является примером реализации методов Монтер-Карло, с другой это пример имитационного моделирования, т.к. использовался естественный процесс изменения положения агентов в пространстве.
Обзор алгоритма
В рамках метода Монте-Карло поставленная в разделе 2 задача сводится к оценке частоты реализации соответствующих событий. Так вероятность
Будем рассматривать процесс заполнения

После каждого этапа пересаживаний будем подсчитывать количество
Частота обнаружения ячеек с двумя агентами (вероятность
Можно предложить несколько вариантов реализации этапа пересаживания:
- Scenario 1 почти полная имитация: для пересаживания случайно выбираются 50% агентов и каждый из них пересаживается в одну из 18 соседних клеток выбранных случайно;
- Scenario 2 пересаживаются все агенты и каждый из них пересаживается в одну из 18 соседних клеток выбранную случайно;
- Scenario 3 для пересаживания случайно выбираются 50% агентов и каждый из них пересаживается в случайную ячейку пространства;
- Scenario 4 пересаживаются все агенты и каждый из них пересаживается в случайную ячейку пространства;
Scenario 1 назван почти полной имитацией процессов в настоящем мире Dракош, так как доля агентов совершающих действия Step 1/2 различна для каждого такта времени и может иметь любое значение а не только 50%. Ниже мы посмотрим, как результаты вычислений зависят от сценария пересаживания агентов.
Процесс начинаем со случайной расстановки агентов по ячейкам. Затем в цикле
- Пересаживаем агентов согласно выбранного сценария;
- Определяем количество заполненных ячеек
и соответствующие значения
и
.
Возникает вопрос: сколько же раз необходимо повторить цикл? С ростом
На графике ниже показана зависимость вероятности
Горизонтальной чертой отмечено точное значение, рассчитанное способом, описанным в разделе 3.

По результатам расчетов получены следующие значения неопределенности (надежность 99.7%):
- Для
неопределенность составила 0.0014 (0.24%);
- Для
0.00047 (0.080%);
- Для
0.00015 (0.025%);
- Для
0.000066 (0.011%);
Т.к. время вычислений растет пропорционально количеству испытаний (для
Результаты Монте-Карло
Давайте посмотрим, как зависит рассчитываемые значения вероятностей от используемого сценария пересаживания агентов. Ниже показаны результаты вычислений для разных сценариев, при

Как можно видеть отличия между рассчитанными значениями вероятностей с использованием разных сценариев намного меньше неопределенности этих результатов. Это позволяет заключить что более точная имитации мира Dракош привела бы к такому же результату. Поэтому не зачем тратить время.
5. Что получилось?
На графике ниже показаны зависимость вероятностей

Совпадение результатов, полученных разными методами, позволяет надеяться, что серьезных ошибок не было допущено.
Теперь можем посмотреть, как выглядит успешность действий Step 1/2, Attack и Sex при полностью случайном (неосознанном) поведении агентов. А точнее на вторые слагаемые в выражениях для индексов осознанности, которые как раз и задают базовый уровень успешности, обусловленный вероятностями обнаружения агентов в пространстве:
Для однозначности примем что частота действия Sex

В дальнейшем мы сравним полученные результаты с успешностью Тупиков (раса агентов без нейронной сети) и полноценных Dракош и посмотрим, как будет расти успешность в результате развития нейронной сети агентов.
- Ландо С.К. Лекции о производящих функциях. 2007. 144 c.
- Соболь И.М. Численные методы Монте-Карло. 1973. 312 с.
- Бусленко Н.П., Шрейдер Ю.А. Метод статистических испытаний (Монте-Карло) и его реализация на цифровых вычислительных машинах. 1961. 226 с.