論文

EM アルゴ リズムによる MMPP パラメータ推定について
山田 雅裕†

中川 健治†† a)

Parameter Estimation for MMPP by EM Algorithm
Masahiro YAMADA† and Kenji NAKAGAWA†† a)

あらまし MMPP（ Markov Modulated Poisson Process ）は，ATM ネットワークにおける画像や音声等のバースト性をもつセルトラヒックに対する近似モデルとして用いられる．そして，MMPP のパラメータである
到着率，状態遷移率を正確に決定することは ATM ネットワークでの有効なトラヒック制御を可能とする．本研
究では，EM アルゴ リズムを用いた MMPP パラメータ推定問題について検討する．従来の研究では，連続時間
MMPP への EM アルゴ リズムの直接の適用 [1] や，MMBP（ Markov Modulated Bernoulli Process ）にモデル変換しての EM アルゴ リズムの適用 [2] があるが ，これらはいずれも計算量が多く，実行時間が長いという問題がある．本研究では複雑な計算を避け，高速なパラメータ推定を行うために，新たなモデルとして MMGP
（ Markov Modulated Geometric Process ）を提案し ，MMPP から MMGP へモデル変換して EM アルゴ リズ
ムを適用する方法を提案する．シミュレーションの結果，パラメータ推定精度は従来法と同程度で，計算時間が従来法 [1] の約 7 分の 1，従来法 [2] の約 115 分の 1 以下となることを確認した．
キーワード

MMPP，EM アルゴ リズム，パラメータ推定，MMGP，セル廃棄率

1. ま え が き

MMPP のモーメントが一致するようにパラメータをMMPP（ Markov Modulated Poisson Process ）は，推定精度が悪いという点が文献 [2] で指摘されている．

決定する方法である．しかし ，この方法はパラメータ

複数の状態間を連続時間マルコフ連鎖に従って状態

そこで，精度の高いパラメータ推定を行うために，ゆが遷移し ，各々の状態でポアソン過程に従って事象がう度関数をもととし た EM アルゴ リズムによるパラ

到着するモデルである．到着率の高い状態と低い状

メータ推定法が提案されている [1], [2]．

態を考え ることによりバースト 性を表現できること

Rydén [1] では，MMPP へ EM アルゴ リズムを直から，ATM（ Asynchronous Transfer Mode ）ネット

接的に適用することによりパラメータ推定を行う方法

ワークにおける画像や音声等のバースト性をもつセルが提案されている．この方法は，事象の到着時間間隔

トラヒックに対する近似モデルとして用いられる．そを観測データとし ，連続時間マルコフ連鎖の状態遷移

して，MMPP のパラメータである到着率，状態遷移を非観測データとして用い，観測データと非観測デー率を正確に決定することは ATM ネットワークでの有タを合わせて完全データと考える．そして完全データ

効なトラヒック制御を可能とする．本研究では，EMの対数ゆう度関数を立て，EM アルゴ リズムを適用し（ Expectation Maximization ）アルゴ リズムを用いたてパラメータ推定を行う．しかし MMPP は連続時間

MMPP パラメータ推定問題について検討する．
過去の研究では ，モー メント マッチング 法による確率過程であるため，直接的な適用には行列の固有値の計算や数値積分等複雑な数値計算を伴うという問題

MMPP パラメータ推定法が Heﬀes ら [3] により提案点がある．

されている．これは，入力トラヒックのモーメントと

ま た Deng ら [2] で は ，モ デ ル を MMPP か ら
MMBP（ Markov Modulated Bernoulli Process ）へ
変換し ，EM アルゴ リズムを適用するパラメータ推定法が提案されている．ここで MMBP とは，事象の到†株式会社エフ・デ ィー・シー，土浦市
FDC Co.,Ltd., Tsuchiura-shi, 300–1173 Japan

††

長岡技術科学大学，長岡市
Department of Electrical Engineering, Nagaoka Univercityof Technology, Nagaoka-shi, 940–2188 Japana) E-mail: nakagawa@nagaokaut.ac.jp

688

着がベルヌーイ過程，状態遷移が離散時間マルコフ連鎖に従うモデルである．この方法では，まず MMPP

電子情報通信学会論文誌 B Vol. J84–B No. 4

pp. 688–696 2001 年 4 月

論文／EM アルゴ リズムによる MMPP パラメータ推定についてから MMBP への観測データの変換を行う．微小時間h を定め，h により到着時間間隔の時間軸を離散化し，
各 h ご とに 2 値データを生成する．すなわち，長さh のタイムスロット中に事象の到着があれば 1，事象の到着がなければ 0 を生成する．その結果，到着時間
間隔系列は 2 値データ系列に変換され，EM アルゴ リ

をいう．

2
−σ11
6
6
6 σ21
Σ=6
.
6
4 ..

い，両者を合わせて完全データと考え，完全データの

σ1R

−σ22
..
.

···
..
.

σ2R
..
.

σR2

···

−σRR

R
X

ただし ，σii =

3
7
7
7
7
7
5

(1)

σij である．

j=1,j=i

行った後，2 値データ系列を観測データとし ，離散時
間マルコフ連鎖の状態遷移を非観測デ ータとし て用

···

σR1

ズムの観測データとなる．次に MMPP から MMBP
へのパラメータの変換を行う．これらのモデル変換を

σ12

3. EM アルゴリズム

対数ゆ う度関数を立て，EM アルゴ リズムを適用しEM アルゴ リズムは，観測できないデータを含む場て MMBP パラメータ推定を行う．EM アルゴ リズム合のモデルパラメータ推定を，最ゆう推定法によってでのパラメータ推定終了後，もとの MMPP へモデル行う繰返しアルゴ リズムで，統計学では古くから知らを逆変換する．しかしこの方法において，各タイムス

れている手法であり，Dempster ら [4] によって一般的

ロットごとに観測データ（ 1 or 0 ）をもつため，観測に定式化がなされた．EM アルゴ リズムは，観測デー

データ量はもとの到着時間間隔データ量と比較して非タのみを用いて最ゆう推定を行うことが困難である場

常に多くなり，メモリ消費量が多い．また，各タイム

合に，観測データと非観測データを合わせて完全デー

スロットごとに状態遷移が起こるかど うかを調べるた

タと考え，完全データの最ゆう推定を間接的に利用しめ，計算量が多いという問題点がある．

て観測データに基づく最ゆう推定を行おうとするもの

本研究ではこれらの問題を解決し ，高速かつ高精である．

度なパラ メータ推定を行うために ，新たなモデ ルと

EM アルゴ リズムは ，E-step（ Expectation-step ）し て MMGP（ Markov Modulated Geometric Pro-と M-step（ Maximization-step ）と呼ばれる二つのス

cess ）を提案し ，MMPP から MMGP にモデル変換

テップを繰り返すことによりパラメータを逐次更新す

して EM アルゴ リズムを適用する方法を提案する．シる方法である．実行回数 p = 1, 2, · · · に対して次の二ミュレーションによる従来法 [1], [2] との比較を行い，つのステップを繰り返す．

提案法の有効性を示す．

2. MMPP モデル
MMPP は，複数の状態間を連続時間マルコフ連鎖に従って状態が遷移し ，各々の状態でポアソン過程に従って事象が到着するモデルである．また MMPP は，【 E-step 】 観測データ y とパラメータ (p) が与えられたという条件のもとで，完全データの対数ゆう度関数 log f (x| ) の非観測データについての条件付き期待値を求める．つまり，Q( | (p) ) = E [ log f (x| )|y ,  (p) ](2)出力系列が与えられても状態遷移系列が一意に決まらを計算する．

ないモデルであり，HMM（ Hidden Markov Model ）【 M-step 】 E-step で求めた対数ゆう度関数の条件の一例である．到着率の高い状態と低い状態を考える

付き期待値を最大化するパラメータ  = (p+1) を求ことによりバースト性を表現できることから，ATM

める．つまり，

ネットワークにおける画像や音声等のバースト性をもつセルト ラヒックに対する近似モデルとし て用いられる．

 (p+1) = arg max Q(| (p) )(3)を計算する．

R 状態 MMPP は，状態 i から j へ状態遷移率 σij

MMPP パラメータ推定への EM アルゴ リズムの適をもつ連続時間マルコフ連鎖に従って状態が遷移し ，用において，事象の到着時間間隔を観測データとして各状態 i において到着率 λi のポアソン過程に従って用い，連続時間マルコフ連鎖の状態遷移を非観測デー

事象が到着するモデルである．ここで，状態遷移率とタとして用いる．MMPP パラメータ推定は観測でき

は式 (1) のような無限小生成行列 Σ の ij 要素のことないデータ，すなわち連続時間マルコフ連鎖の状態遷689

電子情報通信学会論文誌 2001/4 Vol. J84–B No. 4

移，を含むモデルパラメータ推定問題であるため，パ

従来法 [2] ではデータ数は T /h となる．後に示す例でラメータ推定手法には EM アルゴ リズムが有効である．

は a = 0.55，h = 10−3 程度の値なので，従来法 [2]4. EM アルゴリズムによる MMPP
パラメータ推定法の提案のデータ数は非常に大きいものとなる．
また，従来法 [2] では 1 タイムスロットごとに状態遷移が起こるかど うかを調べるため，計算量が多くな

本研究では，複雑な計算を避け，高速なパラメータ

るが，提案法では事象の到着時のみ状態遷移を考える

推定を行うために，新たに MMGP（ Markov Modu-

ため，状態遷移回数が少なくてすみ，計算量をおさえlated Geometric Process ）モデルを提案し ，計算の複雑な連続時間 MMPP から 計算の容易な離散時間
MMGP へのモデルの変換を行い，MMGP に EM ア

ることができる．

ルゴ リズムを適用し てパラ メータ推定を行う．まずMMGP モデル，MMPP から MMGP へのモデル変換について示し，その後 EM アルゴ リズムを用いたパ
ラメータ推定法について示す．

4. 1 MMGP モデル
MMGP は，状態 i から j へ状態遷移確率 rij の離散時間マルコフ連鎖に従って状態が遷移し ，各状態 iにおいて 1 タイムスロット中の事象到着確率 bi の幾

4. 3 パラメータ変換
MMPP と MMGP のパラメータ変換を次式で示す．
bi = λi h for i = 1, · · · , R
σij
rij =
for i, j = 1, · · · , R, i = j
λi

(4)
(5)

ここで，bi は MMGP の状態 i で 1 タイムスロット中に事象が到着する確率，rij は MMGP において事象の到着時に状態 i から j への状態遷移確率を表し，λiは MMPP の状態 i における到着率，σij は MMPPの状態 i から j への状態遷移率を表す．式 (4) の両辺

何分布に従って事象が到着するモデルである．MMGPは，1 タイムスロット中に事象が到着する確率，式 (5)モデルでは，状態遷移が起こるのは事象の到着時のみの両辺は，一つの到着当りの状態 i から j への遷移確である．2 状態 MMGP の状態遷移図を図 1 に示す．

4. 2 観測データ変換
MMPP の状態数 R は既知と仮定し，与えられた観測データ系列，すなわち事象の到着時間間隔の系列をyN
1 = {y1 , y2 , · · · , yN } と示す．まず，長さ h のタイムスロットを定め，時間軸をタイムスロットによって離散化する．そして，到着時間間隔をタイムスロットで数えた系列 lN．
1 = {l1 , l2 , · · · , lN } を生成する（ 図 2 ）提案法では離散化後の観測デ ータ量は 離散化前の観測データ量と変わらず N 個である．しかし従来法 [2]では，タイムスロットごとに観測データ，すなわち 2
値データをもつため，離散化前の観測データ量よりも増加する．観測期間を T とし ，平均到着間隔を a とすると，観測期間中の観測データ量 N の平均値は約率である．また，rii = 1 −RXrij とする．

j=1,j=i

まず，MMGP パラメータ初期値を決定する際に式(4)，(5) のパラメータ変換式を用いる．具体的には，与えられた到着時間間隔の観測値から，k-means クラ
スタリング法 [2] を用いて MMPP パラメータの初期値を推定し ，得られた MMPP パラメータ初期値を式(4)，(5) に代入して，MMGP パラメータの初期値を決定する．k-means クラスタリング法は，推定精度は低いが，非常に高速なパラメータ推定を行えるという特徴があり，MMPP パラメータの初期値決定法として用いる．式 (5) において，σij が λi よりも大きい，すなわち rij > 1 となるような初期値を与えたとしてT /a である．これが提案法のデータ数である．一方，

図 1 2 状態 MMGP モデル
Fig. 1 2-state MMGP model.

690

Fig. 2図 2 観測データ変換
Transformation of observed data.

論文／EM アルゴ リズムによる MMPP パラメータ推定について

も，EM アルゴ リズムでは rij を，後述の式 (12) に+よって求めるため，1 回の繰返し以降 rij ≤ 1 として

N
R X
X

γn (i)[(ln − 1) log{1 − bi } + log bi ]i=1 n=1

推定される．ただし ，rij > 1 のときには，rii < 0 と(8)なってしまうため，この場合には rii の符号を反転したものを初期値とする．また，λi = 0 という初期値がここで，添字 (p) は繰返し回数が p 回目であることを得られた場合には，式 (5) の右辺に注意して，λi = 0示す．また，となる適当な値を初期値とする．

EM アルゴ リズムで計算を行った後，MMGP パラ
メータ推定値 bi ，rij を式 (4)，(5) に代入し，MMPP
パラメータ推定値 λi ，σij へと逆変換する．
4. 4 EM アルゴリズムによるパラメータ推定
R 状態 MMGP において，観測データ系列 lN
1 から
EM アルゴ リズムを適用してパラメータ bi ，rij を推定する．MMGP パラメータ推定への EM アルゴ リズ
ムの適用において，非観測データとして離散時間マル= (s1 , s2 , · · · , sN ) を用いる．

コフ連鎖の状態遷移

ここで，sk は k 番目の状態を表す．観測データと非観測データを合わせて完全データと考える．
【 E-step 】 観 測デ ー タ l

N
1

b

(p)

=

r

(p)
(bi )i=1,··· ,R ， (p)と 現 在 の パ ラ メー タ(p)

= (rij )i,j=1,··· ,R が 与(p)ξn (i, j) = P {sn = i, sn+1 = j|lN
, r (p) }1 ,b(p)
γn (i) = P {sn = i|lN
, r (p) }1 ,bである．ξn (i, j)，γn (i) はそれぞれ forward-backward
アルゴ リズム [8] より計算される．
【 M-step 】 Q(b, r|b(p) , r(p) ) をパラメータ bi ，rij
について最大化する．

∂Q(b, r|b(p) , r(p) )
=0
∂bi

(9)

∂Q(b, r|b(p) , r(p) )
=0
∂rij

(10)

より，
N
X

えられたときの完全データの対数ゆう度関数の非観測
デ ータ，すなわち離散時間マルコフ連鎖の状態遷移についての条件付き期待値 Q(b, r |b(p) , r (p) ) を求

(p+1)

bi

=

める．

l

|b, r} =N
Y(1 − bsn )となり，両辺の対数をとって[(ln − 1) log(1 − bsn )n=1

Q(b, r|b(p) , r(p) )

=

=

R XX
X

|b , r lb ,rN

j=1 i=1 n=1=

ξn (i, j)

n=1
N−1

X

γn (i)

fori, j = 1, · · · , R(12)

パラメータから MMPP パラメータへ逆変換する．
提案法は従来法 [1] で必要とする行列の固有値や数値積分等の複雑な計算を必要とせず，forward-backward
アルゴ リズムによる簡単な計算のみで実行が可能であ

}| N
1 ,

(p)

(p)

]

N
(p)
log P {lN
, r(p) }1 , |b, r }P { |l1 , b
R

(p+1)

rij

パラメータが収束するまで E-step，M-step を繰りとなる．したがって E-step は次式の計算となる．

l

(11)

返す．パラメータ収束後，式 (4)，(5) を用いて MMGP

(7)

X

i = 1, · · · , R

n=1+ log bsn + log rsn sn+1 ]= E[log P { N
1 ,

X

N−1

bsn rsn sn+1
(6)

N
X

γn (i)ln

for
ln −1

n=1

log P {lN
1 , | b, r } =X

n=1

完全データのゆう度関数は

P{ N
1 ,

γn (i)n=1
N

ξn (i, j) log rij

る．また，提案法は従来法 [2] よりも観測データ量や，状態遷移が起こるかど うかを調べる回数が少ないため，少ない計算量，メモリ消費量で実行が可能である．

5. シミュレーション結果と考察
MMPP パラメータ推定のシミュレーション結果に691

電子情報通信学会論文誌 2001/4 Vol. J84–B No. 4

ついて述べる．シミュレーションでは，まずパラメー

上の事象の到着がないように小さな値を用いなければ

タ真値を定めることにより真の MMPP を定め，そのならないが，あまりに小さな値の h を用いると，従来

MMPP から発生する事象の到着時間間隔が与えられ

法 [2] の観測データ量が膨大に増えてしまうため，計算ているとする．この事象の到着時間間隔を観測データが現実的でない．そこで推定精度が変わらない範囲で，として用い，提案法，Rydén による MMPP パラメー

できるだけ大きな h を用いる．case 1 では h = 0.001

タ推定法（従来法 [1] ），Deng による MMPP パラメーとした．

タ推定法（従来法 [2] ）を用いてパラメータ推定を行っパラメータ推定結果を表 2 に示す．推定結果を比較た．そして推定精度，実行時間の比較を行った．推定すると，提案法の λ1 の推定値に推定誤差があるもの

精度の比較には，MMPP パラメータ推定結果をトラの，ほぼ真値に近い値に収束している．従来法 [1]，従ヒック制御に応用するという観点から，MC（ Monte

来法 [2] についてもほぼ 真値に近い値に収束している

Carlo ）法によるセル廃棄率シミュレーションを行い，ことが確認できる．
次に MC 法を用いたセル廃棄率シミュレーションに

セル廃棄率特性での比較を行った．

5. 1 シミュレーションにおける EM アルゴリズムの条件

よるセル廃棄率特性での推定精度の評価を行う．真値及びそれぞれの方法による推定値をサービ ス時間で正シミュレ ーションにおける EM アルゴ リズムの実

規化し ，1 セルのサービ ス時間を単位時間とする．こ行の条件について述べる．まず，EM アルゴ リズムはこで，正規化に用いたサービ ス時間を 0.2 とする．そ繰返しアルゴ リズムであるため，パラメータの初期値して，その正規化し た真値，推定値より各々MMPPを必要とする．今回のシミュレーションでは初期値をを生成し ，セル廃棄率シミュレーションを行った．こk-means クラスタリング法 [2] を用いて決定した．
また，EM アルゴ リズムの終了条件を以下のようにこでシミュレーション回数を 1 × 106 とする．また真定める．

こで λ̄ は次式より得られる．

|b(p) − b(p+1) |b(p)|r (p) − r(p+1) |r(p)≤ 10−5(13)≤ 10−5(14)値，各推定値における平均到着率 λ̄ を表 3 に示す．こλ̄ =λ1 σ21 + λ2 σ12σ12 + σ21(15)

セル廃棄率シミュレーション結果を図 3 に示す．図 3において，TRUE は真値，PROP，CONV1，CONV2

これら二つの式を満たすときに EM アルゴ リズムを終はそれぞれ，提案法，従来法 [1]，従来法 [2] による推了する．

定値を表す．結果を比較すると，三つの方法ともセル

5. 2 シミュレーション結果
5. 2. 1 case 1
シミュレーション条件を表 1 に示す．状態数は既知とし ，パラメータ真値を表 1 中に示す値とし て真の
MMPP を定め，その MMPP から発生する事象の到着
時間間隔 10,000 個を観測データとして用いて MMPP

廃棄率特性は真値に近い値が得られていることがわかパラメータ推定を行った．提案法，従来法 [2] で用いるタイムスロット h は，1 タイムスロット中に二つ以
表 1 シミュレーション条件（ case 1 ）Table 1 Condition of simulation (case 1).
状態数
観測データ数
提案法，従来法 [2] で用いる
タイムスロット h
λ1
λ2
σ12
σ21

692

2
10,000
0.001
10.0
1.0
1.0
0.1

る．そして，提案法は安全側の評価となっているのに対し ，従来法 [1]，従来法 [2] は危険側の評価となっている．
次に実行時間の比較を行う．シミュレーションは並列

Table 2表 2 パラメータ推定値の比較（ case 1 ）Comparison of estimated parameter values(case 1).

λ1
λ2
σ12
σ21

真値 提案法 従来法 [1] 従来法 [2]10.0 10.65
9.77
9.71
1.0
1.098
1.005
1.004
1.0
0.961
0.942
0.941
0.1 0.0847
0.104
0.104

表 3 真値，各推定値における平均到着率 λ̄（ case 1 ）Table 3 True value and estimated values for themean arrival rate λ̄ (case 1).
平均到着率 λ̄

真値
0.364

提案法 従来法 [1] 従来法 [2]
0.374
0.375
0.374

論文／EM アルゴ リズムによる MMPP パラメータ推定について

計算機 SGI Origin 2000，8 GB メモリ，100 GB ディ

し ，提案法では 0.091 と相対的に大きめの値に収束しスク，16 CPUs という環境で行った．実行時間の比較ている．提案法では事象の到着時にのみ状態遷移が起を表 4 に示す．結果を比較すると，提案法は従来法 [1]

こると考えているため，到着率の真値が 0 に近いと，の実行時間の約 1/7，従来法 [2] の実行時間の約 1/115

状態遷移が起こった時間を正しく推定するのが困難とで推定を行っており，提案法が従来法と比較して非常なる．このことは式 (5) の右辺の分母に λi があるこに高速にパラメータ推定を実行していることが確認でとからもわかる．そして，到着した事象がど ちらの状

きる．

態から発生した事象なのかという推定を正確に行うこ5. 2. 2 case 2とが従来法よりも難しい．そのため状態 2 の到着率がシミュレ ーション 条件を表 5 に示す．case 2 では

λ2 = 0，すなわち状態 2 に事象が到着しないという場

0 の場合の推定では，状態 1 での到着を状態 2 での到着であると誤って判断してしまうため，λ2 の推定に誤合のパラメータ推定を行う．

差が生じている．しかし ，以下に示すセル廃棄率でのパラメータ推定結果を表 6 に示す．表 6 の λ2 の推定精度の評価では，提案法のセル廃棄率特性は真値推定結果を比較すると，従来法 [1] が 0.000077，従来と近い値が得られているため，λ2 の推定誤差は，セル

法 [2] が 0.000071 と 0 に近い小さな値であるのに 対

廃棄率特性に大きな影響を与えていないことがわかる．
そこで，セル廃棄率による評価を行う．サービ ス時
間を 0.4 とし，シミュレーション回数を 1 × 106 とする．また，真値，各推定値における平均到着率 λ̄ を
表 7 に示す．
セル廃棄率シミュレーション結果を図 4 に示す．結果を比較すると，提案法，従来法 [1]，従来法 [2] のセル廃棄率特性は，真値に対してやや危険側の評価であるが ，真値と近い特性となっている．このことから，提案法の λ2 の推定での推定誤差は，セル廃棄率特性に大きな影響を与えていないことがわかる．
次に実行時間の比較を行う．実行時間の比較を表 8に示す．結果を比較すると，提案法は従来法 [1] の実行
時間の約 1/8，従来法 [2] の実行時間の約 1/160 で推定を行っており，提案法が従来法と比較して非常に高速にパラメータ推定を実行していることが確認できる．

図 3 MC シ ミュレ ーション に よ る セル 廃 棄 率 推 定 値（ case 1 ）Fig. 3 Estimated values for cell loss rate by MC
simulation (case 1).
表 4 実行時間の比較（ case 1 ）Table 4 Comparison of computation time (case 1).
実行時間

提案法 従来法 [1]9s
65 s

従来法 [2]
1,051 s

表 5 シミュレーション条件（ case 2 ）Table 5 Condition of simulation (case 2).
状態数2
観測データ数
10,000
タイムスロット h 0.001
λ1
10.0
λ2
0.0
σ12
1.0
σ21
0.1

式 (5) では，λi ≥ σij である必要があるが，以上のように λi < σij の場合でも，MMGP に基づく推定

Table 6表 6 パラメータ推定値の比較（ case 2 ）Comparison of estimated parameter values(case 2).

λ1
λ2
σ12
σ21

真値 提案法 従来法 [1]10.0 10.71
9.72
0.0
0.091 0.000077
1.0
1.065
1.033
0.1
0.083
0.101

従来法 [2]9.67
0.000071
1.033
0.102

表 7 真値，各推定値における平均到着率 λ̄（ case 2 ）Table 7 True value and estimated values for themean arrival rate λ̄ (case 2).
平均到着率 λ̄

真値
0.364

提案法 従来法 [1] 従来法 [2]
0.347
0.348
0.347

693

電子情報通信学会論文誌 2001/4 Vol. J84–B No. 4は可能で，真値に近いセル廃棄率特性が得られること

次にセル廃棄率による評価を行う．サービ ス時間を0.01 とし ，シミュレーション回数を 1 × 106 とする．

を確認した．

5. 2. 3 case 3

また，真値，各推定値における平均到着率 λ̄ を表 11

シミュレーション条件を表 9 に示す．case 3 では 3に示す．ここで，3 状態 MMPP の平均到着率 λ̄ は次状態でのパラメータ推定を行う．パラメータ推定結果
を表 10 に示す．
表 10 の推定結果を比較すると，提案法は λ1 ，λ3 ，式で求められる．

λ̄ =λ 1 σ2 σ3 + λ 2 σ1 σ3 + λ 3 σ1 σ2
σ1 σ2 + σ2 σ3 + σ1 σ3σ31 で推定誤差がみられ，また従来法では σ13 で推定誤差が大きくなっている．提案法，従来法ともに推定誤差があるものの，ほぼ真値に近い値に収束している．

3
X

ここで，σi =(16)

σij である．

j=1,j=i

セル廃棄率シミュレーション結果を図 5 に示す．結果を比較すると，提案法のセル廃棄率特性は真値に対して安全側の評価となっているのに対し ，従来法 [1]，従来法 [2] のセル廃棄率特性は，真値に対して危険側の評価となっている．
次に実行時間の比較を行う．実行時間の比較を表 12に示す．パラメータ数が多いため，すべてのパラメー
タが収束条件を満たすまでに多くの繰返し回数を必要とし ，それに伴って 2 状態のときよりも実行時間が長くなっている．しかし ，提案法は従来法 [1] の実行時
間の約 1/40，従来法 [2] の実行時間の約 1/70 で推定を行っており，提案法が従来法と比較して非常に高速にパラメータ推定を実行していることが確認できる．
上記 case 1 ∼ 3 に加えて，更に表 13 に示すパラメータをもつ case 4 ∼ 6 の 2 状態 MMPP について考察する．case 1 ∼ 6 における実行時間と推定精度の図 4 MC シ ミュレ ーション に よ る セル 廃 棄 率 推 定 値（ case 2 ）Fig. 4 Estimated values for cell loss rate by MC
simulation (case 2).

Table 8表 8 実行時間の比較（ case 2 ）Comparison of computation time (case 2).

実行時間

提案法 従来法 [1]
10 s
85 s

従来法 [2]
1,652 s

表 9 シミュレーション条件（ case 3 ）Table 9 Condition of simulation (case 3).
状態数3
観測データ数
10,000
タイムスロット h 0.0001
λ1
100.0
λ2
50.0
λ3
10.0
σ12
5.0
σ13
3.0
σ21
3.0
σ23
2.0
σ31
2.0
σ32
1.0

694

比較を表 14 に示す．表 14 において，時間比 1 及び時間比 2 は，それぞれ，提案法の実行時間の従来法 [1]及び従来法 [2] の実行時間に対する比を表す．また，設表 10 パラメータ推定値の比較（ case 3 ）Table 10 Comparison of estimated parameter values(case 3).
λ1
λ2
λ3
σ12
σ13
σ21
σ23
σ31
σ32

真値
100.0
50.0
10.0
5.0
3.0
3.0
2.0
2.0
1.0

提案法
112.32
53.44
12.77
7.75
2.65
2.42
2.1
1.35
1.39

従来法 [1]102.71
48.83
9.96
7.25
2.22
3.19
1.9
1.75
1.2

従来法 [2]102.41
48.66
9.92
7.33
2.21
3.22
1.91
1.76
1.19

表 11 真値，各推定値における平均到着率 λ̄（ case 3 ）Table 11 True value and estimated values for themean arrival rate λ̄ (case 3).
平均到着率 λ̄

真値
0.392

提案法 従来法 [1] 従来法 [2]
0.401
0.371
0.369

論文／EM アルゴ リズムによる MMPP パラメータ推定について

ズムによる状態滞在確率と状態遷移確率の計算である．観測データ数を N とし ，MMPP の状態数を Rとすると forward-backward アルゴ リズムの計算量はO(N R2 ) である．提案法及び従来法 [2] の観測データ
数を，それぞれ ，N, N  とすると，4. 2 で述べたよ
うに，典型的な場合に N  = 500N である．このことから，従来法 [2] が提案法よりも実行時間がかかることがわかる．状態数 R に関する計算量は提案法，従来法 [2] ともに同程度である．一方，従来法 [1] において最も計算時間がかかるのは ，E-step での状態遷移
回数と状態占有時間の期待値を計算するための数値積
分と行列計算である．観測時間を T ，MMPP の状態
数を R とすると，この部分の計算量は O(T R3 ) である．提案法の観測データ数を N とすると，典型的な
図 5 MC シ ミュレ ーション に よ る セル 廃 棄 率 推 定 値（ case 3 ）Fig. 5 Estimated values for cell loss rate by MC
simulation (case 3).

Table 12

表 12 実行時間の比較（ case 3 ）Comparison of computation time (case 3).

実行時間

提案法 従来法 [1]51 s
2,092 s

従来法 [2]
3,666 s

表 13 シミュレーション条件（ case 4 ∼ 6 ）Table 13 Condition of simulation (case 4 ∼ 6).
λ1
λ2
σ12
σ21

case 4
15.0
2.0
1.0
0.1

case 5
10.0
2.0
3.0
0.5

case 6
100.0
10.0
10.0
1.0

場合に T = 2N 程度であるが，毎回の E-step において複雑な被積分関数に対する数値積分の計算に時間がかかり，結果として表 14 の時間比となっている．状
態数 R に関する提案法の計算量は従来法 [1], [2] の計算量以下である．
また，設計誤差を見ると，6 例中 5 例において安全側であることがわかる．また，危険側となった case 2
でも図 4 からわかるように，従来法 [1], [2] は提案法よりも更に危険側の評価を与えている．また，表 14より，設計誤差の絶対値は状態数 R に依存しないと考えられる．

6. む

す び

EM アルゴ リズムによる MMPP パラメータ推定に
表 14 実行時間と設計誤差の比較
Table 14 Comparison of computation time anddesign error.

case 1case 2case 3case 4case 5case 6

時間比 1 時間比 2 設計誤差
1/7
1/115
10％
1/8
1/160
−7％
1/40
1/70
11％
1/11
1/85
7％
1/24
1/75
11％
1/7
1/93
15％

おいて，新たに MMGP モデルを提案し ，MMPP から MMGP へのモデル変換を行って EM アルゴ リズムを適用するパラメータ推定法を提案した．シミュレー
ションによる従来法 [1]，従来法 [2] との比較を行った
結果，提案法の実行速度は従来法よりも高速であり，また，セル廃棄率に対する推定精度は従来法と同程度であることを確認した．
セル廃棄率特性に関して，本論文で示した例において提案法の特性は真値に対して安全側の評価となる場

計誤差はセル廃棄率 10−3 におけるキュー長 q の提案合が多かった．提案法では，事象の到着時にのみ状態法による推定値と真値との誤差を示す．すなわち，提遷移が起こると考えているため，状態遷移が起こった

案法による推定値を q̂ とし真値を q0 とすると，設計時間の推定に誤差が生じ る．実験によると，その誤差

誤差 = (q̂ − q0 )/q0 × 100[％] である．

によって到着率が大きめの値として推定され，その結計算の実行時間比についてはどの場合も 7 分の 1 以

果，セル廃棄率特性では，真値のセル廃棄率特性より下となっている．提案法及び 従来法 [2] において最もも安全側の評価となることが多くなったと考えられる．

計算時間がかかる部分は forward-backward アルゴ リ

実行時間に関して，提案法は，従来法 [1] で必要と695

電子情報通信学会論文誌 2001/4 Vol. J84–B No. 4

する行列の固有値や数値積分等の複雑な計算を行う必山田雅裕 （ 正員）

要がなく，また従来法 [2] よりも少ない計算量，メモ平 10 長岡技科大・工・電気・電子シス

リ消費量で計算を行うことができるため，実行時間をテム卒．平 12 同大大学院修士課程了．同短縮できたものと考えられる．

年（ 株）エフ・デ ィー・シー入社．

文
[1]

献

T. Rydén, “An EM algorithm for estimation in
Markov modulated Poisson processes,” Computational Statistics & Data Analysis, vol.21, pp.431–447,
1996.

[2]

昭 55 東工大・理・数学卒．昭 60 同大大

Markov modulated Poisson processes via the EM

学院博士課程満期退学．昭 60 NTT 研究所入社．平 4 長岡技科大・工・助教授．待Syst., vol.1, pp.321–338, 1993.
H. Heﬀes and D.M. Lucantoni, “A Markov Modu-

ち行列理論，ネットワーク特性評価，大偏差理論等の研究に従事．理博．IEEE ，情

lated Characterization of Packetized Voice and Data

報理論とその応用学会，日本数学会各会員．

Traﬃc and Related Statistical Multiplexer Performance,” IEEE J. Selected Areas Commun., vol.sac-4,
no.6, pp.856–868, 1986.
[4]

A.P. Dempster, N.M. Laird, and D.B. Rubin, “Maximum Likelihood from Incomplete data via the EM
Algorithm,” J. Royal Statistical Society, vol. B39,
pp.1–22, 1977.

[5]赤穂昭太郎，“EM アルゴ リズムの幾何学，” 情報処理，vol.37, no.1, pp.43–51, 1996.

[6]

宮川雅巳，“EM アルゴ リズムとその周辺，” 応用統計学，[7]

M.A. Tanner, Tools for Statistical Inference; Meth-vol.16, no.1, pp.1–21, 1987.
ods for the Exploration of Posterior Distributions andLikelihood Function, Springer, 1996.
[8]

L.R. Rabiner, “A Tutorial on Hidden Markov Models and Selected Applications in Speech Recognition,”
Proc. IEEE, vol.77, no.2, pp.257–286, 1989.

[9]中川聖一，確率モデルによる音声認識，電子情報通信学会，
1988.

[10]

K.S.

Meier-Hellstern,“A

ﬁttingalgorithm

forMarkov-modulated Poisson processes having two arrival rates,” European J. Operational Research,
vol.29, pp.370–377, 1987.
[11]

W. Fischer and K. Meier-Hellstern, “The Markovmodulated Poisson process (MMPP) cookbook,” Performance Evaluation, vol.18, pp.149–171, 1992.

[12]

小 川 耕 司 ，中 川 健 治 ，“MMPP/D/1 キュー イン グ におけ る最適 IS シ ミュレ ーション 分布，” 信学論（ B-I ），vol.J80-B-I, no.2, pp.64–73, Feb. 1997.

[13]

山田雅裕，中川健治，“EM アルゴ リズムによる MMPP
パラメータ推定，” 情報理論とその応用シンポジウム講演
論文集，pp.889–892, Nov. 1999.

[14]

山田雅裕，“EM アルゴ リズムによる MMPP パラメータ推定について，” 長岡技科大大学院修士論文，2000.
（ 平成 12 年 3 月 10 日受付，9 月 22 日再受付）

696

健治 （ 正員）

L. Deng and J.W. Mark, “Parameter estimation foralgorithm with time descretization,” Telecommun.

[3]

中川