シミュ V ーションソフトと統計的側面
若山邦紘
川
111111111
目11川
刷
11川
11111
川111
川
111
川
附
11川
11111
川
1111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111②システムの状態を記述している変数の値が離徴的に1. はじめに変化する，今回のシミュレージョンソフト特集は OR ソフト特集ものであることを認識しておいていただきたい.状態変の第 2 弾として企画されたものである.これに筆者が O化をもたらす出来事を事象とかイベントと呼ぶが，事象R 誌編集委員長のときに，森戸晋さんにオーガナイザはある時刻に瞬間的に起こり，その瞬間に状態が変化すーをお願いしたのであるが，そのつけがまわり，る，という前提でそデルは記述されるのである.

ミイラ

とりがミイラになった次第である.

シミュレージョンモデル，それも特に離散系モデルを

2. なぜ乱数か?

計算機のプログラムとして記述しようとするとき，このたとえば，ほかほか弁当屋に客がランダムにやってきプログラムが一般の数値計算用プログラムとかなり異なて，何人分かの弁当を注文し貿って帰るとする.ある客った形式のプログラムとなるため，1950年代の後半にシの到着から次の客の到着までの時間間隔が一定ではなくミュレーションプログラムが簡単に書けるように GASP

ランダムだとすると，その時間間隔をランダムに決めるという離散系シミュレーションのための FORTRANためのランダムメカニズムが必要になる.また何個買うサブルーチンが Pritsker によって発表され，まもなくかを決めるのにもランダムメカニズムが L 、る.注文をしといった本格的な離散系のシミて金を払い弁当を受け取るまでの所要時間もランダムメュレーション言語が，それぞれ IBM ， RAND で開発さカユズムで決めたい.このようなランダム要因を含んだGPSS , SIMSCRIPTれた.しみじみと「オブジェクト指向言語が新しいなん

シミュレーション分析に乱数を使用するのは無作為にデていったって，ータを抽出し，煽りの少ない結果を推定しようとする標30数年前に SIMCSCRIPT はオブジェグト指向だったんだもんねえ Jとは，当時 SIMCSCRIPで大規模なシミュレーションを実施していた小池将貴さ

本調査の応用である.
具体的にはランダムメカニズムをそれぞれの要因につん{元三菱電機，現筑波技術短期大学)の言である.小いて，実績データをもとにして確率分布を決定したり，池さんは古い恋人を想い出すようにまだ惚れ込んでいる要因の構造から推察される理論的分布を当てはめたりのです.ただし SIMSCRCIPT は玄人好みであるのでし，これらの分布に従う乱数を発生させて，ランダムに素人にはあまり利用されてこなかったのが実情である.

起こる現象を模擬してゆくので、ある.

このマユュアルに，シミュレーション言語にはまず問弁当慶にやってくる客の到着時間間隔が指数分布に従題をモデル化するための世界観があることが重要だとさっているとか，その客が買ってゆく品物の数がポアソンれ，分布に従っているとか，そのサーピス時聞がアーラン分①シミュレーショ γ の進行を制御するためのタイミン
グルーチ γ布に従っていると古代、った場合，シミュレーションの実行中にそれらの分布に従う乱数を引きながら，事象の発②種々の分布に従う乱数の発生

生時刻をコントロールしたりシステムの状態を表わして③結果の統計と報告書作成

L 、る変数の値を変化させる.現在のところ不確定要素のの機能を備えていなければならないとある.
また離散系の基本的コンセプトが，①ジステムの状態変化が離散的な時刻に起こり，シミュレーションにおけるランダムメカニズムは乱数が基本である.

そこでシミュレーションソフトは各種の乱数を使いやすい関数の形で用意している.

わかやま

くにひろ法政大学工学部

千 184 小金井市梶野町 3 ー 7-2

588 (10)

比較的新しい SLAMnと WITNESS で利用できる確率分布を表 1 に示す.ちなみに SLAMn は FORTRAN をベースに， WITNE© 日本オペレーションズ・リサーチ学会. 無断複写・複製・転載を禁ず.

オベレーションズ・リサーチ

表 1

用意されている乱数発生用関数

確率分布の種類をプログラムをわざわざ書かなくても自動的に集計してSLAMn~WI!NESSくれるような仕組みが組み込まれている.実際，シミュレーションのプログラムを自分で書くと半分以上がこのI
山]一様分布(実数) I(轍) I。

。

。

。

三角分布。

。

指数分布。

。

アーラン分布。

。

ガンマ分布。

変えての実験となるとその組合せば膨大になる.むかしワイプル分布。

1 セットのシミュレーションが終わると何百ページかのベータ分布。

正規分布。

。に，結果を生のデータそのままで見せるのではなく，誰対数正規分布。

。

が見てもわかりやすい形のレポートにまとめる機能が備。

わっているとありがたい.この稿ではふれないが出力形

臥 1] 一様分布(実数)

ような結果の統計を出力するためのステートメントになってしまう.シミュレーションソフトはモデルの基本構造を組み立てるだけで，このような煩雑さからユーザー。

を解放してくれるのである.
ランダムな実験は，やるたびに結果が異なるであろうから，より高い精度を求めるならば多数回のシミュレーションが必要であろう.また，いろいろなパラメータをアウトプット用紙が出てくるのをよく見たものである.
これを見て結果を解析するのは大変な仕事となる.さら二項分布式のひとつとしてシミュレーションの実行経過を視覚的

ポアソン分布。

。

ユーザー定義の離散型分布|。

。

ソフトが時流となっている.ワークステーションの普及。

によりモデルの動きの可視化が実用的になってきた.に訴えるアユメーション機能を含んだシミュレージョン

ユーザー定義の連続型分布|

4.

SS は C をベースにした言語で、ある.

しかし「乱数を用いてシミュレーションを実行したがシミュレーションの観測データのバラツキ

結果がどうもおかしい.乱数は信用できるのか J といっ誰しもが観測データをもとにして何らかの推定値を求た問いかけが多い.これは確かに乱数がおかしい場合もめるためにシミュレーションを行なっている.街単な Mあるが，実験のやり方が悪い場合も考えられる.

川!/ I のシミュレーションから平均待ち時間を次のよう

伏見 [2J には 16 ピット程度のパソコンでは発生されるに計算してみた.

乱数の性質に欠陥があることが解説され，これに対して①初期状態、は空の状態から始める.

M系列を用いた乱数発生法がよいことが述べられてい②客の待ち時間を記録し， N 人までの待ち時間を平均る.関根，高橋，若山 [4J にはいろいろな乱数発生アルゴリズムが詳しく述べられている.いずれにしても，32ピットで演算する計算機なら，通常のシミュレーションでは普通の乱数プログラムの基本として用いられている合同法による乱数に問題はないと考えてよい.

3.

シミュレーションソフトの結果の統計機能する.

③N を 5 万人まで増やしてゆきながら平均待ち時間を計算する.

④乱数の初期値を変えて実験をくりかえす.
その結果が図 1 である.太線はこの実験から得られた95%信頼区間を示す.つまり「平均待ち時間はこの幅の中に入るだろう」といえば， 95% はあたること(ひっくり返せば当たらない確率が 5%) を意味する.信頼区間

シミュレーションソフト，とりわけ離散系シミュレーの幅が広すぎて満足で、きなければ人数をもっと増やさなションソフトは待ち行列のシミュレータであると言いきければならない.実験回数を増やせないとすれば，あたっても叱られまい.そこで，システムの中で発生する待る確率を犠牲にした信頼区間を示すしかない.

ち行列に関する滞留数と滞留時間に関するヒストグラム

図 2 はこの M川1/1 の理論的解析からわかっている待や平均値，分散，最大値，最小値などの基本的な統計量

ち時間の平均 μ と分散 σ2 をもっ母集団から互いに独立

1993 年 11 月号© 日本オペレーションズ・リサーチ学会. 無断複写・複製・転載を禁ず.

(11)587

なデータ Xh X2t … ， XN を抽出し，

120

この実験を

100

内川

N_l~t=ln川U

S2=~L: ~=1 (Xt‑X)2h
町

平均待ち時間

X=(X1+ X 2+ ・・・ +xN)/N

UAHVτAHυ
ge4
2

くりかえした場合のグラフである.周知のようとしたとき，と S2 は μ と σZ の不偏推定量であ。

るから， 95%信頼区聞は1。

X::!
::2s/‑IN234

実行の長き(万人)図 1と表わされる.つまり，信頼区間の幅はデータ数のルートに反比例する.別の言い方をすれば，平均待ち時間の変化の様子

120

推定値の精度を現在より 1 桁あげようとするな

100

らば，現在の 100 倍の実験が必要となる.これ

Z80をルート法則という.この事実から，シミュレ: 60ーションであまり精度にこだわると莫大な計算

信 40

時間を覚悟しなければならない.

20

グラフを一見すれば，平均待ち時間の信頼区

ハυ。リ間の幅は独立なデータの平均の場合に比べては1るかに広い.じつはわれわれのシミュレーションにおいてi 番目の客の待ち時間を叩h サー

ビス時聞を Sj ，2345実行の長き(万人)

図 2

独立なサンプルの場合

tt.i+1 を次の客との到着間隔と=，2 とし，自己相関のある待ち時間を Wh W 2 ， … ， WN とすると，
叩 t+1=max (世}i+Si-tt. i+h 0)

すると，自己相関関数は r(k)=COV(Wi ， 却i+k)/，2 と書という関係があるので，自分の前の客の待ち時間が長いくことができる.そして，平均値は独立な場合と同じくと自分の待ち時間も長くなることが予想され，

事実 Wjと叩 i+1 には強い正の相関がある(図 3 参照)

このよう

üi=(w ， +wg+ …+叩 N)/Nによって推定できるが，分散の推定量は ，2/N とはならに 1 つの時系列のデータの聞に相関があるとき，この時ず，系列には自己相関があるという.叩4 と叩 i-k 相関をこの附)=ヲト{ 1 叫 ~:~(1 ーす )r (k)}実験のデータから計算すると図 4 のようなグラフになる.
このカーブを自己相関関数と呼ぴ ， E(wtl= μ ，と表わされる . r(k)=O ならば独立な場合そのものであV(wtlるが，正の自己相闘がある場合は川島 )>0 であり， {}の中の値を g と書くと相関のない場合に N個のデータから200得られるものと同じ精度を実現するためには ， gN 個の180-セ160:/.
:

140120VVi+1官

100~

80e、2::~，作'.

0.
8.
，・ー

相 0.6 H思1\'・・之.~・.
・~:叩数 0.4t
¥

0.
2。

。

1、k、、、200S

Wi (日)
図 3

引き続く客の待ち時間の相関関係

588(12)

図 48001,000

M/M/l の自己相関関数© 日本オペレーションズ・リサーチ学会. 無断複写・複製・転載を禁ず.

オベレーションズ・リサーチ

叫(川岨)

表 2

相関による標本数への影響

混雑率pg0.8

75.
090.
9330.
850.
951268.
48
図 E

データが必要となる.待ち行列のシミュレーショ
/結果から得られた g の値は表 2 のような結果に再生的過程を利用したサンプリング

1~なる.混雑がはげしい場合にはその影響の大きさに驚く.
われわれは標本調査におけるランダムザンプリ

変ングの考えをとり入れ，乱数を各所に使ったので係数あるが，実行の途中でシステムの状態が一度偏っ事.恥~一一一一司、+・今、-+.・4ド4・4・岬・_.~

動 0.1ー回』『‘
司崎角・IJ_..

たところにいってしまうと，しばらくはその影響が続き，観測データは独立性を失ってしまうとい

面扇...・・陶・‘司剛・渇句"‘・・ぬ岡崎・M・・“・・"隆司・z0.
0104

うまことに厄介なもの，それが確率的なシステムシミュレーションの宿命といえる.

図 85. 独立なデータの抽出法

8121620

サンプル数(万)再生過程を利用したサンプリング方法の精度(A が自己相関関数による方法，B が再生的過程を利用した方法)
このようなシミュレーションのデータの観測法に工夫をこらして，独立なデータを得るための考えについてふ

1 回の長い実験を分割したものと考えてよい.ただし，れよう.

第 1 ブロックの終了状態と第 2 プロックから得られるデ

簡単な方法でプロック形成法と L 、う方法がある.自己ータは相関をもってしまうことになるので，独立と見な相関の影響を打ち消すために観測データをまとめてプロせるようにするには，ブロックの大きさを十分に大きくックにして，それぞれのプロックの平均値をあらたなデとるか，プロック聞のデータを捨てるなどしてプロックータと見れば，これらは独立な観測値となる.ブロック

閥の独立性をできるだけ保つようにするなどが考えられの作り方にはいろいろ考えられているが，いちばん簡単る.しかし，全体の実験の長さが一定だとした場合，ウな方法は，乱数の初期値を変えてその実験ごとに平均待ォームアップもブロック間のデータを捨てることもやめ

ち時聞を計算したものをデータすれば，これらのデータて，プロック聞の独立性が保てるようなできるだけ小さは互いに独立になる.しかし，この実験では窓口が空のいブロックに分けてできるだけ数多くのデータを得るの状態からシミュレーションを始めているから却。 =0 であが賢明な策である.

る.したがって，はじめのうちの待ち時間の値は加。 =0

プロック形成による方法では，自己相関の高い系列にの影響を受けているはずだから小さい値に偏ったもので対してはプロックサイズを大きくしなければならないかあろうことが予想される.そこで，この影響が無視できら独立なデータの個数が少なくなるために推定値の分るくらいまでの待ち時間データは捨てて，それ以後の観散を小さくしようとした当初の目的が果たせなくなる.

測値から平均待ち時間を計算するのがよいといわれてい待ち行列はしばしば空になる.この時の待ち時間は Oる.これをシステムのウォームアップという.それぞれである.これ以前の状態はこれ以後の状態には影響を与の実験でウォームアップが必要で捨てるデータ量が馬鹿

えな L 、から，待ち時間ゼロを境にブロックを形成すれば各プロックの挙動は互いに独立である.これを「再生過にならない.

そこで回目の実験の終了時の状態は定常状態から程を利用したサンプリング jと呼ぶ(図 5 参照).

こののランダムなサンプルと考え，これを 2 回目の実験の初場合，ブロックサイズは異なるから，平均値の推定には期状態とすればウォームアップの無駄が省ける.これは加重平均を計算する.

1993 年 11 月号© 日本オペレーションズ・リサーチ学会. 無断複写・複製・転載を禁ず.

(1 3) 騎9より推定値の精度が高い.
(万)

80(%)

義 14
な 1240 百 l ・対照1の相関が期待できるかも知れない.しかし，モ30 基 L至宝16デルの中で乱数同士が複合的に作用して，出力2011棒ー

10120.
10.
070.
050.
030の債に負の相関が出てこないことがあるので必ずしもうまくゆくとは限らない . M/M/1 の待線

ち時間シミュレーションでは，精度(変動係数)

図 7{r I， r2' … ， rN} を乱数列とするシミュレーションを実行すると結果には負50 ~ I 図通常1プ。
lレ

グ… ， rN} としたとき ，60 与一一一

310数

もっとも簡単なのは，もとの乱数列が {r"r2'70 サ

叩 i+ ， =max (Wi+Si‑ti.i+"O)

精度の達成のために必要なサンプル数の比較 (p=O.8)

において ， Ui=Si-tü i+l の分布に従う乱数を発面= (b，w ， +b 2 W 2+ ・・・ +bMwM)/Nによって推定できるが，分散の推定量は，V( 弱)=I;正午{ 1+2I; ~~1'( 1 ーす ) r(k)}

生させるようにすれば，乱数を 1 種類しか使わないので対照変量法が効果的に働く.図 7 は通常の 2 種類の乱数を使った場合と比較したものである.

7. おわりにとなる.
図 S は再生過程を利用した場合の精度(変動係数で表わしている)を自己相関関数による方法と比較したものである.この方法は大変に小さい分散の推定値を与える.
同じ精度でよいのなら実験回数は約 1/3 ですむ.しかし，複雑なモデルでは再生点を定義しでも，システムの状態が再生点になかなか戻ってこないことになってしまうので，うまく使える場合が少ないのである.

6.

分散減少法の適用

数多いシミュレーションソフトが出回って，計算機プログラミングを知らない人でもモデルが作れるようになってしまった.だからシミュレーション結果はどんどん出てくる.しかし，シミュレーションの結果から「何がわかるのか J [IJ については，相変わらず何の面倒も見てはくれずユーザーまかせ. I シミュレーション公害が出なければよいが」と危倶するのは筆者ばかりではあるまし、.

M川f/l のシミュレーションにおける対照変量法(負参考文献相関法)の適用を考えてみよう.この方法はモデルへ負の相関をもつような入力をいれたとき，モデルから負の相関をもっ出力が得られるならば，これらの出力の平均[1J 逆瀬川浩孝「シミュレーションで何がわかるのか jOR 誌 32 ，1
1

値は分散が小さくなることを利用したものである.すな[2J

かち，相関をもっ 2 つの確率変数の平均値の分散は[3] 伏見他「特集:乱数とその応用 J OR 誌38 ， 12

伏見正則[出力結果の解析 J OR 誌 35 ，

2

[4J 関根，高橋，苦手山[シミュレーション」日科技連，V( 市)=す(1 +p)

1976となるので ， p<O の場合は独立なデータを平均したもの

570(14)© 日本オペレーションズ・リサーチ学会. 無断複写・複製・転載を禁ず.

オベレーションズ・リサーチ