統計学入門⑤　確率分布と可視化 - Skywill CloudApps 技術ブログ

f:id:monozukuri-bu:20190520215822j:plain

こんにちは、モノづくり部のkanaiです。

統計学を勉強していると色々な確率分布が登場してくるかと思います。
今回はそんな確率分布についてご紹介していきたいと思います。

Python環境

以下の環境で動作確認しています。

Python 3.7
Jupyter 4.4.0
matplotlib 3.0.2
Pandas 0.23.4
seaborn 0.9.0
scipy  1.2.1

なお、Jupyter nobtebook上で実装する前提で進めていきます。

確率分布

起こりうる事象に割りあてている値を「確率変数」と言い、確率変数の各値がとる確率を表したものを「確率分布」と言います。

例えば、サイコロを投げて出る目は{1, 2, 3, 4, 5, 6}のいずれかになるので、これらの値が確率変数になります。
そしてそれぞれの目が出る確率は $\frac{1}{6}$ なので、確率分布と言えます。

つまり、サイコロの目をXとするとこんな感じの確率分布になります。

X	1	2	3	4	5	6
P(X)	$\frac{1}{6}$	$\frac{1}{6}$	$\frac{1}{6}$	$\frac{1}{6}$	$\frac{1}{6}$	$\frac{1}{6}$

また、確率分布には離散型と連続型があります。

サイコロの例では確率変数の値が1, 2, 3,...と離れていて、例えば1.73695...の目が出るということはありえません。
つまり確率変数が離散型となるので「離散型確率分布」であると言えます。
反対に、ある人の身長のように連続的な値に対する確率分布を「連続型確率分布」と言います。

離散型には、「二項分布」「ポアソン分布」「幾何分布」「離散一様分布」
連続型には、「正規分布」「連続一様分布」
といった確率分布があります。

scipy

Pythonではscipyのモジュールを使用することで確率分布を扱えるようになります。
scipyは高度な科学計算を行うためのライブラリです。
numpyでも科学計算は行えますが、scipyでは積分や統計など、より高度な計算を行うことができます。

scipyをインストールする前にまずnumpyをインストールする必要があります。
以下のコマンドでインストールしちゃいましょう。

pip install numpy
pip install scipy

インストールできたらjupyter notebookなどを起動しインポートしてみてください。
エラーが発生していなければ、scipyを使用できるようになっています。

import scipy

基本的な使い方としては

scipy.stats.[確率分布名].[メソッド]

の書き方で統一されているので、一度使い方を理解してしまえばいろいろな確率分布を扱えるようになるかと思います！

ベルヌーイ試行

コインを投げて表と裏のどちらが出るか。
このように何かをやったときに結果が2つにしかならない試行のことを「ベルヌーイ試行」といいます。

ベルヌーイ試行では、2つの結果のうち一方を「成功」とし、確率変数Xがとる値を「1」とします。
もう一方の結果を「失敗」とし、確率変数Xがとる値を「0」とします。

成功する確率を $p ( 0 \leq p \leq 1 )$ とすると、それぞれの確率は次のように表されます。

$成功する確率:\,P(X=1) = p \\\\ 失敗する確率:\,P(X=0) = 1-p$

二項分布

このベルヌーイ試行をn回行って、成功する回数Xが従う確率分布を「二項分布」といいます。
n回のベルヌーイ試行を行いk回成功する確率は次の式より計算できます。

$P(X=k)=_n\mathrm{C}_x{p}^{k}{(1-p)}^{(n-k)}$

例えば、表と裏が出る確率がそれぞれ0.5のコインを10回投げて5回表が出る確率は、0.2くらいになります。

$P(X=5)=_{10}\mathrm{C}_5*{0.5}^{5}*{(1-0.5)}^{(10-5)} = 0.205078125$

今度はこの試行を10回と言わず50回繰り返し、二項分布を可視化してみましょう。

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.stats import binom
%matplotlib inline

p = 0.5
N = 50
k = np.arange(N+1)

fig, ax = plt.subplots(1,1)
ax.plot(k, binom.pmf(k, N, p), 'bo', ms=8)
ax.vlines(k, 0, binom.pmf(k, N, p), colors='b', lw=1, alpha=0.2)
ax.set_ylim((0,0.15))
plt.show()

f:id:monozukuri-bu:20190520222427p:plain

二項分布を見ると、50回中表が出る回数が25回となる確率が一番高くなることがわかります。
また25回より多くなったり、少なくなるにつれて確率もだんだん低くなることも読み取れますね。

ポアソン分布

1年で馬に蹴られて死ぬ兵士の確率分布を知りたい。
これが歴史上初めてポアソン分布が使われた事例だそうです。

このように、ポアソン分布は「単位時間あたりに平均 λ 回起こる現象が、単位時間に k 回起きる確率」を扱う場合に使用されます。

例えば、以下のような事象はすべてポアソン分布に従うことが知られています。

・ある時間の範囲に会社にかかってくる電話の数
・ある人物が書いた文章の100文字中に存在するタイプミスの回数
・道路を歩いたときに発生する交通事故の回数

ポアソン分布は以下の式で表されます。

$P(k) = \frac{λ^k}{k!}e^{-λ}$

ポアソン分布は二項分布で λ=np を一定に保ちつつ、n(試行回数)を非常に大きく、p(確率)を非常に小さな値と近似することで得ることができます。

つまり、P(X=k)を二項分布とすると、次の数式が得られます。

$\lim_{n \to \infty} P(X=k) = \frac{λ^k}{k!}e^{-λ}$

それでは、1時間あたりに平均20回電話がかかってくる会社のポアソン分布を描画してみましょう。

from scipy.stats import poisson
import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inline

lamb = 20
k=np.arange(60)
pmf_pois = poisson.pmf(k,lamb)
plt.bar(k,pmf_pois)
plt.show()

f:id:monozukuri-bu:20190520222801p:plain

例えば1時間に22回電話がかかってくる確率を知りたいときは

poisson.pmf(22,lamb)

を実行してみると、7.7%ほどであることがわかります。

幾何分布

ベルヌーイ試行を繰り返し行い、初めて成功するまでの回数が従う分布を幾何分布と言います。
確率変数 X が幾何分布に従う場合、成功確率が p の試行において k 回目で初めて成功する確率は次の式で計算できます。

$P(X=k)=(1-p)^{k-1}p$

例えばサイコロを投げて5回目で初めて6が出る確率は次の計算より8%ほどになります。

$P(X=5)=(1-1/6)^(5-1) * 1/6 = \frac{625}{7776} = 0.0803...$

それでは、例として成功確率が０.1のベルヌーイ試行を繰り返し実施した時の幾何分布を描画してみましょう。

import numpy as np
from scipy.stats import geom
import matplotlib.pyplot as plt
%matplotlib inline

fig, ax = plt.subplots(1, 1)
p = 0.1
x = np.arange(geom.ppf(0.01, p), geom.ppf(0.99, p))
ax.plot(x, geom.pmf(x, p), 'bo', ms=8, label='geom pmf')
ax.vlines(x, 0, geom.pmf(x, p), colors='b', lw=5, alpha=0.5)

f:id:monozukuri-bu:20190520222920p:plain

正規分布

正規分布は最も代表的な分布の一つです。
物理の実験の測定の誤差やテストの点数などは、ほぼ正規分布に従うと考えられます。
統計的な手法においても、データが正規分布に従うと仮定します。

μを平均、σを分散とすると、正規分布の確率密度関数は次の式で表されます。

$f(x) = \frac{1}{\sqrt{2π}σ}\exp(-\frac{(x-μ)^2}{2σ^2})$

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
import matplotlib.pyplot as plt
 
x = np.arange(-10,10,0.01)
y = norm.pdf(x,0,1)
plt.plot(x,y,color='b')
plt.xlim(-5,5)
plt.show()

f:id:monozukuri-bu:20190520223037p:plain

分布からもわかるかと思いますが、正規分布の代表的な性質を以下に挙げてみます。

・平均値、最頻値、中央値が一致する。
・平均値を中心にして左右対称である。(直線x = μに関して対称)

一様分布

一様分布には離散一様分布と連続一様分布があります。
まず離散一様分布とは、確率変数が離散型である場合にすべての事象が起こる確率が等しい分布のことです。

たとえば、さいころを投げて1から6までの目が出る確率はすべて $\frac{1}{6}$ となり等しいので、これは離散一様分布に従います。

また、ルーレットも離散一様分布に従います。
ルーレットは円盤を等間隔に区切って、転がしたボールがどの区間に落ちるか当てるゲームですね。
この場合も、どの区間にボールが落ちるかは同様に確からしいので、離散一様分布になります。

ではこの区間の数を100区間、1000区間、10000区間・・・とどんどん増やしてみるとどうなるでしょうか。
ボールは0度以上360度未満の角度の円盤のあらゆる場所に止まる可能性があるといえます。
この角度は連続値で、例えばボールを投げたら186.39364209・・・度の位置に止まるという現象も起こりえます。

このように確率変数Xが連続確率変数で、xの値に関わらず確率密度関数が常に一定である確率分布を連続一様分布といいます。

連続一様分布はコンピュータで疑似乱数を発生させる場合などに使われます。
離散一様分布の例として、1 ~ 6までの目が出るサイコロを投げた時の振る舞いを可視化してみます。

import matplotlib.pyplot as plt
%matplotlib inline

roll_options = [1, 2, 3, 4, 5, 6]
roll_prob = 1 / 6
plt.bar(roll_options, roll_prob)
plt.show()

f:id:monozukuri-bu:20190520223354p:plain

今度は numpy.random.uniformで 0 ~ 1 の範囲で一様分布に従うランダムを10,000個生成し、データを描画してみましょう。

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

np.random.seed()
N = 10000
x = np.random.uniform(0.0, 1.0, N)
nbins = 50
plt.hist(x, nbins, normed=True)
plt.show()

f:id:monozukuri-bu:20190520223415p:plain