乱数とシミュレーション(simulation with random numbers)とは、一様で独立な数の列である乱数から試行を作り、くり返して確率や面積を見積もる方法である。正方形に $n$ 点を打って四分円の内側の割合を 4 倍すると、平均 $\pi$、標準偏差 $\sqrt{\pi(4-\pi)}/\sqrt n\approx1.642/\sqrt n$ の推定値になり、誤差は $1/\sqrt n$ の速さでしか減らない(誤差 $0.01$ を約 95% の確率で達成するには約 10 万点)。1 次元の面積なら区分求積のほうが速いが、乱数は数えにくい確率や変数の多い積分に向く。計算機の擬似乱数は線形合同法などの決まった計算で作られ、周期は法 $m$ 以下である。周期が短い、並びに規則が残る、点が一様に散らばらないときは、結果が偏る。
前提知識: 大数の法則とさいころの平均, 標本平均の分散, 中心極限定理(高校数学), 確率密度関数と連続型確率変数
数学 B では、さいころ・乱数表・計算機の 乱数 を使って試行をくり返し、起こった割合から確率を見積もる方法を学ぶ。数えるのが難しい確率でも、試行を何度もくり返せば、相対度数は確率に近づく(大数の法則とさいころの平均)。このような実験を シミュレーション という。まず、計算機で実際に試してみる。
さいころを 10 回振るとき、同じ目が 3 回以上続けて出る(たとえば $\cdots4,4,4\cdots$)確率を $P$ とする。10 回の目の並び $6^{10}=60466176$ 通りを数えるのは大変なので、計算機の乱数で「10 回振る」を $N$ 回くり返し、3 回以上続いた回の割合を調べた(乱数の種を固定して作った)。
| くり返しの回数 $N$ | $100$ | $1000$ | $10000$ | $100000$ |
|---|---|---|---|---|
| 3 回以上続いた回数 | $19$ | $172$ | $1798$ | $17986$ |
| 割合 | $0.1900$ | $0.1720$ | $0.1798$ | $0.1799$ |
「直前の何回が同じ目か」で場合分けする漸化式を使って計算機で正確に計算すると、$P=\frac{1827071}{10077696}=0.18130\ldots$ である。割合は $N$ を増やすと $0.18$ のあたりに落ち着くが、$N=100000$ でも小数第 3 位はまだずれている。
$0$ 以上 $1$ 未満の乱数を 2 つ引いて点 $(x,y)$ を作ると、点は一辺 $1$ の正方形の中のでたらめな点になる。点が四分円 $x^2+y^2\le1$ の内側に入る確率は、四分円の面積 $\frac\pi4$ に等しい。そこで、$n$ 点のうち内側に入った点の数を $m$ として、$\pi$ を $4\times\frac mn$ で見積もる。
計算機で 1000 点を打ち、最初の 100 点だけを数えると $m=71$ で $4\times\frac{71}{100}=2.84$、1000 点全部では $m=810$ で $4\times\frac{810}{1000}=3.24$ であった(図 1、図 2。乱数の種を固定して作った)。どちらも $\pi=3.14159\ldots$ から $0.1$ 以上ずれている。
一辺 1 の正方形に打った 100 点。四分円の内側の点を青、外側の点を橙で示す。
同じ乱数で打った 1000 点。内側の点の割合が四分円の面積 π/4 に近いことを見る図。
1000 点も打ったのに、円周率が 1 桁しか合わない。ここで次の問いが出てくる。
| 高校の計算 | この記事での見方 | ボックス |
|---|---|---|
| 乱数さいころ・乱数表 | 独立で一様な乱数の列 | def-rns-random |
| 相対度数で確率を見積もる | 標本平均。平均は確率、標準偏差は $\frac{\sqrt{p(1-p)}}{\sqrt n}$ | thm-rns-main |
| 標本平均の分散 $\frac{\sigma^2}n$ | 誤差は $\frac1{\sqrt n}$ の速さでしか減らない | thm-rns-main、cor-rns-trials |
| 余りの計算・周期 | 擬似乱数の作り方(線形合同法) | def-rns-lcg、prop-rns-period |
乱数さいころ(正二十面体に $0$ から $9$ を 2 回ずつ書いたもの)を振ると (i) が得られる。乱数表は、そのような数字を並べた表である。計算機の乱数は (ii) を目指して作られるが、実際には決まった計算で作る数の列であり、擬似乱数 という(節「擬似乱数の作り方:線形合同法」)。GS06 §1.1・§2.1 も、計算機の乱数を「区間に入る確率が区間の長さに等しい」ように選ばれる数として導入し、実際には決まった手順で作られる数であることを注意している。
一様乱数から、いろいろな確率の試行を作れる。
$U$ を区間 $[0,1)$ の一様乱数とする。
計算機が $U=0.7312$ を返したとする。
ex-rns-pi-intro の見積もりを、確率変数として調べる。点を打つたびに、内側なら $1$、外側なら $0$ となる量を考えるのがこつである。
$(X_1,Y_1),(X_2,Y_2),\dots,(X_n,Y_n)$ を、$2n$ 個の独立な区間 $[0,1)$ の一様乱数から作った点とし、$X_i^2+Y_i^2\le1$ なら $I_i=1$、そうでなければ $I_i=0$ とする。$\hat p=\frac{I_1+\cdots+I_n}{n}$、$\hat\pi=4\hat p$ とすると
$$
E(\hat\pi)=\pi,\qquad \hat\pi\ \text{の標準偏差}=\frac{\sqrt{\pi(4-\pi)}}{\sqrt n}=\frac{1.6421\ldots}{\sqrt n}
$$
である。したがって、誤差の目安である標準偏差を $\frac1{10}$ にするには、$n$ を $100$ 倍にしなければならない。
方針:$I_i$ は確率 $p=\frac\pi4$ で $1$ をとる量なので、平均と分散が $p$ と $p(1-p)$ になる。$\hat p$ はその標本平均なので、標本平均の分散 の主定理を使う。
段 1($I_i$ の平均と分散)。点 $(X_i,Y_i)$ が四分円に入る確率は、四分円の面積 $p=\frac\pi4$ である(確率密度関数と連続型確率変数 と同じく、正方形の中の図形に入る確率はその面積に等しいとする)。$I_i$ は確率 $p$ で $1$、確率 $1-p$ で $0$ なので
$$
E(I_i)=1\cdot p+0\cdot(1-p)=p,\qquad E(I_i^2)=1^2\cdot p+0^2\cdot(1-p)=p,\qquad V(I_i)=p-p^2=p(1-p)
$$
である。
段 2(独立性)。$I_i$ は $X_i,Y_i$ だけから決まり、$2n$ 個の乱数は独立なので、$I_1,\dots,I_n$ は独立である(本記事では、独立な確率変数を別々に使って作った確率変数どうしが独立になることを認めて使う)。
段 3(標本平均)。標本平均の分散 の主定理から $E(\hat p)=p$、$V(\hat p)=\frac{p(1-p)}n$ である。
段 4(4 倍する)。$\hat\pi=4\hat p$ なので、1 次式の変換の公式(確率変数の期待値と分散)から $E(\hat\pi)=4p=\pi$、$V(\hat\pi)=16\cdot\frac{p(1-p)}n$ である。$16p(1-p)=16\cdot\frac\pi4\cdot\frac{4-\pi}4=\pi(4-\pi)$ なので、標準偏差は $\frac{\sqrt{\pi(4-\pi)}}{\sqrt n}$ である。$\sqrt{\pi(4-\pi)}=\sqrt{2.6967\ldots}=1.6421\ldots$ である。
段 5(100 倍)。標準偏差は $\frac1{\sqrt n}$ に比例するので、$n$ を $100$ 倍にすると $\frac1{\sqrt{100}}=\frac1{10}$ 倍になる。$\square$
thm-rns-main の標準偏差は、$n=100$ で $0.164$、$n=10000$ で $0.0164$、$n=1000000$ で $0.00164$ である。計算機で(ex-rns-pi-intro とは別の乱数の種で)実際に打った結果は次のとおりであった。
| 点の数 $n$ | $100$ | $10000$ | $1000000$ |
|---|---|---|---|
| 内側の点の数 $m$ | $78$ | $7895$ | $785658$ |
| $\hat\pi=\frac{4m}{n}$ | $3.12000$ | $3.15800$ | $3.14263$ |
| 誤差 $\hat\pi-\pi$ | $-0.0216$ | $+0.0164$ | $+0.0010$ |
| 標準偏差 | $0.164$ | $0.0164$ | $0.00164$ |
$n=10000$ の誤差は標準偏差とほぼ同じ大きさで、$100$ 万点打っても誤差は $0.001$ 程度で、小数第 2 位までしか合わない。ex-rns-pi-intro の $2.84$($n=100$)と $3.24$($n=1000$)は、どちらも標準偏差の約 $2$ 倍ずれた場合で、珍しいが起こりうる大きさである。
乱数による π の推定の誤差(点)と理論値(実線)、右端の区分求積の誤差(破線)を、横軸・縦軸とも対数目盛で描いた図。乱数の誤差は傾き −1/2、区分求積は傾き −1 の直線に沿って減ることを見る。
図 3 は、各 $n$ で推定を $200$ 回くり返し、誤差の大きさ(誤差の 2 乗の平均の平方根)を点で描いたものである。点は理論値 $\frac{1.642}{\sqrt n}$ の直線によく乗っている。両対数のグラフで傾きが $-\frac12$ であることが、「$n$ を $100$ 倍にして誤差が $\frac1{10}$」を表している。
必要な点の数を見積もるには、中心極限定理(高校数学) を使う。
$\varepsilon>0$ とする。$n$ が大きいとき、$\lvert\hat\pi-\pi\rvert\le\varepsilon$ となる確率を約 $0.95$ にするには、およそ
$$
n\ge\Bigl(\frac{1.96\times1.6422}{\varepsilon}\Bigr)^2
$$
であればよい。これは中心極限定理による近似であり、保証ではない。確率 $0.95$ 以上を保証する数は、Chebyshev の不等式から $n\ge\frac{\pi(4-\pi)}{0.05\,\varepsilon^2}$ である。
$\sigma_n=\frac{\sqrt{\pi(4-\pi)}}{\sqrt n}$ とする。
段 1(近似)。$\hat p$ は独立で同じ分布に従う $I_i$ の標本平均なので、中心極限定理(高校数学) の標本平均の正規近似から、$n$ が大きいとき $P\bigl(\lvert\hat\pi-\pi\rvert\le1.96\sigma_n\bigr)$ は約 $0.95$ である($\lvert\hat\pi-\pi\rvert\le1.96\sigma_n$ は $\lvert\hat p-p\rvert\le1.96\frac{\sigma_n}4$ と同じで、$\frac{\sigma_n}4$ は $\hat p$ の標準偏差である)。$1.96\sigma_n\le\varepsilon$ を $n$ について解くと、$\sqrt n\ge\frac{1.96\sqrt{\pi(4-\pi)}}{\varepsilon}$ から主張の式を得る。
段 2(保証)。Chebyshev の不等式 $P(\lvert Z-E(Z)\rvert\ge\varepsilon)\le\frac{V(Z)}{\varepsilon^2}$(大数の法則とさいころの平均 で証明)を $Z=\hat\pi$ に使うと、$P(\lvert\hat\pi-\pi\rvert\ge\varepsilon)\le\frac{\pi(4-\pi)}{n\varepsilon^2}$ である。右辺が $0.05$ 以下なら、$\lvert\hat\pi-\pi\rvert<\varepsilon$ の確率は $0.95$ 以上である。$\frac{\pi(4-\pi)}{n\varepsilon^2}\le0.05$ を解くと $n\ge\frac{\pi(4-\pi)}{0.05\varepsilon^2}$ である。$\square$
$\frac\pi4$ は四分円の面積なので、$\frac\pi4=\int_0^1\sqrt{1-x^2}\,dx$ である。乱数を使わずに、この積分を区分求積で計算することもできる(区分求積と積分の定義)。
$f(x)=\sqrt{1-x^2}$ は $0\le x\le1$ で減少する。区間 $[0,1]$ を $n$ 等分し、$x_k=\frac kn$ とする。小区間 $x_{k-1}\le x\le x_k$ では $f(x_k)\le f(x)\le f(x_{k-1})$ なので、各辺を積分して足すと
$$
R_n:=\frac1n\sum_{k=1}^{n}f(x_k)\le\frac\pi4\le\frac1n\sum_{k=1}^{n}f(x_{k-1})=:L_n
$$
である(等号は成り立たない)。$L_n-R_n=\frac{f(0)-f(1)}n=\frac1n$ である。
$n=100$ で計算すると(計算機による)、$4R_{100}=3.12042$、$4L_{100}=3.16042$ で、$3.1204<\pi<3.1605$ が 確実に 言える。はさむ区間の幅は $\frac4n=0.04$ である。同じ $100$ 個の点を乱数で使うと、ex-rns-sd のとおり標準偏差が $0.164$ もあり、しかも「確実に」とは言えない。
| 区分求積(右端・左端) | 乱数による推定 | |
|---|---|---|
| 使う点 | 等間隔の $n$ 個 | でたらめな $n$ 個 |
| 誤差の減り方 | $\frac1n$ に比例 | $\frac1{\sqrt n}$ に比例(確率的) |
| $n=100$ での誤差の大きさ | $0.04$ 以下(はさんで保証) | 標準偏差 $0.164$ |
| 誤差の保証 | 不等式で確実 | 確率でしか言えない |
| 次元が上がると | 格子の点が $n^d$ 個必要 | 誤差の式は次元 $d$ によらない |
図 3 の破線は、右端の区分求積 $4R_n$ の誤差である。傾き $-1$ の直線に沿って、乱数の誤差(傾き $-\frac12$)よりずっと速く減る。1 次元の面積なら区分求積のほうがよい。乱数が役に立つのは、正確に計算するのが難しい確率(ex-rns-runs)や、次元の高い積分(rem-rns-mc)のときである。円周率を評価する では、誤差が不等式で保証された $\pi$ の評価をまとめている。
計算機は、決まった計算で「乱数らしく見える」数の列を作る。最も簡単な方法の 1 つが次の方法である。
正の整数 $m$ と、$0$ 以上 $m$ 未満の整数 $a$、$c$、$x_0$ を決め、
$$
x_{k+1}=(a\,x_k+c\ \text{を}\ m\ \text{で割った余り})\qquad(k=0,1,2,\dots)
$$
で数の列 $x_0,x_1,x_2,\dots$ を作る。$u_k=\frac{x_k}{m}$ を $[0,1)$ の擬似乱数として使う。$x_0$ を 種(シード)という。
同じ種からは、いつも同じ列ができる。シミュレーションの結果を他の人が確かめられるように、使った種を記録しておく。
$m=16$、$a=5$、$c=3$、$x_0=1$ とすると
$$
x_1=5\cdot1+3=8,\quad x_2=(5\cdot8+3)-32=11,\quad x_3=(5\cdot11+3)-48=10,\ \dots
$$
と続き、列は
$$
1,\ 8,\ 11,\ 10,\ 5,\ 12,\ 15,\ 14,\ 9,\ 0,\ 3,\ 2,\ 13,\ 4,\ 7,\ 6,\ 1,\ 8,\ \dots
$$
となる。$x_{16}=1=x_0$ で、$0$ から $15$ までの $16$ 個の数がちょうど 1 回ずつ現れてから、同じ並びをくり返す。
def-rns-lcg の列には、$0\le i< j\le m$ で $x_i=x_j$ となる $i$、$j$ がある。このとき $k\ge i$ のすべての $k$ で $x_{k+(j-i)}=x_k$ が成り立つ。つまり列は途中から周期 $j-i$($m$ 以下)でくり返す。
段 1(同じ値が現れる)。$x_0,x_1,\dots,x_m$ は $m+1$ 個の数で、どれも $0,1,\dots,m-1$ の $m$ 通りのどれかである。鳩の巣原理 により、$0\le i< j\le m$ で $x_i=x_j$ となるものがある。
段 2(くり返す)。$x_{k+1}$ は $x_k$ だけから決まる。$x_i=x_j$ なら $x_{i+1}=x_{j+1}$、それから $x_{i+2}=x_{j+2}$、…と、帰納法で $t\ge0$ のすべてについて $x_{i+t}=x_{j+t}$ である。$k=i+t$ とおけば主張の式になる。$j-i\le m-0=m$ である。$\square$
周期は $m$ を超えられない。ex-rns-lcg16 は周期がちょうど $m$ になる場合である。$c=0$ の場合は $x_k=(a^kx_0\ \text{を}\ m\ \text{で割った余り})$ となるので、周期は「$a$ を何乗すると $m$ で割った余りが $1$ に戻るか」($a$ と $x_0$ がどちらも $m$ と互いに素のとき)で決まる。
周期が $m$ でも、よい乱数とは限らない。ex-rns-lcg16 と同じ形の $x_{k+1}=(5x_k+1)$ を $256$ で割った余り(周期 $256$)で、連続する 2 つの値の組 $(u_k,u_{k+1})$ を点として描くと、図 4 のように $5$ 本の平行な直線の上に並ぶ。$u_{k+1}$ は $5u_k+\frac1{256}$ から整数を引いた値なので、傾き $5$ の直線にしか乗れないからである。周期のとても長い擬似乱数(図 5)では、このような並びは見えない。
線形合同法 x_{k+1} = (5x_k + 1) mod 256 の連続する 2 つの値の組。256 点が傾き 5 の 5 本の直線の上に並ぶことを見る図。
周期のとても長い擬似乱数の連続する 2 つの値の組 256 点。正方形の中にばらばらに散らばることを見る図。
シミュレーションの結果が信用できるのは、使う乱数が def-rns-random の性質(一様で独立)に近く、それを正しく使っているときである。条件が崩れた例を並べる。
| 外す条件 | 反例 | 成り立たなくなること |
|---|---|---|
| 周期が十分長い | ex-rns-counter-short:周期 $4$ の線形合同法でさいころを作る | どの目も確率 $\frac16$ で出る |
| 並びが独立 | ex-rns-counter-parity:線形合同法の値の偶奇で硬貨を作る | 2 回続けて表が出る確率 $\frac14$ |
| 点が正方形全体に一様 | ex-rns-counter-diagonal:同じ乱数を $x$ にも $y$ にも使う | $E(\hat\pi)=\pi$ |
ex-rns-lcg-mult の 2($m=16$、$a=3$、$c=0$、$x_0=1$)の列 $1,3,9,11,\dots$ から、prop-rns-transform の方法でさいころの目 $\lfloor6u_k\rfloor+1$($u_k=\frac{x_k}{16}$)を作ると
$$
\Bigl\lfloor\frac{6}{16}\Bigr\rfloor+1=1,\quad\Bigl\lfloor\frac{18}{16}\Bigr\rfloor+1=2,\quad\Bigl\lfloor\frac{54}{16}\Bigr\rfloor+1=4,\quad\Bigl\lfloor\frac{66}{16}\Bigr\rfloor+1=5
$$
がくり返され、$3$ と $6$ の目は何回振っても出ない。何万回くり返しても、$3$ の目の相対度数は $0$ のままで、$\frac16$ に近づかない。
崩れている条件は、周期が十分長いことである。その結果、値が $[0,1)$ に一様に散らばらない。周期 $4$ の列は $4$ つの値しかとらないので、一様乱数の代わりにならない。
ex-rns-lcg16 の列 $1,8,11,10,5,12,\dots$ は周期 $16$ で、$0$ から $15$ をちょうど 1 回ずつとる。しかし偶奇を見ると、奇、偶、奇、偶、…と交互になっている。実際、$a=5$、$c=3$ はどちらも奇数なので、$x_k$ が奇数なら $5x_k+3$ は偶数、偶数なら奇数で、$16$ で割った余りの偶奇もそれと同じである。
そこで「$x_k$ が奇数なら表」として硬貨を作ると、表と裏が必ず交互に出る。表と裏の相対度数はちょうど $\frac12$ ずつで、1 回ずつ見れば公平に見えるが、「2 回続けて表」は決して起こらない(本当に独立な公平な硬貨なら確率 $\frac14$)。
崩れている条件は、並びの独立性である。周期が最大で値が一様に現れても、並び方に規則が残る。図 4 の直線の並びも同じ種類の現象である。
ex-rns-pi-intro で、乱数を節約しようとして、1 つの乱数 $U$ を $x$ にも $y$ にも使い、点を $(U,U)$ としたとする。点は正方形の対角線の上にしか来ない。$U^2+U^2\le1$ は $U\le\frac1{\sqrt2}$ と同じなので、内側に入る確率は $\frac1{\sqrt2}=0.7071\ldots$ で、$4\hat p$ の平均は $\frac4{\sqrt2}=2\sqrt2=2.828\ldots$ になる。点をいくら増やしても $\pi$ には近づかない。
崩れている条件は、$x$ 座標と $y$ 座標が独立であることで、そのために点が正方形全体に一様に散らばらない。thm-rns-main の証明の段 1(点が四分円に入る確率は面積 $\frac\pi4$)が成り立たなくなる。
$[0,1)$ の一様乱数 $U_1,\dots,U_n$ と関数 $g$ について、$\frac1n\sum_{i=1}^ng(U_i)$ の平均は $\int_0^1g(x)\,dx$ で、標準偏差は $\frac{\sigma_g}{\sqrt n}$($\sigma_g$ は $g(U)$ の標準偏差)である。thm-rns-main は、同じことを 2 変数の関数 $g(x,y)$(四分円の内側なら $4$、外側なら $0$)と正方形 $[0,1)^2$ で行った場合にあたる。このように乱数で積分を見積もる方法を Monte Carlo 法 という(Monte Carlo法)。
$d$ 個の変数の関数を $d$ 次元の立方体 $[0,1)^d$ で積分するとき、区分求積のように各辺を $k$ 等分した格子の点を使うと、点は $k^d$ 個必要になる。$d=10$ で各辺を $10$ 等分するだけで $10^{10}$ 個である。一方、乱数による見積もりの標準偏差 $\frac{\sigma_g}{\sqrt n}$ の式には $d$ が現れない。$d$ 個の乱数で 1 つの点を作り、それを $n$ 個使えば、次元によらず誤差は $\frac1{\sqrt n}$ の速さで減る。このため、物理学や金融の計算など、変数の多い積分や確率の計算で乱数が使われる。
ただし ex-rns-riemann のように 1 次元なら、区分求積(や、さらに誤差の小さい台形公式・中点公式)のほうがずっと速い。
平行線を等間隔に引いた床に針を落とし、線と交わる割合から $\pi$ を見積もる方法もある(Buffonの針)。GS06 §2.1 は、乱数で正方形に点を打って面積を見積もる方法(本記事の ex-rns-pi-intro と同じ考え方)に続けて、Buffon の針を扱っている。同書は、確率を見積もるシミュレーションの誤差が、およそ 95% の割合で $\frac1{\sqrt n}$ 以下になることも注意している。これは cor-rns-trials と同じ考え方で、1 回ごとの結果が $1$ か $0$ のときの標準偏差 $\sqrt{p(1-p)}$ が $\frac12$ 以下であることから出る。
正方形に $10000$ 点を打ったところ、四分円の内側に $7890$ 点が入った。推定値 $\hat\pi$ を求め、$\pi$ との差が thm-rns-main の標準偏差の何倍かを求めよ。この結果は、乱数に問題があることを疑わせるか。
$\hat\pi=4\times\frac{7890}{10000}=3.156$ である。$\pi$ との差は $3.156-3.14159=0.01441$ で、標準偏差 $\frac{1.6422}{\sqrt{10000}}=0.016422$ の $0.88$ 倍である。cor-rns-trials の近似では、差が標準偏差の $1.96$ 倍以内に収まる確率が約 $0.95$ なので、$0.88$ 倍はふつうに起こる大きさで、これだけで乱数を疑う理由にはならない。
$m=9$、$a=4$、$c=1$、$x_0=0$ の線形合同法の列を、同じ値に戻るまで書き出し、周期を求めよ。
$x_1=1$、$x_2=5$、$x_3=21-18=3$、$x_4=13-9=4$、$x_5=17-9=8$、$x_6=33-27=6$、$x_7=25-18=7$、$x_8=29-27=2$、$x_9=9-9=0$ で、$x_9=x_0$ となる。列は $0,1,5,3,4,8,6,7,2,0,\dots$ で、$0$ から $8$ までの $9$ 個の値がちょうど 1 回ずつ現れ、周期は $9$($m$ に等しい)である。
Mathpediaは寄付と、参考文献の書籍リンク(Amazonアソシエイト)の紹介料で運営されています。 支援について / 寄付する