x回目に初めて当たる確率：幾何(first sucess, geometric)分布メモ

x回目に初めて当たる確率分布である幾何(first sucess, geometric)分布のメモ

幾何分布とは

幾何分布は，その事象が成功する確率$p$として，確率変数$x$と置いたとき，$x$回試行してちょうど$x$回目に成功する確率を表現する確率分布になる．

ガチャの当たる確率などで試行回数がどのくらい必要になるかなどが計算できる．

$$ x \sim \mathrm{Geom}_{x}\left[ p \right] $$

幾何分布の確率質量関数pmf

$x$回目までの各回数での確率を考えると求めることができる． $x$回目に初めて成功するので，$x$回目の確率は$p$, 1回目から$x-1$回目まで各回で$(1-p)$の確率で失敗しているので，失敗の総確率は$(1-p)^{x-1}$になる. よって，確率質量関数は次のようになる．

$\begin{eqnarray} \mathrm{pmf}(x) = (1-p)^{x-1} \cdot p, (x \geq 1) \end{eqnarray}$

分布の特徴：

$x$回目に成功する確率なので，$x-1$回までは連続で失敗している，つまり，xが大きくなればなるほど確率密度は小さくなっていく．
どのタイミング（1回目やx回目）でも成功確率$p$自体は変わらない．:無記憶性

確率の公理を満たすか

等比数列の和で計算できる．

$\begin{eqnarray} \sum_{x=1}^{\infty} (1-p)^{x-1} \cdot p &=& p \sum_{x=1}^{\infty} 1 \cdot (1-p)^{x-1} \\\\ &=& p \lim_{x \to \infty} \frac{1\cdot (1 - ( 1 - p )^{x-1})}{ 1 - ( 1 - p ) } \\\\ &=& p \frac{1 - 0}{p} \\\\ &=& 1 \\\\ \end{eqnarray}$

例：コインを$x$回振って$x$回目に初めて表がでる確率分布

$p=\frac{1}{2}$ のときの幾何分布がこの問題の確率分布に該当する．

ここで，$x=3$の確率とは，2回裏が出て3回目に表が出る確率になる．

$\begin{eqnarray} \mathrm{pmf}(x) = \left( 1-\frac{1}{2} \right)^{x-1}\frac{1}{2} = \left( \frac{1}{2} \right)^{x} \end{eqnarray}$

import scipy as sp
import scipy.stats as stats
import matplotlib.pyplot as plt
%matplotlib inline

plt.figure()

N = 12
x = sp.arange(1, N+1)
g_p = stats.geom(p=1/2).pmf(x)
plt.bar(x, g_p)
plt.xlim(0, N+1)
plt.grid()
plt.xlabel("$x$")
plt.ylabel("$\mathrm{pmf}(x)$")
plt.savefig("coin_geom_pdf.png")
plt.show()

例：サイコロを$x$回振って$x$回目に初めて1がでる確率分布

$p=\frac{1}{6}$ のときの幾何分布がこの問題の確率分布に該当する．

ここで，$x=3$の確率とは，2回1以外が出て3回目に1が出る確率になる．

$\begin{eqnarray} \mathrm{pmf}(x) = \left( 1-\frac{1}{6} \right)^{x-1}\frac{1}{6} = \left( \frac{5}{6} \right)^{x-1}\frac{1}{6} \end{eqnarray}$

import scipy as sp
import scipy.stats as stats
import matplotlib.pyplot as plt
%matplotlib inline

plt.figure()
N = 30
x = sp.arange(1, N+1)
g_p = stats.geom(p=1/6).pmf(x)
plt.bar(x, g_p)
plt.grid()
plt.xlabel("$x$")
plt.ylabel("$\mathrm{pmf}(x)$")
plt.savefig("dice_geom_pdf.png")
plt.show()

累積密度関数cdf

こちらのほうがpmfより実用的である． pmfではx回目に成功する確率であったが， cdfの場合，累積しているので$x$回目までに1回は成功する確率になる．つまり，これを計算できれば成功するまでにどのくらい試行回数が必要になるかがわかる．

計算方法として，全確率1から$x$回目までに1度も成功しない余事象を引けば計算できる．

$\begin{eqnarray} \mathrm{cdf}(x) = 1 - ( 1 - p )^{x} \end{eqnarray}$

例：コインをx回目まで振って1回は表が出る確率

$\begin{eqnarray} \mathrm{cdf}(x) = 1 - ( 1 - \frac{1}{2} )^{x} = 1 - \frac{1}{2^{x}} \end{eqnarray}$

import scipy as sp
import scipy.stats as stats
import matplotlib.pyplot as plt
%matplotlib inline

plt.figure()

N = 12
x = sp.arange(1, N+1)

plt.figure()
g_cdf = stats.geom(p=1/2).cdf(x)
plt.bar(x, g_cdf)
plt.hlines(0.95, 0, N+1, colors="r", label="95%")
plt.xlim(0, N+1)
plt.legend()
plt.grid()
plt.tight_layout()
plt.savefig("coin_cdf.png")
plt.show()

5回振れば表が出る確率が95%を超えることがわかる．

例：サイコロをx回目まで振って1回は1が出る確率

$\begin{eqnarray} \mathrm{cdf}(x) = 1 - ( 1 - \frac{1}{6} )^{x} = 1 - \left( \frac{5}{6} \right)^{x} \end{eqnarray}$

import scipy as sp
import scipy.stats as stats
import matplotlib.pyplot as plt
%matplotlib inline

plt.figure()

N = 30
x = sp.arange(1, N+1)

plt.figure()
g_cdf = stats.geom(p=1/6).cdf(x)
plt.bar(x, g_cdf)
plt.hlines(0.95, 0, N+1, colors="r", label="95%")
plt.xlim(0, N+1)
plt.legend()
plt.grid()
plt.tight_layout()
plt.savefig("coin_cdf.png")
plt.show()

17回振れば1が出る確率が95%を超える．

幾何分布の期待値

$\begin{eqnarray} E\left[ x \right] &=& \sum_{x=1}^{\infty} x \mathrm{pmf}(x) \\\\ &=& \sum_{x=1}^{\infty} x ( 1 - p )^{x - 1}p \\\\ &=& p\sum_{x=1}^{\infty} x ( 1 - p )^{x - 1} \\\\ &=& p\left( 1 + 2 ( 1 - p ) + 3(1-p)^{2} + \cdots + \right) \\\\ \end{eqnarray}$

漸化式を作って解くと，

$\begin{eqnarray} E\left[x \right] &=& \sum_{x=1}^{\infty} x ( 1 - p )^{x - 1}p \\\\ &=& p\cdot 1 + \sum_{x=2}^{\infty} x( 1 - p )^{x - 1}p \\\\ &=& p\cdot 1 + (1-p)\sum_{x=2}^{\infty} (x-1+1)( 1 - p )^{(x - 1)-1}p \\\\ &=& p\cdot 1 + (1-p)\sum_{x=2}^{\infty} (x-1) ( 1 - p )^{(x-1) - 1}p + p(1-p)\sum_{x=2}^{\infty} 1\cdot ( 1 - p )^{x-2} \\\\ &=& p\cdot 1 + (1-p)\sum_{x=1}^{\infty} x ( 1 - p )^{x - 1}p + p(1-p)\frac{1}{1-(1-p)} \\\\ &=& (1-p)E\left[ x \right] + p\cdot 1 + (1-p) \cdot 1 \\\\ \end{eqnarray}$

この期待値の漸化式は，今までの期待値回数分失敗して，最後に成功する期待値$p$回，最後もまた失敗する期待値$1-p$回で幾何分布（試行回数）の期待値が決まることを示している．

これを解くと

$\begin{eqnarray} E\left[x \right] &=& (1-p)E\left[ x \right] + p + (1-p) \\\\ E\left[x \right] - (1-p)E\left[ x \right] &=& 1 \\\\ pE\left[x \right] &=& 1 \\\\ E\left[x \right] &=& \frac{1}{p} \\\\ \end{eqnarray}$

参考：https://mathtrain.jp/kikabunpu

※等差数列と等比数列の積から導出も可能

例：コインを投げてx回目に初めて表が出る試行回数の期待値

$\begin{eqnarray} E\left[ x \right] = \frac{1}{p} = \frac{1}{\frac{1}{2}} = 2 \end{eqnarray}$

そもそも数学的確率におけるコインを投げて表が出る確率$1/2$というのは，試行回数nを$\infty$まで繰り返すと平均的に2回に1回は表になるという意味なので幾何分布の期待値とこの意味は一致する．

import scipy as sp
import scipy.stats as stats
import matplotlib.pyplot as plt
%matplotlib inline

plt.figure()

p = 1/2
N = 12
x = sp.arange(1, N+1)
g_p = stats.geom(p=p).pmf(x)
plt.bar(x, g_p)

E = 1/p
plt.vlines(E, 0, 0.5, colors="r",label="E")

plt.xlim(0, N+1)
plt.xlabel("$x$")
plt.ylabel("$\mathrm{pmf}(x)$")
plt.grid()
plt.legend()
plt.savefig("coin_geom_pdf_E.png")
plt.show()

例：サイコロを投げてx回目に1が出る試行回数の期待値

$\begin{eqnarray} \mathrm{E}\left[ x \right] = \frac{1}{p} = \frac{1}{\frac{1}{6}} = 6 \end{eqnarray}$

import scipy as sp
import scipy.stats as stats
import matplotlib.pyplot as plt
%matplotlib inline

plt.figure()

p = 1/6
N = 30
x = sp.arange(1, N+1)
g_p = stats.geom(p=p).pmf(x)
plt.bar(x, g_p)

E = 1/p
plt.vlines(E, 0, 0.175, colors="r",label="E")

plt.xlim(0, N+1)
plt.xlabel("$x$")
plt.ylabel("$\mathrm{pmf}(x)$")
plt.grid()
plt.legend()
plt.savefig("dice_geom_pdf_E.png")
plt.show()

幾何分布の分散

$\begin{eqnarray} \mathrm{V}\left[ x \right] &=& \sum_{x=1}^{\infty} (x - E\left[ x \right])^{2} \mathrm{pmf}(x) \\\\ &=& E\left[ x^{2} \right] - (E\left[ x \right])^{2} \end{eqnarray}$

$E\left[ x^{2} \right]$を求める．

$\begin{eqnarray} E\left[ x^{2} \right] &=& \sum_{x=1}^{\infty} x^{2} ( 1 - p )^{x - 1}p \\\\ &=& p \sum_{x=1}^{\infty} x^{2} ( 1 - p )^{x - 1} \\\\ &=& p (1 + 2^{2}( 1 - p ) + 3^{2}(1-p)^{2} + \cdots + n^{2}(1-p)^{n-1} ) \end{eqnarray}$

$ (1 - p) E\left[ x^{2} \right] $ を考える．

$\begin{eqnarray} (1-p)E\left[ x^{2} \right] &=& p \sum_{x=1}^{\infty} x^{2} ( 1 - p )^{x} \\\\ &=& p( 1^{2}(1-p) + 2^{2}( 1 - p )^{2} + 3^{2}( 1 - p )^{3} + \cdots + n^{2}(1-p)^{n} ) \end{eqnarray}$

$\begin{eqnarray} E\left[ x^{2} \right] - (1-p)E\left[ x^{2} \right] &=& p\left(1 + (2^{2} - 1^{2})(1-p) + (3^{2} - 2^{2})(1-p)^{2} + \cdots + \left( ( n+1 )^{2} - n^{2}\right)(1-p)^{n} \right) \\\\ &=& p\left(1 + 3 \cdot 1 (1-p) + 5 \cdot 1(1-p)^{2} + \cdots + \left( (n+1+n)(n+1-n) \right)(1-p)^{n} \right) \\\\ &=& p\left(1 + 3 \cdot 1 (1-p) + 5 \cdot 1(1-p)^{2} + \cdots + (2n+1)(1-p)^{n} \right) \\\\ &=& p\sum_{x=0}^{\infty} (2x+1)(1-p)^{x} \\\\ &=& 2p\sum_{x=0}^{\infty} x(1-p)^{x} + p\sum_{x=0}^{\infty} (1-p)^{x} \end{eqnarray}$

1項目 $\sum_{x=1}^{\infty} x(1-p)^{x} = (1-p) + 2(1-p)^{2} + \cdots + n(1-p)^{n}$を求める．

$\begin{eqnarray} (1-p)\sum_{x=1}^{\infty} x(1-p)^{x} &=& (1-p)^{2} + 2(1-p)^{3} + \cdots + (n-1)(1-p)^{n} + n(1-p)^{n+1} \\\\ \sum_{x=1}^{\infty} x(1-p)^{x} - (1-p)\sum_{x=1}^{\infty} x(1-p)^{x} &=& (1-p) + (1-p)^{2} + (1-p)^{3} + \cdots + (1-p)^{n} \\\\ &=& \frac{(1-p)}{1 - (1-p)} \\\\ p\sum_{x=1}^{\infty} x(1-p)^{x} &=& \frac{(1-p)}{1 - (1-p)} \\\\ \sum_{x=1}^{\infty} x(1-p)^{x} &=& \frac{1-p}{p^{2}} \end{eqnarray}$

よって，

$\begin{eqnarray} \mathrm{E}\left[ x^{2} \right] - (1-p)E\left[ x^{2} \right] &=& 2p \frac{1-p}{p^{2}} + p\frac{1}{1-(1-p)} \\\\ &=& \frac{2(1-p)}{p} + p\frac{1}{p} \\\\ p E\left[ x^{2} \right] &=& \frac{2(1-p)}{p} + 1 \\\\ E\left[ x^{2} \right] &=& \frac{2(1-p)}{p^{2}} + \frac{1}{p} \\\\ &=& \frac{2-p}{p^{2}} \\\\ \end{eqnarray}$

よって，分散は，

$\begin{eqnarray} \mathrm{V}\left[ x \right] &=& E\left[ x^{2} \right] - (E\left[ x \right])^{2} \\\\ &=& \frac{2-p}{p^{2}} - \left( \frac{1}{p} \right)^{2} \\\\ &=& \frac{1 - p}{p^{2}} \\\\ \end{eqnarray}$

ルートを取った標準偏差を考えると，

$\begin{eqnarray} \sigma = \pm \sqrt{ \mathrm{V}\left[ x \right] } = \pm \frac{\sqrt{1-p}}{p} \end{eqnarray}$

例：コインを投げてx回目に表が出る試行回数の分散

$\begin{eqnarray} \mathrm{V}\left[ x \right] &=& \frac{1-p}{p^{2}} \\\\ &=& \frac{1-\frac{1}{2}}{\frac{1}{2}^{2}} \\\\ &=& 2 \\\\ \end{eqnarray}$

標準偏差

$\begin{eqnarray} \sigma = \pm \sqrt{ \mathrm{V}\left[ x \right] } &=& \pm \frac{\sqrt{1-p}}{p} \\\\ &=& \pm \sqrt{2} \\\\ \end{eqnarray}$

$\begin{eqnarray} E \pm \sigma &=& 2 \pm \sqrt{ 2 } &=& \left( 0.586, 3.414 \right) \\\\ \end{eqnarray}$

$E - \sigma, E+\sigma$ 範囲内の確率は，

$\begin{eqnarray} P(E-\sigma \leq x \leq E+\sigma ) = P(0.586 \leq x \leq 3.414) &=& \mathrm{cdf}(3) - \mathrm{cdf}(0) \\\\ &=& 1-(1-\frac{1}{2})^{3} - 0 \\\\ &=& 0.875 \\\\ \end{eqnarray}$

pmfの図にこの区間を可視化すると，

import scipy as sp
import scipy.stats as stats
import matplotlib.pyplot as plt
%matplotlib inline

p = 1/2
N = 12
g_dist = stats.geom(p=p)

E = 1/p
V = (1-p)/p**2
sigma = sp.sqrt(V)
print("V", V, "sigma", sigma)
print("E-s", E-sigma, "E+s", E+sigma)

l = sp.floor(E-sigma)
h = sp.floor(E+sigma)
ra = sp.arange(l, h+1)
p_ra = g_dist.cdf(h) - g_dist.cdf(l)
print(f"p({l+1}<= x <={h}):{p_ra}")

plt.figure()
x = sp.arange(1, N+1)
plt.bar(x, g_dist.pmf(x))

plt.vlines(E, 0, 0.5, colors="r",label="E")
plt.vlines([E-sigma, E+sigma], 0, 0.5, colors="g",label="$E\pm \sigma$")

plt.xlim(0, N+1)
plt.xlabel("$x$")
plt.ylabel("$\mathrm{pmf}(x)$")
plt.grid()
plt.legend()
plt.savefig("coin_geom_pdf_E_std.png")
plt.show()

例：サイコロを投げてx回目に1が出る試行回数の分散

$\begin{eqnarray} \mathrm{V}\left[ x \right] &=& \frac{1-p}{p^{2}} \\\\ &=& \frac{1-\frac{1}{6}}{\frac{1}{6}^{2}} \\\\ &=& 30 \\\\ \end{eqnarray}$

標準偏差

$\begin{eqnarray} \sigma = \pm \sqrt{ \mathrm{V}\left[ x \right] } &=& \pm \frac{\sqrt{1-p}}{p} \\\\ &=& \pm \sqrt{30} = 5.477\\\\ \end{eqnarray}$

$\begin{eqnarray} E \pm \sigma &=& 6 \pm \sqrt{ 30 } &=& \left( 0.52, 11.48 \right) \\\\ \end{eqnarray}$

$E-\sigma, E+\sigma$ 範囲内の確率は，

$\begin{eqnarray} P(E-\sigma \leq x \leq E+\sigma ) = P(0.52 \leq x \leq 11.48) &=& \mathrm{cdf}(11) - \mathrm{cdf}(0) \\\\ &=& 1-(1-\frac{1}{6})^{11} - 0 \\\\ &=& 0.865 \\\\ \end{eqnarray}$

pmfの図にこの区間を可視化すると，

import scipy as sp
import scipy.stats as stats
import matplotlib.pyplot as plt
%matplotlib inline

p = 1/2
N = 12
g_dist = stats.geom(p=p)

E = 1/p
V = (1-p)/p**2
sigma = sp.sqrt(V)
print("V", V, "sigma", sigma)
print("E-s", E-sigma, "E+s", E+sigma)

l = sp.floor(E-sigma)
h = sp.floor(E+sigma)
ra = sp.arange(l, h+1)
p_ra = g_dist.cdf(h) - g_dist.cdf(l)
print(f"p({l+1}<= x <={h}):{p_ra}")

plt.figure()
x = sp.arange(1, N+1)
plt.bar(x, g_dist.pmf(x))

plt.vlines(E, 0, 0.5, colors="r",label="E")
plt.vlines([E-sigma, E+sigma], 0, 0.5, colors="g",label="$E\pm \sigma$")

plt.xlim(0, N+1)
plt.xlabel("$x$")
plt.ylabel("$\mathrm{pmf}(x)$")
plt.grid()
plt.legend()
plt.savefig("coin_geom_pdf_E_std.png")
plt.show()

幾何分布のパラメータ推定

参考：幾何分布の最尤推定について解説 - 理数アラカルト -

最尤推定でパラメータ$p$を求める．

観測データが${ x_{1}, \cdots x_{N} }$のN個の場合の尤度は，

$\begin{eqnarray} L(p) &=& \Pi_{i=1}^{N} (1-p)^{x_{i}-1}p \\\\ &=& p^{N} \Pi_{i=1}^{N} (1-p)^{x_{i}-1} \end{eqnarray}$

対数尤度LL(p)は

$\begin{eqnarray} LL(p) &=& \sum_{i=1}^{N} \log{\left( (1-p)^{x_{i}-1}p \right) } \\\\ &=& N\log{p} + \log{(1-p)} \sum_{i=1}^{N} (x_{i}-1) \end{eqnarray}$

微分して0となる点が最大のpになるので，

$\begin{eqnarray} \frac{\partial LL(p)}{\partial p} &=& 0 \\\\ \frac{N}{p} - \frac{\sum_{i=1}^{N} (x_{i}-1)}{1-p} &=& 0 \\\\ (1-p)N - p\sum_{i=1}^{N} (x_{i}-1) &=& 0 \\\\ (1-p)N - p(\sum_{i=1}^{N} x_{i} - N) &=& 0 \\\\ N - p\sum_{i=1}^{N} x_{i} &=& 0 \\\\ p\sum_{i=1}^{N} x_{i} &=& N \\\\ p &=& \frac{N}{\sum_{i=1}^{N} x_{i}} \\\\ &=& \frac{1}{ \frac{1}{N} \sum_{i=1}^{N} x_{i}} \\\\ &=& \frac{1}{ \bar{x} } \\\\ \end{eqnarray}$