コインを20000回投げて11000回以上表が出る確率

以前こちらで考察したコイン投げ試行。
もう3年以上前という現実に、切なくなります。

このような（自分にとって）直感的でない確率の問題は好きなので、もう一度考察してみます。

まずは前回と同じ、2000回投げて1100回以上表が出る確率です。
初めてこの問題を考えたとき、直感的にはそこそこの確率で起こりそうに感じましたが、
実際には0.00038%という非常に小さい値でした。
もう一度順を追って計算してみます。

まず、
コインを投げて表が出る確率
$$ p=\frac{1}{2} $$

表が出る回数の平均
$$ \mu = 2000 \times \frac{1}{2} = 1000 $$

標準偏差
$$ \sigma = \sqrt{2000\times\frac{1}{2}\times(1-\frac{1}{2})} = \sqrt{500} = 10\sqrt{5}$$

であり、コイン投げの二項分布は試行回数が増えれば正規分布に近似することを考えると、確率変数 $X$ は $ N(1000, 500) $ に従います。

次に、$Z = \frac{X-1000}{10\sqrt{5}}$ による標準化を行います。

表が出る回数が1100回以上2000回以下となる確率は、

\begin{align}
P(1100 \le X \le 2000) &= P\left(\frac{1100-1000}{10\sqrt{5}} \le Z \le \frac{2000-1000}{10\sqrt{5}}\right) \\
&= P(2\sqrt{5} \le Z \le 20\sqrt{5}) \\
&= P(0 \le Z \le 20\sqrt{5}) – P(0 \le Z \le 2\sqrt{5}) \\
&= \displaystyle \int_{0}^{20\sqrt{5}} \frac{1}{\sqrt{2\pi}} \exp{\left(-\frac{z^2}{2}\right)} dz – \int_{0}^{2\sqrt{5}} \frac{1}{\sqrt{2\pi}} \exp{\left(-\frac{z^2}{2}\right)} dz \end{align}

となります。
前回は電卓を叩きましたが、今回は次のような単純なプログラムを作成し、計算させてみました。

	public static void main(String[] args) {
		double step = 1e-6;
		double p = integralND(0, 20*Math.sqrt(5), step) - integralND(0, 2*Math.sqrt(5), step);
		
		System.out.println("p = " + p);
	}
	
	private static double integralND(double a, double b, double step) {
		double res = 0.0;
		double fix = 1.0 / Math.sqrt(2*Math.PI);
		
		for(double x = a; x <= b; x+=step) {
			res += fix * Math.exp(-x*x/2.0) * step;
		}
		return res;
	}

また、より厳密な解を力技で求める方法も試してみます。
今回求めたい確率は、次の式で表すことができます。

$$ p = \displaystyle \sum_{i=1100}^{2000}\frac{ _{2000}C_{i}}{2^{2000}} = \sum_{i=0}^{900}\frac{ _{2000}C_{i}}{2^{2000}}$$

かなり大きな数になってしまいますが、任意の精度を扱えるBigInteger, BigDecimalクラスを用いることで計算可能でした。

	public static void main(String[] args) {
		int N = 2000;
		int p = 900;
		BigInteger denom = bigPow(2, N);
		BigInteger numer = BigInteger.valueOf(0);
		BigDecimal denomD, numerD;

		for (int i = 0; i <= p; i++) {
			numer = numer.add(bigCmb(N, i));
		}

		denomD = new BigDecimal(denom);
		numerD = new BigDecimal(numer);

		System.out.println("p = " + numerD.divide(denomD).setScale(50, BigDecimal.ROUND_DOWN));
	}

	private static BigInteger bigPow(int a, int b) {
		BigInteger res = BigInteger.ONE;

		for (int i = 0; i < b; i++) {
			res = res.multiply(BigInteger.valueOf(a));
		}
		return res;
	}

	private static BigInteger bigCmb(int N, int m) {
		BigInteger res = BigInteger.ONE;
		int rep = (int) Math.min(N - m, m);

		for (int i = 0; i < rep; i++) {
			res = res.multiply(BigInteger.valueOf(N - i));
		}
		for (int i = 0; i < rep; i++) {
			res = res.divide(BigInteger.valueOf(rep - i));
		}
		return res;
	}

前者の結果は

　p = 0.00000387211 = 0.00038711%

後者の結果は

　p = 0.0000042285 = 0.00042285%

スケールは同じですが、数値は結構ずれる結果になってしまいました。
この試行回数であればかなり精度良く正規分布に近似されているはずなので、
もう少し近い値になってくれることを期待していました・・・。
どちらかの方法（あるいは両方）に間違いがあることと思いますが、どこに間違いがあるのかはわからない現状です。

原因解明は一旦置いておいて、他のケースの結果を載せていきます。

前者の結果は

　　p = 0.32735 = 32.73%

後者の結果は

　　p = 0.41190 = 41.19%

起こりそうです。
大きく値がずれているのは、試行回数が少なくて正規分布に近似できないのが原因と思われます。
後者の結果が正しいはずです。

前者の結果は

　　p = 0.078649 = 7.8649%

後者の結果は

　　p = 0.089482 = 8.9482%

大分低くなりました。しかし値がずれます。
正規分布を使った方法に間違いがあるような気がしてきました。

前者の結果は

　　p = 0.00000387211 = 0.00038711%

後者の結果は

　　p = 0.0000042285 = 0.00042285%

二度目の結果掲載です。
やはり、私個人的にはこの確率が最も直感的でないように感じます。
センスのある人なら騙されないことと思いますが、
なんとなく 1% くらいはあるような気がしてなりません。

これは、前者のやり方では計算不可能でした。
後者の結果は

　　p = 9.7704E-46 = 0.0000000000000000000000000000000000000000000001%

とんでもなく小さな結果となります。
まず起こりえない確率と言えるでしょう。

コインを N 回投げて 0.55N 回以上表が出る確率を求めるこの問題の結果は、
私にとっての直感的な結果と大きく異なり、なかなか衝撃的でした。

しかし、私に確率的センスが無いだけで
「こんなの直感的にわかるよ」と思われる方が多数派なのかもしれません。
少し気になります。