論 文
ATM バーストセルトラヒックの待ち時間特性に対する
IS シミュレーション
上田 健志†

中川 健治††IS Simulation for the Waiting Time of ATM Burst Cell Traﬃc
Takeshi UEDA† and Kenji NAKAGAWA††

あらまし ATM 交換機内のバッファサイズを設計するために，様々なトラヒックに対するセル遅延特性を調べる必要がある．本研究では，セル遅延特性，特にセル遅延の α 値をシミュレーションによって求めることを目的とする．しかし ，従来法であるモンテカルロ（ MC： Monte Carlo ）法では計算時間がかかること，及び疑似乱数の精度の限界のため，非常に小さい α に対するセル遅延特性を求めることができない．そこで，本研究ではインポータンスサンプ リング（ IS：Importance Sampling ）法を適用し ，高速に精度良く 10−9 程度以下のセル遅延特性を得ることを目的とする．特に，ランダムセルトラヒックに対する待ち行列モデルである Geo/D/1 や
バーストセルトラヒックに対する 2 状態 MMGP/D/1 待ち行列の定常状態における待ち時間特性に IS 法を適用し ，large deviation 理論に基づいて近似的な最適シミュレーション分布を得る．そして得られたシミュレーショ
ン分布によって実際に IS シミュレーションを行い，MC シミュレーションと比較し高速性と 95％信頼区間の点で優れていることを示す．更に理論値，近似値との比較によって IS 推定値の精度の良さを示す．
キーワード

シミュレーション，IS 法，Geo/D/1，2 状態 MMGP/D/1

1. ま え が きProcess ），及び バースト セルト ラヒックのモデ ルである 2 状態マルコフ変調幾何過程（ MMGP：Markov

ATM 網ではセルの 統計多重のため 各ノード でセModulated Geometric Process ）を 入力とし ，一定

ルがバッファに 蓄積され ，そこで 遅延が 生じ る．一

サービ ス時間，単一サーバの待ち行列システム，すな

つのノード で許容され る遅延 W の α 値，すなわち

わち，Geo/D/1 及び 2 状態 MMGP/D/1 に対して IS

P (W > w) = α となる w の値を調べることが重要となる．しかし ，α が小さいと一般的なモンテカルロ（ MC：Monte Carlo ）シミュレーション法ではセル遅度良く得ることを目的とする．

法を適用し ，10−9 オーダの待ち時間特性を高速に精これらのモデル Geo/D/1 と 2 状態 MMGP/D/1 に延特性を得ることは非常に困難である．そこで本研究対する待ち時間解析は，他の方法，例えば行列解析等

では，非常に小さい α に対する遅延特性を得るための理論的な方法によっても可能であるが，本研究はそに ，インポータン スサンプ リング（ IS： Importanceのような理論的解析が困難なより複雑なトラヒックモ

Sampling ）法を適用する．IS 法は非常に小さい確率

デルへの IS 法の拡張のための基礎的な検討である．

を高速に求めるために開発されたシミュレーション手法である．
本研究では ，ATM 網に おけ るランダ ムセルト ラ
ヒックの モデ ル で あ る 幾 何 過 程（ Geo： Geometric† NEC 航空宇宙システム株式会社，横浜市

2. IS シミュレーションの概要
イン ポ ー タン ス サン プ リン グ（ IS： ImportanceSampling ）シ ミュレ ーション 法は ，通常のモンテカ
ルロ（ MC： Monte Carlo ）シミュレーション法では

NEC Aerospace Systems Limited, Yokohama-shi, 222–0033

得られ な い 小さい 確 率を 高 速に 精 度良 く 得るため

Japanに 開発された方法で ，待ち行列におけ るブ ロック確†† 長岡技術科学大学，長岡市
Department of Electrical Engineering, Nagaoka Universityof Technology, Nagaoka-shi, 940–2188 Japan

1782

率 [2], [3], [10] や通信路符号の誤り率 [6] の推定など 多くの工学的問題に対して適用されている．

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

pp. 1782–1790 1999 年 10 月

論文／ ATM バーストセルトラヒックの待ち時間特性に対する IS シミュレーション

MC 法ではある確率 p を推定する場合に，推定値の分散を同程度にすれば IS 法は MC 法よりシミュレー

標準偏差を 0.1p 以下にするために必要なシミュレー

ション回数が少なくてすむ．これが IS 法が MC 法よション回数は 100/p 以上である．したがって，p が小り高速な理由である．

さいときには時間的な限界のため MC 法は実行不可能となる．また，文献 [5] でも示されているように，正すそ

規分布の 裾の確率を推定するというような簡単な場合
でも疑似乱数の精度の限界のため MC 法では小さい確3. 入力ト ラヒックモデルと待ち時間
シミュレーション
本研究では，ATM におけるランダムセルトラヒックを表現できる幾何過程（ Geo：Geometric Process ）モ率の推定はできない．
以上のような限界を克服するため IS 法が適用され

デルとバーストセルトラヒックを表現できるマルコフ変る．そこで，IS 法の原理について，正規分布の例を用

調幾何過程（ MMGP：Markov Modulated Geomtric
Process ）モデルを考える．
3. 1 幾何過程に従う到着間隔
1 タイムスロットにおけるセルの到着率を q としてセルとセルの到着間隔が τ となる確率 P (τ ) がいて以下に述べる（ 文献 [5] 参照）．

2. 1 正規分布における IS 法の例
X を平均 0, 分散 1 の正規分布 P0 に従う確率変数とし ，a > 0 に対して確率 P0 (X > a) の推定値をシ
ミュレーションで求めるとする（ 図 1 参照）．
a の値が大きいと {X > a} は希少事象となり，通常の MC 法では推定値が得られない．そこで，正規分布の平均値を少しずらし て平均 µ，分散 1 の正規分
布 Pµ を考え，この Pµ に従う確率変数 X̃ の標本値
x̃1 , . . . , x̃n を発生させる．図 1 より，x̃i > a となる

P (τ ) = q(1 − q)τ

(τ = 0, 1, 2, . . .)(2)と表されるとき，到着間隔は到着率 q の幾何過程に従うという．

確率は大きいので，前節で述べた限界が克服され る．

3. 2 2 状態 MMGP に従う到着間隔
2 状態 MMGP は，到着率 q0 の幾何分布に従う到
着間隔でセルが到着する状態 S = 0 と，到着率 q1 のそして，平均値の変化による頻度の変化分を重み関数

幾何分布による状態 S = 1 とをマルコフ的に遷移す

P0 (x̃i )/Pµ (x̃i ) によって補正して，IS 推定値として
1 
P0 (x̃i )δa (x˜i )NPµ (x̃i )N

P̂IS (X > a) =る過程である．状態 S = 0 から S = 1 への状態遷移
確率を r01 , S = 1 から S = 0 へは r10 とする．シ

(1)

i=1

を得る．ただし ，δa (x) は x > a のとき 1, x <=aのとき 0 となる関数である．P̂IS (X > a) は P0 (X > a)の不偏推定量である [1]．平均値 µ をうまく選ぶと IS
推定値の分散を MC 推定値の分散より小さくすることができる．シミュレーション回数を N とすると，分散は N −1 に比例するので，IS 推定値と MC 推定値のミュレーションでは，あるタイムスロットにセルが到着したとすると，そこでマルコフ連鎖に従って新たに状態を決定し ，その状態が S = i のとき到着率 qi の幾何分布によって次のセルの到着時刻を決め，以下これを繰り返す．

3. 3 Geo/D/1 における待ち時間の基本漸化式
Geo/D/1 待ち行列において，n 番目のセルの待ち
時間を W (n) とし ，セルの到着間隔を表す到着率 qの幾何分布に従う確率変数を τ とすると，次の漸化式が成り立つ．

W (n + 1) = max(0, W (n) + 1 − τ )(3)

式 (3) より W (n) はマルコフ連鎖をなすことがわかる．

W (n + 1) − W (n) = y とおくと，W (n) = m という条件のもとでの y の条件付き確率 Pm (y) は m > 0のとき

Pm (y) =図 1 正規分布に対する IS シミュレーション
Fig. 1 IS simulation for normal distribution.

q(1 − q)1−y(1 − q)m+1(−m < y <= 1)
(y = −m)

(4)

となる．
1783

電子情報通信学会論文誌 ’99/10 Vol. J82–B No. 10

3. 4 2 状態 MMGP/D/1 における待ち時間の基
本漸化式

待ち行列モデルに IS 法を適用したときの待ち時間
特性 P (W > w) の推定値と重み関数を以下に示す．22 状態 MMGP/D/1 待ち行列において，n 番目のセルの到着時点での状態を S(n)，到着率 qi の幾何分布

状態 MMGP/D/1 の P (W > w) に対する IS 推定値

P̂IS (W > w) は，RC を用いて，に従う確率変数を τi , i = 0, 1 とすると，次の漸化式

1 
δw (W (n))Gk (n)A

W (n + 1) = max(0, W (n) + 1 − τS(n) )(5)P̂IS (W > w) =S(n) は マ ル コ フ 連 鎖 に 従 うか ら ，式 (5) よ り 組(W (n), S(n)) は マルコフ連鎖をなすことが わか る．
W (n + 1) − W (n) = y, S(n + 1) = j とおくと ，W (n) = m, S(n) = i という条件のもとでの (y, j) の条件付き確率 Pm,i (y, j) は，m > 0 のときPm,i (y, j) =qi (1 − qi )1−y rij(1 − qi )m+1 rij(−m < y <= 1)
(y = −m)
(6)

Nk

A

が成り立つ．

k=1 n=1Nk

1 Gk (n)
BB

k=1 n=1(8)となる．ただし ，A，B はシミュレーション中に発生した RC の数，Nk ，Nk は k 番目の RC 中に発生したセル個数，Gk (n) は重み関数である．式 (8) の重み関
数 Gk (n) は，Gk (n)

 p(W ∗ (c), S ∗ (c)|W ∗ (c − 1), S ∗ (c − 1))

n−1

=

となる．

3. 5 待ち時間特性 P (W > w) の推定値

c=1

W を定常状態における待ち時間，すなわち W =limn→∞ W (n) とする．W がある待ち時間 w を超える確率の MC 法による推定値 P̂ (W > w) は再生サイ
クル（ RC：Regenerative Cycle ）の概念 [3] を用いて
表すと，

1 
δw (W (n))
AA

P̂ (W > w) =Nk

k=1 n=1
B

1  
Nk
B

k

k

k

k

p∗ (Wk∗ (c), Sk∗ (c)|Wk∗ (c − 1), Sk∗ (c − 1))(9)となる．ここで ，分母 p∗ (Wk∗ (c), Sk∗ (c)|Wk∗ (c − 1),

Sk∗ (c − 1)) は ，シ ミュレ ーション 分布に 従って発生し た c − 1 番目のセルに対する待ち時間と状態の組(Wk∗ (c − 1), Sk∗ (c − 1)) が c 番目の (Wk∗ (c), Sk∗ (c))
へ遷移するシ ミュレ ーション 分布による遷移確率を表し ，分子 p(Wk∗ (c), Sk∗ (c)|Wk∗ (c − 1), Sk∗ (c − 1)) は

(7)

k=1

となる．ただし ，δw (W (n)) は W (n) > w のとき 1，

W (n) <
= w のとき 0 となる関数を表す．A，B はシミュレーション中に発生した RC の数を表し，Nk ，Nkは k 番目の RC 中で発生したセル個数を表す．

4. P (W > w) に対する IS 推定値
大きな w の値に対し て {W > w} は希少事象とシ ミュレ ーション 分布に 従って発生し た c − 1 番目のセルに対する (Wk∗ (c − 1), Sk∗ (c − 1)) が c 番目の(Wk∗ (c), Sk∗ (c)) へ遷移するもとの分布による確率を表す．
具体的に，Geo/D/1 及び 2 状態 MMGP/D/1 待ち
時間の重み関数は以下のとおりになる．
まず，Geo/D/1 待ち時間の重み関数 Gk (n) は，

 q(1 − q)τ ∗ (c)n−1

Gk (n) =c=1q ∗ (1 − q ∗ )τ ∗ (c)(10)

な る．そこで ，待ち行列に おけ る IS 法では ，事象となる．ただし，q ∗ はシミュレーション分布の到着率{W > w} の発生頻度を高めるために入力トラヒック

であり，τ ∗ (c) は c 番目と c + 1 番目の到着間隔を表の到着率を増加させてシミュレーションを行い，最後す到着率 q ∗ の幾何分布を表す．ここで，シミュレーに重み関数で補正し て P (W > w) に対する推定値

ション分布の到着率を q ∗ としたが，IS 法では q ∗ をを得る．シミュレーションで用いる分布をシミュレー

ど のように決定するかが重要となる．q ∗ の最適値が，ション分布という．

次章で最適シミュレーション分布によって与えられる．

1784

論文／ ATM バーストセルトラヒックの待ち時間特性に対する IS シミュレーション

次に ，2 状 態 MMGP/D/1 待 ち 時 間の 重み 関 数
Gk (n) は，∗
 qS(c) (1 − qS(c) )τS(c)rS(c)S(c+1)

n−1

Gk (n) =∗∗c=1 qS(c) (1 − qS(c) )

τ∗
S(c)

∗
rS(c)S(c+1)(11)

方法には ，最も一般的な従来法とし て以下の二つがある．
（ 1 ） 事前シミュレーション法による導出 [10]（ 2 ） 理論的な導出 [1], [4]これらに加えて，本研究では最適なシミュレーション
分布を求める第 3 番目の手法として，となる．ただし ，S(c) は c 番目のセルに対する状態（ 3 ） 近似解析による導出∗
を表し，qi∗ , rij
はシミュレーション分布の到着率及び遷移確率である．また，τi∗ は到着率 qi∗ の幾何分布を∗表す．これら qi∗ , rijの最適値が，最適シミュレーショ

を提案する．

ン分布によって与えられる．

上記
（1）
〜
（3）
の手法について順に説明する．説明のために Geo/D/1 の待ち時間を例にする．

5. 1 事前シミュレーション法による導出

4. 1 ダイナミックインポータンスサンプリング
待ち時間に対する IS シミュレーションでは到着率

ション分布で IS シミュレーションを行い，得られたを増加させてシミュレーションを行うために RC にお推定値の平均が 安定し ている領域で分散が 最小になけるセル個数が増加する．したがって，RC 当りのシる部分の到着率を最適到着率 q ∗ として決定する方法

ミュレーション時間も増加してしまう．そこで，IS 法

である [10]．実際に事前シミュレーション法を用いて事前シミュレーション法とは，いろいろなシミュレー

ではシミュレーション時間を減少させるためにダ イナ

Geo/D/1 待ち時間の最適シミュレーション分布を決定

ミックインポータンスサンプ リング（ DIS：Dynamic

した例を示す．図 2 は，到着率 q = 0.4 の Geo/D/1

Importance Sampling ）という手法を用いて高速化を行う [3]．
DIS は，IS 法で各 RC において十分な待ち時間事象を得た後に到着率を下げて待ち時間を 0 にすることでRC を強制的に終了させる手法である．ここで，十分

待ち時間における事前シミュレーションの例である．

れぞれ対数表示したものである．シミュレーション分

な待ち時間事象とは，ある待ち時間 w を超えるセル布の到着率は 0.45 から 0.9 まで 0.0045 刻みで変化さ

図 2 の横軸はシミュレーション分布の到着率，左縦軸は待ち時間 W が待ち時間 50 を超える確率 P (W > 50)の推定値 100 個の標本平均，右縦軸は標本分散を，そ個数が 15 から 30 程度発生した事象をいう．この DIS

せ，その刻みごとに P (W > 50) の推定値の標本平均という手法を導入することで，RC 当りのシミュレーと標本分散を示した．1 標本当りの発生 RC 数は 1,000

ション時間が減少しシミュレーションの高速化が図れである．図 2 において，P (W > 50) の推定値の平均る．本研究では，すべての IS シミュレーションに DISが安定している領域で分散が最小の部分の q ∗ は 0.6法を用いている．

5. 最適シミュレーション分布
IS 法でシミュレーションに用いる分布をシミュレー
ション分布という．このとき重要なことは良いシミュ
レ ーション 分布を使用することである．良いシ ミュレーション分布とは，IS 法によって得られる推定値の分散が小さい分布である．待ち行列に IS 法を適用する場合，勝手なシミュレーション分布を使用すると，
重み関数で補正し ても正し い値を得ることができな
い [5], [10]．IS 推定値の分散が最小となるシミュレー
ション分布を最適シミュレーション分布といい，この最適シミュレーション分布で IS シミュレーションを行うことで正しい推定値を得ることができる．
以下に最適シミュレーション分布を求めるいくつかの手法をあげる．最適シミュレーション分布を求める

図 2 シミュレーション分布の到着率 q∗ に対する IS 推定
値の標本平均と標本分散（ Geo/D/1 ）Fig. 2 Sample average and sample variance of IS estimate vs. arrival rate of simulation distribution(Geo/D/1).

1785

電子情報通信学会論文誌 ’99/10 Vol. J82–B No. 10である．したがって，この 0.6 が最適な到着率となる．

e2θ qeθ − (1 − q)=

しかし ，事前シミュレーションもシミュレーションであるため一つのシミュレーション結果を得るのに大+変な時間を必要とする．また，図 2 を見てわかるよ
うに ，分散が 最小の所でも少し 幅があるために 正確eθ − 1e−mθ (1 − q)m+2 (15)eθ − (1 − q)となる．そこで，式 (13) に基づいて，な最適シミュレーション部分を見極めることが困難となる．したがって，もう少し簡単に正確な最適シミュ
レーション分布を求める必要がある．

5. 2 理論的な導出
IS シミュレーションに large deviation 理論 [1], [4]を適用し IS 推定値の分散を最小にする最適シミュレー
ション分布を得る．文献 [1], [2] によれば以下の結果がMm (θ) = 1(16)の解 θ = θ∗ を用いて，式 (14) より最適なシミュレー
ション分布が∗eθ y Pm (y)(17)と与えられる．

得られる．
一般に，次の Markov 連鎖モデル {xn } を考える．

xn+1 = xn + V (xn , ξn )(n = 0, 1, . . . , )(12)

実際に数値的に得られた θ∗ の値を表 1 に示す．表 1を見てわかるように，理論的解析法では最適な θ∗ は待ち時間 m に依存する．つまり，待ち時間 m ごとに

ここで，ξn は i.i.d.（ independent and identically dis-

式 (16) を数値的に解く必要があり，非常に計算量がtributed；独立同分布 ）確率変数で，V は xn と ξn

多くなってしまう．したがって，更に簡単に最適な θ∗の deterministic な 関数で あ る．xn = x と い う条を求める必要がある．

件の もとで V (xn , ξn ) = y とな る条件付き 確率を
Px (y) ≡ P (V (xn , ξn ) = y|xn = x) とし ，Px (y) に関するモーメント母関数を Mx (θ) =eθy Px (y) とy布を求める近似解析を提案する．この近似解析は，理論的解析法に 近似を入れ ることで ，最適シ ミュレ ー

すると，確率 P (xn > y) に対する IS シミュレーショ
ンにおける最適シミュレーション分布は，

Mx (θ) = 1

5. 3 近似解析による導出
本研究では高速かつ簡単な最適シミュレーション分

ション分布を近似的に求める方法である．

5. 3. 1 Geo/D/1 におけ る最適シミュレ ーション
(13)

分布の到着率に対する近似値
式 (15) より，を満たす解 θ = θ∗ によって∗P (V (xn , ξn∗ ) = y|xn = x) = eθ y Px (y)



(14)

で決まる ξn∗ によって与えられる．
以上の理論的な導出法を Geo/D/1 待ち時間シミュ(1 − q)m+1
1
dMm (θ) 
=2− +

dθ
q
q
θ=0

(18)

となる．Geo/D/1 待ち行列の安定条件は q < 1/2 でレーションに当てはめて説明する．まず，Geo/D/1 の待ち時間において，式 (3) より W (n) がマルコフ連鎖をなすことがわかるので，上記の理論的な導出法が適

表 1 理論解析による Geo/D/1 の最適な θ ∗ の値
Table 1 The optimal θ ∗ of Geo/D/1 by theoreticalanalysis.

用可能である．
まず，m > 0 のとき Pm (y) のモーメント 母関数

Mm (θ) は式 (4) より，
Mm (θ) =1
eθy Pm (y)y=−m= e−mθ (1 − q)m+1+1

y=−m+11786

eθy q(1 − q)1−y

m=1
m=2
m=3
m=4
m=5
m=6
m=7
.
.
.
m = 50.
.
.

q = 0.1
2.092
2.188
2.196
2.197
2.197
2.197
2.197
.
.
.
2.197
.
.
.

q = 0.2
1.163
1.351
1.380
1.385
1.386
1.386
1.386
.
.
.
1.386
.
.
.

q = 0.3
0.491
0.771
0.827
0.841
0.846
0.847
0.847
.
.
.
0.847
.
.
.

q = 0.4
0.000
0.271
0.361
0.389
0.399
0.403
0.404
.
.
.
0.405
.
.
.

論文／ ATM バーストセルトラヒックの待ち時間特性に対する IS シミュレーションあるから（ 後述の式 (37) 参照），式 (18) よりm>

まず，式 (6) より m > 0 のとき Pm,i (y, j) のモー
メント母関数 Mm,i (θ1 , θ2 ) は，log(1 − 2q)−1
log(1 − q)(19)

Mm,i (θ1 , θ2 )



(26)


1

1

のとき dMm (θ)/dθ|θ=0 < 0 となる．Mm (θ)|θ=0 = 1=であるから，式 (19) を満足する m に対して，方程式

Mm (θ) = 1 は解 θ = θ∗ > 0 をもつことがわかる．したがって，式 (15) より，m が大きいとき

rij eθ2 (j−i)

j=0



y=−m



1

=

rij e

θ2 (j−i)

j=0

∗

e2θ q
Mm (θ )  θ∗e − (1 − q)
∗

e2θ q

eθ∗ − (1 − q)(27)と計算される．m が十分大きいとき，式 (20) と同様

1

(21)
Mm,i (θ1 , θ2 ) (22)∗となる．Pm(y) を式 (4) の形の分布として，その到着率を q ∗ とすると，
∗
Pm
(y) =

e2θ1 qirij eθ2 (j−i) θ
e 1 − (1 − qi )
(28)

1−q
q



1

j=0

となり，式 (21) より∗

に式 (27) は∗eθ e2θ1 qieθ1 − (1 − qi )eθ1 − 1e−mθ1 (1 − qi )m+2
+ θ
1e − (1 − qi )(20)と近似される．そこで式 (16)，(20) より，eθ1 y Pm,i (y, j)q ∗ (1 − q ∗ )1−y
(1 − q ∗ )m+1

する方程式
1

j=0

(−m < y <= 1)
(y = −m)e2θ1 qirij eθ2 (j−i) θ=1e 1 − (1 − qi )(i = 0, 1)
(29)

(23)
となるが，そこで

q∗ = 1 − qと近似され る．式 (13)，(28) に基づいて θ1 , θ2 に関(24)とおくと，式 (4)，(22)〜 (24) より，m が大きいとき，計算によりを 以下のよ うに 解く．行列 R = (Rij ) とベクトルz = (zi ) を
Rij = rije2θ1 qiθe 1 − (1 − qi )zi = eθ2 i(i = 0, 1)(i, j = 0, 1)(31)によって定義すると，式 (29) はRz = z∗∗Pm(y)  eθ y Pm (y)

(−m <
=y<
= 1)

(25)

となることが容易に確かめられる．したがって，式 (17)

(30)

(32)

と書ける．式 (32) より R は固有値 1 をもつので，固有方程式 det(R − E) = 0（ E は 2 × 2 単位行列）より解 θ1 = θ1∗ が得られる．また，z がその固有ベクト∗
より，Pm(y)，すなわち式 (24) の到着率 q ∗ = 1 − q をルであることから解 θ2 = θ2∗ が zi の定数倍を除いてもつ分布が近似的に最適シミュレーション分布となる．

得られる（ 文献 [4], [5] 参照）．得られた解 θ1∗ , θ2∗ より

このように，Geo/D/1 待ち時間に近似的解析手法
を用いると事前シミュレーションや数値的解析なしで，最適な到着率を求めることができる．

5. 3. 2 2 状態 MMGP/D/1 に おけ る最適シ ミュ
レーション分布の到着率と遷移確率に対す

最適シミュレーション分布のパラメータ∗qi∗ = 1 − e−θ1 (1 − qi )
∗
2θ1

∗
e qi
∗
rij
= eθ2 (j−i) θ∗rij
e 1 − (1 − qi )ミュレーション分布を近似的に求める．

(33)(i, j = 0, 1)
(34)

る近似値
次に，2 状態 MMGP/D/1 の IS 法における最適シ(i = 0, 1)∗が得られる．式 (33)，(34) により得られた qi∗ , rijの数値例を表 2 に示す．
1787

電子情報通信学会論文誌 ’99/10 Vol. J82–B No. 10
表2近似解析による 2 状態 MMGP/D/1 の最適シミュ
レーション分布の到着率と遷移確率の近似値
Table 2 Approximating values for the optimal arrival
rate and the optimal transition probabilitiesof 2-state MMGP/D/1.
もとの到着率と遷移確率(q0 , q1 ; r01 , r10 )(0.10, 0.40; 0.90, 0.10)(0.10, 0.40; 0.20, 0.20)(0.10, 0.40; 0.50, 0.50)(0.10, 0.40; 0.50, 0.80)(0.10, 0.40; 0.10, 0.90)

最適な到着率と遷移確率の近似値∗∗(q0∗ , q1∗ ; r01, r10)(0.49, 0.66; 0.96, 0.04)(0.59, 0.72; 0.70, 0.03)(0.71, 0.81; 0.77, 0.22)(0.79, 0.86; 0.73, 0.60)(0.87, 0.91; 0.30, 0.70)

図 3 Geo/D/1 待ち時間特性 P (W > w)
Fig. 3 Waiting time probability P (W > w) of
Geo/D/1.

6. 待ち時間特性の理論値または近似値
6. 1 Geo/D/1 待ち時間特性の理論値
式 (3) よ り Geo/D/1 待 ち 時 間 の 状 態 遷 移 確 率

P = P (W (n + 1)|W (n)) は，



1−q

q

0

0

..
.

q(1 − q)q(1 − q)2..
.

q
q(1 − q)
..
.

0
q
..
.

 (1 − q)2
P =  (1 − q)3


···



···

···

..
.

(35)
となる．そして，方程式

pP = p(36)

を解いて待ち時間の定常分布 p = (p0 , p1 , p2 , . . .) を得る．定常分布 p より待ち時間特性 P (W > w) として

P (W > w) =

q
1−q

w+1
(37)

を得る．

case1
case2
case3
case4
case5

(q0 , q1 ; r01 , r10 )(0.10, 0.40; 0.90, 0.10)(0.10, 0.40; 0.20, 0.20)(0.10, 0.40; 0.50, 0.50)(0.10, 0.40; 0.50, 0.80)(0.10, 0.40; 0.10, 0.90)

図 4 2 状態 MMGP/D/1 待ち時間特性 P (W > w)
Fig. 4 Waiting time probability P (W > w) of 2-state
MMGP/D/1.

6. 2 2 状態 MMGP/D/1 待ち時間特性の近似値
文献 [4] における近似値の導出法を適用して，2 状態
MMGP/D/1 待ち時間特性 P (W > w) の近似値 P̃wである．比較のため，Geo/D/1 では式 (24) で示したを以下の手順で導出することができる．

到着率を用いて行った IS 法，MC 法，及び理論値を示

式 (29) を解いて得られた θ1 = θ1∗ によって，近似値

P̃w = e

∗
−θ1
(w+1)

(w = 0, 1, 2, . . .)

(38)

を得る．

7. シミュレーション結果と評価
7. 1 IS 法 を 用 い た Geo/D/1 及 び 2 状 態
MMGP/D/1 待ち時間特性の評価
図 3，図 4 は ，そ れ ぞ れ Geo/D/1 と 2 状 態
1788

MMGP/D/1 待ち時間特性シ ミュレ ーションの結果

し ，2 状態 MMGP/D/1 では式 (33)，(34) で示した到着率と遷移確率を用いて行った IS 法，MC 法，及び式 (38) の近似値を示す．また，MC 法ではシミュレー
ション回数を 106 とし，IS 法ではシミュレーション回
数を 103 としてシミュレーションを行った．
図 3，図 4 を見ると，MC 法では P (W > w) の推定値が 10−6 程度しか得られていないのに対して，IS法では本研究の目的である 10−9 オーダの P (W > w)の推定値も得ることができる．このように，IS 法を用

論文／ ATM バーストセルトラヒックの待ち時間特性に対する IS シミュレーション
表3

MC 法と IS 法の標本平均，標本分散とシミュレー
ション 時 間の 比 較（ 2 状 態 MMGP/D/1，(q0 =0.1, q1 = 0.4 ; r01 = 0.9, r10 = 0.1), P (W >15) の場合）
Table 3 Comparison of sample average, sample variance, and simulation time between MC
and IS methods (in the case of 2-state
MMGP/D/1, (q0 = 0.1, q1 = 0.4 ; r01 =0.9, r10 = 0.1), P (W > 15)).
MC

標本平均
標本分散
RC 数
1RC の平均シミュレー
ション時間（ 秒）
シミュレーション
時間比

ミュレーションと標本数 100 の結果から標本平均と標本分散を計算して，95％信頼区間長を求めた．表 4 より，IS 法は MC 法より信頼区間長が 2 けた小さいことがわかる．つまり，IS 法は MC 法より真値に対するばらつきが少ない推定値を与える．

8. む

す び

IS

1.26 × 10−4

1.19 × 10−4

2.54 × 10−10
1062.40 × 10−10300

る IS 法の有 効性を 示すために ，Geo/D/1，2 状態

3.17 × 10−5

3.16 × 10−4法を適用した．その結果，動画像伝送等で必要とされ

334

1

る 10−9 オーダの待ち時間特性を高速に得ることがで

本 研 究 で は ，待 ち 時 間シ ミュレ ーション に 対 す

MMGP/D/1 待ち行列モデルの待ち時間に対して IS

きた．
表 4 MC 法 と IS 法 の 95％ 信 頼 区 間 の 比 較（ 2 状 態
MMGP/D/1，(q0 = 0.1, q1 = 0.4 ; r01 =0.9, r10 = 0.1), P (W > 15)，RC 数=30,000の場合）
Table 4 Comparison of 95％ conﬁdence intervals between MC and IS methods (in the caseof 2-state MMGP/D/1, (q0 = 0.1, q1 =0.4 ; r01 = 0.9, r10 = 0.1), P (W > 15),
number of RC = 30, 000).
MC 法

95％信頼区間長
3.55 × 10−4

IS 法

−6

3.61 × 10

また，IS 法を待ち行列に用いる場合，最適シミュ
レーション分布を求めることが重要となるが，最適なシミュレーション分布を求める方法として，本研究では近似解析法を提案した．この近似解析法を使用することで，シミュレーションや多くの数値解析を行うことなく最適シミュレーション分布が求められた．本来の最適シミュレーション分布は現時点での状態 m あ
るいは (m, i) に依存するが，m を大きくしていくと一定の分布に収束するので，その収束した分布を最適
分布の代わりに使用することで計算量を大幅に削減できる，というのが提案した近似解析法のポイントであいて低い確率を求めることができることがわかる．

7. 2 高 速 性る．したがって，このような収束の特徴をもつ他の待ち行列モデルへの広い応用が考えられる．

次に，IS 法が MC 法より高速であるということを示すためにシミュレーション時間比を求めた．

文[1]

表 3 は ，2 状 態 MMGP/D/1 待 ち 時 間 に お け

1990.
[2]

tic Algorithms,” IEEE Trans. Automat. Contr.,

レーション時間の比較結果である．標本平均，標本分

vol.AC-28, no.9, pp.907–920, Sept. 1983
[3]

eters for Eﬃcient Simulation of Communication Net-

求めるための標本数は 100 とした．結果として表 3 をworks,” IEEE/ACM Trans. Networking, vol.1, no.3,

見てわかるように，IS 法と MC 法の標本分散を同程度

pp.293–305, June 1993.
[4]

vol.E80-A, no.11, pp.2238–2244, Nov. 1997.
[5]

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

表 4 は ，2 状 態 MMGP/D/1 に お け る (q0 =0.1, q1 = 0.4, r01 = 0.9, r10 = 0.1), P (W > 15)の場合の待ち時間特性に対する 95％信頼区間長を評価したものである．信頼区間は RC 数 30,000 回のシ

K. Nakagawa, “The Importance Sampling Simulation of MMPP/D/1 Queueing,” IEICE Trans. Fund.,

7. 3 95％信頼区間長
次に，IS 法と MC 法の 95％信頼区間長を評価する．

M. Devetsikiotis and J.K. Townsend, “Statistical Optimization of Dynamic Importance Sampling Param-

る RC 数を調節した．ここで，標本平均と標本分散をにしたとき，IS 法は MC 法より約 334 倍高速である．

M. Cottrell, J-C. Fort, and G. Malgouyres, “Large
Deviations and Rare Events in the Study of Stochas-

本平均，標本分散を同程度にそろえたときのシ ミュ
散を同程度になるようにシミュレーション中に発生す

J.A. Bucklew, Large Deviation Techniques in Decision, Simulation, and Estimation, Wiley, New York,る (q0 = 0.1, q1 = 0.4 ; r01 = 0.9, r10 = 0.1),

P (W > 15) の 場合に おいて ，IS 法と MC 法の 標

献[6]

T. Sakai and H. Ogiwara, Importance Sampling for
TCM Scheme over Non-Gaussian Noise Channel,”
IEICE, Trans. Fund., vol.E78-A, no.9, Sept. 1995.

1789電子情報通信学会論文誌 ’99/10 Vol. J82–B No. 10[7]鈴木信吾 ，中川健治 ，“バックプレッシャ制御を 用いた
ATM 交換機の IS 法に よ る特性評価 ，” 信学論（ B-I ），vol.J81-B-I, no.6, pp.371–377, June 1998.

[8]

上田健志，中川健治，“待ち行列シミュレーションへの IS

[9]

法の適用，
” 平 8 信学信越支大，pp.215–216, Oct. 1996.
上田健志，“IS 法を用いた ATM 交換機における待ち時間[10]

Q. Wan and V.S. Frost, “Eﬃcient Estimation of Cell

特性の解析，” 長岡技科大大学院修士論文，1998.
Blocking Probability for ATM Systems,” IEEE/ACM
Trans. Networking, vol.1, no.2, pp.230–235, April
1993.

（ 平成 10 年 7 月 2 日受付，11 年 1 月 6 日再受付）

上田健志 （ 正員）

平 8 長岡技科大・工・電気・電子システ
ム卒．平 10 同大大学院修士課程了．同年
NEC 航空宇宙システム（ 株）入社．

中川

健治 （ 正員）

昭 55 東工大・理・数学卒．昭 60 同大大学院博士課程満期退学．昭 60 NTT 研究所入社．平 4 長岡技術科学大学工学部助教
授．統計的情報理論，待ち行列理論など の研究に従事．理博．IEEE ，情報理論とその応用学会，日本数学会各会員．

1790