[現代数理統計学の基礎 with python]3章：代表的な確率分布

はじめに

　このシリーズは久保川達也氏著「現代数理統計学の基礎」を簡潔に説明し、その内容をpythonで実装してみよう、というものです。実装にはGoogle Colaboratory(以下：Colab)を使いました。ご指摘があればコメント欄に書いて下さるとありがたいです。
　説明が必要だと思ったところにだけ触れて、アウトプットできればいいよねくらいのスタンスで書くので本の内容をちゃんと全部理解したい方にとっては向いていないかもしれないです。
　式番号や命題・定義のインデックスは本に合わせて書いた場合、この記事では番号がとびとびになることがありますので予めご了承ください。

３章の概観

　変数を与えると確率を求められる関数が確率分布でしたね。様々な種類の確率分布にはそれぞれ特徴があり、使いどころがあります。想定する確率分布を間違えると予想を誤ることになるので、それぞれの確率分布にどんな特徴があるかを知っておくことは重要です。
　確率分布の期待値や分散を前章の確率生成関数や積率生成関数、特性関数を用いて求めることができますが、覚えてしまってよいと思います。使ってるうちに覚えるかもしれません。章の最後にスタインの等式とスターリング公式に触れています。記事で紹介していない確率分布も、ググればたくさん出てきます。
　確率生成関数などを用いた命題の証明は、別の機会に「確率生成関数・積率生成関数・特性関数」の記事を書きますのでそのときに合わせて紹介したいと思います。

離散確率分布

$ $
２章で期待値と分散を扱いましたが、期待値と分散の関係には触れていませんでした。$E[X]=\mu$とし、$$V(X)=E[(X-\mu)^2]=E[X^2-2\mu X+\mu^2]=E[X^2]- (E[X])^2$$
$$=E[X(X-1)]+E[X]-(E[X])^2 $$
という関係が導かれます。この関係式は今後たくさん出てくるので覚えておいてください。
　本で紹介している離散確率分布（変数$X$が離散型の確率分布）は以下の通りです。
・離散一様分布
・２項分布
・ポアソン分布
・幾何分布
・負の二項分布
・超幾何分布
では掻い摘んで見ていきましょう。

２項分布

$ $　
　２項分布の前に、ベルヌーイ試行について説明しておきます。本の表現を引用させてもらいますと、

ベルヌーイ試行とは、$p$の確率で'成功',$1-p$の確率で'失敗'する実験のを行うことをいい、確率変数$X$は'成功'のとき$1$,'失敗'のとき$0$をとる。

　２項分布とは、このベルヌーイ試行を独立に（前の試行は次の試行に影響与えない）$n$回行ったとき、変数$X$を「'成功'の回数」とした分布です。$k$回成功,$n-k$回失敗する確率は以下の式で表されます（'成功'、'失敗'は分かりやすく言っているだけで'病気にかかる'、'かからない'など二項対立しているものであればいい）。
$$P(k)={}_nC_kp^k(1-p)^{n-k}, \ k=0,1,2,...,n $$
　なぜ${}_nC_k$があるのかというと、$n$回の試行は独立に行われるため、$n$回のうち好きに$k$回の'成功'を選んでよいからです。試行回数と確率がそれぞれ$n,p$である二項分布(binominal distribution)を$Bin(n,p)$で表します。この本では、最後までこういう表記法で書かれているので慣れるしかありません。

例として、コインを30回投げた時と1000回投げた時の表が出る回数の確率分布を描いてみます。

ポアソン分布

　ポアソン分布は二項分布の特殊な例で、「希に起こる現象」を「大量に観測（試行）」できるとき（例：１日に起こる交通事故の件数の分布）には、二項分布ではなくポアソン分布を用います。つまり、先ほどの二項分布で$$n \to \infty,\ p \to 0$$の極限をとるとポアソン分布に収束します。なおポアソン分布自体の確率分布の式もありますが、この記事では割愛します。$np= \lambda$としたとき、ポアソン分布は$Po(\lambda)$と表します。例えば$n=10,p=0.1$の場合$\lambda=1$（だいたい１０回に１回起こる）になります。

では、pythonで二項分布とポアソン分布を確認してみましょう。

%matplotlib inline
import matplotlib.pyplot as plt
from scipy.special import comb#組み合わせを計算する関数
import pandas as pd

#二項分布のグラフ描画
def Bin(n,p,x_min,x_max,np):
  prob = pd.Series([comb(float(n),k)*p**k*(1-p)**(float(n)-k) for k in range(0,n+1)]) #それぞれのkでの確率を計算
  plt.bar(prob.index,prob,label=np)#棒グラフ（ｙ値,ｘ値）
  plt.xlim(x_min,x_max)
  plt.legend()
  plt.show()

Bin(1000,0.5,0,30,"n=30,p=0.5")#コインを30回
Bin(10000,0.5,4500,5500,"n=1000,p=0.5")#コインを1000回
Bin(40000,0.00007,0,15,"n=40000,p=0.00007")#nを大きく、pを小さくしてみる

これを実行すると、以下の三つのグラフが得られます。

　どうでしょう、同じ関数ですけど少し歪んでポアソン分布のようなものが描けましたね。

残りの三つの離散確率分布もそれぞれ独特の考え方がありますが、離散確率変数$X$が何を表すのかをしっかり意識しておけば読み進められると思います。

#連続分布
　本で紹介している連続分布は以下の通りです
・一様分布
・正規分布
・ガンマ分布、カイ二乗分布
・指数分布、ハザード分布
・ベータ分布
では、ここも掻い摘んで見ていきましょう。

正規分布

$ $　
　正規分布は平均を中心として対象な形をしていて扱いやすく、最も重要な確率分布です。
　確率変数$X$が平均$\mu,$分散$\sigma^2$の正規分布に従うとき、$X$の確率密度関数は$$f_X(x|\mu,\sigma^2)=\frac{1}{\sqrt{2\pi}\sigma}\exp(-\frac{(x-\mu)^2}{2\sigma^2})$$で与えられ、この分布を$\mathcal{N}(\mu,\sigma^2)$で表します。これを標準化した$\mathcal{N}(0,1)$を標準正規分布といい、$$\phi(z)=\frac{1}{\sqrt{2\pi}}\exp(-\frac{z^2}{2})$$と書きます（前の記事にグラフがあります）。
　標準正規分布の累積分布関数(積分値=確率)を$$\Phi(z)=\int_{-\infty}^z \phi(t) dt$$で表し、後の章で出てくる仮説検定や信頼区間を扱うときに使われます。

ガンマ分布、カイ二乗分布

$　$
　ガンマ分布の特殊な場合としてカイ二乗分布がありますが、統計学において重要なのはカイ二乗分布の方です。後の章で出てきますが、カイ二乗分布は母分散の区間推定、適合度検定、独立性の検定などを行う際に使われます。カイ二乗分布は、ガンマ関数を使って表す式よりも4章、5章で出てくる性質のほうが重要なので、ここではカイ二乗分布の形だけ描いておきます。自由度$n$のカイ二乗分布を$\chi_n^2$で表します。自由度というのも、次章以降で説明した方が理解が深まりますので割愛します。

import matplotlib.pyplot as plt
import numpy as np
from scipy import stats

x1 = np.arange(0,15,0.1)
y1 = stats.chi2.pdf(x=x1,df=1)#df=degree of freedom(自由度)です
y2 = stats.chi2.pdf(x=x1,df=2)
y3 = stats.chi2.pdf(x=x1,df=3)
y4 = stats.chi2.pdf(x=x1,df=5)
y5 = stats.chi2.pdf(x=x1,df=10)
y6 = stats.chi2.pdf(x=x1,df=12)

plt.figure(figsize=(7,5))
plt.plot(x1,y1, label='n=1')
plt.plot(x1,y2, label='n=2')
plt.plot(x1,y3, label='n=3')
plt.plot(x1,y4, label='n=5')
plt.plot(x1,y5, label='n=10')
plt.plot(x1,y6, label='n=12')

plt.ylim(0,0.7); plt.xlim(0,15)
plt.legend()
plt.show()

これを実行すると、以下のようになります。

指数分布、ハザード分布

$ $　
　指数分布の確率密度関数は以下の式で与えられ,$Ex(\lambda)$と表されます。$$f_X(x|\lambda)=\lambda e^{- \lambda x},\ x>0$$
　指数分布、ハザード分布は生存時間や機械が故障するまでの期間などの分布として用いられ、確率変数$X$は時間・期間を示すことが多いです。
指数分布の期待値、分散は$$E[X]=\frac{1}{\lambda}=\theta$$$$V(X)=\frac{1}{\lambda^2}$$
　であり、$\theta=2$の指数分布は自由度$n=2$のカイ二乗分布と一致します(上図)。
　$P(X>s)$は時間$s$を超えて生存する確率を表しており、$P(X>s)=1-F_X(s)=e^{-\lambda s}$となります（１は全確率、二項目は時間$s$までの累積分布関数です）。時間$s$生存したという条件のもとでさらに時間$t$を超えて生存する条件付確率は、$$P(X\ge s+t|X\ge s)=...=\frac{e^{-\lambda (s+t)}}{e^{-\lambda s}}=e^{^\lambda t}=P(X\ge t)$$となり、時間$s$生存していたという条件に依存しないことが分かります(条件付確率を思い出して、導いてみてください)。この性質を無記憶性といい、離散確率分布である幾何分布でも成り立ちます。今までラッキーなことが一度もなかった方も、明日幸運が訪れるかもしれません。

ベータ分布

　ベータ分布は確率変数$X$が区間（０，１）上に値をとり、その確率密度関数が$$f_X(x|a,b)=\frac{1}{B(a,b)}x^{a-1}(1-x)^{b-1}$$
で与えられ、$Beta(a,b)$で表します。$B(a,b)$はベータ関数で、
$$B(a,b)=\int_{0}^1 x^{a-1}(1-x)^{b-1} dx$$
です。確率密度関数を積分すると１になるのが簡単に確認できますね。ベータ関数は、6章のベイズ法というところで出てきます。

　少ししか紹介しませんでしたが、3章は以上です。ありがとうございました。

参考文献

「現代数理統計学の基礎」久保川達也氏著