CC BY
33
УДК 519.711.3
СТАТИСТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ПОПУЛЯЦИЙ ВЗАИМОДЕЙСТВУЮЩИХ ЧАСТИЦ С ПРОИЗВОЛЬНЫМ РАСПРЕДЕЛЕНИЕМ ВРЕМЕНИ
ЖИЗНИ
Б.Ю. Пичугин, И.В. Перцев
We consider a stochastic model of interacting populations of the particles with arbitrary distributions of the life spans. The Monte-Carlo approach for the simulation the population dynamics is presented here.
Введение
При вероятностном описании динамики популяций широко применяются марковские случайные процессы как с дискретным, так и с непрерывным временем. Один из наиболее распространенных подходов опирается на случайный процесс рождения и гибели. Этот подход использует систему дифференциальных уравнений Колмогорова, решение которой в некоторых случаях можно выразить в явной форме [1]. Другой подход опирается на ветвящиеся случайные процессы, которые позволяют описать размножение и превращение частиц, а также их взаимодействие [2], [3]. В определенных случаях нахождение вероятностных характеристик численностей популяций является затруднительным или невозможным из-за сложности используемых уравнений. В этих случаях для изучения популяций можно применить метод Монте-Карло [4], [5].
При изучении динамики популяций в задачах демографии, экологии, эпидемиологии и др. приходится рассматривать взаимодействующие популяции частиц, в которых важную роль играет как общая численность частиц, так и их возрастной состав. Если предположить, что продолжительность времени жизни частиц описывается экспоненциальным распределением, то для построения модели можно использовать случайный процесс рождения и гибели, либо ветвящийся случайный процесс с взаимодействием частиц. Если же распределение времени жизни частиц отлично от экспоненциального, то можно применять аппарат полумарковеких случайных процессов [6]. Вместе с тем, если интенсивности взаимодействий зависят от численностей популяций сложным образом, а распределение времени жизни частиц не экспоненциально, то применение указанных выше подходов является, по видимому, невозможным.
© 2001 Б.Ю. Пичугин, Н.В. Перцев
E-mail: pichugin@lib.omskreg.ru, pertsev@univer.omsk.su Омский государственный университет
Предлагаемая имитационная модель предназначена для исследования динамики развития популяций взаимодействующих частиц с произвольным распределением времени жизни. При реализации этой модели на ЭВМ использовались специальным образом оптимизированные структуры представления данных, позволяющие сократить время вычислений и объем используемой памяти,
1. Описание модели
1.1. Популяции
Рассмотрим набор популяций Ai, А2,..., Am. Под популяцией будем понимать класс частиц, объединенных по какому-либо признаку, причем каждая участвующая в модели частица может входить только в одну популяцию. Будем считать, что все существующие в один момент времени частицы отдельно взятой популяции равновероятно принимают участие во взаимодействиях. Также для частиц одной популяции зафиксируем распределения на интервале (0; +оо) времени их жизни: Li - для частиц, рожденных в момент времени t > 0, и L° - для частиц, уже существующих в момент времени / = 0. 1 ^ ^ ///. Примем,
что распределения L, и L° не изменяются со временем и не учитывают влияния других частиц и взаимодействий, в которые частица будет вступать в процессе своего жизненного цикла. Кроме того, полагаем, что Д(0; dt) —>• 0 и L°(0; dt) —>• 0 при dt —>• 0, 1 ^ % С т.
В момент I> 0 рождения частицы х популяции Д (далее будем писать х Е А^ закрепим за ней время ее смерти вследствие старения, d{x) = tx+£x, где распределение случайной величины £х совпадает с 1 ^ ^ ш. Для частиц,
существующих в момент времени / = 0. полагаем <1{х) = £х, где случайная величина £х имеет распределение L°. Случайные величины £х и £у для любых двух частиц х ъ у считаем независимыми.
Также будем считать, что каждая частица i 6 Aj, 1 Д ^ т, в момент своей естественной смерти d(x) уничтожает и порождает rjij частиц популяции Aj, 1 С j ^ т. Так как уничтожение и рождение происходят в один и тот же момент времени, то для определенности примем, что уничтожение частиц предшествует рождению, любая частица популяции Aj может быть уничтожена с равной вероятностью, и если популяция Aj перед уничтожением содержала менее чем £ij частиц, то уничтожаются все частицы. Здесь = (бъ Сг2; ■ ■ ■, 6™) и di = {dm dvi-> • • • ■> dim), 1 ^ г ^ m, - независимые целочисленные случайные вектора с неотрицательными компонентами, конечными множествами возможных значений и постоянными распределениями. Пару этих векторов будем называть результатом смерти вследствие старения.
Заметим, что если частица х Е .1,. 1 ^ ^ ///. гибнет не в результате ста-
рения, то есть раньше предписанного времени d(x), то такая гибель не будет сопровождаться рождением новых и уничтожением существующих частиц.
Численность популяции Д в момент времени t будем обозначать через Xi(t), а численности всех популяций через вектор x(t) = (xi(t), аДД - - -, xm(t)). И, наконец, каждой популяции сопоставим уникальное имя (конечный набор сим-
волов), которое будем использовать для обращения к этой популяции в программной реализации модели,
1.2. Взаимодействия
Под взаимодействием будем понимать любой процесс, приводящий к изменению численное ! и популяций, для которого распределение времени между этими изменениями носит экспоненциальный характер и зависит только от численности популяций.
Будем считать, что за период своей жизни частицы могут участвовать в
и различных типах взаимодействий l\. l-j./„. Каждый тип взаимодействий
Ik-, 1 ^ к ^ п, будем описывать интенсивностью - неотрицательной функцией qk(x), х G Zm, и результатом - парой независимых целочисленных случайных векторов 7fc = (Tfeb Tfe2, ■ ■ ■, Tfem) и Рк = (Д ъ Д 2, ■ ■ ■, А™) с неотрицательными компонентами, конечными множествами возможных значений и постоянными распределениями.
Пусть за интервал времени [t, t+dt) нет гибели частиц в результате старения. Тогда вероятность осуществления за это время ровно одного взаимодействия будем полагать равной Q(x(t))dt + o(dt), где Q{x{t)) = Ylk=i Qk(x(t))- Вероятность возникновения за это время более одного взаимодействия примем равной o(dt).
Определим случайную величину ф(1) как продолжительность времени до первого осуществления одного из взаимодействий начиная с момента t. Предположим, что за интервал времени [t,t + ф(1)) численности популяций не менялись, а в момент времени t + ф(1) произошло взаимодействие. Тогда случайная величина ф(1) имеет экспоненциальное распределение
Р{ф(ф > s} = e-Qix(t))s, О 0.
Номер типа осуществившегося в момент времени t + ф(1) взаимодействия обозначим через к(1) и примем, что к(1) есть случайная величина с распределением
В результате осуществления взаимодействия типа 1др в каждой популяции Д, 1 ^ г ^ т,, в момент времени t + ф(1) будет уничтожено yK{t)i и родится /Зкдр частиц. Здесь так же, как мы определяли результат смерти вследствие старения, предполагается, что уничтожение частиц предшествует их рождению; любая частица популяции А{ может быть уничтожена с равной вероятностью, и если популяция А{ перед уничтожением содержала менее, чем yKp)i частиц, то уничтожаются все частицы,
2. Описание алгоритма моделирования
2.1. Входные параметры имитационной модели
Как видно из описания модели, к ее параметрам относятся: количество популяций ///: количество типов взаимодействий и: распределения времени жизни
частиц Li и L°, 1 С i ^ ш; интенсивности взаимодействий цДт), 1 С k С Щ распределения результатов смерти вследствие старения (Д;?Д, 1 ^ г ^ т, и результатов взаимодействий (ук',/3к), 1 ^ к С Щ начальные численности популяций х° = аДО), 1 С i С, т. Параметры тип вычисляются автоматически при описании популяций и типов взаимодействий,
2.1.1. Распределения времени жизни частиц
На данном этапе разработки имитационной модели проблема удобного представления произвольного распределения времени жизни частиц не имеет эффективного решения. Таким образом, при описании конкретной модели пользователю необходимо для каждой популяции указать пару функций, которые будут моделировать случайные величины £х с распределениями Li и L°. В дальнейшем планируется создать динамическую библиотеку распределений, редактирование которой будет доступно пользователю,
2.1.2. Интенсивности взаимодействий
Интенсивность взаимодействия - это неотрицательная функция целочисленных аргументов. Данная функция должна обладать содержательным смыслом, отражающим специфику изучаемых объектов моделирования. Поэтому в качестве интенсивностей будем брать формулы из класса Q, построенного по перечисленным ниже правилам,
1. q(xi,x2, ■ ■ ■ ,хт) = 1 е Q.
Эта интенсивность описывает взаимодействие, не зависящее от состояний модели, что соответствует регулярному внешнему воздействию на модель, например, притоку частиц,
2. q(xi,x2, ■ ■ ■, хт) = Xi е Q, 1 ^ i С га.
Такая интенсивность соответствует взаимодействию, которое может быть инициировано каждой частицей популяции Д самостоятельно. Примером такого взаимодействия может служить производство частицами популяции Д частиц популяции Д,
3. (/ е Q ^ к * q е Q. к е II?/ .
Здесь * соответствует операции умножения, а взаимодействие, описываемое интенсивностью к * q, осуществляется в к раз чаще взаимодействия, описываемого интенсивностью q.
4. qi & Q, q-< e Q=> tp — q-> L Q,
Интенсивность щ + q2 описывает взаимодействие, возникающее при осуществлении любого из взаимодействий, описываемых интенсивностями щ и q2.
5. qi е Q, q2 е Q => qi * q2 е Q.
Здесь * - это операция умножения, и в случае, когда щ = ащ q2 = хj.
i Ф j, функция qi * q2 = х\ * описывает парное взаимодействие друг с другом частиц из популяций А{ и . 1. 1 sC i.j С /и. Последовательно применяя эту операцию, можно описывать взаимодействия трех, четырех и более частиц различных популяций,
6, q Е Q =>■ com,b(q,k) = О/' е Q. к G N (построена без использования правила 3),
Если q = Xjt, то интенсивность comb(q, к) описывает взаимодействие между собой к частиц одной популяции Ai, 1 ^ г *6, т.
7. Ч\ е Q. 42 е Q =>• mm(gi,g2) G Q, mm(gi,§2) = 42^ mm(?i,?2) G Q.
При помощи функций min и imin можно описывать поведение частиц в условиях потребления ресурса. Например, пусть популяция Аг - потребитель, а популяция А2 выступает в роли ресурса. Тогда взаимодействие с интенсивностью qi(x) = х\ * min(x2, к), к е N, осуществляется в момент потребления частицы-ресурса популяции .С частицей-потребителем популяции Ai. Если .г2 > к. то частицы-потребители не испытывают недостатка в ресурсе, и потребление протекает с интенсивностью, зависящей только от свойств частиц-потребителей, В противном случае частицы-потребители испытывают дефицит ресурса и интенсивность потребления снижается пропорционально объему ресурса. При помощи интенсивности 42{х) = Xi * imin(x2,k) можно описать взаимодействие, которое возникает при дефиците ресурса (например, это может быть гибель частиц-потребителей), В этом случае с уменьшением Х2 интенсивность этого взаимодействия будет возрастать.
Индукцией по длине формулы можно показать, что каждая функция q е Q определена и неотрицательна на множестве Z™,
Пусть интенсивность взаимодействий типа I/.. 1 б. к б. н. представляется в виде формулы q е Q, Тогда, заменяя в формуле q все символы ащ 1 ^ г *6, т, на имена соответствующих популяций, получим строку символов, которая и подается на вход программной реализации модели при описании типа взаимодействий
2.1.3. Распределения результатов смерти вследствие старения и результатов взаимодействий
Представление распределений результатов смерти вследствие старения //,). 1 ^ г *6, ш, и результатов взаимодействий (у*,; Д), 1 к <6, п, принятое в имитационной модели, покажем на примере случайной величины ф. Остальные распределения результатов описываются аналогично.
Совокупность, состоящую из кг частиц популяции Аг, к2 частиц популяции ,12...к,„ частиц популяции Ат, будем обозначать с помощью мультииндекса
Ак = Akl * А*2 * ... * Ак™, Щ е Z+, 1 «С j ^ гп.
Пусть случайный вектор Д принимает г различных значений щ, щ, ..., ur с вероятностями pi, Р2, ■ ■ ■, рг-, соответственно, тогда распределение этого случайного вектора будем записывать следующим образом:
(Ci) = Pi* Аш + Р2* AU2 + ■ ■ ■ +pr * AUt .
Если Д = 0 с вероятностью 1, то его распределение будем записывать (Д) = 0, Далее, заменив в выражении (Д) все символы Д, 1 Д i С ш, на имена соответствующих популяций, а операцию взятия верхнего индекса - символом " " ". получим строку символов, которая и подается на вход программной реализации модели в качестве описания распределения случайной величины Д,
2.1.4. Начальные численности популяций
Начальные численности популяций х° = яД0), 1 Д г Д т, - неотрицательные целые числа, которые задают число частиц популяций в начальный момент времени t = 0 без учета их возраста. Возраст этих частиц вычисляется автоматически в соответствии с распределением Ь°.
2.2. Алгоритм моделирования
Для описания алгоритма будем использовать следующие переменные:
/ - текущее время;
t,4 - время первой гибели вследствие старения, то есть t,4 - это наименьшая из величин d(x) по всем существующим в момент времени t частицам х (см. пункт 1.1.);
// = t+ip(t) - время осуществления очередного взаимодействия (см. пункт 1.2.); Qk = 4k(x(t)), 1 Д k С п, - значения интенсивностей взаимодействий в момент времени t (см. пункт 1.2.);
Q = Q(x(t)) = Щ + Ц2 + - ■ - + <1п - суммарная интенсивность всех взаимодействий (см. пункт 1.2.).
0. Инициализация
0.1. Для всех популяций Д, 1 Д i С ш, сгенерировать х° частиц, для каждой из которых вычислить <1 (х} = £х, где распределение случайной величины £х совпадает с распределением Ь°.
0.2. Инициализировать переменные: t := 0;
t,A '■ = min{d(a:) : х G Д, 1 Д % С т};
Чк := Чк(х°), 1 Д к С щ Q '■= 4l + 42 + ■ ■ ■ + Чп;
1. Итерационная часть
1.1. Если Q = 0, то // := +оо,
иначе tj := t — ^ln(a), где а - равномерно распределенная на (0; 1) случайная величина.
1,2, Если t,4 ^ 1л, mo
1.2.1, Для всех частиц х е Д, 1 ^ г ^ т, таких, что d(x) = Д, вычислить величины и г]х, распределения которых равны Д) и (щ) соответственно,
1.2.2, Вычислить общее число уничтожаемых £, и рождаемых // частиц по формулам:
1.2.3, Уничтожить все частицы, для которых d(x) = Д,
1.2.4, Уничтожить частицы в популяциях в соответствии с вектором £ и соглашениями, принятыми в пункте 1,1,
1.2.5, Породить частицы в популяциях в соответствии с вектором г] и соглашениями, принятыми в пункте 1,1,
1.2.6, Для каждой рожденной частицы я: 6 4; определить время ее гибели в результате старения d(x) = / | + Д, где распределение случайной величины £х совпадает с Д, 1 ^ г ^ т,.
1.2.7, Изменить текущее время / := /
1.3, Иначе (если Д > tj)
1.3.1, Вычислить значение к случайной величины к(1) в соответствии с распределением, определенным в пункте 1,2,
1.3.2, Вычислить величины у и /3, распределения которых равны (ук) и (/Зк) соответственно,
1.3.3, Уничтожить частицы в популяциях в соответствии с вектором у и соглашениями, принятыми в пункте 1,2,
1.3.4, Породить частицы в популяциях в соответствии с вектором /3 и соглашениями, принятыми в пункте 1,2,
1.3.5, Для каждой рожденной частицы х е Д определить время ее гибели в результате старения <1(х) = tj + £х, где распределение случайной величины £х совпадает с Д, 1 ^ г ^ т,.
1.3.6, Изменить текущее время / := //.
1.4, Пересчитать время до первой гибели в результате старения
Д : = min{d(a:) : х е Д, 1 ^ ^ т}.
1.5, Пересчитать интенсивности ц*. := цДаД)), 1 ^ к ^ п.
1.6, Пересчитать суммарную интенсивность Q := щ + q2 + ... + qn.
1.7, Перейти к следующей итерации.
Отметим трудности, возникшие при реализации приведенного алгоритма. Во-первых, при больших численностях популяций перебор всего массива частиц для вычисления времени до первой гибели в результате старения / | становится трудоемким и приводит к значительному увеличению времени вычислений.
d{x)=t,A
d{x)=t,A
Во-вторых, при большом количестве типов взаимодействий п трудоемким становится моделирование случайной величины /,•(/). Указанные трудности эффективно решаются путем использования древовидных структур, что приводит к логарифмической зависимости числа операций от количества частиц или взаимодействий соответственно,
В завершение укажем, что для моделирования равномерно распределенной на (0; 1) случайной величины применяется датчик псевдослучайных чисел М, В, Антипова [7] с периодом последовательности ~ 2 • 1036,
3. Примеры использования и тестирования имитационной модели
3.1. Исследование модели распространения эпидемии
В качестве примера рассмотрим дифференциальную модель Мэя-Андерсона распространения эпидемии с приобретением иммунитета [8]
{х = а(х + у + z) — bx — [3 ху + у z, у = (Зху- (a + b + v)y, z = vy — (b + 7 )z.
Здесь x = x(t) - число чувствительных к болезни частиц, у = y(t) - число заболевших, z = z(t) — число частиц, которые приобрели иммунитет после выздоровления, Все коэффициенты этой модели - неотрицательные действительные числа и интерпретируются следующим образом: а - интенсивность рождения новых частиц;
b - интенсивность смерти частиц в результате воздействия внешних факторов; /3 - интенсивность заражения в случае контакта чувствительной частицы с зараженной;
у - интенсивность потери иммунитета;
а - интенсивность гибели частицы из-за болезни;
v - интенсивность выздоровления и приобретения иммунитета.
Известно, что в приведенной модели существует устойчивое ненулевое положение равновесия (х0, уа, z0), которое возникает при определенных соотношениях на ее коэффициенты. Такое решение указывает на возможность ограничения численности популяций за счет распространения в ней некоторого заболевания, допускающего смертельный исход. Настоящий пример представляет собой одну из вероятностных модификаций модели Мэя-Андерсона, в которой опускается предположение о том, что распределение времени жизни частиц экспоненциально, Целью проводимых расчетов являлось нахождение решений имитационной модели, которые соответствовали бы указанному выше устойчивому решению дифференциальной модели.
Перепишем модель Мэя-Андерсона в терминах представленной имитационной модели. Рассмотрим три популяции X. Y и Z соответственно: чувствительных, заболевших и иммунных частиц. Имя первой популяции пусть будет "X",
имя второй - "Y", имя третьей - "Z", Распределения времени жизни частиц этих популяций имеют вид
L\[s] +оо) = ЬЦз; +оо) = e^bs,
L2[s; +оо) = L^[s] +оо) = e^(a+h+v^s,
L3[s; +oo) = L°[s; ±oo) = e-(h+^s, О 0.
В момент гибели в результате старения частицы популяции X ничего не происходит, значит (£х) = 0, (г]i) = 0,
Если в результате старения погибнет частица популяции ) . то она е вероятностью р = v/(a + b + v) порождает одну частицу популяции Z, то сеть выздоравливает и приобретает иммунитет. Полагаем тогда, что {£,2) = 0, (щ) = p*Zl.
В случае гибели в результате старения частицы популяции Z она с вероятностью q = 7/(6 + 7) порождает одну частицу популяции X. то есть теряет иммунитет. Поэтому принимаем, что (£3) = 0, {//.-,) = q * X1.
В этой модели существует два типа взаимодействий Д и /•_>. Интенсивность взаимодействий типа Д равна qi(x,y,z) = a * (х + у + z) (здесь х, у ж z -численности соответствующих популяций), а результат (ух) = 0, (/Д) = X1 -этот тип взаимодействий описывает процесс появления новых частиц, У взаимодействий типа Д интенсивность равна q2(x,y,z) = /3 * х * у, а результат (72) = X1, (/Д) = У1. Тип взаимодействий !■> описывает процесс заражения чувствительных частиц.
Рассмотрим дифференциальную модель Мэя-Андерсона с параметрами о = 5, b = 1, a = 8, /3 = 0, 01, 7 = 1, v = 1, В этом случае устойчивое положение равновесия существует и равно (x0,y0,z0) = (1000,2000,1000),
При запуске имитационной модели с этими параметрами начальные численности популяций х° = т(0), у° = у(0) и z° = z(0) выбирались: х° - из множества {500,1000,10000,100000}, у° - из множества {50,100, 2000,10000}, а : из
множества {0,100,1000, 2000}, Оценка математических ожиданий численностей x(t), y(t) и z(t) в моменты времени / = 10. / = 15 и / = 20. произведенная по 64 реализациям с доверительной вероятностью 0,95, приведена в таблице 1,
Таблица 1,
t x(t) y(t) z(t)
10 993, 27 ±45, 23 1996,19 ±92, 00 1002, 78 ±46,42
15 1000, 86 ±45, 64 1990, 67 ±91, 38 998, 34 ±45, 86
20 1002, 84 ±46,12 1983, 39 ±91, 62 979, 66 ±45,18
Полученные результаты указывают на наличие у модели устойчивого состояния, причем состояние (1000, 2000,1000) попадает в оценочные интервалы.
Все оценки, встречающиеся в статье далее, также будут производиться с доверительной вероятностью 0, 95,
Предположим далее, что для некоторой болезни накоплен статистический материал о длительности ее протекания и о длительности действия иммунитета к ней. Пусть эти данные описываются рядами {lyi, Д2, ■ ■ ■, 1Ут} и {hi, Ьл, ■ ■ ■, hs} соответственно. Тогда в качестве L-, и 1.Л можно взять дискретные распределения такие, что L2{lyk} = l/r, и L%{lzk} = 1/s, 1 ^ k ^ s. Так
например, можно положить L-,{/,, }• = 1, где ly = (a ± гД1, a 1.2{/~ }• = 1, где lz = 7-1, В качестве L\ и Щ можно взять равномерные на интервалах (0',1У) и (0; lz) распределения. Новые распределения Ь2 и L3 не учитывают гибели частиц в результате воздействия внешних факторов. Это можно исправить, добавив в модель еще два типа взаимодействий
h ■ qs(x, y,z) = b*y, (73) = Y\ (fa) = 0, h : qi(x,y,z) = b* z, (74) = Z1, {(34) = 0,
и изменив вероятности p = v/(a + b + v), q = 7/(6 + 7) на p = v/(a + v) и q = 1, Результаты оценок математических ожиданий численностей популяций для модели с новыми параметрами приведены в таблице 2,
Таблица 2,
t x(t) y(t) z(t)
10 15 20 945, 91 ±43, 37 3316, 97 ± 156, 77 1993, 83 ± 95, 91 953, 73 ±43,40 3253, 64 ± 149,11 1953, 84 ± 89, 63 955, 72 ±43, 37 3242, 81 ± 149, 53 1939, 98 ± 88, 52
По этим данным также можно сделать вывод о наличии стационарного состояния модели, и это состояние отличается от стационарного состояния модели Мэя-Андерсона, Таким образом, замена распределения времени жизни частиц может привести к значительному изменению предельных численностей популяций,
3.2. Тест 1
Имеется две популяции способных Ai и неспособных А2 к воспроизводству частиц, Распределения времени жизни частиц этих популяций имеют вид
L\[s] +00) = Ll[s-,+ 00) = e-k2S,
L2[s; +00) = L^ls] +00) = e-fe3S, s ^ 0,
Результаты гибели вследствие старения задаются соотношениями
(Ci) = 0, (m) =
(С2) = 0, Ы = 0-
За воспроизводство новых частиц отвечает один тип взаимодействий
h : qi(xi,x2) = ki* Xi, {yl}=0, {/3i)=X1.
Для данной модели известно, что при /7 = />-;! и х., = 0 математическое ожидание суммарной численности частиц в момент времени / имеет вид
E(t) = E(xi{t) + x2{t)) = X°l (2A;ie(fel-fe2)* - k2e^klt) .
lk\ — /С2
Оценка математического ожидания и его действительное значение, вычисленное по указанной формуле, приведены в таблицах 3 и 4 для различных параметров модели (усреднение производилось по результатам 1000 испытаний).
Таблица 3, .г, =10 Ад = 0.2 Ад = 0. 1
t E(t) E(t) (оценка)
10 20 35,79 35, 738 ±0,894 98,48 98,136 ±3,031
Таблица 4, .г, = 100 Ад = 0. 1 Ад = 0,15
t E(t) E(t) (оценка)
10 132,25 132,81 ± 1,511
20 106,55 106, 52 ± 1,883
40 48,639 47, 96 ± 1,632
3.3. Тест 2
Рассматриваются пять популяций . 1|. А2, Л3. А4 и Ч5, в каждой из которых объединяются частицы, находящиеся в одном «состоянии». Родившись, частица должна последовательно пройти все пять «состояний» и в конце своего жизненного цикла породить две новые частицы, если не погибнет ранее.
Распределения времени жизни частиц рассматриваемых популяций и результаты гибели вследствие старения имеют вид
Li[s; +00) = Ll[s; +00) = e^kl s, (СО = 0, (Vi) = Ah
L2[s;+oo) = Ц[в;+оо) = e-k2S, (6) = 0, (m) = A\,
L3{r3} = 1, LI - равномерное на (0; т3), (Сз) = 0, {щ} = А\,
L4[s; +00) = L\[s] +00) = e^ki s, (£4) = 0, \m) = Ah
L5{t5} = 1, L5-равномерное на (0;T5), (£5) = 0, (щ) = A(,
где s > 0, r3 > 0, т5 > 0, h > 0, к2 > 0, А:4 > 0.
Частицы популяции . 1;! могут также погибать и в результате воздействия внешних факторов:
h ■ qi(x) = h * х3, (71) = 0, </5i) = 0.
Если значения параметров Ад и т3 таковы, что еГкгТг ^ 1/2, то в этом случае все популяции вырождаются с вероятностью 1. Сопоставление данного утверждения с результатами модельных расчетов происходило следующим образом. Фиксировались параметры к\ = 1, к2 = 2, Ад = 0,1, Ад = 1, т5 = 6, а параметр т3 уменьшался, начиная с т3 = 10, и для каждого т3 оценивалось математическое ожидание Е(Т) времени Т до вырождения всех популяций. Результаты этого теста приведены в таблице 5. Также производилась оценка вероятности p(t) вырождения популяций к моменту времени /. Оценки приведены в таблицах 6 и 7. Для всех оценок производилось 1000 испытаний, каждое из которых стартовало с численностей х° = (10,0,0,0,0).
Таблица 6. т3 = 10 t p(t) (оценка)
~40 0, 092 ±0,0183
60 0, 251 ± 0,0274
80 0,448 ±0,0315
100 0, 637 ± 0,0304 120 0, 765 ± 0,0268 160 0, 897 ±0,0192 180 0, 937 ±0,0154 200 0, 958 ±0,0127
Таблица 7. т3 = 8
t p{t) (оценка)
100 0, 246 ± 0, 0272
200 0,585 ±0,0312
300 0, 779 ± 0, 0262
350 0,844 ±0,0229
400 0,897 ±0,0192
Таблица 5.
n E{T) (оценка)
10 110, 691 ±4, 570
9 141, 047 ±6, 080
8 205,486 ± 10,035
7,5 295, 755 ± 16,808
7,2 437,156 ± 28, 662
7 788, 025 ± 86, 303
Из этих таблиц видно, что с приближением величины еНкзП к пороговому значению 1/2, которому соответствует тз рз 6,931, время до вырождения популяций растет, а вероятность раннего вырождения уменьшается,
В заключение отметим, что предложенная имитационная модель является развитием стохастической модели, построенной е использованием точечных распределений [9,10],
Литература
1. Бартлетт М.С. Введение в теорию случайных процессов. М.: Иностранная литература, 1958.
2. Севастьянов Б.А. Ветвящиеся процессы. М.: Наука, 1971.
3. Севастьянов Б.А., Калинкин А.В. Ветвящиеся случайные процессы с взаимодействием частиц // Доклады АН СССР. 1982. Т.264, JY» 2. С.306-308.
4. Ермаков С.М., Михайлов Г.А. Курс статистического моделирования. М.: Наука, 1976.
5. Соболь II.M. Численные методы Монте-Карло. М.: Наука, 1973.
6. Тихонов В.И., Миронов М.А. Марковские процессы. М.: Советское радио, 1977.
7. Пригарин С.М. Введение в численное моделирование случайных процессов и полей. Часть II. Новосибирск: НГУ, 1999.
8. Anderson R.M., Mav R.M. Population biology of infectious disease. Part Iff Nature. 1979. V.280. P.361-367.
9. Перцев Н.В. Вероятностная модель динамики взаимодействующих частиц с ограниченным временем жизни // Математические структуры и моделирование. Омск. ОмГУ. 1998. Вып.1. С.60-71.
10. Харрис Е. Теория ветвящихся случайных процессов. М.: Мир, 1966.
