和が１を超えるまで加算すべき回数の期待値

# 本日は

ある種の期待値の計算をやってみた話です。

問題はこちら：

(0,1) 区間上で値を取る一様乱数を逐次発生させ(発生の確率は互いに独立)、それらを加算していきます。加算は累積した値が１を超えるまで行い続けるものとします。さて、その加算回数の平均値はいくらでしょう？

すなわち、

```math
x_1, x_2, x_3, x_4,\dots 
```

を一様乱数の列として

```math
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でもうごくとおもいます。

```julia:calc.jl
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()
```

```shell-session
$ 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)$ 上の一様分布に従う確率変数とします。

```math
X_1+X_2+\dots+X_{n-1} < 1,\quad X_1+X_2+\dots+X_{n} \geq 1
```

となる $n$ を値に取る確率変数を $N$ とします。 
冒頭の問題は $N$ の期待値 $E[N]$ を計算することです:

```math
E[N] = \sum_{n=1}^{\infty}nP(N=n).
```

ここで、$P(N=n)$ は $N=n$ となる確率です。この計算を次のように進めていきます。

```math
\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の箇所について：

和の順序を入れ替えて

```math
P(N=1)+P(N=2)+P(N=3)+P(N=4)+\dots = P(N\geq 1)=1,
```
```math
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 の箇所について:

これは　

```math
P(X_1+\dots+X_n<1)=\frac{1}{n!}\quad\mathrm{for}\  n \geq 1
```

が成り立つからです。これはコーシーの反復積分の公式と同様に計算を進めると示すことができます。

うーん。$e$ になるなんておもしろいですね。ちなみにこの事実は最近知りました。これって有名な定理ですか？

# 最後に

並列して計算させるコードをPython＋NumbaとJuliaのコードを紹介して終わろうと思います。

```julia:napier.jl
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()

```

```python:napier.py
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による探索プログラミング―基礎から遺伝的アルゴリズムまで](https://www.amazon.co.jp/C%E3%81%AB%E3%82%88%E3%82%8B%E6%8E%A2%E7%B4%A2%E3%83%97%E3%83%AD%E3%82%B0%E3%83%A9%E3%83%9F%E3%83%B3%E3%82%B0%E2%80%95%E5%9F%BA%E7%A4%8E%E3%81%8B%E3%82%89%E9%81%BA%E4%BC%9D%E7%9A%84%E3%82%A2%E3%83%AB%E3%82%B4%E3%83%AA%E3%82%BA%E3%83%A0%E3%81%BE%E3%81%A7-%E4%BC%8A%E5%BA%AD-%E6%96%89%E5%BF%97/dp/4274067300) で知りました。

## その他いくつかの情報

- [反復積分に関するコーシーの公式](https://ja.wikipedia.org/wiki/%E5%8F%8D%E5%BE%A9%E7%A9%8D%E5%88%86%E3%81%AB%E9%96%A2%E3%81%99%E3%82%8B%E3%82%B3%E3%83%BC%E3%82%B7%E3%83%BC%E3%81%AE%E5%85%AC%E5%BC%8F)

- [What is the expected number of random numbers (generated uniformly) such that their sum of numbers exceeds one?](https://mathoverflow.net/questions/141942/what-is-the-expected-number-of-random-numbers-generated-uniformly-such-that-th)

- [Uniform Sum Distribution](http://mathworld.wolfram.com/UniformSumDistribution.html)
- [Irwin–Hall distribution](https://en.wikipedia.org/wiki/Irwin%E2%80%93Hall_distribution)