乱数とシミュレーション

同義語:乱数によるシミュレーションsimulation with random numbers

概要

乱数とシミュレーション(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$ 以下である。周期が短い、並びに規則が残る、点が一様に散らばらないときは、結果が偏る。

$$\newcommand{C}[0]{\mathbb{C}} \newcommand{div}[0]{\mathbin{÷}} \newcommand{N}[0]{\mathbb{N}} \newcommand{Q}[0]{\mathbb{Q}} \newcommand{R}[0]{\mathbb{R}} \newcommand{Z}[0]{\mathbb{Z}} $$

前提知識: 大数の法則とさいころの平均, 標本平均の分散, 中心極限定理(高校数学), 確率密度関数と連続型確率変数

高校での出発点:乱数で確率を実験する

数学 B では、さいころ・乱数表・計算機の 乱数 を使って試行をくり返し、起こった割合から確率を見積もる方法を学ぶ。数えるのが難しい確率でも、試行を何度もくり返せば、相対度数は確率に近づく(大数の法則とさいころの平均)。このような実験を シミュレーション という。まず、計算機で実際に試してみる。

さいころで同じ目が 3 回続く確率

さいころを 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 点。四分円の内側の点を青、外側の点を橙で示す。 一辺 1 の正方形に打った 100 点。四分円の内側の点を青、外側の点を橙で示す。
同じ乱数で打った 1000 点。内側の点の割合が四分円の面積 π/4 に近いことを見る図。 同じ乱数で打った 1000 点。内側の点の割合が四分円の面積 π/4 に近いことを見る図。

1000 点も打ったのに、円周率が 1 桁しか合わない。ここで次の問いが出てくる。

  1. 計算機の「乱数」とは何で、それからさいころやくじをどう作るのか。→ def-rns-random、prop-rns-transform
  2. $4\times\frac mn$ はどのくらいの誤差で $\pi$ に近いのか。$n$ を増やすと誤差はどう減るのか。→ thm-rns-main
  3. 誤差を $0.01$ にするには、何点打てばよいのか。→ cor-rns-trials
  4. 区分求積で面積を計算する方法と比べて、どちらがよいのか。→ ex-rns-riemann
  5. 計算機はどうやって乱数を作るのか。どんな乱数ではいけないのか。→ def-rns-lcg、ex-rns-counter-short、ex-rns-counter-parity
    高校の計算この記事での見方ボックス
    乱数さいころ・乱数表独立で一様な乱数の列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

乱数とその使い方

乱数
  1. $0,1,\dots,9$ の数字の列 $d_1,d_2,\dots$ で、各 $d_k$ がどの数字も確率 $\frac1{10}$ でとり、互いに独立なものを 一様乱数の列(数字の乱数)という。
  2. $0$ 以上 $1$ 未満の実数の列 $U_1,U_2,\dots$ で、各 $U_k$ が区間 $[0,1)$ 上の一様分布に従い($0\le a\le b\le1$ なら $P(a\le U_k< b)=b-a$)、互いに独立なものを、区間 $[0,1)$ の一様乱数の列 という。

乱数さいころ(正二十面体に $0$ から $9$ を 2 回ずつ書いたもの)を振ると (i) が得られる。乱数表は、そのような数字を並べた表である。計算機の乱数は (ii) を目指して作られるが、実際には決まった計算で作る数の列であり、擬似乱数 という(節「擬似乱数の作り方:線形合同法」)。GS06 §1.1・§2.1 も、計算機の乱数を「区間に入る確率が区間の長さに等しい」ように選ばれる数として導入し、実際には決まった手順で作られる数であることを注意している。
一様乱数から、いろいろな確率の試行を作れる。

一様乱数から試行を作る

$U$ を区間 $[0,1)$ の一様乱数とする。

  1. $0\le p\le1$ について、事象「$U< p$」の確率は $p$ である。
  2. 正の整数 $k$ について、$\lfloor kU\rfloor+1$ は $1,2,\dots,k$ のどれも確率 $\frac1k$ でとる。ここで $\lfloor x\rfloor$ は $x$ 以下の最大の整数(ガウス記号 $[x]$ と同じ)である。
区間の長さで数える
  1. def-rns-random の (ii) で $a=0$、$b=p$ とすると $P(0\le U< p)=p$ で、$U\ge0$ はいつも成り立つので $P(U< p)=p$ である。
  2. $j\in\{1,\dots,k\}$ とする。$\lfloor kU\rfloor+1=j$ は $\lfloor kU\rfloor=j-1$、つまり $j-1\le kU< j$ と同じで、両辺を $k$ で割ると $\frac{j-1}k\le U<\frac jk$ と同じである。(ii) からこの確率は区間の長さ $\frac jk-\frac{j-1}k=\frac1k$ である。$0\le U<1$ なので $0\le kU< k$ で、$\lfloor kU\rfloor+1$ は $1$ から $k$ の外の値をとらない。$\square$
乱数でさいころ・硬貨・くじを作る

計算機が $U=0.7312$ を返したとする。

  1. さいころ:$\lfloor6\times0.7312\rfloor+1=\lfloor4.3872\rfloor+1=5$ で、目は $5$。
  2. 硬貨:$U<0.5$ なら表、そうでなければ裏とすると、$0.7312\ge0.5$ なので裏。
  3. 中心極限定理(高校数学) のくじ($0$ 円 $0.7$、$100$ 円 $0.2$、$1000$ 円 $0.1$):$U<0.7$ なら $0$ 円、$0.7\le U<0.9$ なら $100$ 円、$0.9\le U$ なら $1000$ 円とする。それぞれの区間の長さは $0.7$、$0.2$、$0.1$ なので、def-rns-random の (ii) から確率はそのとおりになる。$U=0.7312$ なら $100$ 円である。
    数字の乱数を使うときは、たとえば 4 けたずつ区切って $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$ 倍にしなければならない。

1 か 0 の量の平均と分散から

方針:$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$

$n$ ごとの標準偏差と実際の値

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 の直線に沿って減ることを見る。 乱数による π の推定の誤差(点)と理論値(実線)、右端の区分求積の誤差(破線)を、横軸・縦軸とも対数目盛で描いた図。乱数の誤差は傾き −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}$ である。

中心極限定理と Chebyshev の不等式

$\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$

誤差 0.01 と 0.001
  1. $\varepsilon=0.01$:近似では $n\ge\bigl(\frac{1.96\times1.6422}{0.01}\bigr)^2=321.9^2=103600$ 程度、つまり約 $10$ 万点である。Chebyshev の不等式による保証では $\frac{\pi(4-\pi)}{0.05\times0.0001}=\frac{2.69676\ldots}{0.000005}=539353.2\ldots$ から $539354$ 点、約 $54$ 万点が必要になる。
  2. $\varepsilon=0.001$:近似では約 $1036$ 万点である。誤差を $\frac1{10}$ にするために、点の数は $100$ 倍になる。
    小数第 3 位まで正しい $\pi$ を得るのに 1000 万点以上が必要で、乱数による推定は $\pi$ の計算方法としては効率がよくない。

区分求積との比較

$\frac\pi4$ は四分円の面積なので、$\frac\pi4=\int_0^1\sqrt{1-x^2}\,dx$ である。乱数を使わずに、この積分を区分求積で計算することもできる(区分求積と積分の定義)。

区分求積で $\pi$ をはさむ

$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$ の線形合同法

$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$ と互いに素のとき)で決まる。

$c=0$ の場合と冪の余り
  1. $m=7$、$a=3$、$c=0$、$x_0=1$:$1,3,2,6,4,5,1,\dots$ で周期 $6$ である。これは $3^k$ を $7$ で割った余りの列で、$3^6$ で初めて余りが $1$ に戻ることに対応する(冪の余りの周期と元の位数)。
  2. $m=16$、$a=3$、$c=0$、$x_0=1$:$1,3,9,11,1,\dots$ で周期は $4$ しかない。$16$ 個の値のうち $4$ 個しか現れない。

周期が $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 本の直線の上に並ぶことを見る図。 線形合同法 x_{k+1} = (5x_k + 1) mod 256 の連続する 2 つの値の組。256 点が傾き 5 の 5 本の直線の上に並ぶことを見る図。
周期のとても長い擬似乱数の連続する 2 つの値の組 256 点。正方形の中にばらばらに散らばることを見る図。 周期のとても長い擬似乱数の連続する 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$)が成り立たなくなる。

大学数学で見る

乱数による積分(Monte Carlo 法)

$[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 次元なら、区分求積(や、さらに誤差の小さい台形公式・中点公式)のほうがずっと速い。

Buffon の針

平行線を等間隔に引いた床に針を落とし、線と交わる割合から $\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$ に等しい)である。

さらに先へ

  • 実際の計算機では、周期が非常に長く、並び方の規則が見つかりにくいように設計された擬似乱数が使われている。擬似乱数の「よさ」を統計的に調べる検定も研究されている(擬似乱数)。
  • 仮説検定(高校数学) で、帰無仮説のもとでの分布が式で求めにくいときは、乱数でその分布を作り、棄却域に入る割合を見積もることができる。
  • 乱数の代わりに、規則的だが偏りの少ない点の列を使って積分の誤差を $\frac1{\sqrt n}$ より速く減らす方法(準モンテカルロ法)もある。

関連項目

参考文献

Mathpediaは寄付と、参考文献の書籍リンク(Amazonアソシエイト)の紹介料で運営されています。 支援について / 寄付する