はじめに
確率分布間の関係の理解を目的として、ベルヌーイ分布を起点に、様々な確率分布を生成する関数をPythonで実装し、まとめました。本記事は以下の方向けです。
- 基本的な確率分布は一通り勉強したことがある
- Pythonを触ったことがある
- 確率分布間の関係を理解したい
実装したPythonの関数は以下の通りです。bernolli、binomなどが関数名です。
例えば、一様分布は次の順にたどることで生成できます。
ベルヌーイ分布 -> 幾何分布 -> 指数分布 -> ガンマ分布 -> ベータ分布 -> 一様分布
関数名や確率分布の定義はscipy.statsを踏襲しています。
それでは、これから各関数を順に説明していきます。
ベルヌーイ分布
まず、出発点のベルヌーイ分布に従うサンプルを生成します。ベルヌーイ分布とは、「当たりが出る確率が$p$のくじを1回引いたときの当たり回数」に従う確率分布です。定義からも分かる通り、とりうる値は0か1です。確率関数は下式で与えられます。
f(x;p) = p^x (1-p)^{1-x}\ \ (x=0,1)
Python実装に進みます。あらかじめnumpyとscipy.statsをインポートしておきます。
import numpy as np
from scipy import stats
確率が$p$のベルヌーイ分布に従うサンプルをsampleSize個発生させる関数は下の通りです。
def bernoulli(p, sampleSize):
return stats.bernoulli.rvs(p, size=sampleSize)
試しにこの関数から20個のサンプルを生成してみます。
bernoulli(0.2, 20)
> array([1, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1])
設定どおり、約0.2の確率で1が出現しています。
scipy.statsを使うのはこれが最後で、以後このbernoulli関数を元に確率分布を生成します。
ベルヌーイ分布から二項分布
二項分布とは、「当たりが出る確率が$p$のくじを$n$回引いた時の当たり回数」に従う確率分布です。確率関数は下式で与えられます。
f(x;n,p) = \left(
\begin{matrix}
n\\
x
\end{matrix}
\right) p^x(1-p)^{n-x}\ \ (x=0,1,\cdots,n)
二項分布はベルヌーイ分布における「くじを引く回数」を$n$倍にしたものです。sampleSize個の「試行回数$n$、確率$p$の二項分布から生成されるサンプル」は、$n\times$sampleSize個の「確率$p$のベルヌーイ分布から生成されるサンプル」から作ることができます。実装は下の通りです。
def binom(n, p, sampleSize):
x = bernoulli(p, n * sampleSize)
x = x.reshape([n, sampleSize])
return np.sum(x, axis=0)
二項分布からポアソン分布
ポアソン分布とは「単位時間あたり平均$\mu$個の当たりが出るくじにおいて、単位時間に出る当たりの数」に従う確率分布です。確率関数は下式で与えられます。
f(x;\mu) = \frac{\mu^x}{x!}\exp(-\mu)\ \ (x=0,1,\cdots)
下の通り、ポアソン分布は二項分布とよく似ています。
| 二項分布 | ポアソン分布 | |
|---|---|---|
| 何を見る？ | 当たり回数 | 当たり回数 |
| どういう条件で | 試行回数一定 | 時間を一定 |
試行回数$n$、確率$\frac{\mu}{n}$の二項分布の$n$を$n\to\infty$とすると、平均$\mu$のポアソン分布に収束します。実装は下の通りです。
def poisson(mu, sampleSize):
n = 1000 # 十分大きい数
p = mu / n
return binom(n, p, sampleSize)
$n$は理論的には大きいければ大きいほど良いですが、計算量が増大します。
ポアソン分布に従うサンプルはbinom関数から生成しました。そのbinom関数はbernoulli関数が元になっています。つまり、ベルヌーイ分布からポアソン分布を生成したことになります。今後もこのような方法をとっていきます。
ベルヌーイ分布から幾何分布
幾何分布とは、「当たりが出る確率が$p$のくじを当たるまで引いたときの試行回数」に従う確率分布です。確率関数は下式で与えられます。
f(x;p)= p(1-p)^{x-1}\ \ (x=1,2,\cdots)
確率$p$の幾何分布から生成されるサンプルは、確率$p$のベルヌーイ分布から作ることができます。実装は下の通りです。
def geom(p, sampleSize):
x = []
for _ in range(sampleSize):
t = 1
while bernoulli(1 - p, 1):
t += 1
x.append(t)
return np.array(x)
def geom(p, sampleSize):
n = 1000
x = []
s = np.array([-1], dtype=int)
while len(x) < sampleSize:
x_ = bernoulli(p, n) # ベルヌーイ分布に従うサンプルを同時にnサンプル生成する。
x_ = np.where(x_ == 1)[0]
x_ = np.concatenate([s, x_])
x.extend(x_[1:] - x_[:-1])
s = np.array([x_[-1] - n])
return np.array(x[:sampleSize])
2つの関数は同じ処理をしていますが、前者のほうがわかりやすく、後者のほうが速いと思います。
幾何分布から負の二項分布
負の二項分布とは、「当たりが出る確率が$p$のくじを$n$回当たるまで引いたときの、外れ回数」に従う確率分布です。確率関数は下式で与えられます。
f(x;n,p) = \left(
\begin{matrix}
n+x-1\\
x
\end{matrix}
\right) p^n(1-p)^{x}\ \ (x=0,1,\cdots,n)
幾何分布との本質的な違いは、「当たりが何回出るまでくじを引くか」のみです。下の通り、負の二項分布は幾何分布から作ることができます。
def nbinom(n, p, sampleSize):
x = geom(p, sampleSize * n) - 1
x = x.reshape([sampleSize, n])
return np.sum(x, axis=1)
幾何分布から指数分布
指数分布とは、「単位時間あたり平均で$\frac{1}{\lambda}$個の当たりが出るくじにおいて、1回当たりが出るまでの時間」に従う確率分布です。確率密度関数は下式で与えられます。
f(x,\lambda) = \lambda\exp(-\lambda x)
指数分布は次のように幾何分布とよく似ています。
| 幾何分布 | 指数分布 | |
|---|---|---|
| 何を見る？ | 1回当たりが出るまでの試行数 | 1回当たりが出るまでの時間 |
| どういうくじ？ | 試行毎の当たり確率が一定 | 単位時間あたりの平均当たり回数が一定 |
時間や試行によらず、当たり確率が変わらない性質を無記憶性と呼びます。幾何分布も指数分布もこの性質を持っています。
幾何分布において$p\to0$とすることで、指数分布を作ることができます。
def expon(lam, sampleSize):
p = 0.01 # 十分小さい数
x = geom(p, sampleSize)
return x * p * lam
$p$は理論的には小さければ小さいほど良いですが、計算量が増大します。
指数分布からポアソン分布
ポアソン分布は、二項分布だけでなく指数分布からも作ることができます。指数分布とポアソン分布のいけない関係のスライドが、非常に分かりやすいです。スライドにはRの実装が載っていますが、ここではPythonでの実装を示します。
def poisson_(mu, sampleSize):
x = []
while len(x) < sampleSize:
t = 0
n = -1
while t <= 1:
t += expon(1 / mu, 1)
n += 1
x.append(n)
return np.array(x)
指数分布からガンマ分布
ガンマ分布とは、「単位時間あたり平均で$\frac{1}{\lambda}$個当たりが出るくじにおいて、$\alpha$回当たりが出るまでの時間」に従う確率分布です。確率密度関数は下式で与えられます。
f(x;\alpha,\lambda)=\frac{\lambda^\alpha}{\Gamma(\alpha)}x^{\alpha-1}\exp(-\lambda x)\ \ (x\ge0)
ただし、$\Gamma(\alpha)$は下式で与えられます。
\Gamma(\alpha)=\int_0^\infty t^{\alpha-1}\exp(-t)dt
指数分布との本質的な違いは、「当たりが何個出るまでの時間を見るか」のみです。ガンマ分布は指数分布から作ることができます。
def gamma(alpha, lam, sampleSize):
x = expon(lam, sampleSize * alpha)
x = x.reshape([sampleSize, alpha])
return np.sum(x, axis=1)
尚、上記gamma関数の引数$\alpha$は自然数である必要があります。一般のガンマ分布の$\alpha$は正の数であればなんでもOKなのですが、自然数以外の$\alpha$だとベルヌーイ分布から作れないですし、説明も難しいので、ここでは自然数に限定しました。
負の二項分布からガンマ分布
ガンマ分布は負の二項分布からも作ることができます。次のようにガンマ分布と負の二項分布はよく似ています。
| 負の二項分布 | ガンマ分布 | |
|---|---|---|
| 何を見る？ | $n$回当たりが出るまでの試行数 | $\alpha$回当たりが出るまでの時間 |
| どういうくじ？ | 試行毎の当たり確率が一定 | 単位時間あたりの平均当たり回数が一定 |
負の二項分布からガンマ分布を作る方法は、幾何分布から指数分布を作る方法と全く同じです。
def gamma_(alpha, lam, sampleSize):
p = 0.01 # 十分小さい数
x = nbinom(alpha, p, sampleSize)
return x * p * lam
ここでも、引数$\alpha$は自然数である必要があります。
二項分布から正規分布
どのような分布も無限に足し合わせると正規分布になります（中心極限定理）。私が初めて作った記事「確率分布の漸近的性質をPythonで確かめる」に例があります。
二項分布はベルヌーイ分布の$n$個の足し合わせでした。$n\to\infty$とすることで二項分布は正規分布に従います。平均$\mu$、分散$\sigma^2$の正規分布の確率密度関数は下式で与えられます。
f(x;\mu,\sigma^2)=\frac{1}{\sqrt{2\pi\sigma^2}}\exp\biggl\{-\frac{(x-\mu)^2}{2\sigma^2}\biggr\}
平均$\mu$、標準偏差$\sigma$の正規分布に従うサンプルは、次の方法で生成することができます。
def norm(mu, sigma, sampleSize):
n = 1000 # 十分大きい数
p = 0.5
x = binom(n, p, sampleSize)
sd = (n * p * (1 - p)) ** 0.5
x = (x - n * p) / sd * sigma + mu
return x
正規分布からカイ二乗分布
くじのような説明はできませんが、仮説検定でよく出てくるのが$\chi^2$分布です。$X_1,X_2,\cdots,X_{df}$が、独立に標準正規分布（平均0、分散1の正規分布）に従うとき、$X_1^2+X_2^2+\cdots+X_{df}^2$を自由度$df$の$\chi^2$分布と言います。実装は下の通りです。
def chi2(df, sampleSize):
x = norm(0, 1, sampleSize * df) ** 2
x = x.reshape([sampleSize, df])
return np.sum(x, axis=1)
カイ二乗分布からガンマ分布
これも、くじのような説明できませんが、自由度1の$\chi^2$分布は、$\alpha=1/2$、$\lambda=1/2$のガンマ分布に一致します。ガンマ分布の再生成という性質より、$\alpha$が$\frac{1}{2}$の$n$倍（$n$は自然数）の場合のガンマ分布を生成できます。再生性については、最後に説明します。
def gamma__(alpha, lam, sampleSize):
df = int(np.round(alpha * 2))
x = chi2(1, sampleSize * df)
x = x.reshape([sampleSize, df])
x = x * lam / 2
return np.sum(x, axis=1)
この関数の引数$\alpha$は$\frac{1}{2}$の$n$倍（$n$は自然数）である必要があります。
ガンマ分布からベータ分布
これもくじのような説明ができませんが、「パラメータが$(\alpha,\lambda)$のガンマ分布」と「パラメータが$(\beta,\lambda)$のガンマ分布」をもとに、パラメータが$(\alpha,\beta)$のベータ分布を生成することができます。
ベータ分布の密度関数は下式で与えられます。
f(x;\alpha,\beta)=x^{\alpha-1}(1-x)^{\beta-1}\ \ (0\le x\le 1)
取りうる範囲が0～1であるため、ベイズ統計学において、確率を表すパラメータの事前分布としてよく用いられます。
def beta(alpha, beta, sampleSize):
x1 = gamma(alpha, 1, sampleSize)
x2 = gamma(beta, 1, sampleSize)
return x1 / (x1 + x2)
ベータ分布から一様分布
ベータ分布の密度関数から明らかなとおり、$\alpha,\beta＝1$の場合、ベータ分布は一様分布と一致します。一様分布に従うサンプルは、ベータ分布から作ることができます。
def uniform(sampleSize):
return beta(1, 1, sampleSize)
まとめ
互いの関係が一目瞭然です。
ここで、再生性について説明します。再生性とは分布の足し算を、パラメータの足し算に置き換えることができる性質のことです。例えば、
- 試行回数が$n$、確率$p$の二項分布をX
- 試行回数が$m$、確率$p$の二項分布をY
とすると、$X+Y$は、試行回数が$(n+m)$、確率$p$の二項分布に従います。1.は、$n$個のベルヌーイ分布の和、2.は$m$個のベルヌーイ分布の和に等しいので、$X+Y$は$(n+m)$個のベルヌーイ分布の和、すなわち、試行回数が$(n+m)$の二項分布と一致するからです。
同じ理論で、上図の複数化の矢印の先の分布には、必ず再生性の性質があります（二項分布、負の二項分布、ガンマ分布、$\chi^2$分布）。そして、二項分布の極限である、正規分布とポアソン分布にも当然、再生性があります。
最後にPython実験
ベルヌーイ分布 -> 幾何分布 -> 指数分布 -> ガンマ分布 -> ベータ分布 -> 一様分布 -> ベルヌーイ分布を試してみます。
あらかじめプロットのための関数を定義しておきます。
import matplotlib.pyplot as plt
import seaborn as sns
sns.set(style='whitegrid')
def plot_descreteDitribution(sample, truePf, xmin, xmax):
sampleSize = len(sample)
x = np.arange(xmin, xmax + 1)
pf = np.array([np.sum(sample == i) for i in x]) / sampleSize # 実験で求めた分布
# 実験で求めた分布を青色でプロット
plt.plot(x - 0.1, pf, 'bo', ms=8)
plt.vlines(x - 0.1, 0, pf, colors='b', lw=5, alpha=0.5)
# 真の分布を赤色でプロット
plt.plot(x + 0.1, truePf(x), 'ro', ms=8)
plt.vlines(x + 0.1, 0, truePf(x), colors='r', lw=5, alpha=0.5)
def plot_continuousDitribution(sample, truePf, xmin, xmax):
sampleSize = len(sample)
x = np.linspace(xmin, xmax, 100)
# 実験で求めた分布を青色でプロット
th = np.linspace(xmin, xmax, 30)
hi = np.array([np.sum(np.logical_and(th[i] < sample, sample <= th[i + 1])) for i in range(30 - 1)]) / (sampleSize * (th[1] - th[0]))
plt.bar((th[:-1] + th[1:]) / 2, hi, width=(th[1] - th[0]))
# 真の分布を赤色でプロット
plt.plot(x, truePf(x), 'r', linewidth=4)
plt.xlim(xmin, xmax)
幾何分布(from ベルヌーイ分布)
p = 0.2
xmin = 1
xmax = 20
sample = geom(p, sampleSize)
truePf = lambda x: stats.geom.pmf(x, p)
plot_descreteDitribution(sample, truePf, xmin, xmax)
- 結果
赤色がscipy.statsで作った分布、青色が今回作成した関数から作った分布です。尚、sampleSizeは10000としました。以後も同様です。
指数分布(from ベルヌーイ分布 -> 幾何分布)
lam = 2
xmin = 0
xmax = 10
sample = expon(lam, sampleSize)
truePf = lambda x: stats.expon.pdf(x, scale=lam)
plot_continuousDitribution(sample, truePf, xmin, xmax)
- 結果
ガンマ分布(from ベルヌーイ分布 -> 幾何分布 -> 指数分布)
alpha = 3
lam = 2
xmin = 0
xmax = 20
sample = gamma(alpha, lam, sampleSize)
truePf = lambda x: stats.gamma.pdf(x, alpha, scale=lam)
plot_continuousDitribution(sample, truePf, xmin, xmax)
- 結果
ベータ分布(from ベルヌーイ分布 -> ... -> ガンマ分布)
alpha = 3
bet = 2
xmin = 0
xmax = 1
sample = beta(alpha, bet, sampleSize)
truePf = lambda x: stats.beta.pdf(x, alpha, bet)
plot_continuousDitribution(sample, truePf, xmin, xmax)
- 結果
一様分布(from ベルヌーイ分布 -> ... -> ベータ分布)
xmin = 0
xmax = 1
sample = uniform(sampleSize)
truePf = lambda x: stats.uniform.pdf(x)
plot_continuousDitribution(sample, truePf, xmin, xmax)
- 結果
ベルヌーイ分布(from ベルヌーイ分布 -> ... -> 一様分布)！？
p = 0.2
xmin = -2
xmax = 4
def bernoulli_(p, sampleSize):
x = uniform(sampleSize)
return (x < p).astype(int)
sample = bernoulli_(p, sampleSize)
truePf = lambda x: stats.bernoulli.pmf(x, p)
plot_descreteDitribution(sample, truePf, xmin, xmax)
- 結果
上手くいったことがわかります。