close

Вход

Забыли?

вход по аккаунту

?

Vorobiev

код для вставкиСкачать
Министерство образования и науки российской федерации
Государственное образовательное учреждение
высшего профессионального образования
Санкт-Петербургский государственный университет
аэрокосмического приборостроения
С. Н. Воробьев, Н. В. Гирина,
И. В. Лазарев, Л. А. Осипов
СТАТИСТИЧЕСКОЕ МОДЕЛИРОВАНИЕ
ИНФОРМАЦИОННЫХ СИСТЕМ
Часть II
Учебное пособие
Допущено Учебно-методическим объединением вузов
по университетскому политехническому образованию
в качестве учебного пособия для студентов высших учебных
заведений, обучающихся по направлению подготовки
230400 Информационные системы и технологии
Санкт-Петербург
2011
УДК 519.7
ББК 22.18
В75
Рецензенты:
доктор технических наук, профессор М. О. Колбанев;
доктор технических наук, профессор С. А. Яковлев
Утверждено
редакционно-издательским советом университета
в качестве учебного пособия
Воробьев, С. Н.
В75Статистическое моделирование информационных систем:
учеб. пособие. Ч. II / С. Н. Воробьев, Н. В. Гирина, И. В. Лазарев, Л. А. Осипов. – СПб.: ГУАП, 2011. – 64 с.: ил.
ISBN 978-5-8088-0615-3
Вводятся марковские модели гауссовых последовательностей.
Теория гауссовых марковских последовательностей первого и высших порядков излагается на уровне основ теории вероятностей и математической
статистики, а также соотношений линейной алгебры. Устанавливаются необходимость и достаточность свойства многодиагональности матрицы точности последовательности. Разрабатываются алгоритмы аппроксимации немарковских процессов марковскими последовательностями минимального
порядка и генерирования марковских последовательностей.
Применение марковских моделей иллюстрируется примерами расчетов
и моделирования плотностей распределения времени пересечения случайными траекториями неслучайных уровней и оценивания временного положения
импульсных сигналов.
Плотность распределения времени пересечения марковского процесса
первого порядка с произвольным дифференцируемым уровнем обобщается на
общий случай стационарного гауссова процесса и достаточно крутого уровня.
УДК 519.7
ББК 22.18
ISBN 978-5-8088-0615-3
©Санкт-Петербургский
государственный университет
аэрокосмического приборостроения
(ГУАП), 2011
©С. Н. Воробьев, Н. В. Гирина, И. В. Лазарев,
Л. А. Осипов, 2011
Введение
В множестве случайных процессов x(t) выделяют класс марковских процессов [1]. Аппарат марковских процессов может оказаться
весьма полезным при разработке информационных систем, так как
марковские модели широко применяются специалистами различных научно – прикладных направлений: в автоматическом регулировании, в расчетах надежности и массовом обслуживании, в биологии и медицине, и др.
Марковские процессы имеют следующее отличительное свойство. Пусть заданы n моментов времени t1 < t2 < ... < tn. Если последнее значение процесса xn = x(tn) зависит только от k предыдущих
xn – 1, xn – 2, ..., xn – k, процесс x(t) называется марковским процессом
k-го порядка. Математическая запись этого условия:
f (xn xn−1, xn−2,..., x1 ) = f (xn xn−1, xn−2,..., xn−k )
– условная плотность распределения включает только k предыдущих значений. Процессы первого порядка (k = 1) обычно называют
марковскими – без указания порядка.
Марковскими процессами называют процессы непрерывного
времени, мгновенные значения которых принимают континуум
значений; цепями Маркова – процессы, дискретизированные во
времени и по амплитуде. Марковские процессы и цепи описываются специальными математическими средствами (уравнениями Фоккера – Планка – Колмогорова, уравнениями Маркова) [1]. В настоящем пособии рассматриваются марковские последовательности X –
процессы, дискретизированные во времени с интервалом Δ. Мгновенное значение xi = x(ti) – отсчет процесса – непрерывная случайная
величина с плотностью распределения f(xi), i = 1, ..., n, n – длина
вектора отсчетов X. Такое описание мгновенного значения является
хорошей моделью, например, в системе MATLAB, где вычисления
выполняются с двойной точностью: под мантиссу и порядок отводятся 8 байт памяти, так что диапазон представления чисел (по модулю) в формате LONG E [2]:
realmin = 2,225073858507202 × 10–308 –
– 1.797693134862316 × 10308 = realmax.
Если f(xi) = f(x), то X – стационарная последовательность (вектор
отсчетов стационарного процесса). В радиоэлектронике сигнал на
входе принято описывать суммой
3
z(t) = x(t) + s(t) ∈ Ν(s(t), R (τ)),
x(t) – шум (стационарный нормальный процесс с функцией корреляции R(τ), s(t) – информационный сигнал. Соответственно вектор
отсчетов
Z = X + S ∈ Ν(S, B),
B – корреляционная матрица, а составляющие вектора распределены с плотностью
 1

f (X) = (2π)−n/2 (det B)−1/2 exp − (X − M)T B−1 (X − M). (1)
 2

Параметры n-мерной нормальной плотности (1) – корреляционная матрица B отсчетов шума и вектор математических ожиданийотсчетов сигнала MT = ST = [s1 s2 ... sn]. Корреляционная матрица –
матрица корреляционных моментов kij = M[xixj], i, j = 1, ..., n, которые задаются соответствующими отсчетами функции корреляции.
Например, значения функции корреляции стационарного нормального марковского процесса (первого порядка)
R (τ) = σ2 exp(−α τ ), (2)
взятые с интервалом Δ = 0.5 при α = 1
R (τ) = σ2 [1.0000 0.6065 0.3679 0.2231 0.1353],
задают матрицу
1.0000

0.6065

B = σ2 0.3679
 0.2231

0.1353

0.6065
1.0000
0.6065
0.3679
0.2231
0.3679
0.6065
1.0000
0.6065
0.3679
0.2231
0.3679
0.6065
1.0000
0.6065
0.1353

0.2231

0.3679 .
0.6065
1.0000 
Следует заметить, что по корреляционной матрице невозможно
увидеть, является ли она матрицей марковской или немарковской
последовательности.
Матрица
D = B−1 – обратная корреляционной, называется матрицей точности [3].
4
(3)
Из n-мерной нормальной плотности (1) можно получить условную плотность распределения нескольких первых или последних
составляющих вектора X при фиксировании остальных его составляющих.
Если вектор X разбить на две части:
T
T
XT = [X1T ; XT
2 ], X1 = [x1,..., xk ], X2 = [xk+1,..., xn ],
то этому разбиению соответствуют разбиение вектора математических ожиданий
T
T
MT = [M1T ; MT
2 ], M1 = [m1,..., mk ], M2 = [xk+1,..., xn ]
и разбиение корреляционной матрицы на блоки
 B11 B12 
, (4)
BX = 
 B21 B22 


B11, B22 – корреляционные матрицы первых k и последних n – k отсчетов. При этом условная плотность распределения f(X2|X1) записывается в виде (1) с условной корреляционной матрицей
Bn−k = B22 − CB12 (5)
и вектором условных математических ожиданий
Mn−k = M2 + C(X1 − M1 ), (6)
−1
C = B21B11
.
(7)
Математическое ожидание M[X] можно положить равным нулю,
что позволяет считать последовательность стационарной и упрощает формулу (6):
Mn−k = CX1.
Соотношения (5) – (7) описывают любые нормальные последовательности, в том числе марковские.
5
1. Гауссовы марковские последовательности первого порядка
1.1. Матрица точности стационарной марковской последовательности
Корреляционная матрица B вектора Xn1 (размерностью n отсчетов), полученного дискретизацией стационарного марковского процесса первого порядка имеет вид [4]
 1
β β2


 β
1 β

2 2
Bnn = σ β
β 1


 ... ... ...

βn−1 ... β2

... βn−1 

... βn−2 

... βn−3  , 
... ... 

β
1 
(1.1)
β = k(∆) – коэффициент корреляции между соседними значениями
вектора Xn1, –1 ≤ β ≤ 1. Частный случай: β = e–α|∆| – гауссов процесс
с экспоненциальной функцией корреляции (2) – марковский пер­
вого порядка согласно первой теореме Дж. Дуба [1]. Действи­
тельно, пусть матрица (1.1) подвергается такому разбиению (4), что
B22 = b22 = σ2 (выделяются последняя строка и столбец), то матрица
B11 имеет размерность (n – 1) × (n – 1), матрица B21 – размерность
1 × (n – 1), матрица B12 – размерность (n – 1) × 1. Матрица, обратная
блоку B11, имеет вид
1
−β
0
0
... 0 


−β 1 + b2
−β
0
... 0 



σ−2  0
−β
−β
1 + b2
... 0 
−1

, B11 =

1 − β2  0
−β
0
1 + b2 −β 0 


...
...
...
... ...
 ...


 0
−β 1 
0
0
...
(1.2)
−1
B11 = I.
в чем можно убедиться, перемножая B11 и B11
:×произведение
первой строки на первый столбец равно
6
1
 
−β
  1 − β2
1 
n−2  
0 =
= 1,
β
β
1
...



   1 − β2
1 − β2 
 ... 
0
 
первой строки на первый столбец
 −β 


1 + β2 


3
3
1 
n−2  
 = −β + β + β − β + ... = 0
1
...
β
χ
−
β


 
1 − β2 
1 − β2

 ... 


 0 
и т. д., т. е.
−1
B11
×B11 = I.
Важно подчеркнуть свойство трехдиагональности матрицы (1.2).
Условное математическое ожидание (6) для марковской после­
довательности первого порядка – число (k = n – 1). Блок B21 =
= σ2[βn – 1 βn – 2 ... β], так что матрица (7)
−1
C = B21B11
= [0 0 ... 0 β ], (1.3)
здесь становится вектор-строкой и имеет один последний ненулевой
элемент, равный b, предыдущие n – 2 элемента равны нулю. Зависимость условного математического ожидания от одного предыдущего значения (последнего) есть одно из определений марковского процесса первого порядка [1].
Условная корреляционная матрица (5) также вырождается в число – условную дисперсию последнего значения xn при фиксированных предыдущих значениях.
Вследствие симметрии корреляционной матрицы
f (xn x1,..., xn−1 ) = f (xn xn−1 ),
f (x1 x2,..., xn ) = f (x1 x2 )
– первое значение марковского процесса зависит только от одного
следующего.
7
Пример 1.1. Корреляционные матрицы (1.1) пяти отсчетов стационарного нормального марковского процесса с единичной дисперсией и коэффициентом корреляции β = ±0.7
 1
±0.7

 ±0.7
1

±0.7
B =  0.49
±0.343
0.49

 0.2401 ±0.343

0.49 ±0.343 0.2401 

0.49
±0.7
±0.343

1
0.49  .
±0.7
1
±0.7
±0.7 
0.49
1 
±0.7
Разбиению вектора XT = [x1 x2 x3 x4| x5] соответствуют матрицы (4)
 1
 0.2401 
±0.7 0.49 ±0.343




 ±0.7
±0.343
1
±0.7
0.49 
 , B12 = 

B11 = 
 0.49  ,
±0.7
1
±0.7 
 0.49


±0.343 0.49 ±0.7
 ±0.7 
1 



B21 = [0.2401 ±0.343 0.49 ±0.7 ], B22 = 1.
Матрицы, обратные матрицам B и B11, трехдиагональны:
 1.9608 ±1.3725

0
0
0


±1.3725 2.9216 ±1.3725

0
0


,
0
0
±1.3725 2.9216 ±1.3725
D = 



0
0
1
3725
2
9216
1
3725
±
.
.
±
.



0
0
0
±1.3725 1.9608 

 1.9608 ∓1.3725

0
0




∓1.3725 2.9216 ∓1.3725
0
−1

.
B11 = 

∓
∓
0
1
.
3725
2
.
9216
1
.
3725



∓1.3725 1.9608 
0
0

Вектор (1.3)
−1
C = B21B11
= [0 0 0 ±0.7 ]
задает условное математическое ожидание (6): при M[X] = 0 (n = 5,
k = 4)
8
 x1 
 
 x2 
mx5 x1, x2, x3, x4 = C   = ±0.7x4
 x3 
x 
 4 
– условное математическое ожидание значения x5 зависит только от
одного предыдущего значения x4, причем равно произведению случайного x4 и коэффициента корреляции β. Условная дисперсия последнего значения равна
σ2x5 x1, x2, x3, x4 = σ2x5 − CB12 = 1 − (0.7)2 = 0.51.
При разбиении вектора XT = [x1|x2 x3 x4 x5] матрица B11 становится матрицей B22, матрица B12 – матрицей B21, так что
mx1 x2, x3, x4 , x5 =± 0.7x2, σ2x1 x2, x3, x4 , x5 = 0.51.
Для двух последних значений xn – 1, xn при фиксированных
x1, ..., xn – 2 блоки
 1
β β2


 β
1 β

2 2
B11 = σ β
β 1


...
...
...


βn−3 ... β2

... βn−3 

... βn−4 

... βn−5  ,

... ... 

β
1 
βn−2 βn−3 ... β 
 , B = BT , B = σ2 1 β .
B21 = σ2 
 12
21
22
β 1
 βn−1 βn−2 ... β2 


Матрица (7)
0 0 ... 0 β 

C = 
2
0 0 ... 0 β 
определяет зависимость вектора условных математических ожиданий (6) двух последних значений от одного предыдущего фиксированного xn – 2:
9
 βxn−2 
. Μ  xn , xn−1 xn−2, xn−3,..., x1  =  2

β
x
n−2 

(1.4)
Условная корреляционная матрица (5)
1 − β2 β − β3 

β 
 = σ2 (1 − β2 )1
.
Bn−k = B22 − CB12 = σ2 

β 1 + β2 
β − β3 1 − β4 
В условиях примера 1.1 матрица
 1.9608

−1 
B11
=  ∓1.3725


0

∓1.3725
0

2.9216 ∓1.3725
∓1.3725 1.9608 
по-прежнему трехдиагональна, матрица
0 0 ±0.7

C=
0 0 0.49 


определяет условное математическое ожидание (1.4)
±0.7x3 
.
Μ  x5, x4 x3, x2, x1  = 
 0.49x3 


Условная корреляционная матрица при σ2 = 1
 0.51 0.357 
.
B2 = B22 x3, x2, x1 = 
0.357 0.7599


В общем случае при фиксированных первых k значениях последовательности (n-k) × k – матрица C
 0 0 ... 0
β 


 0 0 ... 0
2 
β


(1.5)
C=

... ... ... ... ... 


 0 0 ... 0 βn−k 


определяет зависимость всех последующих значений {xk+1, xk+2, ...}
последовательности X от последнего фиксированного значения xk,
причем зависимость ослабевает с ростом номера значения (реально
|β| < 1). Условная корреляционная матрица имеет общий вид
10
 1 − β2
β(1 − β2 ) β2 (1 − β2 ) β3 (1 − β2 ) ...



 β(1 − β2 )
1 − β4
β(1 − β4 ) β2 (1 − β4 ) ...


B2 = σ2 β2 (1− β2 ) β(1− β4 )
1 − β6
β(1− β6 ) ... , (1.6)

 3

β (1− β2 ) β2 (1− β4 ) β(1− β6 )
1 − β8
...




...
...
...
...
...


свидетельствующий о том, что условная дисперсия c ростом номера
значения растет, приближаясь к σ2, коэффициенты корреляции
приближаются к значениям β, β2, β3, ..., как в исходной корреляционной матрице (1.1). Эти свойства условных математического ожидания и корреляционной матрицы описывают переходной процесс –
развитие марковской траектории, начинающейся в фиксированной
точке xk, в которой дисперсия равна нулю.
Полученные ранее результаты для марковских последовательностей первого порядка – частный случай общей теории случайных
процессов. В теории случайных процессов рассматриваются в том
числе условные процессы, начинающиеся в заданной точке x0. С течением времени условные процессы приближаются к безусловным
процессам. Развитие условного процесса от точки x0 описывается
условной функцией корреляции
R (τ | x0 ) = R (τ) 1 − R 2 (τ). (1.7)
и условным математическим ожиданием M | x0.
m(t) | x0 = x0 R (t)+ m(t),
где m(t) – математическое ожидание безусловного процесса. При
m(t) = 0
m(t) x0 = x0 R (t). (1.8)
Дискретизация функций (1.6) и (1.7) дает корреляционную матрицу B|x0 и математическое ожидание M[X|x0] условной последовательности X|x0. Нетрудно заметить, что переход от безусловной корреляционной матрицы марковского процесса (1.1) к условной матрице (1.5), а также формирование условного математического ожидания с помощью матрицы (1.5) соответствуют общим правилам
(1.7) и (1.8).
Пример 1.2. Зафиксировано значение xk = 2.5 марковской последовательности с нулевым средним и функцией корреляции
11
R (τ )= exp(−α τ ), α = −ln (β), β = 0.7, α = 0.3567;
интервал дискретизации ∆ = 1. Далее в этом примере обозначено
x0 = xk. Моделируются значения N = 1000 условных траекторий,
описывающихся формулами (1.7) и (1.8). Программа моделирования
N=1000
x0=2.5
del=1
T=15
a=-log(0.7) t=0:del:T;
r=exp(-a*t) n=length(t) %
%
%
%
%
число реализаций
начальная точка x0
интервал дискретизации
длительность реализаций
показатель экспоненты
% функция корреляции
% число отсчётов
for i=1:n
for j=1:n
b(i,j)=r(abs(i-j)+1)*sqrt((1-r(i)^2)*(1-r(j)^2));
end
end
B=b;
% условная корреляционная матрица
[u,v]=eig(B);
x=randn(N,n);
BX=cov(x);
[ux,vx]=eig(BX);
A=u*v^(1/2)*u‘*ux*vx^(-1/2)*ux‘; % оператор генерирования
y=A*x‘;
Y=y‘;
m=x0*r
% условное м.о.
for i=1:N
MX(i,:)=m(:);
end
z=Y+MX;
% условные марковские траектории
Z=z(1:10,:) % десять их реализаций
z1=3*(1-r.^2)+m
% верхняя граница
z2=-3*(1-r.^2)+m
% нижняя граница
plot(t,Z,’k’,t,z1,’--’,t,z2,’--’,t,m,t,zeros(1,length(t)))
Десять траекторий, развивающихся из точки x0 вокруг услов­
ного математического ожидания (жирный пунктир) показаны на
рис. 1.1. Линиями 1 и 2 показаны границы
12
9
Y
U
s
s
Рис. 1.1. Условные траектории, σ2 = 1
z1 = m(t | x0 ) + 3σ (t | x0 ), z2 = m(t | x0 ) − 3σ (t | x0 ),
соответствующие правилу “шести сигма”.
Следует отметить, что марковская последовательность первого
порядка может описываться функцией корреляции
R (τ ) = σ2exp(−α τ )cos πτ (1.9)
при интервале дискретизации ∆ = 1. Функция (1.9), показанная на
рис. 1.2 штрихами, отличается от функции (2) тем, что она знакопеременна. Матрица точности и в этом случае трехдиагональна.
3
s
s
s
T
Рис. 1.2. Функции корреляции
13
9
Y
U
s
s
s
Рис. 1.3. Условные траектории, σ2 = 1
Условная корреляционная матрица такая же, как в предыдущем
случае.
Условное математическое ожидание также знакопеременно
(рис. 1.3).
1.2. Общие соотношения для марковских последовательностей первого порядка
Матрица (7) по сути является матрицей связи наблюдаемых значений с предыдущими. Необходимый и достаточный признак марковской последовательности первого порядка – вид матрицы связи
(1.5) – элементы всех ее столбцов кроме последнего равны нулю.
Одно из фундаментальных свойств марковских последовательностей: любая подпоследовательность является марковской [1]. Следовательно, блок B11 так же, как матрица B, есть корреляционная матрица марковской последовательности первого порядка (укороченной по сравнению с исходной). Марковская последовательность описывается обратной матрицей D = B–1 (матрицей точности) вида (1.2),
−1
подпоследовательность – матрицей B11
. Матрица точности марковской последовательности первого порядка и матрицы точности под−1
последовательностей (матрицы B11
) имеют характерную форму –
они трехдиагональны. Переход от описания последовательности
корреляционной матрицей к описанию матрицей точности позволяет сформулировать необходимый признак принадлежности стацио14
нарной гауссовой последовательности классу марковских первого
порядка [4]: трехдиагональность матрицы точности.
Трехдиагональность матрицы точности есть также достаточный
признак марковости первого порядка. Если матрицы B и D разбиты
на блоки
 B11 B12 


 , D =  D11 D12  ,
B = 
2 
 D21 d22 


 B21 σ 
соответствующие условной плотности f(xn| x1, ..., xn – 1), их произведение
 B11D11 + B12D21 B11D12 + B12d22 


 = I =  I11 I12  , (1.10)
BD = 

2
2


B21D12 + σ d22 
I21 1 
 B21D11 + σ D21
I21, I12 – строка и столбец нулей; I11 – единичная (n – 1)-матрица.
−1
Матрица связи (7) C = B21 B11
,. так что
B21 = CB11. (1.11)
Условная дисперсия (5)
σ2c = σ2 | x1,..., xn−1 = σ2 −CB12,
так что
CB12 = σ2 − σ2c . (1.12)
Одно из соотношений (1.10)
B11D11 + B12D21 = I
определяет равенство
B11D11 = I − B12D21. (1.13)
Другое соотношение (1.10)
B21D11 + σ2D21 = 0
С учетом (1.11), (1.13) и (1.12) записывается
B21D11 + σ2D21 = CB11D11 + σ2D21 = C(I − B12D21 ) + σ2D21 =
= C − (σ2 − σ2c )D21 + σ2D21 = C + σ2c D21 = 0,
т. е.
C = −σ2c D21. (1.14)
15
Если матрица точности D трехдиагональна, матрица D21 =
= [0 ... 0 β],
C = −σ2c D21 = [0 ... 0 β ]
– отличие от нуля одного последнего элемента в векторе связи определяет первый порядок последовательности: свойство трехдиагональности матрицы точности является также достаточным признаком марковской последовательности первого порядка.
Свойство трехдиагональности распространяется на нестационарные марковские последовательностиc [5]. Поскольку случайная
траектория может иметь произвольное среднее, его можно положить равным нулю, тогда нестационарность сводится к переменной
дисперсии отсчетов, т. е. корреляционная матрица нестационарной
марковской последовательности приобретает вид

σ12
βσ1σ2 β2σ1σ3
...
βn−1σ1σn 



 βσ2σ1
σ22
βσ2σ3
...
βn−2σ2σn 


 , σi ≠ σ j .
2
n−3
Bí =  β2σ σ
βσ
σ
σ
...
β
σ
σ

3 1
3 2
3
3 n


...
...
...
...
...




−
n
1
2
β σ σ

...
...
βσn σn−1
σn
n 1


Соответствующая матрица точности трехдиагональна с элементами главной диагонали
{σ1−2
2
2
(1 + β2 )σ−
... (1 + β2 )σ−
2
n−1
2
2
σ−
n } / (1 − β )
и элементами следующих диагоналей
1 −1
2
dij = −bσ−
i σ j / (1 − b ),
i, j – номера строк и столбцов. Например, пусть корреляционная матрица B стационарной последовательности трансформируется в матрицу Bн последовательности с переменной дисперсией:
1 β


β 1
B=
 2 β
β
 3
β β2

 1
b2 β3 


2 
 2β
β β
 , Bí = 

 2
1 β
 3β

 3
4β
β 1 

матрица точности трехдиагональна:
16
2β
4
6β
8β2
4β3 

6β 8β2 
;
9 12β 

12β 16 
3β2
 1
−β / 2
0
0 


−β / 2 (1 + β2 ) / 4

−
β
/
6
0
1 

D=

.
2
1 − β2  0
(1 + β ) / 9 −β / 12
−β / 6


 0
0
−β / 12
1 / 16 

Соотношения (1.10) позволяют получить еще один результат. Из
равенства
B21D12 + σ2d22 = 1
следует σ2 =
1 − B21D12
.
d22
Условная дисперсия (5) значения xn равна
−1
1
B D + d22B21B11
B12
−1
σ2c = σ2 − B21B11
B12 =
− 21 12
. (1.15)
d22
d22
Из равенства
B11D12 + d22B12 = I12 = 0
следует
−1
D12 = −d22B11
B12. (1.16)
Подстановка (1.16) в (1.15) дает соотношение для условной дисперсии
σ2c =
−1
1
−d B B−1B + d22B21B11
B12
1
− 22 21 11 12
=
d22
d22
d22
(1.17)
– она обратна последнему элементу матрицы точности. В справедливости формулы (1.17) можно убедиться по примеру 1.1: значение
d22 = 1.9608, условная дисперсия σc2 = 1/d22 = 0.51.
Этот результат справедлив также для нестационарных марковских последовательностей. Так, в предыдущем примере для β = 1/2
матрица связи
C = [0 0 2 / 3], d22 = 1 / 12,
2
условная дисперсия σc2 = d–1
22 = σ – CB12 = 16 – 4 = 12.
17
2. Гауссовы марковские последовательности высших порядков
2.1. Матрица точности
Марковская последовательность k-го порядка определяется как
последовательность с условной плотностью распределения
f (x p+1,..., xn | x1,..., x p )= f (x p+1,..., xn | x p−k+1,..., x p )
– последние значения зависят от k предыдущих. Такая зависимость
задает общий вид матрицы связи (7), (1.5): в ней отличен от нуля
хотя бы один элемент каждого из k последних столбцов [4,5]:
0

0

C=
...

0

c1 p 

... 0
c2( p−k+1)
...
c2 p 
. ... ...
...
...
... 

... 0 c(n− p)( p−k+1) ... c(n−k)p 

...
0
c1( p−k+1)
...
(2.1)
Например, при наблюдении семи значений процесса третьего порядка с нулевыми средними условное среднее (7) шестого и седьмого значений равно
 c13x3 + c14 x4 + c15x5 
,
Μ  x6 , x7 | x1,..., x5  = CX = 
 c23x3 + c24 x4 + c25x5 


так как матрица связи (2.1) имеет вид
0 0 c13
C=
0 0 c23

c14
c24
c15 
.
c25 
При умножении (1.10) блочных матриц DB = I одно из четырех
произведений
D21B12 + D22B22 = I,
I – (n – p) × (n – p) – единичная матрица (p первых значений фиксированы). Отсюда
1
B22 = D−
22 (I − D21B12 ),
1
−1
D−
22 = B22 + D22 D21B12 , (2.2)
следовательно, для матрицы связи и условной корреляционной матрицы Bc имеют место соотношения
18
1
C = −D−
22 D21, (2.3)
1
Bc = D−
22 , (2.4)
так как тогда условная корреляционная матрица (5)
Bc = B22 − CB12
и матрица (2.2) равны. Таким образом, если элементы m первых
1
столбцов матрицы D21 равны нулю, а элементы матрицы D−
22 отличны от нуля, матрица (2.3) имеет вид матрицы D21: m первых
столбцов нулевые. Равенство (2.4) обобщает результат (1.17).
На основании этих результатов можно сформулировать утверждение о том, что матрица точности марковской последовательности
k-го порядка имеет 2k + 1 главных (ненулевых) диагонали [4, 5]
(остальные элементы матрицы точности равны нулю).
Необходимость 2k + 1-диагональности матрицы точности имеет
смысл утверждения о том, что матрица точности марковской последовательности k-го порядка 2k + 1-диагональна. Матрица связи последовательности k-го порядка имеет вид (2.1), а это возможно, если
ее блок D21 при k < p имеет p – k первых нулевых и k ненулевых
столбцов, т. е. для 2k + 1-диагональных матриц D. Например, условно изображенную матрицу точности (звездочками обозначены ненулевые элементы) с семью ненулевыми диагоналями можно разбить
на блоки следующим образом (k = 3, p = 5):
∗ ∗ ∗ ∗ 0

∗ ∗ ∗ ∗ ∗

∗ ∗ ∗ ∗ ∗

D =  ∗ ∗ ∗ ∗ ∗
0 ∗ ∗ ∗ ∗


0 0 ∗ ∗ ∗

0 0 0 ∗ ∗
Форма произведения (2.3)
0
0
∗
∗
∗
∗
∗
0

0

0
0 0 ∗ ∗ ∗
.
∗  , D21 = 
0 0 0 ∗ ∗



∗

∗

∗ 
∗ ∗ 0 0 ∗ ∗ ∗ 0 0 ∗ ∗ ∗
1


=

C = −D−
22 D21 = 

 

∗ ∗ 0 0 0 ∗ ∗ 0 0 ∗ ∗ ∗
определяет вид матрицы связи двух значений марковской последовательности третьего порядка с пятью предыдущими (два первых
столбца – нулевые).
19
Смысл достаточности 2k + 1-диагональности матрицы точности
в том, что если матрица точности 2k + 1-диагональна, то D – матрица точности последовательности k-го порядка. Достаточность можно доказать следующим образом.
1
Матрица связи (2.3) C = −D−
22 D21 = −B c D21, условная корреляционная матрица Bc, как любая корреляционная матрица, – теоретически неотрицательно определенная: detBc ≥ 0. Практически матрица
Bc – положительно определенная матрица: detBc > 0. Пусть элементы
первого столбца матрицы C равны нулю (фиксированы первые p из n
значений последовательности k-го порядка; rij – элементы матрицы
–Bc). Произведение –BcD21 для первого столбца в развернутом виде
r( p+1)( p+1)d( p+1)1 + r( p+1)( p+2)d( p+2)1 + ... + r( p+1)ndn1 = c11 = 0,
r( p+2)( p+1)d( p+1)1 + r( p+2)( p+2)d( p+2)1 + ... + r( p+2)ndn1 = c21 = 0,
...............................................................................
rn( p+1)d( p+1)1 + rn( p+2)d( p+2)1 + ... + rnndn1 = c(n−k)1 = 0.
(2.5)
При заданных rij (2.5) – система однородных линейных уравнений относительно dm1, m = n – p. Определитель detBc > 0 не равен
нулю, следовательно, система (2.5) имеет нулевое решение. Система
уравнений вида (2.5) существует для второго, третьего и т. д. нулевого столбца матрицы связи. Таким образом, столбец матрицы связи (2.1) равен нулю только в том случае, если равен нулю соответствующий столбец матрицы D21, т. е., если матрица точности
2k + 1-диагональна.
Пример 2.1. Матрица точности определяется корреляционной
матрицей, но не наоборот. Матрица точности не может быть произвольной. Попытка задать
 1.00
0.70 −0.30 0.10 


 0.70
1.00
0.70 −0.30


D=

−
0
.
30
0
.
70
1
.
00
0
.
70


 0.10 −0.30 0.70
1.00 

равносильна заданию матрицы
 0.2916
0.5680 −0.8019 0.7026 


 0.5680 −0.2948 0.9381 −0.8019
−1

,
b=D =

−
−
0
.
8019
0
.
9381
0
.
2948
0
.
5680


 0.7026 −0.8019 0.5680
0.2916 

20
 0.9162
0.4089

 0.4089
0.6640

−0.2476 0.1153
B = D−1 = 
−0.2917 −0.2321
 0.0147 −0.1396


0.0147
 0.1516
−0.2476 −0.2917 0.0147
0.1516 

0.11523 −0.2321 −0.1396 0.0147 

0.6777
0.2667 −0.2321 −0.2917
0.2667
0.6777
0.1153 −0.2476
−0.2321 0.1153
0.6640
0.4089 

−0.2917 −0.2476 0.4089
0.9162 
которая не является корреляционной.
Пример 2.2. В общем случае корректно заданная матрица точности соответствует нестационарной последовательности. Так, добавление пары диагоналей к матрице D последовательности первого
порядка из примера 1.1 задает матрицу точности второго порядка
 1.9608 −1.3725 0.9500

0
0
0


−1.3725 2.9216 −1.3725 0.9500

0
0


 0.9500 −1.3725 2.9216 −1.3725 0.9500

0

,
D=

0
0
9500
1
3725
2
9216
1
3725
0
9500
.
.
.
.
.
−
−



0
0
0.9500 −1.3725 2.9216 −1.3725



0
0
0
0.9500 −1.3725 1.9608 

соответствующую корреляционной матрице
 0.9162
0.4089 −0.2476 −0.2917 0.0147
0.1516 


 0.4089
0.6640
0.1153 −0.2321 −0.1396 0.0147 



0.1153
0.6777
0.2667 −0.2321 −0.2917
−1 −0.2476
B=D =
.
0.6777
0.1153 −0.2476
−0.2917 −0.2321 0.2667
 0.0147 −0.1396 −0.2321 0.1153
0.6640
0.4089 



0.0147 −0.2917 −0.2476 0.4089
0.9162 
 0.1516
Свойство 2k + 1-диагональности матрицы точности распространяется и на случай, когда последними ненулевыми являются k – е
диагонали, а среди предыдущих диагоналей встречаются нулевые
(главная диагональ имеет номер k = 0). Последовательности с такими матрицами точности по-прежнему имеют порядок k.
Пример 2.3. Отсчеты процесса с функцией корреляции R(τ) =
= e–|τ|cosπτ, взятые с интервалом, образуют марковскую последовательность второго порядка. Действительно, n = 6 отсчетов описываются корреляционной матрицей
21
 1.0000

0
0
0.11353
0
−0.3679



0
1.0000
0
0
0.1353 
−0.3679


−0.3679

0
1.0000
0
0
−0.3679


B=

0
0
3679
0
1
0000
0
0
3679
.
.
.
−
−


 0.1353

0
0
1.0000
0
−0.3679




0
0
1353
0
0
3679
0
1
0000
.
−
.
.


и пяти – диагональной матрицей точности
1.1565
0
0.4255
0
0
0 


 0
1.1565
0
0.4255
0
0 


0.4255
0
1.3130
0
0.4255
0 

D=
.
0.4255
0
1.3130
0
0.4255
 0
 0
0
0.4255
0
1.1565
0 



0
0
0.4255
0
1.1565
 0
При фиксировании первых трех значений последовательности
матрица связи (2.3) и условная корреляционная матрица (2.4)
0 −0.3679

 0.8647
0
0
0.3181




0
0
0.8647
0  .
C = 0
−0.3679 , Bc = 
0 0.1353

−0.3181
0
0
0.9817



Условные математические ожидания равны
m4 x1,..., x4 = m4 − 0.3679x2,
m5 x1,..., x4 = m5 − 0.3679x3,
m6 x1,..., x4 = m6 + 0.1353x3.
Этот результат предсказуем корреляционными свойствами последовательности: ее значения зависят через одно.
Анализ четырех составляющих блочного произведения (1.10)
DB = I приводит также к следующим соотношениям для матрицы
связи
−1
−1 −1
C = −B22D21D11
B11 = −B22D21 (I − B12D21 )
(2.5)
и вектору связи последнего значения xn при фиксированных предыдущих [4]
22
C =−
1 
D21
−1
0 ... 0 d p−k+1 ... d p  , B22 = dnn
. (2.6)
=−
dnn
dnn 
Таким образом, матрица связи и условная корреляционная матрица могут рассчитываться с использованием только исходной
корреляционной матрицы B по формулам (5) и (7), только матрицы
точности D по формулам (2.3), (2.4) и (2.6) или по смешанным формулам (2.5).
2.2. Аппроксимация немарковского гауссова процесса марковской последовательностью
Марковское приближение отсчетов X немарковского гауссова
процесса с корреляционной матрицей BX реализуется приведением
его матрицы точности DX к многодиагональному виду D заменой
нулями всех элементов за исключением элементов 2k + 1 главных
диагоналей. Матрица D при условии положительной определенности матрицы B = D–1 есть матрица точности марковской последовательности k-го порядка (в общем случае нестационарной). Если обнуляются элементы dij ≈ 0, можно ожидать, что погрешность такого
представления невелика, т. е. B ≈ BX. На базе корреляционной матрицы B можно построить программу – генератор марковской последовательности XM, близкой к последовательности X по корреляционным свойствам. Условные матрицы (16), (17) позволяют генерировать траектории XM точка за точкой, парами, триадами и т. д.
точек [4,5].
Пример 2.4. Корреляционные матрицы и матрицы точности отсчетов X1 и X2 процессов с функциями корреляции, показанными
на рис. 2.1, а, 2.1, б (кривые 1),
−τ
R1 (τ ) = 0.5e
1.0000

 0.4872

0.2516
B1 = 
 0.1365
0.0768


0.0444
−2 τ
+ 0.5e
0.4872
1.0000
0.4872
0.2516
0.1365
0.0768
−τ
, R2 (τ ) = 2e
0.2516
0.4872
1.0000
0.4872
0.2516
0.1365
−2 τ
−e
0.1365
0.2516
0.4872
1.0000
0.4872
0.2516
, ∆ = 0.5, n = 5
0.0768
0.1365
0.2516
0.4872
1.0000
0.4872
0.0444

0.0768

0.1365
,
0.2516
0.4872

1.0000 
23
3
¸
¹
3
U
U
Рис. 2.1. Аппроксимация марковскими последовательностями
1.0000

0.8452

0.6004
B2 = 
0.3965
0.2524


0.1574
0.8452
1.0000
0.8452
0.6004
0.3965
0.2524
0.6004
0.8452
1.0000
0.8452
0.6004
0.3965
0.3965
0.6004
0.8452
1.0000
0.8452
0.6004
0.2524
0.3965
0.6004
0.8452
1.0000
0.8452
0.1574

0.2524

0.3965
,
0.6004
0.8452

1.0000 
 1.3119

−0.6269

−0.0186
D1 = 
−0.0092
−0.0046


−0.0030
−0.6269
1.6114
−0.6180
−0.0142
−0.0071
−0.0046
−0.0186
−0.6180
1.6117
−0.6179
−0.0142
−0.0092
−0.0092
−0.0142
−0.6179
1.6117
−0.6180
−0.0186
−0.0046
−0.0071
−0.0142
−0.6180
1.6114
−0.6269
 4.2059

−5.1482

 2.2235
D2 = 
−0.5542
 0.1350


−0.0258
−5.1482
10.5074
−7.8691
2.8985
−0.7058
0.1350
2.2235
−7.8691
11.6785
−8.1443
2.8985
−0.5542
−0.5542
2.8985
−8.1443
11.6785
−7.8691
2.2235
0.1350 −0.0258

−0.7058 0.1350 

2.8985 −0.5542
.
−7.8691 2.2235 
10.5074 −5.1482

−5.1482 4.2059 
24
−0.0030

−0.0046

−0.0092
,
−0.0186
−0.6269

1.3119 
Элементы d16, d15, d14, d13; d26, d25, d24; d36, d35; d45 матрицы D1
по абсолютной величине значительно меньше предыдущих, следовательно, можно обнулить их:
 1.3119 −0.6269

0
0
0
0


−0.6269 1.6114 −0.6180

0
0
0




0
0
0
−0.6180 1.6117 −0.6179


D1′ = 
.
0
0
0
6179
1
6117
0
6180
0
.
.
.
−
−



0
0
0
−0.6180 1.6114 −0.6269



0
0
0
0
−0.6269 1.3119 

При этом исходная корреляционная матрица воспроизводится
с относительно небольшими погрешностями (рис. 2.1, а кривая 2):
0.9855

 0.4671

0.2183
B1′ = 
 0.1021
 0.0481


0.0230
0.4671
0.9775
0.4567
0.2137
0.1007
0.0481
0.2183
0.4567
0.9695
0.4536
0.2137
0.1021
0.1021
0.2137
0.4536
0.9695
0.4567
0.2183
0.0481
0.1007
0.2137
0.4567
0.9775
0.4671
0.0230

0.0481

0.1021
.
0.2183
0.4671

1.0000 
Немарковская последовательность с функцией корреляции R1(τ)
удовлетворительно аппроксимируется марковской первого порядка.
В матрице D2 можно обнулить только элемент d16, аппроксимируя немарковскую последовательность с функцией корреляции
R2(τ) марковской четвертого порядка (рис. 2.1, б кривая 2). Попытка
уменьшить порядок аппроксимации до третьего, обнуляя элементы
d16, d15; d26, приводит к недопустимым погрешностям (рис. 2.1, б
кривая 3).
Необходимо заметить, что реализация конечной протяженности
n немарковской последовательности фактически есть реализация
марковской последовательности n – 1-го порядка, так как последнее
значение зависит от n – 1 предыдущих. На этом основании можно
было бы не заниматься аппроксимацией, однако, ее цель – минимизация порядка, т. е. максимальное упрощение модели.
25
3. Генерирование марковских траекторий
Рассматриваются генераторы двух типов: генераторы множества
реализаций и генераторы полубесконечных траекторий. Множество
N реализаций гауссовых траекторий протяженностью n значений
предназначено для моделирования задач математической статистики. Оно описывается корреляционной матрицей Bn и вектором средних значений. Полубесконечная траектория генерируется в единственном экземпляре. Она имитирует реализацию эргодического
гауссова процесса произвольной длительности.
Понятия произвольной и конечной протяженности траекторий
условны. Например, в системе MATLAB без существенных затрат
времени можно сформировать массив N = 1000 траекторий протяженностью n = 1000 значений каждая. Иными словами, каждая
траектория имитирует достаточно протяженную последовательность, которая может относиться к полубесконечной (в зависимости
от конкретной задачи).
Марковские нормальные последовательности, разумеется, можно генерировать общими методами, разработанными для нормальных процессов [6, 7]. Однако некоторые особенности марковских последовательностей позволяют расширить арсенал методов генерирования.
3.1. Генерирование множества реализаций
Генераторы нормальных последовательностей обычно пытаются
построить как линейные системы, преобразующие векторы белого
шума X ∈ Ν(0, I) в векторы Y ∈ Ν(0, B) с заданной корреляционной
матрицей B ≠ I. Такие системы называют окрашива-ющими. Поскольку программ – генераторов отсчетов строго белого шума не существует, вместо белого шума приходится использовать окрашенный – результат работы функций типа RANDN [2]. Следовательно,
неизбежны погрешности генерирования.
Пример 3.1. Фрагмент MATLAB – программы
N=1000; n=6
X=randn(N,n);
BX=cov(X)
формирует (1000 × 6) – массив X с корреляционной матрицей
26
 0.8902 −0.0528 0.0462
0.0078
0.0171 −0.0019


−0.0528 1.0635
0.0025
0.0408 −0.0166 −0.0012


 0.0462
0.0462
1.0502 −0.0150 −0.0588 0.0171 

BX = 
,
0.0408 −0.0150 0.9826
0.0070
0.0374 
 0.0078
 0.0171 −0.0166 −0.0588 0.0070
1.0532 −0.0151



0.0374 −0.0151 1.0627 
−0.0019 −0.0112 0.0171
отличающейся от единичной, т. е. формируется массив окрашенного шума.
Преобразование вектора X в вектор Y задается общим векторным
соотношением
Y = AX, (3.1)
в котором квадратная матрица A называется оператором системы.
Если X – вектор окрашенного шума, соотношение (3.1) описывает
процедуру перекрашивания – преобразования вектора X ∈ Ν (0, B X )
в вектор Y ∈ Ν (0, BY ). Задача синтеза перекрашивающей системы –
поиска оператора A по заданным матрицам BX и BY – решается с помощью их сингулярного разложения [6, 7]: оператор
1/2
1/2 T
" #1Y/2#
6Y ,1Y/26T
Y 6X , X 6X , X
(3.2)
где UX, UY – собственные векторы, ΛX, ΛY – собственные значения
матриц BX и BY. Оператор (3.2) описывает генерирование как стационарных, так и нестационарных последовательностей.
Пример 3.2. Марковская последовательность четвертого порядка, аппроксимирующая немарковский процесс из примера 2.4, нестационарна, ее корреляционная матрица
0.9926

0.8359

 0.5883
B = 
0.3798
0.2298


 0.1312
0.8359
0.9896
0.8331
0.5852
0.3768
0.2298
0.5883
0.8331
0.9882
0.8323
0.5852
0.3798
0.3798
0.5852
0.8323
0.9882
0.8331
0.5883
0.2298
0.3768
0.5852
0.8331
0.9896
0.8359
0.1312

0.2298

0.3798
.
0.5883
0.8359

0.9926
Собственные векторы
27
 0.1530
0.3070
0.4463
0.5424 −0.5267

−0.4016 −0.5675 −0.3812 0.0643 −0.4397

 0.5615
0.2894 −0.3943 −0.4491 −0.1711
U = 
0.3943 −0.4491 0.1711
−0.5615 0.2894
 0.4016 −0.5675 0.3812
0.0643
0.4397


0
1530
.
−
0
.
3070
−
0
.
4463
0
.
5424
0.5267

0.3341

0.4170

0.4632
,
0.4632
0.4170

0.3341
собственные значения
¨0.0354
0
0
0
0
0 ·
©
¸
© 0
0.0665
0
0
0
0 ¸
©
¸
© 0
¸
0
0
1596
0
0
0
.
¸;
, ©©
0
0
0.4484
0
0 ¸¸
© 0
© 0
0
0
0
1.4352
0 ¸¸
©
©
¸
0
0
0
0
3.7957¸¹
©ª 0
корреляционная матрица B воспроизводится как B’ = UΛUT с погрешностью
−0.0777 −0.0111

−0.0111 0.0222


0
−14 −0.0222
×
∆B = 10
0
0333
0
0333
.
.

 0.0278
0.0222


 0.0139 −0.0111
−0.0222
0.0111
−0.0111
0.0111
−0.0111
−0.0111
0.0333 0.0278
0.0028 

0.0333 0.0222 −0.0111

0.0111 0.0111 −0.0111
,
0.0333 0.0333
0.0333 
0
−0.0888 −0.0111

0.0222 −0.0111 0.1110 
соответствующей машинной точности.
Пример 3.3. Стационарный гауссов процесс с функцией кор­
реляции
π
R (τ ) = exp(− τ )cos τ
2
на интервале T = 4 аппроксимируется марковской последовательностью с интервалом дискретизации Δ = 1/2. Подобрать минимальный
порядок аппроксимации.
В матрице точности размерностью 9 × 9 последовательно обнуляются угловые элементы
28
d19 ;
d18 , d19 ; d29 ;
.........................
d13,..., d19 ; ... ;d69
3
для аппроксимации после- довательностями восьмого, седьмого, ..., второго порядка. Функции корреля- ции показаны на рис. 3.1,
U
кривые: 1 – исходная; 2, 3,
4 – аппроксимация после- s
довательностями четвертоРис. 3.1. Функции корреляции
го, третьего и второго порядка. По-видимому, следует остановиться на аппроксимации третьего порядка. Аппроксимирующая последовательность задается
в строгом смысле нестационарной – дисперсии значений
σˆ 2 =[0.9914;0.9879;0.9886;0.9847;0.9781;0.9847;0.9886;0.9879;0.9914];
практически же такую последовательность можно считать стационарной.
Несколько реализаций Y, изображенных функцией PLOT показаны на рис. 3.2. Оценки дисперсии значений, полученные по массиву N = 1000 реализаций,
σˆ 2 =[0.9914;0.9879;0.9886;0.9847;0.9781;0.9847;0.9886;0.9879;0.9914]
равны заданным.
:
U
s
s
s
Рис. 3.2. Марковские траектории третьего порядка
29
3.2. Генерирование полубесконечных последовательностей
Полубесконечная марковская последовательность произвольного порядка как нами отмечалась ранее, может генерироваться число за числом, парами чисел, триадами, квартетами и т. д. на базе
матрицы связи C (2.3) и условной корреляционной матрицы Bc (2.4).
Генерируются последовательности с нулевыми средними и единичной дисперсией.
3.2.1. Генерирование число за числом
Матрица связи очередного генерируемого числа в последовательности k-го порядка имеет структуру
C = [0 ... 0 ∗ ... ∗]
– последние k значения отличны от нуля, предыдущие – равны
нулю. Следовательно, матрица связи записывается без нулевых элементов:
C = [c1 ... ck ].
Дисперсия генерируемого числа в установившемся режиме (1.17)
σ2c =
1
,
d22
d22 – последний диагональный элемент матрицы точности. Установившийся режим наступает с генерирования k + 1-го числа, первые
k чисел формируются последовательно в соответствии с формулами
(5) – (7).
Таким образом, процедура генерирования последовательности X
включает:
1) генерирование первого значения x1 ∈ Ν (0,1);
2) по корреляционной матрице двух значений
 b11
B=
 b21

b12 

b22 
находятся коэффициент (7) C = b21/b11 и с. к. о. (5) σ2 = b22 − cb12 ,
что позволяет получить второе значение x2 ∈ Ν (cx1, σ2 );
30
3) из корреляционной матрицы трех значений
 b11 b12 b13 


B =  b21 b22 b23 
b

 31 b32 b33 
−1
находятся C = [c1 c2] = B21 B11
и σ3 = b33 −CB12 , генерируется
x3 ∈ Ν (c1x1 + c2x2, σ3 );
4) этап 3 повторяется с наращиваемыми корреляционными матрицами, и процесс заканчивается с матрицей
 b11 ... b1k 


B =  ... ... ... 
b

 k1 ... bkk 
−1
получением вектора связи (2.3) C = [x1 ... xk – 1] = B21 B11
= –D21/dkk
k−1

и с. к. о. σk, генерируется последнее xk ∈ Ν  ∑ ci xi , σk  из группы

 i=1
чисел, от которой зависит следующее число;
5) начинается установившийся режим: следующие значения
xk+1, xk+2, ..., xk+j, ... генерируются с постоянными матрицей связи
k
C и дисперсией σ2c (2.6); средние значения Μ [xk+ j ] = ∑ ci xi+ j . Генеi=1
рируется последовательность с нулевыми средними, так как процедура линейная.
Завершающий этап – оценка корреляционных свойств полученной последовательности. Классические оценки корреляционных
моментов центрированных чисел
rˆk =
N−k+1
1
∑ xi xi+k−1, k = 1, 2,..., N − k +1 i
(3.3)
строго говоря, применимы к эргодическим процессам, а генерируемая марковская последовательность четвертого порядка нестационарна. Тем не менее, нестационарность выражена слабо, и оценку
(3.3) можно рассматривать, как приближенную.
Пример 3.4. В примере 2.4 марковской последовательностью
с четвертого порядка с интервалом дискретизации ∆ = 0.5 аппрокси31
мируется процесс с функцией корреляции R2(τ) = 2e–|τ| – e–2|τ|. Корреляционная матрица аппроксимирующей последовательности
0.9927

0.8360

0.5884
B = 
0.3800
0.2303


0.1317
0.8360
0.9897
0.8333
0.5853
0.3771
0.2303
0.5884
0.8333
0.9883
0.8324
0.5853
0.3800
0.3800
0.5853
0.8324
0.9883
0.8333
0.5884
0.2303
0.3771
0.5853
0.8333
0.9897
0.8360
0.1317

0.2303

0.3800
.
0.5884
0.8360

0.9927
Программа генерирования
N=10000 % число реализаций
del=1/2 % интервал дискретизации
T=2.5 % длительность реализаций
a=2
b=4/2*pi
t=0:del:T;
r=2*exp(-t)-exp(-2*t) % функция корреляции
n=length(t) % число отсчётов
for i=1:n
for j=1:n
b(i,j)=r(abs(i-j)+1);
end
end
BB=b % исходная корреляционная матрица
D=inv(BB) % матрица точности
D(1,6)=0
D(6,1)=0
B=inv(D % аппроксимирующая корреляционная матрица
% корреляционные матрицы первых значений
B3=[B(1,1)
B(2,1)
B4=[B(1,1)
B(2,1)
B(3,1)
B5=[B(1,1)
B(2,1)
32
B(1,2)
B(2,2)]
B(1,2)
B(2,2)
B(3,2)
B(1,2)
B(2,2)
B(1,3)
B(2,3)
B(3,3)]
B(1,3) B(1,4)
B(2,3) B(2,4)
B(3,1)
B(4,1)
B(3,2)
B(4,2)
B(3,3) B(3,4)
B(4,3) B(4,4)]
x(1)=randn*sqrt(B(1,1))
% первое значение последовательности
c=B(2,1)/B(1,1)
si=sqrt(B(2,2)-c*B(1,2))
x(2)=c*x(1)+randn*si % второе значение последовательности
C=[B(3,1) B(3,2)]*inv(B3)
si=sqrt(B(3,3)-C*[B(1,3);B(2,3)])
x(3)=C(1)*x(1)+C(2)*x(2)+randn*si % третье значение последовательности
C=[B(4,1) B(4,2) B(4,3)]*inv(B4)
si=sqrt(B(4,4)-C*[B(1,4);B(2,4);B(3,4)])
x(4)=C(1)*x(1)+C(2)*x(2)+C(3)*x(3)+randn*si
% четвёртое значение
%--------------------------------------------------------- % установившийся режим генерирования
C=[B(5,1) B(5,2) B(5,3) B(5,4)]*inv(B5)
% матрица связи
s=sqrt(B(5,5)-C*[B(1,5);B(2,5);B(3,5);B(4,5)])% с.к.о.
X=[x(1);x(2);x(3);x(4)]
for i=5:N
x(i)=C*X+s*randn; % следующее значение
for j=1:4
X(j)=x(i-4+j); % сдвиг предыдущих четырёх
end
end
plot(x) % марковская последовательность, рис. 3.3
pause
for i=1:n
s=0;
for j=1:N-i+1
s=s+x(j)*x(j+i-1);
end
rr(i)=s/(N-i+1); % оценка функции корреляции последовательности
end
plot(t,b(1,:),’r’,t,rr)
33
9
U
s
s
s
Рис. 3.3. Марковская траектория четвертого порядка, N = 10000
3
/
T
Рис. 3.4. Корреляционные моменты
Марковская последовательность показана на рис. 3.3.
На рис. 3.4 кривая 1 показаны девять первых заданных корреляционных моментов марковской последовательности, на рис. 3.4 кривая 2 – соответствующие оценки (3.3).
3.2.2. Генерирование группами чисел
Генерирование марковская последовательности k-го порядка
группами по m чисел отличается от генерирования число за числом
в установившемся режиме. Если первые k чисел получены (вектор
34
Xk), то группа из m следующих чисел (вектор Xm) формируется в соответствии с m- мерной условной нормальной плотностью
−m/2
f (Xm | Xk )= (2π )
 1
 2


1

(det Bc )−1/2 exp− (Xm − Μm )T B−
c (Xm − Μ m ),
где Bc – корреляционная матрица (2.4); вектор средних Mm = CXk;
C – матрица связи (2.3). Должна использоваться матрица точности
D размерностью k + m.
Пример 3.5. Установившийся режим генерирования парами чисел в предыдущем примере базируется на матрице точности, разбитой на блоки:
  4.2059

 −5.1482

  2.2235

D =  
 −0.5542

 0.1350
 
  0

−5.1482 2.2235 −0.5542

10.5074 −7.8691 2.8985 

−7.8691 11.6785 −8.1443
2.8985 −8.1443 11.6785 
−0.7058 2.8985 −7.8691

0.1350 −0.5542 2.2235 
 0.1350

0


−0.7058 0.1350  


 2.8985 −0.5542 

.


−7.8691 2.2235  

 10.5074 −5.1482 


−5.1482 4.2059  


Матрица связи
−0.0321 0.1285 −0.5279 1.2239
1

,
C = −D−
22 D21 = 

−0.0393 0.1253 −0.5144 0.9695
корреляционная матрица
0.2378 0.2910
1

.
Bc = D−
22 = 

0.2910 0.5940
Сингулярное разложение позволяет записать оператор (3.2) системы, окрашива-ющей два значения белого шума в соответствии
с матрицей Bc:
0.4159 0.2545
. A = B1c/2 = Uc Λ1c/2UT
c = 

0.2545 0.7275
(3.4)
Таким образом, к сформированным первым четырем значениям
{x1, x2, x3, x4} добавляются два следующих: x5 и x6 с дисперсиями,
задаваемыми матрицей Bc, и средними
m5 = −0.0321x1 + 0.1285x2 − 0.5279x3 + 1.2239x4 ,
35
m6 = −0.0393x1 + 0.1253x2 − 0.5144x3 + 0.9695x4 .
Далее процесс повторяется с множеством {x3, x4, x5, x6}.
Следует отметить, что оператор (3.4) записан в усеченном вариан1/2
те по сравнению с (3.2): отсутствует декоррелирующий блок B−
X .
Поскольку функция RANDN генерирует не строго белый шум, такое упрощение оператора может привести к дополнительным погрешностям.
Фрагмент программы (установившийся режим):
X=[x(1);x(2);x(3);x(4)];
% блоки матрицы точности
D11=[D(1,1) D(1,2) D(1,3)
D(2,1) D(2,2) D(2,3)
D(3,1) D(3,2) D(3,3)
D(4,1) D(4,2) D(4,3)
D21=[D(5,1) D(5,2) D(5,3)
D(6,1) D(6,2) D(6,3)
D22=[D(5,5) D(5,6)
D(6,5) D(6,6)];
C=-inv(D22)*D21
BC=inv(D22) [u,v]=eig(BC)
A=u*v^(1/2)*u’ D(1,4)
D(2,4)
D(3,4)
D(4,4)];
D(5,4)
D(6,4)];
% матрица связи
% корреляционная матрица
% оператор окрашивания
for i=5:2:N
xx=randn(2,1);
y=A*xx+C*X;
x(i)=y(1);
x(i+1)=y(2);
X=[x(i-2);x(i-1);x(i);x(i+1)];
end
% значения белого шума
% два следующих значения
% сдвиг значений
Удовлетворительное качество генерирования парами чисел показывает рис. 3.5.
Генерирование триадами, квартетами и т. д. чисел реализуется
по той же схеме с выделением блока D22 соответствующей размерности: 3 × 3,4 × 4 и т. д. и сдвигом формируемой последовательности
на три, четыре и т. д. числа.
36
3
/
T
Рис. 3.5. Корреляционные моменты
Указанные способы генерирования не имеют видимых взаимных
преимуществ. Можно говорить лишь о разнообразии арсенала генерирования марковских траекторий.
37
4. Примеры применения марковских моделей при расчете информационных систем
Марковские модели широко применяются для описания состояний информационных систем и сигналов [1]. При этом решаются
важные практические задачи, называющиеся задачами о достижении границ [1, 8, 9]. К ним относятся, например, задачи о срыве автоматического слежения в дальнометрии, фазовой автоподстройке
частоты в приемниках (ФАПЧ) и т. п.
Описание марковскими моделями помехи (шума) может оказаться полезным при решении задачи пересечения случайными процессами неслучайных уровней, имеющей множество приложений в информационных системах.
Одно из приложений – измерение интервалов между импульсными сигналами, оцениваемых разностью моментов пересечения заданных уровней фронтами сигналов. Решение этой задачи базируется на плотности распределения времени пересечения tu уровня
u(t) (фронта сигнала) случайным процессом x(t), описывающим шумовую составляющую.
Простую инженерную методику оценки плотности f(tu) можно
получить в пространстве дискретного времени переходом от процессов к последовательностям X, U с использованием марковских моделей и геометрического распределения.
Пусть последовательности X, U формируются на промежутке [0,
T] дискретизацией процессов x(t), u(t) с интервалом ∆. Если начальные значения x0 < U0, первое пересечение возможно снизу вверх.
Первое пересечение может произойти в k-м интервале дискретизации с вероятностью
k−1





f (tk ) = p 
(xk > uk )∩ ∩ (xi < ui ), k = 1, 2,... 



i=0


(4.1)
Если положить интервал ∆ = 1, а время первого пересечения tu
фиксировать как номер k первого выполнения неравенства xi < ui,
i = 0,1, ..., k – 1, значения f(tk) задают плотность распределения дискретного времени первого пересечения.
Если последовательность X некоррелирована (дискретный белый шум), вероятности (1) рассчитываются как произведения
k−1
f (tk ) = pk ∏ (1 − pi ), pi = p(xi > ui ). i=0
38
(4.2)
В случае коррелированной последовательности непосредственный расчет вероятностей (4.1) проблематичен даже для гауссовых
последовательностей.
4.1. Пересечение постоянного уровня
Время первого пересечения постоянного уровня U = u дискретным белым шумом распределено по геометрическому закону [3]:
k
f (tk ) = (1 − p0 ) p0 , p0 = p {x > u}= 1 − Φ (u / σ), (4.3)
соответствующему произведению (4.2) при pi = p0; Ф(λ) – интеграл
вероятности. Плотность (4.3) обобщается на слабокоррелированные
марковские последовательности конечного порядка n, задаваемые
условием τ0 << T; τ0 – интервал корреляции последовательности X.
Замена немарковской последовательности марковской позволяет
при расчете вероятностей ограничить число учитываемых составляющих. Рассчитываются вероятности пересечения уровня u первыми n значениями:
p0 = p {x > u}= 1 − Φ (u / σ),
u
pk =
∫
−∞
u ∞
... ∫
∫ f (x0,..., xk )dx0...dxk, k = 1,..., n −1.
(4.4)
−∞ u
При k ≥ n нормированные значения геометрической плотности [10]
n−1
k−n
fˆ(tk ) ≈ pk (1 − P ) = (1 − P ) pn (1 − pn )
, P = ∑ pi , (4.5)
i=0
описывают время первого пересечения уровня u.
Пример 4.1. Процесс x (t)∈ Ν (0, R (τ )) на отрезке [0, 30] дискретизируется с интервалом Δ = 0.5; функция корреляции (рис. 4.1, а кривая 1)


α
R (τ ) = exp(−ατ )cos βτ + sin βτ, α = 1 / 4, β = π / 2.


β

Аппроксимация процесса x(t) марковской последовательностью
n-го порядка базируется на приведении матрицы точности к 2n +
+ 1-диагональному виду. Функция корреляции аппроксимирующей
последовательности X порядка n = 5 показана на рис. 4.1, а (марков39
¸
3
¹
I
V
U
s
º
IQ
U
»
IQ
V
V
U
U
Рис. 4.1. Гистограммы и вероятности
ский процесс менее инерционен). На рис. 4.1, б – гистограммы времени первого пересечения уровня u = 1 исходным процессом
(N = 5000 траекторий) и марковской последовательностью X (пунктир). Вероятности (4.4) получены статистическим моделированием
пересечений последовательностью X; для уровня u = 1 в одном из
экспериментов они приняли значения p0 = 0.1618, p1 = 0.0702,
p2 = 0.0786, p3 = 0.0730, p4 = 0.0700, p5 = 0.0588, так что в расчетах
вероятностей (4.5) положены pn = 0.0588, 1 – P = 0.5464. На рис. 4.1,
в, г показаны гистограммы времени первого пересечения уровней
u = 1 и u = 2 исходным процессом и оценки вероятностей (4.5), изображенные пунктирными линиями.
4.2. Пересечение переменного уровня
Формула (4.1) по сути обобщает геометрическое распределение:
если вероятности события pi в независимых экспериментах различны, вероятность первого его наступления в k-м эксперименте равна
k−1
Pk = pk ∏ (1 − pi ), k = 1, 2,... (4.6)
i=1
Вероятность первого пересечения уровня U сверху вниз гауссовой последовательностью X с нулевым средним в i-м интервале
40
pi =
1
2π
ui
 x2 
u 
exp− 2 dx = Φ  i .
σ
 2σ 
−∞
∫
Если пересекающий процесс аппроксимировать марковской последовательностью n-го порядка, обобщенное геометрическое распределение (4.6) можно упростить [10]:
k−1
Pk = pk ∏ (1 − pi ), k ≤ n + 1,
i=1
Pk = pk
k−1
∏
i=k−n
(1− pi ), k ≥ n + 2. (4.7)
При этом первые n + 1 вероятностей рассчитываются по формуле (4.6), следующие вероятности – по формуле (4.7), сохраняющей произведение n сомножителей, что сокращает объем вычи­
слений.
Пример 4.2. Уровень u(t) = 6Ф((t-4)/2)-3, имитирующий передний
фронт импульсного сигнала, пересекаемый сверху вниз траекторией процесса x (t)∈ Ν (0, R (τ )), R(τ) = exp(–ατ)cosβτ, α = 1/2, β = π/2
показан на рис. 4.2, а. Процесс x(t) дискретизируется с интервалом
∆ = 1/2 и аппроксимируется марковской последовательностью пятого порядка. Вероятности pi из формул (4.6) и (4.7) равны pi = Ф(ui).
¸
¹
VY
IQ
U
s
s
U
Рис. 4.2. Пересечение уровня, имитирующего фронт сигнала
41
Гистограмма времени первого пересечения уровня U последовательностью X показана на рис. 4.2, б (штриховая диаграмма). Там же
приведены визуально неразличимые вероятности (4.6) и (4.7), изображенные функцией PLOT в виде непрерывной кривой.
4.3. Плотность распределения времени пересечения заданного уровня марковским процессом
Гауссов марковский процесс (первого порядка), начинающийся
в точке x0 (x(t = 0) = x0) имеет среднее
m(t) = m(t | x0 )= x0 exp(−αt)
и дисперсию
σ2 (t) = σ2 (t | x0 )= σ2 (1 − exp(−2αt)).
Пусть пересекающий уровень постоянен и равен u. Если x0 > u,
а в момент t > 0 значение процесса x(t) < u, время первого пересечения (сверху вниз) t0 < t. Вероятность такого события
p(x (t)< u)=
u
∫
−∞
 (x − m(t))2 
 u − m(t)
1

.
exp −
dx = Φ 
2


 σ (t) 
2πσ (t)
2σ (t) 


В соответствии с определением эта вероятность есть условная
функция распределения времени пересечения уровня u, возможно
не единственного:
 u − m(t0 )
. (4.8)
F (t0 | x0 )= Φ 
 σ (t0 ) 
Условная плотность распределения
 (u − m(t ))


1  u − m(t0 )′
0 
−

f (t0 | x0 )= F ′(t0 | x0 )=
. (4.9)
 exp 
2


2π  σ(t0 ) 
σ
2
t
(
)


0




Функция (4.8) обобщается [11] на случай переменного уровня u(t)
его аппроксимацией ступенчатой функцией (рис. 4.3)
v(t) = ui = u(ti−1 ),
где ti – 1 ≤ t < ti, i = 1, ..., k, индекс i = k соответствует общей длительности T.
42
VY
Y
VJ
$
Js J
5
J
Рис. 4.3. Аппроксимация уровня ступенчатой функцией
Первое пересечение может произойти в одном из k интервалов ∆
с вероятностью
 u − mi 
 u − mi−1 
 − Φ  i
,
pi = Φ  i
 σi−1 
 σi 
mi = m(ti), σi = σ(ti). Эти события несовместны, поэтому вероятность
пересечения равна
k
p = p {t0 ∈ (0, T )}= ∑ pi . (4.10)
i=1
При малом a с использованием разложения интеграла вероятности в степенной ряд [12]
3
5

1  z3 z5
Φ(z)−Φ(z+a)=
z− + −...−(z+a)+ (z+a) −(z+a) +...=

2π  6 40
6
40

 z2 
a 
z2
z4
z6
a
1−
=−
+
−
+...+δ(a)=−
exp− +δ(a). (4.11)

2π  2.1! 4.2! 8.3!
2π
 2 
На левой границе ступеньки ui


 u − mi−1 
 u − mi

 u − mi + mi′ ∆ 
 ≈ Φ  i
Φ  i
+ bi ,  = Φ  i
 σi − σi′ ∆ 

 σi−1 

 σ1



m ′σ −(mi − ui )σi′
bi = i i 2
∆.
σi − σi σi′ ∆
(4.12)
Подстановка (4.11) и (4.12) в (4.10) позволяет предельным переходом Δ → 0 получить условную функцию распределения времени пересечения уровня u(t):
43
F (t0 | x0 ) =lim∆→0 P =−
=−
 u −m 2 
m ′σ +(ui −mi )σi
 ( i
ν
i) 
lim∆→0 ∑ i i
exp
−
∆ =
2
2
2π
2σi 
σi

i
2
t0

 (u(t)−m(t)) 
m ′(t)σ(t)+(u(t)−m(t))σ ′(t)
ν
−
exp

dt,

2π ∫
2σ2(t) 
σ2(t)

0
коэффициент ν = 1/F(∞|x0) – нормирующий. Условная плотность
распределения
f (t0 | x0 )= F ′(t0 | x0 )=
=−
2

 (u(t)− m(t)) 
ν  m ′(t) (u(t)− m(t))σ ′(t)
dt. (4.13)
+
−
exp





2π  σ (t)
σ2 (t)
2σ2 (t) 


В установившемся режиме σ2(t) = σ2, σ’(t) = 0. Тогда плотность
(4.13) становится безусловной и записывается
f (t0 ) = ν
2

 (u(t0 )− m(t0 )) 
. exp−
2


2πσ
2
σ


m ′(t0 )
(4.14)
Абсолютное значение производной m’(t0) в (4.14) применено для
унификации решения: оно применимо для пересечения сверху вниз
и снизу вверх. Далее используется значение дисперсии σ2 = 1.
Следует отметить, что плотность (4.14) является частным случаем плотности
 u2 (t )
u′(t0 )

0 
f (t0 ) = ±
exp−
,

2 
2π

полученной в первой части пособия (разд. 3, формула (3.9)) для общего случая (стационарного нормального процесса). Это естественно, так как модель гауссова марковского процесса – частный случай
общей нормальной модели.
4.4. Оценивание времени прихода импульсного сигнала
Пусть на входе апериодического звена с весовой функцией
h (t) = α exp(−αt)
прямоугольный сигнал длительностью T, маскируемый белым шумом. Сигнал на выходе
44
t


 Aα ∫ exp(−ατ )dτ = A (1 − exp(−αt)), t ≤ T,
s(t) = 
0

 A (1 − exp(−αt)− exp(−α (t − T ))), t > T,

(4.15)
A – амплитуда сигнала. Функция корреляции шума x(t) на выходе
в установившемся режиме
∞
∞
0
0
R (τ ) = ∫ h (t)h (t + τ )dt = α2 exp(−ατ )∫ exp(−2t)dt =
α2
exp(−ατ ).
2
Коэффициент α2/2 можно положить равным единице, тогда дисперсия шума σ2 = 1, отношение сигнал – шум d = A/σ = A. Функция
корреляции
R (τ )= exp(−ατ ) (4.16)
показывает, что гауссов шум – марковский процесс первого порядка.
Время прихода сигнала измеряется методом временной фиксации [13, 14] как момент пересечения переднего или заднего фронта
сигнала (суммы z(t) = x(t) + s(t)) с заданным уровнем u:
в среднем
x (t)+ s(t) = u, (4.17)
Μ  x (t) + u = s(t), s(t) = u. (4.18)
Сигнал (4.15) длительностью T = 5, амплитудой A = 5 с коэффициентом a = 2 и a = 10 и траектории шума x(t) с функцией корреляции (4.16) показан на рис. 4.4.
Сигналы различаются крутизной фронтов. Фронты сигнала (без
шума) пересекают уровень u = A/2 в моменты t1 = 0.3466 (снизу
вверх) и t2 = 5.3466 (сверху вниз) при α = 2, t1 = 0.0693 и t2 = 5.0693,
которые характеризуют время прихода сигнала. Можно отметить,
что при немалом отношении сигнал – шум d = 5 могут появляться
ложные пересечения уровня шумом (точки τ1, τ2). Достоверные измерения времени прихода импульсных сигналов возможны при отношении сигнал – шум d > 5.
Траектории z(t) = x(t) + s(t) тех же сигналов с другими реализациями шума показаны на рис. 4.5. Моменты пересечения становятся случайными – время прихода измеряется со случайными погрешностями. Более пологие фронты случайного сигнала могут пересекать уровень неоднократно, вероятность этого растет с умень45
TY
B
V
U
U
TY
U T T
B
V
U
s
Рис. 4.4. Сигнал, траектории шума, u = A/2, T = 5
[
B
U
U
[
U B
U
U
U
Рис. 4.5. Сигналы, маскируемые шумом
шенем крутизны фронтов. Для того чтобы избежать неоднозначности,
время прихода принято фиксировать по первому пересечению.
Уравнение (4.17) , записанное в виде
v (t) = x (t)+ s(t) − u = 0, (4.19)
задает моменты t1 и t2 как точки пересечения траекторий ν(t) с нулевым уровнем (рис. 4.6, а).
46
W
¸
U
U
s
¹
JOE
U
U
s
Рис. 4.6. Эквивалентные траектории
Можно рассчитывать и моделировать плотность распределения
(4.14) времени пересечения f(t1) переднего и f(t2) заднего фронта сигнала (передним фронтом считается сигнал (4.15) на интервале t ≤ T,
задним – на интервале t > T). Но в силу симметрии фронтов и стационарности шума достаточно получить одну из плотностей – другая будет такой же, сдвинутой на время T.
Среднее значение переднего фронта в (4.14) – сигнал
s(t) = A (1 − exp(−αt))= m(t),
производная
m ′(t) = αA exp(−αt),
так что
2


1  


2 


A
1
−
exp
−
α
t
−
(
)
(
)


0

α

2  

f (t0 ) =
exp  A exp(−αt0 )−
. (4.20)


2
2π












Нормированная плотность (4.20) показана на рис. 4.18: 1 – α = 10,
2 – α = 2.
Как и следовало ожидать, более крутой фронт пересекается
с меньшей дисперсией.
Моделирование пересечений включает:
– исключение траекторий, не пересекающих уровень ни разу;
– вычисление индикаторной функции (рис. 4.6, б)
47
G
E
Рис. 4.7. Плотности распределения
U
ind = diff (sign(v)); (4.21)
– построение гистограммы времени первого пересечения уровня.
Индикатор отмечает положительными пиками точки пересечения снизу вверх, сверху вниз – отрицательными. Гистограмма сопоставляется с плотностью распределения.
Для моделирования задач пересечений разработан пакет специальных файл – функций [13, разд. 2]. Без их использования программа вычисления сигнала, переднего фронта, плотности распределения
(4.20), гистограммы времени первого пересечения может иметь вид
N=1000
T=5
U=5
a=2
de=0.05
nn=T/de
t=0:de:2*T;
n=length(t)
s=1-exp(-a*t);
z=a*exp(-a*t);
%
%
%
%
%
длительность сигнала на входе
амплитуда сигнала
показатель экспоненты
интервал дискретизации
число отсчётов сигнала на входе
% общее число отсчётов
% переходная функция
% её производная
ss=zeros(1,n);
zz=ss;
ss(nn+1:n)=s(1:nn+1);
zz(nn+1:n)=z(1:nn+1);
s=U*(s-ss);
48
% сигнал (4.15)
z=U*(z-zz);
% производная сигнала
% ------------------------% генератор шума
r=exp(-a*t);
% функция корреляции
for i=1:n
for j=1:n
b(i,j)=r(abs(i-j)+1);
end
end
B=b;
% корреляционная матрица
[u,v]=eig(B);
x=randn(N,n);
bx=cov(x);
[ux,vx]=eig(bx);
A=u*v^(1/2)*u’*ux*vx^(-1/2)*ux’;
y=A*x’;
Y=y’;
% траектории шума
%-------------------------------------------% передний фронт сигнала
t1=0:de:T-de
n1=length(t1)
s1(1:n1)=s(1:n1)
% передний фронт
z1(1:n1)=z(1:n1)
% его производная
Y1(1:N,1:n1)=Y(1:N,1:n1); % траектории шума, маскирующего передний фронт
f1=z1.*normpdf(t1,s1-U/2,1)% расчёт плотности (4.20)
ss1=trapz(t1,f1)
f1=f1/ss1
% нормированная плотность распределения
% ------------------------------------------for i=1:N
for j=1:n1
Y1(i,j)=Y1(i,j)+s1(j)-U/2;
% функция (4.19)
end
end
for i=1:N
YY(i,:)=diff(sign(Y1(i,:)));
end
% индикатор (4.21)
for i=1:N
49
kl(i)=length(find(YY(i,:)>0));
end
% число положительных выбросов
lk=find(kl==0)
YY(lk,:)=[]; % исключение траекторий, не пересекающих уровень
[N1,n]=size(YY)
for i=1:N1
k(i)=min(find(YY(i,:)>0)); % номера отсчётов первого пересечения
end
N2=length(find(kl==1)) % число траекторий, пересекающих уровень один раз
pause
TT=k*de
tt=0:de:T-2*de;
h=hist(TT,tt)/N1/de
% моменты первого пересечения
% гистограмма
stem(tt-de/2,h)
hold on
plot(t1,f1)
m=trapz(t1,t1.*f1)
% м.о. первого пересечения
si=trapz(t1,(t1-m).^2.*f1) % с.к.о. первого пересечения На рис. 4.8 и 4.9 показаны плотности распределения (4.20) и гистограммы времени первого пересечения при α = 2 и нескольких
GI
GI
B
E
/
B
E
/
U
U
Рис. 4.8. Плотности распределения, гистограммы; N = 1000, ∆ = 0.05
50
GI
GI
B
E
/
B
E
/
U
U
Рис. 4.9. Плотности распределения, гистограммы; N = 1000, ∆ = 0.02
значениях отношения сигнал – шум d. Число N1 – количество траекторий, пересекающих уровень u = A/2 один раз (моделировались
пересечения N = 1000 траекториями).
При недостаточном отношении сигнал – шум (d < 10) значительное число траекторий пересекают уровень неоднократно, поэтому
плотность распределения моментов всех пересечений отличается от
гистограммы времени первого пересечения. Только при увеличении
отношения сигнал – шум до d = 11 почти все траектории пересекают
уровень один раз, и гистограмма хорошо согласуется с плотностью.
Аналогичные характеристики пересечения траекторий с более
крутым фронтом – при α = 10 показаны на рис. 4.10.
GI
GI
B
E
/
B
E
/
U
U
Рис. 4.10. Плотности распределения, гистограммы; N = 1000, Δ = 0.02
51
Таблица 4.1
α
d
2, t1 = 0.3466
5
7
9
10, t1 = 0.0693
11
5
7
9
11
kr
5
7
9
11
25
35
45
55
mt
σt
0.618
0.557
0.462
0.409
0.087
0.078
0.074
0.072
0.245
0.205
0.107
0.065
4.10–3
2.10–3 7. 10–4 7. 10–5
Удовлетворительное согласование плотности и гистограммы наблюдается уже при отношении сигнал – шум d = 7, когда 90 % траекторий пересекают уровень один раз.
Результат очевиден: решающее значение имеет крутизна фронта
сигнала. Крутизна сигнала (4.20) в точке t1 = ln2/α (рис. 4.4) равна
(табл. 4.1)
αd
(4.22)
kr = s′(t = t1 ) =
.
2
Крутизна (4.22) при α = 10 в пять раз больше, чем при α = 2, поэтому удовлетворительные результаты достигаются при меньшем отношении сигнал – шум. Например, при d = 7 смещение оценки времени первого пересечения и с.к.о. оценки незначительны: ∆t = mt –
t1 = 0.009, σ ≈ 0.002.
Таким образом, можно сделать вывод о том, что при дисперсии
шума σ2 = 1 крутизна фронта сигнала kr ≥ 30 достаточна для оценивания времени прихода сигнала с приемлемой погрешностью.
52
5. Обобщение марковской модели пересечений
Плотность распределения (4.14) моментов пересечений гауссова
марковского процесса с уровнем u(t)
f (t0 ) = ν
2

 (u(t0 )− m(t0 )) 

exp−
2


2πσ
2
σ


m ′(t0 )
распространяется на общий случай пересечения уровня стационарным гауссовым процессом с нулевым средним и единичной дисперсией [13, 14]:
 u2 (t)
u′(t)

. f (t0 ) = ±
exp−
(5.1)
2 
2π

Это обобщение обосновывается следующим образом.
Пересечение процесса x(t) с уровнем u(t) в момент времени t0 можно трактовать как преобразование случайного значения x(t0) в случайное значение t0. Преобразование случайных величин – типовая
задача теории вероятностей: плотность распределения величины
y = j(x)
f (y) = ψ ′(y) fx (y), (5.2)
где Ψ = j–1 – обратная функция преобразования (x = Ψ(y)); fx(y) –
плотность распределения f(x), в которой случайная x заменена на
Ψ(y). В такой постановке уровень u(t) – обратная функция Ψ преобразования случайной величины x ∈ Ν (0,1) в случайное время t0, что
и следует из сопоставления плотностей (5.1) и (5.2).
Таким образом, если u(t) – достаточно крутой фронт импульсного
сигнала, описываемый дифференцируемой функцией u(t), то (5.1) –
плотность распределения его временного положения.
Пример 5.1. На входе ФНЧ (фильтра нижних частот) Баттерворта
четвертого порядка [15] с частотой среза F = 4 Гц белый шум с единичной дисперсией и прямоугольный сигнал длительностью T = 1 с.
Весовая функция фильтра
h (t) = 23.2196{exp(−23.2196t)(cos 9.6179t + 2.4142 sin 9.6179t)−
−exp(−9.6179t)(cos 23.2196t + 0.4142 sin 23.2196t)}.
Фрагмент программы вычисляет весовую функцию и сигнал на
выходе фильтра на (рис. 5.1):
53
I
U
4
LS
s
s
U
U
Рис. 5.1. Весовая функция, сигнал на выходе ФНЧ, крутизна сигнала
T=1
del=0.025
t=0:del:T;
n=length(t)
h=23.2196*((exp(-23.2196*t).*(cos(9.6179*t)+2.4142*sin(9.6179*t))-...
exp(-9.6179*t).*(cos(23.2196*t)+0.4142*sin(23.2196*t)))); subplot(3,1,1),plot(t,h,t,zeros(1,n))
s=rectpuls(t,2)
% сигнал на входе
S=conv(s,h)*del
% сигнал на выходе
kr=diff(S)/del % крутизна
tt=0:del:2*T
ttt=0:del:2*T-del
subplot(3,1,2),plot(tt,S,tt,zeros(1,length(tt)))
subplot(3,1,3),plot(ttt,kr,ttt,zeros(1,length(ttt)))
Крутизна фронта kr ≈ 10 в окрестности s ≈ A/2 при отношении
сигнал – шум d = 1. Увеличение отношения сигнал – шум до рабочих значений d > 5 приведет к kr > 50, что в соответствии с выводом
предыдущего раздела гарантирует приемлемость погрешностей
оценивания временного положения сигнала. Действительно, в [13,
разд. 6] приведен пример с ФНЧ Баттерворта, имеющим частоту
среза F = 2 Гц, что уменьшает крутизну фронта в два раза. Тем не
менее, расчетная плотность (4.23) согласуется с гистограммой моделирования пересечений достаточно хорошо даже при d = 5 (рис. 5.2, б).
54
Передний фронт и несколько траекторий шума показан на рис. 5.2, а.
Оценка с. к. о. σˆ t = 0.05.
Пример 5.2 [13]. На входе колебательного звена белый шум с единичной дисперсией и прямоугольный сигнал длительностью T = 1
с амплитудой A = 5. Нормированная весовая функция звена
h (t) =
2σ
α α2 + β2 exp(−αt)sin (βt),
β
(
)
α = 5, β = 2π (собственная частота f0 ≈ 1 Гц). Функция корреляции
шума


α
(5.3)
R (τ ) = exp(−α τ )cos βτ + sin β τ  
β

определяет сигнал на выходе фильтра (рис. 5.3, а)

y(t) ïðè 0 ≤ t < T,
s(t) = 



y(t)− y(t − T ) ïðè t ≥ T, y(t) = A (1 − R (t)).
(5.4)
На рис. 5.3, б показана также производная сигнала
y ′(t) ïðè 0 ≤ t < T,

s′ (t) = 


′
′

y (t)− y (t − T ) ïðè t ≥ T,
y ′(t) = A
α2 + β2
sin βt exp(−αt).
β
¹
IG
¸
49
V
(5.5)
U
U
Рис. 5.2. Пересечение переднего фронта сигнала, d = 5
55
Эти и последующие расчеты примера выполняются программой
Ampl=5
% амплитуда сигнала
a=5;b=2*pi
dt=0.01
T=1
% длительность сигнала
t=0:dt:T;
% ось абсцисс для переднего фронта
tt=0:dt:2*T;
% ось абсцисс для всего сигнала
ttt=T:dt:2*T;
% ось абсцисс для заднего фронта
n=length(t)
nn=length(tt)
r=exp(-a*t).*(cos(b*t)+a/b*sin(b*t));
R=exp(-a*tt).*(cos(b*tt)+a/b*sin(b*tt)); % функция корреляции (5.3)
s=Ampl*(1-r);
% передний фронт сигнала
S=Ampl*(1-R);
ss=zeros(1,nn);
ss(1:nn)=S(1:nn);
ss(n+1:nn)=ss(n+1:nn)-ss(1:n-1)% сигнал (5.4)
s1=ss(n:nn);
% задний фронт сигнала
% производная сигнала (5.3)
spr1=exp(-a*tt).*sin(b*tt)
spr2=-exp(-a*t).*sin(b*t)
spr3=zeros(1,nn);
spr3(n+1:nn)=spr2(1:n-1)
spr=Ampl*(a^2+b^2)/b*(spr1+spr3)
% производная сигнала (5.5)
subplot(2,1,1),plot(tt,ss)
subplot(2,1,2),plot(tt,spr)
pause
% рис. 5.3
[n,Y,db,B]= gener(r,N,t)
YY=Y+Ampl/2;
size(YY)
% генератор траекторий
% траектории шума с м.о. Ampl/2
% рис. 5.4-а
subplot(1,2,1),plot(t,s,t,YY(34,:),t,YY(321,:),t,YY(789,:),t,zeros(1,n
),t,2.5*ones(1,n))
% пересечение переднего фронта
56
% f - плотность распределения (5.1) времени первого пересечения сверху вниз
f=Ampl*(a^2+b^2)/b/sqrt(2*pi)*sin(b*t).*exp(-a*t-(s-Ampl/2).^2/2)
YY=Y+Ampl/2;
s=Ampl*(1-r);
% траектории шума с м.о. Ampl/2
% передний фронт сигнала
[IN,YY,N,NN]=ind(s,YY,dt);
% индикатор пересечений
dtt=0.05;T1=0;T2=1
[tt,h,NN,tc]=cro1(IN,dt,dtt,T1,T2)
% гистограмма времени первого
пересечения
[th,ff,b]=inthist(T1,T2,tt,dtt,h)
% интерполяция гистограммы
bb=b’
subplot(1,2,2),stem(tt,h/dtt) % рис. 5.4-б
hold
plot(t,f,’r’,th,ff)
mt=mean(tc)
% оценка среднего
si2=var(tc)
% оценка дисперсии
Время прихода сигнала s(t) оценивается методом временной фиксации по уровню u = A/2. Сигнал и его производная показаны на
рис. 5.3. На рис. 5.4, а показано пересечение трех траекторий стационарного шума с функцией корреляции (5.3) и переднего фронта
сигнала (5.4). Экспериментальные результаты (гистограмма) приведены на рис. 6.2, б: гистограмма, плотность распределения (5.1) –
кривая 1, интерполированная гистограмма ([13], разд. 2.2, файл –
T
5
"
U
LS
s
U
Рис. 5.3. Сигнал, производная сигнала
57
¸
¹
49
IG
V
U
s
U
Рис. 5.4. Пересечение переднего фронта сигнала, d = 5
функция INTHIST) – кривая 2. Распределение времени пересечения
вследствие асимметрии плотности отлично от нормального. Из
N = 1000 траекторий в типичном эксперименте пересеклись с фронтом NN = 995. Крутизна фронта сигнала kr ≈ 20 оказывается достаточной для приемлемых погрешностей оценивания времени прихода: типичные оценки с. к. о. σˆ t = 0.06.
В заключение следует подчеркнуть независимость плотности
распределения (5.1) от корреляционных свойств шума, т. е. от его
инерционности [13, разд. 3.2]. Этот результат аналогичен числу пересечений уровня: достаточно крутой фронт сигнала пересекается
практически один раз, точно так же инерционность шума пренебрежима при пересечении крутого фронта.
58
Заключение
В учебном пособии на уровне простых соотношений линейной
алгебры изложена теория гауссовых марковских последовательностей первого и высших порядков. Установлено характерное свойство
последовательности k-го порядка – 2k + 1-диагональность матрицы
точности. Получены расчетные соотношения для условных математического ожидания и корреляционной матрицы n значений последовательности k-го порядка при l фиксированных значениях (l ≥ k).
Разработаны алгоритмы аппроксимации немарковского процесса
марковской последовательностью минимального порядка и генерирования марковских последовательностей.
Приведены примеры применения марковских моделей при решении задачи пересечений случайных процессов с неслучайными
уровнями, в частности, для оценивания временного положения импульсных сигналов. Получена плотность распределения времени
пересечения марковского процесса первого порядка с произвольным дифференцируемым уровнем.
Плотность распределения времени пересечения обобщается на
общий случай стационарного гауссова процесса и достаточно крутого уровня (фронта импульсного сигнала).
59
Библиографический список
1. Тихонов В. И., Миронов М. А. Марковские процессы. – М.: Сов.
радио, 1977. – 488 с.
2. Ануфриев И. и др. MATLAB 7. – СПб.: БХВ – Петербург, 2005. –
1104 с.
3. Королюк В. С. и др. Справочник по теории вероятностей и математической статистике. – М.: Наука, 1985. – 640 с.
4. Осипов Л. А., Воробьева Ю. Г. Генерирование гауссовых марковских последовательностей // Информационно – управляющие системы. 2006. № 4 (23). С. 4–9.
5. Воробьев С. Н., Гирина Н. В., Осипов Л. А. Гауссовы марковские
последовательности. Изв. вузов. Приборостроение. 2011. Т. 54. № 1.
С. 23–31.
6. Воробьев С. Н., Осипов Л. А. Линейные системы. Расчет и моделирование / СПбГУАП. СПб., 2004. – 122 с.
7. Воробьев С. Н. Эффективное обнаружение детерминированных
сигналов / СПбГУАП. СПб., 2003. – 139 с.
8. Тихонов В. И. Нелинейные преобразования случайных процессов. – М.: Радио и связь, 1986. – 296 с.
9. Тихонов В. И., Хименко В. И. Проблема пересечений уровней
случайными процессами. Радиофизические приложения. Радиотехника и электроника, 1998. Т. 43. № 5. С. 501–523.
10. Воробьев С. Н., Гирина Н. В. Пересечение стационарных гауссовых последовательностей с неслучайными уровнями // Информационно – управляющие системы, 2009. № 3. С. 7–12.
11. Воробьев С. Н. Пересечение гауссовым марковским процессом
детерминированного уровня // Информационно – управляющие системы. 2004. № 2. С. 16–20.
12. Двайт Г. Б. Таблицы интегралов. – М.: Наука, 1983. – 181 с.
13. Воробьев С. Н., Гирина Н. В., Лазарев И. В., Осипов Л. А. Статистическое моделирование информационных систем, Ч. 1 / ГУАП.
СПб., 2010. – 152 с.
60
14. Воробьев С. Н., Гирина Н. В., Лазарев И. В. Оценивание временного положения импульсного сигнала // Информационно – управляющие системы. 2010. № 4. С. 9–14.
15. Сергиенко А. Б. Цифровая обработка сигналов. – СПб.: Питер,
2002. – 606 с.
61
Содержание
Введение......................................................................... 1. Гауссовы марковские последовательности первого порядка. 1.1. Матрица точности стационарной марковской последовательности................................................................. 1.2. Общие соотношения для марковских последовательностей первого порядка................................................. 2. Гауссовы марковские последовательности высших порядков.
2.1. Матрица точности.................................................. 2.2. Аппроксимация немарковского гауссова процесса марковской последовательностью........................................ 3. Генерирование марковских траекторий............................ 3.1. Генерирование множества реализаций...................... 3.2. Генерирование полубесконечных последовательностей.
3.2.1. Генерирование число за числом.......................... 3.2.2. Генерирование группами чисел.......................... 4. Примеры применения марковских моделей при расчете
информационных систем................................................... 4.1. Пересечение постоянного уровня.............................. 4.2. Пересечение переменного уровня.............................. 4.3. Плотность распределения времени пересечения заданного уровня марковским процессом................................ 4.4. Оценивание времени прихода импульсного сигнала..... 5. Обобщение марковской модели пересечений...................... Заключение..................................................................... Библиографический список............................................... 62
3
6
6
14
18
18
23
26
26
30
30
34
38
39
40
42
44
53
59
60
Учебное издание
Воробьев Станислав Николаевич,
Гирина Наталия Владимировна,
Лазарев Игорь Владимирович,
Осипов Леонид Андроникович
Статистическое моделирование
информационных систем
Часть II
Учебное пособие
Редактор В. П. Зуева
Верстальщик С. В. Барашкова
Сдано в набор 23.03.11. Подписано в печать 18.05.11. Формат 60 × 84 1/16.
Бумага офсетная. Усл. печ. л. 3,5. Уч.-изд. л. 3,8.
Тираж 150 экз. Заказ № 218.
Редакционно-издательский центр ГУАП
190000, Санкт-Петербург, Б. Морская ул., 67
Для заметок
Документ
Категория
Без категории
Просмотров
0
Размер файла
2 219 Кб
Теги
vorobiev
1/--страниц
Пожаловаться на содержимое документа