УДК 519.21
А. В. Мастихин
РЕШЕНИЕ СТАЦИОНАРНОГО ПЕРВОГО УРАВНЕНИЯ КОЛМОГОРОВА ДЛЯ МАРКОВСКОГО ПРОЦЕССА ЭПИДЕМИИ, РАЗВИВАЮЩЕЙСЯ ПО СХЕМЕ Ti + T2 ^ Ti + T3, Ti + T3 ^ Ti, T1 ^ 0
Для трехмерного марковского процесса специального вида решено стационарное первое уравнение Колмогорова для переходных вероятностей. Получено интегральное представление для производящей функции финальных вероятностей. Найдены асимптотики для математического ожидания и дисперсии финального распределения и установлены предельные теоремы.
Определение марковского процесса. На множестве состояний
N3 = {а = («1,^2, а3), а^ а2, а3 = 0,1, 2,... } рассматривается однородный по времени марковский процесс £(£) = (£1(£),£2(£),£3(£)), £ € [0, то), с переходными вероятностями Р^вГвз^ (£) = Р{£(£) = = (в1,в2,в3) | £(0) = (а1,а2,а3)}. Пусть при Д£ ^ 0 переходные вероятности имеют вид (А1 > 0, Л2 > 0, А3 > 0)
P
t1;^ а 3+i)(Ai) = Aia^At + o(At),
P
tS^-vA) = A2a: a3At + o(At),
P
fc^i 3)(At) = Asa: At + o(At),
P
(<а11,а221а33))(At) = 1 - (Aiaia2 + A2aia3 + Asai)At + o(At).
Определим производящие функции (|s1| < 1, |s2| < 1, |s3| < 1)
oo
Fa (t; si, s2, s3) - ^bÄ^f (t)Se1 s22 s33 .
в1,в2,вз=0
Вторая (прямая) система дифференциальных уравнений Колмогорова для переходных вероятностей равносильна уравнениию в частных производных [2, 6]
dFa (t; s) w , ö2Fa (t; s)
-—-- Aii S1S3 - sis^l --+
dt 4 7 ösiös2
, , д2Fa (t; s) Л ( , ÖFa (t; s) (1)
+ Msi - + A3i1 - si)^sr~, (1)
с начальными условиями
Fa(0; si,s2,s3) = sa1 sa2sa3.
Введем экспоненциальную (двойную) производящую функцию [4, 5]
_ ~а1 -а2 ~а3
Т(¿; ^2, гз; 8Ь 82, вз) = > —:—:—г(¿; §1, 82, вз).
а: ,а2,аз=0
Первая (обратная) система дифференциальных уравнений Колмогорова для переходных вероятностей процесса £(£) приводится к виду
дТ = Л г г ( 92Т д2Т \ + Л г г ( дТ _ д2Т \ +
+ Лзг^Т - , Т(0, г1, ¿2, гз; 81, 82, ) = е5121+5222+5323. (2)
Интерпретация процесса. Задача о финальных вероятностях.
Марковский процесс £(£) может быть интерпретирован как модель эпидемии [1,2, 6], а именно как модель распространения инфекции с двумя стадиями заболевания. Процессу соответствует схема взаимодействий [2, 6]:
Т1 + Т2 ^ Т1 + Тз, Т1 + Тз ^ Т1, Т1 ^ 0, (3)
где частицы типа Т1 — зараженные особи (источники инфекции); частицы типа Т2 — здоровые особи (восприимчивые к инфекции, не имевшие пока контактов с зараженными); частицы типа Тз — особи, имевшие один контакт с зараженными (ставшие носителями инфекции). Здоровая особь после двух контактов с зараженным удаляется из популяции. Состояние (а1, а2, аз) означает наличие а1 частиц типа Т1, а2 частиц типа Т2 и аз частиц типа Тз. Через случайное время ть с распределением вероятностей Р{т1 < ¿} = е-а1а2А:^, пара частиц типа Т1 и типа Т2 взаимодействует и независимо от других частиц превращается в частицу типа Т1 и частицу типа Тз .Процесс переходит в состояние, соответствующее вектору (а1, а2 — 1, аз + 1). Через случайное время т2, с распределением вероятностей Р{т2 < ¿} = е-а1азА24, пара частиц типа Т1 и типа Тз взаимодействует и независимо от других частиц превращается в частицу типа Т1. Процесс переходит в состояние, соответствующее вектору (а1, а2, аз — 1). Кроме того, через случайное время тз, с распределением вероятностей Р{тз < ¿} = е-а:Аз4, частица типа Т1 умирает и процесс переходит в состояние, соответствующее вектору (а1 — 1, а2, аз). Случайные величины т1, т2 и тз независимы, в состоянии (а1, а2, аз) процесс находится случайное время т = шт{т1, т2, тз}.
Предложенный марковский процесс рассмотрен в работе [2] в случае Л1 = Л2 = 1,Лз = В [2] методом преобразования Лапласа решено
уравнение (1). Определим финальные вероятности для поглощающих состояний (0,72,7э), 72,7э = 0,1, 2,...,
те
(«1 ,а2|аз) = 1- р (а1,а2|а3)(.) „(а1,а2,аэ) = 1
%72,7З) = I™ Р(°.72,73) 3(0,72,7З)
72,73 =0
Для финальных вероятностей в [2] получено соотношение
а2 «2-72 а2
EV „(а1,а2,аз)„72 „73 _ V^ ^__(е _ i\ Y2V
0 ^ д(0^2,^з) „2 _ (72 + ^)Y2+«1 («2 - 72) ^ 1
72 =0 73 =0 72 =0 4 х 4 '
72 lo Ч-2 { (72 + (—-1)
-££... £ 1—гг^
10=011=0 га1-1=0 а1 1
Однако это выражение малопригодно для исследования асимптотических свойств рассматриваемого марковского процесса.
В настоящей работе процесс £(£) исследуется предложенным в работах [4, 5] методом экспоненциальной производящей функции. Получены интегральное представление для производящей функции финальных вероятностей, асимптотики для математического ожидания и дисперсии финального распределения при а2 ^ то и предельные теоремы. Такие теоремы "порогового" типа позволяют определить пороговую численность инфицированных особей, при превышении которой принято говорить о начале эпидемии [2, 6]. Согласно работе [2] исследован случай А1 = Л2 = 1, Л3 = случай А1 = А2 будет рассмотрен в работе, готовящейся к печати.
Стационарное первое уравнение Колмогорова. Вводим производящую функцию финальных вероятностей
Ф(а1,а2,аз) („2, „з) _ £ fe^^ „Г, I „2 | < 1, |„з| < 1,
72,73 =0
и двойную производящую функцию
а1 а2 аз
$(Zl,Z2,Z3i S2,S3)_ ^ 1 | 2 | 3 | Ф(аьа2,аз)(„2, S3).
aila2i«31
а1,а2,аз=0
Аналитичность функции Ф(а1,а2,а3)(^2,53) при |з2| < 1, |з3| < 1 устанавливаем, исходя из неравенства |Ф(а1,а2,а3)(^2, в3)| < 1 в рассматриваемой области.
Аналитичность функции Ф(г1, -2, -3; 82, 83) при любых г1, -2, -3 следует из неравенства
|Ф(-1,-2,-3; 82,^3)1 <
< Е к1|а1 ^!2У3 |Ф(в1'«а'«з)(*2,*3)|< е|21|+|22|+|231. ^ п а1!а2!а3!
а1,а 2,а з=0
Из первого уравнения (2) аналогично теореме 1.4 работы [5] получаем, что функция Ф(г1, -2, -3; 82, 83) удовлетворяет стационарному уравнению
( д2Ф д2Ф \ (дФ д2Ф \ ( дФ \
-1-Ч д-д-3 - д-д-2) + -1Ч д-1- д-д-3) + КФ - д-") =0,
или, после преобразований,
д2Ф д2Ф дФ
(-2 - -3) д-д-3- -2 д-д-2- - -3) д-1+=0 (4) Получим граничные условия. Пусть = 0. Очевидно, что
9(0,а2,а3) = 1 и в = 0 пРи «2 = ^2 или а3 = 73. СлеДовaтельно,
ГС -"2-а3
Ф(0,-2,-3; 82,83)= У^ 2 3 Ф(0'О2'Оз) (®2,83) =
^ п «2^3! 4
«2 -ОД
\ ' -2 -3 8а2 8аз = ^222+3323 «2!«3! 2 3 .
Пусть -2 = 0, -3 = 0. Для состояний (а1,0, 0), а1 = 0,1, 2 ..., финальные вероятности определяются следующим образом: 0'°э)0) = 1 и = 0 при 72 = 0 или 73 = 0. Следовательно,
-«1 ^ -а1
Ф(-1, 0, 0; 82,83) = V Ф(аь0 ,0)(82, 83) = V = е21.
' а1! а1!
а1=0 1 а1=0 1
так, уравнение (4) рассматривается при условиях
Ф(-1, 0, 0; 81,82,83) = е21, Ф(0, -2, -3; 81,82,83) = ев222+в323. (5)
В работе получено решение задачи (4), (5), удовлетворяющее условию аналитичности. Вопросы существования и единственности решений
уравнений вида (4) не являются целью настоящей работы, отметим лишь, что они сложны и мало изучены.
Замена переменных. Функция Римана. Поставленная задача решения уравнения (4) с тремя переменными и условиями (5) сводится к граничной задаче решения гиперболического уравнения с двумя переменными.
Рассмотрим замену переменной х = г1, у = е-2з/22, ( = г2е2з/22. Тогда г1 = х, г2 = £у, г3 = — £у 1п у. После вычислений получим уравнение в частных производных для функции Ф(х, у) = Ф(х, £у, —Су1пу; «2,^3):
Фху +(- + С 1п ^ Фх — -ф = 0, (6)
уу
с условиями ф(х, 0) = ех, ф(0, у) = е^у(в2-в31пу). Запишем уравнение (6) в виде
дх (Ф у+(у+< 1п у)Ф) — уФ=0-
Решая обыкновенное дифференциальное уравнение с разделяющимися переменными ( )
Ф у + + с 1п у)ф = 0,
приходим к замене ф(х, у) = и(х, у)у-^е-у^(1пу-1). Подставляя последнее выражение в уравнение (6), получаем для функции и(х, у) уравнение
мЖу — = 0 (7)
У
с граничными условиями
и(х, 0) = 0, и(0,у) = у^((в2-1)-(в3-1)1пу). (8)
Полученная задача Гурса (7),(8) решается методом Римана [7, 8]. В общем случае для гиперболического уравнения Ь(м) = иху + + аих + + си = 0 решение граничной задачи и(х, у0)=ф(х), и(хо,у) = ^(у), ф(х0) = ^(у0), дается формулой [7]
и(х, у) = Л(х, у0; х, у)ф(х) + Л(х0, у; х, у)^(у) — Л(х0, у0; х, у)ф(х0)+
гх г д 1
у0)Л(^ у0;x, у) — ^ту0; х у)
+
J xo
dt
"У
+у
д
a(x0, t)R(x0, t; x, y) - ^R(x0, t x, y) ^(t)dt, (9)
где Л(х0, у0; х, у) — функция Римана. По определению функция Римана Л(х0, у0; х, у) является решением сопряженного уравнения Ь*(м) =
= — (ам)ж — (Ьи)у + си = 0 относительно переменных х и у. В случае уравнения (7) рассматриваемое уравнение совпадает с сопряженным. Кроме того, требуется выполнение следующих соотношений:
Д(хо,уо; Хо,у) = Дж(хо,уо; х,уо) =
= Дхо (хо,у; х,у) = Дуо (х,уо; х,у) = 0,
Д(х, у; х, у) = Д(хо, уо; хо, уо) = 1. Лемма 1. Функция Римана для уравнения (7) имеет вид
Д(хо, уо; х, у) = 7о (—и(х — хо)1п у^, (10)
где 7о — функция Бесселя нулевого порядка.
Доказательство. Функцию Римана для уравнения Ь(м) = 0 ищем в виде ряда [9]
т->/- \ ^ VJ(x,y)(x - x0)j(y - yo)J R(xo, yo; x, y) = -j-.
j=0 j 'j' Нулевой коэффициент, следуя [9], находим из соотношения
(x , y)
lnv0(x,y) = J (adX + bdY) = 0,
(xo,yo)
таким образом, v0 (x,y) = 1. И нтеграл берется по отрезку, соединяю -щему точки (x0, y0) и (x, y); в полярных координатах имеем x = x0 + + r cos y = y0 + r sin 9 для граничных точек и X = x0 + s cos 9, Y = y0 + s sin 9 для внутренних точек из этого отрезка, где 9 = const, s £ [0, r]. Рекуррентная формула для вычисления коэффициентов ряда получается при подстановке ряда в сопря енное уравнение и имеет вид [9]
r
V (x,y) = -j J sj-1L*(Vj-i(X, Y))ds. 0
Вычисление v1(x, y), v2(x, y) составляет основание индукции по j:
( r ( r
L*(Vo) = -- Vl = f --ds = d(y0 + s sin 9) -
y' r j y r J (y0 + s sin 0) sin в
0 0
r
- ln(y0+s sin в)
r sin в
- - (ln(y0+r sin в) -ln У0) = —-— ln —;
0
r sin в y - y0 y0
L>i) =
u u2 y u2 y0 + s sin 0
=--v1 = —--— In — =----- In-,
y (y — yo)y yo s sin %o + s sin У0
r
v = í___ln yo +s sin 0ds = ln2 y_
2 r2 J s(yo + s sin 0) sin 0 yo (y - yo)2 yo '
o
y
Предположим, что коэффициент v = --— lnj — вычислен.
(y - yo)j yo
Определим vj+i:
) = _ lnj ^ = _ j bJ ((yo + s sin 0)/yo)
j (y - yo)jy yo (yo + ssin0)(s sin0)j '
v = _ j + 1 /-sjj lnj((yo + s sin0)/yo) ds = ^ lnj+i y_
Vj+1 rj+V sj(yo + s sin0)(ssin0)j s (y - yo)j+ n yo ' o
Окончательно, функция Римана имеет вид
Rxo,*; x,y) = Í: Vj (x,y><x -, X f <y - y°'j =
j=o j 'j'
= £ ln(y/yo)(x - xo))j = Jo^-^x - xo)ln у-)' j=o j 'j' V y°
Лемма доказана.
Интегральное представление для производящей функции
(s2,s3).
Теорема 1. Производящая функция финальных вероятностей имеет вид
г ™
$(ai,a2,aa)(s2,s3) = 3x^V^X
"i а2 г "1 аз
X e-x((s2 - 1) + x(s3 - 1)) + 1 e-x(s3 - 1) + 1 dx, (11)
где ^ > 0, a1 = 0.
Доказательство. Для уравнения (7) рассмотрим решение uo(x,y) задачи Гурса (xo > 0, yo > 0) с граничными условиями
uo(xo,y)= ^(y), uo(x,yo) = 0(x),
где ^(у) = (у — уо)^ес(у-у°)((в2-1)-(в3-1)1п(у-у°», ф(х) = 0. Формула (9) приобретает вид
у
[ д
мо(х,у)= ^(у) — д^ ^^ (хо, х,уМс)^.
уо
Используя метод интегрирования по частям и учитывая равенства
Д(хо,уо; х, у) = 7о(0) = 1, получим
у
у Г д-^(С)
ио(х,у) = ^(у) — ^(¿)Д(хо,£; х,у) + / Д(хо,С; х,у)^С =
у° J дс
У0
y
= f dm
J dt
y0
R(x0, t; x, y)dt.
Рассмотрим предел при x0 ^ 0, y0 ^ 0, получим
y y
u(x,y) = f R(0,t; x,y)dt = f ^e-Zi((s2-1)-(s3-1)1n 4)x
0
x
^ + tZ ((S2 - 1) - (S3 - 1)(ln t +1))
X
х 7о( 2а/^х 1п (12)
Возвращаясь к переменным г1, ¿2, и функции Ф(г1, ¿2, ¿3; 82, 83) = = м(г1, егз/г2 )е^2з/22+22+23, получаем выражение для экспоненциальной производящей функции
Ф(^1,^2,^з; 82,83) =
/z2+Z2ez3/z2т((s2-1)-(s3 —1) 1nt)+Z2+Z3 ^
0
X
^+TZ2eZ3/z2 ((S2-1)-(S3- 1)(ln T+1))] Jo(ln T + Z3))dT. После замены t = e-z3/z2-x получаем
Ф(^1, Z2, Z3; S2, S3) =
e-^x+22(e-x(s2-1)+x(s3-1)+1)+23(e-x(s3-1)+1)x
OG
0
x((s2-1) + (x-1)(s3-1)))+z3e x(s3-1) Jo(dx.
(13)
Совершенный вше предельный переход нуждается в обосновании; однако непосредственной подстановкой выражения (13) в уравнение (4) и проверкой условий (5) убеждаемся, что (13) является решением задачи (4), (5). Единственность устанавливается, исходя из аналитичности решения (ср. [9, 5]).
Учитывая определение производящей функции ¿2, г3; 52, 53) и разложения в ряды
е* - V 7 (.)-V (^ (У2)*'
е — 2.. п, ТУТ! ,
п=0 ^=0 '
запишем разложение подынтегральной функции в ряд по г1 и г3 . Имеем
52, 5з) -
сю
„ а1!а1! азЧ
«1=0 аз=0 о
„ + ^2(е-ж((82 - 1) + (х - 1)(5з - 1)))(е-х(5з - 1) + 1)аз +
X
+ аз(е-ж(вз - 1))(e-x(s3 - 1) + 1)аз-1
dx.
алее, разбивая интеграл на сумму двух интегралов соответственно слагаемым, стопим в квадратных скобках, и производя интегрирование по частям во втором интеграле, получаем
оо
zal za3 „ai
^^^3)= £ „1,(0,1 - 1)!o3i X
ai,a3=0 4 '
ai-1e-^x+z2(e x(s2 —1)+x(s3 —1)+1) (g-x (S3 - 1) + 1)a3dx.
x
Представляя экспоненту под интегралом в виде ряда по г2, получаем утверждение теоремы. Теорема доказана.
Замечание. При начальном состоянии процесса (а1, 0, а3) имеем марковский процесс эпидемии Вейса со схемой взаимодействий [3]
Т1 + Тз ^ Т1, Т1 ^ 0.
Соответственно, при а2 — 0 выражение (11) совпадает с производящей функцией финальных вероятностей процесса эпидемии Вейса [5].
Асимптотические свойства финального распределения. Для
марковского процесса £(£) частицы типов Т2 и Т3 называются финальными [10]. Обозначим п2а1,а2,аз), Пза1'а2'аз) случайные числа финальных частиц, которые останутся после того, как процесс выродится, т.е. не останется частиц типа Т1. Совместное вероятностное распределение величин п2а1'а2 аз), п3а1'а2'аз) определяется производящей функцией (13). Одномерные распределения финальных вероятностей случайных величин п2 и п3 задаются производящими функциями
вероятностей Ф(а1,а2,аз)(82,1) и Ф(а1,а2,а3) (1, 83). Для факториальных моментов из формулы (11) при а2 ^ то получаем
К, (а1,а2,аз )_ дФ(а1,«2,аз) (1, 1)_ ( » \а М2 = -я- = а2 I -—7 I ,
д82 V +1/
м (а1,а2,аз) = дФ(«1,«2,аз)(1, 1) «1 «2 ( ^ Ча1 д2ф(а1,а2 ,аз)(1, 1) , 1А( ^ У1
-021-= а2(а2 — 1Ч^г) ,
д 2 Ф(а1
(1, 1) «2 («2 — 1)«1(«1 + 1) ( ^ Ча1
ds2 у + 2 + 2
pw (01,02,03) _ д Ф(ах,а2,03) (1, 1) | дФ(а1,а2,а3) (1, 1)
= я„2 +
>2 dS2
(дФ(а1,а2,а3)(1, 1)\2 ( ч2 Г( ^ V1 ( V Ч2а1 ]
(—^—) ~(а2) Кут^) - (ул)
а1 - 1 ( у \а1 а1 ( у \2а1
Dn3"1,"2,a3) „ «1(«2)2[-01-i^f01 - -JLJ)
13 П ' Цу + 2)2Vу + 27 (у + 1)2Vу + 1/
Теорема 2. Пусть x G [0,1]. Тогда
(п(а1 ,а2 ,а3) а1-1 ( in x)n
lim pi %- < Л = x^ ^ (-у 1ПХ)
«2^^ I а2 I ' п!
4 ' п=о
Доказательство. Преобразование Лапласа от функции распределения случайной величины п2а1,а2,аз)/а2 имеет вид [11]
ч (а1,а2 ,аз) /
М(б-Ап2 1 2 з/а2 ) = Ф(а1,а2,аз) (б-"7"2 , 1) =
/7а1 Г™ У I
Х
(«1 - 1)! Л)
(\ а2
е-ж((е-Л/а2 - 1) + 1)j dx.
При а2 ^ то получаем
$(ai,a2,a3)(e-Va2 , 1) ~ ^ ^J Xai1 - ^g-*) 2 dx
1 0 2
(«1 - 1)^0
После замены y = g-x имеем
„ai Г1
„ 1 (-ln y)ai-y-1e-Äy dy.
. ai ^
„ ' xai-1g-^xg-Äe-xdx.
(«1 - 1)! Л
Полученное выражение является преобразованием Лапласа плотности
иа1
распределения f (у) — --у)а1-1у^-1 для случайной вели-
(а 1 1)у
чины, распределенной на отрезке [0, 1]. По теореме непрерывности [11] имеем сходимость функций распределений
Л- < х — f (у)йу.
«2 J ,/0
алее вычисляем
Т^Г [Х(- 1пуГ-У-Чу — х^ ^^^.
(«1-1)! л п=0 п!
Теорема доказана.
Теорема 3. Пусть х Е [0,1]. Тогда
{„(а1,а2,аз) /• X
П3- < х — Л(у)^у,
«2 J У0
где преобразование Лапласа от плотности распределения вероятностей f 1 (у) при Л > 0 имеет вид
Гс „«1 Гс
е-АуЛ(у)ф — —„„-гт ха1-1 е-^Хе-Ахе хЖх
(«1 - 1)Uo
0
Теорема 4. Пусть xi, ж2 G [0,1]. Тогда
{(а1,а2,аз) п(а1,а2,аз)
П2- < xi, - < ж2 }> =
а2 а2
ЛЖ2
/ / /2(Уl,У2)dУldУ2, '0 J 0
где двойное преобразование Лапласа [12] от двумерной плотности распределения вероятностей ^(у^у2) имеет вид (Л1 > 0, Л2 > 0)
/ / е-Л1У1-Л2 y2 /2(yi,y2)dyidy2 = 'о J0
,а1 /'•те
ж
на _
--а1-1е-^же-(Л1+Л2ж)е x
(«1 - 1)^0
Доказательство теорем 3 и 4 аналогично доказательству теоремы 2. Явный вид для плотностей вероятностей ^(у^ и ^(уь у2) может быть найден в виде рядов по специальным функциям.
те />те
СПИСОК ЛИТЕРАТУРЫ
1. Эпидемии процесс // Математическая энциклопедия. Т. 5. - М.: Советская энциклопедия, 1985. - Кол. 1008.
2. G a n i J. Approaches to the Modelling of Aids // Lecture notes in biomathematics. V. 86. Stochastic processes in epidemic theory. - Heidelberg: Springer, 1990. - P. 145154.
3. Weiss G. Onthe spread of epidemics by carries//Biometrics. - 1965.-V. 21, № 2. -P. 481-490.
4. Севастьянов Б. А., Калинкин А. В. Ветвящиеся случайные процессы с взаимодействием частиц // Докл. АН СССР. - 1982. - Т. 264, вып. 2. - C. 306-308.
5. Калинкин А. В. Финальные вероятности ветвящегося процесса с взаимодействием частиц и процесс эпидемии // Теория вероятн. и ее примен. - 1998. -Т. 43, вып. 4. - С. 773-780.
6. К а л и н к и н А. В. Марковские ветвящиеся процессы с взаимодействием // УМН. - 2002. - Т. 57, вып. 2. - С. 23-84.
7. Бицадзе А. В., Калиниченко Д. Ф. Сборник задач по уравнениям математической физики. - М.: Наука, 1985. - 312 с.
8. Курант Р. Уравнения с частными производными. - М.: Мир, 1964.
9. C o p s o n E. T. Partial Differential Equation. - Cambridge: Cambridge Univ. Press, 1975. -280 p.
10. Севастьянов Б. А. Ветвящиеся процессы. - М.: Наука, 1971. - 436 c.
11. Феллер В. Введение в теорию вероятностей и ее приложения. Т. 2. - М.: Мир, 1984.
12. Д и т к и н В. А., Прудников А. П. Операционное исчисление по двум переменным и его приложения. - М.: Физматгиз, 1958.
13. Мастихин А. В. Функция Римана для стационарного уравнения марковской эпидемии // Обозрение прикл. промышл. матем. Cер. вероятн. и статист. - 2003. - T. 10, № 2. - C. 502.
Статья поступила в редакцию 17.02.2005