本日の講義資料

1
2015年6月4日
エネルギー工学演習II
乱数とモンテカルロ
― 第１回：乱数の発生 ―
渡辺幸信
2
エネルギー工学演習IIのスケジュール
１．コンピュータ演習－１
２．コンピュータ演習－２
３．コンピュータ演習－３
４．コンピュータ演習－４
５．コンピュータ演習－５
６．コンピュータ演習－６
７．コンピュータ演習－７
８．コンピュータ演習－８
９．英作文－１
7/2
１０．英作文－２
7/9
１１．コンピュータ演習－９
１２．コンピュータ演習－１０
１３．英作文－３
7/16
１４．英作文－４
7/23
4/16
4/23
4/30
5/14
5/21
5/28
6/4
6/11
山本
山本
6/18
6/25
山本
山本
山口
山口
松清
山口
山口
山口
渡辺
渡辺
渡辺
渡辺
変更
本演習のＨＰ
http://enep.ence.kyushu-u.ac.jp/lectures06.html
3
１．はじめに
• モンテカルロ法（Monte Carlo method）は乱数
（random number）を用いた数値解析技法の総称。
• カジノで有名なモナコの地名モンテカルロに由来し
て命名された。
• 乱数とはある指定された確率分布をもつ数列を意
味する。
4
モンテカルロ法の対象
• モンテカルロ法の対象となる問題は、確率的な事象とそ
うでない決定論的事象に大別できる。
• 確率的な事象は広範囲の問題に含まれており、モンテ
カルロ法は自然科学、工学、医学、経済・社会科学等の
様々分野におけるコンピュータ・シミュレーションで広く
用いられている。
– 例えば、原子炉内のランダムな中性子の運動を、その素過程
である原子核衝突の断面積（反応確率を意味する）を用いて乱
数で表現し、炉内の中性子挙動を調べること（中性子拡散現象
のシミュレーション）に使われる。
– 決定論的事象への典型的な適用例として多次元空間での定
積分の数値計算がある。
5
PHITS計算結果の例
137Cs
から放出された100,000個の光子の挙動を模擬
個々の放射線挙動を模擬することにより，全体的な挙動（平均値）を導出
6
JQMD (JAERI Quantum Molecular Dynamics) モデル
• 原子核を核子の集合体と仮定して，全ての核子間力を数値解析で解く手法
• 入射放射線の核種・エネルギーとターゲット核種の情報から核反応の終状態
を予測することができる
Iwase et al
Production
yields
of residual
nuclides
Double
differential
cross
section of
neutron
56Fe
K. Niita et al., Phys. Rev. C52 (1995) 2620
800 MeV/u on 208Pb
7
中性子発生装置周辺の電離量計算
Energy
Deposition
tt == 201000
410nsec
100
nsec
nsec
nsec
Be reflector
Fe reflector
Moderators
Proton
Hg target
Moderators
水銀ターゲット付近の電離量計算
Hg target
水銀ターゲット付近の幾何形状
M. Harada et al. J. Nucl. Material 343, 197 (2005)
8
本演習の内容
• 組み込み関数を使った一様乱数の発生（6/4）
• 一様乱数を用いたサンプリング手法（6/4）
• 数値積分への応用（6/11）
– pの計算
– モンテカルロ法による定積分
9
２．乱数
• さいころやルーレットは簡単な乱数発生器の１つ
• 乱数表や物理乱数源（ランダム物理現象を利用した電子管
の雑音、放射性物質の崩壊など）
• 一般にモンテカルロ計算では算術乱数が使われる。この乱
数を発生させる方法が、計算機内の簡単な四則演算で行わ
れるので“算術的”であり、ランダムでない方法で近似的にラ
ンダムな数列を作るので、“疑似”乱数と呼ばれる。
乱数の中でも基本となるのは、
[0,1]の区間に一様に分布する一様乱数
10
演習１： [0,1]の一様乱数発生
11
乱数発生用組み込みサブルーチンrandom_numberの使い方
10個の一様乱数を発生して画面上に出力
program random_uni
implicit none
real r
integer i
do i = 1, 10
call random_number(r)
write(*,*) r
enddo
end program random_uni
Call 文でサブルーチンを
呼び出す。結果は引数rに入る。
一様乱数の出力結果
3.9208680E-07
2.5480442E-02
0.3525161
0.6669145
0.9630555
0.8382882
0.3353550
0.9153272
0.7958636
0.8326931
何度実行しても同じ結果
が出力される。
毎回、違う乱数列が
ほしい場合
乱数の初期値設定サブ
ルーチン random_seed
（毎回、異なった乱数列を
発生できる）
何度か実行して結果が変
化することを確認すべし
program random_uni
implicit none
real r
integer i
call random_seed()
do i = 1, 10
call random_number(r)
write(*,*) r
enddo
end program random_uni
12
出席チェックタイム
下記のアドレスへメール発信してください。
[email protected]
件名：学生番号と氏名の併記
1TE11**** 九州太郎
講義の印象、感想、要望等のコメントがあれば、
記入してください。
13
14
一様乱数の頻度分布作成
• 頻度分布をヒストグラム（柱状グラ
フ）で表すプログラム（hist_random）
を用いて乱数の一様性をテストする
• ここでは[0,1]を１０個のビン（すなわ
ち、0～0.1, 0.1～0.2, …., 0.9～1.0）
に分割する。発生した乱数がどの
ビンに入るか演算して、そのビンの
個数を１つ増やしていく。
• 最終的に同じビンに入った個数を
集計してヒストグラム表示でグラフ
化する。
０
１０個のビン
１
本演習のＨＰ
http://enep.ence.kyushu-u.ac.jp/lectures06.html
15
一様乱数の頻度分布作成
program hist_random
implicit none
real r, dr, rmax
integer nb, nmax, ir, i, hist(0:10)
dr=0.1
! ヒストグラムのビン幅
rmax=1.0
! ヒストグラム横軸の最大値（上限値）
nb=rmax/dr
! ヒストグラムのビン数 = 上限値(rmax)/ビン幅(dr)
nmax = 10000
! 乱数の発生数
do i =0,nb
! ヒストグラムデータの初期化注）do-loopを使う代わりに、
hist(i)=0
! hist(:)=0 と一行に書いて一括代入も可
enddo
call random_seed() ! 乱数の初期値設定（毎回、異なった乱数列を発生可）
do i= 1, nmax
call random_number(r) ! 乱数１個発生
ir= r / dr
! ヒストグラムのビン番号irの計算
hist(ir)=hist(ir)+1
! ヒストグラムのビンに一つ加える
enddo
do i=0,nb
write(*,*) i,hist(i)
! ヒストグラムの出力
enddo
end program hist_random
16
リダイレクション機能
集計結果をグラフ化するために、リダイレクション機能を使っ
て結果をファイル（randhist.dはファイル名）に保存する。
./a.out > randhist.d
write(*,*) : 標準出力（スクリーン上）
出力先をファイルに変更する方法
open文を使わずにファイル出力可（便利）
注） a.outはFortranコンパイラifort で作成した実行形式ファイ
ルの名称。-oオプションをつけないと, a.outが生成される。
gnuplotを用いてヒストグラム表示でグラフ化
gnuplot > plot [] [0:] ‘randhist.d’ with step
ｘ軸はそのまま、y軸の始まりを0からに変更
17
一様乱数の頻度分布の例
nmaxを変化させて計算した結果をレポートには示して、考察すること。
nmax=100
nmax=10000
nmax=1000
nmax=100000
18
一様でない乱数の発生方法
ー直接法－
19
p(x)dx : 確率変数 x が [x,x+dx] にある確率とする。
p(x)は確率密度関数、区間[a,b]で規格化されているとする。

b
a
確率分布関数
（累積分布関数）
p( )d  1
p(x)
x
P( x)   p( )d  r
a
1
x  P (r )
a
b x
1
P(x)
r
単調増加関数
[0,1]の一様乱数 r より x が一意に決まる。
0
a
x
b
一様でない乱数の発生方法
ー指数関数分布－
確率密度関数
1
 exp(  x /  ) ( x  0)
p( x)   

0
( x  0)

(3-3)に代入して、確率分布関数 P(x) を求める。
x
1
0

r  P( x)  
exp(  )d  1  exp( x)
ここで逆変換(3-4)を行って
x   ln(1  r )
[0,1]の一様乱数 r から指数関数の確率分布に従う確率変数 x を
サンプリングすることができる。
20
21
演習課題（２）
演習課題(1)のプログラムを参考にして、上記の指数関数分布
（λ=1 とせよ）に対するサンプリングのプログラムを作成せよ。
•
•
•
•
nmax=10000
ビン幅は0.1
x軸の上限を10.0
縦軸 linear表示とlog表示の頻度分布
ヒストグラム
自然対数（底がe）
y=log(x)
ヒント： gnuplotのコマンド
set logscale y
22
演習課題（３）発展問題
それぞれサイコロの目が出る確率が表１で与えられる”サイバー
サイコロ”をコンピュータ上で製作し、試行回数10, 100, 1000とした
場合の各目が出る頻度分布（ヒストグラム）を作成せよ。
(3-3) 式
サイの目i
1
2
3
4
5
6
確率pi
0.1
0.2
0.3
0.1
0.2
0.1
1
x
P( x)   p( )d  r
a
0.8
r=0.75
0.6
i
Pi
Pi   pk  r
k 0
ヒント： if文を使うことになる。
0.4
0.2
0
1
2
3
4
5
サイコロの目
6
レポート提出の注意点
23
• 電子ファイル形式のレポートをメール添付で提出する際、空メールで
なく、教員に対する通常のメール（手紙）の書き方をすること。（宛先、
提出する旨の簡単なメール文章、発信者の氏名等）
• Gnuplotで作成した図ファイル（PNG形式）はワードの挿入機能（各自
で使い方を調べる）を使ってワードファイル内に図として取り込むこと
ができるので、レポートは図も含めた１つのファイルにすること。
• レポートファイル名：学生番号＋氏名
• レポートをファイルとして提出する際も、通常の紙ベースのレポートと
同じく、氏名や学生番号、提出日等の情報をレポートの始めに記載
するのが望ましい。
• レポートの書き方の参考書、例えば
木下是雄「理系のための作文技術」（中公新書）
を読んで、執筆スキル（技巧）を磨くこと。（講義HP参照）

Download Report