月: 2021年7月

前回の記事で名前だけ登場した指数分布についてついでに整理しておきます。
参考: 幾何分布の無記憶性について

指数分布は幾何分布の連続分布版のような確率分布です。
(古さに関係なく一定確率で壊れる機械について)機械が故障するまでの時間や、
(単位時間あたり一定確率で発生する災害について)災害が発生するまでの時間など、
一定確率で発生する何かしらの事象が、次に発生するまでの時間が従う分布です。

数学的には次のように定義されます。
パラメーター$\lambda > 0$に対して、次の確率密度関数を持つ分布を指数分布と呼び、$Exp(\lambda)$と書きます。
$$
f(x;\lambda) = \left\{
\begin{align}
&\lambda e^{-\lambda x} \quad & (x \geq 0)\\
&0 \quad & (x < 0)
\end{align}
\right.
$$
期待値は$\frac{1}{\lambda}$、分散は$\frac{1}{\lambda^2}$です。

モーメント母関数を使うと簡単に導出できますので見ておきましょう。
まず、モーメント母関数は$t<\lambda$の範囲で次のように定義されます。
$$
\begin{align}
M_X(t) &= E(e^{tX})\\
&= \int_{0}^{\infty} e^{tx}\lambda e^{-\lambda x} dx\\
&= \lambda \int_{0}^{\infty} e^{(t-\lambda)x} dx\\
&= \frac{\lambda}{\lambda -t}.
\end{align}
$$

これの微分は簡単ですね。1回微分と2回微分はそれぞれ次のようになります。
$$
\begin{align}
\frac{d}{dt}M_X(t) &= \frac{\lambda}{(\lambda-t)^2}\\
\frac{d^2}{dt^2}M_X(t) &= \frac{2\lambda}{(\lambda-t)^3}.
\end{align}
$$

これを使うと期待値と分散は次のように計算できます。
$$
\begin{align}
E(X) &= \left.\frac{d}{dt}M_X(t)\right|_{t=0}\\
&= \frac{\lambda}{(\lambda-0)^2}\\
&= \frac{1}{\lambda}.
\end{align}
$$
$$
\begin{align}
E(X^2) &= \left.\frac{d^2}{dt^2}M_X(t)\right|_{t=0}\\
&= \frac{2\lambda}{(\lambda-0)^3}\\
&= \frac{2}{\lambda^2}
\end{align}
$$
より、
$$
\begin{align}
V(X) &= E(X^2) – E(X)^2\\
&= \frac{2}{(\lambda-0)^2} – \left(\frac{1}{\lambda}\right)^2\\
&= \frac{1}{\lambda^2}.
\end{align}
$$

前回の記事でも触れました通り、指数分布は無記憶性を持つ連続分布です。
$x_1, x_2 \geq 0$に対して、$P(X\geq x_1+x_2|X\geq x_1) = P(X\geq x_2)$が成り立ちます。
実際、$x>0$とすると、 $P(X\geq x)=e^{-\lambda x}$ ですから、
$$
\begin{align}
P(X\geq x_1+x_2|X\geq x_1) &= \frac{P(X\geq x_1+x_2)}{P(X\geq x_1)}\\
&= \frac{e^{-\lambda (x_1+x_2)}}{e^{-\lambda x_1}}\\
&= e^{-\lambda x_2}\\
&= P(X\geq x_2)
\end{align}
$$
となります。

ここ数回の記事で幾何分布に関連する話を取り上げているので、ついでに幾何分布が持つ無記憶性という性質について紹介します。
これは条件付き確率を用いて、次の数式で表される性質です。

$$
P(X > m+n|X > m) = P(X > ​n) \quad \text{ただし}m, n\geq 0.
$$

まず、幾何分布について上の数式が成り立つとを確認しておきましょう。
$P(X=k) = p(1-p)^{k-1}$ ですから、
$$
\begin{align}
P(X>n) &= \sum_{k=n+1}^{\infty}p(1-p)^{k-1}\\
&= p\cdot\frac{(1-p)^{n}}{1-(1-p)}\\
&= (1-p)^{n}
\end{align}
$$
となります。

よって、
$$
\begin{align}
P(X> m+n|X > m) &= \frac{P(X> m+n \land X > m)}{P(X > m)}\\
&= \frac{P(X > m+n)}{P(X > m)}\\
&= \frac{(1-p)^{m+n}}{(1-p)^{m}}\\
&= (1-p)^{n}\\
&= P(X>n)
\end{align}
$$
となり、幾何分布が冒頭の数式を満たすことが示されました。

これはどういうことか説明します。
幾何分布は確率$p$で成功する独立な試行を、初めて成功するまで繰り返すときに要した回数の分布ですから、
$P(X>n)$というのは、初めて成功するまでに$n+1$回以上かかる確率、言い換えると初めて成功するまでに$n$回以上失敗する確率になります。
これに対して、$P(X > m+n|X > m)$はどういうことかというと、成功するまでに$m+1$回以上かかる、つまりすでに$m$回失敗したという条件のもとで、
成功するのに$m+n+1$回以上かかる、つまり追加で$n$回以上失敗し成功するまでに$n+1$回以上かかる確率を意味します。

この二つが等しいということはどういうことかというと、
成功するまでに$n$回以上失敗する確率は、今の時点で何回失敗しているかという事実に全く影響を受けないということです。

例えば、$1/20$の確率で当たりが出るクジで、連続して20回ハズレを引くと、
そろそろ当たるんじゃないかなという気がしてくる人も多いと思うのですが、
そんなことは全くなく、この先あたりを引くまでにかかる回数の期待値は全く変わってないということを示しています。

この無記憶性は、離散分布の中では幾何分布だけが持つ性質です。
(逆にいうと、離散分布で、無記憶性を持っていたらそれは幾何分布だと言えます。)
このほか、連続分布まで範囲を広げると、指数分布が幾何分布同様に無記憶性を持ちます。

前回の記事で、コンプガチャの期待値と分散を求めましたが、いまいち自信がなかったのでシミュレーションしてみました。
参考: 全種類の景品を集めるのに必要な回数の期待値

おさらいしておくと、$n$種類の景品があるクジを景品が全種類揃うまで引く回数は、
期待値が$n\sum_{k=1}^{n}\frac{1}{k}$, 分散が$n\sum_{k=1}^{n-1}\frac{k}{(n-k)^2}$です。

実際にそのような結果になるのか、仮に$n=20$として、プログラムで繰り返し実行してみましょう。
ちなみに、$n=20$の場合の期待値と分散は次ようになります。

import numpy as np


n = 20
# 期待値
print(sum(n/np.arange(1, n+1)))
# 71.95479314287363

# 分散
print(sum([n*k/(n-k)**2 for k in range(1, n)]))
# 566.5105044223357

シミュレーションに使うために、景品が全種類揃うまでクジを引く関数を実装します。

def complete_gacha(n):
    # 揃ったアイテムの配列
    item_list = []

    # n種類揃うまでクジを引く
    while len(set(item_list)) < n:
        item_list.append(np.random.randint(n))

    return item_list

ためしに$n=5$で実行すると次のような結果が得られます。

print(complete_gacha(5))
# [4, 2, 0, 2, 4, 1, 3]

それでは、この関数を100000回実行し、回数(=返された配列の長さ)のリストを作って、平均値と不変分散を出してみましよう。

result_list = np.array([len(complete_gacha(20)) for _ in range(100000)])

# 期待値
print(result_list.mean())
# 71.86858

# 不偏分散
print(result_list.var(ddof=1))
# 566.9716985005849

試行回数がかなり大きいのもあって、理論値にかなり近い結果が得られましたね。
どうやら前回の記事の結果は一応正しそうです。

いわゆるコンプガチャの問題です。

$n$種類の景品があるクジにおいて、全ての景品を揃えるためには何回程度引けば良いのか(=期待値)を考えていたところ、
うまく解けたので記事にすることにしました。

まず改めて問題の前提条件を整理しておきます。
– クジには$n$種類の景品がある。
– どの景品も当たる確率は等しく$1/n$である。
– クジは無限にあり、過去の景品は将来のクジの確率に影響しない。
– 全集類の景品を最低1回引くまでクジを引き続け、その回数の期待値を求める。

3番目の条件は重要です。要するに景品Aがなかなか出なかったからといって、そのあと景品Aを引ける確率が上がったりしないということです。

当初、場合分けを色々考えてアプローチしていたのですが、次のように考えるとすんなり解けました

この問題は、全ての景品が揃うまでのクジの回数を次のように分解して考えます。

全ての景品が揃うまでの回数 =
1種類目の景品を引くまでの回数
+ 1種類持っている状態から2種類目の景品を引くまでの回数
+ 2種類持っている状態から3種類目の景品を引くまでの回数
・・・
+ $n-1$種類持っている状態から$n$種類目の景品を引くまでの回数

すると、期待値の線形性から次のようになります。
$$
E(\text{全ての景品が揃うまでの回数}) = \sum_{k=1}^{n}E(\text{k-1種類持っている状態からk種類目の景品を引くまでの回数})
$$

あとは、$k-1$種類持っている状態から$k$種類目の景品を引くまでの回数の期待値がわかれば良いです。
ここで、$k-1$種類持っているということは、持っていない景品は$n-k+1$種類であり、
これは、$\frac{n-k+1}{n}$の確率で当たり(=まだ持ってない景品)を引けるクジと考えることができます。
そして、そのあたりを引くまでの回数は、$p=\frac{n-k+1}{n}$の幾何分布に従うので、
前回の記事で見た通り、その期待値は$\frac{1}{p}=\frac{n}{n-k+1}$となります。

よって、
$$
E(\text{k-1種類持っている状態からk種類目の景品を引くまでの回数}) = \frac{n}{n-k+1}
$$
ですから、
$$
E(\text{全ての景品が揃うまでの回数}) = \sum_{k=1}^{n}\frac{n}{n-k+1} = \frac{n}{n} + \frac{n}{n-1} + \cdots + \frac{n}{1}
$$
となり、和の順番を入れ替えて$n$で括ると、
$$
\begin{align}
E(\text{全ての景品が揃うまでの回数}) &= n\left(1+\frac{1}{2}+\frac{1}{3}+\cdots \frac{1}{n}\right)\\
&= n\sum_{k=1}^{n}\frac{1}{k}
\end{align}
$$
となります。
シンプルで美しい結果になりましたね。

ついでにですが、分散も求めておきましょう。
$X$と$Y$が独立の時、期待値同様に分散も$V(X+Y)=V(X)+V(Y)$と分解できることに注意すると以下のようになります。
$$
V(\text{全ての景品が揃うまでの回数}) = \sum_{k=1}^{n}V(\text{k-1種類持っている状態からk種類目の景品を引くまでの回数}).
$$
パラメーターが$p$の幾何分布の分散は、$\frac{1-p}{p^2}$ですから、$p=\frac{n-k+1}{n}$を代入すると、
$$
\begin{align}
V(\text{k-1種類持っている状態からk種類目の景品を引くまでの回数}) &= \frac{1-\frac{n-k+1}{n}}{\left(\frac{n-k+1}{n}\right)^2}\\
&=\frac{n(k-1)}{(n-k+1)^2}
\end{align}
$$
となります。
よって、求めたい分散は、
$$
\begin{align}
V(\text{全ての景品が揃うまでの回数}) &= \sum_{k=1}^{n}\frac{n(k-1)}{(n-k+1)^2}\\
&= n\sum_{k=1}^{n-1}\frac{k}{(n-k)^2}
\end{align}
$$
となります。

この次の記事で、幾何分布の性質(期待値)を使いたいのでおさらいしておきます。

おさらい:
成功確率が$p$である独立なベルヌーイ試行を繰り返す時、初めて成功するまでの試行回数$X$が従う確率分布を幾何分布と言います。
(本やサイトによっては、初めて成功する回数ではなく、初めて成功するまでの失敗回数で定義することもあります。
負の二項分布との関係が明確になったり、台が0始まりになったりするので実は個人的にはそちらの方が好みです。)

確率質量関数は次のようになります。
$$
P(X=k) = (1-p)^{k-1}p \quad (k = 1, 2, 3, \cdots).
$$
確率$1-p$で発生する失敗を$k-1$回続けた後に、確率$p$で発生する成功を1回、と考えれば自明ですね。

期待値は$\frac{1}{p}$、分散は$\frac{(1-p)}{p^2}$です。

モーメント母関数は次の式になります。
$$
M_X(t) = \frac{pe^t}{1-(1-p)e^t} \qquad(t<-\log{(1-p)}). $$ 期待値と分散はこのモーメント母関数から算出することもできるのですが、見ての通り、この関数の微分、2回微分を計算していくのは結構手間です。 なので、幾何分布に関しては、期待値と分散は直接計算するのがおすすめです。 (といってもこれもそこそこトリッキーなことをするのですが。) では、期待値から導出していきましょう。まず期待値の定義です。 $$ E(X) = \sum_{k=1}^{\infty}k(1-p)^{k-1}p. $$ ここで、$\frac{1}{1-y}$のマクローリン展開を考えます。 $$ \frac{1}{1-y} = 1+y+y^2+\cdots = \sum_{k=0}^{\infty}y^k. $$ これを両辺$y$で微分すると次の式になります。(しれっと微分と極限の順序交換をしています。数学科の学生さんなどはこの辺りも厳密に議論することをお勧めします。) $$ \frac{1}{(1-y)^2} = \sum_{k=1}^{\infty}ky^{k-1}. $$ この式に、$y=1-p$を代入すると次のようになります。 $$ \sum_{k=1}^{\infty}k(1-p)^{k-1} = \frac{1}{(1-1+p)^2} = \frac{1}{p^2}. $$ 両辺に$p$をかけることで、 $$ E(X) = \frac{1}{p} $$ が導かれました。 成功率が$p=\frac{1}{n}$のベルヌーイ試行は、平均$\frac{1}{p}=n$回で成功する、と考えると直感とよくあいますね。 続いて、分散$V(X)$を導出するために、$E(X^2)$を計算していきましょう。 先ほどのマクローリン展開の微分の式から始めます。 $$ \frac{1}{(1-y)^2} = \sum_{k=1}^{\infty}ky^{k-1}. $$ この両辺に、$y$を掛けます。 $$ \sum_{k=1}^{\infty}ky^{k} = \frac{y}{(1-y)^2} $$ そして、この両辺をもう一回$y$で微分します。 $$ \sum_{k=1}^{\infty}k^2y^{k-1} = \frac{1+y}{(1-y)^3}. $$ $y=1-p$とすると、 $$ \sum_{k=1}^{\infty}k^2(1-p)^{k-1} = \frac{2-p}{p^3}. $$ 両辺に$p$をかけて、 $$ \sum_{k=1}^{\infty}k^2(1-p)^{k-1}p = \frac{2-p}{p^2}. $$ この左辺は$E(X^2)$ですね。 よって、分散$V(X)$は、次のように求まります。 $$ \begin{align} V(X) &=E(X^2) - E(X)^2\\ &=\frac{2-p}{p^2} - \left(\frac{1}{p}\right)^2\\ &=\frac{1-p}{p^2}. \end{align} $$

前回の記事で二項分布の期待値と分散を直接計算したわけですが、記事中でも述べている通り、
二項分布の期待値や分散を導出するのはモーメント母関数を使った方が楽です。
マイナーな方法だけ紹介しているというのも変なので、この記事で二項分布のモーメント母関数について紹介します。

早速ですが、確率関数は
$$P[X=k] = {}_{n}\mathrm{C}_{k}p^k(1-p)^{n-k}$$
なので、モーメント母関数は次のようになります。
$$
\begin{align}
M_X(t) &= E(e^{tX})\\
&= \sum_{k=0}^{n} e^{tk}{}_{n}\mathrm{C}_{k}p^k(1-p)^{n-k}\\
&= \sum_{k=0}^{n} {}_{n}\mathrm{C}_{k} (pe^{t})^k(1-p)^{n-k}\\
&= (pe^t + 1 -p)^n.
\end{align}
$$
最後の行の式変形は二項定理を使いました。

モーメント母関数を$t$で1回微分すると次式になります。
$$
\frac{d}{dt} M_X(t) = npe^t(pe^t + 1 -p)^{n-1}.
$$
$t=0$を代入することで期待値が得られます。
$$
\begin{align}
E(X) &= \left.\frac{d}{dt} M_X(t)\right|_{t=0}\\
&=np(p+1-p)\\
&=np.
\end{align}
$$
前回の記事の直接計算するのに比べて若干楽なのが感じられると思います。

続いて分散です。モーメント母関数の2回微分は次のようになります。
$$
\frac{d^2}{dt^2} M_X(t) = npe^t(pe^t + 1 -p)^{n-1} + n(n-1)p^2e^{2t}(pe^t + 1 -p)^{n-2}.
$$
正確には$n=1$の場合と$n\geq2$の場合でそれぞれ計算しないといけないのですが、結局どちらの場合も上の式で表されることが証明できます。

ここから二項分布の2次のモーメントが次のように計算できます。
$$
\begin{align}
E(X^2) &= \left.\frac{d^2}{dt^2} M_X(t)\right|_{t=0}\\
&= np(p+1-p)^{n-1} + n(n-1)p^2(p+1-p)^{n-2}\\
&= np +n^2p^2-np^2
\end{align}
$$

よって、二項分布の分散は次のように導出されます。
$$
\begin{align}
V(X) &= E(X^2)-E(X)^2\\
&= np +n^2p^2-np^2 – (np)^2\\
&= np(1-p).
\end{align}
$$

分散に関しては、直接計算するに比べてモーメント母関数を使った方がはるかに楽に導出できましたね。

おさらい:
成功確率が$p$のベルヌーイ試行を独立に$n$回行った時の成功回数を確率変数とする分布を二項分布と呼び、$B(n, p)$と書きます。
確率関数は次の式になります。
$$
P[X=k] = {}_{n}\mathrm{C}_{k}p^k(1-p)^{n-k}.
$$

期待値$E(X)$と分散$V(X)$は次の式で表されることが知られています。
$$
\begin{align}
E(X) &= np.\\
V(X) &= np(1-p).
\end{align}
$$

色々なテキストを見ると期待値や分散の導出はモーメント母関数を使われているのをよく見かけます。
最近復習と計算の練習を兼ねて、これらをモーメント母関数を使わずに定義から直接算出してみたところ、思ったより手こずったので記事に残すことにしました。

式変形の途中で二項係数の次の性質を使いますので注意してみてください。
$x\geq1$の時、$x\cdot{}_{n}\mathrm{C}_{n} = n\cdot {}_{n-1}\mathrm{C}_{x-1}$です。

証明 $x\geq1$とすると、
$$
\begin{align}
x\cdot{}_{n}\mathrm{C}_{x} &= x\frac{n!}{x!(n-x)!}\\
& = n\frac{(n-1)!}{(x-1)!((n-1) – (x-1))!}\\
& = n \cdot {}_{n-1}\mathrm{C}_{x-1}.
\end{align}
$$

それでは、本題に戻って期待値$E(X)$から算出していきます。
$$
\begin{align}
E(X) &= \sum_{x=0}^{n}x\cdot{}_{n}\mathrm{C}_{x}p^{x}(1-p)^{n-x}\\
&= \sum_{x=1}^{n}x\cdot{}_{n}\mathrm{C}_{x}p^{x}(1-p)^{n-x}\\
&= \sum_{x=1}^{n}n\cdot{}_{n-1}\mathrm{C}_{x-1}p^{x}(1-p)^{n-x} \qquad& (\text{冒頭の二項係数の性質から})\\
&= np\sum_{x=1}^{n}{}_{n-1}\mathrm{C}_{x-1}p^{x-1}(1-p)^{n-x}\\
&= np\sum_{x=0}^{n-1}{}_{n-1}\mathrm{C}_{x}p^{x}(1-p)^{n-1-x} \qquad&(\text{x-1をxに置き換え})\\
&= np\{p+(1-p)\}^{n-1}\\
&= np.
\end{align}
$$

以上で、$B(n, p)$の期待値が$np$であることが証明できました。
つぎは分散$V(X)$ですが、$V(X)=E(X^2)-E(X)^2$を利用して算出するので、$E(X^2)$を計算していきます。
$$
\begin{align}
E(X^2) &= \sum_{x=0}^{n}x^2\cdot{}_{n}\mathrm{C}_{x}p^{x}(1-p)^{n-x}\\
&= \sum_{x=1}^{n}x^2\cdot{}_{n}\mathrm{C}_{x}p^{x}(1-p)^{n-x}\\
&= \sum_{x=1}^{n}xn\cdot{}_{n-1}\mathrm{C}_{x-1}p^{x}(1-p)^{n-x}\\
\end{align}
$$
ここで、$\sum$の中の最初の$x$を、$x=(x-1)+1$と変形して、2項にわけます。
$$
\begin{align}
E(X^2) &= \sum_{x=1}^{n}(x-1)n\cdot{}_{n-1}\mathrm{C}_{x-1}p^{x}(1-p)^{n-x}+\sum_{x=1}^{n}n\cdot{}_{n-1}\mathrm{C}_{x-1}p^{x}(1-p)^{n-x}\\
&= np\sum_{x=1}^{n}(x-1)\cdot{}_{n-1}\mathrm{C}_{x-1}p^{x-1}(1-p)^{n-x}+np\sum_{x=1}^{n}{}_{n-1}\mathrm{C}_{x-1}p^{x-1}(1-p)^{n-x}\\
&= np\sum_{x=0}^{n-1}x\cdot{}_{n-1}\mathrm{C}_{x}p^{x}(1-p)^{n-1-x}+np\sum_{x=0}^{n-1}{}_{n-1}\mathrm{C}_{x}p^{x}(1-p)^{n-1-x}
\end{align}
$$
ここで、$\sum_{x=0}^{n-1}x\cdot{}_{n-1}\mathrm{C}_{x}p^{x}(1-p)^{n-1-x}$は$B(n-1, p)$の期待値なので、$(n-1)p$です。
さらに、$\sum_{x=0}^{n-1}{}_{n-1}\mathrm{C}_{x}p^{x}(1-p)^{n-1-x}$は$B(n-1, p)$の確率関数の全体の和なので$1$になります。
よって、
$$
E(X^2) = n(n-1)p^2 + np
$$
となります。

あとはこれを使って、
$$
\begin{align}
V(X) &= E(X^2) – E(X)^2\\
&= n^2p^2-np^2+np-(np)^2\\
&= np(1-p)
\end{align}
$$
が導出されました。