何回コインを投げれば連続10回表が出るんだろう(1)

http://blogs.wankuma.com/episteme/archive/2007/12/13/112757.aspx


この分布を擬似乱数で出すのは難しい、擬似乱数の質を調べるのによい題材、ということだと思う。


しかし、それにはまず正しい分布を知ることが必要。
というか、私にはそちらを出すほうが面白い。


m回コインを投げてはじめてn回連続表になる確率をp(m,n)、
m回コインを投げて現在n回連続表である確率をq(m,n)とする。


m > nのとき、裏のあと表がn回続くのがq(m,n)だから、2-n-1
他の場合とまとめると、

 q(m,n) = \begin{cases} 0 & (m < n) \\ 2^{-n} & (m = n) \\ 2^{-n-1} & (m > n) \end{cases}

現在n回連続表であるとき、はじめてn回連続になった場合と、それ以前にn回連続になったことがある場合がある。それは、m > 2nのとき。m - n - 1回までは任意で、その次が裏、そのあとn回表となる。だから、nからm - n - 1回コインを投げてはじめてn回になった可能性がある。これを式にすると、

 q(m,n) = \begin{cases} (p(n,n) + \cdots + p(m-n-1,n))2^{-n-1} + p(m,n) & (m > 2n) \\ p(m,n) & (n \le m \le 2n) \\ 0 & (m < n) \end{cases}

この漸化式から各pと期待値を求めることができる。


var n = 10;
var N = 100000;
var p = [ ];

p[n] = 1;
for(var i = 0; i < n; i++) p[n] /= 2;

var E = n * p[n]; // 期待値
for(var m = n + 1; m <= n * 2; m++) {
p[m] = p[n] / 2;
E += m * p[m];
}

var sum = p[n];
for(var m = n * 2 + 1; m <= N; m++) {
p[m] = p[n+1] * (1 - sum);
E += m * p[m];
sum += p[m-n];
}

pのグラフは次のようになる。

そして、期待値は2046.00000016722となった。
10万で打ち切っているが、そのあとはどうなるのだろう。

 r(m) = p(n, n) + \cdots + p(m, n)

とすると、漸化式を微分で近似できて、

 \frac{dr}{dm} = 2^{-11}(1 - r)

これを解いて積分すると、残りの和は、
p(N,n) = 5.50E-17より、

 (2^{11}N + 2^{22})p(N,n) = 1.15E-8

となり、無視できる。
よって、期待値は約2046になる。
ただ、しょせん浮動小数点演算なので誤差がある。