確率変数$X$と$Y$の期待値について、以下が成り立つらしい。(期待値の線形性)
$E(X+Y) = E(X) + E(Y)$
試行の平均を見てみる
確率変数の前に、実際に試行した結果の平均を見てみましょう。
2つのサイコロd1とd2を同時に振った目の合計について考えます。
import random
def gen_dice(values=[1,2,3,4,5,6], prob=None):
def dice():
return random.choices(values, prob)[0]
return dice
import pandas as pd
d1 = gen_dice()
d2 = gen_dice()
results = [[d1(), d2()] for i in range(10000)]
df = pd.DataFrame(results, columns=["d1", "d2"])
df["d1 + d2"] = df.d1 + df.d2
df.describe()
上の実行結果のmean(平均)に着目すると、d1とd2の合計がd1 + d2の値になってますね。
これはなんてことはない、d1の出目の数列($\{d1_n\}$)とd2の出目の数列($\{d2_n\}$)を足し合わせてnで割る(平均を出す)のと、それぞれの数列で平均を出して足すのが同じになるという$\sum$の基本的な計算によります。
$$
\begin{aligned}
&\frac{1}{n}\sum_{k}^{n}(d1_k + d2_k) \\
&= \frac{1}{n}(\sum_{k}^{n}d1_k + \sum_{k}^{n}d2_k ) \\
&= \frac{1}{n}\sum_{k}^{n}d1_k + \frac{1}{n}\sum_{k}^{n}d2_k \\
\end{aligned}
$$
証明をコードで確認する
上の具体例から、確率変数の平均(期待値)についても大体同じような話になりそうな気がしなくもない雰囲気を感じつつ証明を見てみると、以下のようになります。
参照：和の期待値は期待値の和 | 高校数学の美しい物語
$\sum$内の$\Bigg(\Bigg)$は筆者が追加しました。
$$
\begin{align}
&E(X+Y)\\
&=\sum_{x}{\Bigg(\sum_{y}{\Bigg((x+y)P(X=x, Y=y)\Bigg)}\Bigg)} &\quad\\
&=\sum_{x}{\Bigg(\sum_{y}{\Bigg(xP(X=x, Y=y)\Bigg)}\Bigg)} +
\sum_{x}{\Bigg(\sum_{y}{\Bigg(yP(X=x, Y=y)\Bigg)}\Bigg)} &\quad\\
&=\sum_{x}{\Bigg(x\sum_{y}{\Bigg(P(X=x, Y=y)\Bigg)}\Bigg)} +
\sum_{y}{\Bigg(y\sum_{x}{\Bigg(P(X=x, Y=y)\Bigg)}\Bigg)} &\quad\\
&=\sum_{x}{\Bigg(xP(X=x)\Bigg)} + \sum_{y}{\Bigg(yP(Y=y)\Bigg)} &\quad\\
&=E(X) + E(Y) &\quad\\
\end{align}\\
$$
なるほど、分からん。
理解の助けになると思って最初に具体例を出した意味が全くなかったですが、気を取り直して上の証明をコードで考えてみます。
式(2)：素朴な実装
$\sum$が2つ並んでいるの初めて見ましたが、ちょっとプログラミングをかじった人間であれば直感的に「二重ループ」と気づくでしょう。
例によってサイコロを2つ同時に振る場合で考えると、ループを$6 \times 6 = 36$回回しながら、$出目の合計(x+y) \times 確率(P(X=x, Y=y) = \frac{1}{36})$を足し合わせていくという作業を表しています。
素朴に実装すると以下のようになります。
def exp_sum(x_values, x_prob, y_values, y_prob):
exp = 0
for (xv, xp) in zip(x_values, x_prob):
for (yv, yp) in zip(y_values, y_prob):
exp += (xv + yv) * xp * yp
return exp
exp_sum([1,2,3,4,5,6], [1/6,1/6,1/6,1/6,1/6,1/6], [1,2,3,4,5,6], [1/6,1/6,1/6,1/6,1/6,1/6])
式(3)(4)：$\sum$計算
式(3)(4)ではその次のステップに進むための基礎的な変形を行っています。
コードを修正するとこんな感じ。
def exp_sum2(x_values, x_prob, y_values, y_prob):
exp = 0
for (xv, xp) in zip(x_values, x_prob):
prob = 0
for (yv, yp) in zip(y_values, y_prob):
prob += xp * yp
exp += xv * prob
for (yv, yp) in zip(y_values, y_prob):
prob = 0
for (xv, xp) in zip(x_values, x_prob):
prob += xp * yp
exp += yv * prob
return exp
exp_sum2([1,2,3,4,5,6], [1/6,1/6,1/6,1/6,1/6,1/6], [1,2,3,4,5,6], [1/6,1/6,1/6,1/6,1/6,1/6])
式(5)：証明のキモ
唐突に式(5)が出てくるとびっくりするかもしれません。前述の参考ページでは「周辺化」として事前に紹介されていますが、コードでも丁寧に見てみます。
まず内側のループでの処理を以下のように修正します。
def exp_sum3(x_values, x_prob, y_values, y_prob):
exp = 0
for (xv, xp) in zip(x_values, x_prob):
prob = 0
for (yv, yp) in zip(y_values, y_prob):
# prob += xp * yp
prob += yp
# exp += xv * prob
exp += xv * xp * prob
for (yv, yp) in zip(y_values, y_prob):
prob = 0
for (xv, xp) in zip(x_values, x_prob):
# prob += xp * yp
prob += xp
# exp += yv * prob
exp += yv * yp * prob
return exp
exp_sum3([1,2,3,4,5,6], [1/6,1/6,1/6,1/6,1/6,1/6], [1,2,3,4,5,6], [1/6,1/6,1/6,1/6,1/6,1/6])
定数を掛けて足し合わせるのと、足し合わせた値に定数を掛けるのは同じという訳です。
証明の式を追加すれば以下のようになります。
$$
\begin{aligned}
&\sum_{x}{\Bigg(x\sum_{y}{\Bigg(P(X=x, Y=y)\Bigg)}\Bigg)} + \sum_{y}{\Bigg(y\sum_{x}{\Bigg(P(X=x, Y=y)\Bigg)}\Bigg)} &\quad &(4)\\
&= \sum_{x}{\Bigg(xP(X=x)\sum_{y}{\Bigg(P(Y=y)\Bigg)}\Bigg)} + \sum_{y}{\Bigg(yP(Y=y)\sum_{x}{\Bigg(P(X=x)\Bigg)}\Bigg)} &\quad &(4.5)\\
\end{aligned}\\
$$
ここで、内側のループは$X, Y$それぞれの全ての場合の確率を足し合わせたものなので1になります。
するとコードは以下のように修正できます。(式(5)に対応)
def exp_sum4(x_values, x_prob, y_values, y_prob):
exp = 0
for (xv, xp) in zip(x_values, x_prob):
# yの全ての確率を足し合わせると1
# prob = 0
# for (yv, yp) in zip(y_values, y_prob):
# prob += yp
exp += xv * xp
for (yv, yp) in zip(y_values, y_prob):
# xの全ての確率を足し合わせると1
# prob = 0
# for (xv, xp) in zip(x_values, x_prob):
# prob += xp
exp += yv * yp
return exp
exp_sum4([1,2,3,4,5,6], [1/6,1/6,1/6,1/6,1/6,1/6], [1,2,3,4,5,6], [1/6,1/6,1/6,1/6,1/6,1/6])
最終的に残ったコードを見ると、結局$X, Y$の期待値を足しているだけということで証明完了です。
