統計数理（2006）
第 54 巻 第 2 号 425–444
c 2006 統計数理研究所


特集「予測と発見」［原著論文］

ベイズ型時系列モデルによる成分分解を用いたネットワークトラヒックの定量的分析法
1

生駒 哲一 ・八名 和夫

2（受付 2006 年 1 月 16 日；改訂 2006 年 8 月 8 日）

要旨インターネットが広く普及した現在，ネットワークトラヒックを分析してこれをネットワー
ク管理に役立てることが重要となってきている．その分析の一方法として，ネットワークトラ
ヒックをトレンド，1 週間周期の変動，1 日周期の変動，および有色雑音成分に分解するベイズ
型時系列モデルを用いる方法を提案する．提案法で用いるモデルは，これらの成分ごとのシステム方程式と，各成分の和をネットワークトラヒックとする観測方程式とから成る線形ガウス型状態空間モデルで表現される．カルマンフィルタによる状態推定により各成分が求まる．モデルに含まれる分散や AR 係数のパラメータは，尤度の最大化により定める．またモデル間の比較や各成分の要不要の決定は，AIC（赤池情報量規準）に基づいて行う．提案法によるネット
ワークトラヒックの分析例として，法政大学小金井キャンパスのダイアルアップアクセス，および法政大学と SINET を結ぶネットワークに対して成分分解を行った結果について報告する．
キーワード： ネットワークトラヒック，成分分解，ベイズ型時系列モデル，状態空間
モデル，カルマンフィルタ．
1.

はじめに

現在，インターネットは広く普及し，我々の生活に欠かせないものとなりつつある．インター
ネットを支える技術も日々進歩し続けており，そこを流れるトラヒックも変化し多様になってきている．一方，企業や大学等の組織においては，組織内ネットワークを管理運営する必要がある．そのためには，トラヒックを解析してその特性の把握や将来の予測をし，適切な設備更新やあるいは動的に経路を振り分ける等の対応をとることになる．
図 1 に，ネットワークトラヒックの一例を示す．これは，法政大学小金井キャンパスに設置されたダイアルアップアクセスポイントにおいて，5 分ごとに下り（外部から法政大学へ）のト
ラヒック量（単位：bit/sec）を記録したものである．期間は 2001 年 5 月 28 日（月）から同年 7 月
8 日（日）の 6 週間である．これは前期授業が定常的に行われている期間である．また，この期間に祝日はない．この期間のうちの最初の 1 週間を図 2 に示す．図 2 は，2001 年 5 月 28 日の0 時から始まっている．
このトラヒックをみてまず判ることは，顕著な 1 日周期の変動があることである．この 1 日
周期の変動は，トラヒックが電話回線を経由したアクセスである事から，電話料金が格安にな
1
2

九州工業大学 工学部電気工学科：〒804–8550 福岡県北九州市戸畑区仙水町 1 番 1 号
法政大学 工学部電子情報学科：〒184–8584 東京都小金井市梶野町 3–7–2

426

統計数理 第 54 巻

図 1.

第 2 号 2006

法政大学小金井キャンパスのダイアルアップアクセストラヒック：2001 年 5 月 28 日（月）∼7 月 8 日（日）（6 週間）．

図 2. 法政大学小金井キャンパスのダイアルアップアクセストラヒックの第 1 週目．

る深夜にアクセスが集中し，昼間は電話料金が高い事や利用者の多くが大学に居てアクセスが少ない事等が，その要因として考えられる．また，日によって，アクセスの多い少ないの違いが見られる．全体的な傾向としても，1∼2 週目は多いトラヒックが，その後半ではやや少なめになってきているのも判る．そして最終週の中頃にて，トラヒックの一時的な上昇もみられる．
これらの変動傾向は，人間の視覚によって定性的に解ることである．これらを定量的に捉えるためには，何らかのモデルに基づく分析を行う必要がある．
ネットワークトラヒック時系列の統計的モデル化に関する研究は，Leland et al.（1994）, Paxon
and Floyd（1995）による報告以来数多くあり，大別すると，現象的なモデル化と構成的なモデル化に分類される．現象的なモデル化においては，トラヒック時系列を単一の時系列とみなし，その自己相似性，長相関性などが分析され，それらの特徴を有する FARIMA モデル（Resnick,1997）等が提案されている．しかしながら，本論文で取り上げるようにトラヒック時系列を特徴的なマクロな成分に分解して定量的に解析した例はない．構成的なモデル化（Willinger et al.,
1997; Yang and Petropulu, 2001）においては，トラヒックを要素となる成分に分解するが，パケット単位のミクロな成分分解であり，本論文で取り扱う成分分解とは異なる．
本研究では，上述のような特徴をもつネットワークトラヒックを解析の対象とする．このようなネットワークトラヒックを，トレンド，1 週間周期の変動，1 日周期の変動等の成分に分解する事を目的とする．そのために，時系列をトレンド，循環成分などに分解するベイズ型時系
列モデル（Kitagawa and Gersch, 1984）に基づく，新しい分析手法を提案する．
さて，Kitagawa and Gersch（1984）で提案されているベイズ型時系列モデルは，時系列デー
タをいくつかの特徴的な成分に分解することができる．成分としては，トレンド，季節性，曜日効果，AR 成分などがある．このモデルは経済データにおける季節調整を念頭において開発

ベイズ型時系列モデルによる成分分解を用いたネットワークトラヒックの定量的分析法

427

されたもので，典型的には月次データや日次データを扱い，季節性は 12ヶ月や 1 週間を周期とする成分を意図している．また曜日効果は，月次データの際に使われる成分で，各月に含まれる同一曜日の日数の違いを表すものである．
ネットワークトラヒックに対しては，単に Kitagawa and Gersch（1984）で挙げられている季節調整モデルをそのまま適用すれば良い訳ではない（高橋 他, 2001）．その理由は，一つにはトラヒックの持つ周期性が単一ではなく，1 日周期の変動と 1 週間周期の変動という，二つの周期性を持つことが挙げられる．またもう一つの理由としては，トレンド自体も，日内の短いスパ
ンでのトレンドと，数日に渡るような長期的なトレンドとがあり得るからである．さらに，これら以外の成分として，有色の雑音成分があると考えられる．よって，ネットワークトラヒッ
ク分析のためには，以上の成分を含むように Kitagawa and Gersch（1984）のモデルを拡張する必要がある．
関連する類似の研究としては，まず月次データに対する季節調整として，例えば電話料金収入の系列について，料金改定の効果項を入れたモデル（阿部・上田, 1985）や，さらに曜日効果も考慮したモデル（Yada et al., 1994）がある．これらは単一の周期（12ヶ月）成分とトレンド成分（月に関する階差）のみを扱っており，本研究で課題としている二つの周期性やトレンドは考慮していない．また，日次データの系列に関する成分分解としては，電話の通信トラヒックに対して，電話料金改定の効果を考慮した研究（矢田・北川, 1994）があり，そこでは料金改定をモデル構造の変化点として考え，観測ノイズ分散を大きくしてロバストな推定を試みたり，非ガウス型フィルタを用いた構造変化の推定を試みている．これら三つの研究のうち，最初の二つは電話料金収入の系列に対するモデル化である．三つめはトラヒック系列のモデル化ではあるが，電話料金改正に伴う収入の変化に主眼を置いたものである．また，データ系列のサンプ
リング間隔も月次や日次と長く，時々刻々変化するトラヒックを扱うには不十分である．
以上の事柄や背景をふまえて，本論文では，Kitagawa and Gersch（1984）で提案されているベイズ型時系列モデルに基づき，上述したネットワークトラヒックの性質を考慮した新しい分析手法を提案する．提案手法により，法政大学のダイアルアップアクセス（清水 他, 2003; Ikoma
et al., 2003; Shimizu et al., 2003）および SINET のトラヒック（今津 他, 2004）の分析を行う．得
られた結果を吟味し，トラヒックの性質について考察する．
2.

モデル

1 時間をサンプリングタイムとするトラヒックデータを想定する．トラヒック量（単位：kbit/sec）の対数をとったものを解析対象とする．対数を取る理由は，トラヒックデータは正の値を取る
ので，これを実数域に変換する為である．これらを満たすために，前処理として，まずトラヒッ
クの原系列（5 分サンプリング）12 期の平均を計算して 1 時間サンプリングの系列を得，その後対数をとるものとする．得られた系列を図 3 に示す．以下では，時刻（単位：hour）を整数 k で表し，解析に用いるトラヒックデータを y(k) と表すことにする．
2.1 トレンド成分
成分の階差がほぼ 0 となる成分を「トレンド成分」と呼ぶものとする．ただし階差は日を単位として取るものとする．これは，例えば階差が 1 の場合には，前日の値と本日の値との差がほぼ 0 となるものをトレンド成分とすることになる．以下では特に，トレンドの階差が 2 の場
合について説明する．他の階差についても，同様に定義することができる．
トレンド成分は，システム方程式(2.1)

xT (k) = F T (k)xT (k − 1) + GT (k)vT (k)428

統計数理 第 54 巻

第 2 号 2006

図 3. ダイアルアップアクセストラヒックの対数系列．

に従うものとして定義される．ここで vT (k) はシステムノイズで，正規分布 N (0, τT2 ) に独立に従うものとする．式（2.1）の各要素は，まず(2.2)

xT (k) = [ T1 (k) T2 (k) ]は，トレンド成分を時間順にその階差の日数分だけ並べた状態ベクトルである．次に行列 F T (k)は，時刻 k により 2 通りの値をとり，k = 0(mod 24) の時には



(2.3)

2
F T (k) =
1

−1
0



とする．ここで行列の 1 行目は階差 2 を取るための係数であり，2 行目以下は要素のシフトの役割を果たす．k = 0(mod 24) の時には F T (k) は単位行列とする．そしてベクトル GT (k) は，
時刻 k により 2 通りの値をとり，k = 0(mod 24) の時には
(2.4)



GT (k) = [ 1 0 ]であり，k = 0(mod 24) の時には GT (k) は零ベクトルとする．
行列 F T (k) および GT (k) を，上述のように時間変化するものとはせず，時刻 k によらず一定で式（2.3），（2.4）の値を取るものとすると，階差を時間（hour）としたトレンド成分を定義することができる．このトレンド成分の状態ベクトルを xt (k) = [ t1 (k) t2 (k) ] ，システムノイズ
を vt (k) ∼ N (0, τt2 )，上述の行列を F t ，Gt と表すことにする．必要に応じて，このような短いスパンのトレンドを用いることができる．
2.2 曜日成分
7 日を周期とする周期的変動成分を「曜日成分」と呼ぶものとする．曜日成分は，システム
方程式
(2.5)

xD (k) = F D (k)xD (k − 1) + GD (k)vD (k)2に従うものとして定義される．ここで vD (k) はシステムノイズで，正規分布 N (0, τD) に独立に従うものとする．式（2.5）の各要素は，まず(2.6)

xD (k) = [ D1 (k) D2 (k) · · · D6 (k) ]は，曜日成分を 6 日分並べた状態ベクトルである．次に行列 F D (k) は，時刻 k により 2 通り

ベイズ型時系列モデルによる成分分解を用いたネットワークトラヒックの定量的分析法の値をとり，k = 0(mod 24) の時には

(2.7)


−1 −1 · · ·

 1


1
F D (k) = 

..

.


−1

1

429











0

であり，k = 0(mod 24) の時には F D (k) は単位行列とする．そしてベクトル GD (k) は，時刻 k
により 2 通りの値をとり，k = 0(mod 24) の時には
(2.8)



GD (k) = [ 1 0 · · · 0 ]であり，k = 0(mod 24) の時には GD (k) は零ベクトルとする．
なお，祝日を含む週については，次のような扱いが可能である．祝日のトラヒック量を日曜のトラヒック量と同じと考えることにする．式（2.7）の行列は，曜日成分のベクトルに適用すると，曜日成分のベクトルが 1 日分シフトするという性質がある．よって，翌日が祝日の場合には，式（2.7）を複数回適用（行列のべき乗を用いる）し，曜日成分のベクトルの第 1 要素が日曜日になるようにすればよい．ただし，どの要素が日曜に対応するかを把握しておく必要がある．
2.3 日内変動成分
24 時間を周期とする周期的変動成分を「日内変動成分」と呼ぶものとする．日内変動成分
は，システム方程式(2.9)

xs (k) = F s xs (k − 1) + Gs vs (k)に従うものとして定義される．ここで vs (k) はシステムノイズで，正規分布 N (0, τs2 ) に独立に従うものとする．式（2.9）の各要素は，まず(2.10)

xs (k) = [ s1 (k) s2 (k) · · · s23 (k) ]は，日内変動成分を 23 時間分並べた状態ベクトルである．次に行列 F s は(2.11)

−1 −1 · · ·


 1

1
Fs =

..

.


−1

1











0

であり，ベクトル Gs は(2.12)



Gs = [ 1 0 · · · 0 ]とする．
2.4 AR 成分
有色の雑音成分を表すために，自己回帰過程に従う成分を考える．p 次の自己回帰過程に従う成分を「AR 成分」と呼ぶものとする．AR 成分は，システム方程式(2.13)xa (k) = F a xa (k − 1) + Ga va (k)430

統計数理 第 54 巻

第 2 号 2006に従うものとして定義される．ここで va (k) はシステムノイズで，正規分布 N (0, τa2 ) に独立に従うものとする．式（2.13）の各要素は，まずxa (k) = [ c1 (k) c2 (k) · · · cp (k) ](2.14)

は，AR 成分をその次数の分だけ並べた状態ベクトルである．次に行列 F a は，
a1 a2 · · ·


1


1
Fa =

..

.


(2.15)

ap

1











0

で，ここで a1 , a2 , . . . , ap は自己回帰係数である．そしてベクトル Ga は，


Ga = [ 1 0 · · · 0 ](2.16)とする．
2.5 観測方程式
各成分の和をネットワークトラヒックとする．これを観測方程式(2.17)

yk = T1 (k) + D1 (k) + s1 (k) + c1 (k) + wk

で表す．ここで wk は観測ノイズで，正規分布 N (0, σ 2 ) に独立に従うものとする．
2.6 状態空間表現
これまで定義してきたシステム方程式，観測方程式をまとめて，次の状態空間表現を得る．
xk = F k xk−1 + Gk v kyk = hxk+ wk(2.18)

ここで状態ベクトル xk は，各成分の状態ベクトル xT (k)，xD (k)，xs (k)，および xa (k) を並べたものであり，またシステムノイズベクトル v k は，各成分のシステムノイズ vT (k)，vD (k)，vs (k)，および va (k) を並べた列ベクトルである．状態遷移行列 F k は，各成分の状態遷移行列
F T (k)，F D (k)，F s ，および F a を対角要素に持つブロック対角行列とする．また行列 Gk は，各成分の列ベクトル GT (k)，GD (k)，Gs ，および Ga を対角要素に持つブロック対角行列である．行ベクトル h は，式（2.17）が成り立つように定義するものとする．
このように状態空間表現されたモデルを，状態空間モデルと呼ぶ．式（2.18）の状態空間モデ
ルは，全てのノイズが正規分布に従い，また線形式であるので，線形ガウス型状態空間モデルのクラスに属する．
3.

モデルの推定3.1 状態推定
状態推定とは，状態空間モデルにおいて，観測の系列 Yk = (y1 , y2 , . . . , yk ) が与えられた下での状態の条件付分布を求めることである．状態推定は，与えられる観測系列の時刻（表記 Yk のk）と推定すべき状態の時刻との関係に応じて，
「予測」，
「ろ波」，
「平滑化」の 3 つに分類することができる．
ろ波とは，観測系列の時刻と状態の時刻が等しい場合をいう．例えば，現在までの全ての観測データが与えられた下で，現在の状態を推定することは，ろ波である．予測とは，状態の時

ベイズ型時系列モデルによる成分分解を用いたネットワークトラヒックの定量的分析法

431

刻が，観測系列の時刻よりも未来になっている場合である．平滑化とは，予測の逆の状況で，状態の時刻が，観測系列の時刻よりも過去になっている場合である．つまり平滑化は，現在までの観測データを使って，過去の状態を推定していることに相当する．求める条件付分布は，
時刻 k の状態の推定に関しては，1 期先予測が p(xk |Yk−1 )，ろ波が p(xk |Yk )，そしてラグ L > 0の平滑化が p(xk |Yk+L ) となる．
本論文で用いるモデルは，式（2.18）のように，線形ガウス型状態空間モデルである．線形ガ
ウス型状態空間モデルは，状態の分布がガウス分布になるので，その平均ベクトルと共分散行列を求めるカルマンフィルタにより，状態推定が可能である．また平滑化についても，カルマ
ンフィルタと同様に，平均ベクトルと共分散行列を求めるアルゴリズムが提案されている．以
下にこれらのアルゴリズムを簡潔に示す．
3.2 カルマンフィルタ
状態の初期分布を適切に設定し，これを p(x0 ) と表すものとする．初期分布 p(x0 ) はガウス
分布で，その平均ベクトルが x0 ，共分散行列が V 0 であるものとする．
状態推定のための平均ベクトル，共分散行列を次のように初期化する．
(3.1)

x0|0 = x0
V 0|0 = V 0

ここで左辺の添え字は，縦線の前が状態の時刻，後が観測系列の時刻を表す．
時刻を k = 1 から開始し，これを 1 時刻ずつ進めながら，以下で示す 1 期先予測とろ波の手
続きを交互に繰り返す．
1 期先予測(3.2)

xk|k−1 = F k xk−1|k−1
V k|k−1 = F k V k−1|k−1 F k + Gk QGk

ここで正方行列 Q は，システムノイズベクトル v k の共分散行列である．
観測の 1 期先予測 p(yk |Yk−1 ) もガウス分布に従い，その平均 ȳk|k−1 と分散 Σk|k−1 は(3.3)

ȳk|k−1 = hxk|k−1
Σk|k−1 = hV k|k−1 h + σ 2

となる．
ろ波
(3.4)

xk|k = xk|k−1 + kk ek
V k|k = (I − kk h) V k|k−1

ここで，I は単位行列，ek は観測の 1 期先予測（平均値）の誤差(3.5)

ek = yk − ȳk|k−1

であり，また kk はカルマンゲインと呼ばれ，(3.6)と計算される．



kk = V k|k−1 h Σk|k−1

−1

432

統計数理 第 54 巻

第 2 号 2006

3.3 平滑化
本論文では，固定区間平滑化と呼ばれる平滑化法を用いる．この方法は，観測データ系列を固定して，全ての時刻の状態の平滑化推定を求める．そのためには，観測データ系列の全ての時刻について，上で述べたろ波と 1 期先予測を予め計算し記憶しておく必要がある．
全観測データ系列を YN = (y1 , y2 , . . . , yN ) とする．固定区間平滑化は，データ系列の末尾から開始し，データ系列の先頭方向に 1 時刻ずつ遡りながら，次の手続きを繰り返し適用する．
(3.7)

xk|N = xk|k + Ak xk+1|N − xk+1|k

V k|N = V k|k + Ak V k+1|N − V k+1|k Ak

ここで行列 Ak は，次の通り計算する．


Ak = V k|k F k+1 V −1
k+1|k

(3.8)

3.4 パラメータの推定
上述のモデルで状態推定を行う際には，観測ノイズの分散 σ 2 とシステムノイズの分散 τ 2 =22
τT , τD, τs2 , τa2 ，及び AR 係数 a = (a1 , a2 , . . . , ap ) を与える必要がある．これらのパラメータをまとめて θ と表すものとする．
パラメータ θ は，尤度を最大化することにより求める．最大化には，数値的最適化手法（グリッドサーチ，準ニュートン法等）を用いる．尤度関数の対数を取ると，(3.9)l(θ) ≡ log p(YN ; θ)Nlog p(yk |Yk−1 ; θ)

=
k=1

=−
−

1
N
log 2π −22Nlog Σk|k−1
k=1

N
yk − ȳk|k−1
1
2 k=1
Σk|k−12となる．最後の式の各要素は，カルマンフィルタの 1 期先予測の過程（式（3.3））で得られる．
ここで扱っているように観測が 1 変量の場合には，観測ノイズの分散を 1 としたモデルの尤度を最大化することで，パラメータを推定することができる（北川, 1993）．具体的には，観測ノイズの分散を 1 としたモデルを想定し，これにシステムノイズ分散 τ̃ 2 を与えカルマンフィ
ルタを適用する．その結果から，元のモデルの観測ノイズ分散の推定値 σ̂ 2 及びシステムノイ
ズ分散 τ 2 = σ̂ 2 τ̃ 2 と，対数尤度関数値とを得ることができる．この手順では，AR 係数 a の扱
いについては，変更は生じない．
観測ノイズの分散を 1 としたモデルに基づく，観測の 1 期先予測の平均と分散を，それぞれȳ˜k|k−1 ，Σ̃k|k−1 と表記する．ここで，ȳk|k−1 = ȳ˜k|k−1 および Σk|k−1 = σ 2 Σ̃k|k−1 となることが導かれる（σ 2 は元のモデルの観測ノイズ分散である）．これらを式（3.9）に代入し，偏微分して0 とおくと，元のモデルの観測ノイズ分散 σ 2 の推定値
(3.10)1σ̂ 2 =
NN

yk − ȳ˜k|k−1

k=1

Σ̃k|k−12を得る．これを再び式（3.9）に代入すると，対数尤度関数は(3.11)l(τ̃ 2 , a) = −N
log 2πσ̂ 22

ベイズ型時系列モデルによる成分分解を用いたネットワークトラヒックの定量的分析法

433

N

−

1
N
log Σ̃k|k−1 −2 k=12となる．これを最大化して，パラメータ τ̂ 2 ，â を求める．
3.5 モデル選択
異なるモデルが存在してそれらの比較をする場合や，モデルの成分の要不要を決定する場
合，もしくはモデル中のトレンド階差や AR 次数を決定する場合には，AIC（赤池情報量規準）（Akaike, 1974）(3.12)AIC = −2l(θ̂) + 2 × (θの自由度)を用いる．候補モデル中で AIC の値の最も小さいモデルを，最も良いモデルとして選択する．
AIC は，パラメータの自由度が比較するモデル間で異なる事を考慮した規準である．これは，解析に用いるモデルと真のモデルとの分布（観測値の予測分布）の相違をカルバック ライブラー
情報量により計るものとし，適切な仮定の下でその近似値を与えるものである．
4.

トラヒックデータの解析

ネットワークトラヒックの解析結果について報告する．解析は，提案法で示したモデルの他に，成分の適切な組み合わせから成るモデルについても行った．成分としては，トレンド成分（日についての階差 T ，および，時間についての階差 t），曜日成分 D，日内変動成分 s，そして AR 成分 c を用いた．トレンド成分の階差は 1 または 2 とし，AR 成分の次数は 1 または 2とした．トレンド成分の階差が異なる場合や，AR 成分の次数が異なる場合は，別のモデルとみなすものとする．
パラメータの推定では，分散パラメータの数値的最適化法として，2 段階のグリッドサーチを用いた．1 段目では，分散パラメータの指数部分を探索し，2 段目では 1 桁の有効数字を探索した．また AR 係数の推定もグリッドサーチで行い，偏自己相関係数（PARCOR）を用いて，PARCOR を (−1, 1) の範囲で 0.05 刻みでグリッドを構成し探索した．分散パラメータのグリッ
ドサーチと PARCOR のグリッドサーチは個別に交互に行い，両グリッドにて最適解が見つかり収束するまで探索を行う方法を採った．ただし分散パラメータには最小値（ここでは 10−10 ）を設定し，これより小さな値については探索しないものとした．
初期分布の平均 x0 と共分散行列 V 0 については，その適切な値の求め方はいくつか知られている（1．パラメータとして扱い，尤度最大化により求める．2．適当な初期分布にてフィル
タリングと平滑化を一旦行い，その平滑化推定値を初期分布に用いる．等）．ここでは簡便な方法として，次の通り値を与えることにした．まず，成分は互いに独立と仮定した．次に各成分内では，成分を構成する状態ベクトルの要素ごと独立に正規分布に従うものとした．各成分で用いた正規分布は次の通りである．トレンドについては，平均をデータ系列の先頭 24 時刻（すなわち 24 時間＝ 1 日）の算術平均とし，分散を 1 とした．曜日成分については平均 0 分散
1 とし，AR 成分については平均 0 分散 1 とした．日内変動成分については，尤度に対して比較的大きな影響を与えることが実験的に分かったので，それを考慮し，平均 0 分散 τs20 と設定
し，τs20 としては適宜 10 や 102 を与えるものとした．
以下，法政大学小金井キャンパスのダイアルアップアクセストラヒックと，法政大学と SINET
を結ぶ SINET トップトラヒックについて，それぞれの解析結果を述べる．
4.1 ダイアルアップアクセストラヒック
図 1（第 1 週目は図 2）の小金井キャンパスのダイアルアップアクセスのトラヒック（5 分間隔

434

統計数理 第 54 巻

図 4.

表 1.

表 2.

第 2 号 2006

ダイアルアップアクセストラヒック対数系列の第 1 週目．

ダイアルアップアクセストラヒックの最大対数尤度 l(θ̂) と AIC．

ダイアルアップアクセストラヒック：T（1）+D+s+c（2）モデルにて推定されたパラメータ．

で計測）について，12 期の平均を計算しその対数を取った系列が図 3 である．この系列を解析対象とし（その第 1 週目を図 4 に示す），提案法のモデル及びその他のモデルを適用した．ここで，日内変動の初期分散は τs20 = 10 とした．
各モデルの最大対数尤度 l(θ̂) と AIC を，表 1 に示す．表中のモデル名の箇所にて，括弧内の数字は，トレンドの階差，および AR 次数を表す．表 1 を見ると，提案法のモデル族が他のモデルと比べて良い AIC 値であることが判る．特に，AR 成分を持つモデルと持たないモデルとで，AIC 値に大きな違いがあることが見てとれる．これは，トラヒックが有色雑音成分を含

ベイズ型時系列モデルによる成分分解を用いたネットワークトラヒックの定量的分析法

図 5.

ダイアルアップアクセストラヒック：トレンド成分の推定結果（平滑化の平均：実線）と解析系列（破線）．

図 6.

図 7.

ダイアルアップアクセストラヒック：曜日成分の推定結果（平滑化の平均）．

ダイアルアップアクセストラヒック：日内変動成分の推定結果（平滑化の平均）．

図 8.

ダイアルアップアクセストラヒック：AR 成分の推定結果（平滑化の平均）．

435

436

統計数理 第 54 巻

第 2 号 2006

図 9.

ダイアルアップアクセストラヒック：推定の結果（平滑化の平均），得られた残差．

図 10.

ダイアルアップアクセストラヒック：推定の結果（平滑化の平均）得られた残差のヒス
トグラム．

むことを示唆していると考えられる．提案法のモデル族のうち，トレンド成分（日についての階差 ）T の階差が 1 のものが良く，そのうち AR 次数が 2 のモデルが微小ながら良い結果であることが判る．しかしこの微小な差は有意ではないと考えられるので，状況によっては AR 次
数が 1 のモデルを選択することもあり得る．
この最も良いモデル T（1）+D+s+c（2）にて得られたパラメータ値を表 2 に，推定された成分（平均）を図 5，6，7，および 8 に示す．また，推定された成分（平均）の和と，トラヒック系
列との差（残差）を，図 9 に示す．
以下，推定結果について考察する．まず，推定されたトレンド成分（図 5）は，ほぼ一定値となっている．これは，トレンド項システムノイズ分散 τT2 が，探索グリッドの最小値 10−10 となっていることと符合する．すなわち，本トラヒックデータに関しては，トレンド成分は確定的な（確率的要素を含まない）一定値とみなせることがわかる．
次に曜日成分（図 6）についてみると，大半の場合，木曜日が最もトラヒックが多く，次にト
ラヒックが多いのは日曜日である．水曜土曜はその次，月曜火曜がそれに続くトラヒック量を持ち，金曜日が最も少ないことがわかる．
日内変動成分（図 7）についてみると，深夜 0 時の少し前から徐々にトラヒック量が増え始め，
その後減ってゆくという変動パターンが得られている．これは，ダイアルアップアクセスが 23時以降の電話料金が格安になる時間帯に集中していることを，成分として抽出したものと考えられる．
そして，AR 成分（図 8）をみると，他の 3 つの成分では表しきれなかった微細な変動を，成ベイズ型時系列モデルによる成分分解を用いたネットワークトラヒックの定量的分析法

437

分として抽出していることがわかる．特に，最終週中盤の一時的なトラヒック上昇において，その変動が顕著に抽出されている．提案法のモデルでは，推定された AR 成分が，数時間程度の短期的なトレンドを表すものと考えることができる．
最後に，残差（図 9）が正規分布に従うかどうかを調べた．残差のヒストグラムを図 10 に示す．図 10 をみると，やや左右非対称な単峰の分布であることが判る．残差の平均は −0.0056，
標準偏差は σ = 0.1468 である．−4σ の以下の範囲にも残差値が 7 点存在している．これら 7 点
を外れ値とみなしてカイ二乗適合度検定を行ったところ，値 Q = 65.10 が得られた．ヒストグ
ラムの階級数 10，正規分布のパラメータ数 2 より，カイ二乗分布の自由度は (10 − 1) − 2 = 7である．自由度 7 のカイ二乗分布の上側 5％点は 14.067，上側 1％点は 18.475 であることから，得られた Q の値では残差が正規分布に従うとは言えない．
残差が正規分布に従わない原因を考察してみると，提案法で用いたモデルには，トレンド成
分（T 成分），曜日成分（D 成分）といった，日付変更時に値が切り換わる成分を持ち，これが原因の一つではないかと考えられる．すなわち，日付変更時にこれら成分の値が急に変化する為，
AR 成分が日付変更時にレベルシフトを伴うこととなり，また残差成分も日付変更時に大きな値を取ることとなる．実際，図 8 の AR 成分見ると，日付変更時にレベルシフトが生じており，
図 9 の残差成分を見ると，日付変更時に大きな値を取っていることがわかる．別の原因としては，線形ガウスモデルを適用する為に原データに適用している対数変換が考えられる．対数変換ではなく，別の変換，例えば Box-Cox 変換（Box and Tiao, 1973）などを用いることで，これを改善できる可能性がある．
4.2 SINET トップトラヒック
図 11 に，法政大学と SINET とを結ぶネットワークのトラヒックを示す．これもダイアル
アップアクセストラヒック同様，5 分ごとに下り（外部から法政大学へ）のトラヒック量（単位：bit/sec）を記録したものであり，期間は 2001 年 5 月 28 日（月）0 時から同年 7 月 8 日（日）の 6週間である．この期間のうちの最初の 1 週間を図 12 に示す．
乱れの多いデータ系列ではあるが，顕著な 1 日周期の変動と，1 週間周期の変動パターンを見ることができる．この 1 日周期の変動は，利用者の多い昼間にアクセスが集中し，夜間は少なくなっている．また 1 週間の変動パターンについては，日曜日が顕著に少ないトラヒックであることがわかる．全体的な傾向としても，2∼3 週目でトラヒックが一旦少なめになり，その後回復してやや多くなっているのが判る．
この SINET トップのトラヒック系列を，ダイアルアップアクセストラヒック同様，12 期の平均を計算した後対数を取り，得られた系列を解析の対象とする．解析対象の系列と，その最初の 1 週間を，それぞれ図 13，図 14 に示す．
この解析系列に，提案法で示したモデルと，適切な成分の組み合わせから成るモデルを適用した．ここで，日内変動の初期分散としては，T+D+s+c モデル族（提案法のモデル族）については τs20 = 102 とし，その他のモデルについては τs20 = 10 とした．
得られた最大対数尤度と AIC の値を，モデルごとにそれぞれ，表 3 に示す．表を見ると，提案法のモデル族が良い AIC 値を得ている．特に，ダイアルアップアクセストラヒックと同様に，AR 成分を持つモデルと持たないモデルとで AIC 値に大きな違いがあり，トラヒックが有色雑音成分を含むことを示唆していると考えられる．提案法のモデル族の中でも，トレンド階差 1 のモデルが良い AIC 値となっていて，そのうち，AR 次数が 2 のモデル T（1）+D+s+c（2）
が最も良い AIC 値を得ている．最良の T（1）+D+s+c（2）モデルについて，2 段階のグリッド

438

統計数理 第 54 巻

図 11.

第 2 号 2006

法政大学 SINET トップトラヒック：2001 年 5 月 28 日（月）∼7 月 8 日（日）（6 週間）．

図 12.

法政大学 SINET トップトラヒックの第 1 週目．

図 13.

図 14.

SINET トップトラヒックの対数系列．

SINET トップトラヒック対数系列の第 1 週目．

ベイズ型時系列モデルによる成分分解を用いたネットワークトラヒックの定量的分析法
表 3.

439

SINET トップトラヒックの最大対数尤度 l(θ̂) と AIC．

表 4. SINET トップトラヒック：T（1）+D+s+c（2）モデルにて推定されたパラメータ．

図 15.

SINET トップトラヒック：トレンド成分の推定結果（平滑化の平均）と解析系列．

サーチにより得られたパラメータの値を表 4 に示す．
T（1）+D+s+c（2）モデルについて，推定された各成分（平均）を図 15，16，17，および 18 に示す．また，推定された成分の和と，トラヒック系列との残差を図 19 に示す．
以下，結果について考察する．まず，推定されたトレンド成分（図 15）は，ややなだらかではあるが，2∼3 週目でやや減少し，その後増加に転じるという，トラヒック系列の長いスパンでの傾向を抽出している．次に曜日成分（図 16）は，日曜日が最もトラヒックが少なく，続いて土曜日が少なくなっていることがわかる．その他の曜日は週によって異なるが，週の前半で

440

統計数理 第 54 巻

図 16.

図 17.

SINET トップトラヒック：曜日成分の推定結果（平滑化の平均）．

SINET トップトラヒック：日内変動成分の推定結果（平滑化の平均）．

図 18.

図 19.

第 2 号 2006

SINET トップトラヒック：AR 成分の推定結果（平滑化の平均）．

SINET トップトラヒック：推定の結果（平滑化の平均），得られた残差．

ベイズ型時系列モデルによる成分分解を用いたネットワークトラヒックの定量的分析法

図 20.

441

SINET トップトラヒック：推定の結果（平滑化の平均）得られた残差のヒストグラム．

はトラヒックが多く，後半では少ない傾向が見てとれる．日内変動成分（図 17）についてみると，昼間トラヒック量が多く夜少ないという変動パターンが抽出されている．そして，AR 成分（図 18）をみると，微細な変動を成分として抽出しているが，特に日曜日のトラヒック量減少をこの成分でも表していることがわかる．
残差（図 19）が正規分布に従うかどうか調べた．残差のヒストグラムを図 20 に示す．やや左右非対称な単峰の分布をしている．残差の平均は 3.7837 × 10−4 ，標準偏差は σ = 0.0139 である．±4σ の範囲外に残差値が 4 点存在しており，これら 4 点を外れ値とみなして，カイ二乗適合度検定を行った．その結果，階級数 10 に対して Q = 53.92 が得られ，正規分布に従うとは言えないという結論となった．
5.

おわりに

Kitagawa and Gersch（1984）で提案されているベイズ型時系列モデルに基づき，ネットワー
クトラヒックの系列をトレンド，1 週間周期の変動，1 日周期の変動，および有色雑音成分に分解する新しい分析手法を提案した．提案法を，法政大学のダイアルアップアクセス，および
SINET トップのトラヒックに適用し，成分分解の結果を得，結果について若干の考察を行った．
トラヒックへの適用に際しては，提案法で示したモデルの持つ成分を適切に組み合わせをした
モデルをいくつか挙げ，それらモデルと提案法のモデルとを AIC により比較した．比較の結果，提案法のモデルが他のモデルよりも良いことが判った．
本手法は，ネットワークトラヒックの定量的な成分分解を可能とするものである．これには
次のような意義がある．まず，分析手法の妥当性が，情報量規準（AIC）の意味において保証されていることが挙げられる．各種の分析手法が工夫され提案されてはいるが，どれが妥当な手法であるか決めるのは容易なことではない．本論文で提案するような，統計モデルの枠組みでの分析手法では，情報量規準による手法の比較が可能であり，その意味で妥当な分析法を採用できることになる．本論文では，成分の組み合わせに関して，この妥当性を得ている．次に，成分分解が自動的かつ定量的に行えることが，利点として挙げられる．扱っているトラヒック
データは，目視であっても，1 日の周期的変動や 1 週間の周期的変動傾向が見て取れるものではある．しかし，それは人間の感覚に基づくものであって，数値化するのには別の客観的合理的方法が必要である．また，人手による手作業でそのような成分分解を行うのは現実的ではない．提案手法はこれを自動的に行えるので，ネットワークトラヒックに対する予測や制御を自442

統計数理 第 54 巻

第 2 号 2006

律的自動的に行うシステムの開発に利用できるというメリットがある．そして，成分分解を継続的に実施することで，長期にわたる定量的データが蓄積でき，ネットワークの更新設計などに際し有益な指針となり得る点も強調しておく．
今後の課題としては，以下の項目が挙げられる；
1）提案法で示したモデルは，カルマンフィルタを適用可能とする為，線形ガウス型状態空間モデルの枠内で，モデルが立てられている．しかし，実際のトラヒックの性質は，必ずしも線形性やガウス性を満たすとは限らない．トラヒックの性質を考慮した非線形項や，非ガウスなノイズ項の導入により，より良いモデルが開発できる可能性がある．その際，状態推定には，
最近盛んに研究されている非ガウス型フィルタ（Doucet et al., 2001）を用いることになる．
2）本論文で行ったネットワークトラヒックの分析では，残差が白色ガウス雑音にはならな
かった．その原因としては，トラヒック原系列に対して行う変換や，日付変更時のレベルシフトなどが考えられることを，4 章にて考察した．これを改善するため，今後の課題として，まずは，トラヒック原系列に対して行う変換としてより適切なもの（例えば Box-Cox 変換）を用いることが挙げられる．次に，日付変更時のレベルシフトに対応できるモデルを採用することが考えられる．現時点では，カルマンフィルタの枠組みで推定を行う為，モデルを線形ガウスに限定していたが，たとえば裾の重い分布を持つ非ガウスモデルに変更することで，値の急激な変化にも対応可能なモデルが構築できる可能性がある．
3）本論文では，トラヒック系列を 1 時間サンプリングタイムの系列に変換した上で，ベイズ型時系列モデルを適用している．その理由は，提案手法のモデルの枠組みでは，あまり長い周期の日内変動を記述しづらい事による．より具体的には，モデル中の日内変動成分の周期 24時間を表すために状態ベクトルのうち 23 次元を使用しているが，仮にサンプリングタイムが原系列のまま 5 分であるとすると，日内変動を表すのに（24 × 12 − 1）次元のベクトルが状態に現れることになる．このように高次元の状態ベクトルは，状態推定に困難を生じると考えられる．この問題点を解消するアイデアが望まれる．
4）トラヒックの短期予測が可能であれば，予測に基づいたルーティングの逐次修正などを可能としトラヒック制御の立場から興味深い．そのためには，2）で述べたように，原系列のサンプリングタイム 5 分のままでモデル化を行うのが望ましい．

謝辞本研究は，統計数理研究所共同研究プログラム（16–共研–0010）による．トラヒックデータの整理と予備的な解析は，法政大学大学院の高橋兼基君（平成 14 年修士課程修了．現在，日本テ
レコム株式会社勤務），清水威臣君（平成 16 年修士課程修了．現在，株式会社電通勤務），今津
和之君（平成 17 年修士課程修了．現在，日本電気株式会社勤務）による．

参 考 文 献阿部威郎，上田 徹（1985）. 状態空間モデルを用いた電話収入予測，電子通信学会論文誌（A）, J68-A（5）,
437–443.
Akaike, H.（1974）. A new look at the statistical model identiﬁcation, IEEE Transactions on Automatic Control, AC-19, 716–723.
Box, G. E. P. and Tiao, G. C.（1973）. Bayesian Inference in Statistical Analysis, Addison-Wesley.
（2001）. Sequential Monte Carlo Methods inDoucet, A., de Freitas, N. and Gordon, N. J.（eds.）
Practice, Springer, New York.

ベイズ型時系列モデルによる成分分解を用いたネットワークトラヒックの定量的分析法

443

Ikoma, N., Shimizu, T., Imazu, K. and Yana, K.（2003）. Decomposition of network traﬃc of dialup
access into trend and periodic components based on nonstationary Bayesian model, Proceedingsof SICE Annual Conference 2003, Fukui, Japan, 359–362.

今津和之，清水威臣，八名和夫，生駒哲一（2004）. トレンド，周期変動および AR の特徴を抽出する状
態空間モデルによる SINET トラヒック系列の成分分解，第 20 回ファジィシステムシンポジウ
ム講演論文集，545–548.
北川源四郎（1993）. 『FORTRAN77 時系列解析プログラミング』, 214–217, 岩波書店，東京．
Kitagawa, G. and Gersch, W.（1984）. A smoothness prior state space approach to the modelingof time series with trend and seasonality, Journal of the American Statistical Association,（386）, 378–389.
79
Leland, W. E., Taqqu, M. S., Willinger, W. and Wilson, D. V.（1994）. On the self-similar nature of
Ethernet traﬃc, IEEE/ACM Transactions on Networking, 2, 103–115.
Paxon, V. and Floyd, S.（1995）. Wide-area traﬃc: The failure of Poisson modeling, IEEE/ACM
Transactions on Networking, 3, 226–244.
（5）, 1805–1848.
Resnick, S. I.（1997）. Heavy tail modeling and teletraﬁc data, Annals of Statistics, 25
清水威臣，今津和之，生駒哲一，八名和夫（2003）. 非定常べイズモデルによるネットワークトラヒックの成分分解，2003 年春季電子情報通信学会総合大会講演論文集, A–4–47, 136–136.
Shimizu, T., Imazu, K., Ikoma, N. and Yana, K.（2003）. Nonstationary Bayesian time series analysis
for network traﬃc by decomposition model into trend, periodic and AR components, Proceedings of Science of Modeling
The 30th Anniversary of the Information Criterion（AIC） ,Yokohama, Japan, Dec. 14–17, 408–409.
高橋兼基，生駒哲一，八名和夫（2001）. 非定常ベイズモデルによるネットワークトラヒックの成分分
解，第 17 回ファジィシステムシンポジウム講演論文集, 153–156.
Willinger, W., Taqqu, M. S., Scherman, R. and Wilson, D. V.（1997）. Self-similarity through highvariability: Statistical analysis of Ethernets LAN traﬃc at the source level, IEEE/ACM Transactions on Networking, 5, 71–86.
矢田 健，北川源四郎（1994）. 競合環境下における通信トラヒックデータの分析，第 62 回日本統計学会講演報告集, 34–35.
Yada, T., Kitawaga, G. and Inoue, M.（1994）. Prediction of telephone revenue by using a Kalmanﬁlter, Proceedings of the 14th International Teletraﬃc Congress（ITC）, 14, 1165–1174.
Yang, X. and Petropulu, A. P.（2001）. The extended alternating fractal renewal process for modelingtraﬃc in high-speed communication networks, IEEE Transactions on Signal Processing, 49,
1349–1363.

444

Proceedings of the Institute of Statistical Mathematics Vol. 54, No. 2, 425–444(2006)

Quantitative Analysis Method of Network Traﬃc
by Component Decomposition Using Bayesian Time Series Model
Norikazu Ikoma1 and Kazuo Yana2
1 Faculty of Engineering, Kyushu Institute of Technology
2 Faculty of Engineering, Hosei University

The ever-increasing expansion of the Internet constantly increases the importance ofanalyzing network traﬃc and making use of those results. To analyze network traﬃc inthis study, we used a Bayesian time series model, which decomposes network traﬃc into
trends, weekly and daily cycles, and colored noise components. The model is represented ina linear-Gaussian state space model that consists of system equations for each component
and an observation equation, in which sums of the components yields the observed network
traﬃc. State estimation by Kalman ﬁlter is used to obtain each component. Variance andAR parameters in the model are obtained to maximize the likelihood. Then a comparison
among candidate models and selection of the components proceed based on the Akaike
Information Criterion (AIC). To provide an example of an analysis using this proposed
method, we report on the decomposition results of traﬃc data of dialup access at Koganei
Campus and those at the top network domain of Hosei University, Tokyo which connects
to the SINET.

Key words: Network traﬃc, decomposition, Bayesian time series model, state space model, Kalmanﬁlter.