前回のあらすじ
赤玉４つ、黒玉３つ、白玉５つの合計１２個の玉が入った袋から、一つずつ玉を取り出して一列に並べる、という試行についての問題。入試問題では「赤玉が隣合わない確率$p_0$」を計算することになっているが、まずは「玉を取り出して並べる」シミュレーターをpythonで書いてみた（aka-dama.py)。このプログラムの動作確認として、赤玉が最初の１、２番目に引き当てられる確率（$p_\text{XX}=1/11$)を数値実験で確認した。このとき、大数の法則を考慮するとサンプリングを6000回以上やると次第に理論値に近づくことも確認した（収束は遅かったが）。
今回は、同じシミュレータを使って、いよいよ入試問題で問われている「赤玉が隣合わない確率$p_0$」の数値実験に挑む。
予想
まずは$p_0$が５０％より大きいのか小さいのか予想を立ててみよう。前回の記事でもちょっと触れているが、$p_\text{XX}=1/11$の確認に気を取られ過ぎて、考察が途中で終わってしまった（反省）。
まず簡単に計算できる例として最初の1,2回目が赤玉である確率を計算し、すでに$p_{12}=1/11$という結果を「理論的」にも数値実験的にも得ている。同じ考え方を$p_{n,n+1}$について適応してみよう（つまり$n$回目と$n+1$回目が連続して赤になる確率）。
事象A (連続する2つの赤玉が最初の赤玉の場合）
$ p_{12} $と違うのは1回目から$n-1$回目の間に「散発的に」赤玉が出てしまう事象を考慮しなくてはならないという点である。もし$n-1$回目までに赤玉がまったく出なければ、赤玉は４つ丸々残っている。一方、$n$回目の試行の際、袋には非赤玉は$8-n$個残っている。全体としては$(8-n)+4=12-n$個残っている。したがって、$n\ge 2$に対して \begin{equation}p_{n,n+1}=\frac{4}{13-n}\frac{3}{12-n}q_{n-1}\end{equation} となる。$q_{n-1}$は$n-1$回目までに赤玉がまったく出ない確率であり、 \begin{equation}q_{n-1} = \frac{8}{12}\frac{7}{11}\cdots\frac{10-n}{14-n}\end{equation} という形にまとめることができる。上の2つの式をまとめ、階乗の部分をうまく処理すると \begin{equation} p_{n,n+1} = \frac{(11-n)(10-n)}{11\cdot 10\cdot 9}\end{equation}という形になる。 $n\ge 2$という前提で考察してきたが、$n=1$を上の式に代入すると上手い具合に$p_{12}=1/11$が再現できていることがわかる。 したがって、上の式は$n\ge 1$で成立する。また、非赤玉の総数は８個なので、$n=9$が上限であることもわかる。 したがって、$n=1,2,\cdots,9$についての和をとれば、「連続赤玉の前に赤玉が一つも出ない」という事象の確率$p_A$が得られる。 \begin{equation}p_A = \sum_{n=1}^9 p_{n,n+1} = \frac{1}{11\cdot 10\cdot 9}\sum_{n}(11-n)(10-n)\end{equation}
かっこの中を展開すると \begin{equation} 11\cdot 10 - 21n + n^2\end{equation}となるから、初項に関しては \begin{equation}\frac{1}{11\cdot 10\cdot 9}\sum_{n=1}^9 11\cdot 10 = \frac{11\cdot 10}{11\cdot 10\cdot 9}\sum_{n=1}^9 1 = \frac{11\cdot 10}{11\cdot 10\cdot 9}\cdot 9 = 1\end{equation}
第二項は和の公式を用いて \begin{equation}-\frac{21}{11\cdot 10\cdot 9}\sum_{n=1}^9 n = -\frac{3\cdot 7}{11\cdot 10 \cdot 9}\frac{9\cdot (9+1)}{2} = -\frac{21}{22} \end{equation} となる。第３項も同様に和の公式を用いて \begin{equation}\sum_{n=1}^9 n^2 = \frac{9(9+1)(2\cdot 9 + 1)}{6} = \frac{19}{66}\end{equation} を得る。まとめると、 \begin{equation}p_A = 1-\frac{21}{22} + \frac{19}{66} = \frac{66-63+19}{66} = \frac{22}{66} = \frac{1}{3}\end{equation}となり、意外にも綺麗な形に落ち着いた。
もちろん、問題となっている事象の補事象としてはこれだけではまだまだ不足である。とはいえ「不足の程度」というものがあるので、なんとなくこの結果から全体の見通しができないかと期待するわけだが、1/3というのはいささか微妙な値である。果たして赤玉が連続しない方が起きやすいのか、それとも起きにくいのか？残念ながら、なんとも言えないのが現状である。
事象B（連続する赤玉の前に１度だけ赤玉が出ている場合）
しかたないので、もう一つのパターンを考えてみよう。二つの赤玉連続の前に単発で赤玉が出てしまう場合である。我々の記号を使えば、XO@@XXOOXO@Oといったような場合である。XXの連続の前にXが一つだけ出てくる事象である。上の場合の考え方に則れば、一般的な形を求めることはそれほど難しくない。
n回目とn+1回目に赤玉(X)が出るとする。$n-1$回目までに赤玉は１つ出てしまっているので、赤玉の残りは$n$回目の試行の際は３つである。つまり、このタイミングで赤玉が連続する確率は\begin{equation}\frac{3}{13-n}\cdot\frac{2}{12-n}\end{equation}である。$n-1$回目までに赤玉が1回だけでる確率として、まずは赤玉が$k$回目に出る確率$q_{n-1}^{(k)}$を考えてみる。ただし$1\le k \le n-2$である。$k=n-1$の場合を除外する理由は、このケースは3回赤玉が連続するケースとなってしまい、それは$p_{n-1,n}$に寄与するべき事象になるからである（数えすぎになってしまう）。 \begin{equation} q_{n-1}^{(k)} = \frac{8}{12}\cdots \frac{4}{12-k}\cdots \frac{11-n}{14-n}\end{equation}となる。赤玉が出るタイミングが$k=1$から$k=n-2$まであるものの、$q^{(k)}$の中身（つまり、分子と分母に現れる整数の組み合わせ）は変わらないから、 \begin{equation}q'_{n-1} = \sum_{k=1}^{n-2}q_{n-1}^{(k)} = (n-2)\frac{8\cdot 7 \cdots (11-n)\cdot 4}{12\cdot 11 \cdots (14-n)}\end{equation}となる。 したがって、全体の確率は（分子分母に階乗$(11-n)!$をかけて整理すると） \begin{equation}p_{n,n+1}=\frac{3}{13-n}\cdot\frac{2}{12-n} q'_{n-1} = \frac{4!8!}{12!}(11-n)(n-2) = \frac{13n-n^2-22}{5\cdot 9 \cdot 11}\end{equation}となる。具体例を使うなどして簡単に考察すれば、$n=3,4,\cdots, 9$の場合に有効であることが確認できる。したがって、この場合の全確率は \begin{equation}P_B = \sum_{n=3}^{9}p_{n,n+1}=\frac{14\cdot 8}{55\cdot 9}\end{equation}となる。
事象AとBの合計
先に計算しておいた$p_A$を合わせると \begin{equation} p_A + p_B = \frac{1}{3} + \frac{14\cdot 8}{55\cdot 9} = \frac{277}{495} \simeq 0.56\cdots >\frac{1}{2} \end{equation}となり、ようやくここにきて「赤玉が連続した方が確率は大きい」という「結果」が得られた。
赤玉が連続する「数え残し」の事象は、赤玉が2つ連続する前に「散発的に2回赤玉が出る」事象である。たとえば、 XOX@@@XXOOOO といったような例である。これを数え尽くせば「補事象」が出尽くすがなかなか面倒くさそうな雰囲気である。というのは、前半部分の赤玉は連続してはいけないからである。このあたりの場合分けが多いと「数え間違い」が発生し、間違いということになるだろう。
もちろん、ここまでの計算があっているとは限らない。ということで、いつものように数値計算して確かめてみよう。
event Aの数値計算
pythonで作ったSequence generatorを使って数値実験を行う。前回は11回連続で赤白黒のタマのsequenceを生成し、それを１セットとした。各セット毎に確率を計算し、その確率の平均値を複数回のセットについてとって最終結果とした。500セットの計算で30秒かかったから、5500個のsequence生成に30秒もかかったことになる。これはなかなかイライラする。
なんとか計算時間を短縮しようと思い、今回はSequence generatorで一気にsequenceを多数生成するやり方に変更してみた。11万個生成で10秒、110万個で108秒だったので、列の生成だけなら、この方法の方が圧倒的に早い。問題は判定にどれだけ時間がかかるかである。まずはA事象の$p_{2,3}=4/55\simeq 0.07272\cdots$を確かめてみた。判定にはawkを使った。
time aka-dama.py |awk '{if( substr($2,1,1) != "X" && substr($2, 2, 2) == "XX" ) print $1,$2;}'|wc -l
11万個の場合は0.2秒、110万個の場合は3秒だった。やはり、コンピュータは同じことを繰り返すのが得意なので、一度やり出したら手続きをあまり細切れにしない方が速いようだ。
結果の方は、11万個の場合は$8102/110000 = 0.07365...$, 110万個の場合は$80181/1100000 = 0.07289....$であった。おそらく解析的な考察とほぼ同じ値、しかも大数の法則が効いている感じに見えるから、まあこれで大丈夫だろう。
$p_{3,4}= 28/495\simeq 0.0565656...$の方も実験してみた。awkのスクリプトを少しだけ複雑にするだけで済む。
time aka-dama.py |awk '{if( substr($2,1,1) != "X" && substr($2,2,1) != "X" && substr($2, 2, 2) == "XX" ) print $1,$2;}'|wc -l
結果は、11万個の場合が$6307/110000=0.05733...$、110万個の場合が$62223/1100000=0.0565663....$となり、かなり良い精度で考察で得た確率と一致した。
事象Aの確率の総和を求めるために、$n=1, \cdots, 9$について$p_{n,n+1}$の和を数値的に計算してみた。110万個の場合、
1 99663 2 80181 3 62223 4 46364 5 33466 6 22072 7 13345 8 6650 9 2211
という結果が得られた。総和をとって110万で割ると、 \begin{equation}366175/1100000=0.33288\cdots\end{equation}となり、考察で得た$p_A=1/3=0.3333...$に非常に近い結果が得られた！数え落としや勘違いはなさそうである。
