How to become an author
Search
Write a publication
Pull to refresh

Часть 2. Dракоши. Комбинаторика мира Dракош

Reading time14 min
Views1.9K
Mathematics*Matlab*Artificial Intelligence
В этой части я расскажу как рассчитывал вероятности обнаружения агентов в ячейках в зависимости от плотности заселения пространства Dракошами. А также о том, как можно наблюдать что Dракоши становятся более сознательными жильцами своей вселенной.

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

image

1. Зачем нам вероятности?
2. Постановка задачи
3. Подсчет способов расстановки агентов
4. Моделирование вероятностей методом Монте-Карло
5. Что получилось?

1. Зачем нам вероятности?


Как уследить за Dракошами


Прежде чем начать запускать моделирование надо бы понять, а как мы будем «смотреть» на мир Dракош? Ведь их будут тысячи, каждый будет заниматься своими делами, куда-то идти, что-то кушать и т.д.…

Поэтому в самом начале пришлось задумался о параметрах, которые будут описывать текущее состояние моделируемой системы и показывать динамику происходящих в ней изменений. Получилось всего 24 величины, за которыми интересно проследить, помимо самого пространства с распределением агентов и пищи. Среди этих параметров есть вполне очевидные: количество агентов, количество пищи в пространстве, рождаемость и смертность… Но о них более подробно будет рассказано в третей части. Сейчас же я больше остановлюсь на специальной группе параметров, которые позволят оценить степень «осознанности» действий агентов. Другими словами, отличить случайное поведение от «осознанного».

Напомню, агенты могут совершать шесть разных действий: Pass, Step 1, Step 2, Eat, Attack, Sex. Для каждого из них можно рассчитывать долю агентов, которые его совершает в течении каждого такта времени: $f_{Pass}$, $f_{Stp 1}$, $f_{Stp 2}$, $f_{Eat}$, $f_{Att}$, $f_{Sex}$. Эти шесть параметров позволят понять, что делают Dракоши.

Все действия кроме Pass могут быть успешными или неуспешными, в зависимости от сложившейся ситуации в текущей и в соседних ячейках, а также от желаний соседа по занимаемой ячейке. Для отслеживания успешности используются четыре параметра: $L_{Stp}$, $L_{Eat}$, $L_{Att}$, $L_{Sex}$. Они показывают какая доля из совершенных действий была успешной. При этом оба действия Step 1 и Step 2 учитываются в одном параметре $L_{Stp}$.

Индексы осознанности или зачем нам вероятности


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

Так успешность действий Step 1/2 определяется наличием в ячейке, которую агент выбрал для перемещения, свободного места (всего в ячейках есть по два места для агентов). И будет зависеть как от правильности выбора агента, так и от количества агентов в пространстве. Текущее количество агентов удобно отслеживать с помощью параметра $RDr$ – плотность агентов в пространстве, которая равна отношению количества агентов к количеству ячеек в пространстве. Когда агентов в пространстве мало, плотность близка к нулю и вероятность того, что в выбранной агентом ячейке будет свободное место близка к единице. Поэтому несмотря на наличие или отсутствие «осознанности» у агентов при низкой плотности агентов ($RDr≈0$) успешность действий Step 1/2 будет высокой ($L_{Stp}≈1$). По мере увеличения количества агентов вероятность обнаружить в выбранной случайно ячейке свободное место снижается вместе с успешностью $L_{Stp}$, если только выбор ячейки для совершения шага не будет делаться осознано. Когда количество агентов будет близко к максимальной емкости пространства ($RDr≈2$) вероятность найти свободное место станет близкой к нулю ($L_{Stp}≈0$). Можно ввести параметр осознанности действий Step 1/2, который рассчитаем, как разность между текущей успешностью $L_{Stp}$ и вероятностью обнаружить свободное место в ячейке $(1-W_2)$:

$I_{Stp}=L_{Stp}-(1-W_2)$


Здесь $W_2$ – обозначает вероятность присутствия в выбранной ячейке двух агентов.

Таким же образом успешность действий Attack и Sex будет зависеть от вероятности ($W_1$) что агент обнаружит в ячейке, которую он занимает, еще одного агента.

Осознанность совершения атаки можно определить как превышение успешности $L_{Att}$ над вероятностью $W_1$:

$I_{Att}=L_{Att}-W_1$


А вот успешность $L_{Sex}$ зависит не только от вероятности найти партнера в ячейке, но и от вероятности что партнер тоже совершает в данный момент действие Sex. Но даже если партнер оказался сговорчивым, для успешного полового акта необходимо что бы в окружающем пространстве нашлось место для малыша. Поэтому вероятность случайного успеха будет произведением вероятностей трех независимых событий: $W_1$ – в ячейке есть потенциальный партнер, $f_{Sex}$ – партнер так же совершает действие Sex, $1-W_2$ – в ячейке, которую выбрала мать для размещения ребенка, есть свободное место. Поэтому осознанность брачного поведения предлагаю определить таким образом:

$I_{Sex}=L_{Sex}-f_{Sex}W_1(1-W_2)$



2. Постановка задачи


Для вычисления введенных выше индексов осознанности необходимо научится рассчитывать две вероятности: $W_1$ и $W_2$. Очевидно, что они зависят от плотности агентов в пространстве, а также от равномерности их распределения по пространству. Но для упрощения расчетов в дальнейшем будем предполагать равномерное распределения агентов по пространству. То есть будем считать, что эти вероятности не зависят от положения агента в пространстве.

Сформулируем задачу следующим образом.
Имеется пространство из $S$ ячеек. В одной ячейке может находится не более $e$ агентов. В этих ячейках располагаются $N$ агентов ($N\leqslant eS$). Необходимо определить вероятность того, что в выбранной ячейке с одним агентом будет находится еще один агент – $W_1$. И вероятность того, что в выбранной ячейке будет $e$ агентов – $W_2$.

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

Во втором подходе использовалось понятие частотной вероятности. Расчеты выполнялись с применением метода Монте-Карло. Вероятно, есть и другие способы вычислить эти вероятности, и возможно даже более простые. Буду благодарен за наводку на таковые.

3. Подсчет способов расстановки агентов


Вывод рекуррентной формулы


Предложенный ниже подход к расчету вероятностей позволяет решить поставленную задачу даже в более общей постановке, когда вместимости ячеек различны. Зададим вместимости ячеек пространства в виде последовательности $e_i$, где $i$ – номер ячейки (от $S$ до $1$). Для удобства вычислений перенумеруем ячейки таким образом что бы у выбранной для определения вероятностей ячейки был номер $S$.

Поставленную задачу можно свести к подсчёту количества способов расставить агентов по ячейкам пространства. Тогда вероятность $W_1$ – это отношение количества способов расстановки агентов, при которых в ячейку $S$ попадает хотя бы один агент, к полному количеству способов расставить $N-1$ агентов по $S$ ячейкам пространства. Для мира Dракош вместимость всех ячеек равна двум. Но при расчете $W_1$ вместимость ячейки $S$ принимается $e_S=1$, т.к. изначально известно, что в этой ячейке уже есть один агент (см. определение вероятности $W_1$ выше).

Вероятность $W_2$ – это отношение количества способов расстановки агентов, при которых в ячейку $S$ попадает предельное для этой ячейки количество агентов – $e_S$, к полному количеству способов расставить $N$ агентов по $S$ ячейкам пространства.

Обозначим через $F_i(n)$ количество способов расставить $n$ агентов по ячейкам с $i$ по первую. Так же для вычислений нам понадобится знать количество свободных мест в ячейках с $i$ по первую – $\sigma_i$:

$\sigma_i=\sum_{k=i}^{1}e_k$


Рассмотрим процесс, при котором мы последовательно перебираем все ячейки пространства (от $S$ до $1$) и на каждом шаге помещаем в очередную ячейку – $i$ некоторое количество агентов – $j$, пока не расставим все $n$ неразмещенных агентов или не исчерпаем все ячейки. В ячейку $i$ можно поместить от $0$ до $e_i$ агентов. Но если агентов осталось меньше максимальной вместимости ячейки, то в текущую ячейку можем поместить не более $n$ агентов. Таким образом $j=0..\min(n,e_i)$.

Количество способов заполнить ячейку агентами – это сочетание без повторений из $n$ по $j$:

$C_n^j=\frac{n!}{j!(n-j)!}$


Поместив $j$ агентов в ячейку $i$, останется еще $n-j$ агентов для расстановки в ячейки с $i-1$ по $1$. И каждому $j$ будет соответствовать $F_{i-1}(n-j)$ способов расстановки оставшихся агентов. Произведение $C_n^jF_{i-1}(n-j)$ – даёт нам количество способов разместить $j$ агентов в ячейке $i$ и $n-j$ агентов по остальным ячейкам. Суммируя по всем допустимым значениям $j$, получим искомое количество способов расстановки:

$F_i(n)=\sum_{j=0}^{\min(n,e_i)}{C_n^jF_{i-1}(n-j)}$


Для того что бы воспользоваться полученной рекуррентной формулой необходимо задать граничные условия:

  • Если $n=0$, т.е. у нас не осталось агентов которых можно поместить в ячейку, тогда у нас есть всего один способ выполнить расстановку агентов – никого не помещать в ячейку:

    $F_i(0)=1$

  • Для последней ячейки $i=1$ так же есть только один способ расставить агентов – оставить все $n$ агентов в ней, но только если им хватит места:

    $F_1(n)=1,\qquad n\leq\sigma_1$

  • Если места не хватает для размещения $n$ агентов в ячейках с $i$ по $1$, то способов заполнить ячейку нету:

    $F_i(n)=0,\qquad n>\sigma_i$


Теперь можем записать выражения для вероятности $W_1$:

$W_1=\frac{F'_S(N-1)}{F_S(N-1)}$


где $F’_S(N-1)$ – это количество способов расстановки агентов, при которых в ячейку $S$ попадает один и более агентов:

$F'_S(N-1)=\sum_{j=1}^{\min(N-1,e_S)}{C_{N-1}^jF_{S-1}(N-1-j)}$


Обратите внимание, здесь начинаем суммировать с единицы, т.к. нас интересуют те расстановки, при которых в ячейку $S$ попадает хотя бы один агент.

Выражения для вероятности $W_2$:

$W_2=\frac{F''_S(N)}{F_S(N)}$


где $F’’_S(N)$ – это количество способов расстановки агентов, при которых в ячейку $S$ попадает предельное для этой ячейки количество агентов:

$F''_S(N)=C_N^{e_S}F_{S-1}(N-e_S)$


Пример вычисления вероятностей для S=5
Рассмотрим пример вычисления количества способов расстановки агентов по пяти ячейкам ($S=5$), вместимость которых равна двум ($e=2$).

Результаты вычислений сведены в таблицу:



Упрощение и переход приведенному виду


Количество расстановок очень быстро растет и уже для $S=50$ чисел двойной точности не хватит. Даже если воспользоваться пакетами символьных вычислений с произвольной точностью (например, Symbolic Math Toolbox и Variable-precision arithmetic в среде MatLab) останется проблема переполнения стека вызовов рекуррентной функции при больших размерах пространства. Поэтому полученные выше выражения необходимо преобразовать. А еще лучше вывести формулу, с помощью которой сразу вычисляется $F_i(n)$ без необходимости считать предыдущие значения.

Первое что можно сделать это вынести из-под знака суммы $n!$.

Но прежде введем обозначение приведенного количества расстановок $n$ агентов по $i$ ячейкам, которое в $n!$ раз меньше:

$f_{i,n}=\sum_{j=0}^{\min(n,e_i)}{\frac{1}{j!(n-j)!}F_{i-1}(n-j)}=\sum_{j=0}^{\min(n,e_i)}{\frac{1}{j!}f_{i-1,n-j}}$


Маленький алгебраический трюк

$F_i(n)=\sum_{j=0}^{\min(n,e_i)}{\frac{n!}{j!(n-j)!}F_{i-1}(n-j)}=n!\sum_{j=0}^{\min(n,e_i)}{\frac{1}{j!(n-j)!}F_{i-1}(n-j)}$


Полное количество способов расстановки теперь будет:

$F_i(n)=n!f_{i,n}$


Аналогично для $F_{i-1}(n-j)=(n-j)!f_{i-1,n-j}$, и после подстановки в $f_{i,n}$ получим:

$f_{i,n}=\sum_{j=0}^{\min(n,e_i)}{\frac{1}{j!}f_{i-1,n-j}}$



И можем сократить выражения для вероятностей на факториал:

$W_1=\frac{F'_S(N-1)}{F_S(N-1)}=\frac{(N-1)!f'_{S,N-1}}{(N-1)!f_{S,N-1}}=\frac{f'_{S,N-1}}{f_{S,N-1}} \\ W_2=\frac{F''_S(N)}{F_S(N)}=\frac{N!f''_{S,N}}{N!f_{S,N}}=\frac{f''_{S,N}}{f_{S,N}}$


Где

$f'_{i,n}=\sum_{j=1}^{\min(n,e_i)}{\frac{1}{j!}f_{i-1,n-j}} \\ f''_{i,n}=\frac{1}{e_i!}f_{i-1,n-e_i}$


После введения новых обозначения требуется переписать выражения для граничных условий. Для дальнейших выкладок удобнее представить граничные условия в следующем виде:

  • $f_{0,0}=1$, при $i=0$, $n=0$;
  • $f_{i,0}=f_{i-1,0}$, при $i>0$, $n=0$;
  • $f_{i,1}=f_{i-1,1}+f_{i-1,0}$, при $i>0$, $n=1$;
  • $f_{i,n}=0$, при $n>\sigma_i$.

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

Пример вычисления вероятностей для S = 5 в приведенном виде
Самое время проверить что ничего не испорчено нашими упрощениями.


Следующим шагом предлагаю избавится аж от трех видов выражений для приведенного количества расстановок: $f_{i,n}$, $f’_{i,n}$, $f’’_{i,n}$. Для этого достаточно выразить вероятности через $f_{S-1,n}$. Тут нам придется вспомнить что вместимость ячеек в мире Dракош везде $e_i=2$, кроме выбранной ячейки в которой одно место уже занято и поэтому $e_S=1$. Дальнейшие решение поставленной задачи в общей постановке (когда вместимость ячеек произвольна) существенно усложнит выкладки в и так затянувшемся погружении в мир комбинаторики.

Избавление
Распишем $f_{S,N-1}$, $f_{S,N}$, $f’_{S,N-1}$, $f’’_{S,N}$ согласно их основным формулам и подставим в выражения для $W_1$ и $W_2$:

$W_1=\frac{F'_S(N-1)}{F_S(N-1)}=\frac{f'_{S,N-1}}{f_{S,N-1}}=\frac{f_{S-1,N-2}}{f_{S-1,N-1}+f_{S-1,N-2}}=\frac{1}{1+\frac{f_{S-1,N-1}}{f_{S-1,N-2}}}$


$W_2=\frac{F''_S(N)}{F_S(N)}=\frac{f''_{S,N}}{f_{S,N}}=\frac{\frac{1}{2}f_{S-1,N-2}}{f_{S-1,N}+f_{S-1,N-1}+\frac{1}{2}f_{S-1,N-2}}=\frac{1}{1+2\frac{f_{S-1,N}}{f_{S-1,N-2}}+2\frac{f_{S-1,N-1}}{f_{S-1,N-2}}}$



Готовые выражения для вероятностей теперь зависят только от приведенного количества способов расстановки агентов одного вида – $f_{i,n}$:

$W_1=\frac{1}{1+\frac{f_{S-1,N-1}}{f_{S-1,N-2}}}$


$W_2=\frac{1}{1+2\frac{f_{S-1,N}}{f_{S-1,N-2}}+2\frac{f_{S-1,N-1}}{f_{S-1,N-2}}}$


Построение производящей функции


Осталось научится рассчитывать элементы последовательности $f_{i,n}$ для произвольно больших значений $i$ и $n$. Порывшись в учебниках по комбинаторике я нашёл довольно популярный и мощный метод для решения всяческих комбинаторных вопросов – метод производящих функций.

Хорошим пособием для начинающих комбинаторов по производящим функциям можно назвать учебник Ландо С.К. [3]. Еще один ресурс типа справочник по производящим функциям. И еще вот этот пост на хабре. Но наиболее полезными для меня были видео лекции Игоря Клейна по производящим функциям, выложенные тут.

Но в перечисленных ресурсах в основном речь ведется о производящих функциях для одномерных последовательностей, с одним индексом. А наша последовательность двумерная, с двумя индексами. Производящие функции от двух переменных лишь упоминаются. Поэтому пришлось рискнуть и попробовать вывести производящую функцию для нашей двумерной последовательности на основе одномерных примеров.

Пару слов о том, что такое производящая функция последовательности
Это формальная замена последовательности $a_0$, $a_1$, $a_2$, $a_3$, … степенным рядом:

$\langle a_0,a_1,a_2,a_3,...\rangle \Leftrightarrow G(x)=\sum_{i=0}^\infty{a_ix^i=a_0+a_1x+a_2x^2+a_3x^3+\dots}$


Идея метода в том, что получив каким-то из способов выражение для производящей функции $G(x)$ можно затем разложить его в степенной ряд $\sum_{i=0}^\infty{a_ix^i}$ и тем самым, например, получить выражение для общего члена последовательности $a_i$ в замкнутом виде (не требующем вычисления предыдущих членов последовательности). Это как раз то, что нам требуется.

Суммируя то, что имеется на данный момент по нашей последовательности приведенных количеств расстановок запишем рекуррентную формулу для $f_{i,n}$ совместно с граничными условиями.

$\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.$


Наша последовательность $f_{i,n}$ зависит от двух индексов: $i$ – количество ячеек в пространстве для расстановки агентов, $i=0\dots S$; $n$ – количество агентов для расстановки в ячейках, $n=0\dots N$. Поэтому и производящая функция ей нужна двумерная – $\mathscr{F}(x,y)$.

$\mathscr{F}(x,y)={ f_{0,0}x^0y^0+f_{0,1}x^0y^1+f_{0,2}x^0y^2+f_{0,3}x^0y^3+f_{0,4}x^0y^4+\dotsb \\ f_{1,0}x^1y^0+f_{1,1}x^1y^1+f_{1,2}x^1y^2+f_{1,3}x^1y^3+f_{1,4}x^1y^4+\dotsb \\ f_{2,0}x^2y^0+f_{2,1}x^2y^1+f_{2,2}x^2y^2+f_{2,3}x^2y^3+f_{2,4}x^2y^4+\dotsb \\ \ldots }$


Используя специальный технический примем можно из рекуррентного выражения для элементов последовательности получить выражение для производящей функции. У меня получилось вот это:

$\mathscr{F}(x,y)=\frac{1}{1-x \left(1+y+\frac{1}{2}y^2 \right)}$


Тот самый технический прием
Суть приема в том, чтобы выписать первые элементы последовательности в столбик согласно рекуррентной формуле и умножить каждую строку на соответствующий множитель $x^iy^n$.

$\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}$


Просуммируем полученные выражения в столбик следующим образом. Отдельно суммируем левый от равно столбец и получим собственно производящую функцию — $\mathscr{F}(x,y)$

Суммируя синие слагаемые справа от равно, получим:

$\color{Blue}{1+f_{0,0}x^1y^0+f_{1,0}x^2y^0+f_{1,1}x^2y^1+f_{1,2}x^2y^2+\cdots}$


Замете, что начиная со второго слагаемого можно вынести за скобки множитель $x$ и там останется наша производящая функция:

$\color{Blue}{1+x \left( f_{0,0}x^0y^0+f_{1,0}x^1y^0+f_{1,1}x^1y^1+f_{1,2}x^1y^2+\cdots \right) = 1 + x \mathscr{F}(x,y)}$


Суммируя зеленые слагаемые, получим:

$\color{Green}{ \begin{array}{c} f_{0,0}x^1y^1+f_{1,0}x^2y^1+f_{1,1}x^2y^2+f_{1,2}x^2y^3+\cdots = \\ = xy \left( f_{0,0}x^0y^0+f_{1,0}x^1y^0+f_{1,1}x^1y^1+f_{1,2}x^1y^2+\cdots \right) = \\ = xy \mathscr{F}(x,y) \end{array} }$


Суммируя красные слагаемые, получим:

$\color{Red}{ \begin{array}{c} \frac{1}{2}f_{0,0}x^1y^2+\frac{1}{2}f_{1,0}x^2y^2+\frac{1}{2}f_{1,1}x^2y^3+\frac{1}{2}f_{1,2}x^2y^4+\cdots = \\ = \frac{1}{2}xy^2 \left( f_{0,0}x^0y^0+f_{1,0}x^1y^0+f_{1,1}x^1y^1+f_{1,2}x^1y^2+\cdots \right) = \\ = \frac{1}{2}xy^2 \mathscr{F}(x,y) \end{array} }$


После выполненных суммирований получим уравнение относительно искомой функции:

$\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)}$


Решая которое относительно $\mathscr{F}(x,y)$ получим выражение для производящей функции.

Теперь дело за малым, разложить полученную производящую функцию в ряд по степеням $xy$. Собственно, произвести элементы последовательности:

$\mathscr{F}(x,y)=\sum_{i=0}^{\infty}{\left( 1+y+\frac{1}{2}y^2 \right)^i x^i}=\sum_{i=0}^{\infty}{\sum_{k=0}^i{\sum_{m=0}^k{\frac{1}{2^m}C_i^kC_k^m y^{k+m}}}}$


Процесс разложения в ряд
Для начала обозначим:

$\begin{array}{c} Y=1+y+\frac{1}{2}y^2=1+z \\ z=y \left( 1+\frac{1}{2}y \right)=y(1+p) \\ p=\frac{1}{2}y \end{array}$


С учетом этих обозначений производящая функция запишется:

$\mathscr{F}(x,y)=\frac{1}{1-Yx}$


Вот тут лежит аж семь объяснений того почему $\frac{1}{1-x}$ это степенной ряд $\sum_{i=0}^{\infty}x^i=1+x+x^2+x^3+\cdots$. С учетом этого наша функция раскладывается в ряд по степеням $x$ следующим образом:

$\mathscr{F}(x,y)=\sum_{i=0}^{\infty}Y^ix^i$


Теперь разберемся с $Y^i$.

$Y^i=(1+z)^i$


А это уже бином Ньютона:

$(1+z)^i=\sum_{k=0}^iC_i^kz^k$


Аналогично для $z^k$:

$z^k=y^k(1+p)^k=y^k\sum_{m=0}^k{C_k^mp^m}=y^k\sum_{m=0}^k{\frac{1}{2^m}C_k^my^m}=\sum_{m=0}^k{\frac{1}{2^m}C_k^my^{k+m}}$


Осталось подставить полученное выражение для $Y^i$ в производящую функцию:

$\mathscr{F}(x,y)=\sum_{i=0}^{\infty}Y^ix^i=\sum_{i=0}^{\infty}{\sum_{k=0}^i{\sum_{m=0}^k{\frac{1}{2^m}C_i^kC_k^m y^{k+m}}}}$



Теперь все что нам нужно это выделить коэффициент при $x^iy^n$, это и будет искомое выражение для элемента последовательности $f_{i,n}$:

$[x^iy^n]\mathscr{F}(x,y)=f_{i,n}$


Для этого необходимо сгруппировать все слагаемые, у которых степень $y$ — $k+m$ равна $n$. Пределы суммирования при этом будут следующими: при $m=0$, $k = n$; при $m = k$, $k = \frac{n}{2}$.

В результате получим два выражения, одно для чётного $n$:

$\eta=\frac{n}{2},\qquad f_{i,n}=\sum_{k=\eta}^{2\eta}{\frac{1}{2^{2\eta-k}}C_i^kC_k^{2\eta-k}}$


и одно для нечетного $n$:

$\eta=\frac{n+1}{2},\qquad f_{i,n}=\sum_{k=\eta}^{2\eta-1}{\frac{1}{2^{2\eta-k-1}}C_i^kC_k^{2\eta-k-1}}$


Полученные выражения уже можно попробовать вычислить, но так как биноминальные коэффициенты, которые входят в них, всё же растут довольно быстро, для расчетов необходимо воспользовался инструментами символьных вычислений MatLab.

Пример реализации в MatLab
Вычисление $f_{i,n}$.

function out = f(i,n)
%% Символьное вычисление f(i,n)

nu = ceil(n/2); % нижний предел суммирования

syms k

if 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

Вычисление $W_1$.

function out = W1sym(ii,N)
%% Символьное вычисление вероятности обнаружить второго агента в ячейке пространства

out = 1/(1 + f(ii-1,N-1)/f(ii-1,N-2));

Вычисление $W_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) используется для установления статистических характеристик других случайных величин или для приближённой оценки некоторого значения $a$. Во втором случае придумывается некоторая случайная величина $\xi$, такая что ее мат. ожидание равно искомой величине $a$. Тогда рассчитав $N$ независимых значений случайной величины $\xi_1$, $\xi_2$, …, $\xi_N$ определяют величину $a$ как среднее арифметическое выборки $\xi_1$, …, $\xi_N$ [2].

Для установления статистических характеристик случайных величин, например, для определения вероятности $p$ случайного события $А$, можно реализовать другой вариант метода Монте-Карло. Когда проводится моделирование такой случайной величины $\delta_i$, которая равна $1$ при реализации случайного события $А$ на испытании $i$, или равна $0$ если события $А$ не наступило. При проведении $N$ независимых испытаний частота события $А$, равная $\frac{1}{N}\sum_{i=1}^N\delta_i$, стремится к искомой вероятности $p$ по мере роста количества испытаний $N$ [3]. Если для моделирования случайной величины $\delta_i$ в модели полностью или частично воспроизводится изучаемый естественный процесс, то такой подход так же называют – имитационным моделированием.

На сегодня метод имитационного моделирования развился в самостоятельное направление методов исследований. Но для решения поставленной в разделе 2 задачи использовался подход, который с одной стороны является примером реализации методов Монтер-Карло, с другой это пример имитационного моделирования, т.к. использовался «естественный» процесс изменения положения агентов в пространстве.

Обзор алгоритма


В рамках метода Монте-Карло поставленная в разделе 2 задача сводится к оценке частоты реализации соответствующих событий. Так вероятность $W_1$ – это частота обнаружения агентами соседей по ячейке. А вероятность $W_2$ – это частота обнаружения ячеек с двумя агентами.

Будем рассматривать процесс заполнения $S$ ячеек $N$ агентами, так что в одной ячейке может быть не более $e=2$ агентов. На каждом шаге будем пересаживать агентов в новые ячейки, и если при попытке пересадить агента окажется что в выбранной ячейке нет свободного места, то оставляем агента на старом месте.



После каждого этапа пересаживаний будем подсчитывать количество $C_2$ полностью заполненных ячеек (с двумя агентами).

$W_1=\frac{2C_2}{N}$


Частота обнаружения ячеек с двумя агентами (вероятность $W_2$) это просто:

$W_2=\frac{C_2}{S}$


Можно предложить несколько вариантов реализации этапа пересаживания:

  • Scenario 1 – почти полная имитация: для пересаживания случайно выбираются 50% агентов и каждый из них пересаживается в одну из 18 соседних клеток выбранных случайно;
  • Scenario 2 – пересаживаются все агенты и каждый из них пересаживается в одну из 18 соседних клеток выбранную случайно;
  • Scenario 3 – для пересаживания случайно выбираются 50% агентов и каждый из них пересаживается в случайную ячейку пространства;
  • Scenario 4 – пересаживаются все агенты и каждый из них пересаживается в случайную ячейку пространства;

Scenario 1 назван почти полной имитацией процессов в настоящем мире Dракош, так как доля агентов совершающих действия Step 1/2 различна для каждого такта времени и может иметь любое значение а не только 50%. Ниже мы посмотрим, как результаты вычислений зависят от сценария пересаживания агентов.

Процесс начинаем со случайной расстановки агентов по ячейкам. Затем в цикле $n$ раз повторяются два этапа:

  1. Пересаживаем агентов согласно выбранного сценария;
  2. Определяем количество заполненных ячеек – $С_2$ и соответствующие значения $W_1$ и $W_2$.

Возникает вопрос: сколько же раз необходимо повторить цикл? С ростом $n$ неопределенность вычислений убывает пропорционально $1/\sqrt{n}$. Т.е. для увеличения точности в 10 раз необходимо увеличить количество испытаний в 100 раз. Метод Монте-Карло похож на натурные физические эксперименты в том плане что точность результатов можно определить только после получения результата эксперимента. Точнее случайную составляющую неопределенности. И так же как в натурном эксперименте ею можно управлять изменяя количество измерений.

На графике ниже показана зависимость вероятности $W_1$ для $S=10000$ и $RDr=1$ ($N=10000$) от количества испытаний. Для каждой точки показаны планки неопределенности, которая рассчитывалась на основе стандартного отклонения полученных выборок и коэффициента охвата $k=3$ (надежность 99.7%).

$u_{W_1}=k\sqrt{\frac{\sum_{i=1}^n{(w_{1_i}-\overline{w_1}})^2}{n(n-1)}}$


Горизонтальной чертой отмечено точное значение, рассчитанное способом, описанным в разделе 3.



По результатам расчетов получены следующие значения неопределенности (надежность 99.7%):

  • Для $n = 100$ неопределенность составила 0.0014 (0.24%);
  • Для $n = 1000$ – 0.00047 (0.080%);
  • Для $n = 10000$ – 0.00015 (0.025%);
  • Для $n = 50000$ – 0.000066 (0.011%);

Т.к. время вычислений растет пропорционально количеству испытаний (для $n = 50000$ составило около 6.5 часов), а желания сильно оптимизировать код у меня не было, решено было для последующих экспериментов использовать по 1000 испытаний.

Результаты Монте-Карло


Давайте посмотрим, как зависит рассчитываемые значения вероятностей от используемого сценария пересаживания агентов. Ниже показаны результаты вычислений для разных сценариев, при $S = 10000$ и $RDr = 1$ ($N = 10000$) и количестве испытаний по $n=1000$.



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

5. Что получилось?


На графике ниже показаны зависимость вероятностей $W_1$ и $W_2$ от плотности агентов в пространстве. Синяя и красная области — это массивы данных, насчитанные методом Монте-Карло (по 1000 итераций для каждого значения $N$). Черные линии — это точные значения, вычисленные методом, описанным в разделе 3.



Совпадение результатов, полученных разными методами, позволяет надеяться, что серьезных ошибок не было допущено.

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

$\begin{array}{1} I_{Stp}=L_{Stp}-(1-W_2) \\ I_{Att}=L_{Att}-W_1 \\ I_{Sex}=L_{Sex}-f_{Sex}W_1(1-W_2) \end{array}$


Для однозначности примем что частота действия Sex – $f_{Sex}=1$ (т.е. если партнер был замечен в ячейке, то он принимает участие в размножении 18+ content).



В дальнейшем мы сравним полученные результаты с успешностью Тупиков (раса агентов без нейронной сети) и полноценных Dракош и посмотрим, как будет расти успешность в результате развития нейронной сети агентов.

  1. Ландо С.К. Лекции о производящих функциях. 2007. 144 c.
  2. Соболь И.М. Численные методы Монте-Карло. 1973. 312 с.
  3. Бусленко Н.П., Шрейдер Ю.А. Метод статистических испытаний (Монте-Карло) и его реализация на цифровых вычислительных машинах. 1961. 226 с.
Tags:
Hubs:
+3
Comments0
16
Karma
0
Rating
@TRSloz

Теплофизика

Send message

Articles

Show the best of all time
Unlock dark mode

Your account

Sections

Information

Services

Support
© 2006–2024, Habr