炉物理プログラム演習：乱数の発生と確率分布の計算
ある変数を考えたとき、それがとりうる値に対して確率が与えられている場合には、その変数を確率変数という。可算集合 {x1 , x2 , ...} の中の値をとる確率変数は離散型と言われ、例えばサイコロの出目が挙げられるであろう。一方、確率変数 X のとる値が関数 f (x) によって∫ b
P (a ≤ X ≤ b) =f (x)dx(1)aのように表される場合、X は連続型と言われる。
計算機上でこのような確率変数を扱う場合、一般的には擬似乱数と呼ばれるものを利用する 1 。計算機での乱数の発生法には様々なものがあるが、ここでは最も基本的と言えるであろう C 言語の RAND 関数を利用する。RAND 関数は[0,RAND MAX] の範囲の整数値を返すので、[0,1] の範囲の実数の乱数を必要とするときには、RAND 関数の戻り値をRAND MAX で除してやればよいであろう。
RAND 関数は標準ライブラリ stdlib.h をインクルードすることにより利用でき、rand() もしくは random() メソッドにより [0,RAND MAX] の範囲の整数値が返される。以下に RAND 関数を用いた乱数の発生例を示す。

Listing 1: RAND 関数を用いた乱数の発生例
123456789101112131415#include < s t d l i b . h>#include <i o s t r e a m >using namespace s t d ;i n t main ( ){srand ( 1 0 ) ;cout<<RAND MAX<<” \n” ;cout<<rand()<<” \n” ;cout<<random()<<” \n” ;};return 0 ;2 行目は出力のための cout メソッドを利用するためのライブラリのインクルードを示し、4 行目は std::の記法を省略するためのものである。また、8 行目の srand メソッドは発生させる乱数列を変えるためのものであり、この引数を変えなければ実行のたびに同一の結果が得られ、これを変更することにより異なる結果を得ることができる（試してみるとよい）。
現在、非常に使い勝手が良く、周期の長い擬似乱数発生機能が C++に複数実装されているようであるが、今回は演習が目的なので、この RAND 関数を利用する課題を提示する。
はじめに、乱数の発生方法の演習として、以下の課題に取り組んでもらいたい。


問題１：乱数を用いて円周率を計算する。これは一辺の長さが 1 の正方形と直径が 1 の円の面積の比が 1:π/4 であることを利用する。すなわち、直径が 1 の円が内接する一辺の長さが 1 の正方形内にランダムに位置を決め、それが円内に存在する確率を計算すればよい。乱数の発生回数と評価された円周率の関係を求め、乱数の発生回数の増加に伴い、円周率の推定値が厳密値に近づいていくことを確認せよ。


次に、乱数を用いてある確率分布に従う確率変数の標本群を求め、それがどの程度、想定している確率変数を再現するか評価してみよう。


問題２：サイコロの出目の頻度分布を作成せよ。試行回数を 10、100、1,000、10,000 と増やしていったときに、得られた頻度分布が厳密な確率分布に一致していく様子を観察せよ。

今度は、複数のサイコロを同時に振ったときの、その出目の平均値の確率分布を考えよう。サイコロが１個のときは一様分布になることが容易に想像できるが、同時に振るサイコロの数が増えるに従って、この確率分布はどのように変化していくであろうか。


問題３：複数のサイコロの出目の平均値の頻度分布を作成せよ。試行回数は 10,000 回とする。なお、同時に振るサイコロの個数を 2 から 5 まで増加させ、頻度分布がどのように変化していくか観察せよ。
1 ここで「擬似」としているのは、これが厳密に「乱れた（ランダムな）」変数ではなく、超長期的にはある種の周期性を持つことが理由である。

1

工学の問題でしばしば現れるのは、確率変数が正規分布に従うというものである。正規分布に従う確率変数についてのランダムな標本を得たい場合、擬似乱数は [0,1] の一様分布で発生するので、これを正規分布のものに変換する必要がある。このあたりの考え方については文献 [2] に詳しく記述されており、ここではその概要を述べる。
確率変数 X が [a,b] の範囲で非ゼロの値をとる確率密度関数 f (x) に従うとしたとき、X が [a,x] の範囲の値をとる確率 F (x) は以下のように定義できる。

∫ x
F (x) =f (x′ )dx′(2)a

この F (x) を累積分布関数と呼ぶ。f (x) は確率密度を示す関数であるが、F (x) は確率そのものを示す関数となる。定義から明らかなように、ここで示した例では F (a) = 0、F (b) = 1 であり F (x) は単調増加関数となる。
このように確率密度関数 f (x) が与えられたときにこれに従う確率変数（の標本）を乱数によりランダムに決定（抽出）する場合は、以下の方法により行う。まずは [0,1] の一様分布に従う確率変数からランダムに値を抽出する（これを ξ とする）。そして、以下の式から、f (x) に従う確率変数 X の標本 x̃ を決める。

ξ = F (x̃)(3)

x̃ = F −1 (ξ)(4)

この式は、F の逆関数 F −1 を用いてと書ける。
一様分布に従う乱数を用いた正規分布に従う乱数の発生方法としては Box-Muller 法が広く知られている。これは以下の定理に基づく。


確率変数 X 及び Y が互いに独立で、ともに (0,1) 上での一様分布に従うものとする。このとき

Z1Z2√−2 log X cos 2πY,√=−2 log X sin 2πY=(5)(6)で定義される Z1 、Z2 は、平均 0、分散 1 の標準正規分布に従う互いに独立な確率変数となる。


2

任意の平均 µ、分散 σ の正規分布に従う確率変数の標本 z̃ を得たいときには、上述の方法で得た標準正規分布に従う確率変数の標本 z をz̃ = zσ + µ(7)のように変換すればよい。


問題４：平均 5.38、標準偏差 0.55 の正規分布に従う確率変数を考える。乱数を用いてこの確率変数の標本を複数個作成し、それらに基づいて平均値、分散、歪度、尖度を計算せよ。そして、標本数をパラメータとして、これら統計量推定値を図示し、標本数が増加するにつれてこれら統計量がもとの確率分布（母分布）のものに一致することを示せ。

なお、上記の統計量推定値は、I 個の標本 xi を用いて以下のように計算される。

x̄=I
∑xi ,

(8)

i=1I
1 ∑2(xi − x̄) ,I − 1 i=1)3I (∑xi − x̄I,(I − 1)(I − 2)σV= σ2 =Sk=(9)

(10)

I=1

Ku=I4
∑
3(I − 1)2
I(I + 1)(xi − x̄)−(I − 1)(I − 2)(I − 3) i=1
σ4(I − 2)(I − 3)(11)

ここで、x̄、V 、σ 、Sk、Ku はそれぞれ平均、分散、標準偏差、歪度そして尖度の推定値に対応する。なお、ここで「推定値」としているのは、これらは平均 5.38、標準偏差 0.55 の正規分布（従って、歪度と尖度はともにゼロである）の母分布から取り出した有限個の標本に基づいて推定した統計量であるためである。従って、これらの統計量推定値は母分布の統計量とは一致しない（採取した標本数が無限大であれば一致する）。

2なお、歪度、尖度については、標本分散の計算式に基づいた以下のような誤った式を用いてしまう場合があるので注意が必要である。

1 ∑(I − 1)
ISk=(=)3,(12)I
41 ∑ (xi − x̄)−3(I − 1) i=1
σ4(13)

I=1Kuxi − x̄σ問題４において、この誤った式を用いた場合の結果も示すとよいであろう。
乱数を用いて抽出された確率変数の標本から推定した平均値や分散といった統計量は、それ自身が確率変数であり、それらに対しても平均、分散等の統計量が定義される。そこで、次は、標本から推定した平均値（標本平均）と分散（標本不偏分散）の不確かさ（分散）について考えよう 2 。
ここで、確率変数 X 、Y について、その線形結合 aX + bY の分散 V (aX + bY ) は以下のように計算されることを記載しておく。

V (aX + bY ) = a2 V (X) + b2 V (Y )(14)

なお、V (X)、V (Y ) はそれぞれ X 、Y の分散を示す。
確率変数 X の I 個の標本 xi から得られる標本平均 x̄ =I∑

xi /I の分散 V (x̄) は以下のように求められる。

i=1( I)II∑
1∑
1∑
1V (x̄) = V
xi /I = 2V (xi ) = 2V (X) = V (X)II
I
i=1
i=1i=1(15)また、I 個の標本 xi から得られる標本不偏分散 s2 については、確率変数 X が正規分布に従うときには、Z = (I − 1)I2∑
s2(xi − x̄)=V (X)V (X)
i=1(16)で定義される確率変数 Z が自由度 I − 1 のカイ二乗分布に従うことを利用すれば求めることが出来る。自由度 I − 1 のカV (X)
イ二乗分布の平均は I − 1、分散は 2(I − 1) になることから、s2 =Z より、標本不偏分散 s2 の平均値として V (X)
I −1
22 {V (X)}が、分散としてが、それぞれ得られる 3 。
I −1
問題５：問題４で求めた平均値（標本平均）、分散（標本不偏分散）について、それらの不確かさ（標準偏差）を上述の方法に基づいて評価し、それら標本から推定した統計量（問題４で求めた平均値、分散）ともとの確率分布（母分布）の統計量（X の平均値 E(X) と分散 V (X)）とのずれが、推定した統計量の不確かさと整合していることを確認せよ。この評価はいくつかの異なる標本数で行えばよい。なお、標本平均、標本不偏分散の分散を求める際にはもとの確率分布の分散 V (X) が必要となる。そこで、もとの確率分布の分散が分かっている（分散を 0.552 ）とした場合と、標本不偏分散で代用した場合とで結果を比較せよ。


問題６：着目する確率変数が正規分布に従うことが分かっていない場合は、標本不偏分散の不確かさを上述のように求めることは出来ない。このような場合は、有限個の標本に基づく標本不偏分散の計算を何度も繰り返し行うことにより、標本不偏分散の頻度分布を得て、そこから標本不偏分散の平均値（期待値）と分散を推定すればよい。問題４で求めた平均値と分散について、100 回の同様の計算を行い、その結果から各々の不確かさを評価せよ。
√
なお、標本不偏分散から計算される標準偏差 s = s2 の分散 V (s) は以下の式を用いて計算できる [4]。

(2V (s ) =ds2ds)22

V (s) = (2s) V (s) = 4s2 V (s)2 ここでの記述は文献 [3, 4] に拠る。

(3 V (s2 ) = V

V (X)ZI −1)=

(V (X))22 (V (X))2V (Z) =2(I − 1)
I −1

3(17)

複数の確率変数についてランダムに標本を抽出する場合、それぞれが独立、すなわち相関係数がゼロであるならば、それぞれの確率変数について上記の方法で標本を計算すればよい。一方、確率変数間に相関がある場合には、それに対する配慮が必要である 4 。

L 個の確率変数 Fl (l = 1, ..., L) を考え、それぞれについてその確率分布に基づいて I 個の標本を得たとし、Fl の i 番目の標本を fli と書く。Fl と Fl′ の標本共分散 mll′ は以下のように定義できる。
I∑(fli − f¯lmll′ = i=1)( i)fl′ − f¯l′(18)I −1

ここで f¯l は Fl の標本平均を示す。
この複数の確率変数をまとめてベクトルで f と表記し、i 番目の標本及び標本平均をそれぞれ f i 、f̄ と表記すると、fの共分散行列 M は以下のように記述できる。

()()T
M = F − F̄ F − F̄
ここで、F、F̄ はいずれも L × I の行列で、

(F =

f1···

f2fI)

(F̄ =f̄f̄···)f̄

f11

f12······

fL1

f22..
.
fL2

 1
 f2=
 ..
 .

(19)

f¯1
 ¯
 f2=
 ..
 .
f¯L

f¯1
f¯2..
.
¯
fL············

f1I

f2I .. ,. (20)

fLI


f¯1
f¯2 .. ,. f¯L(21)である。
さて、次に f の線形変換 f ′ を以下のように定義する。

f ′ = Pf
式 (19) の両辺について、左から P、右から PT を乗じると以下の式を得る。
()()T
PMPT = PF − PF̄ F − F̄ PT()( ())T= PF − PF̄ P F − F̄()()T= PF − PF̄ PF − PF̄()()T= F′ − F̄′ F′ − F̄′ = M′

(22)

(23)

つまり、行列 PMPT は f ′ の共分散行列 M′ となることが分かる。
ここで、線形変換行列 P を適切に選ぶことにより、M′ を対角行列にすることが出来れば、すなわちパラメータ f ′ がそれぞれ独立関係にあるように出来れば、各々のパラメータ fl′ について前述の方法で標本を求め、以下の式で f に変換すればよいことが分かる。

f = P−1 f ′(24)次に、行列 M′ を対角行列にするための方法として、ここでは特異値分解を用いた方法を説明する。
行列 M に特異値分解を施すと以下の式が得られる。

M = UDVT = UDUT

(25)

ここで、M は対称行列であることから、U = V を用いている。また、D は対角行列である。U はユニタリ行列である
ため、UUT = UT U = I より、以下の式が得られる。

UT MU = D

(26)

4 例えば、２つの確率変数 A、B の間に強い正の相関がある場合、A の標本としてその平均値よりも大きい値が得られた場合、B の標本としてもその平均値より大きい値が得られるべきであると考えられる。

4

この式と式 (23) を比較することにより、線形変換行列として P = UT を用いることにより、それぞれのパラメータが独立となる確率変数 f ′ を定義できることが分かる。f ′ の標本が得られた後は、f = Uf ′ により、もとの確率変数に戻せばよい。
以上の記述に基づいた、複数の確率変数についてランダムに標本を抽出する方法は以下のようになるであろう。

1. 共分散行列 M の対角化を行い、行列 U を得る。
2. f を線形変換した f ′ について、その平均値 f̄ ′ と共分散行列 D に基づいて標本を抽出する。ここで、f ′ はそれぞれが独立な確率変数であるので、それぞれ独立に標本抽出が可能である。
3. 得られた f ′ についての標本から、目的とする確率変数に関する標本を f = Uf ′ の関係より求める。
なお、一般的には以下のような方法が採用されている。共分散行列の対角化の式を以下のように記述する。
()T()T
M = UDUT = UD1/2 D1/2 UT = UD1/2 D1/2 UT = UD1/2 UD1/2
= AAT

(27)

( )−1
()T
ここで、A−1 = D−1/2 UT 、 AT
= UD−1/2 = A−1 なので、()TA−1 M A−1 = I(28)が成り立つ。これより、f ′ = A−1 f なる変換を行った場合、f ′ の共分散行列は I となることが分かる。従って、標準正規分布に従う確率変数の標本を f ′ としたときに、それに対応する元の確率変数の標本は以下の式で求められる。
()
f = Af ′ = UD1/2 f ′(29)

この方法では、前述の方法のように線形変換後の確率変数の平均値 f̄ ′ を求めておく必要がないなどの利点がある。


問題７：平均 5.38、標準偏差 0.55 の正規分布に従う 2 つの確率変数を考える。両者の相関を 0.8 として 100、1,000、10,000 個の標本を作成し、それら標本群から共分散行列を求め、標本数の増加に伴いもとの（母分布の）共分散行列の再現性が向上することを確認せよ a 。
a なお、標準偏差がいずれも σ であり相関がない２つの確率変数 x、y について、x′ = x、y ′ = Cx +√1 − C 2 y という変換を行って得られる x′ 、y ′ を考えると、これらは標準偏差が σ で相関が C となる。これについては実際に線形変換行列 P を(P=1C√ 01 − C2)(30)として、x′ 、y ′ の共分散行列を計算し、そのような結果が得られること確認してみるとよい。問題６はこのようにしても解くことができるが、ここでは解説で述べられているように、もとの共分散行列の固有ベクトルから計算する方法を用いてもらいたい。



共分散行列とは、対角成分が分散、非対角成分が共分散に対応する行列である。共分散行列 M の要素 Mii は fi の分散になる。また fi と fj の相関係数 Cij は共分散行列を用いて以下のように定義される。

MijMij√Cij = √=σi σjMii Mjj(31)

ここで σi は fi の標準偏差を示す。
また共分散行列のような対称正方行列について特異値分解を行う場合、それはこの行列の固有値を求める操作と等価となる。すなわち、共分散行列 M の固有値を λi 、固有関数を xi とすると以下の方程式が成り立つ。

Mxi = λi xi(32)

ここで、X = (x1 x2 · · · xL )、D = diag (λ1 , λ2 , · · · , λL ) とすると

MX = XD(33)

XT MX = D(34)が得られる。X−1 = XT よりとなり、線形変換により対角成分がゼロの共分散行列が得られる。問題７では行列のサイズが小さいため固有値と固有ベクトルを解析的に得ることができることに留意すること。
最後に、確率変数の対数が正規分布に従う対数正規分布の場合について、標本の抽出方法を以下にまとめておく。こ
こで、確率変数 x の対数 z = ln(x) が正規分布に従うものとする。

5

1. x の平均値 µx と共分散 Vx から、z の平均値 µz と共分散 Vz を計算する。
2. z について、µz 、Vz に基づいて標本の抽出を行う。
3. x = exp(z) より、z の標本に対応する x の標本を得る。

参考文献[1] 東京大学教養学部統計学教室編、「基礎統計学 I、統計学入門」、東京大学出版会.
[2] 長家康展、「モンテカルロ計算の基礎理論及び実験解析への適用」、第 38 回炉物理夏期セミナーテキスト、(2006).
[3] 山本章夫、「不確かさ評価の基礎」、第 44 回炉物理夏期セミナーテキスト、(2012).
[4] 遠藤知弘、「ランダムサンプリング法で得られた不確かさの統計誤差」、私信.

6