確率や期待値を求める問題は、数学やプログラミングにおいて中級者（？）向けの典型問題であると言えます。
単発の確率であれば簡単な計算で求めることができますが、試行が複数回にわたる場合、すぐに状況が複雑化します。また、そのような傾向は期待値の問題では更に著しくなります。それらがどのように難しくなるか、それを解析的に（エレガントに）または力任せに（エレファントに）解くにはどうすれば良いか、という点について論じます。
なお、エレガント、エレファントという呼称についてはこの記事を参考にしています。
エレガントな解法、エレファントな解法 〜モンテカルロ法を添えて〜
コイン
半分の確率で表か裏を出すコインがある。このコインが 3 回連続で表になる確率はいくらか。
これはすぐに計算できます。$(\frac{1}{2})^3=0.125$ です。
では、これが次のような問題になった場合どうでしょうか。
半分の確率で表か裏を出すコインがある。このコインを 10 回投げた時に 3 回連続で表になることがある確率はいくらか。
状態数は $2^{10}=1,024$ 通りなので、手計算の全探索は厳しいです。場合分けして考えるのも、実はかなり面倒です。
こういうのはプログラムに頑張ってもらいます。
from itertools import product
result = [0, 1]
okay = 0
total = 0
for res in product(result, repeat=10):
total += 1
for i in range(8):
if res[i:i+3] == (1, 1, 1):
okay += 1
break
print(okay/total) # 0.5078125
という訳で、答えは 0.5078125 でした。
確率が不均等な場合
先の問題は、コインを 1 回投げるごとに枝分かれする状態が等しく 2 つだったから計算しやすかったですが。これが例えば「表が出る確率が $\frac{1}{3}$ のコイン」だったらどうでしょうか。
$\frac{1}{3}$ の確率で表を出すコインがある。このコインを 10 回投げた時に 3 回連続で表になることがある確率はいくらか。
この場合、上のようなプログラムでは対応不可能です。表と裏の結果を等価な状態として扱っても良いのは、確率が 0.5 である時に限られるからです。
実際に乱数を使う（いちばん力任せ）
一番力任せな方法を考えます。
実際に乱数を発生させて大量の試行を行うことです。
最初の問題ではこれを使うまでもありませんでしたが、今はこれに頼ってみましょう。
from random import random
def trial() -> int:
res = [0 for i in range(10)]
for i in range(10):
res[i] = random() < 1/3
for i in range(8):
if all(res[i:i+3]) == True:
return 1
return 0
num_of_trial = 100000
results = [0 for i in range(num_of_trial)]
for i in range(num_of_trial):
results[i] = trial()
print(sum(results)/num_of_trial) # 0.20435 この結果は実行ごとに変動します
とりあえず 0.20435 という結果が得られました。これは乱数による結果なので、この結果が正しいかどうかはとりあえず保留します。
状態の遷移を考える
さて、厳密解を求めるには、プログラムの力を借りるにしても工夫が必要です。
ここで発想を変えて、「状態が変わる」という考え方を持つ必要があります。
つまり、
- 0 回連続した状態
- 1 回連続した状態
- 2 回連続した状態
- 3 回連続した状態
という 4 つの状態を考えて、それらがどのように確率的遷移をするかを考えます。連続する限りは 1 つずつステップを刻むことができますが、1 回でも裏を出した時点で振り出しに戻り「0 回連続した状態」に戻ります。ただし、一回でも 3 連続を達成すれば後はどうしようが自由です。
図で表すとこんな感じになります。
これをプログラムで実現します。
int main()
{
vector<vector<double>> dp(11, vector<double>(4));
double p = double(1)/3;
dp[0][0] = 1; // 0 回目の試行にいく確率は 100 %
rep(i, 10)
{
rep(j, 3)
{
dp[i + 1][j + 1] += dp[i][j] * p; // 表が出たら連続数が増加
dp[i + 1][0] += dp[i][j] * (1 - p); // 裏が出たらふりだしに戻る
}
dp[i + 1][3] += dp[i][3]; // 一回 3 回連続すれば後は安泰
}
cout << dp[10][3] << endl; // 0.202561
}
0.202561 という結果が得られました。これは厳密解となります。
やはり乱数方式の結果は、正解にある程度近い数値が出たものの、$100,000$ 回程度の結果では変数の揺らぎを無視できないほど受けてしまいます。試行回数を増やせばこの差はもちろん縮まりますが、計算時間がそれだけ大きくなります。求める精度にもよりますが、十分な精度を得るまでに現実的な処理時間を超えてしまうケースが多いでしょう。
ちなみに、これらの状態数がどのように遷移するかを詳細に示した結果が以下になります。
// p = 0.5
1 0 0 0
0.5 0.5 0 0
0.5 0.25 0.25 0
0.5 0.25 0.125 0.125
0.4375 0.25 0.125 0.1875
0.40625 0.21875 0.125 0.25
0.375 0.203125 0.109375 0.3125
0.34375 0.1875 0.101562 0.367188
0.316406 0.171875 0.09375 0.417969
0.291016 0.158203 0.0859375 0.464844
0.267578 0.145508 0.0791016 0.507812
// p = 0.3
1 0 0 0
0.666667 0.333333 0 0
0.666667 0.222222 0.111111 0
0.666667 0.222222 0.0740741 0.037037
0.641975 0.222222 0.0740741 0.0617284
0.625514 0.213992 0.0740741 0.0864198
0.609053 0.208505 0.0713306 0.111111
0.592593 0.203018 0.0695016 0.134888
0.576741 0.197531 0.0676726 0.158055
0.561297 0.192247 0.0658436 0.180613
0.546258 0.187099 0.0640824 0.202561
p = 0.5（つまり表と裏が当確率で出る状態）で、上と同じ結論 0.507182 が得られていることがわかります。
つまり、状態で分けて管理することにより確率の計算を解析的に解くことができることがわかります。
期待値（基礎）
さて、いよいよ期待値の問題に移ります。
$p$ の確率で表を出すコインがある。表が 1 回出るのに必要な試行の期待値はいくらか。
ネタバレをすると答えは $\frac{1}{p}$ です。例えば $p=\frac{1}{3}$ だったら、3 回コイントスをしたら 1 回は出るだろうというのは自然な感覚です。しかし、これを厳密に証明するためにはどうすればよいでしょうか。
期待値の定義から考えると、これは 1 回目で出る確率、2 回目で出る確率、...、無限回目で出る確率をすべて求めて、それらと回数との積を足し上げていく必要があります。つまり、無限回の計算が必要で、原理上は永遠に終わりません。期待値のナイーブな計算が難しい根本的な原因はここにあります。
とは言っても、実は無限回足したからといって和が無限になるとは限りません。100 回振って表が 1 回も出ない、などというのは天文学的に少ない確率なので、確率×試行回数の累計は何らかの値に収束するであろうという予測はすぐに立ちます。ここは頑張って解析的に考えてみます。
求める期待値を $E$ とします。
\displaylines{
E = p + 2(1-p)p + 3(1-p)^2p + ... \\
E = p(1 + 2(1-p) + 3(1-p)^2 + ...) \\
E = p\sum_{n=1}^\infty(n(1-p)^{n-1})
}
シグマの中身は等比等差数列という、高校数学だとけっこう難問になるやつですが、大人は微分でズルできることを知っているのでズルします。（この微分テクは特殊な場合のみ適用可能ですが、これは大丈夫です）
\displaylines{
E = p\sum_{n=1}^\infty(n(1-p)^{n-1}) \\
E = p(-\sum_{n=0}^\infty(1-p)^n)'\\
E = p(-(\frac{1}{1-(1-p)}))' \\
E = p(-\frac{1}{p})' \\
E = p(\frac{1}{p^2}) = \frac{1}{p} }
繰り返すと、「成功するまで繰り返す」回数の期待値は、成功確率の逆数になります。言葉にすると自明な感じですけど、解析的に証明できることも重要です。なお、この定理（？）はコンプガチャ問題でも良く引き合いに出されます。
期待値（発展）
次は少し難しくなります。
$\frac{1}{3}$ の確率で表を出すコインがある。表が 3 回連続で出るまでに必要な試行数の期待値を求めよ。
これをナイーブに解析的に解くのはかなり大変です。
とりあえず乱数で力押しします。
from random import random
def trial():
okay = 0
total = 0
while okay < 3:
total += 1
if random() < 1/3:
okay += 1
else:
okay = 0
return total
num_of_trial = 100000
results = [0 for i in range(num_of_trial)]
for i in range(num_of_trial):
results[i] = trial()
print(sum(results)/num_of_trial) # 38.94518
$38.94518$ という結果が得られましたが、当然これは誤差が大きいと予想されます。
次に、上の問題でやったように、$S_0,S_1,S_2,S_3$ に分けて管理します。
ここで、コンプガチャ問題のように「初期状態から遷移に必要な期待値をそれぞれ計算して足し合わせる」ということをやろうとしてもうまくいきません。なぜかというと、$S_0 \rightarrow S_1 \rightarrow S_0 \rightarrow ...$ のようにループする経路が存在するからです。
しかし、各試行において $S_3$ に初めて至る確率はわかるので、それに試行数をかけたものを累積していくことを考えます。
問題の設定により、$S_3$ から出ている矢印が消えている点に注意です。
このような遷移を書きだせば、「初めて3回連続した時の試行数×確率」を足していくことにより、期待値を求めることができます。
int main()
{
int num_of_trial = 1000;
double p = double(1) / 3;
vector<vector<double>> dp(num_of_trial+1, vector<double>(4));
dp[0][0] = 1;
rep(i, num_of_trial)
{
rep(j, 3)
{
dp[i + 1][j + 1] += dp[i][j] * p; // 表が出たら連続数が増加
dp[i + 1][0] += dp[i][j] * (1 - p); // 裏が出たらふりだしに戻る
}
}
cout << dp << endl;
double ans = 0;
rep(i, num_of_trial+1)
{
ans += dp[i][3] * i;
}
cout << ans << endl; // 39
}
$39$ という結果が得られました。ここで気になる点として、これは厳密解なのか？ があります。原理上は無限回足し合わせないと厳密解は出ません。例えば $1,000$ 回目で試行を打ち切った時は、$1,001$ 回目で初めて 3 回連続する可能性を無視しているので、これは厳密解よりも小さい可能性があります。
しかし、結論を言えば実用上問題ない値と考えてよいです。なぜかというと、後半に行けば行くほど確率は累乗スケールで小さくなっていくので、その寄与度は無視できる程小さくなります。問題の設定によりますが、$1,000$ 回程度の試行で十分な精度を得られると思います。
つまり、期待値のような、原理上、無限回の試行が要求される問題であったとしても、状態を分けて管理することにより現実的な時間で力押しで計算することができます。エレガントな考察とエレファントな計算力の組み合わせです。
期待値（厳密ver）
ところが、競プロの問題には、そのような力押しさえ許さない設定が存在します。出力を整数ではなく、mod 分数の形で答えさせる問題です。これは、小数で答えさせると誤差が出るので、分数の形で答えさせるものです。ただし、分数だと出力が煩雑なので、逆元の形で答えさせるものがほとんどですが。
（例：1/3 は mod7 においては 5、なぜならば 3 をかけると 1 になるため）
この形式だと、先のような力押しが通用しません。
ここで、より解析的な考えを進めて試行数の期待値を考える必要があります。
先の状態図と同じように、期待値の試行数の状態図を考えます。$i$ 回連続した状態で、最終状態に進むための期待値を $E_i$ とします。図としては同じですが、遷移は以下のようになります。
- $E_3 = 0$
- $E_2 = (E_3 + 1)\frac{1}{3} + (E_0 + 1)\frac{2}{3}$
- $E_1 = (E_2 + 1)\frac{1}{3} + (E_0 + 1)\frac{2}{3}$
- $E_0 = (E_1 + 1)\frac{1}{3} + (E_0 + 1)\frac{2}{3}$
「$+1$」がついているのは、「一回のステップで次の状態にいった後、その状態からの期待値のステップを要する」というニュアンスです。
これらは 4 変数の連立一次方程式になり、まともに解くと面倒くさいですが、大人は線形代数でズルできることを知っているのでズルします。
\displaylines{
E_3 = 0 \\
-2E_0 + 3E_2 - E_3 = 3 \\
-2E_0 + 3E_1 - E_2 = 3 \\
E_0 - E_1 = 3\\
\begin{bmatrix}
1 & -1 & 0 & 0\\
-2 & 3 & -1 & 0\\
-2 & 0 & 3 & -1\\
0 & 0 & 0 & 1 \\
\end{bmatrix}
\begin{bmatrix}
E_0 \\
E_1 \\
E_2 \\
E_3 \\
\end{bmatrix}
=
\begin{bmatrix}
3 \\
3 \\
3 \\
0 \\
\end{bmatrix}
\\
\begin{bmatrix}
E_0 \\
E_1 \\
E_2 \\
E_3 \\
\end{bmatrix}
=
\begin{bmatrix}
1 & -1 & 0 & 0\\
-2 & 3 & -1 & 0\\
-2 & 0 & 3 & -1\\
0 & 0 & 0 & 1 \\
\end{bmatrix}^{-1}
\begin{bmatrix}
3 \\
3 \\
3 \\
0 \\
\end{bmatrix}
\\
\begin{bmatrix}
E_0 \\
E_1 \\
E_2 \\
E_3 \\
\end{bmatrix}
=
\begin{bmatrix}
9 & 3 & 1 & 1\\
8 & 3 & 1 & 1\\
6 & 2 & 1 & 1\\
0 & 0 & 0 & 1 \\
\end{bmatrix}
\begin{bmatrix}
3 \\
3 \\
3 \\
0 \\
\end{bmatrix} \\
\begin{bmatrix}
E_0 \\
E_1 \\
E_2 \\
E_3 \\
\end{bmatrix}
=
\begin{bmatrix}
39 \\
36 \\
27 \\
0 \\
\end{bmatrix} }
こうして、求める期待値 $E_0 = 39$ が厳密に求まりました。
確率や期待値との向き合い方（まとめ）
- 全探索できるなら全探索
- 「表、裏、表......」等を列挙できる場合
- 状態に分けて遷移を考えることにより高速化
- 「$n$ 回表が連続して出た時」等のようにまとめた状態を考える
- 期待値は、状態遷移を考えた上で累積計算することにより高精度で推定可能
- 期待値の厳密解を求める時は、試行数の期待値の状態遷移を考えることが有効
- 乱数シミュレーションはすべての場合において有効だが、精度は低い
- 検算にはいいかもしれない