pyhtonで学ぶ統計1 離散型確率分布

メモ書きになりますのでご了承ください。

統計をpython使って勉強していこうと思います。

## 一様分布(descrete uniform distribution)
$N$を正の整数とします。離散型確率変数$X$が、**離散一様分布**に従うとき確率関数は以下のように表せます。
$$
P(X=x|N)=\frac{1}{N},\ x=1,2,\cdots,N
$$

#### 期待値
$$
E[X]=\frac{1}{N}\sum_{x=1}^Nx=\frac{1}{N}・\frac{N(N+1)}{2}=\frac{N+1}{2}
$$

#### 分散
$$
E[X^2]=\frac{1}{N}\sum_{x=1}^Nx^2=\frac{1}{N}・\frac{N(N+1)(2N+1)}{6}=\frac{(N+1)(2N+1)}{6}
$$
から、
$$
Var(X)=E[X^2]-\{E[X]\}^2=\frac{(N+1)(2N+1)}{6}-\frac{N+1}{2}=\frac{(N+1)(N-1)}{12}
$$
となります。

pythonで確認していきます。
本来でなれば一様分布にはscipy.statsにuniformという関数がありますが、連続分布になってしまうのでここではnumpyを使用します。
$N=10$として、一様分布でサンプリングを行い確率分布を作ってみます。

```python
from numpy.random import randint
import matplotlib.pyplot as plt
import numpy as np

samples = randint(1, 11, 1000)
x, count = np.unique(samples, return_counts=True)
p = count / 1000
plt.bar(x=x, height=p)
plt.hlines(0.1, 0, 11, linestyles='--')
plt.xlim(0,11)
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/e8a541d2-bcf9-f57b-de83-2075dc56f9e6.png)

どの数値においても$1/N=1/10=0.1$付近となっていることが分かります。
次に、平均・分散の計算を行います。まず、確率変数とその確率から計算します。

```python
Ex = np.sum(p*x)
Ex2 = np.sum(p*x**2)
Varx = np.sum(p*(x-Ex)**2)

print('E[X]:         ', Ex)
print('E[X2]:        ', Ex2)
print('Var[X]:       ', Varx)
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/af7bdd1c-2583-895a-97f2-670fb78195df.png)

上記の式から算出します。

```python
print('E[X]:         ', (10+1)/2)
print('E[X2]:        ', (10+1)*(2*10+1)/6)
print('Var[X]:       ', (10+1)*(2*10+1)/6 - ((10+1)/2)**2)
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/8dcac83b-cb7c-8267-fb43-5130e4e93aa4.png)

今後のために、確率変数とその確率から平均・分散の計算をする関数を定義しておきます。

```python
def calc_mean_var(x ,p):
    Ex = np.sum(p*x)
    Ex2 = np.sum(p*x**2)
    Varx = np.sum(p*(x-Ex)**2)

    print('E[X]:         ', Ex)
    print('Var[X]:       ', Varx)
```

## ベルヌーイ分布(Bernoulli distribution)
ベルヌーイ試行とは、$p$の確率で成功、$1-p$の確率で失敗する実験を行うことをいいます。
確率偏巣$X$は成功のとき1、失敗のときに0をとるものとすると、確率変数の従う確率関数は次のように書けます。

```math
P(X=x|p)=\left\{
\begin{array}\\
p&(x=1)\\
1-p&(x=0)
\end{array}
\right.
```
これを**ベルヌーイ分布**と呼び、$Ber(p)$で表します。
#### 期待値
$$
E[X]=1×p+0×(1-p)=p
$$
#### 分散
$$
Var(X)=(1-p)^2×p+(0-p)^2×(1-p)=p(1-p)
$$

scipy.statsのbernoulliを使用して、分布を確認します。

```python
from scipy.stats import bernoulli

samples = bernoulli.rvs(p=0.2, size=1000)

x, count = np.unique(samples, return_counts=True)
p = count / 1000
plt.bar(x=x, height=p)
plt.hlines(0.2, -0.5, 1.5, linestyles='--')
plt.xlim(-0.5,1.5)
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/9a90bfb1-24cf-9e6f-9368-dc4570fc9b2e.png)


```python
Ex = np.sum(p*x)
Varx = np.sum(p*(x-Ex)**2)

print('E[X]:         ', Ex)
print('Var[X]:       ', Varx)
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/2c8ab4be-24bc-4203-c8a2-f53bdfe94821.png)

統計値の確認をします。左から平均、分散、尖度そして歪度の順番で出力されます。

```python
bernoulli.stats(p=0.2, moments='mvsk')
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/73ec14ce-7af8-4a4e-1e92-dcc09ec9d04f.png)

確率分布はpmfで計算できます。

```python
bernoulli.pmf(k=np.array([0,1]), p=0.2)
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/fcb7f1b5-93fc-5914-3de2-1542139745ed.png)

## 2項分布(binomial distribution)
ベルヌーイ試行を独立に$n$回行ったときの成功の回数の分布が２項分布となります。
$n$回の試行のうち$k$回成功、$n-k$回失敗となる確率は、**2項係数**$_nC_k$を、

```math
_nC_k=\frac{n!}{k!(n-k)!}=\begin{pmatrix}n\\k\end{pmatrix}
```
と表したきに、

```math
P(Y=k|n,p)=\begin{pmatrix}n\\k\end{pmatrix}p^k(1-p)^{n-k},\ k=0,1,2,\cdots,n
```
と書けます。これを**2項分布**と呼び、$Bin(n,p)$で表します。

#### 期待値
```math
\begin{align}
E[X]&=\sum_{k=0}^n\frac{n!}{k!(n-k)!}p^k(1-p)^{n-k}\times k\\
&=n\sum_{k=1}^n\frac{(n-1)!}{(k-1)!\bigl((n-1)-(k-1)\bigr)!}p^k(1-p)^{n-k}\\
&=np\sum_{k=1}^n\frac{(n-1)!}{(k-1)!\bigl((n-1)-(k-1)\bigr)!}p^{(k-1)}(1-p)^{(n-1)-(k-1)}\\
&=np\sum_{k=0}^{n-1}\frac{(n-1)!}{k!\bigl((n-1)-k\bigr)!}p^{k}(1-p)^{(n-1)-k}\\
&=np
\end{align}
```
#### 分散

```math
\begin{align}
E[X^2]&=\sum_{k=0}^n\frac{n!}{k!(n-k)!}p^k(1-p)^{n-k}\times k^2\\
&=\sum_{k=0}^nk(k-1)\frac{n!}{k!(n-k)!}p^k(1-p)^{n-k}+\sum_{k=0}^nk\frac{n!}{k!(n-k)!}p^k(1-p)^{n-k}\\
&=n(n-1)p^2\sum_{k=2}^{n}\frac{(n-2)!}{(k-2)!\bigl((n-2)-(k-2)\bigr)!}p^{k-2}(1-p)^{(n-2)-(k-2)}+np\\
&=n(n-1)p^2+np
\end{align}
```
であるから、

```math
\begin{align}
Var(X)&=E[X^2]-\{E[X]\}^2\\
&=n(n-1)p^2+np-(np)^2\\
&=np(1-p)
\end{align}
```
#### モーメント母関数

```math
\begin{align}
M_X(t)&=E[e^{tX}]\\
&=\sum_{k=0}^n\frac{n!}{k!(n-k)!}p^k(1-p)^{n-k}\times e^{tk}\\
&=\sum_{k=0}^n\frac{n!}{k!(n-k)!}(pe^t)^k(1-p)^{n-k}\\
&=(pe^t+1-p)^n
\end{align}
```
モーメント母関数の1次微分と2次微分はそれぞれ

```math
M'_X(t)=npe^{t}(pe^t+1-p)^{n-1}\\
M''_X(t)=np\bigl\{e^{t}(pe^t+1-p) + e^{t}(n-1)pe^{t} \bigr\}(pe^t+1-p)^{n-2}
```
と計算できます。$t=0$とすると、

```math
E[X]=M'_X(0)=np(p+1-p)^{n-1}=np\\
M''_X(0)=np\times (1 + p(n-1))\\
Var(X)=M''_X(0)-\{M'_X(0)\}^2=n(n-1)p^2
```
と求められる。

成功の確率が0.2であるベルヌーイ試行を10回行ったときの成功の回数の分布を確認してみます。

```python
samples1 = [np.sum(bernoulli.rvs(p=0.2, size=10)) for _ in range(10000)]
samples2 = [np.sum(bernoulli.rvs(p=0.6, size=10)) for _ in range(10000)]


x1, count1 = np.unique(samples1, return_counts=True)
x2, count2 = np.unique(samples2, return_counts=True)

p1 = count1 / 1000
p2 = count2 / 1000

plt.bar(x=x1, height=p1, alpha=.5, label='p=0.2')
plt.bar(x=x2, height=p2, alpha=.5, label='p=0.6')
plt.xlim(-1,11)
plt.legend()
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/c0704b1e-d184-c51a-76ec-6eaafec48f43.png)


もちろん専用の関数が用意されています。scipy.statsのbinomialを使ってサンプリングをして分布を出してみます。
また、binom.pmfで確率変数の値を引数に与えれば、そのときの確率を出力してくれます。
これを使って算出した確率分布も重ねて表示します。

```python
from scipy.stats import binom

samples1 = binom.rvs(p=0.2, n=10, size=1000)
samples2 = binom.rvs(p=0.6, n=10, size=1000)

x1, count1 = np.unique(samples1, return_counts=True)
x2, count2 = np.unique(samples2, return_counts=True)

p1 = count1 / 1000
p2 = count2 / 1000


x = np.arange(0, 10)
f1 = binom.pmf(k=x, n=10, p=0.2)
f2 = binom.pmf(k=x, n=10, p=0.6)

plt.bar(x=x1, height=p1, alpha=.5, label='p=0.2')
plt.bar(x=x2, height=p2, alpha=.5, label='p=0.6')
plt.plot(x, f1, '-o', alpha=.8, label='p=0.2')
plt.plot(x, f2, '-o', alpha=.8, label='p=0.6')


plt.xlim(-1,11)
plt.legend()
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/1de935e0-2529-d7ea-e381-e0c516a2c4b2.png)

サンプリングして得られた確率分布から期待値と分散を計算してみます。

```python
print('p=0.2')
calc_mean_var(x1 ,p1)
print('p=0.6')
calc_mean_var(x2 ,p2)
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/e5f3b1f4-9ade-c60d-2f7f-9844282385b5.png)

上で説明した計算結果に値を代入してみます。

```python
print('p=0.2')
print('np=', 10*0.2)
print('np(1-p)=', 10*0.2*(1-0.2))

print('p=0.6')
print('np=', 10*0.6)
print('np(1-p)=', 10*0.6*(1-0.6))
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/e09542be-4cdc-f248-230e-ddcfe28c54b4.png)

計算しなくてもbinom.statsを使って確認できます。

```python
print('p=0.2:  ',binom.stats(n=10, p=0.2))
print('p=0.6:  ',binom.stats(n=10, p=0.6))
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/f48c6701-f817-3789-97fb-c4e5d95d8cdf.png)


## ポアソン分布(poisson distribution)
稀な現象の大量観測によって発生する現象の個数の分布を表すときポアソン分布が用いられます。
この現象が起こる個数を$X$で表すと、$X$の確率関数は次のように表されます。

```math
P(X=x|\lambda)=\frac{\lambda^k}{k!}e^{-\lambda},\ k=1,2,\cdots
```

この確率分布を**ポアソン分布**と呼び、$Po(\lambda)$と表します。
$\lambda>0$は強度と呼ばれるパラメータであり、回数の平均を表します。
#### 期待値

```math
\begin{align}
E[X]&=\sum_{k=0}^n\frac{\lambda^k}{k!}e^{-\lambda}\times k\\
&=\lambda\sum_{k=1}^n\frac{\lambda^{k-1}}{(k-1)!}e^{-\lambda}\\
&=\lambda
\end{align}
```

#### 分散

```math
\begin{align}
E[X^2]&=\sum_{k=0}^n\frac{\lambda^k}{k!}e^{-\lambda}\times k^2\\
&=\sum_{k=1}^nk(k-1)\frac{\lambda^{k}}{k!}e^{-\lambda}+\sum_{k=0}^nk\frac{\lambda^{k}}{k!}e^{-\lambda}\\
&=\lambda^2\sum_{k=2}^n\frac{\lambda^{k-2}}{(k-2)!}e^{-\lambda}+\lambda\\
&=\lambda^2+\lambda
\end{align}
```
から、

```math
\begin{align}
Var(X)&=E[X^2]-\{E[X]\}^2\\
&=(\lambda^2+\lambda)-\lambda^2\\
&=\lambda
\end{align}
```

#### モーメント母関数

```math
\begin{align}
M_X(t)&=E[e^{tX}]\\
&=\sum_{k=0}^n\frac{\lambda^k}{k!}e^{-\lambda}\times e^{tk}\\
&=e^{-\lambda(1-e^t)}\sum_{k=0}^n\frac{(\lambda e^t)^k}{k!}e^{-\lambda e^t}\\
&=e^{-\lambda(1-e^t)}
\end{align}
```
モーメント母関数の1次微分と2次微分はそれぞれ

```math
M'_X(t)=\lambda e^t e^{-\lambda(1-e^t)}\\
M''_X(t)=\lambda e^t(1 + \lambda e^t)e^{-\lambda(1-e^t)}
```
と計算できます。$t=0$とすると、

```math
E[X]=M'_X(0)=\lambda\\
M''_X(0)=\lambda(1+\lambda)\\
Var(X)=M''_X(0)-\{M'_X(0)\}^2=\lambda
```

と求められる。

scipy.statsのpoissonを使って、サンプリングを行ってみます。

```python
from scipy.stats import poisson

samples1 = poisson.rvs(mu=2, size=1000)
samples2 =poisson.rvs(mu=5, size=1000)
samples3 =poisson.rvs(mu=10, size=1000)

x1, count1 = np.unique(samples1, return_counts=True)
x2, count2 = np.unique(samples2, return_counts=True)
x3, count3 = np.unique(samples3, return_counts=True)

p1 = count1 / 1000
p2 = count2 / 1000
p3 = count3 / 1000


x = np.arange(0, 20)
f1 = poisson.pmf(k=x, mu=2)
f2 = poisson.pmf(k=x, mu=5)
f3 = poisson.pmf(k=x, mu=10)

plt.bar(x=x1, height=p1, alpha=.5, label='mu=1')
plt.bar(x=x2, height=p2, alpha=.5, label='mu=5')
plt.bar(x=x3, height=p3, alpha=.5, label='mu=10')
plt.plot(x, f1, '-o', alpha=.8, label='p=2')
plt.plot(x, f2, '-o', alpha=.8, label='p=5')
plt.plot(x, f3, '-o', alpha=.8, label='p=10')

plt.xlim(-1,21)
plt.legend()
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/b7aca080-3059-1bdd-f1e8-647c1e36188c.png)

```python
print('mu=2')
calc_mean_var(x1 ,p1)
print('mu=5')
calc_mean_var(x2 ,p2)
print('mu=10')
calc_mean_var(x3 ,p3)
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/6c122966-50a4-ae6e-6eee-536f8af6b0e6.png)

```python
print('mu=2:  ', poisson.stats(mu=2))
print('mu=5:  ', poisson.stats(mu=5))
print('mu=10:  ', poisson.stats(mu=10))
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/1a820869-c5f9-991c-7717-20b0b5037dbb.png)

$np=\lambda$を一定とする条件のもとで$n→\infty,p→0$とすると、2項分布はポアソン分布に収束するという性質があります。
ここで$n=100000$として2項分布とポアソン分布の比較を行ってみます。

```python
x = np.arange(0, 20)
n = 100000

g1 = binom.pmf(k=x, n=n, p=2/n)
g2 = binom.pmf(k=x, n=n, p=5/n)
g3 = binom.pmf(k=x, n=n, p=10/n)

plt.bar(x, g1, alpha=.5, label='np=2(binom)')
plt.bar(x, g2, alpha=.5, label='np=5(binom)')
plt.bar(x, g3, alpha=.5, label='np=10(binom)')

plt.plot(x, f1, '-o', alpha=.5, label='p=2(poisson)')
plt.plot(x, f2, '-o', alpha=.5, label='p=5(poisson)')
plt.plot(x, f3, '-o', alpha=.5, label='p=10(poisson)')

plt.xlim(-1,21)
plt.legend()
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/f1c0cd43-f07f-960e-b1bc-283d7028ad1c.png)

2項分布がポアソン分布とほぼ同じ形状になることがわかりました。

## 幾何分布(geometric distribution)
成功確率$p$のベルヌーイ試行を独立に行っていき、初めて成功するまでに要した失敗の回数を$X$とするとき、$X$の分布が幾何分布となります。
$X=k$となる確率は、次のように書けます。

```math
P(X=k|p)=p(1-p)^k,/ k=0,1,2,\cdots
```
これを**幾何分布**と呼び、$Geo(p)$で表します。

#### 期待値

```math
\begin{align}
E[X]&=\sum_{k=0}^np(1-p)^k \times k
\end{align}
```

ここで、

```math
f(x)=\frac{1}{1-x}
```
のマクローリン展開を考えます。
マクローリン展開は、

```math
f(x)=\sum_{k=0}^\infty\frac{f^{(k)}(0)}{k!}x^k
```

と書けるので、まずはこのf$(x)$の微分を考えます。

```math
f'(0)=\frac{1}{(1-x)^2}|_{x=0}=1\\
f''(0)=\frac{2}{(1-x)^3}|_{x=0}=2\\
f'''(0)=\frac{2\times3}{(1-x)^3}|_{x=0}=2\times3
```

となることから、

```math
f^{(k)}(0)=k!
```

と考えられます。
従って、

```math
f(x)=\frac{1}{1-x}=\sum_{k=0}^\infty x^k
```

が成り立ちます($x=0$付近で)。
ここで、$x=1-p$と置くと、

```math
\frac{1}{p}=\sum_{k=0}^\infty (1-p)^k
```

さらに、両辺を$p$で微分すると(期待値の式に合わせるためです)、

```math
\frac{1}{p^2}=\sum_{k=0}^\infty k(1-p)^{k-1}
```

以上より、

```math
\begin{align}
E[X]&=\sum_{k=0}^np(1-p)^k \times k\\
&=p(1-p)\times\sum_{k=0}^n k(1-p)^{k-1}\\
&=\frac{1-p}{p}=\frac{q}{p}
\end{align}
```

#### 分散

```math
\begin{align}
E[X(X-1)]&=\sum_{k=0}^np(1-p)^k\times k(k-1)
\end{align}
```

期待値同様に、

```math
f(x)=\frac{1}{1-x}
```

のマクローリン展開を考えます。
先ほどの結果から、

```math
f(x)=\frac{1}{1-x}=\sum_{k=0}^\infty x^k
```

が成り立ちます。
$x=1-p$と置くと、

```math
\frac{1}{p}=\sum_{k=0}^\infty (1-p)^k
```

さらに、両辺を$p$で微分すると、

```math
\frac{1}{p^2}=\sum_{k=0}^\infty k(1-p)^{k-1}
```

そして、もう一度$p$で微分します。

```math
\frac{2}{p^3}=\sum_{k=0}^\infty k(k-1)(1-p)^{k-2}
```

これを$E[X(X-1)]$の式に代入します。

```math
\begin{align}
E[X(X-1)]&=\sum_{k=0}^np(1-p)^k\times k(k-1)\\
&= p(1-p)^2\sum_{k=0}^n k(k-1)(1-p)^{k-2}\\
&=2\frac{(1-p)^2}{p^2}
\end{align}
```

従って、

```math
\begin{align}
Var(X)&=E[X(X-1)]+E[X]-\{E[X]\}^2\\
&=2\frac{(1-p)^2}{p^2}+\frac{1-p}{p}-\frac{(1-p)^2}{p^2}\\
&=\frac{1-p}{p^2}
\end{align}
```

#### モーメント母関数

```math
\begin{align}
M_X(t)&=E[e^{tX}]\\
&=\sum_{k=0}^np(1-p)^k \times e^{tk}\\
&=\sum_{k=0}^np(e^t-pe^t)^k\\
&=p\frac{1}{1-(1-p)e^t}
\end{align}
```

モーメント母関数の1次微分と2次微分はそれぞれ

```math
M'_X(t)=\frac{p(1-p)e^t}{(1-(1-p)e^t)^2}\\
M''_X(t)=\frac{p(1-p)e^t(1+(1-p)e^t)}{(1-(1-p)e^t)^{3}}
```

と計算できます。$t=0$とすると、

```math
E[X]=M'_X(0)=\frac{1-p}{p}\\
M''_X(0)=\frac{(1-p)(2-p)}{p^2}\\
Var(X)=M''_X(0)-\{M'_X(0)\}^2=\frac{(2-3p+p^2)-(1-2p+p^2)}{p^2}=\frac{1-p}{p^2}
```

scipy.statsのgeomで分布などを確認します。

```python
from scipy.stats import geom

samples1 = geom.rvs(p = 0.2, size=1000)
samples2 = geom.rvs(p = 0.4, size=1000)
samples3 = geom.rvs(p = 0.8, size=1000)

x1, count1 = np.unique(samples1, return_counts=True)
x2, count2 = np.unique(samples2, return_counts=True)
x3, count3 = np.unique(samples3, return_counts=True)

p1 = count1 / 1000
p2 = count2 / 1000
p3 = count3 / 1000

x = np.arange(1, 20)
f1 = geom.pmf(k=x, p=0.2)
f2 = geom.pmf(k=x, p=0.4)
f3 = geom.pmf(k=x, p=0.8)

plt.bar(x=x1, height=p1, alpha=.5, label='p=0.2')
plt.bar(x=x2, height=p2, alpha=.5, label='p=0.4')
plt.bar(x=x3, height=p3, alpha=.5, label='p=0.8')
plt.plot(x, f1, '-o', alpha=.5, label='p=0.2')
plt.plot(x, f2, '-o', alpha=.5, label='p=0.4')
plt.plot(x, f3, '-o', alpha=.5, label='p=0.8')
plt.xlim(-1,21)
plt.legend()
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/7d79fe52-0eda-bd43-41dc-bc3741160d5f.png)

確率分布から平均と分散を求めてみます。

```python
print('p=0.2')
calc_mean_var(x1 ,p1)
print('p=0.4')
calc_mean_var(x2 ,p2)
print('p=0.8')
calc_mean_var(x3 ,p3)
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/64d31310-1fce-d2a5-e6e5-ec4a27d4f1c2.png)

関数を使って求めます。

```python
print('p=0.2:  ', geom.stats(p=0.2))
print('p=0.4:  ', geom.stats(p=0.4))
print('p=0.8:  ', geom.stats(p=0.8))
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/d43b56e0-8560-43f9-6bcb-a395280a37b4.png)

## 負の2項分布(negative binomial)
成功確率が$p$のベルヌーイ分布について、$r$回成功するまでに要した失敗の回数を$X$とするとき、$X$の分布が負の2項分布になります。
合計$r+k$の試行のうち、成功が$r$回、失敗が$k$回となるので、確率分布は次のようになります。

```math
P(X=x|r,p)=\begin{pmatrix}r+k-1\\k\end{pmatrix}p^rq^k,\ k=0,1,2,\cdots
```

これを**負の2項分布**と呼び、$NB(r,p)$と表します。
#### 期待値

```math
\begin{align}
E[X]&=\sum_{k=0}^\infty\begin{pmatrix}r+k-1\\k\end{pmatrix}p^rq^k \times k\\
&=\sum_{k=0}^\infty k\frac{(r+k-1)\times\cdots\times r}{k\times\cdots\times 1}p^rq^k\\
&=r\frac{q}{p}\sum_{k=1}^\infty\frac{(r+k-1)\times\cdots\times (r+1)}{k-1\times\cdots\times 1}p^{r+1}q^{k-1}\\
&=\frac{rq}{p}
\end{align}
```

#### 分散

```math
\begin{align}
E[X(X-1)]&=\sum_{k=0}^\infty\begin{pmatrix}r+k-1\\k\end{pmatrix}p^rq^k \times k(k-1)\\
&=\sum_{k=0}^\infty k(k-1)\frac{(r+k-1)\times\cdots\times r}{k\times\cdots\times 1}p^rq^k\\
&=r(r+1)\frac{q^2}{p^2}\sum_{k=1}^\infty\frac{(r+k-1)\times\cdots\times (r+2)}{k-2\times\cdots\times 1}p^{r+2}q^{k-2}\\
&=\frac{r(r+1)q^2}{p^2}
\end{align}
```

従って、

```math
\begin{align}
Var(X)&=E[X(X-1)]+E[X]-\{E[X]\}^2\\
&=\frac{r(r+1)q^2}{p^2}+\frac{rq}{p}-\frac{r^2q^2}{p^2}\\
&=\frac{rq}{p^2}
\end{align}
```

#### モーメント母関数

```math
\begin{align}
M_X(t)&=E[e^{tX}]\\
&=\sum_{k=0}^\infty\begin{pmatrix}r+k-1\\k\end{pmatrix}p^rq^k \times e^{tk}\\
&=\sum_{k=0}^\infty\begin{pmatrix}r+k-1\\k\end{pmatrix}p^r(qe^t)^k\\
&=\frac{p^r}{(1-qe^t)^r}\sum_{k=0}^\infty\begin{pmatrix}r+k-1\\k\end{pmatrix}(1-qe^t)^r(qe^t)^k\\
&=\frac{p^r}{(1-qe^t)^r}
\end{align}
```

モーメント母関数の1次微分と2次微分はそれぞれ

```math
M'_X(t)=rp^rq(1-qe^t)^{-r-1}(e^t)\\
M''_X(t)=rp^rq((-r-1)(1-qe^t)^{-r-2}(-qe^{2t}) + (1-qe^t)^{-r-1}(e^t))
```

と計算できます。$t=0$とすると、

```math
E[X]=M'_X(0)=\frac{rq}{p}\\
M''_X(0)=r\frac{q(1+rq)}{p^2}\\
Var(X)=M''_X(0)-\{M'_X(0)\}^2=\frac{rq}{p^2}
```

scipy.statsのnbinomで分布などを確認します。

```python
from scipy.stats import nbinom

samples1 = nbinom.rvs(p = 0.2, n=5, size=1000)
samples2 = nbinom.rvs(p = 0.4, n=5, size=1000)
samples3 = nbinom.rvs(p = 0.8, n=5, size=1000)

x1, count1 = np.unique(samples1, return_counts=True)
x2, count2 = np.unique(samples2, return_counts=True)
x3, count3 = np.unique(samples3, return_counts=True)

p1 = count1 / 1000
p2 = count2 / 1000
p3 = count3 / 1000

x = np.arange(0, 40)
f1 = nbinom.pmf(k=x, p=0.2, n=5)
f2 = nbinom.pmf(k=x, p=0.4, n=5)
f3 = nbinom.pmf(k=x, p=0.8, n=5)

plt.bar(x=x1, height=p1, alpha=.5, label='p=0.2')
plt.bar(x=x2, height=p2, alpha=.5, label='p=0.4')
plt.bar(x=x3, height=p3, alpha=.5, label='p=0.8')
plt.plot(x, f1, '-o', alpha=.5, label='p=0.2')
plt.plot(x, f2, '-o', alpha=.5, label='p=0.4')
plt.plot(x, f3, '-o', alpha=.5, label='p=0.8')

plt.xlim(-1,41)
plt.legend()
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/c1e96ad8-6cb6-33f3-f957-ed475d14437a.png)

```python
print('p=0.2 n=5')
calc_mean_var(x1 ,p1)
print('p=0.4 n=5')
calc_mean_var(x2 ,p2)
print('p=0.8 n=5')
calc_mean_var(x3 ,p3)
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/0abadcf5-4539-3d22-9f7b-e60453ec0804.png)

```python
print('p=0.2 n=5:  ', nbinom.stats(p=0.2, n=5))
print('p=0.4 n=5:  ', nbinom.stats(p=0.4, n=5))
print('p=0.8 n=5:  ', nbinom.stats(p=0.8, n=5))
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/034e0704-6962-c142-c7dd-3fcfa420ee25.png)

```python
samples1 = nbinom.rvs(p = 0.4, n=2, size=1000)
samples2 = nbinom.rvs(p = 0.4, n=5, size=1000)
samples3 = nbinom.rvs(p = 0.4, n=10, size=1000)

x1, count1 = np.unique(samples1, return_counts=True)
x2, count2 = np.unique(samples2, return_counts=True)
x3, count3 = np.unique(samples3, return_counts=True)

p1 = count1 / 1000
p2 = count2 / 1000
p3 = count3 / 1000

x = np.arange(0, 40)
f1 = nbinom.pmf(k=x, p=0.4, n=2)
f2 = nbinom.pmf(k=x, p=0.4, n=5)
f3 = nbinom.pmf(k=x, p=0.4, n=10)

plt.bar(x=x1, height=p1, alpha=.5, label='n=2')
plt.bar(x=x2, height=p2, alpha=.5, label='n=5')
plt.bar(x=x3, height=p3, alpha=.5, label='n=10')
plt.plot(x, f1, '-o', alpha=.5, label='n=2')
plt.plot(x, f2, '-o', alpha=.5, label='n=5')
plt.plot(x, f3, '-o', alpha=.5, label='n=10')
plt.xlim(-1,41)
plt.legend()
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/b6da785f-1bce-80bc-bc75-954fe6a09bf4.png)


```python
print('p=0.4 n=2')
calc_mean_var(x1 ,p1)
print('p=0.4 n=5')
calc_mean_var(x2 ,p2)
print('p=0.4 n=10')
calc_mean_var(x3 ,p3)
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/43b2834b-d4a9-4e4e-baa3-c87a6b3761b4.png)

```python
print('p=0.4 n=2:   ', nbinom.stats(p=0.4, n=2))
print('p=0.4 n=5:   ', nbinom.stats(p=0.4, n=5))
print('p=0.4 n=10:  ', nbinom.stats(p=0.4, n=10))
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/d9249682-b4f5-3fca-b90f-264c943a2657.png)


## 超幾何分布(hypergeometric distribution)
$M$個の赤いボールと$N-M$個の白いボールが入っている壺の中から$K$個のボールを無作為に非復元中抽出で抽出したとき、赤いボールが$X$個であったとします。
このとき、$X$は次の確率分布に従います。

```python
P(X=x|N,M,K)=\frac{\begin{pmatrix}M\\x\end{pmatrix}\begin{pmatrix}N-M\\K-x\end{pmatrix}}{\begin{pmatrix}N\\K\end{pmatrix}},\ x=0,1,\cdots,k
```

この分布を**超幾何分布**と呼びます。

#### 期待値

```math
\begin{align}
E[X]&=\sum_{x=0}^\infty\frac{\begin{pmatrix}M\\x\end{pmatrix}\begin{pmatrix}N-M\\K-x\end{pmatrix}}{\begin{pmatrix}N\\K\end{pmatrix}} \times x\\
&=\sum_{x=1}^\infty x\frac{\frac{M!}{x!(M-x)!}\times\frac{(N-M)!}{(K-x)!(N-M-K+x)!}}{\frac{N!}{K!(N-K)!}}\\
&=\sum_{x=1}^\infty K\frac{M}{N}\frac{\frac{(M-1)!}{(x-1)!((M-1)-(x-1))!}\times\frac{(N-M)!}{(K-x)!(N-M-K+x)!}}{\frac{(N-1)!}{(K-1)!((N-1)-(K-1))!}}\\
&= K\frac{M}{N}\sum_{x=1}^\infty \frac{\frac{(M-1)!}{(x-1)!((M-1)-(x-1))!}\times\frac{((N-1)-(M-1))!}{((K-1)-(x-1))!((N-1)-(M-1)-(K-1)+(x-1))!}}{\frac{(N-1)!}{(K-1)!((N-1)-(K-1))!}}\\
&=Kp\ (p=M/N)
\end{align}
```

#### 分散

```math
\begin{align}
E[X(X-1)]&=\sum_{x=0}^\infty\frac{\begin{pmatrix}M\\x\end{pmatrix}\begin{pmatrix}N-M\\K-x\end{pmatrix}}{\begin{pmatrix}N\\K\end{pmatrix}} \times x(x-1)\\
&=\sum_{x=1}^\infty x(x-1)\frac{\frac{M!}{x!(M-x)!}\times\frac{(N-M)!}{(K-x)!(N-M-K+x)!}}{\frac{N!}{K!(N-K)!}}\\
&=K(K-1)\frac{M(M-1)}{N(N-1)}
\end{align}
```

従って、

```math
\begin{align}
Var(X)&=E[X(X-1)]+E[X]-\{E[X]\}^2\\
&=K(K-1)\frac{M(M-1)}{N(N-1)}+K\frac{M}{N}-K^2\frac{M^2}{N^2}\\
&=K(K-1)\frac{p(M-1)}{N-1}+Kp-K^2p^2\\
&=\frac{Kp}{N-1}(-K-Np+N+Kp)\\
&=\frac{N-K}{N-1}Kp(1-p)
\end{align}
```

scipy.statsのhypergeomを使って分布の確認していきます。

```python
from scipy.stats import hypergeom

samples1 = hypergeom.rvs(M=100, n=10, N=10, size=1000)
samples2 = hypergeom.rvs(M=100, n=10, N=50, size=1000)
samples3 = hypergeom.rvs(M=100, n=10, N=90, size=1000)

x1, count1 = np.unique(samples1, return_counts=True)
x2, count2 = np.unique(samples2, return_counts=True)
x3, count3 = np.unique(samples3, return_counts=True)

p1 = count1 / 1000
p2 = count2 / 1000
p3 = count3 / 1000

x = np.arange(0, 40)
f1 = hypergeom.pmf(k=x, M=100, n=10, N=10)
f2 = hypergeom.pmf(k=x, M=100, n=10, N=50)
f3 = hypergeom.pmf(k=x, M=100, n=10, N=90)

plt.bar(x=x1, height=p1, alpha=.5, label='M=100 n=10 N=10')
plt.bar(x=x2, height=p2, alpha=.5, label='M=100 n=10 N=50')
plt.bar(x=x3, height=p3, alpha=.5, label='M=100 n=10 N=90')
plt.plot(x, f1, '-o', alpha=.5, label='M=100 n=10 N=10')
plt.plot(x, f2, '-o', alpha=.5, label='M=100 n=10 N=50')
plt.plot(x, f3, '-o', alpha=.5, label='M=100 n=10 N=90')

plt.xlim(-1,11)
plt.ylim(0,0.6)
plt.legend(loc='upper center')
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/232c5e0e-6dbe-0569-9a1e-f2a26172540b.png)


```python
print('M=100 n=10 N=10')
calc_mean_var(x1 ,p1)
print('M=100 n=10 N=50')
calc_mean_var(x2 ,p2)
print('M=100 n=10 N=90')
calc_mean_var(x3 ,p3)
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/1c3b2d56-e0d0-0f31-5614-46506c898de8.png)

```python
print('M=100 n=10 N=10:   ', hypergeom.stats(M=100, n=10, N=10))
print('M=100 n=10 N=50:   ', hypergeom.stats(M=100, n=10, N=50))
print('M=100 n=10 N=90:  ', hypergeom.stats(M=100, n=10, N=90))
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/7016c0af-4d05-b5d4-b3b5-fac8ca8b421d.png)


```python
samples1 = hypergeom.rvs(M=100, n=20, N=20, size=1000)
samples2 = hypergeom.rvs(M=100, n=50, N=20, size=1000)
samples3 = hypergeom.rvs(M=100, n=80, N=20, size=1000)

x1, count1 = np.unique(samples1, return_counts=True)
x2, count2 = np.unique(samples2, return_counts=True)
x3, count3 = np.unique(samples3, return_counts=True)

p1 = count1 / 1000
p2 = count2 / 1000
p3 = count3 / 1000

x = np.arange(0, 40)
f1 = hypergeom.pmf(k=x, M=100, n=20, N=20)
f2 = hypergeom.pmf(k=x, M=100, n=50, N=20)
f3 = hypergeom.pmf(k=x, M=100, n=80, N=20)

plt.bar(x=x1, height=p1, alpha=.5, label='M=100 n=20 N=20')
plt.bar(x=x2, height=p2, alpha=.5, label='M=100 n=50 N=20')
plt.bar(x=x3, height=p3, alpha=.5, label='M=100 n=80 N=20')
plt.plot(x, f1, '-o', alpha=.5, label='M=100 n=20 N=20')
plt.plot(x, f2, '-o', alpha=.5, label='M=100 n=50 N=20')
plt.plot(x, f3, '-o', alpha=.5, label='M=100 n=80 N=20')

plt.xlim(-1,21)
plt.ylim(0,0.6)
plt.legend(loc='upper center')
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/c6613599-884c-ee47-2791-b4bf1821c669.png)


```python
print('M=100 n=20 N=20')
calc_mean_var(x1 ,p1)
print('M=100 n=50 N=20')
calc_mean_var(x2 ,p2)
print('M=100 n=80 N=20')
calc_mean_var(x3 ,p3)
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/93980397-ac67-b581-a68f-3028ec023d6b.png)


```python
print('M=100 n=20 N=20:   ', hypergeom.stats(M=100, n=20, N=20))
print('M=100 n=50 N=20:   ', hypergeom.stats(M=100, n=50, N=20))
print('M=100 n=80 N=20:  ', hypergeom.stats(M=100, n=80, N=20))
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/3f603532-bf93-c039-d1a5-f3d43d324276.png)

### 参考

https://www.amazon.co.jp/%E7%8F%BE%E4%BB%A3%E6%95%B0%E7%90%86%E7%B5%B1%E8%A8%88%E5%AD%A6%E3%81%AE%E5%9F%BA%E7%A4%8E-%E5%85%B1%E7%AB%8B%E8%AC%9B%E5%BA%A7-%E6%95%B0%E5%AD%A6%E3%81%AE%E9%AD%85%E5%8A%9B-%E4%B9%85%E4%BF%9D%E5%B7%9D-%E9%81%94%E4%B9%9F/dp/4320111664/ref=pd_bxgy_img_1/357-1227394-8103025?pd_rd_w=SNde5&pf_rd_p=d8f6e0ab-48ef-4eca-99d5-60d97e927468&pf_rd_r=5VV11NNQY83GR759QRKF&pd_rd_r=1505e9a1-de31-4a38-9d1b-43e8b4d805c9&pd_rd_wg=IY0ds&pd_rd_i=4320111664&psc=1