4) процедуру повторяют необходимое число раз (п. 1-3), получают последовательность значений m1,m2 ,m3,... случайной величины

X B(n, p). Описанный метод называется методом браковки. В вычисли-

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

Моделирование случайной величины, распределенной по закону Пуассона. Распределение Пуассона P(λ) - предельная форма биномиаль-

ного распределения B(n, p) при n → ∞,	p → 0, np → λ , т.е.
lim Cnm pm (1− p)n−m =		λ me−λ	.		(3.4.2)
n→∞,		m!
np→λ
Алгоритм моделирования случайной величины, распределенной по
закону Пуассона, учитывая (3.4.2), может быть, например, таков:
1) задать λ и выбрать n такое, чтобы вероятность p = λ				n	была дос-
таточно малой (p < 0.01);

2)получить последовательность значений r1, r2 ,..., rn случайной величины R R[0,1];

3)для каждого числа ri , i =1,2,...,n проверить выполнение неравен-

ства ri < p. Если оно выполняется, то Xi =1, в противном случае считают

Xi = 0 ;

	n
4) вычислить	Y = ∑Xi .	Это и есть значение случайной величины
	i=1

YP(λ);

5)п. 1-4 повторить требуемое число раз.

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

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

Метод обратной функции основан на теореме Смирнова .

Теорема 3.1. (Смирнова). Если X удовлетворяет уравнению

X

∫dFX (t) = α, т.е.

−∞

Николай Васильевич Смирнов (1900-1966) - советский математик.

68

	X = FX−1(α),	(3.5.1)
где α - величина, распределенная равномерно на [0,1], то X распре-
делено по закону FX (t).
Действительно,	введем случайную величину α = FX (t), обратим вни-
мание на то, что так как 0 ≤ FX ≤1 , то и 0 ≤ α≤1. Найдем функцию рас-
пределения случайной величины α:
	0,	z ≤ 0,
Fα(z) = P(α< z) = P{FX (t)< z}= P{t < FX−1(z)}= FX {FX−1(z)}= α, 0 ≤ α≤1,
	1,	z ≥1.
Таким образом,	α- случайная величина, имеющая равномерное рас-
пределение на отрезке [0,1].

Порядок действия при моделировании конкретного распределения методом обратной функции следующий:

1)разыгрывается реализация α равномерной случайной величины

αR[0,1];

2)решается уравнение F(x)= α. Его решение X = F −1(α) дает слу-

чайную величину X с заданным законом распределения F(x). Здесь F −1 функция обратная к F(x). Во многих случаях удается найти явное выра-

жение для F −1(α). Тогда моделирование происходит наиболее просто. Недостатком описанного метода являются аналитические трудности

при вычислении F −1 . В «чистом виде» метод обратной функции используется редко на практике, так как для многих распределений, например, нормального, даже F(x) не выражается через элементарные функции, а

табулирование F −1(α) существенно усложняет моделирование. На практике метод обратной функции дополняют аппроксимацией F(x) или соче-

тают с другими методами.

Метод суммирования основан на центральной предельной теореме (ЦПТ). Например, в простейшем случае ЦПТ может быть сформулирована так.

Теорема 3.2. Если случайные величины X1, X 2 ,..., X n независимы,

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

X1 + X 2 +... + X n неограниченно приближается к нормальному.

Практически оказывается, что для получения хорошего приближения к нормальному распределению достаточно сравнительно небольшого чис-

69

ла слагаемых.		Пусть	r1, r2 ,..., rn - независимые		случайные	величины
ri R[0,1].	Обозначим через Y			сумму этих величин: Y = r1 + r2 +... + rn ,
тогда Mri	= 0.5	и Dri	=1 12, i	=1,2,..., n. Отсюда	MY = 0.5n,	DY = n 12.

При достаточно большом n по ЦПТ можно считать, что Y имеет нормальный закон распределения с математическим ожиданием MY = 0.5n и дисперсией DY = n12 , т.е.

Y N (0.5n,			n 12 ).				(3.5.2)
Перейдем к стандартной нормально распределенной случайной вели-
чинеU = Y − MY = (Y − 0.5n) 12		= (Y − 0.5n) 6			, U N (0,1).		Например,
DY	n			3n
при n =12
	n		6		12
				= MX			(3.5.3)
X = MX + σXU = MX + σX	∑ri − 0.5n		3n		+ σX ∑ri − 6	.
i=1					i=1
Тогда X N(MX ,σX ),	а ri R[0,1].		Таким образом, имея двенадцать
значений случайной величины		R , получаем значение нормальной слу-

чайной величины X , имея следующие двенадцать значений R , получаем следующее значение X и так далее.

Для других распределений основной принцип метода суммирования и его опора на ЦПТ остаются без изменения. Моделирующие формулы типа (3.5.3) для каждого конкретного распределения разные.

Рассмотрим еще один пример моделирования гамма-распределения.

		f (x) =	1				−	x		1
	Его плотность				xαe β , x > 0. Положим λ =						, тогда
			Γ(α +1)β α+1							β
			f (x) =	λ α+1			xαe−λx , x > 0.				(3.5.4)
				Γ(α+1)

Предлагаемый способ моделирования гамма-распределения основывается на следующей теореме.

Теорема 3.3. Если y1, y2 ,...yn - независимые стандартные экспо-

n

ненциально распределенные случайные величины, то x = ∑yn имеет i=1

гамма-распределение с параметром α = n −1.

70

Если y - экспоненциально распределенная случайная величина, то

f (y) = λe−λ y , y ≥ 0, λ - параметр экспоненциального распределения. Для

моделирования экспоненциального распределения применим метод обратной функции.

y

где γ R[0,1].

F(y)= λ ∫e−λt dt =1−e−λy ,

F(y)= γ 1−e−λy

= γ,

0

1 − γ = e−λy , − λ y = ln(1 − γ) = ln γ ,

так как

если γ R[0,1],

то

и

(1 − γ) R[0,1].

Итак,

y = −

1

ln γ

(3.5.5)

λ

формула, моделирующая показательное распределение. По теореме 3.3

α+1

α+1

1

1

α+1

1

α+1

x = ∑ yi = ∑

∑ln γi = −

∏

(3.5.6)

−

λ

ln γi = −

λ

ln

γi .

i=1

i=1

λ i=1

i=1

Формула (3.5.6) справедлива, если

α- целое. В случае дробного

α

формулу (3.5.6) модернизируют следующим образом:

1

N

γ

+ x

, α

+1 = N + 0.5,

x = −

ln

λ

∏

i

(3.5.7)

i=1

x

= −

1

ln γN +1 cos

2

(2π γN +1 ).

λ

При 0 < α+1 <1 существуют модификации рассмотренного метода. Для значений параметра (α+1) > 1 , можно в комбинации с указанной техникой использовать свойство аддитивности гамма-распределения.

3.6. Лабораторная работа № 3. Моделирование некоторых распределений с помощью базовых случайных величин в пакете MATHCAD

В качестве базового распределения для моделирования многих случайных величин используется стандартное равномерное распределение R[0,1]. Само равномерное распределение может быть получено несколькими способами с помощью линейных, нелинейных и смешанных формул

метода сравнений (см. подразд. 3.3). Рассмотрим формулу
xn+1 = (bxn +θ)mod P -	(3.6.1)
71

частный случай формулы (3.3.3), вырабатывающую последовательность целых чисел, равномерно распределенных в области {0,1,2,...P −1}. Эта

формула используется в программе URAND (Universal RANDom number generator) [22]. Как уже упоминалось ранее, пусковые константы b,θ и P

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

Для формулы (3.6.1) в [11] следующим образом подытожены теоре- тико-числовые ограничения на выбор b,θ и P :

							(b)mod8 = 5,
						1)
					2)		P	< b < P − P ,	(3.6.2)
						100
	θ		1		1
		≈		−		3	≈ 0.21132, θ- нечетноечисло.
3)
	P		2		6

Следует заметить, что формула (3.6.1) использует целые числа, которые в разных ЭВМ имеют разную длину и хранятся по-разному. При выполнении арифметических операций с целыми числами большое значение имеют особенности машинной арифметики [26]. В пакете MATHCAD реально приходится оперировать с вещественными числами, поэтому приведенная ниже программа Urand, взятая из [22] и адаптированная в пакете MATHCAD, нуждается в подробном тестировании, которое автором не проводилось.

Программа реализует формулу (3.6.1) для P = 231 и предназначена для вычисления псевдослучайного числа, имитирующего реализацию случайной величины со стандартным равномерным распределением R[0,1].

ORIGIN :=1

iy := -1023 t := runif (100,0,1)

72

k := Urand(iy) k = 0.48016216

k1:= Urand(k ) k1 = 0.399884106

d := MassUrand(iy,100)

Приведем в заключение еще одну программу, вычисляющую псевдослучайные числа со стандартным равномерным распределением, и входящую в типовое математическое обеспечение ЭВМ фирмы IBM и называе-

мую Randu.

73

ix: =19510 u := Randu(ix,1) u = 0.595 u1:= Randu(ix,100)

Программа Randu реализует ту же формулу (3.6.1), где b = 65539, θ = 0 и P = 231 . Анализ свойств получаемых ею псевдослучайных чисел приведен в работе [22]. Там показано, что в последовательности, вырабатываемой программой Randu, наблюдается крайне высокая корреляция между тремя подряд идущими случайными числами.

Итак, при моделировании стандартного равномерного распределения в пакете MATHCAD при проведении студенческих лабораторных работ можно пользоваться двумя встроенными функциями пакета runif(m,a,b) и rnd(x), а также двумя приведенными подпрограммами.

Другие распределения из списка лабораторной работы № 2 моделируются чаще всего двумя самыми распространенными способами: методом обратной функции и методом суммирования. Далее приведен список распределений и моделирующий каждое конкретное распределение алгоритм.

1. Биномиальное распределение. Метод браковки, описанный в подразд. 3.4 – стандартный способ имитационного моделирования дис-

74

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

n,p – за-

i:=0

r:=rnd

r<p

да x:=x+1

i:=i+1

i<n

нет

X:=x

дано

x:=0

нет

да

- задают n и p по схеме Бернулли;

- получают очередное стандартное равномерно распределенное число

r ;

количество чисел r , которые меньше p , есть случайное число B(n, p), распределенное биномиально.

Использование геометрического распределения. Если p мало, то ме-

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

номиальное случайное число B(n, p)= k −1, где k - минимальное число,

k

такое, что ∑yi > n. yi - геометрически распределенное случайное число.

i=1

Блок-схема этого метода такова:

с:=1/ln(1-p)

да

n,p зада-

r:=rnd

[

]

+1

s:=s+y

X:=x

но

s:=0

s>n

y:= c

lnr

x:=0

нет

x:=x+1

2. Геометрическое распределение. Случайные числа G(p) получа-

ются из случайных чисел, распределенных равномерно на [0,1],						с помо-
щью соотношения
		ln r
xi	=	i			,	(3.6.3)
	ln(1 −		p)

где [...]- целая часть числа. Алгоритм основан на следующей теореме. 75

Теорема 3.4. Если r R[0,1] - базовая случайная величина, то случайная величина x = [ln rln(1 − p)], где [x]- целая часть x , имеет рас-

пределение P(X = x) = p(1 − p)x , x = 0,1,2,...

Доказательство

								0,	x ≤ 0,
	Так как для r R[0,1] по определению								< x <1, то
					F (x) = P(r < x) = x, 0
						r			x ≥1,
								1,
		ln r
	P(X = x) = P x ≤			< x +1 = P{x ln(1 − p) ≥ ln r > (x +1)ln(1 − p)}=
		ln(1 − p)
						− p)x +1 < r ≤ (1 − p)x }=
	= P{ln(1 − p)x +1 < ln r ≤ ln(1 − p)x }= P{(1
	= (1 − p)x − (1 − p)x +1 = (1 − p)x (1 −1 + p) = p(1 − p)x = pq x .

3. Распределение Пуассона. Первый способ моделирования описан в подразд. 3.4 и употребляется когда p мало. Его блок-схема такова.

λ=np,	x:=0	r:=rnd				r:=r-p			r<0		да
n,p – за-											X:=x
дано
											нет
			p :=	pλ				x:=x+1
				x

1) Задают n , p и λ = np так, чтобы n → ∞ , p → 0 ( p < 0.01) .

2)Получают очередное стандартное равномерно распределенное чис-

ло r .

3)Количество чисел r , меньших p , есть случайное число P(λ), рас-

пределенное по закону Пуассона.

Второй способ основан на связи пуассоновского с показательным и эрланговским распределениями и базируется на теореме 3.5.

Теорема 3.5. Случайная величина x , определенная соотношением

	N +1		распределена по закону Пуассона.
x = min N : ∏ rk < e−λ,		N = 0,1,2,...
	k =1

Блок-схема алгоритма выглядит следующим образом.

76

λ – задано	x:=0,							да
	а:=1,	r:=rnd	a:=a*r			a<p
								X:=x
	p:=e-λ
								нет
					x:=x+1

4. Экспоненциальное распределение. Самый распространенный ме-

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

x	= −	1	ln r ,	r R[0,1].	(3.6.4)
i		λ i		i

5.Классическое распределение Вейбулла. Распределение Вейбулла

спараметром α =1 совпадает с экспоненциальным распределением со

средним b = (1λ)1c , поэтому моделируется также методом обратной

функции. Действительно, если	F(x)=1	−e−λxα = r, r R[0,1], то
1 − r = e−λxα , r = e−λxα , − λxα = ln r,	xα = −1 λln r, x = (−1 λln r)1 α .

Таким образом, моделирующая формула имеет следующий вид:

	x	i	= (−1 λln r )1 α ,					r R[0,1].		(3.6.5)
						i		i
6. Распределение Парето. Функция плотности вероятности этого
распределения равна			f (x) = α x0 (x0					x)α+1,	x > x0 , где x0 - параметр по-
ложения, левая	граница				области			возможных		значений (x0 > 0),
α- параметр	формы					(α > 0).			Функция	распределения
F(x) =1 − (x0 x)α,	x > x0 легко обращается, в результате получается мо-
делирующая формула							(1 r )1 α , r R[0,1].
			x	i	= x	0				(3.6.6)
							i	i

7. Распределение Эрланга. Распределение Эрланга ( α +1)-го порядка – это гамма-распределение с целым параметром α (см. формулу (2.5.1)). Графики и формулы оценивания параметров для гаммараспределения применимы также и для распределения Эрланга. Распределение моделируется методом суммирования по теореме 3.3. Основная моделирующая формула имеет вид

77

	α+1		(3.6.7)
xi = −βln	∏ri , ri R[0,1].
i=1

8. Гамма-распределение. I. В том случае, когда α- целое число, моделирование происходит по формуле (3.6.7).

II. Если α- не целое, может быть использован следующий алгоритм. Сначала для −1 < α < 0 выберем r1,r2 и r3 - три независимые случайные величины, равномерно распределенные на отрезке [0,1]. Положим

S	= r1 (α+1)	, S	2	= r1 α	. Если	S	+ S	2	>1 , возьмем вместо r ,r другую
1	1			2		1			1	2

пару таких же случайных величин. Так будем поступать до тех пор, пока не получим S1 + S2 ≤1 . В этом случае случайная величина

x = βS1 ln r3 (S1 + S2 ) будет иметь гамма-распределение с параметрами α и β . Блок-схема этого алгоритма имеет вид

r1 := rnd

α,β - задано r1 := rnd r1 := rnd

	Вычисле-	да		βS	1	ln r
	ние	S1 + S2 ≤1	X :=			3
				S1 + S2
	S1 , S1

нет

Наконец,

для параметра формы α > 0

псевдослучайные числа,

под-

m

чиненные гамма-распределению, дает формула y = x −βln

∏ri , где слу-

i=1

чайные числа x

имеют гамма-распределение

с

параметрами

β и

α1 = α −[α],

m = [α +1]. Здесь [...]- целая часть числа.

Итак,

в

случае нецелого α

моделирующие формулы для гамма-

распределения выглядят следующим образом.

для −1 < α < 0

S

= r1 α+1, S

2

= r1 α,

1

1

2

если S + S

2

≤ 1,

x

i

=

βS1 ln r3

,

1

S1 + S2

(3.6.8)

m

для α1 > 0

yi

= xi − βln ∏ri

,

α

= α − [α],

m

i =1

1

= [α +1], r R[0,1].

i

9. Распределение Коши. Функция плотности вероятности этого распределения приведена в лабораторной работе № 2 (см. формулу (2.7.3)).

78

Функция распределения имеет вид F(x) = 12 + π1 arctg x −λ μ . Методом об-

ращения из нее легко получается моделирующая формула

				1		+ μ, ri R[0,1].
	xi	= λtg π ri −					(3.6.9)
				2

10. Нормальное распределение. Наиболее употребительный метод моделирования нормального распределения – метод суммирования. Моде-

	12
лирующая формула (3.5.3)
	xi = mx + σx ∑ri − 6	, ri R[0,1], дает слу-
	i=1

чайную величину, распределенную нормально с математическим ожиданием M [X ]= mx и дисперсией D[X ]= σ2 . Для моделирования стандартной нормальной случайной величины применяется формула

			12
			xi = ∑ri − 6 .			(3.6.10)
			i=1
11. Логарифмически нормальное (логнормальное) распределение.
Логнормальное распределение имеет				плотность	вероятности	равную
f (x) =	1	e−ln 2 (x m)	2σ2 , где m -	параметр масштаба (медиана), σ-
xσ	2π		случайная величина X		распределена	по лог-
параметр	формы. Если

нормальному закону с параметрами m,σ , то случайная величина Y = ln X подчиняется нормальному закону распределения с математическим ожиданием μ = ln m и средним квадратическим отклонением σ.

Так как X = eY , то моделирующая формула, очевидно, имеет вид

x	= mexp			12	− 6			R[0,1].	(3.6.11)
			σ	r		, r
i				∑ i			i
			i=1

12. χ2 - распределение подробно описано в подразд. 2.1. Для генерирования случайных чисел используются несколько алгоритмов. Самый распространенный алгоритм опирается на определение χ2 - распределения и реализуется формулой

n

xi = ∑u2j , u j N (0,1). (3.6.12)

j =1

Величины xi , полученные по формуле (3.6.12), имеют χ2 - распреде-

ление с n степенями свободы.

79

13.t - распределение Стьюдента. Распределение Стьюдента связано

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

оно совпадает с распределением Коши с параметром положения μ = 0 и

параметром масштаба λ =1. При n → ∞ распределение Стьюдента сходится к стандартному нормальному распределению. При произвольном n для генерирования случайных чисел используется статистика (2.2.1). То-

гда tn = u χ2n n , u N (0,1),		12				χn2	n
		u = ∑ri			− 6,		= ∑ui2 , ri R[0,1]. Мо-
		i =1					i =1
делирующая формула выглядит следующим образом:
	ti	=		ui		,
				n
		1			2
				∑u j
			n				(3.6.13)
				j =1
	12
ui = ∑rk − 6, rk R[0,1].

k=1

14.F -распределение. Как и в двух предыдущих случаях моделиро-

вание псевдослучайных чисел основано на определении основной статистики F -распределения. Случайная величина, имеющая распределение

Фишера, связана с независимыми случайными величинами										χn2	и χn2
										1	2
следующим соотношением	x	n1	,n2	= χ2	n	2	(χ2		n ) . Тогда		модели-
				n			n	2	1
				1

рующая формула имеет вид

	xn1,n2
	12
ui	= ∑rj
	j =1

	n
	n2 ∑1	ui2
=	i =1
	n
	n1 ∑2 ui2		(3.6.14)

i=1

−6, rj R[0,1].

15. Логистическое распределение. Функция плотности вероятности

этого распределения равна f (x) = ( e(x −μ) λ ) , где μ- параметр поло-

λ1 + e(x −μ) λ 2

жения (− ∞ < μ < ∞), а λ - параметр масштаба (λ > 0). Функция распреде80