表題の質問を生徒から受けました。簡単だと思っていたら、意外と手こずりました。
イメージしやすいように具体的な例で言うと、3個の箱に4個のボールをでたらめに入れていく場合、すべての箱に1個以上のボールが入る確率を求めよということです。この場合は、m=3, n=4です。このくらいのm、nの数なら、まあ、なんとか解けるかと思います。ただし、一般のm, nについて定式化するのは難しいようです。
以下、3通りの方法で説明したいと思います。
方法１
乱数を使ってコンピューターシミュレーションしてみましょう。これは簡単です。
乱数によりボールをm個の箱のうち、どの箱に入れるかを決定します。これを次々にn回繰り返します。
そしてすべての箱に入ったからどうかをチェックします。
この作業を、例えば10万回繰り返して、そのうち何回がすべての箱に入ったかを求めます。
これで、確率の近似解を得ることができます。
Pythonプログラムは以下のようになります。
# ボールn個、箱m個のとき シミュレーション
import numpy as np
repeat = 100000 # number of trials
count = 0 # success counter
n = 9 # number of balls
m = 3 # number of boxes
for _ in range(repeat):
box = np.zeros(m)
for _ in range(n):
#box[np.random.randint(0, m)] += 1
box[np.random.randint(0, m)] = 1
if(np.sum(box) == m):
count += 1
print(count/repeat)
方法２
Pythonであれば、順列、組み合わせを簡単に扱うことができます。今回は重複あり順列を使います。ボールn個を箱m個に入れるすべてのケースをリストにします。リストの各要素には選ばれた箱が入っています。なので、リストの要素一個ずつについて、すべての箱が入っているかどうかをチェックすれば良いわけです。この方法では、正確な確率の値が求まります。ただし、n、mが大きくなると膨大な要素数のリストになり、長い計算時間が必要となります。
# ボールn個、箱m個のとき 重複あり順列で1つずつチェック
import itertools
import numpy as np
from fractions import Fraction
n = 8
m = 7
np.box = [i for i in range(m)]
np.all_case = list(itertools.product(np.box, repeat=n))
#print(np.all_case)
np.in_box_case = np.ones(len(np.all_case))
for i in range(len(np.box)):
np.in_box_case = np.in_box_case * [np.box[i] in np.all_case[j] for j in range(len(np.all_case)) ]
occupy = int(np.sum(np.in_box_case))
total = len(np.all_case)
p = Fraction(occupy, total)
print('occupy: {0}, total: {1}, p: {2}, p: {3}'.format(occupy, total, p, round(float(p), 3)))
方法３
求める確率を箱の数mとボールの個数nを使った式で書ければ一番良いわけですが、これがなかなか難しいです。
以下のプログラムではm=7まで対応しています。
考え方は、空き箱ができる場合を求めて、ボールの入れ方の総数m**nから引くことにより、空き箱のない場合の数を求めます。
空き箱を1つ以上作る場合の数は、残りの箱(m-1)個にボールを入れていくので
(m-1)**n
となります。ところがこの場合の数には空き箱2個以上の場合も重複しているので重複分を引かなければなりません。そのため、空き箱2個以上の場合の数を計算しますが、この計算には空き箱3個以上の場合が重複しているのでこの重複分を引かなければなりません。という具合に入れ子になっているわけです。ここがややこしいところです。
# ボールn個、箱m個のとき 重複順列と組み合わせを使った数式による計算 箱7個まで
import math
from fractions import Fraction
def c(n, r):
return math.factorial(n) // (math.factorial(n - r) * math.factorial(r))
n = 8 # number of balls
m = 7 # number of boxes
if(m==2):
double_count = 0
elif(m==3):
double_count = int((m-2)**n * c(m-1, 1)/2)
elif(m==4):
double_count = int((m-2)**n*c(m-1, 1)/2
- (m-3)**n)
elif(m==5):
double_count = int((m-2)**n*c(m-1, 1)/2
- ((m-3)**n*c(m-1, 2)/3 - (m-4)**n))
elif(m==6):
double_count = int((m-2)**n*c(m-1, 1)/2
- ((m-3)**n*c(m-1, 2)/3 - ((m-4)**n*c(m-1, 3)/4 - (m-5)**n)))
elif(m==7):
double_count = int((m-2)**n*c(m-1, 1)/2
- ((m-3)**n*c(m-1, 2)/3 - ((m-4)**n*c(m-1, 3)/4 - ((m-5)**n*c(m-1, 4)/5 - (m-6)**n))))
else:
m =2
double_count = 0
print('m out of range. m is set to 2')
total = m**n
not_occupy = ((m-1)**n - double_count)*m
occupy = total - not_occupy
p = Fraction(occupy, total)
print('occupy: {0}, total: {1}, p: {2}, p: {3}'.format(occupy, total, p, round(float(p), 3)))
どの方法が良いかはケースバイケースですが、m、nが大きいときには最初の数値シミュレーションが良いように思います。
