はじめに
確率分布間の関係の理解を目的として、ベルヌーイ分布を起点に、様々な確率分布を生成する関数をPythonで実装し、まとめました。本記事は以下の方向けです。
- 基本的な確率分布は一通り勉強したことがある
- Pythonを触ったことがある
- 確率分布間の関係を理解したい
実装したPythonの関数は以下の通りです。bernolli、binomなどが関数名です。
例えば、一様分布は次の順にたどることで生成できます。
ベルヌーイ分布 -> 幾何分布 -> 指数分布 -> ガンマ分布 -> ベータ分布 -> 一様分布
関数名や確率分布の定義はscipy.statsを踏襲しています。
それでは、これから各関数を順に説明していきます。
ベルヌーイ分布
まず、出発点のベルヌーイ分布に従うサンプルを生成します。ベルヌーイ分布とは、「当たりが出る確率が$p$のくじを1回引いたときの当たり回数」に従う確率分布です。定義からも分かる通り、とりうる値は0か1です。確率関数は下式で与えられます。
f(x;p) = p^x (1-p)^{1-x}\ \ (x=0,1)
ベルヌーイ分布の平均は$p$、分散は$p(1-p)$です。
[確率関数からの導出]
\begin{align}
E(X) &= 0\times(1-p) + 1\times p = p \\
V(X) &= E[(X-p)^2] \\ &= p^2 \times (1-p) + (1-p)^2 \times p \\&= p(1-p)
\end{align}
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)
二項分布の平均は$np$、分散は$np(1-p)$です。
[確率関数からの導出]
\begin{align}
E(X) &= \sum_{x=0}^n x\left(\begin{matrix}
n\\
x
\end{matrix}\right)p^x(1-p)^{n-x} \\
&= np\sum_{x=1}^n \left(\begin{matrix}
n-1\\
x-1
\end{matrix}\right)p^{x-1}(1-p)^{n-x} \\
&= np\sum_{y=0}^{n-1} \left(\begin{matrix}
n-1\\
y
\end{matrix}\right)p^{y}(1-p)^{n-1-y} \\
&= np \\
E[X(X-1)] &= \sum_{x=0}^n x(x-1) \left(\begin{matrix}
n\\
x
\end{matrix}\right)p^x(1-p)^{n-x} \\
&=
n(n-1)p^2\sum_{x=2}^n \left(\begin{matrix}
n - 2\\
x - 2
\end{matrix}\right)p^{x-2}(1-p)^{n-x} \\
&=
n(n-1)p^2\sum_{y=0}^{n-2} \left(\begin{matrix}
n - 2\\
y
\end{matrix}\right)p^{y}(1-p)^{n-2-y} \\
&=n(n-1)p^2 \\
V(X) &= E(X^2)-E(X)^2 \\ &= E[X(X-1)] + E(X) - E(X)^2 \\ &= n(n-1)p^2 + np - (np)^2 \\ &= np(1-p)
\end{align}
二項分布はベルヌーイ分布における「くじを引く回数」を$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)
ポアソン分布の平均、分散はともに$\mu$です。
[確率関数からの導出]
\begin{align}
E(X) &= \sum_{x=0}^\infty x\cdot\frac{\mu^x}{x!}\exp(-\mu) \\
&= \sum_{x=1}^\infty \cdot\frac{\mu\cdot\mu^{x-1}}{(x-1)!}\exp(-\mu) \\
&=\mu
\\
E[X(X-1)] &= \sum_{x=0}^{\infty} x(x-1) \cdot \frac{\mu^x}{x!}\exp(-\mu) \\
&= \sum_{x=2}^\infty \cdot\frac{\mu^2\cdot\mu^{x-2}}{(x-2)!}\exp(-\mu) \\
&=\mu^2
\\
V(X) &= E(X^2)-E(X)^2 \\
&= E[X(X-1)] + E(X) - E(X)^2 \\
&= \mu^2 + \mu - \mu^2 \\
&= \mu
\end{align}
下の通り、ポアソン分布は二項分布とよく似ています。
| 二項分布 | ポアソン分布 | |
|---|---|---|
| 何を見る？ | 当たり回数 | 当たり回数 |
| どういう条件で | 試行回数一定 | 時間を一定 |
試行回数$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)
幾何分布の平均は$\frac{1}{p}$、分散は$\frac{1-p}{p^2}$です。
[確率関数からの導出]
E(X) = \sum_{x=1}^\infty x p(1-p)^{x-1}
これを$S$とおく。
\begin{align}
S &= \sum_{x=1}^\infty x p(1-p)^{x-1} \cdots(1)\\
&= \sum_{x=0}^\infty (x+1) p(1-p)^x \\
&= p+\sum_{x=1}^\infty (x+1) p(1-p)^x \cdots(2)\\
\end{align}
$(1-p)\times(1)$をすると、
(1-p)S = \sum_{x=1}^\infty x p(1-p)^x \cdots(3)\\
$(2)-(3)$をすると
pS=p+ \sum_{x=1}^\infty p(1-p)^x = 1
よって
E(X)=S=\frac{1}{p}
次に、$E(X^2)$を求める
E(X^2)=\sum_{x=1}^\infty x^2 p(1-p)^{x-1}
これを$T$とおくと
\begin{align}
T&=\sum_{x=1}^\infty x^2 p(1-p)^{x-1} \cdots(4)\\
&=\sum_{x=0}^\infty (x+1)^2 p(1-p)^x \\
&=p+\sum_{x=1}^\infty (x+1)^2 p(1-p)^x \cdots(5)
\end{align}
$(1-p)\times(4)$をすると、
(1-p)T=\sum_{x=1}^\infty x^2 p(1-p)^x \cdots(6)
$(5)-(6)$をすると、
pT = p+\sum_{x=1}^\infty(2x+1)p(1-p)^x
$U=\sum_{x=1}^\infty(2x+1)p(1-p)^x$とおくと、
\begin{align}
U&=\sum_{x=1}^\infty(2x+1)p(1-p)^x \cdots(7)\\
&=\sum_{x=0}^\infty(2x+3)p(1-p)^{x+1} \\
&=3p(1-p)+\sum_{x=1}^\infty(2x+3)p(1-p)^{x+1} \cdots(8) \\
\end{align}
$(1-p)\times(7)$をすると、
(1-p)U=\sum_{x=1}^\infty(2x+1)p(1-p)^{x+1} \cdots(9)
$(8)-(9)$をすると、
\begin{align}
pU&=3p(1-p)+\sum_{x=1}^\infty 2p(1-p)^{x+1}\\
&=3p(1-p)+2(1-p)^2\\
&=-p^2-p+2
\end{align}
以上より、
\begin{align}
E(X^2)&=T=\frac{p+S}{p}=\frac{2-p}{p^2} \\
V(X)&=E(X^2)-E(X)^2= \frac{2-p}{p^2}-\frac{1}{p^2}=\frac{1-p}{p^2}
\end{align}
確率$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)
負の二項分布の平均は$\frac{n(1－p)}{p}$、分散は$\frac{n(1-p)}{p^2}$です。
[確率関数からの導出]
\begin{align}
E(X) &= \sum_{x=0}^\infty x
\left(\begin{matrix}n+x-1\\x\end{matrix}\right)
p^n(1-p)^{x}\\
&= n\sum_{x=1}^\infty \left(\begin{matrix}n+x-1\\x-1\end{matrix}\right)p^n(1-p)^{x}\\
&= n\sum_{y=0}^\infty \left(\begin{matrix}n+y\\y\end{matrix}\right)p^n(1-p)^{y+1}\\
&= \frac{n(1-p)}{p}\sum_{y=0}^\infty \left(\begin{matrix}n+y\\y\end{matrix}\right)p^{n+1}(1-p)^{y}\\
&=\frac{n(1-p)}{p} \\
E[X(X-1)] &= \sum_{x=0}^\infty x(x-1)
\left(\begin{matrix}n+x-1\\x\end{matrix}\right)
p^n(1-p)^{x}\\
&=n(n+1)\sum_{x=2}^\infty
\left(\begin{matrix}n+x-1\\x-2\end{matrix}\right)
p^n(1-p)^{x}\\
&=n(n+1)\frac{(1-p)^2}{p^2}\sum_{y=0}^\infty
\left(\begin{matrix}n+y+1\\y\end{matrix}\right)
p^{n+2}(1-p)^{y}\\
&=\frac{n(n+1)(1-p)^2}{p^2}\\
V(X) &= E(X^2)-E(X)^2 \\
&= E[X(X-1)] + E(X) - E(X)^2 \\
&= \frac{n(n+1)(1-p)^2}{p^2} + \frac{n(1-p)}{p} - \frac{n^2(1-p)^2}{p^2} \\
&= \frac{n(1-p)}{p^2}
\end{align}
「負の二項分布」という名前の由来は、確率関数を変形すると二項分布に類似した形になるからです。
[確率関数の変形]
\begin{align}
\left(\begin{matrix}n+x-1\\x\end{matrix}\right) p^n(1-p)^{x}&=\frac{(n+x-1)(n+x-2)\cdots(n+1)n}{x!}p^n(1-p)^x\\
&=\frac{-n(-n-1)\cdots(-n-x+1)}{x!}(-1)^xp^n(1-p)^x \\
&=\frac{-n(-n-1)\cdots(-n-x+1)}{x!}p^n(p-1)^x\\
&=\left(\begin{matrix}-n\\x\end{matrix}\right)p^n(p-1)^x
\end{align}
幾何分布との本質的な違いは、「当たりが何回出るまでくじを引くか」のみです。下の通り、負の二項分布は幾何分布から作ることができます。
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)
指数分布の平均は$\frac{1}{\lambda}$、分散は$\frac{1}{\lambda^2}$です。
[確率密度関数からの導出]
\begin{align}
E(X)
&=\int_0^\infty \lambda x \exp(-\lambda x)dx \\
&= - \int_0^\infty x \{ \exp(-\lambda x)\}^{'} dx \\
&= \int_0^\infty \exp(-\lambda x)dx \\
&= \frac{1}{\lambda} \\
E(X^2)
&=\int_0^\infty \lambda x^2 \exp(-\lambda x)dx \\
&= - \int_0^\infty x^2 \{ \exp(-\lambda x)\}^{'} dx \\
&= \int_0^\infty 2x \exp(-\lambda x) dx \\
&= \frac{2}{\lambda} E(X) \\
&= \frac{2}{\lambda^2} \\
V(X)
& = E(X^2) - E(X)^2\\
& = \frac{2}{\lambda^2} - \frac{1}{\lambda^2}\\
&= \frac{1}{\lambda^2}
\end{align}
指数分布は次のように幾何分布とよく似ています。
| 幾何分布 | 指数分布 | |
|---|---|---|
| 何を見る？ | 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
ガンマ分布の平均は$\frac{\alpha}{\lambda}$、分散は$\frac{\alpha}{\lambda^2}$です。
[確率密度関数からの導出]
\begin{align}
E(X)
&=\int_0^\infty \frac{\lambda^\alpha}{\Gamma(\alpha)}x^{\alpha}\exp(-\lambda x)\ dx \\
&=\frac{\Gamma(\alpha+1)}{\lambda\Gamma(\alpha)}\int_0^\infty \frac{\lambda^{\alpha+1}}{\Gamma(\alpha+1)}x^{\alpha}\exp(-\lambda x)\ dx \\
&=\frac{\alpha}{\lambda} \\
E(X^2)
&=\int_0^\infty \frac{\lambda^\alpha}{\Gamma(\alpha)}x^{\alpha+1}\exp(-\lambda x)\ dx \\
&=\frac{\Gamma(\alpha+2)}{\lambda^2\Gamma(\alpha)}\int_0^\infty \frac{\lambda^{\alpha+2}}{\Gamma(\alpha+2)}x^{\alpha}\exp(-\lambda x)\ dx \\
&=\frac{\alpha(\alpha+1)}{\lambda^2}\\
V(X)&=E(X^2)-E(X)^2\\
&=\frac{\alpha}{\lambda^2}
\end{align}
指数分布との本質的な違いは、「当たりが何個出るまでの時間を見るか」のみです。ガンマ分布は指数分布から作ることができます。
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$とすることで二項分布は正規分布に従います。正規分布の確率密度関数は下式で与えられます。
f(x;\mu,\sigma^2)=\frac{1}{\sqrt{2\pi\sigma^2}}\exp\biggl\{-\frac{(x-\mu)^2}{2\sigma^2}\biggr\}
正規分布の平均は$\mu$、分散は$\sigma^2$です。
[確率密度関数からの導出]
\begin{align}
E(X)&=\int_{-\infty}^{\infty}\frac{x}{\sqrt{2\pi\sigma^2}}\exp\biggl\{-\frac{(x-\mu)^2}{2\sigma^2}\biggr\}\ dx \\
&=\int_{-\infty}^{\infty}\frac{y+\mu}{\sqrt{2\pi\sigma^2}}\exp\biggl\{-\frac{y^2}{2\sigma^2}\biggr\}\ dy\\
&=\int_{-\infty}^{\infty}\frac{y}{\sqrt{2\pi\sigma^2}}\exp\biggl\{-\frac{y^2}{2\sigma^2}\biggr\}\ dy + \mu\\
&=\mu \\
V(X)&=\int_{-\infty}^{\infty}\frac{(x-\mu)^2}{\sqrt{2\pi\sigma^2}}\exp\biggl\{-\frac{(x-\mu)^2}{2\sigma^2}\biggr\}\ dx \\
&=\int_{-\infty}^{\infty}\frac{y^2}{\sqrt{2\pi\sigma^2}}\exp\biggl\{-\frac{y^2}{2\sigma^2}\biggr\}\ dx \\
&=-\sigma^2 \int_{-\infty}^{\infty}\frac{y}{\sqrt{2\pi\sigma^2}}\Biggl[\exp\biggl\{-\frac{y^2}{2\sigma^2}\biggr\}\Biggr]^{'}\ dx \\
&=\sigma^2 \int_{-\infty}^{\infty}\frac{1}{\sqrt{2\pi\sigma^2}}\exp\biggl\{-\frac{y^2}{2\sigma^2}\biggr\}\ dx \\
&=\sigma^2
\end{align}
平均$\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_{k}$が、独立に標準正規分布（平均0、分散1の正規分布）に従うとき、$X_1^2+X_2^2+\cdots+X_{k}^2$を自由度$k$の$\chi^2$分布と言います。
カイ二乗分布の確率密度関数は下式で与えられます。
f(x;k)=\frac{1}{2^{k/2}\Gamma(k/2)}x^{k/2-1}e^{-x/2}
カイ二乗分布の平均は$k$、分散は$2k$です。
[確率密度関数からの導出]
\begin{align}
E(X)&=\int_{0}^{\infty}\frac{x}{2^{k/2}\Gamma(k/2)}x^{k/2 \ -1}e^{-x/2} \ dx\\
&=\frac{2\Gamma(k/2 \ +1)}{\Gamma(k/2)}\int_{0}^{\infty}\frac{1}{2^{k/2 \ +1}\Gamma(k/2 \ +1)}x^{k/2}e^{-x/2} \ dx\\
&=2\times\frac{k}{2}=k\\
E(X^2)&=\int_{0}^{\infty}\frac{x^2}{2^{k/2}\Gamma(k/2)}x^{k/2 \ -1}e^{-x/2} \ dx\\
&=\frac{2^2 \ \Gamma(k/2 \ +2)}{\Gamma(k/2)}\int_{0}^{\infty}\frac{1}{2^{k/2 \ +2}\Gamma(k/2 \ +2)}x^{k/2 \ +1}e^{-x/2} \ dx\\
&=2^2\times\Bigl(\frac{k}{2}+1\Bigr)\frac{k}{2}=k(k+2)\\
V(X)&=E(X^2)-E(X)^2=2k
\end{align}
カイ二乗分布に従うサンプルは、次の方法で生成することができます。
def chi2(k, sampleSize):
x = norm(0, 1, sampleSize * k) ** 2
x = x.reshape([sampleSize, k])
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)=\frac{x^{\alpha-1}(1-x)^{\beta-1}}{B(\alpha, \beta)}\ \ (0\le x\le 1)
ただし、$B(\alpha, \beta)$は下式で与えられます。
B(\alpha, \beta)=\int_0^1 t^{\alpha-1}(1-t)^{\beta-1} dt
ベータ分布の平均は$\frac{\alpha}{\alpha+\beta}$、分散は$\frac{\alpha\beta}{(\alpha+\beta)^2(\alpha+\beta+1)}$です。
[確率密度関数からの導出]
\begin{align}
E(X)
&=\int_0^1 \frac{x^{\alpha}(1-x)^{\beta-1}}{B(\alpha, \beta)} dx\\
&=\int_0^1 \frac{x^{(\alpha + 1)-1}(1-x)^{\beta-1}}{B(\alpha + 1, \beta)} dx \ \ \frac{{B(\alpha + 1, \beta)}}{{B(\alpha , \beta)}}\\
&=\frac{{B(\alpha + 1, \beta)}}{{B(\alpha , \beta)}}\\
&=\frac{\alpha}{\alpha+\beta}\\
E(X^2)
&=\int_0^1 \frac{x^{\alpha+1}(1-x)^{\beta-1}}{B(\alpha, \beta)} dx\\
&=\int_0^1 \frac{x^{(\alpha + 2)-1}(1-x)^{\beta-1}}{B(\alpha + 2, \beta)} dx \ \ \frac{{B(\alpha + 2, \beta)}}{{B(\alpha , \beta)}}\\
&=\frac{{B(\alpha + 2, \beta)}}{{B(\alpha , \beta)}}\\
&=\frac{\alpha+1}{\alpha+\beta+1}\frac{{B(\alpha + 1, \beta)}}{{B(\alpha , \beta)}}\\
&=\frac{\alpha(\alpha+1)}{(\alpha+\beta)(\alpha+\beta+1)}\\
V(X)&=E(X^2)-E(X)^2\\
&=\frac{\alpha\beta}{(\alpha+\beta)^2(\alpha+\beta+1)}
\end{align}
def beta(alpha, beta, sampleSize):
x1 = gamma(alpha, 1, sampleSize)
x2 = gamma(beta, 1, sampleSize)
return x1 / (x1 + x2)
ベータ分布から一様分布
区間$[a,b]$の一様分布の密度関数は下式で与えられます。
f(x;a,b)=\frac{1}{b-a}\ \ (b\le x\le a)
区間 $[a,b]$の一様分布の平均は$\frac{a+b}{2}$、分散は$\frac{(b-a)^2}{12}$です。
[確率関数からの導出]
\begin{align}
E(X) &= \int_a^b \frac{x}{b-a} dx = \frac{a+b}{2}\\
E(X^2) &= \int_a^b \frac{x^2}{b-a} dx = \frac{a^2+ab+b^2}{3}\\
V(X) &= E[X^2] - E[X]^2 = \frac{(b-a)^2}{12}
\end{align}
ベータ分布の密度関数から明らかなとおり、$\alpha,\beta＝1$の場合、ベータ分布は区間[0,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)
- 結果
上手くいったことがわかります。