東京大学2023数学[2]-(1) : 赤玉確率の問題 - 複天一流：どんな手を使ってでも問題解決を図るブログ

確率の問題に挑戦
12個の玉を袋から取り出し並べるプログラム
- 小学生バージョン
- 高校生バージョン
大数の法則
最初に赤玉が2回連続する確率

確率の問題に挑戦

東京大学の数学の入試問題に挑戦してみたい。といっても、前回と同じように正攻法ではなく、搦め手である。使えるものは全て使って解決を目指す。

今回の問題は確率である。コンピュータを使った解決法にはうってつけである。シミュレーションとか数値実験とか、そういう類の試みからまずは始めてみたい。一方で、正攻法で解くのはなかなか厄介である。そもそもこんな問題が軽々解けるような人間だったら、今頃こんなブログを書いているはずはないのである。

とはいえ、確率の問題というのは、日頃からかなり練習を積んでおかないと、案外意外なところでの数え落としがあるものである。だからこそコンピューターの愚直さというものが我々には必要なのである。あの藤井７冠でも日頃はコンピュータに教えを請うて、数え落としについて修行を積んでいるそうである。どんな優秀な人間でもミスって数え落とすものであることは覚えておくべきであろう。特に原発管理者の諸君には声を大にして呼びかけたい。2011年の大地震が来るまで「絶対に私は事故らない」と断言していた原発関係の工学者たちは、ここ10年ほどは少しは反省していたように見えたのだが、最近だんだん昔のようになってきたのをちょっと感じる。今回のような確率の問題は彼らの十八番なんだろうか？「私は絶対に数え落とさない」とか今でも彼らは断言できるものなのであろうか？

ちなみに私はよく数え落としをやってしまう。原発管理者になったとすれば、「あっ！」とか言った直後に世界を破滅させる自信がある。

さて問題文を見てみよう。

袋から赤玉、黒玉、白玉を取り出すという古典的な博打である。連続して赤玉を出さない確率を求めよ、というが５０％より高いのか、それとも低いのか、まずは簡単に考えてみよう。

袋の中に玉は全部で12こ入っている。赤玉は4つしか入ってない。全体の1/3である。最初に赤玉を引き当てる確率がこれにあたる。２回目に赤玉を引き当てる確率は3/11となる。したがって、１回目と２回目に連続して赤玉を引き当てる確率は \begin{equation} \frac{1}{3}\cdot\frac{3}{11} = \frac{1}{11}\end{equation}であろう。私は私の数え上げの能力を信じていないので、数え落としはないかマシンに確認させてみたい。

12個の玉を袋から取り出し並べるプログラム

まずは黒玉３つ、赤玉４つ、そして白玉５つが入った袋から、一つずつ取り出し横一列に並べるプログラムを書いてみよう。

重要なのは乱数の生成である。完全な乱数の生成はかなり難しい問題だが、今回の問題については普通程度の乱数で多分十分であろう。答えがわかっている関数についてモンテカルロ積分をやってみれば、どの程度ランダムな乱数になっているかはある程度把握できる（が、その確認は別の機会にやることにする）。

今回もpythonを利用しよう。pythonにおける乱数生成の方法は多種あるらしいが、もっとも標準で簡単なものがrandomライブラリだそうである。プログラムの文頭で読み込めば使えるようになる。

import random

乱数を生成する関数はrandom()なので次のようにして呼び出す。

random.random()

この乱数生成関数は0から1までの浮動小数の形で乱数を生成する。

我々が欲しいのは1から12までの自然数の乱数である。pythonでは配列などは０からスタートすることを考慮して、1から12の代わりに、0から11までの自然数を出力する関数が欲しいところである。その目的のためには、普通次のようにする。

 int( random.random() * 12 )

intというのは整数化のための関数で、少数点以下を切り捨てる機能を持つ。mathとか呼ばずにそのままデフォルトで利用できる。

ASCIIで玉の色を表現することにする。Xを赤玉、@を黒玉、Oを白玉とみなす。

pythonでプログラムを書くと次のような感じになった。

#!/usr/local/bin/python3

# Aka-dama 4, Kuro-dama 3, Shiro-dama 5
#
# ver.0 (Oct.6, 2023): how to use the random-function


import random
NMAX_DM = 12
NMAX_TR = 110


for i in range(NMAX_TR):
    dmbox = [1,1,1,1,1,1,1,1,1,1,1,1]
    while True:    
        x = int( random.random() * NMAX_DM ) 
        if dmbox[x] == 1:
            dmbox[x] = 0

            if x <= 3:
                print("X",end="")
            elif x >= 4 and x <= 6:
                print("@",end="")
            elif x >= 7:
                print("O",end="")
        sum = 0
        for i in range(12):
            sum += dmbox[i]
        if sum == 0:
            print("")
            break

このプログラムでは、赤玉が"X"である。白玉、黒玉はそれぞれ"@"と"O"で表現する。デフォルトの設定では、袋から取り出すという事象(NMAX_TR)を110回繰り返させている。もし私の考察が正しく、先頭の２つがXXとなる確率が1/11であるならば、この数値実験では10回程度の事象が引っかかるはずである。

とはいえ、いきなり110回の計算をやって、そのうちの10回分程度の事象が該当するかどうかを（目でみて）調べるのはきつい作業なので、いろいろと便利な道具をかき集めて使ってみることにする。

まずは「お試し計算」である。デバッグの目的も兼ね（プログラムの動作確認）、手計算で確認できる程度の試行回数に落として実行する。しばらくはNMAX_TR=11でやってみる。

袋からの取り出しを11回行う試行をプログラムにやらせたとき、予想としては1回程度の「ヒット」があるはずである。実行結果は、例えば次のようになった。

aka-dama.py | cat -n

     1  X@X@OOO@XOOX
     2  XOOXX@OOO@X@
     3  @OOX@OXXO@OX
     4  O@OX@X@OXOOX
     5  OX@XXOOOX@@O
     6  XO@OXOO@XOX@
     7  OOX@@OXXOO@X
     8  OOXO@X@X@XOO
     9  XXXO@@O@OOOX
    10  OOX@XO@X@OOX
    11  XXOXX@O@OO@O

9番と11番の2回が該当する事象である。この実験の場合は2/11の確率となる。

該当する事項がいくつあるかを識別するためにawkを使ってみる。任意の文字の切り出しにはsubstr関数を使うと良いらしい。

今やりたいのは最初とその次の文字がXかどうかである。どうやって実装するかというと、

x1 = substr($2, 1, 2);

とやるそうである。今考えているpythonのプログラムでは１行に2つの項目の出力があり、最初が通し番号（整数）で、2つ目が玉の色を表すX,@,Oからなる１２文字の文字列である。substrの対象は2つ目の項目（文字列）の方だから、awkでは$2を指定し、そこに処理をかける。先頭から２文字取り出すので1,2という指定となる。もし3文字目から４文字分取り出すときは3,4となる。

pythonは先頭の項目を「０番目」として数える定義であったが、awkでは「１番目」として数え始める定義である。

早速awkを組み込んでみよう。

aka-dama.py | cat -n | awk '{if( substr($2,1,2) == "XX" ) print $1,$2;}'

上の試行結果に対して実施した結果が以下である。

9 XXXO@@O@OOOX
11 XXOXX@O@OO@O

ちゃんと該当する行が印字された。

何回か上のコマンドを繰り返してみる。

(最初の試行ではnull return、つまり該当するものなし）

（2回目）
1 XXOO@O@OOXX@
2 XXOO@XOOOX@@
7 XX@OXO@O@OOX
8 XXO@OOX@@OXO

（3回目）
null

（4回目）
null

（5回目）
8 XXO@XOOX@OO@

（6回目-8回目）
null

（9回目）
4 XXO@XOO@@XOO
9 XX@X@OOOXO@O

（10回目）
7 XXOO@X@OOX@O

であった。1回のセットにつき11回の試行を行う計算を、10回セット行なった場合の確率計算方法はいろいろあるが、まずは下準備をしておこう。

まずは各セットでヒットした回数を書き出す。

セット番号	1	2	3	4	5	6	7	8	9	10
ヒット数	0	4	0	0	1	0	0	0	2	1

次に、この表をヒストグラムに変換

ヒット数	0	1	2	3	4
頻度	6	2	0	0	1

各ヒット数に相当する確率は、当然ながら0, 1/11, 2/11, 3/11, 4/11 となり、期待値は1/11になるはずである（....）。

小学生バージョン

まずは小学生バージョンで計算してみよう。トータルのヒット数は0回x６＋1回x２＋2回x０＋3回x０＋4回x１= 6回。1回のセットで11回の試行を行い、１０セット行なったので、トータルの試行回数は11x10=110回。したがって、赤玉が最初に連続する平均の確率は（この実験に関しては） \begin{equation} \frac{6}{110} = \frac{3}{55} \simeq 0.054\cdots < 1/11 (\simeq 0.09\cdots)\end{equation}という結果となった。

ただし、それぞれのセットにおける確率は0, 4/11, 0, 0, 1/11, 0, 0, 0, 2/11, 1/11となっている。お目当の1/11は2回だけ発生している。とはいえ、再頻出は「ヒット０回」であり、望まれる「ヒット1回」とは「誤差１」に過ぎない。小学生バージョンで確率を計算すると、予想される値と大きく違ってしったのは、この「０回」の存在である。「惜しい」結果なのに０という値を持つがゆえに、確率の計算から弾かれてしまうのだ。

そこでヒット数の期待値を計算してみることにする。10セットやって、0が6回、1が2回、4が1回だから、(1x2 + 4x1)/10 = 0.6である。切り捨てれば０だが、四捨五入すれば１であるから、$\frac{\langle 1\rangle}{11}$となって期待通りの結果になる。０と１の差が小さいようで大きいのが問題なのである。

高校生バージョン

高校数学になると各セットにおける確率$p_i$と頻度$n_i$を使って平均の確率や期待値を計算する別の公式を習う。平均の確率は\begin{equation} \langle p \rangle = \frac{\sum_{i=0}^N n_i \cdot i}{\sum_{i=0}^N 11}= \frac{1}{\sum_{i=0}^N 1} \sum_{i=0}^N \left(\frac{n_i \cdot i}{11}\right) = \frac{\sum_{i=0}^N n_i p_i}{N}\end{equation}

最後の表式における$n_i/N$は「確率の重み」とか「確率分布」と呼ばれる。$w_i=n_i/N$と書くと、 \begin{equation} \langle p \rangle = \sum_{i=0}^N w_ip_i\end{equation}となる。

我々の数値実験において、期待されるヒストグラムは、ヒット数が1回の項目に大きなピークが来て、その周辺にわずかに揺らぐというものである。しかし、試行回数が少ないため、データの揺らぎが大きくなってしまっている。ピークの位置がわずかに「左」にずれたがために、1回→0回にピークの位置が移動してしまったのだ。たまたま逆に揺らいで1回→2回となれば、もう少し理想に近い形になったかもしれないが、たまたま今回はそうならなかった。

試行回数を増やして確率分布を予想されるものに近づける必要がある。この性質を「大数の法則」という。証明は確率統計の教科書に（たしか）載っていたと思う（.....汗）。

ということで、セット数を一気に増やして、大数の法則が成立するようにしてみよう。

大数の法則

pythonのプログラムはいじらず（つまりNMAX_TR=11のまま）、コマンドの実行回数を100にしてみよう。シェルスクリプトで実装する。今回はbashを使ってみる。

まず、与えられたセットの中で該当する事象が何回起きたかを判別するためにwcを利用する。行数を数えることでヒット数を判定するのだ。

aka-dama.py| awk '{if(substr($1,1,2) == "XX") print $1}'| wc

     2       2      26

表示は（１）行数（２）語数（３）文字数である。ヒット数が一つあると、玉の識別文字１２文字＋改行で１３文字になっている。

この命令をコマンドに100回打ち込むのは、かなり厳しい「修行」でしかない。そこで新たにスクリプトを書くことにする。

#!/bin/bash                                                                                                                           
NMAX_TR=100
i=0
while [ 1 ]
do
    aka-dama.py | cat -n |awk '{if(substr($2,1,2) == "XX") print $1,$2}'|wc
    i=$((i+1))
    if [ $i -ge ${NMAX_TR} ]; then
        break
    fi
done

これで自動的に100回計算してくれるようになった。このシェルスクリプトはre-aka-dama.shと名付ける。

ヒストグラムも自動的に生成してしまおう。awkを噛ませて実現しよう。ただ、ちょっと長めの命令になるのでスクリプト化することにする。

BEGIN{n0=0; n1=0;n2=0;n3=0;n4=0; n5=0;}

{
    if($1 == 0){
        ++n0;
    }else if($1 == 1){
        ++n1;
    }else if($1 == 2){
        ++n2;
    }else if($1 == 3){
        ++n3;
    }else if($1 == 4){
        ++n4;
    }else{
        ++n5;
    }
}

END{print n0,n1,n2,n3,n4,n5;}

switch() -- caseを使おうと思ったがエラーが出るので、どうもawkでは使えないようである。このスクリプトはhistgr.awkと名付ける。

re-aka-dama.sh| awk -f histgr.awk 

34 41 21 3 1 0

となった。なかなか良さそうである。

ちょっと時間はかかるが、200回の繰り返しをやってみると

76 74 36 12 2 0

となったが、また「ヒット0回」が最頻出になってしまい、確率計算の観点からは微妙な感じである。

一気に500回に増やしてみた。私の環境では30秒もかかったが、次の結果を得た。

178 181 103 31 7 0

さらに倍の1000回と三倍の3000回に挑戦である（それぞれ43秒と125秒だった）。

331 395 204 57 13 0
1035 1167 581 184 27 6

最後の６という数字が不味い感じがするが、これは5/11の確率の回数だと思って計算する（この後でプログラムを改良し５文字連続の事象も扱えるようにした）。さらに追加として6000回の計算をやってみたが、4分14秒かかった。

2126 2270 1167 355 66 16 0

さて、平均の確率を計算してみよう。

回数	100	200	500	1000	3000	6000
確率	0.16636	0.08636	0.09236	0.09327	0.0906	0.09111

となった。

$1/11 = 0.090909\cdots$なので、次第に大数の法則が効き始めているように見えなくもないが、収束は悪そうである。これは乱数の系列が偏っているせいなのかもしれない。

最初に赤玉が2回連続する確率

ということで、大数の法則を込みの上で、我が考察「赤玉が最初に2回連続する確率は1/11」については、数値実験により「まあ証明できた」といってよいのではないだろうか？数学の証明に対して「まあいいか」という態度をとるとはけしからん、と思う人も多いだろうが、ここは「理解したものがち」をテーマにやっているので、なりふり構わず全ての道具を使って「我が理解」に向かって突進するのである。

次回は、数値計算によって、入試問題の解答の候補を挙げてみたいと思う。（つづく）