1. 下ネタサイコロ
今は販売されていないようですが、ちんかまんかサイコロと呼ばれるものがあるらしいです。飲み会などで使われるサイコロですが、これを応用して空白のサイコロを5個用意して適当に文字を書き、お⚪︎んちんが出やすいようなサイコロを設計したいわけです。
ただし、
- 全部の面が「お」を1つ
- 全部の面が「ち」を2つ
- 全部の面が「ん」を2つ
と準備するのは確率が1になってつまらないので、いくつか制約をつけながら設計しようと思います。
2. 確率の基礎
導入として、「お」、「ち」、「ん」の3文字が2回ずつ出現するサイコロが5個あったときに、これらのサイコロを同時に投げたとき、「お⚪︎んちん」が出る確率を求めてみましょう。一つのサイコロに着目したとき
$$
P\left( \text{お} \right) = \frac{2}{6} \quad P\left( \text{ち} \right) = \frac{2}{6} \quad P\left( \text{ん} \right) = \frac{2}{6}
$$
となります。
ここで、考えたいのは、
- 「お」が出るサイコロを1つ
- 「ち」が出るサイコロを2つ
- 「ん」が出るサイコロを2つ
選べばいいことになります。
この場合の数は、まず「お」が出るサイコロを1つ選び、次に「ち」が出るサイコロを2つ選び、「ん」が出るサイコロを2つ選べばいいことになります。すなわち組み合わせを用いて
$$
_{5}{C} _{1} \times _{4}{C} _{2} \times _{2}{C} _{2}=\frac{5!}{1!4!}\frac{4!}{2!2!} \frac{2!}{2!0!}=\frac{5!}{2!2!}=30
$$
と求めことができます。それぞれの場合の数は上記の通りそれぞれの目が出る確率は$\frac{2}{6}=\frac{1}{3}$ですから、$\left(\frac{1}{3} \right)^{5}$となり、
$$
P\left( \text{おちんちん} \right)= 30 \times \left(\frac{1}{3} \right)^{5} = \frac{10}{81} \simeq 0.12345\cdots
$$
3. 問題設定の作成
ここで、二つの問題を考えてみたいと思います。
- 「お⚪︎んちん」が出やすくなるようにサイコロの目を設計する
- サイコロを振って「おち⚪︎ちん」が出たときに、どのような目の構成だったかを推測する
これらを以下で考えていきます。
4.「お⚪︎んちん」が出やすくなるようにする
4.1. サイコロの目の構成を5個とも同じにする場合
5個とも同じ目とする場合どのように考えればいいでしょうか？なんとなく、「お」を少なめ、2回出現する「ち」と「ん」を多めにすれば良さそうな検討がつきます。すなわち目の数が
$$
\left(\text{お}, \text{ち}, \text{ん} \right)=\left(1,2,3 \right),\left(1,3,2 \right)
$$
とすると確率$P\left( \text{おちんちん} \right)$が最大化しそうな気がします。考察してみます。
今回はサイコロの目が全部同じなので、
$$
P\left( \text{お} \right) = \frac{a}{6} \quad P\left( \text{ち} \right) = \frac{b}{6} \quad P\left( \text{ん} \right) = \frac{c}{6}
$$
とおきます。ただし、
$$
a+b+c=6
$$
となります。2.で考えた設定と同じく
- 「お」が出るサイコロを1つ
- 「ち」が出るサイコロを2つ
- 「ん」が出るサイコロを2つ
と考えれば良いので、サイコロの組み合わせは先ほどと同様に
$$
\frac{5!}{2!2!}=30
$$
ですから
\begin{align}
P\left( \text{おちんちん} \right) &= 30 \times \left(\frac{a}{6} \right)^{1} \left(\frac{b}{6} \right)^{2} \left(\frac{c}{6} \right)^{2} \\
&= 30 \times \frac{1}{6^{5}} \left(6-b-c \right) b^{2}c^{2}
\end{align}
これを最大化することなります。これは定数を無視して
$$
f\left( b,c \right) = b^{2}c^{2}\left(6-b-c \right)
$$
を最大化することになります。ただし、$1 \leq b,c < 5$かつ$b+c \leq 5$の整数という制約があります。さらに$b,c$について対称なので$b \leq c$という制約を設けます。意味的にいうと「⚪︎ちんちん」を構成する「ち」と「ん」の出現回数は等しいので、数学上区別する必要はありません。このくらいであれば、頑張れば紙の上でも議論できますが、コードを書いて考えます。
import numpy as np
values = []
# 制約： b, c は 1〜4, b + c <= 5, b <= c
for b in range(1, 5):
for c in range(1, 5):
if (b <= c) and (b + c <= 5):
# 6面のうち残りが「お」
a = 6 - b - c
# f(b,c) = b^2 c^2 (6 - b - c)
f = (b ** 2) * (c ** 2) * a
values.append((a, b, c, f))
# numpy配列にして f を基準に降順ソート
dtype = [("a", int), ("b", int), ("c", int), ("f", float)]
arr = np.array(values, dtype=dtype)
order = np.argsort(arr["f"])[::-1]
print("=== f(b,c) の値と構成（降順） ===")
for i in order:
a, b, c, f = arr[i]
print(f"(a,b,c)=({a},{b},{c}), f={f}")
この出力は
=== f(b,c) の値と構成（降順） ===
(a,b,c)=(1,2,3), f=36.0
(a,b,c)=(2,2,2), f=32.0
(a,b,c)=(2,1,3), f=18.0
(a,b,c)=(1,1,4), f=16.0
(a,b,c)=(3,1,2), f=12.0
(a,b,c)=(4,1,1), f=4.0
となります。ここで、$b \leq c$の制約を加えたので、この二つは入れ替えてもよいので$\left(a,b,c \right)=\left(1,2,3 \right), \left(1,3,2 \right)$のいずれかになり
\begin{align}
\max{P\left( \text{おちんちん} \right)} &=30 \times \frac{1}{6^{5}} \times 1 \times 2^{2}\times 3^{2} \\
&= \frac{5}{36} \simeq 0.1389
\end{align}
となります。
4.2. 5個のサイコロで異なる目の構成を許す場合
今度は5個のサイコロそれぞれ違う目の構成を許すことにします。一つサイコロに着目したときの確率を
$$
P\left( \text{お} \right) = \frac{a}{6} \quad P\left( \text{ち} \right) = \frac{b}{6} \quad P\left( \text{ん} \right) = \frac{c}{6}
$$
とおきます。ただし$a,b,c >0$という制約を置き、必ず一つのサイコロに各文字1回は出てくるものとします。頑張れば解析的にできなくもないですが、非常に手間なのでいきなりシミュレーションに頼ります。$a,b,c$の組み合わせは
\begin{align}
\left(a,b,c \right) =& \left(1,1,4 \right),\left(1,2,3 \right),\left(1,3,2 \right),\left(1,4,1 \right), \\
& \left(2,1,3 \right),\left(2,2,2 \right),\left(2,3,1 \right), \\
&\left(3,1,2 \right),\left(3,2,1 \right), \\
&\left(4,1,1 \right)
\end{align}
の10通りです。いま、サイコロが5個ですので
$$
10^5
$$
通りですので、コンピューターであれば十分に計算できる計算量です。数も多くないので、簡単のために同じ構成のサイコロが複数あってもよいものとしましょう。
import numpy as np
import itertools
# 1個のサイコロの全候補 (a,b,c)
candidates = [(a,b,6-a-b) for a in range(1,5) for b in range(1,5) if 1 <= 6-a-b <= 4]
best = None
best_p = -1
# 候補のサイコロを5個準備して全パターンを生成してループ
for dice in itertools.product(candidates, repeat=5):
# 各サイコロが出す確率
probs = np.array(dice) / 6 # shape (5,3) → 列順は [お,ち,ん]
# 「おち⚪︎ちん」が出る確率（サイコロそれぞれ独立）
# = sum over 全ての配置 だが組み合わせが固定なので：
# → 各サイコロから「お」1個、「ち」2個、「ん」2個選ぶ確率の総和
# 全サイコロを順に見て、どのサイコロが何を出すか決める
# brute force: 5!/(2!2!) = 30 通りの割当を全探索
p = 0
for assign in set(itertools.permutations(["お","ち","ち","ん","ん"])):
indices = {"お":0, "ち":1, "ん":2}
p += np.prod([probs[i, indices[ch]] for i, ch in enumerate(assign)])
if p > best_p:
best_p = p
best = dice
print("=== 最適構成 ===")
print(best)
print("最大確率 =", best_p)
順を追って解説して行きます。
# 1個のサイコロの全候補 (a,b,c)
candidates = [(a,b,6-a-b) for a in range(1,5) for b in range(1,5) if 1 <= 6-a-b <= 4]
best = None
best_p = -1
# 候補のサイコロを5個準備して全パターンを生成してループ
for dice in itertools.product(candidates, repeat=5):
この部分で1つのサイコロの全候補をまず作り、それを5つ分全組み合わせを作ってfor文を回します。repeat=5で5個分増幅しています。
次に
probs = np.array(dice) / 6
で、$\left[\frac{a}{6},\frac{b}{6},\frac{c}{6} \right]$という配列を作っています。一つのサイコロでそれぞれの目が出る確率となります。これを5行分作っています。各行がそれぞれのサイコロとなります。
for assign in set(itertools.permutations(["お","ち","ち","ん","ん"])):
ここでは5個のサイコロがそれぞれ「お」、「ち」、「ん」のどの組み合わせを担当するかの全パターンを列挙します。
p += np.prod([probs[i, indices[ch]] for i, ch in enumerate(assign)])
probsの行がサイコロ、列が「お」、「ち」、「ん」の順でそのサイコロがその文字を出現する確率になりますので、サイコロを5個準備したときにおちん⚪︎んが確率をそれぞれのパターンで求めてその確率を加えています。
これをすべての「おち⚪︎ちん」の出現パターンで加えて、もし今まで出てきた最大値よりおち⚪︎ちんの出現確率が大きくなる場合はそのパターンを上書きして保持します。
if p > best_p:
best_p = p
best = dice
これの出力結果は以下の通りです。
((1, 1, 4), (1, 4, 1), (1, 1, 4), (1, 4, 1), (4, 1, 1))
最大確率 = 0.21913580246913583
すなわち、
- 「お」担当のサイコロを1つ(お:4、ち:1, ん:1)
- 「ち」担当のサイコロを2つ(お:1、ち:4, ん:1)
- 「ん」担当のサイコロを2つ(お:1、ち:1, ん:4)
と書くとおち⚪︎ちんが出る確率は最大となり、その最大値は
\begin{align}
\max{P\left( \text{おちんちん} \right)} \simeq 0.2191
\end{align}
となります。
5 おち⚪︎ちんから目の構成を推測する
今度は観測結果から、目の構成を推測したいと思います。実際には観測結果そのものよりも構成を推測したい場合が多いのでその手法を紹介します。4.2のケースは少し難しいので、4.1のサイコロが全て同じ目の構成だったと仮定します。そのとき「おち⚪︎ちん」の出現状況により、目の構成を推測することを考えます。このように得られた結果からパラメータ（今回は$a,b,c$。すなわちサイコロの目の構成）を推測する手法を「最尤法」と呼びます。
5.1. サイコロ一回振ってお⚪︎んちんが出た場合
まずサイコロを一回だけ振って「おち⚪︎ちん」が出たとします。
4.1.の議論からそれぞれの目が出る確率は
$$
P\left( \text{お} \right) = \frac{a}{6} \quad P\left( \text{ち} \right) = \frac{b}{6} \quad P\left( \text{ん} \right) = \frac{c}{6}
$$
となります。このときおちんちんが出る確率は
\begin{align}
P\left( \text{おちんちん} \right) &= 30 \times \frac{1}{6^{5}} \left(6-b-c \right) b^{2}c^{2}
\end{align}
となります。このとき最尤法の発想は「おち⚪︎ちん」が最も出やすくなるようにパラメータ($a,b,c$)を選びましょうという発想です。このとき、次のように考えます。
$$
L\left(a, b,c \right) = 30 \times \frac{1}{6^{5}} a b^{2}c^{2}
$$
を考え、この関数を最大化するように$a,b,c$を決めます。この最大化する関数のことを$尤度関数$と呼びます。ここでポイントになっているのは$a,b,c$が変数になっているということです。確率を求めるときはこれらのパラメータはすでに決まっている値となりますが、最尤法ではこれらが変数となっている点です。
このとき
$$
L\left(a, b,c \right) \propto b^{2}c^{2} \left(6-b-c \right)
$$
ですから、先ほど考えた
$$
f\left( b,c \right) = b^{2}c^{2}\left(6-b-c \right)
$$
と同じ関数となり、$\left(a,b,c \right)=\left(1,2,3 \right), \left(1,3,2 \right)$のいずれかになります。
5.2. サイコロ複数回振って複数回お⚪︎んちんが出た場合
ここでサイコロを複数回振っておちんちんが何度か観測されたとします。実際にはおち⚪︎ちん以外にも「おおちちん」などが出る場合もあります。$k$回観測したとして、その結果が最も起こりやすくなるようにサイコロの目の構成を選ぶのが最尤法の考え方です。$i$回目にとある文字が出たとします。これは「お」、「ち」、「ん」の三文字が0個以上合計5個出ることになります。そこで、それぞれの個数を$n_{i\text{お}},n_{i\text{ち}},n_{i\text{ん}}$とおきます。これは
$$
n_{i\text{お}},n_{i\text{ち}},n_{i\text{ん}} \geq 0, \quad n_{i\text{お}}+n_{i\text{ち}}+n_{i\text{ん}} = 5
$$
を満たします。これまでと同じように一つのサイコロの目の構成は
$$
P\left( \text{お} \right) = \frac{a}{6} \quad P\left( \text{ち} \right) = \frac{b}{6} \quad P\left( \text{ん} \right) = \frac{c}{6}
$$
とおきます。このとき、とある文字の出現確率は
\begin{align}
P\left( \text{*} \right) &= _{5}{C} _{n_{i\text{お}}} \times _{5-n_{i\text{お}}}{C} _{n_{i\text{ち}}} \times _{5-n_{i\text{お}}-n_{i\text{ち}}}{C} _{n_{i\text{ん}}} \left( \frac{a}{6} \right)^{n_{i\text{お}}} \left( \frac{b}{6} \right)^{n_{i\text{ち}}} \left( \frac{c}{6} \right)^{n_{i\text{ん}}} \\
&= \frac{5!}{n_{i\text{お}}!\left(5-n_{i\text{お}} \right)!} \frac{\left(5-n_{i\text{お}} \right)!}{n_{i\text{ち}}!\left(5-n_{i\text{お}}-n_{i\text{ち}} \right)!}\frac{\left(5-n_{i\text{お}}-n_{i\text{ち}} \right)!}{n_{i\text{ん}}!0!} \left( \frac{a}{6} \right)^{n_{i\text{お}}} \left( \frac{b}{6} \right)^{n_{i\text{ち}}} \left( \frac{c}{6} \right)^{n_{i\text{ん}}} \\
&= \frac{5!}{n_{i\text{お}}!n_{i\text{ち}}!n_{i\text{ん}}!} \left( \frac{a}{6} \right)^{n_{i\text{お}}} \left( \frac{b}{6} \right)^{n_{i\text{ち}}} \left( \frac{c}{6} \right)^{n_{i\text{ん}}}
\end{align}
となります。ここで
\bf{n}_{i}=
\left(
\begin{array}{c}
n _{i\text{お}} \\
n _{i\text{ち}} \\
n _{i\text{ん}}
\end{array}
\right)
とおくと、このとき$k$回試行したとしたときの尤度は
\begin{align}
L\left(a,b,c \right) &= P\left( \bf{n}_{1} \right)P\left( \bf{n}_{2} \right) \cdots P\left( \bf{n}_{k} \right) \\
&= \prod_{i=1}^{k}{P\left( \bf{n}_{i} \right)}
\end{align}
一般にこの尤度は確率の掛け算のため1より小さく、それの大量の掛け算となるため、アンダーフローが起こる場合などがあるため、尤度に対数を取った対数尤度$l\left( * \right)=\log{L\left(* \right)}$が用いられます。ここで対数尤度
\begin{align}
l\left(a,b,c \right) &= \log{\prod_{i=1}^{k}{P\left( \bf{n}_{i} \right)}} \\
&= \sum_{i=1}^{k}\log{P\left( \bf{n}_{i} \right)} \\
&\propto \sum_{i=1}^{k} \log{\left( \frac{a}{6} \right)^{n_{i\text{お}}} \left( \frac{b}{6} \right)^{n_{i\text{ち}}} \left( \frac{c}{6} \right)^{n_{i\text{ん}}}} \\
&= \log{\left( \frac{a}{6} \right)^{\sum_{i=1}^{k}n_{i\text{お}}} \left( \frac{b}{6} \right)^{\sum_{i=1}^{k}n_{i\text{ち}}} \left( \frac{c}{6} \right)^{\sum_{i=1}^{k}n_{i\text{ん}}}}
\end{align}
これを最大化するように求めます。
\begin{align}
l^{*} \left(a,b,c \right) = \log{\left( \frac{6-b-c}{6} \right)^{\sum_{i=1}^{k}n_{i\text{お}}} \left( \frac{b}{6} \right)^{\sum_{i=1}^{k}n_{i\text{ち}}} \left( \frac{c}{6} \right)^{\sum_{i=1}^{k}n_{i\text{ん}}}}
\end{align}
とおきます。ここで$b$で微分すると
\begin{align}
\frac{\partial l^{*} }{\partial b} &= \frac{\partial}{\partial b}{\left\{ \sum_{i=1}^{k}n_{i\text{お}}\log{\left( 6-b -c \right)} + \sum_{i=1}^{k}n_{i\text{ち}}\log{b}\right\}} \\
&= -\sum_{i=1}^{k}n_{i\text{お}} \frac{1}{6-b-c} + \sum_{i=1}^{k}n_{i\text{ち}} \frac{1}{b}
\end{align}
同様に$c$で微分すると
\begin{align}
\frac{\partial l^{*} }{\partial c} &= -\sum_{i=1}^{k}n_{i\text{お}} \frac{1}{6-b-c} + \sum_{i=1}^{k}n_{i\text{ん}} \frac{1}{c}
\end{align}
となる。
\frac{\partial l^{*} }{\partial b} = \frac{\partial l^{*} }{\partial c} = 0
とおくと、
\frac{b}{c}= \frac{\sum_{i=1}^{k}n_{i\text{ち}}}{\sum_{i=1}^{k}n_{i\text{ん}}}
制約の使い方を変えて対数尤度を微分して同様の操作をすると
\frac{a}{b}= \frac{\sum_{i=1}^{k}n_{i\text{お}}}{\sum_{i=1}^{k}n_{i\text{ち}}}
となることから、
$$
a: b : c = \sum_{i=1}^{k}n_{i\text{お}} : \sum_{i=1}^{k}n_{i\text{ち}} : \sum_{i=1}^{k}n_{i\text{ん}}
$$
これにより、目の数は全体の実験をしたときの、トータルの「お」、「ち」「ん」の比率であることがわかる。
最後にこれを実験してみましょう。
import numpy as np
from collections import Counter
# --- おちんちん判定関数 ----------------------------------
def is_ochinchin(roll):
"""
roll: 配列例 ["お","ち","ん","ち","ん"]
並び替えて「おちんちん」(「お」1・「ち」2・「ん」2) を作れるなら True
"""
cnt = Counter(roll)
return cnt["お"] == 1 and cnt["ち"] == 2 and cnt["ん"] == 2
# --- サイコロ設定 ----------------------------------------
true_a, true_b, true_c = 3, 2, 1 # (お, ち, ん) の出現数
def roll_once(a, b, c, n_dice=5, rng=None):
"""1回の試行で5個のサイコロを振る（同じ構成のサイコロを5個とする）"""
if rng is None:
rng = np.random.default_rng()
faces = np.array(["お"] * a + ["ち"] * b + ["ん"] * c)
return rng.choice(faces, size=n_dice)
def simulate_rolls(a, b, c, n_dice=5, n_trials=1000, seed=0):
"""n_trials 回振って文字カウントと、おちんちん出現数を返す"""
rng = np.random.default_rng(seed)
counts = Counter()
ochinchin_count = 0
for _ in range(n_trials):
r = roll_once(a, b, c, n_dice=n_dice, rng=rng)
counts["お"] += np.sum(r == "お")
counts["ち"] += np.sum(r == "ち")
counts["ん"] += np.sum(r == "ん")
# ここでおちんちん判定
if is_ochinchin(r):
ochinchin_count += 1
return counts, ochinchin_count
def mle_face_counts(counts):
"""観測データから (a,b,c) の比率と a+b+c=6 での推定値を返す"""
ratios = np.array([counts["お"], counts["ち"], counts["ん"]], dtype=float)
abc_hat = 6 * ratios / ratios.sum()
return ratios, abc_hat
# ==== 実験 ====
counts, ochinchin_count = simulate_rolls(true_a, true_b, true_c, n_dice=5, n_trials=10000, seed=42)
ratios, abc_hat = mle_face_counts(counts)
print("[観測された合計個数]")
print(counts)
print("\n[最尤推定された比 (お:ち:ん)]")
print(ratios)
print("\n[a+b+c=6 制約下での推定値]")
print(abc_hat)
print("\n[真のサイコロ構成]")
print((true_a, true_b, true_c))
print("\n[おち⚪︎ちん出現回数]")
print(ochinchin_count)
print("頻度 ≈", ochinchin_count / 10000)
順を追って解説します。
def roll_once(a, b, c, n_dice=5, rng=None):
"""(a,b,c) のサイコロを n_dice 個ふったときの出目（文字列配列）を1回分返す。"""
if rng is None:
rng = np.random.default_rng()
faces = np.array(["お"] * a + ["ち"] * b + ["ん"] * c)
roll = rng.choice(faces, size=n_dice)
return roll
これはサイコロの各目の数を受け取って、一回分試行(5個振る)を一回やる関数です。
def simulate_rolls(a, b, c, n_dice=5, n_trials=1000, seed=0):
"""n_trials 回振って文字カウントと、おちんちん出現数を返す"""
rng = np.random.default_rng(seed)
counts = Counter()
ochinchin_count = 0
for _ in range(n_trials):
r = roll_once(a, b, c, n_dice=n_dice, rng=rng)
counts["お"] += np.sum(r == "お")
counts["ち"] += np.sum(r == "ち")
counts["ん"] += np.sum(r == "ん")
# ここでおちんちん判定
if is_ochinchin(r):
ochinchin_count += 1
return counts, ochinchin_count
実際のシミュレーションをする関数です。サイコロを10000回振って「お」、「ち」、「ん」の出現回数をカウントしています。
def mle_face_counts(counts):
"""観測データから (a,b,c) の比率と a+b+c=6 での推定値を返す"""
ratios = np.array([counts["お"], counts["ち"], counts["ん"]], dtype=float)
abc_hat = 6 * ratios / ratios.sum()
return ratios, abc_hat
ここで、「お」、「ち」、「ん」の出現回数から目の数を推測しています。
この例ではシミュレーションですが、実際には、シミュレーションの部分が実験結果をファイルなどで読み込ませて、そこから同様の計算をして目の数を推測することになります。今回は乱数固定をしていないので、結果は環境依存になりますが、おおよそ次のような出力結果が得られるはずです。
観測された合計個数]
Counter({'お': np.int64(24971), 'ち': np.int64(16730), 'ん': np.int64(8299)})
[最尤推定された比 (お:ち:ん)]
[24971. 16730. 8299.]
[a+b+c=6 制約下での推定値]
[2.99652 2.0076 0.99588]
[真のサイコロ構成]
(3, 2, 1)
[おち⚪︎ちん出現回数]
485
頻度 ≈ 0.0485
6.まとめ
今回はサイコロの目の数から最尤推定について概観しました。実務上は解析的に解けないような関数のパラメータを予想することも多く、その場合は勾配法などを用いて最小値を求めていきます。考え方は同じなので、この記事が少しでも役に立てますと幸いです。