Добавил:
Upload Опубликованный материал нарушает ваши авторские права? Сообщите нам.
Вуз: Предмет: Файл:
Лабораторные работы по теории вер..pdf
Скачиваний:
40
Добавлен:
04.06.2015
Размер:
2.62 Mб
Скачать

Работа 6

СЛУЧАЙНЫЕ ПРОЦЕССЫ С ДИСКРЕТНЫМИ СОСТОЯНИЯМИ И НЕПРЕРЫВНЫМ ВРЕМЕНЕМ.

Решение задач с использованием матрицы интенсивностей. Численное моделирование динамики

распределения вероятностей состояний. Статистическое моделирование процесса с непрерывным временем

Время на выполнение и защиту 2 часа

Цель работы:

1)компьютерное решение задач о процессах с непрерывным временем (динамика распределения вероятностей системы с выходом на стационарный режим);

2)моделирование переходов в цепи Маркова с непрерывным временем по методу Монте-Карло;

3)изучение ряда функций Excel и Mathcad.

Матрица интенсивностей

В ситуации, когда изменения состояний системы могут происходить в любые моменты времени, говорят о случайном процессе с непрерывным временем. Очевидно, описание такого процесса в терминах вероятностей перехода за определённое количество шагов, как это делается для процесса с дискретным временем (работа 5), уже не годится. Действительно, в этом случае говорить о вероятности перехода системы из одного состояния в другое в любой определённый момент времени бессмысленно (она равна нулю, подобно вероятности любого значения непрерывной случайной величины). Можно говорить только о вероятности перехода за промежуток времени t (подразумевается, что за малый промежуток изменение не может произойти более 1 раза).

Рассмотрим систему с дискретными состояниями и непрерывным временем. Пусть в момент t система находилась в состоянии Si . Вероятность пере-

хода системы в состояние S j в течение промежутка времени t , отсчитывая с момента t , обозначим через Pij (t + ∆t) :

Pij (t + ∆t) = P(S(t + ∆t) = S j

 

S(t) = Si ) .

(6.1)

 

В таком случае обозначение Pij (t) будет соответствовать вероятности перехода в течение нулевого промежутка времени или, иначе говоря, вероятности того,

78

что система одновременно (в момент t ) находится и в состоянии Si , и в состоянии S j . Легко понять, что

Pij (t) = P(S(t) = S j

 

S(t) = Si )

0,

i j,

(6.2)

 

=

i = j.

 

 

 

 

 

 

 

 

 

1,

 

 

 

 

 

 

 

 

 

 

 

Плотностью вероятности перехода системы из состояния Si

в состоя-

ние S j (или интенсивностью перехода) называется величина

 

λij

(t) = lim

Pij (t + ∆t) Pij (t)

= Pij(t) .

(6.3)

 

 

 

t0

 

 

 

t

 

 

 

 

 

 

Используя выражение (6.2) для Pij (t), получаем:

 

 

 

 

 

 

P

(t + ∆t)

 

 

 

 

 

 

 

lim

ij

 

 

 

0,

 

i

j,

 

 

 

 

t

 

 

λij

t0

 

 

 

 

 

 

 

(6.4)

(t) =

Pii (t + ∆t) 1

 

 

 

 

 

lim

0, i = j.

 

 

 

 

 

 

 

 

 

t

 

 

 

 

 

 

 

t0

 

 

 

 

 

 

 

 

 

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

Из формулы 6.4

вытекает следующий важный вывод. Если t мало, то

Pij (t + ∆t) λij t при

i j . При совпадении начального и конечного состояний

перехода Pii (t + ∆t) 1+ λii t .

Если λij не зависит от момента времени t , с которого начался промежу-

ток времени t , то процесс называется однородным (аналогичное определение было в работе 5 и для процесса с дискретным временем). Далее мы будем рассматривать именно однородные процессы. В этом случае можно пользоваться так называемой матрицей интенсивностей

 

λ

λ

...

λ

 

 

11

12

...

1k

 

λ21

λ22

λ2k

Λ =

 

 

...

...

 

... ...

 

 

λk1

λk 2

...

 

 

 

λkk

(аналог матрицы переходных вероятностей). Рассмотрим её свойства.

79

1. Диагональные элементы этой матрицы неположительны, остальные элементы – неотрицательны: λij 0 при i = j ; λij 0 при i j (6.4).

2. Сумма элементов каждой строки равна нулю, т.е.

k

 

 

 

 

 

λij

= 0,

i =

1,k

.

(6.5)

j=1

 

 

 

 

 

Пример 17. По матрице интенсивностей

 

3

1

2

 

0

2

2

 

Λ =

 

 

1

0

 

 

 

 

1

построить размеченный граф состояний.

 

 

 

 

 

Рассматривая первую строку матрицы интенсивно-

 

 

 

S2

 

 

1

стей, приходим к выводу, что из состояния S1 возможны

 

S1

 

 

 

 

переходы в состояние S2 (с интенсивностью 1) и в состоя-

 

2

2

 

 

 

 

 

ние S3 (с интенсивностью 2); из состояния S2 возможен

 

 

 

1

 

S3

только переход в состояние S3

(с интенсивностью 2), а из

 

 

 

 

состояния S3 — переход в S1

(с интенсивностью 1). Эти

 

 

 

 

Рис. 6.1 рассуждения приводят нас к графу, изображённому на рис. 6.1.

Системы уравнений Колмогорова

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

пример, имеется три дискретных состояния и P(t) = (P1(t) P2 (t) P3 (t)) – вектор вероятностей того, что в момент времени t система будет находиться в со-

стояниях S1 , S2 , и S3 соответственно, а P(t) = (P1(t) P2(t) P3(t)) – вектор производных этих вероятностей. Тогда система дифференциальных уравнений Колмогорова должна быть записана в виде

(6.6)

P (t) = P(t)Λ.

Так, в условиях примера 17, учитывая правило перемножения матриц, имеем:

P1(t) = −3P1(t) + P3 (t),P2(t) = P1(t) 2P2 (t),

P3(t) = 2P1(t) + 2P2 (t) P3 (t).

80

Полученную систему уравнений Колмогорова следует решать с заданными начальными условиями. Например, если при t = 0 система находилась в состоянии S1 , то начальные условия будут таковы:

P1(0) =1, P2 (0) = 0, P3 (0) = 0.

В системе уравнений Колмогорова всегда одно уравнение являются линейной комбинацией остальных. Поэтому система должна быть дополнена условием нормировки Pi (t) =1 для любого t .

i

Предельный стационарный режим

Условие стационарности имеет вид

 

 

 

 

i =1,k ,

(6.7)

Pi (t) = 0,

или (в матричном виде)

 

 

 

 

 

 

st Λ = 0,

(6.8)

 

P

где Pst = (P1, P2 , ..., Pk ) .

В условиях примера 17 найдём стационарное распределение вероятностей, добавив в систему условие нормировки вероятности:

3P1 + P3 = 0,

P1 2P2 = 0,

2P1 + 2P2 P3 = 0,P1 + P2 + P3 =1.

3-е уравнение является комбинацией двух первых, исключаем его. Решая систему, получаем единственное решение

P1 = 92 , P2 = 19 P3 = 23 .

Для эргодического процесса стационарное распределение вероятностей есть предельное распределение вероятностей при неограниченном возрастании времени:

P = lim P (t), i =1, 2, ..., k .

i

t→∞ i

Условие протекания такого процесса в системе, как уже говорилось в работе 5, следующее: все существенные состояния системы должны быть сообщающимися. Действительно, в примере 17 все состояния являются существен-

81

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

Численное моделирование динамики процесса с непрерывным временем

Безусловно, система уравнений Колмогорова (6.6) может быть решена аналитически, однако это непростая задача, усложняющаяся с ростом числа состояний. Мы рассмотрим здесь очень простое численное решение, которое легко может быть реализовано как в Excel, так и в Mathcad.

Применим аппроксимацию производных по времени:

Pi (t + ∆t) Pi (t)

 

 

 

 

 

 

 

 

i =1,k ,

 

Pi (t)

 

,

 

t

 

 

 

 

 

 

 

 

 

откуда

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

(6.9)

 

 

 

 

 

 

 

Pi (t + ∆t) = Pi (t) + Pi (t)t , i =1,k ,

или в матричном виде P(t + ∆t) = P(t) + P(t)t .

Используя это равенство, при заданном начальном векторе P(0) можно

получить временной ход P(t) с шагом t . На каждом шаге для вычисления вектора производных используем равенство (6.6)

P(t) = P(t)Λ.

Статистическое моделирование случайного процесса с дискретными состояниями и непрерывным временем

Идея состоит в том, чтобы свести процесс с непрерывным временем к процессу с дискретным временем. Будем считать, что время изменяется дискретно с некоторым малым шагом t . Мы знаем, что если t мало, то

Pij (t + ∆t) λij t

(6.10)

при условии, что i j . Напротив, при совпадении начального и конечного состояний перехода,

Pii (t + ∆t) 1+ λii t .

(6.11)

Пример 18 (продолжение примера 17). Задана матрица интенсивностей переходов для процесса с непрерывным временем

3

1

2

 

 

 

 

 

 

Λ =

0

2

2

.

 

 

 

 

 

 

1

0

 

 

 

1

82

Записанные приближённые равенства (6.10) и (6.11) могут быть использованы для того, чтобы перейти от процесса с непрерывным временем, характеризуемого матрицей интенсивностей Λ, к процессу с дискретным временем, харак-

теризуемому матрицей перехода

~

Например,

при t = 0,1 можно получить

P .

матрицу

 

 

 

 

 

 

 

~

 

 

0,7

0,1

0,2

 

 

 

 

0

0,8

0,2

 

,

P

=

 

 

 

 

0,1

0

0,9

 

 

 

 

 

 

 

удовлетворяющую всем необходимым требованиям к матрице перехода. Решение подобной задачи методом Монте-Карло (только с меньшим числом состояний) рассматривалось в работе 5 (задание 5.2) .

Задания для лабораторной работы

В данной работе случайные процессы с дискретными состояниями и непрерывным временем моделируются двумя способами: численным решением системы уравнений Колмогорова и методом Монте-Карло.

Задание 6.1. Марковский процесс с непрерывным временем задан матрицей интенсивностей

2

0

2

 

 

6

8

2

 

Λ =

.

 

3

4

7

 

 

 

В начальный момент времени система находится в состоянии S1 . Требу-

ется исследовать динамику распределения вероятностей системы и получить стационарное распределение вероятностей.

Задание 6.2. Задана матрица интенсивностей переходов для процесса с непрерывным временем (примеры 17, 18):

3

1

2

 

 

 

 

 

 

Λ =

0

2

2

.

 

 

 

 

 

 

1

0

 

 

 

1

Найти стационарное распределение вероятностей по состояниям, разыграв случайный процесс методом Монте-Карло.

Инструкция по выполнению заданий в Excel

Приступим к выполнению задания 6.1. Заданную матрицу интенсивностей введём в диапазон B3:D5. Стационарное распределение вероятностей по

83

состояниям должно удовлетворять матричному уравнению Pst Λ = 0, дополненному условием нормировки вероятности. Для этого нужно из системы одно-

родных уравнений PstΛ =0, где Pst = (P1, P2 , ..., Pk ) , исключить одно (любое) уравнение, заменив его условием

k

Pi =1.

i=1

Врезультате система примет вид Pst A = B , где матрица A есть изменённая матрица Λ (скажем для определённости, что последний столбец заменён еди-

ницами),

а

B = (0, 0, ..., 0, 1)

матрица

размером

1× k . Таким

образом,

Pst

= BA1 . Реализовав этот алгоритм с помощью функций МОБР и МУМНОЖ

подобно тому, как это описано в инструкции по выполнениюзадания 5.1, полу-

чим стационарное распределение вероятностей по состояниям, которое (в про-

стых дробях) выглядит как P1 = 2 / 3, P2 =1/ 9,

P3 = 2 / 9 .

 

 

1

A

B

C

D

E

F

G

H

I

J

 

 

 

 

 

 

 

 

 

 

2

Матрица интенсивностей

Матрица А

 

Матрица А обр

 

3

 

-2

0

2

-2

0

1

-0,16667

0,055556

0,111111

4

 

6

-8

2

6

-8

1

-0,04167

-0,06944

0,111111

5

 

3

4

-7

3

4

1

0,666667

0,111111

0,222222

6

 

 

 

 

Матрица В

Стац распред. вероятностей

7

 

 

 

 

0

0

1

0,666667

0,111111

0,222222

Перейдём к численному моделированию динамики распределения вероятностей системы. В ячейку C7 введём значение t = 0,05. (Отметим здесь, что шаг по времени должен быть достаточно мал). В ячейку A10 введём начальное значение времени (0), а в ячейку A11 — формулу =A10+$C$7. В диапазон

B10:D10 введём начальное распределение вероятностей состояний P(0) . Выделим диапазон E10:G10 и введём в него функцию МУМНОЖ, умножающую

матрицу P(0) на матрицу интенсивностей. (Напомним, что это функция массива, и при её применении потребуется перейти в режим правки F2, как это описано в инструкции по выполнению задания 5.1. Учтите также, что ссылки на массив матрицы интенсивностей должны быть безусловными для последующего автозаполнения таблицы).

В ячейки B11:D11 необходимо ввести формулы (6.9), которые будут вычислять вероятности состояний на новом временном шаге через значения вероятностей и их производных с предыдущего шага. Для ячейки B11 это будет выглядеть как =B10+E10*$C$7, а к С11 и D11 можно применить автозаполнение. Теперь выделим диапазон E10:G10 и выполним автозаполнение в область E11:G11. После этого остаётся лишь выделить строку A11:G11 и выполнить автозаполнение произвольного количества строк в направлении вниз.

84

 

A

B

C

D

E

F

G

H

I

J

1

 

 

 

 

 

 

 

 

 

 

2

 

Матрица интенсивностей

Матрица А

Матрица А обр

3

 

-2

0

2

-2

0

1

-0,1667

0,05556

0,11111

4

 

6

-8

2

6

-8

1

-0,0417

-0,0694

0,11111

5

 

3

4

-7

3

4

1

0,66667

0,11111

0,22222

6

 

 

 

 

Матрица В

Стац распред. вероятностей

7

 

дельта t

0,05

 

0

0

1

0,66667

0,11111

0,22222

8

 

 

 

 

 

 

 

 

 

 

9

t

Вероятности

 

Производные

 

 

 

10

0

1

0

0

-2

0

2

 

 

 

11

0,05

0,9

0

0,1

-1,5

0,4

1,1

 

 

 

 

 

 

Изучая полученные результаты, легко убедиться, что:

1)с ростом времени распределение вероятностей состояний приближается к стационарному распределению, полученному нами ранее в ячейках H7:J7;

2)производные от вероятностей при этом приближаются к нулю (этим критерием и определяется количество строк для автозаполнения вниз).

Приступим к выполнению задания 6.2. В примере 18 было показано, что при введении искусственной дискретизации времени с шагом t = 0,1 из исход-

ной матрицы интенсивностей можно получить матрицу переходных вероятностей

~

 

0,7

0,1

0,2

 

 

0

0,8

0,2

 

P

=

.

 

 

0,1

0

0,9

 

 

 

 

Таким образом, перейдя к процессу с дискретным временем, мы получили задачу, отличающуюся от задания 5.2 лишь числом возможных состояний процесса. Приступим к её решению методом Монте-Карло. В диапазон ячеек B2:D4 введём матрицу интенсивностей, в ячейку D5 – выбранный шаг по времени, в диапазон B6:D8 формулы, вычисляющие переходные вероятности в соответствии с правилами (6.10) и (6.11). Диапазон ячеек A9:C11 заполняется подобно тому, как это было сделано в задании 5.2. Однако формула для генера-

 

 

 

A

B

C

D

1

 

Матрица интенсивностей

 

 

2

 

 

-3

1

2

 

 

 

 

3

 

 

0

-2

2

4

 

 

1

0

-1

 

5

 

Матрица

перехода за время t =

0,1

6

 

 

0,7

0,1

0,2

 

7

 

 

0

0,8

0,2

 

8

 

 

0,1

0

0,9

 

 

 

 

9

 

Время

Сл. число

Состояние

 

 

10

 

0

 

1

 

 

11

 

=A10+$D$5

=СЛЧИС85()

Формула

 

 

 

 

 

 

 

 

(6.12)

 

 

12

 

 

 

 

 

ции состояния процесса на текущем шаге (в ячейке C11) будет, разумеется, несколько более сложной, чем в задании 5.2 (из-за увеличения размерности матрицы перехода):

=ЕСЛИ(C10=1;ЕСЛИ(B11<$B$6;1;ЕСЛИ(B11<$B$6+$C$6;2;3)); ЕСЛИ(C10=2;ЕСЛИ(B11<$B$7;1;ЕСЛИ(B11<$B$7+$C$7;2;3)); (6.12)

ЕСЛИ(B11<$B$8;1;ЕСЛИ(B11<$B$8+$C$8;2;3)))).

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

Дальнейшие действия, по сути, ничем не отличаются от того, что мы делали при решении задания 5.2 (автозаполнение, подсчёт времени, проведённого системой в разных состояниях, оценивание предельных вероятностей). Сравните оценки вероятностей состояний с ранее полученным нами аналитическим решением примера 17

P1 = 92 , P2 = 19 P3 = 23 .

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

Инструкция по выполнению задания в Mathcad

Приступим к выполнению задания 6.1. Присвоим переменной ORIGIN значение, равное 1. Стационарное распределение вероятностей по состояниям

должно удовлетворять матричному уравнению Pst Λ =0 , дополненному условием нормировки вероятности. Для этого нужно из системы однородных уравне-

ний Pst Λ =0 , где Pst = (P1, P2 , ..., Pk ) , исключить одно (любое) уравнение, заменив его условием

k

Pi =1.

i=1

Врезультате система примет вид Pst A = B , где матрица A есть изменённая матрица Λ (скажем для определённости, что последний столбец заменён еди-

ницами), а B = (0, 0, ..., 0, 1) — матрица размером 1× k . Таким образом,

Pst = BA1 .

86

Для реализации этого алгоритма введём заданную в условии матрицу ин-

2

0

2

 

 

 

6

8

2

 

и указанную матрицу B = (0, 0,1) , используя

тенсивностей Λ =

 

 

3

4

7

 

 

 

 

 

панель инструментов Матрицы. Далее матрице A присвоим значения матрицы Λ и, используя дискретные переменные, присвоим элементам последнего столбца значения, равные 1 ( i :=1..3 Ai,3 :=1). Напомним, что знак присваива-

ния вводится клавишей « : », а выражение 1..3 можно ввести последовательно нажимая «1» « ; » «3» или воспользовавшись кнопкой панели Матрицы.

Затем введем выражение Pst := B A1 для вычисления стационарного распределения вероятностей по состояниям. Знак умножения вводится нажатием кла-

виши « * », а знак обратной матрицы щелчком по кнопке панели Матрицы. В результате получим Pst := (0.667, 0.111, 0,222), что является приближённой

записью простых дробей P1 = 2 / 3, P2 =1/ 9, P3 = 2 / 9 .

Перейдём к численному моделированию динамики распределения вероятностей системы.

У нас имеется три дискретных состояния и P(t) = (P1(t) P2 (t) P3 (t)) – вектор вероятностей того, что в момент времени t система будет находиться в

состояниях S1 , S2 , и S3 соответственно, а P(t) = (P1(t) P2(t) P3(t)) – вектор производных этих вероятностей. Тогда система дифференциальных уравнений Колмогорова (6.6) для условий задания 6.1 должна быть записана в виде:

P1(t) = −2P1(t) +6P2 (t) +3P3 (t),P2(t) = −8P2 (t) + 4P3 (t),

P3(t) = 2P1(t) + 2P2 (t) 7P3 (t).

Полученную систему уравнений следует решать с заданными начальными условиями. Поскольку по условию при t = 0 система находилась в состоянии S1 , то начальные условия будут таковы:

P1(0) =1, P2 (0) = 0, P3 (0) = 0.

В Mathcad решить приближенно заданную систему линейных дифференциальных уравнений можно с помощью функции rkfixed(p0, t1, t2, intvls, D) в категории Differential Equation Solving (Решение дифференциальных уравнений) с использованием метода Рунге-Кутта с постоянным шагом.

Здесь, t1, t2 – начальная и конечная точки отрезка интегрирования, параметр intvls задает число строк матрицы результатов (число шагов интегрирования), D – имя вектора-функции D(t,p), содержащей правые части системы дифференциальных уравнений. В нашем случае t1= 0. Отметим, что шаг интегрирования по времени должен быть достаточно мал. Возьмем его равным

87

0,05. Зададим число шагов равным 20. Тогда t2 будет равно 1. Если полученная точность будет недостаточной, можно увеличить число шагов и получить новый результат.

Введем вектор начальных условий

1

p0 := 0

0

Зададим вектор-функцию D правых частей системы уравнений следующим образом:

2 p1 + 6 p2 + 3 p3

D(t ,p) := −8 p2 + 4 p3

2 p1 + 2 p2 7 p3

Для этого нужно:

ввести имя D(t,p) и знак присваивания (двоеточие);

щелкнуть по значку в панели Матрицы. В появившейся диалоговой панели ввести число строк (3) и столбцов (1).

После нажатия кнопки OK открывается поле для ввода элементов векто-

ра. Заполните метки-заполнители соответствующими выражениями. Напомним, что для ввода индексных выражений нужно обязательно нажать клавишу « [ » (левую квадратную скобку). При этом курсор перемещается вниз, и индексные выражения оказываются смещенными относительно имени массива. После завершение ввода индексных выражений нужно нажать клавишу «пробел» или «

».

Обращение к функции rkfixed будет выглядеть следующим образом:

Pst1 := rkfixed(p0 ,0,1,20,D)

Результат получим в виде матрицы, имеющей 4 столбца и 21 строку. В первом столбце указано время, во втором, третьем и четвертом приближенные значения P1(t) , P2 (t) и P3(t) соответственно. Если матрица результата видна не

полностью, щелкните в ячейке, содержащей знак многоточия (в правом нижнем углу). Появится знак прокрутки, с помощью которого можно посмотреть результат до конца.

Изучая полученные результаты, легко убедиться, что с ростом времени распределение вероятностей состояний приближается к стационарному распределению Pst , полученному нами ранее.

Приступим к выполнению задания 6.2. В примере 18 было показано, что при введении искусственной дискретизации времени с шагом t = 0,1 из исход-

88

ной матрицы интенсивностей можно получить матрицу переходных вероятностей

~

0,7

0,1

0,2

 

 

0

0,8

0,2

 

P

=

.

 

 

0,1

0

0,9

 

 

 

 

Таким образом, перейдя к процессу с дискретным временем, мы получили задачу, отличающуюся от задания 5.2 лишь числом возможных состояний процесса. Приступим к её решению методом Монте-Карло.

Для этого, как и в предыдущей работе, составим программу-функцию f(n, sn, р), где n – продолжительность реализации случайного процесса (10000 моментов времени), sn – начальное состояние, р – исходная матрица перехода цепи Маркова.

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

Определим эту формулу таким образом:

– система на i-ом шаге процесса будет находиться в состоянии S1, если выполняются условия:

( sn =1 и a < P11 ) или ( sn = 2 и a < P21 ) или ( sn = 3 и a < P31 );

– система на i-ом шаге процесса будет находиться в состоянии S3, если выполняются условия:

( sn =1 и a > P11 + P12 ) или ( sn = 2 и a > P21 + P22 ) или ( sn = 3 и

a > P31 + P32 ).

В противном случае система будет находиться в состоянии S2. Организуем эту процедуру с помощью оператора if . Знаки логических

отношений: >, <, =, (логическое И), (логическое ИЛИ) вводятся с помощью панели инструментов Булева алгебра.

Определим вектор Т=(T1 ,T2 ,T3 ), элементами которого будет суммарное

время, проведенное системой в состояниях S1, S2 и S2 соответственно. Стационарное предельное распределение вероятностей состояний может быть оценено как отношения T1 , T2 и T3 к общей продолжительности реализации случайного

процесса Т=10000. После завершения описания программы-функции присвоим переменной РР значение этой функции РР:=f(10000,1,р). При этом исходная матрица перехода цепи Маркова р должна быть определена до обращения к программе функции. Напомним, что переменные, обозначенные большими и малыми буквами, воспринимаются Mathcad как разные.

Реализация данной задачи будет выглядеть примерно так:

89

f(n,sn ,p) :=

 

for

i 1..3

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Ti

0

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

for

i 1..n

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

a rnd(1)

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

si 1 if sn

 

1 a < p1,1 sn

 

2 a < p2,1 sn

 

3 a < p3,1

 

 

 

 

 

 

 

 

 

 

 

 

 

si 3 if sn

 

1 a > p1,1 + p1,2 sn

 

2 a > p2,1 + p2,2 sn

 

3 a > p3,1 + p3,2

 

 

 

 

 

 

 

 

 

 

 

 

 

si 2 otherwise

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

sn si

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

T1 T1 + 1 if si

 

1

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

T2 T2 + 1 if si

 

2

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

T3 T3 + 1 otherwise

 

 

 

 

 

 

 

 

 

 

 

 

for

i 1..3

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

T

Ti

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

n

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

i

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

T

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

0.7

0.1

0.3

 

 

 

 

 

 

 

 

 

0.215

 

PP :=

f(10000,1,p)

 

 

 

p :=

0

 

0.8

0.2

 

 

PP

= 0.108

 

 

 

 

 

 

 

 

 

 

 

 

0

0.9

 

 

 

 

 

 

 

 

 

0.677

0.1

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Нажмите несколько раз комбинацию клавиш <Ctrl> + <F9>. При каждом нажатии будет генерироваться новая реализация продолжительностью T =10000, и, следовательно, будут появляться новые значения T1 , T2 и T3 .

Сравните оценки вероятностей состояний с ранее полученным нами аналитическим решением примера 17

P1 = 92 , P2 = 19 P3 = 23 .

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

Контрольные вопросы

1.Какие случайные процессы называют процессами с непрерывным временем?

2.Почему для описания процессов с непрерывным временем применяется не матрица перехода, а матрица интенсивностей?

90