本日は
ある種の期待値の計算をやってみた話です。
問題はこちら：
(0,1) 区間上で値を取る一様乱数を逐次発生させ(発生の確率は互いに独立)、それらを加算していきます。加算は累積した値が１を超えるまで行い続けるものとします。さて、その加算回数の平均値はいくらでしょう？
すなわち、
x_1, x_2, x_3, x_4,\dots
を一様乱数の列として
x_1+x_2+\dots+x_{n-1} < 1,\quad x_1+x_2+\dots+x_{n} \geq 1
となる $n$ の取りうる期待値はいくらでしょうか？
プログラムを書いて計算してみる
Julia で書いてみます。（Julia Version 1.0.0） で動作確認していますが、ふるいJuliaでもうごくとおもいます。
function count_upto_one()
counter = 0
accumulated = 0.0
while true
accumulated += rand()
counter += 1
if accumulated >= 1.0
break
end
end
return counter
end
function main()
n_trial =1e8
total_count=0
for trial in 1:n_trial
total_count+=count_upto_one()
end
println(total_count/n_trial)
end
main()
$ julia calc.jl
2.718199
（以下何回か同じコマンドを実行）
2.71825619
2.71820856
2.71839289
2.71836157
ふーむ。どうやら２から３の間の実数のようです。
ちなみに、JuliaのREPLを起動して
\euler と書いてTABを打つと自然対数の底である ℯ が得られます。
この値は
julia> ℯ
ℯ = 2.7182818284590...
でした。値が近いのは偶然でしょうか？
計算してみましょう
$X_n$ を $n$ 回目に発生させる区間 $I=(0,1)$ 上の一様分布に従う確率変数とします。
X_1+X_2+\dots+X_{n-1} < 1,\quad X_1+X_2+\dots+X_{n} \geq 1
となる $n$ を値に取る確率変数を $N$ とします。
冒頭の問題は $N$ の期待値 $E[N]$ を計算することです:
E[N] = \sum_{n=1}^{\infty}nP(N=n).
ここで、$P(N=n)$ は $N=n$ となる確率です。この計算を次のように進めていきます。
\begin{align}
E[N] &\overset{\mathrm{def}}{=} \sum_{n=1}^{\infty}nP(N=n) \\
&=\quad P(N=1)\\
&\quad +P(N=2)+P(N=2)\\
&\quad +P(N=3)+P(N=3)+P(N=3)\\
&\quad +P(N=4)+P(N=4)+P(N=4)+P(N=4) \\
&\quad +\cdots\\
&\overset{A}{=}1 + \sum_{n=2}^{\infty} P(N\geq n) \\
&\overset{B}{=}1 + \sum_{n=1}^{\infty} P(X_1+\dots+X_n < 1) \\
&\overset{C}{=}\sum_{k=0}^{\infty}\frac{1}{k!} = e.
\end{align}
上記計算過程において、
- Aの箇所について：
和の順序を入れ替えて
P(N=1)+P(N=2)+P(N=3)+P(N=4)+\dots = P(N\geq 1)=1,
P(N=2)+P(N=3)+P(N=4)+\dots = P(N\geq 2)
のように式を変形しています。
-
Bの箇所について:
$P(X_1+\dots+X_n < 1)$ は確率変数の和が１未満を取る確率、すなわち、 $X_1+\dots+X_n<1$ となる確率を表しています。ここでの式変形の気持ちとしては、$n$ 乱数を取得してもまだ１を超えないという状態は目標を達成する試行回数が $n+1$ 以上であることとおなじだよねということです。 -
C の箇所について:
これは
P(X_1+\dots+X_n<1)=\frac{1}{n!}\quad\mathrm{for}\ n \geq 1
が成り立つからです。これはコーシーの反復積分の公式と同様に計算を進めると示すことができます。
うーん。$e$ になるなんておもしろいですね。ちなみにこの事実は最近知りました。これって有名な定理ですか？
最後に
並列して計算させるコードをPython＋NumbaとJuliaのコードを紹介して終わろうと思います。
using Distributed
function count_upto_one()
counter = 0
accumulated = 0.0
while true
accumulated += rand()
counter += 1
if accumulated >= 1.0
break
end
end
return counter
end
function calc_napier_parallel()
n_trial =1e9
total_count = @distributed (+) for i in 1:n_trial
counter = 0
accumulated = 0.0
while true
accumulated += rand()
counter += 1
if accumulated >= 1.0
break
end
end
counter
end
println(total_count/n_trial)
end
addprocs()
@time calc_napier_parallel()
import numpy as np
from numba import jit, prange
@jit(nopython=True)
def count_upto_one():
counter = 0
accumulated = 0.0
while True:
accumulated += np.random.random()
counter += 1
if accumulated >= 1.0:
break
return counter
@jit(nopython=True, parallel=True)
def main():
n_trial = int(1e9)
total_count = 0
for _ in prange(n_trial):
total_count += count_upto_one()
print(total_count / n_trial)
if __name__ == '__main__':
main()
参考文献
元ネタ
- Cによる探索プログラミング―基礎から遺伝的アルゴリズムまで で知りました。