最適化・シミュレーション演習
第8回 シミュレーションの種類
とモンテカルロ法
• 授業サポートページ
http://www.morito.mgmt.waseda.ac.jp/optsim/
• 数理計画による最適化と（離散事象型）シミュレーションに
関する演習を行う．使用するソフトウェアは，AMPL（＋
GurobiまたはCPLEX)，および，Simul8を想定している．
• 演習では，数理計画による最適化やシミュレーションの実践
的能力を身につけることを目指す。履修者は，Cを使用でき
る環境を有するPCを持参すること。受講者は，実験室にて，
演習で使用する C，AMPL，Simul8をダウンロードできる。
1
後半の主な内容（シミュレーション）
• 離散型シミュレーションの演習では，C，および，Simul8を用いて様々な待
ち行列系システム，離散事象システムのモデル化，実験の進め方，簡単
な分散減少法と結果の分析について学ぶ．さらに，シミュレーションを用
いた最適化の演習を行う．
•
第８回 11/20 シミュレーションの種類とモンテカルロシミュレーション
シミュレーションの種類，モンテカルロシミュレーション，一様乱数，C による円周率πの推定
• 第９回 11/27 Simul8によるシミュレーション
M/M/1，並列キャッシュマシンの並び方，レンタルビデオの本数
• 第１０回 12/04 ラベル，リソース，シミュレーション実験の進め方と結果の表示・分析
• 第１１回 12/11 やや複雑なSimul8のモデル，分散減少法
宅配ピザ屋のオーブン容量とバイク台数，分散減少法（円周率πの推定，M/M/1）
• 第１２回 12/18 離散事象型シミュレーションの考え方とCによるM/M/1モデルの実装
事象，事象カレンダ，事象ルーチン，事象制御ロジック，単一サーバ待ち行列（M/M/1)モデ
ル実装の考え方
• 第１３回 2015/01/08 CによるM/M/1モデルの実装 （ディスカッションと講評）
指数乱数など，実装，結果の分析，解析結果との比較
• 第１４回 2015/01/15 シミュレーションを用いた最適化，まとめ
サーバー数の最適化，分散減少法（M/M/1，円周率πの推定）
• 第１５回 未定 自作の最適化モデルとシミュレーションモデル，および，その分析
本日のトピックス
• いろいろな種類のシミュレーション
• 擬似乱数
– 乱数の発生方法
– 一様乱数
– 特定の分布に従う乱数
• モンテカルロ法（例：①定積分，②πの推定）
– 当たり外れのモンテカルロ
– 入門的モンテカルロ
3
いろいろなシミュレーション（その１）
• 連続型(continuous)シミュレーション
微分方程式や差分方程式で表現されるモ
デルのシミュレーション
• 離散型（discrete, または、離散事象型
discrete-event）シミュレーション
待ち行列型のモデルで表現されるシステム
のシミュレーション
– 銀行のキャッシュマシン、生協食堂、コールセ
ンター、病院、工場のショップ（ジョブショップ等、
電話、通信ネットワーク、．．．）
– 混雑現象のシミュレーション
• その他（モンテカルロシミュレーション、エー
ジェントシミュレーション等々）
4
いろいろなシミュレーション（その２）
• 確定的(deterministic)シミュレーション
– 連続型シミュレーションに多い
• 確率的（probabilistic, または、stochastic)シ
ミュレーション
– 離散型シミュレーションのほとんどは確率的シ
ミュレーション
– 計算機上で乱数（擬似乱数）を発生させて確率
的変動を再現する
5
擬似乱数
• シミュレーションで使用する乱数は、一見したと
ころはランダムにみえる、一定の生成方法に基
づく数列
• 再現可能であったり、一定の範囲で制御可能で
あることを活用
• なんらかの生成方法で擬似乱数列を生成し、
（多くの場合）それを一様乱数列に変換
• 一様乱数以外の、特定の分布（たとえば、正規
分布、指数分布、等）にしたがう乱数は、多くの
場合、一様乱数をもとに適当な方法で生成
⇒森戸・逆瀬川 2000
6
擬似乱数の作り方（その１）
• 二乗採中法
• 乗算合同法
𝑥n+1 ＝𝑎𝑥𝑛 （mod 𝑃）
– 𝑥𝑛 は0,1,2, . . . , 𝑃 − 1の𝑃通りしかない
– 𝑥0 , 𝑥1 , 𝑥2 , … と並べていくと、最初の𝑃 +
1個のなかに同じものが存在
– 𝑛番目の𝑥𝑛 と𝑚番目の𝑥𝑚 とが同じだと、
𝑥𝑛+1 と𝑥𝑚+1 、𝑥𝑛+2 と𝑥𝑚+2 、・・・は一致
– 周期性が存在 （周期が長ければ問題ない）
– 適正なパラメータ（𝑎, 𝑃）の選択
7
擬似乱数の作り方（その２）
• 混合合同法
𝑥n+1 ＝𝑎𝑥𝑛 + 𝑐（mod 𝑃）
– 𝑃 = 232 の場合、𝑐は任意の奇数、𝑎 − 1が4の倍数に
なるように選択することによって周期が𝑃となる
float urand (unsigned long *seed) {
static unsigned long a=77520261;
static float fm=.23283064E-9;
*seed *= a;
(*seed)++;
return ((*seed) * fm);
}
//森戸・逆瀬川2000，p.151
• M系列法、．．．
8
一様でない乱数の生成（その１）
関数関係を利用する方法
• 一定の分布にしたがう確率変数と一様分布にし
たがう確率変数の関係を利用する方法
• 正規乱数の場合（その１：＿＿＿定理に基づく方
法）
– [0,1]の連続一様分布を12個加えて正規分布に変換
する → 標準正規変量𝑁(0, 12 )に変換
– 𝑎, 𝑏
𝑏−𝑎
(𝑏−𝑎)2
の連続一様分布：平均＝ ,分散＝
2
12
• 正規乱数の場合（その２：Box-Muller法）
– 𝑢, 𝑣は0-1間の一様乱数のとき,以下の𝑥, 𝑦は独立な標
準正規分布にしたがう乱数
– 𝑥 =sin(2πu) −2 log 𝑣, 𝑦 =cos(2πu) −2 log 𝑣,
9
一様でない乱数の生成（その２）
逆変換法
• 特定の分布の累積分布関数
を𝐹(𝑥)とするとき、 𝐹(𝑥)の
逆関数𝐹 −1 (𝑥)と一様乱数𝑢
から𝐹 −1 𝑢 を求めて乱数と
する方法
• 指数乱数の場合
1
u
𝑭(𝒙)
– 平均𝑐の指数分布の分布関数
𝑥
(− 𝑐 )
𝐹（𝑥）＝1－𝑒
– 𝐹 −1 𝑢 = −𝑐 log(1 − 𝑢)
𝒚 = 𝑭−𝟏 (𝒖)
10
C（またはC++）の乱数
• Cは、rand関数という 0 から RAND_MAX の区間
の疑似乱数を発生する関数を提供
– RAND_MAX は、標準ヘッダファイル stdlib.h の中に定
義あり
• 疑似乱数の初期値設定にはsrand関数を提供
• Rand関数は、混合合同法により乱数生成
• 特定の分布（一様乱数、正規乱数、指数乱数、
等）にしたがう乱数を生成するためには関数を自
作する必要あり
• 「rand関数は、”お手軽に”乱数を生成する方法である。しか
し、rand関数にはいくつかの問題があることが知られている
ために、この方法を研究などに使用することは推奨されない
ということに注意されたい。あくまでも導入のための”お手軽”
な方法である。」
11
C による一様乱数の生成
【書式】
#include <stdlib.h>
int rand(void);
【説明】
0～RAND_MAX の間の疑似乱数を返します。
乱数とは、要は「でたらめな値」ですが、コンピュータ内ではある規則にしたがって「で
たらめな値」を生成しています。これを擬似乱数と言います。
rand を使って発生させた擬似乱数は srand関数 で乱数の発生系列を変更しない限り、
実行のたびに同じ擬似乱数を発生します。下の使用例のプログラムを実行すると、
実行のたびに同じ繰返しで擬似乱数を発生します。
【引数】 なし
【戻り値】
0～RAND_MAX の間の疑似乱数。RAND_MAX は stdlib.h の中でマクロ定義されてい
ます。
12
C における乱数の初期値設定：ｓｒａｎｄ
【書式】
#include <stdlib.h>
void srand(unsigned seed);
【引数】
unsigned seed ： 擬似乱数の発生系列を変更する種。
同じ数値を設定すると同じ繰返しでrand関数
は擬似乱数を発生する。
【戻り値】 なし
「本当に」ランダム（時刻依存）にしたい場合
#include <stdlib.h>
#include <stdio.h>
#include <time.h>
int main(void)
{
int i,j;
/* 乱数系列の変更 */
srand((unsigned) time(NULL));
/* 1～100の擬似乱数を10個ずつ10行発生 */
for (i=1; i<=10; i++) {
for (j=0; j<10; j++) {
printf("%3d ",rand()%100+1);
}
printf("\n");
}
return 0;
}
13
モンテカルロ法
• 乱数を用いてモデルの解を求める技法の総
称で、別名、確率実験
• 実物による実験が不可能とか、解析的扱い
が困難であるような場合に有用
• 一般に相当の計算量を必要とするが、どのよ
うな問題に対してでも適用できるところが強み
• 広義の多重積分の評価と等価
14
基本的なモンテカルロ法
当たり外れ(hit or miss)のモンテカルロ
1
𝑓
0
• 𝐼=
𝑥 𝑑𝑥の評価例
• 𝑥 ∈ [0,1] 上で0 ≤ 𝑓(𝑥) ≤ c とし、
0≤ 𝑦 ≤ 𝑐, 0≤ 𝑥 ≤ 1の矩形領
域を考慮
𝒄
• [0,1]上の一様乱数αと[0,c]上
の一様乱数βの組(α，β)を繰り
返しN組生成
• 点(α，β)が関数𝑓(𝑥)の下側に
ある回数，すなわち，β≤f(α)と 𝟎
なる回数がMであれば，cM/N
が𝐼の推定値
𝟏
15
モンテカルロ法（その２）
入門的(introductory)モンテカルロ
• 𝐼=
1
𝑓
0
𝑥 𝑑𝑥の評価例
• [0,1]上の一様乱数
𝛼𝑡 (t=1,…,N)を生成
• 𝑁
𝑡=1 𝑓( 𝛼𝑡 )/𝑁が𝐼の推定
値
𝒄
𝟎
𝟏
16
授業内演習（宿題８．１）
• πの値を、
1 当たり外れのモンテカルロ
（２）入門的モンテカルロ
で推定しなさい。同じ「サンプルサイズ」に対し
て、どちらの方法の「精度」がよいと考えます
か？
• 授業時間中に少なくとも一方の方法でプログ
ラムを動かせるように！
• 授業中にプログラムができたか否かに関わら
ず宿題8.1を提出すること．
17
宿題８
（プログラム、結果、考察を提示すること）
• 宿題８．１
πの値を(1)当たり外れのモンテカルロ、(2)入門的モンテカ
ルロ、で推定し、どちらの方法が精度のよい推定結果を生
むか調べなさい。さらに、いずれかの方法（両方でもよい）
を想定して、サンプルサイズと結果の精度との間にどのよう
な関係があるか調べて、考察しなさい。
• 宿題８．２
さいころの目をシミュレートするプログラムを作成した上で
実行し、いくつかの基本的な方法を用いて、得られた結果
が妥当と考えられる（または、考えられない）ことを示しなさ
い。解答にあたっては、シミュレーション結果をExcelを用い
て分析してもよい。
• 提出の必要のない課題（後の演習で使用予定）
平均λの指数分布にしたがう乱数を生成する関数を作れ。
18
参考文献
［１］森戸晋，逆瀬川浩孝，システムシミュレー
ション，朝倉書店，2000．
［２］D.P.Kroese, T.Taimre, Z.I.Botev, Handbook of
Monte Carlo Methods, Wiley, 2011. （伏見正則，
逆瀬川浩孝監訳，「モンテカルロ法ハンドブッ
ク」，朝倉出版，2014）
［３］津田孝夫，モンテカルロ法とシミュレーショ
ン＜改訂版＞，培風館，1977.
19
© Copyright 2024 ExpyDoc
