»Sy3=limit(R,omega,-inf,'right') Sy3 =

-5/2*pi

»Dy=Sy2-Sy3

Dy = 5*pi

Задание. Задание определяется и выполняется так, как указано на с. 105 и в табл. 3.

4. МОДЕЛИРОВАНИЕ СЛУЧАЙНЫХ ПРОЦЕССОВ

4.1. Общая характеристика методов моделирования случайных процессов

К настоящему времени разработан ряд методов моделирования случайных процессов ξ(t) с различными заданными характеристи-

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

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

Моделирование случайных процессов сводится к моделированию независимых случайных величин. В свою очередь моделирование независимых случайных величин с любыми законами распределения осуществляется с помощью датчиков случайных чисел.

109

В настоящее время в практике цифрового моделирования случайных процессов довольно часто используется около десятка алгоритмов: метод формирующего фильтра, метод скользящего суммирования, методы канонических и неканонических разложений, методы авторегрессии-скользящего среднего, метод условных распределений, метод Неймана и т.д.

4.2. Метод условных распределений

Рассмотрим случайный процесс ξ(t) с заданной n -мерной плотностью вероятности f (x1,t1, x2 ,t2 ,..., xn ,tn ). Совместная

n -мерная плотность вероятности удовлетворяет следующим условиям:

1)положительности f (x1,t1, x2 ,t2 ,..., xn ,tn )≥ 0 ;

2)нормировки

∞	∞	∞
∫	∫	... ∫ f (x1,t1, x2 ,t2 ,..., xn ,tn )dx1dx2 ...dxn =1;	(4.1)
−∞ −∞		−∞

nраз

3)согласованности

f (x1,t1, x2 ,t2 ,..., xk ,tk )=

∞	∞	∞
∫	∫	... ∫ f (x1,t1x2 ,t2 ,..., xk ,tk , xk+1,tk+1,..., xn ,tn )×
−∞ −∞		−∞
	n−k	раз

×dxk+1dxk+2 ...dxn ,k < n

(см. также формулу (1.11)).

Цифровое моделирование случайного процесса ξ(t), заданно-

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

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

110

В общем виде n → ∞, на практике чем больше n , тем детальнее статистическое описание случайного процесса ξ(t). При за-

данном n случайный процесс ξ(t) представляется случайным вектором ξ1,ξ2 ,...,ξn , и моделирование случайного процесса можно

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

Метод основан на определении условных одномерных распре-

делений для каждого момента ti	случайного процесса
ξ(ti ),i =1,2,...,n . Каждое значение ξ(ti )	моделируется с помощью

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

ного процесса ξ(t1 )

на основе формулы (1.11) определяется безус-

ловная плотность вероятности:

f1 (x1,t1 )

∞

∞

∞

f

(x1,t1, x2 ,t2 ,..., xn ,tn )dx2dx3...dxn . (4.2)

= ∫

∫

... ∫

−∞ −∞

−∞

n−1 раз

Для второго и последующих отсчётов текущей реализации оп-

ределится

условная

одномерная

плотность

вероятности

f

x

2

,t

2

=

f12 (x1,t1, x

2 ,t2 )

,

где

по

формуле

(1.11)

x1

f (x ,t

)

2

,t1

1

1

∞

∞

∞

f (x1,t1, x2 ,t2 ,..., xn ,tn )dx3dx4 ...dxn ,

f12 (x1,t1, x2 ,t2 )= ∫

∫

... ∫

−∞ −∞

−∞

n−2 раз

f123 (x1,t1, x2 ,t2

,t3 )

аналогично

f

x

,t

, x3

и

3

3

=

f

(x ,t

, x

,t

)

3

x1,t1, x2 ,t2

12

2

2

1

1

f123 (x1,t1, x2 ,t2 , x3 ,t3 )

∞

∞

∞

f (x1,t1, x2 ,t2 ,..., xn ,tn )dx4dx5...dxn.

= ∫

∫

... ∫

−∞ −∞

−∞

n−3 раз

Очевидно, что для каждого очередного отсчёта условная плотность вероятности, определяемая таким же образом, будет рав-

на:

f

x

k

,t

k

=

f12...k (x1,t1, x2 ,t2 ,..., xk ,tk )

.

f

(x

,t , x

,t

,..., x

,t

)

k

x1,t1,..., xk−1,tk−1

12...k−1

2

2

k−1

k−1

1

1

111

Во всех этих формулах значения величин, стоящих в знаменателях, определяются на предыдущих шагах.

Приведённые рекуррентные соотношения для условных плотностей используются для построения n датчиков случайных чисел с данными (условными) плотностями вероятностей. Каждый i -й датчик моделирует i -й отсчёт конкретной реализации ξ(ti ) слу-

чайного процесса ξ(t).

Многомерная плотность вероятности f (x1,t1, x2 ,t2 ,..., xn ,tn )

выражается через условные плотности вероятностей следующим образом:

f (x

,t , x

,t

,..., x

,t

)=

f (x

,t ) f

x

2

,t

2

... ×

x1

1

1

2

2

n

n

1

1

1

,t1

(4.3)

x

,t

×

x

,t

, x

,t

,..., x

,t

f

n

n

2

2

n−1

.

1 1

n−1

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

вероятности f (x1,t1, x2 ,t2 ,..., xn ,tn ). Следует иметь в виду, что на практике приходится вычислять интегралы, выражающие условные плотности вероятности f12...k , причём часто эти интегралы не

выражаются в конечном виде и их вычисляют численно. Это один из недостатков метода условных распределений. Второй недостаток – метод требует реализации n датчиков случайных чисел с

несовпадающими распределениями f1, f2 ,..., fn . В общем случае

вид плотностей вероятностей f (xk ,tk x1,t1, x2 ,t2 ,..., xk−1,tk−1 ) может меняться при смене текущей реализации, что ещё более ус-

ложняет практическое применение метода (третий недостаток). Пример 16 [22]. Рассмотрим моделирование случайного про-

цесса с двумерной функцией плотности

f (x	,t	, x	2	t	2	)=	1		(x1t1 )2 + 4(x2t2 )2 + (1 L)2	, L > 0 .
									5
1	1
							4πL		[(x1t1 )2 + (x2t2 )2 + (1 L)2 ]2

Определим безусловную плотность вероятности для первого отсчёта t1:

112

f (x1,t1 )= ∞∫ f (x1,t1, x2 ,t2 )dx2 =

−∞

1

∞

(x t

)2

+ 4(x

2

t

2

)2

+ (1 L)2

=

∫

1 1

dx2 .

5

4πL −∞ [(x1t1 )2 + (x2t2 )2 + (1 L)2 ]2

Неопределённый интеграл этого типа принадлежит к классу

интегралов

вида

2

t

2

+ n

2

Он рационализируется и

∫R t, m

dt .

приводится

к

интегралу

от

тригонометрических функций вида

∫R(sin z,cos z)dz соответствующей подстановкой:

f (x ,t )=

(x t

)+

1

2

= a2

1

∞ a2 + 4(x

2

t

2

)2

dx

=

∫

2

=

5

1 1

1 1

L

4πL −∞

[a2 + (x2t2 )2 ]2

	x2	=		a		tgz,
				t2
	=					dz
	dx2 =		a
						cos2 z
			t2

		π	a2 + 4(atgz)2
	a	2				dz
=		∫						=
			[a2 + (atgz)2 ]25			cos2
	4πLt2 −π2						z

π

π

=

a

2∫

1+ 4tg2 z

dz

=

a

2∫cos z(cos2 z +

2

2

4πLt

a

π

cos

z

4πLt

a

π

2

[1+ tg2 z]2

2

−2

−2

+ 4sin2 z)dz =

a

=

1

=

πLt2a2

πLt2

[(x1t1 )2 + (1 L)2 ]

1

1

1

1

1 Lt1

.

=

=

πt2

(Lx t )2 +1

t1t2

π

x2

+ (1

Lt

)2

1 1

1

1

Полученное

распределение является распределением Коши,

для которого

f (x)=

λ

,

где µ – параметр положе-

π[λ2 + (x −µ)2 ]

ния,

а

λ > 0 –

параметр

рассеивания.

В

нашем случае

µ = 0,

λ =

1

, присутствует ещё множитель

1

.

Lt

t t

1

1 2

113

Воспользуемся методом обращения (обратной функции) для моделирования случайной компоненты ξ(t1 ). Здесь каждая компо-

нента моделируется как скалярная величина [23]. Найдём функцию распределения F(x1,t1 ):

			x														x		1	dx
														1					Lt1		1
	1																1
F(x1,t1 )= ∫					f (x1,t1 )dx1					=							∫									=
																								2
			−∞										πt1t2 −∞					2			1
																		x1		Lt
																			+
																					1
	1 x1				d (Lt x )									1			arctg(Lt x				)		x1
=					1 1					=															=
	πt1t2 −∞∫ 1 + (Lt x )2											πt1t2
																			1 1
					1 1																		−∞
	=		1				1	arctg(Lt x						)	+		1 .
											1 1
			t1t2 π														2

Используя выражение для F(x1, t1 ), получим моделирующий алгоритм для ξ(t1 ) методом обращения. Пусть γ1 R(0,1), тогда, по методу обратной функции,

		γ	1	=	1		1		arctg(Lt x				)+ 1 , γ						t t	2		− 1		=
					t1t2					1	1			2				1 1				2
						π
	1	arctg(Lt				x		),ξ(t		)= x				1		tg						t t			1		.
=												=							π γ				2	−
	π													Lt											2
				1 1				1		1											1 1
														1

Условная плотность вероятности x2 , т.е. второго отсчёта случайного процесса ξ(t2 ), будет равна:

						x			,t											f (x1,t1, x2								,t2 )
					f			2		2						=														=
											x1												f	(x ,t				)
														,t1
																									1	1
																				1			2											1		2
			(x t			)2	+ 4(x						t		)2	+								t t				π x		2	+
																																		Lt
				1 1								2		2						L					1		2		1
	1																																	1
=																																					=
	4πL																											5
												2									2			1		2 2					1
							(x t				)			+ (x				t			)		+
																									L						Lt1
									1 1								2		2

114

														2																2					1				2														2										1					2
							(x t						)				+				4(x						t			)				+									(x t								)						+
																																					L																										L
		t2						1 1																	2			2																			1 1
=																																																																								=
	4																																																				5
																												2												2			1				2 2
																		(x t								)					+ (x					t			)	+
																																													L
																						1 1													2			2
										2										1				2								2												2		[a		2																			2		]
						(x t			)					+													= a							,					t	2			a								+ 4(x											t		)
	=					1 1														L																=																								2				2			5			,
							(Lx1t1 )2 = (La)2 −1																																4				[a2 + (x2t2 )2 ]2
где a2 = (x t						)2						1					2		,a =										[ξ(t					)t		]2							1 2
								+																															+						.
	1 1											L																					1		1								L
Для моделирования ξ(t2 ) найдём F																																											x			2	,t		2
																																																				x1									.
																																																									,t1
								x				2		,t			2																x2				x			2		,t	2
							F																						= ∫						f										x1						dx2 =
																				x1,t1												−∞															,t1
																																														x2 =											a						,
	a2t						x2		4(x									t			)2			+ a2																													t2tgz
																																																		a
																2			2																																								dz
=					2 ∫																												dx2 = dx2 =																																			, =
																														5																																	2
	4						−∞		[a2 + (x2t2 )2 ]2																																								t2								cos										z
																																											z = arctg															x2t2									.
																																																											a
			1																		x					2	t		2																		x					2			t	2											x2
	=							sin arctg																									+ sin3								arctg																													=
			4																		a																										a
																																																																		−∞

= 14 (sin z + sin3 z + 2).

Воспользуемся вновь методом обращения. Пусть γ2 R(0,1),

тогда		γ2	−	1		. Данное кубическое уравнение отно-
	sin3 z +sin z = 4			2

сительно z может быть решено по формулам Кардано. Напомним

115

их: решение кубического уравнения

у3 +3ру + 2q = 0

зависит от

значения дискриминанта D = q2 + p3 . Если D > 0 ,

то имеется один

действительный

и два

комплексных

корня. В

нашем

случае

1

1

= (2γ2 −1)

2

1

sin z = y, p =

, q

= −2 γ2

−

, D

+

> 0 .

Находим един-

3

2

27

1

ственный действительный

корень

y = u + v,

u = (2γ2 −1 +

)3 ,

D

1

x2t2

v = (2γ2 −1−

D )3 ,

т.е.

z = arcsin y ,

arctg

= arcsin y ,

a

y

ay

x

2

t

2

= atg(arcsin y)= atg arcsin

=

.

1− y

2

1

− y

2

Сведём воедино все формулы, решающие поставленную

задачу:

R(0,1),

ξ(t

)=

1

tg

t t

1

,

γ

1

π

γ

2

−

Lt

2

1

1 1

1

2

1

2

γ2 R(0,1),

a =

[ξ(t1 )t1 ]

+

,

L

1

1

(2γ2 −1)2 + p3 ,

p

=

, q = −2 γ

2 −

, D =

3

2

1

1

u = (2γ2 −1 +

)3 , v = (2γ2 −1 −

)3 , y = u + v ,

D

D

ξ(t2 )=

ay

.

1 − y2

Смоделируем теперь средствами пакета Mathcad 10 отсчётов в

моменты t1 = 2

и t2 = 5 разобранного примера.

Это можно сде-

лать, например, следующим образом:

116

117