前回のあらすじ
東京大学2023年の数学入試問題で出題された確率の問題を「ありとあらゆる手段」で理解してしまおうという試みを、ここしばらく行なっている。
問題の概要を再確認しておこう。十二個の玉が布袋のようなものに入っていて、そのうちの4つが赤玉。この袋から一つずつ玉を取り出して一列に並べるとき、赤玉が隣合わない確率を求めたい。
我々の考察では真正面から取り組むのを避け、まずは手始めに（簡単そうに見えた）補事象の「赤玉が隣り合う確率」から考察を開始した。最初の成果として解析的な表現の導出に成功、次いでpythonなどによるプログラミングを用いた数値実験によりその結果も確かめた。ただ数値実験で確かめたのは、"event-A"と名付けた事象（赤玉２つが並んでいるとき、それが最初の赤玉である場合）のみである。
今回は"event-B"、つまり連続する赤玉が最初ではなく、その前に一回だけ赤玉がすでに取り出されたものの連続しなかった場合についての数値実験を行いたい。先に行った数学的な考察の結果、その確率は$112/495=0.226262...$であると思われるが、それが正しいかどうか検算するのである。
その前に計算の高速化の試み
event-Bの数値実験を行う前に, いままでawkでおこなってきた事象判定の高速化に取り組もうと思う。前回の試みでは、110万回の試行の判定にawkでは3秒ほどを要した。これをコンパイラーであるC言語で書き直し、高速化しようと思うのである。C言語を選んだのは、awkの文法がC言語によく似ているからである。FORTRANでもやってみたいのだが、別の機会に譲ることにした。
まずはaka-dama.pyで110万個のsequenceを生成し、aka-dama-sequence.datというファイルに結果を書き込む。そのファイルをリダイレクトして、aka-dama-count.exeと名付けたCプログラムに流し込んで処理をさせる方法を採用した。
C言語ではfgetsなどを利用して外部ファイルを開く方法もあるのだが、ファイル名がプログラム中に固定されてしまうのが面倒である（工夫をすれば避けることはできるが初心者には面倒...）。scanfを用いて外部ファイルからリダイレクトできるようにすると（ファイル名をプログラム中で気にしなくて良いので）便利である。
scanfの返り値はファイルエラーの場合などに利用されるので、調べてみた。man scanfによると、
RETURN VALUES
These functions return the number of input items assigned. This can be fewer than provided for, or even zero, in the event of a matching failure. Zero indicates that, although there was input available, no conversions were assigned; typically this is due to an invalid input character, such as an alphabetic character for a `%d' conver- sion. The value EOF is returned if an input failure occurs before any conversion such as an end-of-file occurs. If an error or end-of-file occurs after conversion has begun, the number of conversions which were successfully completed is returned.
scanfはint型の関数として定義され、成功すれば読み込んだ項目数を返し、失敗したときにはEOFが返るという。「失敗」とはつまりファイルが終わって読み込むものがなくなったとき、という意味であろう。
赤玉、黒玉、白玉からなるsequenceの長さは１２である（12個の玉を１列に並べる）から、文字列の終端を表す\0を含めサイズ13の文字型配列を用意し、そこに格納することにする。
#define NSIZE 13 .... char buffer[NSIZE]; .... while ( scanf("%s", buffer) != EOF ){ ... }
あとは、判定条件をコマゴマ丁寧に書き出すだけ（下のプログラムは主要部分の抜粋のみ）。
.... #define NSEQ 13 #define NDATA 1100000 int main(){ char seq[ NSEQ ]; int nsum[NSEQ], tsum; int i, nc=0; for(i =0; i< NSEQ; ++i){ nsum[i] = 0; } while ( scanf("%s", seq) != EOF ){ if(seq[0] == 'X' && seq[1] == 'X') { ++nsum[0]; }else if ( seq[0] != 'X' && seq[1] == 'X' && seq[2] == 'X'){ ++nsum[1]; }else if ( seq[0] != 'X' && seq[1] != 'X' && seq[2] == 'X' && seq[3] == 'X'){ ++nsum[2]; }else if ( seq[0] != 'X' && seq[1] != 'X' && seq[2] != 'X' && seq[3] == 'X' && seq[4] == 'X'){ ++nsum[3]; }else if ( seq[0] != 'X' && seq[1] != 'X' && seq[2] != 'X' && seq[3] != 'X' && seq[4] == 'X' && seq[5] == 'X'){ ++nsum[4]; }else if ...... ...... ...... }else{ } ++nc; }
gcc -o aka-dama-count.exe aka-dama-count.c
実行ファイルを実行し、aka-dama.pyの結果をリダイレクトで飲み込ませる。
aka-dama.py > aka-dama-sequence.dat; time aka-dama-count.exe < aka-dama-sequence.dat nc, total trial number : 1100000 1100000 (1,2) 99663 0.0906027 (2,3) 80181 0.0728918 (3,4) 62223 0.0565664 (4,5) 46364 0.0421491 (5,6) 33466 0.0304236 (6,7) 22072 0.0200655 (7,8) 13345 0.0121318 (8,9) 6650 0.0060455 (9,10) 2211 0.0020100 ---- Total probability: 0.3328864 real 0m0.140s user 0m0.124s sys 0m0.010s
awkで3秒かかっていた処理が、0.15秒で終わっている...やはりコンパイラは圧倒的に速い！
event-Bの判定プログラム
今度は、赤玉が２つ連続する前に、一度赤玉が単独で出てしまった事象の確率を数値計算する。
上のプログラムでは文字Xに対してif文を延々と書き下して判定させているが、いささか見苦しく、見にくくてデバッグしにくい感じがする。 そこで、文字列そのものを処理するのではなく、数字に置き換え、代数処理ができるよう工夫することを考える。赤玉Xは1に、その他の非赤玉(O,@)は負の値に変換してから別の配列に格納し直す（今回は白玉Oは-10に、黒玉@は-100にした）。
#define NSEQ 13 int main(){ .. .. char seq[NSEQ]; int kseq[NSEQ]; while ( scanf("%s", seq) != EOF ){ for(i=0; i < NSEQ; ++i){ switch(seq[i]){ case 'X': kseq[i] = 1; break; case 'O': kseq[i] = -10; break; case '@': kseq[i] = -100; break; default: kseq[i] = 0; break; } } .. ..
次はいよいよ事象の判定アルゴリズムである。まずn1回目とn1+1回目に赤玉が連続する条件を考える。プログラムではkseqという配列に数字化したsequenceが格納されており、赤玉の場合は1が入ることになっている。したがって赤玉が連続するならkseq[n1]=1, kseq[n1+1]=1であるから、その和は２になる。つまり、
if (kseq[n1] + kseq[n1+1] == 2)
が赤玉が連続する条件である。
event-Aの場合は、この連続赤玉が「最初赤玉」であるから、n1-1回までの配列には全て負の数が入っているはずである。だからといって、2つ目の条件を
sum0 = kseq[n1] + kseq[n1+1]; sum1=0.0; for(i=0; i <= n1-1; ++i){ sum1 += kseq[i]; } if(sum1 < 0 && sum2 == 2){ ... }
とやってはならない！なぜなら、赤玉の場合の配点数1は白玉の-10点、黒玉の-100点に比べると、小さすぎるからである。例えば、n1=4の場合に@XOXX....という列が得られたとしよう。最初の黒玉には-100点、二回目の赤玉には+1点、三回目の白玉には-10点が割り振られるが、その和は-109点となり「負の数」である。つまり、和が負であることは赤玉が出なかったことの印にはなり得ないのである（プログラムの最初のバージョンではこれにより数え間違いが発生してしまった.....）。 そこで利用するのが、赤玉の得点は一桁の数（つまり＋1点）という性質である。上の例でも-109点であるから、100で割ったときの余りが-9となる。つまり、赤玉が出ると一桁に対応する数字が非零になるのである。この性質を使うと、判定は次のように改善される。
sum0 = kseq[n1] + kseq[n1+1]; sum1=0.0; for(i=0; i <= n1-1; ++i){ sum1 += kseq[i]; } sum1 = (sum1 % 100) % 10; if(sum1 == 0 && sum2 == 2){ event-Aが発生したと認識できる }
続けてevent-Bの判定条件を考えよう。
まずn1-1回目が非赤でなくてはならない。もし赤ならば赤玉三回連続となるが、その事象は$p_{n1-1, n1}$に寄与すべきものである。「数えすぎ」になってしまう。そこでkseq[n1-1] < 0という条件を加えることになる。また、赤玉が一つ入ったとき、n1-1回までの得点の和をとると一桁の位が-9になるという性質がある。これを利用する。以上をまとめると、event-Bの判定条件は
sum0 = kseq[n1] + kseq[n1+1]; sum1=0.0; for(i=0; i <= n1-1; ++i){ sum1 += kseq[i]; } sum1 = (sum1 % 100) % 10; if(sum1 == -9 && sum2 == 2 && kseq[n1-1] < 0){ event-Bが発生したと認識できる }
さらにevent-Bはn1=3以上で許される事象であるから、その条件も組み込む必要がある。
ということで最終的なプログラムは次のようになった。
#include <stdio.h> #include <stdlib.h> #define NSEQ 13 #define NDATA 1100000 int main(){ char seq[NSEQ]; int kseq[NSEQ]; int nAsum[NSEQ], nBsum[NSEQ]; int tAsum, tBsum, rrsum,rcsum, ires; int i, j, n1, nc=0; for(i =0; i< NSEQ; ++i){ nAsum[i] = 0; nBsum[i] = 0; kseq[i] = 0; } while ( scanf("%s", seq) != EOF ){ for(i=0; i < NSEQ; ++i){ switch(seq[i]){ case 'X': kseq[i] = 1; break; case 'O': kseq[i] = -10; break; case '@': kseq[i] = -100; break; default: kseq[i] = 0; break; } } for(n1=0; n1 < 9; ++n1){ rrsum = kseq[n1] + kseq[n1+1]; rcsum = 0; for(j=0; j < n1; ++j){ rcsum += kseq[j]; } ires = (rcsum % 100) % 10; // Case A if(rrsum == 2 && ires == 0){ ++nAsum[n1]; } // Case B if(n1 > 1 && rrsum == 2 && ires == -9 && kseq[n1-1] < 0){ ++nBsum[n1]; } } ++nc; } printf("nc, total trial number : %d %d\n\n", nc, NDATA); tAsum = 0; tBsum = 0; for(i=0; i< 9; ++i){ printf("(%d,%d) %d %d %10.7f %10.7f \n",i+1,i+2,nAsum[i], nBsum[i], 1.0*nAsum[i]/NDATA, 1.0*nBsum[i]/NDATA); tAsum += nAsum[i]; tBsum += nBsum[i]; } printf("----\n"); printf("Total probability: %10.7f %10.7f %10.7f \n", 1.0*tAsum/NDATA, 1.0*tBsum/NDATA, 1.0*(tAsum+tBsum)/NDATA); return EXIT_SUCCESS; }
コンパイルして実行した結果が以下である。
nc, total trial number : 1100000 1100000 (1,2) 99663 0 0.0906027 0.0000000 (2,3) 80181 0 0.0728918 0.0000000 (3,4) 62223 17718 0.0565664 0.0161073 (4,5) 46364 30845 0.0421491 0.0280409 (5,6) 33466 39712 0.0304236 0.0361018 (6,7) 22072 44136 0.0200655 0.0401236 (7,8) 13345 44569 0.0121318 0.0405173 (8,9) 6650 40204 0.0060455 0.0365491 (9,10) 2211 31075 0.0020100 0.0282500 ---- Total probability: 0.3328864 0.2256900 0.5585764 real 0m0.674s user 0m0.449s sys 0m0.028s
計算結果は次のようにして解釈する。 まず最初の行は、sequenceは(pythonによって）110万個生成され、（Cプログラムによって）110万個が分析された、ことを表す。 ２行目以降は各々の確率$p_{n,n+1}$についてのデータが印字されている。かっこの中が$(n1,n1+1)$を意味し、次の2つの（大きな）整数がevent-Aとevent-Bに対応する発生数を表し、最後の２つの少数がevent-Aとevent-Bに対応する（このサブ事象の）確率を表す。 (9,10)の後の行で"----"というセパレーターが印字され、各事象についての確率の和がevent-A, event-B、そして両者の和に対して計算され印字される。 real, user, sysは、この計算に要した計算時間である。結果は0.67秒であった。
event-Bの数値実験の値は0.22569であるのに対し、考察によって求めた値は$112/495=0.226262...$つまり、小数位３桁までが一致している。event-Bの各確率についてもなかなかよい精度が出ているし、$p_A+p_B=277/495=0.5595...$の理論値に対し、数値実験は0.55857..という値を返している（小数位３桁まで一致）。
ということで、数値計算も理論計算も、大丈夫そうな感じである！
さて次はいよいよ問題の本丸へと切り込むことにしよう。
