ビンゴの確率その2 - ヒビノキロク

id:nozom:20061229:1167403674の続き。

前回はm個の穴が空いたカードにビンゴが１つ以上ある確率q(m)を求める手前で終わっていた。これを数式で出すのは難しかったので、以下のようなプログラムで全ての場合を数えるという非常に安直な方法をとった。

#!/usr/bin/env gosh

(use srfi-1)
(use util.combinations)
(use util.list)

(define (count-permutation n m)
  (let rec ((n n) (m m) (num 1))
    (if (positive? m)
        (rec (- n 1) (- m 1) (* num n))
        num)))

(define (count-combination n m)
  (let rec ((n n) (m m) (numer 1) (denom 1))
    (if (positive? m)
        (rec (- n 1) (- m 1) (* numer n) (* denom m))
        (/ numer denom))))

(define (bingo-set n)

  (define (pos-to-id x y)
    (+ (* n y) x))

  (define (center-topleft)
    (let loop ((i (/ (- n 1) 2)) (r ()))
      (if (> i 0)
          (loop (- i 1) (cons (pos-to-id (- i) (- i))
                              (cons (pos-to-id i i) r)))
          r)))

  (define (center-topright)
    (let loop ((i (/ (- n 1) 2)) (r ()))
      (if (> i 0)
          (loop (- i 1) (cons (pos-to-id i (- i))
                              (cons (pos-to-id (- i) i) r)))
          r)))

  (define (center-horizontal)
    (let loop ((i (/ (- n 1) 2)) (r ()))
      (if (> i 0)
          (loop (- i 1) (cons (pos-to-id i 0)
                              (cons (pos-to-id (- i) 0) r)))
          r)))

  (define (center-vertical)
    (let loop ((i (/ (- n 1) 2)) (r ()))
      (if (> i 0)
          (loop (- i 1) (cons (pos-to-id 0 i)
                              (cons (pos-to-id 0 (- i)) r)))
          r)))

  (define (horizontal-line y)
    (let loop ((i n) (r ()))
      (if (> i 0)
          (let ((x (- i (/ (+ n 1) 2))))
            (loop (- i 1) (cons (pos-to-id x y) r)))
          r)))

  (define (horizontal-lines)
    (let loop ((i (/ (- n 1) 2)) (r ()))
      (if (> i 0)
          (loop (- i 1) (cons (horizontal-line (- i))
                              (cons (horizontal-line i) r)))
          r)))

  (define (vertical-line x)
    (let loop ((i n) (r ()))
      (if (> i 0)
          (let ((y (- i (/ (+ n 1) 2))))
            (loop (- i 1) (cons (pos-to-id x y) r)))
          r)))

  (define (vertical-lines)
    (let loop ((i (/ (- n 1) 2)) (r ()))
      (if (> i 0)
          (loop (- i 1) (cons (vertical-line (- i))
                              (cons (vertical-line i) r)))
          r)))

  (append (list (center-topleft)
                (center-topright)
                (center-horizontal)
                (center-vertical))
   (horizontal-lines)
   (vertical-lines)))

(define (make-bingo? n)
  ;; (define (cmpfn set1 set2)
  ;;   (lset<= = set2 set1))
  (define (cmpfn set1 set2)
    (or (null? set2)
        (let ((obj (memv (car set2) set1)))
          (and obj
               (cmpfn (cdr obj) (cdr set2))))))
  (let ((s (map sort (bingo-set n))))
    (lambda (set)
      (let ((n (count (cut cmpfn (sort set) <>) s)))
        (and (> n 0) n)))))

(define (bingo-holes n)
  (let loop ((i (/ (- (* n n) 1) 2)) (lst ()))
    (if (zero? i)
        (sort lst)
        (loop (- i 1) (cons i (cons (- i) lst))))))

(define (display-bingo-table n s e)
  (define bingo? (make-bingo? n))
  (let ((holes (bingo-holes n)))
    (define (count-bingo i)
      ;; (count bingo? (combinations holes i))
      (let ((count 0) (ht (make-hash-table)))
        (combinations-for-each (lambda (set)
                                 (let ((k (bingo? set)))
                                   (when k
                                     (set! count (+ count 1))
                                     (hash-table-update! ht k (cut + 1 <>) 0))))
                               holes i)
        (sort (hash-table->alist ht) (lambda (a b) (< (car a) (car b))))))
    (let loop ((i s))
      (if (<= i e)
          (begin
            (display i)
            (display " ")
            (display (cons (count-combination (- (* n n) 1) i)
                           (count-bingo i)))
            (newline)
            (loop (+ i 1)))))))

(define (main args)
  (if (> (length args) 1)
      (let ((n (string->number (list-ref args 1))))
        (display-bingo-table n
                             (if (> (length args) 2)
                                 (string->number (list-ref args 2))
                                 (- n 1))
                             (if (> (length args) 3)
                                 (string->number (list-ref args 3))
                                 (- (* n n) 1))))
      (begin
        (display (list-ref args 0))
        (display " n [s] [e]")
        (newline))))

これを実行すると、m個の穴の組み合わせでビンゴになる場合の数が以下のように得られる(M=24つまり5x5の場合)。これは例えばm=8のとき、全部で735471通りの組み合わせのうち１列ビンゴになるのが27072通り、２列ビンゴになるのが30通りということを表している。

% ./bingo.scm 5
4 (10626 (1 . 4))
5 (42504 (1 . 88))
6 (134596 (1 . 912))
7 (346104 (1 . 5928))
8 (735471 (1 . 27072) (2 . 30))
9 (1307504 (1 . 92016) (2 . 504))
10 (1961256 (1 . 240120) (2 . 3972))
11 (2496144 (1 . 488304) (2 . 19344) (3 . 48))
12 (2704156 (1 . 775944) (2 . 64428) (3 . 728))
13 (2496144 (1 . 955248) (2 . 153096) (3 . 5008) (4 . 8))
14 (1961256 (1 . 891852) (2 . 262344) (3 . 20188) (4 . 236))
15 (1307504 (1 . 607832) (2 . 319792) (3 . 51760) (4 . 2032) (5 . 8))
16 (735471 (1 . 283724) (2 . 266280) (3 . 85696) (4 . 8563) (5 . 180) (6 . 2))
17 (346104 (1 . 81448) (2 . 139856) (3 . 88368) (4 . 20024) (5 . 1336) (6 . 24))
18 (134596 (1 . 11868) (2 . 39824) (3 . 51380) (4 . 25568) (5 . 4564) (6 . 220) (7 . 4))
19 (42504 (1 . 584) (2 . 4536) (3 . 13456) (4 . 15664) (5 . 7128) (6 . 1064) (7 . 48))
20 (10626 (2 . 90) (3 . 922) (4 . 3204) (5 . 4136) (6 . 1982) (7 . 282) (8 . 10))
21 (2024 (4 . 56) (5 . 456) (6 . 872) (7 . 544) (8 . 96))
22 (276 (6 . 8) (7 . 100) (8 . 132) (9 . 36))
23 (24 (9 . 8) (10 . 16))
24 (1 (12 . 1))

この結果を元に、次のような表ができる。

m	${}_{M}C_{m}$	bingo(m)	q(m)	dq(m)
4	10626	4	0.0004	0.0004
5	42504	88	0.0021	0.0017
6	134596	912	0.0068	0.0047
7	346104	5928	0.0171	0.0104
8	735471	27102	0.0368	0.0197
9	1307504	92520	0.0708	0.0339
10	1961256	244092	0.1245	0.0537
11	2496144	507696	0.2034	0.0789
12	2704156	841100	0.3110	0.1076
13	2496144	1113360	0.4460	0.1350
14	1961256	1174620	0.5989	0.1529
15	1307504	981424	0.7506	0.1517
16	735471	644445	0.8762	0.1256
17	346104	331056	0.9565	0.0803
18	134596	133428	0.9913	0.0348
19	42504	42480	0.9994	0.0081
20	10626	10626	1.0000	0.0006
21	2024	2024	1.0000	0.0000
22	276	276	1.0000	0.0000
23	24	24	1.0000	0.0000
24	1	1	1.0000	0.0000

ただし、左から順にm、M個の穴からm個を選ぶ組み合わせの数 ${}_{M}C_{m}$ 、そのうち１つ以上ビンゴになっている場合の数 $\mbox{bingo}(m)$ 、ビンゴが１つ以上ある確率 $q(m)=\mbox{bingo}(m)/{}_{M}C_{m}$ 、最初のビンゴになる確率 $\mbox{d}q(m)=q(m)-q(m-1)$ 。

これでq(m)が求まったので、前回漸化式の形で表したp(n,m)も実際に値を出して、r(n)を計算してみる。
結果だけ書くと、5x5のビンゴで玉の数が100個の場合こんな風になる。

n	r(n)	dr(n)	1/r(n)
4	0.0000	0.0000	980306.250
5	0.0000	0.0000	192060.000
6	0.0000	0.0000	62739.600
7	0.0000	0.0000	26361.176
8	0.0001	0.0000	12927.142
9	0.0001	0.0001	7046.299
10	0.0002	0.0001	4149.572
11	0.0004	0.0001	2592.729
12	0.0006	0.0002	1697.738
13	0.0009	0.0003	1154.864
14	0.0012	0.0004	810.813
15	0.0017	0.0005	584.653
16	0.0023	0.0006	431.317
17	0.0031	0.0008	324.558
18	0.0040	0.0009	248.495
19	0.0052	0.0012	193.197
20	0.0066	0.0014	152.270
21	0.0082	0.0017	121.493
22	0.0102	0.0020	98.018
23	0.0125	0.0023	79.879
24	0.0152	0.0027	65.700
25	0.0183	0.0031	54.498
26	0.0219	0.0036	45.561
27	0.0261	0.0041	38.367
28	0.0307	0.0047	32.528
29	0.0360	0.0053	27.753
30	0.0420	0.0059	23.821
31	0.0486	0.0067	20.560
32	0.0561	0.0074	17.840
33	0.0643	0.0082	15.559
34	0.0733	0.0091	13.634
35	0.0833	0.0100	12.002
36	0.0942	0.0109	10.612
37	0.1061	0.0119	9.423
38	0.1190	0.0129	8.401
39	0.1330	0.0140	7.519
40	0.1480	0.0150	6.756
41	0.1641	0.0161	6.092
42	0.1814	0.0172	5.513
43	0.1997	0.0183	5.007
44	0.2192	0.0194	4.563
45	0.2397	0.0205	4.172
46	0.2613	0.0216	3.827
47	0.2839	0.0226	3.522
48	0.3076	0.0236	3.251
49	0.3322	0.0246	3.011
50	0.3576	0.0254	2.796
51	0.3838	0.0262	2.605
52	0.4108	0.0269	2.434
53	0.4383	0.0276	2.281
54	0.4664	0.0281	2.144
55	0.4948	0.0284	2.021
56	0.5235	0.0287	1.910
57	0.5524	0.0288	1.810
58	0.5812	0.0288	1.721
59	0.6098	0.0287	1.640
60	0.6382	0.0284	1.567
61	0.6661	0.0279	1.501
62	0.6935	0.0273	1.442
63	0.7201	0.0266	1.389
64	0.7459	0.0258	1.341
65	0.7707	0.0248	1.298
66	0.7944	0.0237	1.259
67	0.8169	0.0225	1.224
68	0.8381	0.0212	1.193
69	0.8580	0.0199	1.166
70	0.8764	0.0185	1.141
71	0.8934	0.0170	1.119
72	0.9089	0.0155	1.100
73	0.9230	0.0141	1.083
74	0.9356	0.0126	1.069
75	0.9468	0.0112	1.056
76	0.9566	0.0098	1.045
77	0.9650	0.0085	1.036
78	0.9723	0.0073	1.028
79	0.9784	0.0061	1.022
80	0.9835	0.0051	1.017
81	0.9876	0.0041	1.013
82	0.9909	0.0033	1.009
83	0.9935	0.0026	1.007
84	0.9955	0.0020	1.004
85	0.9970	0.0015	1.003
86	0.9981	0.0011	1.002
87	0.9988	0.0007	1.001
88	0.9993	0.0005	1.001
89	0.9996	0.0003	1.000
90	0.9998	0.0002	1.000
91	0.9999	0.0001	1.000
92	1.0000	0.0001	1.000
93	1.0000	0.0000	1.000
94	1.0000	0.0000	1.000
95	1.0000	0.0000	1.000
96	1.0000	0.0000	1.000
97	1.0000	0.0000	1.000
98	1.0000	0.0000	1.000
99	1.0000	0.0000	1.000
100	1.0000	0.0000	1.000

最初dr(n)の期待値を出してみたらE(dr(n))=54.98となったんだけど、これって結局のところ中心極限定理で1-100の平均近くになるというだけの話で何とも当たり前の結果になってしまった。それで試しにr(n)の逆数を取ってみたところ、今度はなかなか面白い数字が出た。この数字は要するに何分の１の確率でビンゴになるかを表していて、つまり会場にこれだけの人数がいれば１人ぐらいビンゴになる、という風に解釈できる。例えばn=30ぐらいで1/r(n)=25となるので、20〜30人でビンゴをやるとだいたい30回ぐらいで誰かがビンゴになる計算になる。逆に、100万人ぐらい集まればn=4つまり最短でビンゴになる人も出てくるというわけだ。

ちなみに、このページではモンテカルロシミュレーション *1によってビンゴの確率を求めていて、それと比較してこの結果はまずまず妥当なものだと言えるだろう。

今日はここまで。次回は、プログラムで数えて求めたq(m)を後付けで数式によって表現してみる予定(m=12まではできてる）。

似たようなことを考えた人たち

*1:全然関係ないけどモンテカルロシミュレーションから辿ったマルコフ連鎖モンテカルロのキーワードページにある説明がやたら詳しくて驚いた