JP6430360B2 - シミュレーション装置、シミュレーション方法及びシミュレーションプログラム - Google Patents

Description

本発明は、トラヒックシミュレーションにおいて、統計的に有効な数のサンプルが得られた時点で損失率を計算する装置、方法及びプログラムに関する。

例えば、インターネットにおいて、動画配信が盛んに行われており、また、社内ＩＰ（Ｉｎｔｅｒｎｅｔ Ｐｒｏｔｏｃｏｌ）ネットワーク等のプライベートネットワークをはじめとして、研修や訓示や講演等の映像や音声を同時に各社員の受信装置（ＰＣ等）に多地点配信（マルチキャスト）することが普及しつつある。

多地点配信を行うマルチキャストサーバからネットワークを介して送られる映像や音声等のマルチキャストパケットはＵＤＰ（Ｕｓｅｒ Ｄａｔａｇｒａｍ Ｐｒｏｔｏｃｏｌ）であり、ＴＣＰによる再送処理がされないため、映像アプリケーションによるＦＥＣ（Ｆｏｒｗａｒｄ Ｅｒｒｏｒ Ｃｏｒｒｅｃｔｉｏｎ）等の保護が働くとしても、パケット損失の影響が映像や音声の品質に表れてしまうことがある。このため、様々なトラヒック形態やトラヒック量に応じてネットワーク内で発生するパケット損失を、事前に予測しておくことは重要である。

また、これまで様々な近似計算式が提案されているが、トラヒック形態やトラヒック量によってどの程度の誤差があるか不明な場合が大半であるため、実際にはシミュレーションも行ってデータグラムの損失数を求め、それを用いて近似計算式の正しさを検証する必要がある。このように、トラヒックシミュレーションは現在でも重要な手段である。

尚、ここでいうデータグラムとは、ＩＰパケットを始めとする各種通信プロトコルのパケットやイーサネットフレーム、ＡＴＭ（Ａｓｙｎｃｈｒｏｎｏｕｓ Ｔｒａｎｓｆｅｒ Ｍｏｄｅ）セルやＸ．２５パケット等の何らかの情報を転送するための一塊のビット列を指している。

シミュレーションを行ってパケット損失を求める場合、バッファ長が長くなるにつれて、すなわち、損失率が極めて小さくなるにつれて、１個の損失を得られるまでのシミュレーション時間が長くなる。このため、シミュレーション時間を短くするために、できるだけ少ない損失数でデータグラムの損失率を得たいが、データグラムの損失数が何個で統計的に有効か、具体的な方法はこれまで見当たらず、非特許文献１や非特許文献２に書かれているように、経験的にシミュレーションを行っていた。

そこで、シミュレーション時間を短くするために、統計的手法を用いて、できるだけ少ない損失数でデータグラムの損失率を得る方法を考案した（例えば、特許文献１参照。）。

特開２０１３−１９２００５号公報

村田正幸、大崎博之、原井洋明、宮原秀夫、"通信トラヒック理論とその応用〔ＶＩ〕−シミュレーション技法とその適用例−、"信学誌、ｖｏｌ．７８、ｎｏ．３、ｐｐ．２６４−２７０、Ｍａｒ．１９９５． 小川耕司、中川健治、"インポータンスサンプリングシミュレーション法における最適シミュレーション分布、"信学総大、通信（２）、Ｂ−６４６、ｐ．７８、Ｍａｒ．１９９６． Ｒ．Ｓ．バーリントン・Ｄ．Ｃ．メイ著、林知己夫・脇本和昌監訳、「確率・統計解析ハンドブック」、森北出版、１９７５．

特許文献１の方法では、図１のデータグラムの損失数Ｌ_ｊの時にシミュレーションの終了条件を満たしたとすると、開始からそれまでにシミュレーションで発生させたデータグラムの入力数Ｎ_ｉを用いて、
（数１）
Ｒ＝Ｌ_ｊ／Ｎ_ｉ （１）
としてデータグラムの損失率Ｒを求めている。

しかし、本来は任意の入力数Ｎに対して、それまでに発生した損失数Ｌを用いて損失率Ｒを計算する方がより正確な損失率を得られると考えられる。

ここで、次の（Ｌ_ｊ＋１）個目の損失が入力数（Ｎ_ｉ＋ｋ＋１）のときに発生したとする。損失率を計算する入力数を、損失数Ｌ_ｊとは無関係に、無作為に決めるとすると、
入力数Ｎ_ｉのときに損失率を計算する確率Ｐ_ｉ、
入力数Ｎ_ｉ＋１のときに損失率を計算する確率Ｐ_ｉ＋１、
・・・
入力数Ｎ_ｉ＋ｋのときに損失率を計算する確率Ｐ_ｉ＋ｋ、
は全て等しいと考えられる。

上記の場合の損失数は全てＬ_ｊ個であり、これらの平均入力数は（Ｎ_ｉ＋ｋ／２）である。よって、より正確な損失率Ｒは
（数２）
Ｒ＝Ｌ_ｊ／（Ｎ_ｉ＋ｋ／２） （２）
となる。

一般的な長時間のシミュレーションの場合、Ｎ_ｉ≫ｋ／２となるため、
（数３）
Ｌ_ｊ／（Ｎ_ｉ＋ｋ／２）≒Ｌ_ｊ／Ｎ_ｉ （３）
となり、式（１）と式（２）のどちらの式で損失率Ｒを計算しても値は殆ど同じであり、実質的に問題ない。

しかしながら、シミュレーション時間を短くするために、できるだけ少ない損失数で損失率Ｒを求める場合には、Ｎ_ｉは小さい値であるため、より正確な損失率Ｒを得るためには、ｋ／２を無視することはできない。例えば、特許文献１では、早ければＬ_ｊ＝２０でシミュレーションを終了する。

このため、ｋ／２を求める必要があるが、ｋ／２を求めるためには、シミュレーション時間を短くしたいにもかかわらず、さらに損失数が（Ｌ_ｊ＋１）個になるまでシミュレーションを行う必要がある。

そこで、損失数がＬ_ｊ個になった時点で、先に述べたｋ／２を考慮して損失率を計算する方法が必要である。

本発明に係るシミュレーション装置は、
トラヒックシミュレーションにおいて、
データグラムの損失数が設定されたＬ_ｊ個になったときに、
データグラムの損失数が（Ｌ_ｊ−１）個のときのデータグラムの入力数Ｎ_ｍ、及び、損失数が前記Ｌ_ｊ個のときのデータグラムの入力数Ｎ_ｉと前記Ｎ_ｍとの差分値ｋ_ｄ、すなわち、式（４）から得られた値、及び、前記Ｌ_ｊを式（５）に適用するか、または、前記Ｎ_ｍ、及び、前記Ｎ_ｉ、及び、前記Ｌ_ｊを式（５’）に適用することで、データグラムの損失率Ｒを算出する。

本発明に係るシミュレーション装置は、
トラヒックシミュレーションにおいて、
データグラムの損失数が設定されたＬ_ｊ個になったときに、
損失数がＬ_ｊ個のときのデータグラムの入力数Ｎ_ｉ、及び、シミュレーションの開始から損失数がＬ_ｊ個になるまでの“１損失当たりの平均入力数”であるｋ_ａ、すなわち、式（６）から得られた値、及び、前記Ｌ_ｊを式（７）に適用するか、または、前記Ｎ_ｉ、及び、前記Ｌ_ｊを式（７’）に適用することで、データグラムの損失率Ｒを算出する。

本発明に係るシミュレーション方法は、
コンピュータが実行するトラヒックシミュレーション方法であって、
データグラムの損失数が設定されたＬ_ｊ個になったときに、
データグラムの損失数が（Ｌ_ｊ−１）個のときのデータグラムの入力数Ｎ_ｍ、及び、損失数がＬ_ｊ個のときのデータグラムの入力数Ｎ_ｉと前記Ｎ_ｍとの差分値ｋ_ｄ、すなわち、式（４）から得られた値、及び、前記Ｌ_ｊを式（５）に適用するか、または、前記Ｎ_ｍ、及び、前記Ｎ_ｉ、及び、前記Ｌ_ｊを式（５’）に適用することで、データグラムの損失率Ｒを算出するステップ、
を備える。

本発明に係るシミュレーション方法は、
コンピュータが実行するトラヒックシミュレーション方法であって、
データグラムの損失数が設定されたＬ_ｊ個になったときに、
損失数がＬ_ｊ個のときのデータグラムの入力数Ｎ_ｉ、及び、シミュレーションの開始から損失数がＬ_ｊ個になるまでの“１損失当たりの平均入力数”であるｋ_ａ、すなわち、式（６）から得られた値、及び、前記Ｌ_ｊを式（７）に適用するか、または、前記Ｎ_ｉ、及び、前記Ｌ_ｊを式（７’）に適用することで、データグラムの損失率Ｒを算出するステップと、
を備える。

本発明に係るプログラムは、本発明に係るシミュレーション方法に備わる各ステップをコンピュータに実行させるためのシミュレーションプログラムである。

上述のように、特許文献１の方法では、ｋ／２を考慮せずに損失率Ｒを計算していたが、本発明により、ｋ／２を考慮した損失率Ｒを求めることができるため、より正確な損失率を得ることができる。

本発明におけるシミュレーションのデータグラムの入力と損失発生の一例。 本発明の第３の実施の形態におけるシミュレーション装置のブロック図である。 本発明の第３の実施の形態におけるシミュレーションのフローチャートである。 シミュレーションのトラヒックを発生させる２状態のＭＭＰＰ（Ｍａｒｋｏｖ Ｍｏｄｕｌａｔｅｄ Ｐｏｉｓｓｏｎ Ｐｒｏｃｅｓｓ）の遷移図である。

以下、本発明の実施形態について、図面を参照しながら詳細に説明する。なお、本発明は、以下に示す実施形態に限定されるものではない。これらの実施の例は例示に過ぎず、本発明は当業者の知識に基づいて種々の変更、改良を施した形態で実施することができる。なお、本明細書及び図面において符号が同じ構成要素は、相互に同一のものを示すものとする。

［第１の実施の形態］
ｋ／２を考慮して損失率Ｒを求める方法の本発明の第１の実施の形態として、シミュレーション中に、データグラムの損失数が（Ｌ_ｊ−１）個になったときのデータグラムの入力数Ｎ_ｍ個を記録しておき、データグラムの損失数がＬ_ｊ個になったときのデータグラムの入力数Ｎ_ｉ個を用いて、
（数４）
ｋ_ｄ＝Ｎ_ｉ−Ｎ_ｍ （４）
として、
（数５）
Ｒ＝（Ｌ_ｊ−１）／（Ｎ_ｍ＋ｋ_ｄ／２） （５）
または、直接、
Ｒ＝２（Ｌ_ｊ−１）／（Ｎ_ｍ＋Ｎ_ｉ） （５’）
を計算する、ことを行う。

一例として、特許文献１の発明に本発明の内容を適用して、より正確な損失率を得ることができる。上記方法は、コンピュータ上にソフトウェアで実現することができる。

また、上述の動作を行う電子回路をハードウェア上に実装し、シミュレーションを行うことができる。ハードウェアとしては、論理素子を結合させた電子回路であり、ロジックＩＣ（Ｉｎｔｅｇｒａｔｅｄ Ｃｉｒｃｕｉｔ）等を基板上に実装して実現してもよいし、ＦＰＧＡ（Ｆｉｅｌｄ−Ｐｒｏｇｒａｍｍａｂｌｅ Ｇａｔｅ Ａｒｒａｙ）やＰＡＬ（Ｐｒｏｇｒａｍｍａｂｌｅ Ａｒｒａｙ Ｌｏｇｉｃ）等で実現してもよいし、ＬＳＩ（Ｌａｒｇｅ Ｓｃａｌｅ Ｉｎｔｅｇｒａｔｉｏｎ）中に実現してもよい。

尚、出力数＝入力数Ｎ_ｍ−損失数（Ｌ_ｊ−１）であることから、入力数Ｎ_ｍに比べて損失数（Ｌ_ｊ−１）が極めて小さい場合（例えば、損失率＝１０^−６の場合、損失数（Ｌ_ｊ−１）：入力数Ｎ_ｍ＝１：１，０００，０００）には、出力数≒入力数になるため、本実施の形態において、入力数Ｎ_ｍの代わりに損失数が（Ｌ_ｊ−１）個になったときの出力数を用い、入力数Ｎ_ｉの代わりに損失数がＬ_ｊ個になったときの出力数を用いて計算する方法もある。

［第２の実施の形態］
ｋ／２を考慮して損失率Ｒを求める方法の本発明の第２の実施の形態として、ｋに、シミュレーションを開始してからそれまでの“１損失当たりの平均入力数”であるｋ_ａを用いる、即ち、データグラムの損失数がＬ_ｊ個になったときのデータグラムの入力数Ｎ_ｉ個を用いて、
（数６）
ｋ_ａ＝Ｎ_ｉ／Ｌ_ｊ （６）
として、
（数７）
Ｒ＝Ｌ_ｊ／（Ｎ_ｉ＋ｋ_ａ／２） （７）
または、直接、
Ｒ＝２Ｌ_ｊ ^２／（Ｎ_ｉ×（２Ｌ_ｊ＋１）） （７’）
を計算する、ことを行う。

一例として、特許文献１の発明に本発明の内容を適用して、より正確な損失率を得ることができる。上記方法は、コンピュータ上にソフトウェアで実現することができる。

また、上述の動作を行う電子回路をハードウェア上に実装し、シミュレーションを行うことができる。ハードウェアとしては、論理素子を結合させた電子回路であり、ロジックＩＣ等を基板上に実装して実現してもよいし、ＦＰＧＡやＰＡＬ等で実現してもよいし、ＬＳＩ中に実現してもよい。

尚、出力数＝入力数Ｎ_ｉ−損失数Ｌ_ｊであることから、入力数Ｎ_ｉに比べて損失数Ｌ_ｊが極めて小さい場合（例えば、損失率＝１０^−６の場合、損失数Ｌ_ｊ：入力数Ｎ_ｉ＝１：１，０００，０００）には、出力数≒入力数になるため、本実施の形態において、入力数Ｎ_ｉの代わりに損失数がＬ_ｊ個になったときの出力数を用いて計算する方法もある。

［第３の実施の形態］
本実施形態では、特許文献１の発明への第１の実施の形態の適用例について説明する。

図２に、本実施形態に係るシミュレーション装置の一例を示す。本実施形態に係るシミュレーション装置は、記憶部１１と、データグラム生成部１２と、バッファ部１３と、損失検出部１４と、時間間隔測定部１５と、時間間隔判定部１６と、χ^２検定部１７と、損失率算出部１８と、発生回数判定部１９と、を備える。χ^２検定部１７は、実度数計数部２１と、実度数判定部２２と、統計量算出部２３と、統計量判定部２４と、を備える。

本実施形態に係るシミュレーション方法について、一例として、ケンドールの記法（または記号）によるＭＭＰＰ／Ｄ／１／Ｋのシミュレーションで実施する場合の動作について、図３に示すフローチャートを参照しながら説明する。本実施形態に係るシミュレーション方法は、損失検出手順と、時間間隔測定手順と、統計量判定手順と、シミュレーション終了手順を順に有する。損失検出手順では、ステップＳ１０１〜ステップＳ１０３を実行する。時間間隔測定手順では、ステップＳ１０４〜Ｓ１０５を実行する。統計量判定手順では、ステップＳ１０６〜Ｓ１１０を実行する。シミュレーション終了手順では、ステップＳ１１１〜Ｓ１１２を実行する。

本発明を使用する場合には、以下に述べる動作をシミュレーション中に行う。
ステップＳ１００では、シミュレーションに先立ち、記憶部１１が、設定された初期値を記憶する（詳しくは後述）。
ステップＳ１０１では、データグラム生成部１２が、データグラムを生成してトラヒックを発生させる。本実施形態の場合、データグラムは固定長であり、以下、これを単位長とする。データグラムの生成方法は、任意であるが、例えば、コンピュータで擬似乱数を用いて生成する。本実施形態では、データグラム生成部１２は、記憶部１１に設定されたＭＭＰＰパラメータを取得し、当該パラメータに従って、データグラムを生成する。

ステップＳ１０２では、データグラム生成部１２が、生成したデータグラムをバッファ部１３に出力する。これにより、バッファ部１３にデータグラムが書き込まれ、状況が損失検出部１４に通知される。バッファ部１３に書き込まれたデータグラムは、本実施形態の場合、１単位時間に１単位長（＝１データーグラム）が読み出される。尚、バッファ部１３は、生成したデータグラムをバッファに実際に書き込んでもよいが、少なくとも、空き又は使用中のバッファ量を管理でき、バッファに空きがないために損失が発生したことを把握できればよい。

ステップＳ１０３では、損失検出部１４が、バッファ部１３からの通知により、損失が発生したか否かを判定し、損失を検出する。

損失が発生した場合、損失検出部１４は、データグラムの損失数Ｌ_ｊ及びデータグラムの入力数Ｎ_ｉを記憶部１１に記憶する。データグラムの入力数Ｎ_ｉは、バッファ部１３に書き込まれたデータグラムの数である。また損失が発生した場合、損失検出部１４は、損失が発生した時刻を時間間隔測定部１５へ出力し、ステップＳ１０４へ移行する。一方、損失が発生しない場合、損失検出部１４は、次のデータグラムを生成する旨をデータグラム生成部１２へ通知し、ステップＳ１０１へ移行する。
ステップＳ１０４では、時間間隔測定部１５が、損失検出部１４から取得した時刻を用いて、損失が発生した時間間隔ｔを測定する。

ステップＳ１０５では、時間間隔判定部１６が、時間間隔測定部１５の測定した時間間隔ｔが閾値ＴＨ以上であるか否かを判定する。これにより、損失の発生時間間隔により、ランダムな損失（ランダムに発生する損失）と引き続く損失（ランダムに発生する損失に引き続く損失）に分類する。時間間隔ｔが閾値ＴＨ以上である場合、すなわちランダムな損失である場合、時間間隔判定部１６は、時間間隔ｔを実度数計数部２１へ通知してステップＳ１０６へ移行し、以下χ^２検定を行う。一方、ステップＳ１０５において時間間隔ｔが閾値ＴＨ未満である場合、すなわちランダムな損失でない場合、時間間隔判定部１６は、次のデータグラムを生成する旨をデータグラム生成部１２へ通知し、ステップＳ１０１のデータグラム生成に戻る。

ステップＳ１０６では、実度数計数部２１が、通知されたランダムな損失の発生時間間隔を記憶し、ランダムな損失の発生回数ｍを１つインクリメントする。ここで、本実施形態における「ｍ」は損失数Ｌ_ｊに相当する。
ステップＳ１０７では、実度数計数部２１が、それまでに記憶したランダムな損失の発生時間間隔を用いて、平均発生時間間隔を計算する。
次に、実度数計数部２１が、計算した平均発生時間間隔を用いて、χ^２検定のためのＬ個の区間の境目となる（Ｌ−１）個の境界値を計算する。
次に、実度数計数部２１が、記憶しているランダムな損失の発生時間間隔を使って各区間ｉ（ｉ＝１，２，・・・，Ｌ）の実度数ｆ_ｉを計数する。

ステップＳ１０８では、実度数判定部２２が、各区間ｉの実度数ｆ_ｉがすべてｆ_ｍｉｎ以上であるか否かを判定する。
すべての区間において実度数ｆ_ｉがｆ_ｍｉｎ以上である場合、実度数判定部２２は、各区間ｉの実度数ｆ_ｉを統計量算出部２３へ出力してステップＳ１０９へ移行する。
一方、実度数ｆ_ｉが５未満の区間がある場合、実度数判定部２２は、次のデータグラムを生成する旨をデータグラム生成部１２へ通知し、ステップＳ１０１のデータグラム生成に戻る。

ステップＳ１０９では、統計量算出部２３が、各区間ｉの実度数ｆ_ｉからχ^２検定用の統計量χ^２を計算する。

ステップＳ１１０では、統計量判定部２４が、統計量χ^２が臨界値χ_０ ^２よりも大きいか否かを判定する。統計量χ^２≦臨界値χ_０ ^２の場合、統計量判定部２４は、その旨を損失率算出部１８へ通知してステップＳ１１１へ移行する。一方、統計量χ^２＞臨界値χ_０ ^２の場合、統計量判定部２４は、ランダムな損失の発生回数ｍを発生回数判定部１９へ出力してステップＳ１１３へ移行する。

ステップＳ１１１では、損失率算出部１８が、損失検出部１４の検出した全損失数と、生成された全データグラムの数から、損失率Ｒを計算する。このとき、損失率算出部１８は、式（４）を用いて「ｋ_ｄ」を求め、式（５）を用いて損失率Ｒを計算するか、或いは、式（５’）を用いて損失率Ｒを計算するか、或いは、式（６）を用いて「ｋ_ａ」を求め、式（７）を用いて損失率Ｒを計算するか、或いは、式（７’）を用いて損失率Ｒを計算する。
ステップＳ１１２では、損失率算出部１８が、計算結果である損失率を出力する。

ステップＳ１１３では、発生回数判定部１９は、これまで発生したランダムな損失の発生回数ｍが最大発生回数Ｍよりも多いか否かを判定する。
ｍ≦Ｍの場合、発生回数判定部１９は、次のデータグラムを生成する旨をデータグラム生成部１２へ通知し、ステップＳ１０１へ移行してデータグラム生成に戻る。
ｍ＞Ｍの場合、発生回数判定部１９は、ステップＳ１１４を実行し、その旨を出力する。

ＭＭＰＰ／Ｄ／１／Ｋのシミュレーションの例の場合の、ステップＳ１００で設定する初期値、即ち、図３のフローチャートの初期設定部分の各項目について、初期値をどのような値に設定するのか等を、以下に説明する。

（ａ）ＭＭＰＰパラメータ値
ＭＭＰＰは、トラヒック発生形態の一つであるＭａｒｋｏｖ Ｍｏｄｕｌａｔｅｄ Ｐｏｉｓｓｏｎ Ｐｒｏｃｅｓｓのことであり、図４に示すような、ランダムに切り替わる複数の異なる密度でトラヒックが発生するモデルである。ＭＭＰＰは複数の異なる密度のトラヒック発生状態を持つことができるが、図４の最も簡単な２状態を例に説明する。
・λ_０とλ_１は、状態０と状態１のそれぞれの状態で発生するデータグラムの平均到着率であり、通常、０以上１以下の値である。
・γ_０１とγ_１０は、それぞれ、ＭＭＰＰの状態０から状態１への遷移確率と、状態１から状態０への遷移確率であり、０以上１以下の値である。同様に、１−γ_０１と１−γ_１０は、それぞれ、状態０と状態１が続く確率である。

（ｂ）バッファサイズ（Ｋ−１）
発生したＭＭＰＰトラヒックが流入するバッファ部１３の大きさである。

（ｃ）しきい値ＴＨ
損失が発生する度に、その直前の損失との発生時間間隔ｔを計算する。しきい値ＴＨは、事前に設定した発生時間間隔ｔの基準値であり、これと比較することにより、ランダムな損失と引き続く損失に分類する。

（ｄ）臨界値χ_０ ^２
臨界値χ_０ ^２を設定するためには、まず、有意水準αを設定する必要がある。有意水準αは、それまでに得られた標本が指数分布に従っている場合に、標本から得られた統計量χ^２が臨界値χ_０ ^２を超える確率（第１種過誤）がαである、という意味である。本実施形態では、標本にランダムな損失の発生時間間隔を用いる。
αの値は、通常５％または１％が用いられる。本実施形態では、後述する様に、自由度ｄ＝１としており、臨界値χ_０ ^２は統計量χ^２の分布の統計数値表より、
α＝５％の時，χ_０ ^２＝３．８４１
α＝１％の時，χ_０ ^２＝６．６３５
である。

（ｅ）ランダムな損失数の最大値Ｍ
統計量χ^２は損失数が増加すると増減しつつ次第に大きくなる傾向があるため、ランダムな損失数ｍが数百個程度までに終了しない場合には、統計量χ^２≦臨界値χ_０ ^２にならない可能性も考えられる。この場合に対応するために、シミュレーションを終了させるランダムな損失数ｍの最大値Ｍを設定しておく（ランダムな損失数ｍが最大値Ｍを超えたら終了（ステップＳ１１３〜Ｓ１１４））。
なお、Ｌとｆ_ｍｉｎについては後述する。

シミュレーション中に行う動作（ステップＳ１０５）で述べているランダムな損失と引き続く損失について説明する。
空きが大きい状態のバッファ部１３にＭＭＰＰトラヒックが入力することによって稀に損失が発生する場合には、損失はランダムに発生すると考えられる（ランダムな損失）。しかし、一度損失が発生した状態のバッファ部１３にＭＭＰＰトラヒックが入力することによって損失が発生する場合には、バッファ部１３内の空きがほとんど無いため、ゼロまたは極めて近い間隔で損失が発生することが多くなる（引き続く損失）。よって、前述のトラヒックモデルで発生する損失には、ランダムな損失と引き続く損失の２種類の損失が混在すると考えられる。

ランダムな損失の発生時間間隔は指数分布に従うと考えられ、ＭＭＰＰ／Ｄ／１／ＫのＤの出力時間すなわちバッファ部１３からのデータグラムの読み出し時間を１［単位時間］とすると、通常、発生時間間隔は数十万［単位時間］になり、発生時間間隔の存在確率は薄く広く分布する。一方、引き続く損失の発生時間間隔はゼロまたはゼロに近い値、例えば１００［単位時間］以下になることから、厳密ではないが、しきい値ＴＨにより分離できる。

シミュレーション中に行う動作（ステップＳ１０６〜Ｓ１１０）のχ^２検定について概要を説明する。
シミュレーションで発生するランダムな損失の発生時間間隔（以下、ランダム損失発生時間間隔と記す）が指数分布に従っているかどうかの一致度合い（適合度）を数値的に評価し判定する方法として、χ^２検定がある。χ^２検定は、標本の分布が母集団（理論）分布にどの程度一致しているかの検定によく用いられる方法である。ランダムに発生する事象の生起間隔は指数分布になることが知られており、このような場合は母集団の分布として指数分布を想定する。これ以外に、対象となる事象により、二項分布や正規分布など様々な分布が想定される。

χ^２検定では、まず、標本と母集団の分布の違いを統計量χ^２として求める。ここでは、標本はランダム損失発生時間間隔ｗ_ｊ（ｊ＝１，２，…，ｎ）であり、その集合をＷ_ｎとする。また、母集団分布は指数分布である。統計量χ^２は、一般に標本の分布が母集団分布に近づくほど値が小さくなる傾向がある。

次に、統計量χ^２と予め求められている臨界値χ_０ ^２を比較し、統計量χ^２が臨界値χ_０ ^２以下であれば、標本は有意水準α％で母集団分布であるとされる。具体的には、シミュレーションで順次ランダム損失発生時間間隔ｗ_ｊを得る度にその集合Ｗ_ｎに対してχ^２検定を行い、指数分布であると判定された場合には、その時点の損失数とバッファ部１３への入力数から損失率を計算し、シミュレーションを終了する。尚、臨界値χ_０ ^２は、自由度ｄと有意水準α％を用いて、複雑な計算により求めることができるが、よく使われる値は非特許文献３などに統計数値表として示されている。

シミュレーション中に行う動作（ステップＳ１０７）のχ^２検定の具体的な実装方法（４つの区間や境目となる３つの境界値など）について説明する。

・母集団分布のパラメータ
本実施形態では、母集団分布を指数分布とするが、そのパラメータ（平均と分散）は不明である。
そこで、標本からランダムな損失の平均発生時間間隔１／τを算出し、これを用いて母集団分布を推定する。指数分布では、平均値が１／τの時、分散は１／τ^２である。そこで、新たな“ランダムな損失”の発生による新たな発生時間間隔のデータを得られる度に平均値を計算し、平均値から分散を計算し、母集団分布のパラメータとする。

・χ^２検定の区間（カテゴリー）数Ｌ
χ^２検定では、シミュレーションで得られたランダムな損失の各発生時間間隔ｗ_ｊをその値により複数の区間の中のいずれかに分類し、実度数を計数する。このため、区間の数と各区間の範囲を決定する必要がある。
区間の数は、実際の検定では数個から数十個が使われているようであるが、計算時間の短縮及び計算負荷の軽減のためにできるだけ少ないことが要求される。
χ^２検定では、自由度をｄ、区間数をＬ、標本から最尤推定したパラメータ数をｓとするとき、
（数８）
ｄ＝Ｌ−１−ｓ≧１ （８）
である。上述したように、母集団分布のパラメータは標本から推定した平均及び分散であるためｓ＝２である。ゆえに、式（８）から
（数９）
Ｌ≧４ （９）
となる。よって、区間数Ｌ＝４とする。

・区間の最小実度数ｆ_ｍｉｎ
χ^２検定では、各区間の標本数も自由に設定できる。ただし、χ^２検定が成立する標本数は一般には５以上とされているため、本実施形態では、各区間の最小実度数ｆ_ｍｉｎを５とする。これにより、本実施形態に必要な最低損失数は２０である。

・各区間の範囲
標本数が少ない場合、標本の値には偏りが生じていることが多い。しかし、標本数が多くなると指数分布に一致することが示されているため、指数分布を想定して各区間の範囲を決定することにする。
χ^２検定を行うための区間数が４で、各区間の最小実度数ｆ_ｍｉｎが５と等しいことから、区間の範囲は、指数分布における４区間の各累積確率分布が等しくなる、すなわち、０．２５になる範囲とすることが妥当と考えられる。

４区間を指数分布（ｗ≧０）の左側から区間ｉ（ｉ＝１，…，４）と呼ぶことにし、区間ｉと区間ｉ＋１との境界の値を境界値ｙ_ｉ（但し、ｉ＝１，２，３）と呼ぶことにする。各境界値ｙ_ｉは、式（８）から、
（数１０）
０．２５×ｉ＝１−ｅｘｐ（−τｙ_ｉ） （１０）
これを解いて、
（数１１）
ｙ_ｉ＝−ｌｏｇ_ｅ（１−０．２５×ｉ）／τ （１１）
となる。

・χ^２検定の統計量χ^２
χ^２検定を行うために、それまでに得られた標本の各ランダム損失発生時間間隔ｗ_ｊを、式（１１）で得られた３つの境界値を用いて４つの区間に分類する。もし、いずれかの区間の実度数ｆ_ｉ（＝損失数）が最小実度数ｆ_ｍｉｎである５未満のときは、再度、次のランダム損失発生までシミュレーションを進める（ステップＳ１０８）。

そして、全区間の実度数ｆ_ｉが最小実度数ｆ_ｍｉｎである５以上になった場合には、以下の式（１２）の統計量χ^２を計算する（ステップＳ１０９）。

ただし、Ｆ_ｉは理論度数（実数）で、本実施形態では、

で表される。

このようにして求められた統計量χ^２を用いて、χ^２検定の臨界値χ_０ ^２と比較する（ステップＳ１１０）。

もし、統計量χ^２＞臨界値χ_０ ^２の場合は、それまでに得られた標本値が有意水準αで統計的に指数分布に従っていないと判定されたことになるので、再度、次のランダム損失発生までシミュレーションを進める。

もし、統計量χ^２≦臨界値χ_０ ^２となった場合には、ランダム損失については指数分布に従っていると判定されたことになり、損失率の計算（ステップＳ１１１）に進む。ステップＳ１１１では、それまでの全損失数を全入力数で割って損失率を計算し、シミュレーションを終了する。

なお、前述の動作を行う電子回路をハードウェア上に実装し、シミュレーションを行うことも可能である。ハードウェアとしては、論理素子を結合させた電子回路であり、ロジックＩＣ等を基板上に実装して実現してもよいし、ＦＰＧＡやＰＡＬ等で実現してもよいし、カスタムＬＳＩ中に実現してもよい。また、本発明の装置は、コンピュータとプログラムによっても実現でき、プログラムを記録媒体に記録することも、ネットワークを通して提供することも可能である。

以上、ＭＭＰＰ／Ｄ／１／Ｋのシミュレーションを行う場合について、本発明に係るシミュレーション方法、シミュレーション装置及びシミュレーションプログラムについて説明したが、本発明はＭＭＰＰ／Ｄ／１／Ｋのシミュレーションに限定されない。例えば、データグラムはＭＭＰＰ／Ｄ／１／Ｋ以外のＭ／Ｍ／１／Ｋ、Ｍ／Ｇ／１／Ｋ等のトラヒックシミュレーションにおいて用いられる任意の送受信データを用いることも可能である。尚、Ｋが有限の値の場合は有限バッファであり、Ｋを無限の値とする場合は無限バッファを意味する（無限バッファの場合は、Ｋが省略されて記される場合もある）。本発明は、ＩＰパケット、イーサネット（登録商標）フレーム、ＡＴＭセル、Ｘ．２５パケットなどの任意のデータグラムを使用した任意の通信システムのシミュレーションに適用され得る。

インターネットにおけるコンサート等のライブインターネット配信や、社内ＩＰネットワーク等のプライベートネットワークにおける研修や訓示や講演等の映像や音声を各社員の受信装置（ＰＣ等）に多地点配信（ＵＤＰを用いたマルチキャスト）する状況において、ネットワーク設計として、配信前に様々なトラヒック形態やトラヒック量における損失率を予測する必要がある。

そのような場合に、シミュレーションは現在でも有効な方法の一つであり、本発明の実施例で述べた統計的手法を用いてシミュレーションを終了する方法では、統計的に有効と判断された時点で損失率を計算することができるため、短時間にシミュレーションを終了させることができる。また、得られた損失率について、統計的妥当性を持たせることができる。

本発明は、上記の短時間でシミュレーションを終了させるような場合に、より正確な損失率を得られる計算方法を提供する。

１１：記憶部
１２：データグラム生成部
１３：バッファ部
１４：損失検出部
１５：時間間隔測定部
１６：時間間隔判定部
１７：χ^２検定部
１８：損失率算出部
１９：発生回数判定部
２１：実度数計数部
２２：実度数判定部
２３：統計量算出部
２４：統計量判定部

Claims (5)

1. トラヒックシミュレーションにおいて、
   データグラムの損失数が設定されたＬ_ｊ個になったときに、
   データグラムの損失数が（Ｌ_ｊ−１）個のときのデータグラムの入力数Ｎ_ｍ、及び、損失数が前記Ｌ_ｊ個のときのデータグラムの入力数Ｎ_ｉと前記Ｎ_ｍとの差分値ｋ_ｄ、及び、前記Ｌ_ｊを式（Ｃ１１）に適用するか、または、前記Ｎ_ｍ、及び、前記Ｎ_ｉ、及び、前記Ｌ_ｊを式（Ｃ１２）に適用することで、データグラムの損失率Ｒを算出する、
   シミュレーション装置。
   （数Ｃ１１）
   Ｒ＝（Ｌ_ｊ−１）／（Ｎ_ｍ＋ｋ_ｄ／２） （Ｃ１１）
   （数Ｃ１２）
   Ｒ＝２（Ｌ_ｊ−１）／（Ｎ_ｍ＋Ｎ_ｉ） （Ｃ１２）
2. トラヒックシミュレーションにおいて、
   データグラムの損失数が設定されたＬ_ｊ個になったときに、
   損失数がＬ_ｊ個のときのデータグラムの入力数Ｎ_ｉ、及び、シミュレーションの開始から損失数がＬ_ｊ個になるまでの“１損失当たりの平均入力数”であるｋ_ａ、及び、前記Ｌ_ｊを式（Ｃ２１）に適用するか、または、前記Ｎ_ｉ、及び、前記Ｌ_ｊを式（Ｃ２２）に適用することで、データグラムの損失率Ｒを算出する、
   シミュレーション装置。
   （数Ｃ２１）
   Ｒ＝Ｌ_ｊ／（Ｎ_ｉ＋ｋ_ａ／２） （Ｃ２１）
   （数Ｃ２２）
   Ｒ＝２Ｌ_ｊ ^２／（Ｎ_ｉ×（２Ｌ_ｊ＋１）） （Ｃ２２）
3. コンピュータが実行するトラヒックシミュレーション方法であって、
   データグラムの損失数が設定されたＬ_ｊ個になったときに、
   データグラムの損失数が（Ｌ_ｊ−１）個のときのデータグラムの入力数Ｎ_ｍ、及び、損失数がＬ_ｊ個のときのデータグラムの入力数Ｎ_ｉと前記Ｎ_ｍとの差分値ｋ_ｄ、及び、前記Ｌ_ｊを式（Ｃ３１）に適用するか、または、前記Ｎ_ｍ、及び、前記Ｎ_ｉ、及び、前記Ｌ_ｊを式（Ｃ３２）に適用することで、データグラムの損失率Ｒを算出する、
   シミュレーション方法。
   （数Ｃ３１）
   Ｒ＝（Ｌ_ｊ−１）／（Ｎ_ｍ＋ｋ_ｄ／２） （Ｃ３１）
   （数Ｃ３２）
   Ｒ＝２（Ｌ_ｊ−１）／（Ｎ_ｍ＋Ｎ_ｉ） （Ｃ３２）
4. コンピュータが実行するトラヒックシミュレーション方法であって、
   データグラムの損失数が設定されたＬ_ｊ個になったときに、
   損失数がＬ_ｊ個のときのデータグラムの入力数Ｎ_ｉ、及び、シミュレーションの開始から損失数がＬ_ｊ個になるまでの“１損失当たりの平均入力数”であるｋ_ａ、及び、前記Ｌ_ｊを式（Ｃ４１）に適用するか、または、前記Ｎ_ｉ、及び、前記Ｌ_ｊを式（Ｃ４２）に適用することで、データグラムの損失率Ｒを算出する、
   シミュレーション方法。
   （数Ｃ４１）
   Ｒ＝Ｌ_ｊ／（Ｎ_ｉ＋ｋ_ａ／２） （Ｃ４１）
   （数Ｃ４２）
   Ｒ＝２Ｌ_ｊ ^２／（Ｎ_ｉ×（２Ｌ_ｊ＋１）） （Ｃ４２）
5. 請求項３又は４に記載の各動作をコンピュータに実行させるためのシミュレーションプログラム。

Images