close

Вход

Забыли?

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

?

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

код для вставкиСкачать
ВЕСТНИК ТОМСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА
2011
Управление, вычислительная техника и информатика
№ 4(17)
ОБРАБОТКА ИНФОРМАЦИИ
УДК 519.21
А.М. Горцев, В.Л. Зуевич
ОПТИМАЛЬНАЯ ОЦЕНКА ПАРАМЕТРОВ
АСИНХРОННОГО ДВАЖДЫ СТОХАСТИЧЕСКОГО ПОТОКА СОБЫТИЙ
С ПРОИЗВОЛЬНЫМ ЧИСЛОМ СОСТОЯНИЙ
Решена задача оптимальной оценки неизвестных параметров асинхронного
дважды стохастического потока с произвольным (конечным) числом состояний. Оценка параметров производится на основе наблюдения за моментами
наступления событий потока. Оценки имеют минимальное среднеквадратическое отклонение от истинных значений параметров потока.
Ключевые слова: асинхронный дважды стохастический поток событий,
оптимальная оценка параметров, апостериорная плотность вектора параметров, цифровые сети интегрального обслуживания.
Важной сферой приложения теории массового обслуживания является проектирование и создание информационно-вычислительных сетей и различных сетей
связи, которые можно объединить единым термином – цифровые сети интегрального обслуживания (ЦСИО). Данная сфера была определена развитием информационных технологий в конце ХХ века. Возникла необходимость в разработке математических моделей потоков событий, адекватно описывающих реальные информационные потоки, функционирующие в ЦСИО. Одними из первых работ в
этом направлении были [1–3]. Отметим, что на практике параметры, характеризующие поток событий, частично либо полностью неизвестны. Кроме того, параметры могут изменяться с течением времени случайным образом, что приводит к
рассмотрению дважды стохастических потоков событий. Поскольку функционирование системы обслуживания непосредственно зависит от параметров входящего потока, важной задачей является оценка в произвольный момент времени его
параметров по наблюдениям за этим потоком. Исследования по оценке параметров дважды стохастических потоков были проведены, например, в работах [4 – 6].
В настоящей статье получен явный вид оценок параметров асинхронного дважды стохастического потока с конечным числом состояний [7]. Оценки оптимальны в смысле минимума среднеквадратического отклонения от истинных значений параметров. Приводится алгоритм оценки параметров асинхронного потока
событий.
1. Постановка задачи
Рассматривается асинхронный дважды стохастический поток событий с произвольным конечным числом состояний (далее асинхронный поток либо просто
поток). Интенсивность потока является кусочно-постоянным случайным процес-
А.М. Горцев, В.Л. Зуевич
26
сом λ(t) с n состояниями: λ1, λ2, … ,λn (λ1 > λ2 > … > λn > 0). Процесс (поток) в
момент времени t находится в i-м состоянии, если λ(t) = λi ( i = 1, n ). В течение
времени пребывания в i-м состоянии поток ведет себя как пуассоновский с интенсивностью λi ( i = 1, n ). Длительность пребывания в i-м состоянии есть экспоненциально распределенная случайная величина с функцией распределения
Fi (τ) = 1 − eαii τ , где αii = −
n
∑
j =1, j ≠ i
αij ( i = 1, n ); αij > 0 ( i, j = 1, n , i ≠ j) – интенсив-
ность перехода процесса λ(t) из состояния i в состояние j, т.е. величины αij образуют матрицу интенсивностей (матрицу инфинитезимальных коэффициентов) переходов между состояниями αij
n
1
. В сделанных предпосылках λ(t) – транзитив-
ный марковский процесс [7].
Значения параметров потока λi, αij ( i, j = 1, n , i ≠ j) неизвестны. Процесс λ(t)
является принципиально ненаблюдаемым. Предполагается, что о потоке известно
только число состояний n и наблюдению доступны только моменты наступления
событий потока t1, t2, … . Необходимо по наблюдениям t1, t2, … оценить параметры потока λi, αij ( i, j = 1, n , i ≠ j) в момент окончания наблюдения за потоком.
Рассматривается стационарный (установившийся) режим функционирования
наблюдаемого потока событий, поэтому переходными процессами на интервале
наблюдения (t0 , t), где t0 – начало наблюдений, t – окончание наблюдений (момент
вынесения решения), пренебрегаем. Без потери общности можно положить t0 = 0.
Обозначим θ = (λ1 ,..., λ n , αij ; i, j = 1, n, i ≠ j ) – вектор неизвестных параметров
потока, θ(t ) – вектор соответствующих оценок параметров в момент времени t.
θk и θk (t ) – k-е компоненты вектора параметров и вектора оценок соответствен-
но. Обозначим N – размерность векторов θ и θ(t ) , N = n2. Обозначим p( θ | t) =
= p( θ | t1, t2, …, tm) – апостериорная плотность распределения вероятностей вектора
параметров θ в момент времени t при условии, что в моменты времени t1, t2, …, tm
(0 < t1 < t2 < … < tm < t) наблюдались события потока. Обозначим Θ = {λ1 > λ2 > …
> λn > 0, αij > 0; i, j = 1, n , i ≠ j} – область значений вектора параметров θ . Будем
использовать оценку θ(t ) , оптимальную в смысле минимума среднеквадратического отклонения от истинного значения вектора θ :
θk (t ) = ∫ θk p (θ | t )d θ ( k = 1, N ),
Θ
где dθ = dθ1dθ2 … dθN, или в векторной форме
θ(t ) = ∫ θp (θ | t )d θ .
(1)
Θ
Выражение (1) дает оптимальную оценку параметров потока в виде апостериорного среднего. Чтобы воспользоваться формулой (1), необходимо найти выражение для апостериорной плотности p( θ | t).
Оптимальная оценка параметров асинхронного дважды стохастического потока событий 27
2. Явный вид апостериорных вероятностей состояний потока
Для нахождения апостериорной плотности p( θ | t), как будет показано в разделе 3 настоящей статьи, необходимо знать ω(λj | t) ( j = 1, n ) – апостериорные вероятности пребывания асинхронного потока в j-м состоянии в произвольный момент
времени t. Явный вид апостериорных вероятностей ω(λj | t) ( j = 0, n ) для любого
момента времени найден в работе [7].
Согласно [7], поведение ω(λj | t) ( j = 1, n ) на полуинтервале [tk , tk+1) (k = 1, 2,
…) между соседними наблюдёнными событиями асинхронного потока и на полуинтервале [t0, t1) между началом наблюдений и наблюдением первого события определяется выражением
n
∑ s jl zl (tk )eω (t −t )
l
ω(λ j | t ) =
k
l =1
n n
( j = 1, n , k = 0, 1, … ),
(2)
∑∑ sil zl (tk )eωl (t −tk )
i =1 l =1
n
n
где tk ≤ t < tk+1, ω(λ j | tk ) = ω(λ j | tk + 0) = ∑ s jl zl (tk ) , zl (tk ) = ∑ sli −1ω(λi | tk + 0) ,
l =1
k = 0, 1, …; ωl ( l = 1, n ) – собственные числа матрицы D =
dll = αll – λl ( i, l = 1, n ); sil – элементы матрицы S = sil
n
1
i =1
n
dil 1
, dil = αli (i ≠ l),
, в которой l-й столбец яв-
ляется собственным вектором, соответствующим собственному числу ωl ;
sli−1 ( l , i = 1, n ) – элементы матрицы S −1 , обратной матрице S. В момент времени tk
(в момент наступления события асинхронного потока) апостериорные вероятности (2) претерпевают разрывы 1-го рода [7] и пересчитываются по формуле
λ j ω(λ j | tk − 0)
( j = 1, n , k = 1, 2, … ),
(3)
ω(λ j | tk + 0) = n
λ
ω
(
λ
|
t
−
0)
∑ i i k
i =1
где ω(λj | tk – 0) вычисляется по формуле (2) в момент t = tk, когда t изменяется в
полуинтервале [tk–1 , tk), соседнем с полуинтервалом [tk , tk+1). В качестве начального значения ω(λj | t0 + 0) = ω(λj | t0 = 0) в (2) выбираются априорные финальные вероятности πj ( j = 1, n ) состояний процесса λ(t) в стационарном режиме (t → ∞),
которые удовлетворяют системе линейных алгебраических уравнений [7]
n
n
i=1
j=1
∑ πi αij = 0 ( j = 1, n − 1 ), ∑ π j = 1 .
(4)
3. Нахождение явного вида апостериорной плотности вероятностей
и оценок параметров потока
Получим рекуррентную формулу для определения апостериорной плотности
вероятностей. Воспользуемся методикой [8]: рассмотрим наблюдения за потоком
через равные достаточно малые промежутки времени длительности ∆t, а затем
устремим ∆t к нулю. Пусть наблюдения за потоком начинаются в момент времени
А.М. Горцев, В.Л. Зуевич
28
t = 0 и время t изменяется дискретно с конечным шагом ∆t: t = k∆t (k = 0,1,…).
Обозначим r(k∆t) – число событий, наблюдённых на полуинтервале времени
[(k − 1)∆t, k∆t) (r(k∆t) = 0,1,…, k = 0,1,…). На полуинтервале [–∆t, 0) наблюдений
не производится, поэтому r(0) можно положить произвольным, например, положим r(0) = 0. Обозначим r (m∆t ) = (r (0), r (∆t ),..., r (m∆t )) – последовательность
значений количества наблюденных событий на временных полуинтервалах
[(k − 1)∆t, k∆t) ( k = 0, m ). Рассмотрим момент времени t, такой, что t = m∆t, t + ∆t =
= (m + 1)∆t. Тогда имеем r(m∆t) = r(t), r((m + 1)∆t) = r(t + ∆t), r (m∆t ) = r (t ) ,
r ((m + 1)∆t ) = r (t + ∆t ) .
Введем p( θ | r (m∆t ) ) = p( θ | r (t ) ) = p( θ | t) – апостериорную плотность вероятностей вектора параметров θ в момент времени t при условии, что на полуинтервалах времени [(k – 1)∆t, k∆t) ( k = 0, m ) наблюдалось r(k∆t) событий потока соответственно. Аналогично введем p( θ | t + ∆t) – апостериорную плотность в момент времени t + ∆t.
Теорема 1. Апостериорная плотность p( θ | t + ∆t) определяется рекуррентной
формулой
n
∑ ω(λ j | t )
j =1
p ( θ | t + ∆t ) = p ( θ | t )
n
(λ j ∆t ) r (t +∆t )
r (t + ∆t )!
ω(λ j | t )
∫ p (θ | t ) ∑
j =1
Θ
e
−λ j ∆t
(λ j ∆t ) r (t +∆t )
r (t + ∆t )!
,
e
−λ j ∆t
(5)
dθ
где ω(λj | t) ( j = 1, n ) – апостериорная вероятность того, что поток в момент времени t (t > t0) находится в j-м состоянии, определяемая формулами (2), (3).
Доказательство. Используя формулу для условной вероятности, находим
p(θ | r (t + ∆t )) =
p(θ, r (t + ∆t ), r (t )) p (r (t + ∆t ) | θ, r (t )) p(θ | r (t )) p (r (t ))
.
=
p(r (t + ∆t ))
p(r (t + ∆t ))
Используя в последней дроби формулу полной вероятности для условной вероятности, получим
p(θ| r (t +∆t )) = p (θ| r (t ))
p (r (t )) n
∑ p(r (t +∆t )| r (t ),θ,λ(t ) =λ j ) p(λ(t ) =λ j | r (t ),θ) . (6)
p(r (t +∆t )) j =1
Рассмотрим в (6) произведение p(r (t + ∆t ) | r (t ), θ, λ (t ) = λ j ) p(λ (t ) = λ j | r (t ), θ) .
Имеем p(λ (t ) = λ j | r (t ), θ) = ω(λ j | t ) ( j = 1, n ), где ω(λj | t) ( j = 1, n ) – апостериорная вероятность того, что поток в момент времени t находится в j-м состоянии,
определяемая в зависимости от момента времени t формулами (2), (3). В силу
предпосылок, количество наблюдённых на полуинтервале [t, t + ∆t) событий не
зависит от последовательности r (t ) наблюдённых до момента t событий, а зависит только от состояния потока в момент времени t, т.е. от значения интенсивности наступления событий потока λ(t) = λj ( j = 1, n ) в момент времени t. Поэтому
имеем
p(r (t + ∆t ) | r (t ), θ, λ (t ) = λ j ) = p (r (t + ∆t ) | θ, λ (t ) = λ j ) . Поскольку асин-
Оптимальная оценка параметров асинхронного дважды стохастического потока событий 29
хронный поток в любом состоянии ведет себя как пуассоновский, получаем
(λ j ∆t ) r (t +∆t ) −λ j ∆t
. Таким образом, находим формулу (6)
p(r (t + ∆t ) | θ, λ (t ) = λ j ) =
e
r (t + ∆t )!
в виде
(λ j ∆t ) r (t +∆t ) −λ j ∆t
p (r (t )) n
(
|
)
.
p(θ | r (t + ∆t )) = p(θ | r (t ))
ω
λ
t
e
∑ j
p (r (t + ∆t )) j =1
r (t + ∆t )!
Находя в последнем выражении множитель
p (r (t ))
из условия нормировки
p(r (t + ∆t ))
∫ p(θ | r (t + ∆t ))d θ = 1 , приходим к (5). Теорема доказана.
Θ
Рассмотрим случай, когда на полуинтервале [t, t + ∆t) нет событий, т.е.
r(t + ∆t) = 0. Это означает, что полуинтервал [t, t + ∆t) находится на временной
оси между моментами наступления соседних событий tk и tk+1 (k = 1, 2, …) либо
между моментом начала наблюдений за потоком t0 и моментом наступления перn
вого события t1. Введем обозначение s (t , θ) = ∑ λ j ω(λ j | t ) .
j =1
Лемма 1. Апостериорная плотность p( θ | t) между моментами наступления событий удовлетворяет интегро-дифференциальному уравнению
⎡
⎤
dp (θ | t )
= − p(θ | t ) ⎢ s (t , θ) − ∫ s (t , θ) p (θ | t )d θ ⎥ (tk < t < tk+1, k = 0, 1, …).
dt
⎣⎢
⎦⎥
Θ
(7)
Доказательство. Поскольку r(t + ∆t) = 0, получаем (5) в виде
n
∑ ω(λ j | t )e
−λ j ∆t
j =1
p ( θ | t + ∆t ) = p ( θ | t )
∫
Θ
n
p(θ | t )∑ ω(λ j | t )e
−λ j ∆t
dθ
j =1
или, разложив экспоненты в числителе и знаменателе в ряд с точностью до o(∆t),
получим
−1
⎡
⎤
p(θ | t + ∆t ) = ⎡⎣ p (θ | t ) − s (t , θ) p (θ | t )∆t + o(∆t ) ⎤⎦ ⎢1 − ∆t ∫ s(t , θ) p (θ | t )d θ + o(∆t ) ⎥ .
⎢⎣
⎥⎦
Θ
Раскладывая второй сомножитель в ряд, получаем
p(θ | t + ∆t ) = p (θ | t ) − ∆ts(t , θ) p (θ | t ) + ∆tp(θ | t ) ∫ s(t , θ) p (θ | t )d θ + o(∆t ) .
Θ
Перенося в последнем равенстве p( θ | t) влево, деля на ∆t и устремляя ∆t к нулю,
приходим к (7). Лемма доказана.
Асинхронный поток обладает свойством ординарности, поскольку в каждом
состоянии ведет себя как пуассоновский. Поэтому вероятность наступления на
полуинтервале [t, t + ∆t) двух и более событий равна o(∆t). Рассмотрим случай,
когда на полуинтервале [t, t + ∆t) наступает одно событие потока в момент времени tk (t < tk < t + ∆t). При этом r(t + ∆t) = 1.
А.М. Горцев, В.Л. Зуевич
30
Лемма 2. Апостериорная плотность p( θ | t) в момент наступления события пересчитывается по формуле
p (θ | tk − 0) s (tk − 0, θ)
p(θ | tk + 0) =
(k = 1, 2, …).
(8)
∫ p(θ | tk − 0) s(tk − 0, θ)d θ
Θ
Доказательство. Поскольку r(t + ∆t) = 1, получаем (5) в следующем виде:
n
∆t ∑ ω(λ j | t )λ j e
−λ j ∆t
j =1
p ( θ | t + ∆t ) = p ( θ | t )
n
∆t ∫ p (θ | t )∑ ω(λ j | t )λ j e
Θ
,
−λ j ∆t
dθ
j =1
или, разложив экспоненты в числителе и знаменателе в ряд с точностью до o(∆t) и
вводя величины ∆t′ и ∆t″ такие, что t = tk – ∆t′, t + ∆t = tk + ∆t″, находим
∆tp (θ | tk − ∆t ′) s (tk − ∆t ′, θ) + o(∆t )
p(θ | tk + ∆t ′′) =
.
∆t ∫ p (θ | tk − ∆t ′) s (tk − ∆t ′, θ)d θ + o(∆t )
Θ
Поделим числитель и знаменатель последней дроби на ∆t, после чего устремим ∆t
к нулю (при этом t = tk – ∆t′ стремится к tk слева, t + ∆t = tk + ∆t″ стремится к tk
справа). В результате предельного перехода приходим к (8). Лемма доказана.
Замечание к лемме 2. В момент времени t0 начала наблюдений за потоком апостериорная плотность p( θ | t0) задается, исходя из априорных данных о параметрах асинхронного потока. Если таких данных нет, можно задать плотность p( θ | t0)
как произведение N плотностей равномерно распределенных случайных величин,
каждая из которых распределена в некотором интервале допустимых значений
для соответствующего параметра потока.
Теорема 2. Поведение апостериорной плотности p( θ | t) на временной оси
определяется интегро-дифференциальным уравнением (7) и формулой пересчета
вероятностей (8), в которых tk ≤ t < tk+1, p( θ | tk) = p( θ | tk + 0) (k = 0,1,…). В момент
времени t0 апостериорная плотность p( θ | t0) задается согласно замечанию к
лемме 2.
Доказательство следует из лемм 1 и 2 путем синхронизации формул (7) и (8).
Теорема доказана.
Теорема 3. Апостериорная плотность вероятностей p( θ | t) вектора параметров
θ на полуинтервале времени [tk, tk+1) (k = 0, 1, …) определяется формулой
⎡ t
⎤
p(θ | tk ) exp ⎢ − ∫ s (τ, θ)d τ ⎥
⎢⎣ tk
⎥⎦
p (θ | t ) =
(tk ≤ t < tk+1, k = 0, 1, …),
⎡ t
⎤
∫ p(θ | tk ) exp ⎢⎢− ∫ s(τ, θ)d τ ⎥⎥ d θ
Θ
⎣ tk
⎦
(9)
где p( θ | tk) = p( θ | tk + 0) вычисляется в момент наступления события tk (k = 1, 2,
…) по формуле (8), а в момент t0 апостериорная плотность p( θ | t0) задается согласно замечанию к лемме 2.
Оптимальная оценка параметров асинхронного дважды стохастического потока событий 31
Доказательство. Преобразуем интегро-дифференциальное уравнение (7) к
следующему виду:
⎡
⎤
dp (θ | t )
= − ⎢ s (t , θ) − ∫ s (t , θ) p (θ | t )d θ ⎥ dt .
p (θ | t )
⎣⎢
⎦⎥
Θ
Интегрируя последнее интегро-дифференциальное уравнение в пределах от tk до t,
получаем
t
⎡
⎤
⎡ p (θ | t ) ⎤
ln ⎢
=
−
⎢ s (τ, θ) − ∫ s (τ, θ) p (θ | τ)d θ ⎥ d τ ,
⎥
∫
⎥⎦
⎣ p (θ | t k ) ⎦
⎣
Θ
tk ⎢
или, проделывая необходимые преобразования, получаем
⎡ t
⎤
⎡t
⎤
p(θ | t ) = p(θ | tk ) exp ⎢ − ∫ s (τ, θ)d τ ⎥ exp ⎢ ∫ ∫ s (τ, θ) p (θ | τ)d θd τ ⎥ .
⎢⎣ tk
⎥⎦
⎢⎣ tk Θ
⎥⎦
Находя последний сомножитель из условия нормировки
∫ p(θ | t )d θ = 1 , приходим
Θ
к (9). Теорема доказана.
Подставляя (9) в (1), получаем явный вид оценок параметров потока для произвольного момента времени наблюдения за потоком. При этом в момент времени
tk+1 (k = 0, 1, …) имеем
θ(tk +1 ) = θ(tk +1 + 0) = ∫ θp (θ | tk +1 + 0)d θ ,
(10)
Θ
где p( θ | tk + 0) вычисляется по формуле (8).
4. Приближенные формулы для расчета оценок параметров потока
Сделаем важное замечание. При расчете вектора оценок параметров θ(t ) по
формулам (1), (8) – (10) возникают существенные сложности реализации алгоритма данного расчета на ЭВМ. Во-первых, интегрирование во всех указанных
формулах ведется по N-мерной области Θ, которая в общем случае может быть не
ограничена. Во-вторых, интегралы являются N-кратными. Поэтому в настоящей
статье приведен приближенный алгоритм расчета оценок θ(t ) , который позволяет
избавиться от указанных сложностей.
Идея приближенного алгоритма состоит в предположении достаточной близости вектора оценок θ(t ) к истинному значению вектора параметров θ , что позволяет воспользоваться разложениями в ряд Тейлора и значительно упростить формулы (1), (8) – (10) для вычислений.
В предположении достаточной близости оценок θ(t ) к начальному значению
θ(tk + 0) (k = 0, 1, …) разложим в интегралах вида ∫ f (⋅, θ) p(θ | t )d θ , входящих в
Θ
состав формул (1), (8) – (10), подынтегральные функции в ряд Тейлора в окрестности точки θ(tk + 0) (k = 0, 1, …), ограничиваясь первыми тремя членами разложения.
А.М. Горцев, В.Л. Зуевич
32
Введем обозначения:
mi (t ) = θi (t ) = ∫ θi p(θ | t )d θ , c jl (t ) = ∫ (θ j − m j (t ))(θl − ml (t )) p (θ | t )d θ ,
Θ
Θ
⎡
⎤
⎡ t
⎤
f (t , tk , θ) = exp ⎢ − ∫ s (τ, θ)d τ ⎥ , fi (t , tk , θ) = θi exp ⎢ − ∫ s (τ, θ)d τ ⎥ ,
⎢⎣ tk
⎥⎦
⎢⎣ tk
⎥⎦
t
(11)
где θi , θi (t ) – элементы векторов θ , θ(t ) , соответственно.
В (11) величины mi (t ) – апостериорные средние параметров θi , i = 1, N . Величины cjl(t) ( j , l = 1, N ) – ковариации компонент вектора параметров θ . Вычисление интегралов в формулах (11) представляет собой сложную задачу в части
построения численного алгоритма на ЭВМ интегрирования по области Θ.
Формула (1), согласно (9), с учетом обозначений (11), запишется в следующем
виде:
⎡ t
⎤
p (θ | tk + 0)exp ⎢− ∫ s (τ, θ)d τ⎥
∫ fi (t ,tk ,θ) p(θ | tk + 0)d θ
⎢⎣ tk
⎥⎦
Θ
mi (t ) = ∫ θi
dθ =
, i = 1, N . (12)
⎡ t
⎤
f (t ,tk , θ) p (θ | tk + 0)d θ
Θ
∫
∫ p(θ | tk + 0)exp ⎢⎢− ∫ s(τ,θ)d τ⎥⎥ d θ Θ
Θ
⎣ tk
⎦
Разложим в (12) функции (11) в ряд Тейлора в окрестности точки
(
= m(tk + 0) , пренебрегая членами порядка малости o ∆θ
2
θ(tk + 0) =
) , ∆θ = θ − m(t
k
+ 0) .
Получим
N
∂f (t , tk , m(tk + 0))
(θ j − m j (tk + 0)) +
∂θ j
j =1
f (t , tk , θ) = f (t , tk , m(tk + 0)) + ∑
+
1
2
∂ 2 f (t , tk , m(tk + 0))
(θ j − m j (tk + 0))(θl − ml (tk + 0)) ,
∑
∂θ j ∂θl
j ,l =1
N
N
∂fi (t , tk , m(tk + 0))
(θ j − m j (tk + 0)) +
∂θ j
j =1
fi (t , tk , θ) = fi (t , tk , m(tk + 0)) + ∑
+
1
2
∂ 2 fi (t , tk , m(tk + 0))
(θ j − m j (tk + 0))(θl − ml (tk + 0)) .
∑
∂θ j ∂θl
j ,l =1
N
Рассмотрим числитель в (12). Учитывая (13) и обозначения (11), имеем
∫ fi (t , tk , θ) p(θ | tk + 0)d θ = fi (t, tk , m(tk + 0)) ∫ p(θ | tk + 0)d θ +
Θ
Θ
N
∂f (t , t , m(tk + 0))
+∑ i k
∫ (θ j − m j (tk + 0)) p(θ | tk + 0)d θ +
∂θ j
j =1
Θ
+
1
2
∂ 2 fi (t , tk , m(tk + 0))
∑
∫ (θ j − m j (tk + 0))(θl − ml (tk + 0)) p(θ | tk + 0)d θ =
∂θ j ∂θl
j ,l =1
Θ
N
(13)
Оптимальная оценка параметров асинхронного дважды стохастического потока событий 33
1
2
= fi (t , tk , m(tk + 0)) +
∂ 2 fi (t , tk , m(tk + 0))
c jl (tk + 0) .
∂θ j ∂θl
j ,l =1
N
∑
(14)
Аналогично найдем знаменатель в (12):
∫
f (t , tk , θ) p(θ | tk + 0)d θ = f (t , tk , m(tk + 0)) +
Θ
1
2
∂ 2 f (t , tk , m(tk + 0))
c jl (tk + 0) .(15)
∂θ j ∂θl
j ,l =1
N
∑
Подставляя (14) и (15) в (12), находим
1
2
∂ 2 fi (t , tk , m(tk + 0))
c jl (tk + 0)
∑
∂θ j ∂θl
j ,l =1
1
f (t , tk , m(tk + 0)) +
2
∂ 2 f (t , tk , m(tk + 0))
c jl (tk + 0)
∑
∂θ j ∂θl
j ,l =1
fi (t , tk , m(tk + 0)) +
mi (t ) =
N
N
, i = 1, N .
(16)
Согласно формуле (16), для оценки параметра θi ( θi = mi (t ) ) в произвольный
момент времени t необходимо знать начальные значения m(tk + 0) , матрицу ковариаций
C (tk + 0) = c jl (tk + 0)
f (t , tk , m(tk + 0)) ,
N
1
,
значения
функций
fi (t , tk , m(tk + 0))
и
i = 1, N . Отметим, что интегрирование при вычислении
fi (t , tk , m(tk + 0)) и f (t , tk , m(tk + 0)) , i = 1, N , производится на конечном отрезке
времени [tk, t], t < tk+1.
Начальные значения m (t0 + 0) и C(t0 + 0) задаются согласно апостериорной
плотности p(θ | t0 + 0) . Если плотность p(θ | t0 + 0) выбрана в соответствии с замечанием к лемме 2, то элементы вектора m (t0 + 0) и элементы матрицы C(t0 + 0)
получаются следующим образом:
mi (t0 + 0) =
c jl (t0 + 0) = δ( j , l )
1
(0)
θ(1)
j −θj
1
θi(1) − θi(0)
θi(1)
∫
θi d θi , i = 1, N ,
θi(0)
θ(1)
j
∫
(θ j − m j (t0 + 0))(θl − ml (t0 + 0))d θ j , j , l = 1, N ,
(17)
θ(0)
j
где ⎡⎣ θi(0) , θi(1) ⎤⎦ – отрезок определения равномерной плотности вероятностей
p(θi | t0 + 0) = 1/(θi(1) − θi(0) ) параметра θi , i = 1, N , δ( j , l ) – символ Кронекера.
Рассмотрим формулу (16) для расчета mi(t) в любой момент времени t на полуинтервалах [tk, tk+1) ( i = 1, N , k = 0, 1, …). Для вычисления величин mi(t) ( i = 1, N ),
во-первых, необходимо получить выражения для векторов m (tk + 0), матриц ковариаций C(tk + 0) при k = 1, 2, … . Bo-вторых, необходимо получить вторые част∂ 2 f (t , tk , m(tk + 0)) ∂ 2 fi (t , tk , m(tk + 0))
,
, i, j , l = 1, N .
ные производные функций
∂θ j ∂θl
∂θ j ∂θl
В силу (11), аналитические формулы вторых частных производных содержат определенные интегралы, поэтому данные производные рассчитываются численно.
А.М. Горцев, В.Л. Зуевич
34
Получим выражения для векторов m (tk + 0) и матриц ковариаций C(tk + 0) при
k = 1, 2, … . Для этого рассмотрим m (tk + 0) (k = 1, 2, …). С учетом выражений
(11) и (8), находим mi(tk + 0) в виде
∫ θi s(tk − 0, θ) p(θ | tk − 0)d θ
Θ
mi (tk + 0) = ∫ θi p (θ | tk + 0)d θ =
,
∫ s(tk − 0, θ) p(θ | tk − 0)d θ
Θ
Θ
i = 1, N , k = 1, 2, … .
Обозначим
si (tk − 0, θ) = θi s (tk − 0, θ)
и
разложим
(18)
s (tk − 0, θ) ,
функции
si (tk − 0, θ) ( i = 1, N ) в ряд Тейлора в окрестности θ = m(tk − 0) с точностью до
(
o ∆θ
2
) , ∆θ = θ − m(t
k
− 0) . Получим
N
∂s (tk − 0, m(tk − 0))
(θ j − m j (tk − 0)) +
∂θ j
j =1
s (tk − 0, θ) = s (tk − 0, m(tk − 0)) + ∑
+
1
2
∂ 2 s (tk − 0, m(tk − 0))
(θ j − m j (tk + 0))(θl − ml (tk + 0)) ,
∂θ j ∂θl
j ,l =1
N
∑
N
∂si (tk − 0, m(tk − 0))
(θ j − m j (tk − 0)) +
∂θ j
j =1
si (tk − 0, θ) = si (tk − 0, m(tk − 0)) + ∑
+
∂ 2 si (tk − 0, m(tk − 0))
(θ j − m j (tk − 0))(θl − ml (tk − 0)) , i = 1, N .
∑
∂θ j ∂θl
j ,l =1
N
1
2
Подставляя полученные выражения в (18), получаем
1
2
∂ 2 si (tk − 0, m(tk − 0))
c jl
∂θ j ∂θl
j ,l =1
1
s (tk − 0, m(tk − 0)) +
2
∂ 2 s (tk − 0, m(tk − 0))
c jl
∑
∂θ j ∂θl
j ,l =1
si (tk − 0, m(tk − 0)) +
mi (tk + 0) =
N
∑
N
,
j , l = 1, N , k = 1, 2, …,
(19)
где cjl = cjl(tk – 0).
Получим выражение для cjl = cjl(tk – 0) в (18) ( j , l = 1, N , k = 1, 2, …). Рассмотрим (11) для cjl(t) на полуинтервале времени [tk–1, tk) (k = 1, 2, …). Учитывая формулу (9) для плотности p(θ | t ) и обозначения (11), находим
c jl (t ) = ∫ (θ j − m j (t ))(θl − ml (t ))
Θ
f (t , tk −1 , θ) p(θ | tk −1 + 0)
∫ f (t , tk −1 , θ) p(θ | tk −1 + 0)d θ
dθ .
(20)
Θ
Обозначим
F jl (t , tk −1 , θ) = (θ j − m j (t ))(θl − ml (t )) f (t , tk −1 , θ) , tk–1 ≤ t < tk, j , l = 1, N .
(21)
Оптимальная оценка параметров асинхронного дважды стохастического потока событий 35
С учетом введенного обозначения, (20) на полуинтервале времени [tk-1, tk) примет
вид
c jl (t ) =
∫ Fjl (t , tk −1 , θ) p(θ | tk −1 + 0)d θ
Θ
∫ f (t , tk −1 , θ) p(θ | tk −1 + 0)d θ
, j , l = 1, N .
(22)
Θ
Подставляя в (22) t = tk – 0, получаем
c jl (tk − 0) =
∫ F jl (tk − 0, tk −1 , θ) p(θ | tk −1 + 0)d θ
Θ
∫
f (tk − 0, tk −1 , θ) p(θ | tk −1 + 0)d θ
, j , l = 1, N .
(23)
Θ
Раскладывая функцию (21) в ряд Тейлора в окрестности точки θ = m(tk −1 + 0) с
(
точностью до o ∆θ
2
) , ∆θ = θ − m(t
k −1
+ 0) , находим
N
∂F jl (t , tk −1 , m(tk −1 + 0))
u =1
∂θu
F jl (t , tk −1 , θ) = F jl (t , tk −1 , m(tk −1 + 0)) + ∑
+
(θu − mu (tk −1 + 0)) +
2
1 N ∂ F jl (t , tk −1 , m(tk −1 + 0))
(θu − mu (tk −1 + 0))(θv − mv (tk −1 + 0)) , j , l = 1, N .
∑
2 u ,v =1
∂θu ∂θv
Раскладывая в ряд Тейлора функцию f (t , tk −1 , θ) , стоящую в знаменателе (22), в
окрестности точки θ = m(tk −1 + 0) и пренебрегая членами порядка малости
(
o ∆θ
∫
2
) , ∆θ = θ − m(t
k −1
+ 0) , получаем
f (t ,tk −1 , θ) p (θ|tk −1 + 0)d θ= f (t ,tk −1 , m(tk −1 + 0)) +
Θ
1 N ∂ 2 f (t ,tk , m(tk −1 + 0))
c jl (tk −1 + 0) .
∑
2 j ,l =1
∂θ j ∂θl
Подставив полученные разложения в (23) и учитывая, что t = tk – 0, получаем (23)
в виде
c jl (tk − 0) =
2
1 N ∂ F jl (tk − 0, tk −1 , m(tk −1 + 0))
F jl (tk − 0, tk −1 , m(tk −1 + 0)) + ∑
cuv
2 u ,v =1
∂θu ∂θv
f (tk − 0, tk −1 , m(tk −1 + 0)) +
1 N ∂ 2 f (tk − 0, tk −1 , m(tk −1 + 0))
cuv
∑
2 u ,v =1
∂θu ∂θv
j , l = 1, N , k = 1, 2, … ,
,
(24)
где cuv = cuv(tk–1 + 0) ( u, v = 1, N ).
Формула (24) дает элементы матрицы ковариаций C(t) в момент времени
t = tk – 0, необходимые для пересчета апостериорных средних mi (tk + 0) ( i = 1, N )
по формуле (19).
Для расчета оценок параметров θi ( t ) = mi (t ) ( i = 1, N ) на полуинтервале времени [tk, tk+1) по формулам (16) необходимо найти элементы cjl(tk + 0) ковариаци-
А.М. Горцев, В.Л. Зуевич
36
онной матрицы C(t) в момент времени t = tk + 0 ( j , l = 1, N ). Имеем
c jl (tk + 0) = ∫ (θ j − m j (tk + 0))(θl − ml (tk + 0)) p(θ | tk + 0)d θ , j , l = 1, N .
Θ
С учетом введенных обозначений и формулы (8) пересчета плотности вероятностей p(θ | t ) в момент времени tk наступления события имеем
c jl (tk + 0) =
∫ (θ j − m j (tk + 0))(θl − ml (tk + 0))s(tk − 0, θ) p(θ | tk − 0)d θ
Θ
.
∫ s(tk − 0, θ) p(θ | tk − 0)d θ
(25)
Θ
Введем обозначение
A jl (tk + 0, θ) = (θ j − m j (tk + 0))(θl − ml (tk + 0)) s (tk − 0, θ) , j , l = 1, N .
С учетом введенного обозначения, формула (25) принимает следующий вид:
c jl (tk + 0) =
∫ Ajl (tk + 0, θ) p(θ | tk − 0)d θ
Θ
∫ s(tk − 0, θ) p(θ | tk − 0)d θ
, j , l = 1, N .
(26)
Θ
Раскладывая функцию
A jl (tk + 0, θ)
(
θ = m(tk − 0) с точностью до o ∆θ
2
в ряд Тейлора в окрестности точки
) , ∆θ = θ − m(t
− 0) , получаем
N
∂A jl (tk + 0, m(tk − 0))
u =1
∂θu
A jl (tk + 0, θ) = A jl (tk + 0, m(tk − 0)) + ∑
+
k
(θu − mu (tk − 0)) +
2
1 N ∂ A jl (tk + 0, m(tk − 0))
(θu − mu (tk − 0))(θv − mv (tk − 0)) , j , l = 1, N .
∑
2 u ,v =1
∂θu ∂θv
Разлагая функцию s (tk − 0, θ) в ряд в точке θ = m(tk − 0) , находим
N
∂s (tk − 0, m(tk − 0))
(θu − mu (tk − 0)) +
∂θu
u =1
s (tk − 0, θ) = s (tk − 0, m(tk − 0)) + ∑
+
1 N ∂ 2 s (tk − 0, m(tk − 0))
(θu − mu (tk − 0))(θv − mv (tk − 0)) , j , l = 1, N .
∑
2 u ,v =1
∂θu ∂θv
Подставляя данные разложения в (26), получаем формулу для пересчета элементов
cjl(t) матрицы ковариации C(t) в момент времени tk наступления события потока:
A jl (tk + 0, m(tk − 0)) +
c jl (tk + 0) =
2
1 N ∂ A jl (tk + 0, m(tk − 0))
cuv (tk − 0)
∑
2 u ,v =1
∂θu ∂θv
1 N ∂ 2 s (tk − 0, m(tk − 0))
s (tk − 0, m(tk − 0)) + ∑
cuv (tk − 0)
2 u ,v =1
∂θu ∂θv
j , l = 1, N .
,
(27)
Оптимальная оценка параметров асинхронного дважды стохастического потока событий 37
4. Алгоритм оценки параметров потока
Полученные в настоящей статье формулы позволяют сформулировать численный алгоритм расчета оценок параметров асинхронного дважды стохастического
потока событий с конечным числом состояний в произвольный момент времени t
наблюдения за потоком:
1) в момент времени t0 = 0 исходя из априорных данных о потоке задаются начальные значения оценок параметров потока θ(t0 ) = m(t0 + 0) и начальная матрица ковариаций С(t0 + 0). В частности, можно считать начальные оценки параметров N независимыми случайными величинами, распределенными равномерно в
интервалах допустимых значений параметров потока; в этом случае m(t0 + 0) и
С(t0 + 0) задаются по формулам (17). После задания m(t0 + 0) по формулам (4)
вычисляются начальные оценки апостериорных вероятностей ω(λj | t0) ( j = 1, n ),
зависящие от m(t0 + 0) .
2) Для k = 0 в любой момент времени t (tk < t < tk+1), где tk+1 – момент наблюдения (k + 1)-го события потока, рассчитываются оценки параметров m(t ) по формулам (16); вероятности ω(λj | t) ( j = 1, n ), входящие в состав формул (16), рассчитываются по формулам (2).
3) В момент времени t = tk+1 наблюдения события асинхронного потока
m(tk +1 − 0) рассчитываются по формулам (16) и C (tk +1 − 0) – по формулам (24).
Вероятности ω(λj | t) ( j = 1, n ) – по формулам (2).
4) В момент времени t = tk+1 по формулам (19) и (27) рассчитываются
m(tk +1 + 0) и C (tk +1 + 0) соответственно. Вероятности ω(λj | tk+1 + 0) ( j = 1, n ) пересчитываются по формулам (3).
5) Шаги 2 – 4 алгоритма повторяются для k = 1 и т.д.
Отметим, что вероятности ω(λj | t) ( j = 1, n ) в любой момент времени t зависят
от значений оценок параметров потока в этот момент времени, так как в каждый
момент времени для расчета апостериорных вероятностей по формулам (2) – (4)
вместо истинных значений параметров потока используются соответствующие
оценки параметров потока m(t ) .
5. Численный эксперимент
Для численного эксперимента по оценке параметров потока на основе алгоритма, приведенного в разделе 4 настоящей статьи, была реализована программа
расчета оценок на языке С++.
Первый этап расчета – имитационное моделирование потока с заданными параметрами.
Второй этап расчета – вычисление оценок параметров, а также апостериорных
вероятностей и элементов матрицы ковариации С по формулам (2) – (4), (16), (17),
(19), (24), (27).
В приведенном ниже примере расчеты были произведены для следующих значений параметров: n =2 (асинхронный поток имеет два состояния), λ1 = 0,5,
λ2 = 0,05, α1 = α12 = 0,08, α2 = α21 = 0,04. При этом начальная плотность распределения параметров потока выбрана равномерной: параметр λ1 равномерно распре-
А.М. Горцев, В.Л. Зуевич
38
делен в интервале [0,337;0,836], параметр λ2 – в интервале [0,021;0,073], параметр
α1 – в интервале [0,04;0,091], параметр α2 – в интервале [0,031;0,048]. Время моделирования T изменяется от 0 до 100 ед. времени.
На рис. 1 – 4 приведены траектории изменения оценок параметров потока в зависимости от времени моделирования T.
λ1(t)
0,7
0,6
0,5
0
20
40
60
80
T
80
T
Рис. 1. Траектория оценки параметра λ1 (t )
λ2(t)
0,052
0,050
0,048
0
20
40
60
Рис. 2. Траектория оценки параметра λ 2 (t )
α1(t)
0,08
0,07
0,06
0
20
40
60
Рис. 3. Траектория оценки параметра α1 (t )
80
T
Оптимальная оценка параметров асинхронного дважды стохастического потока событий 39
α2(t)
0,0400
0,0398
0,0396
0
20
40
60
80
T
Рис. 4. Траектория оценки параметра α 2 (t )
Скачки траекторий соответствуют моментам наступления событий потока – в
эти моменты оценки параметров пересчитываются по формулам (19) и претерпевают разрывы первого рода.
Анализ приведенного в настоящей статье, а также многочисленных проведенных численных экспериментов показывает, что оценки параметров потока
имеют достаточно хорошую стабильность, в том числе в случае равномерного
распределения начальных оценок параметров в интервалах, допускающих отклонение от истинного значения параметра на 100 %. Колебание оценок для
λ1 (t ) происходит в диапазоне от 0,578281 до 0,737368; для λ 2 (t ) – в диапазоне
от 0,04702 до 0,0486; для α1 (t ) – в диапазоне от 0,06279 до 0,06572; для α 2 (t ) –
в диапазоне от 0,0396 до 0,03968. Еще раз подчеркнем, что полученные оценки
обеспечивают минимум среднеквадратического отклонения оценок параметров
от истинных значений.
Заключение
Полученные в статье результаты позволяют находить оптимальные оценки параметров асинхронного дважды стохастического потока событий с произвольным
конечным числом состояний. Оценки оптимальны в смысле минимума среднеквадратического отклонения от истинных значений параметров. Оценивание параметров потока событий необходимо для управления системами и сетями массового обслуживания, поскольку на практике параметры функционирующих в системах и сетях потоков часто неизвестны.
Полученный в статье явный вид оценок параметров позволяет находить оценки с высокой точностью и скоростью без привлечения численных методов. Для
упрощения численной оценки параметров реальных потоков событий предложен
приближенный алгоритм; проведенные эксперименты показывают возможность
оценки параметров асинхронного потока с помощью данного алгоритма в режиме
реального времени.
ЛИТЕРАТУРА
1. Башарин Г.П., Кокотушкин В.А., Наумов В.А. О методе эквивалентных замен расчета
фрагментов сетей связи. Ч.1 // Изв. АН СССР. Техн. кибернетика. 1979. № 6. С. 92−99.
40
А.М. Горцев, В.Л. Зуевич
2. Башарин Г.П., Кокотушкин В.А., Наумов В.А. О методе эквивалентных замен расчета
фрагментов сетей связи. Ч.2 // Изв. АН СССР. Техн. кибернетика. 1980. № 1. С. 55−61.
3. Neuts M.F. A versatile Markov point process // J. Appl. Probab. 1979. V. 16. P. 764−779.
4. Горцев А.М., Нежельская Л.А. Оценка параметров синхронного альтернирующего пуассоновского потока событий методом моментов // Радиотехника. 1995. № 7−8. С. 6−10.
5. Горцев А.М., Шмырин И.С. Оптимальная оценка параметров дважды стохастического
пуассоновского потока событий при наличии ошибок в измерениях моментов наступления событий // Изв. вузов. Физика. 1999. № 4. С. 19−27.
6. Бушланов И.В., Горцев А.М., Нежельская Л.А. Оценка параметров синхронного дважды
стохастического потока событий // Автоматика и телемеханика. 2008. № 9. С. 76−93.
7. Горцев А.М., Зуевич В.Л. Оптимальная оценка состояний асинхронного дважды стохастического потока событий с произвольным числом состояний // Вестник Томского государственного университета. Управление, вычислительная техника и информатика. 2010.
№ 2(11). С. 44–65.
8. Хазен Э.М. Методы оптимальных статистических решений и задачи оптимального
управления. М.: Сов. радио, 1968. 256 с.
Горцев Александр Михайлович
Зуевич Владимир Леонидович
Томский государственный университет
E-mail: amg@fpmk.tsu.ru; zuevichv@ya.ru
Поступила в редакцию 23 августа 2011 г.
1/--страниц
Пожаловаться на содержимое документа