LoginSignup
7
11

More than 3 years have passed since last update.

@g-k

Pythonで学ぶ統計検定2級の確率分布②

Posted at

はじめに

統計検定の勉強をしていると様々な確率分布が出てきますが、数式だけ見ても中々イメージが付きにくいと思います。Pythonで色々パラーメータを動かしながら確率分布を描画してそのイメージをつけていきます。(前回の投稿はこちらです。前回は二項分布とポアソン分布を対象にしました。)

参考

確率分布の説明については下記を参考にしています。

様々な確率分布

本記事では各種数式の導出等の細かな説明は行わず、各分布の形とその分布が意味することのイメージを掴むことを主眼においています。この記事では下記3つの分布を扱います。

  • 幾何分布
  • 指数分布
  • 負の二項分布

幾何分布

「コインを投げたときに表が出るか裏が出るか」のように起こる結果が2つしかない互いに独立した試行(ベルヌーイ試行)が初めて成功するまでの回数$X$が従う分布を幾何分布と呼びます。二項分布が$n$回実施した中で成功する回数が従う分布なので非常に似ていますね。(二項分布については前回の記事をご参照ください。)

  • サイコロを振り続けて初めて1が出るまでの回数
  • コインを投げ続けて初めて表が出るまでの回数

などが幾何分布に従います。

幾何分布の確率質量関数の数式は下記のように表されます。


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

$p$はその試行の成功確率です。

また確率変数$X$が幾何分布に従う時、の期待値$E(X)$と分散$V(X)$は以下のようになります。


E(X) = \frac{1}{p}



V(X) = \frac{1-p}{p^2}

例えばサイコロを振り続けて初めて1が出るまでの試行回数の期待値は$\frac{1}{\frac{1}{6}} = 6$というイメージです。

それでは$p$の値(成功確率)が変化するにつれて幾何分布の形がどのように変化するかを描画していきます。


import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

def comb_(n, k):
    result = math.factorial(n) / (np.math.factorial(n - k) * np.math.factorial(k))
    return result

def geometric_dist(p, k):
    result = ((1 - p) ** (k - 1)) * p
    return result

fig = plt.figure()

def update(a):
    plt.cla()

    x =  np.arange(1, 50, 1)

    y = [geometric_dist(a,  i) for i in x]

    plt.bar(x, y, align="center", width=0.4, color="blue", 
                 alpha=0.5, label="binomial p= " + "{:.1f}".format(a))

    plt.legend()
    plt.ylim(0, 0.3)
    plt.xlim(0, 50)


ani = animation.FuncAnimation(fig,
                              update,
                              interval=1000,
                              frames = np.arange(0.1, 1, 0.1),
                              blit=True)
plt.show()
ani.save('Geometric_distribution.gif', writer='pillow')

Geometric_distribution.gif

こちらが成功確率$10$%の事象から$90$%の事象までの幾何分布です。成功確率が高くなるほど試行回数が少ない段階で成功する確率が高いのは感覚的にも理解できるかと思います。成功確率$90$%の事象で$10$回連続で失敗するのは殆ど奇跡であることがわかります。

また、幾何分布にはある事象が発生する確率は前に発生した結果の影響を受けないという性質があります。コインの表が$5$回連続で出たからといって次に裏が出る確率が高いわけではない、というイメージです。これは無記憶性と呼ばれています。(過去の情報を記憶しないという意味。)

指数分布

指数分布とは単位時間当たりに平均$\lambda$回起こる事象の発生間隔が$x$単位時間である確率を表す確率分布で、下記のような例に用いられます。

  • 災害が発生する間隔
  • 故障時間が一定であるシステムの偶発的な故障が発生する間隔
  • 店においてある客が来てから次の客が来るまでの間隔

指数分布の確率密度関数

指数分布の確率密度関数は下記のように表されます。

\begin{equation}
f(x)=
    \left\{
    \begin{aligned}
          &\lambda \mathrm{e}^{-\lambda x}　&(x\geq0) \\
          &0 &(x<0)\\
    \end{aligned}
    \right.
\end{equation}

また確率変数$X$が指数分布に従う時、の期待値$E(X)$と分散$V(X)$は以下のようになります。


E(X) = \frac{1}{\lambda}



V(X) = \frac{1}{\lambda}

例えば$1$時間(単位時間)あたり$5$回発生する事象の発生間隔の期待値は$\frac{1}{5}時間 = 12分$という感じです。何となくの感覚に合致していると思います。

それでは$\lambda$の値(単位時間あたりに発生する事象の平均回数)が変化させると指数分布の形がどのように変わるかを描画していきます。

例えばお客さんが1時間当たり平均30人来店する店から平均1人しか来店しない店まで、お客さんの来店間隔の時間を表している分布であると思って見てみましょう。

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

fig = plt.figure()

def update(a):
    plt.cla()

    x =  np.arange(0, 1, 0.01)

    y = [expon.pdf(i, scale = 1/a) for i in x]

    plt.plot(x, y, label="expon λ= %d" % a)
    plt.legend()
    plt.ylim(0, 35)
    plt.xlim(0, 1.0)

ani = animation.FuncAnimation(fig,
                              update,
                              interval=500,
                              frames = np.arange(31, 1, -1),
                              blit=True)
plt.show()
ani.save('Exponential_distribution.gif', writer='pillow')

Exponential_distribution.gif

分布の形状が何となく幾何分布と似ているのがわかるでしょうか。実は幾何分布の連続バージョンが指数分布です。（幾何分布において試行回数を増やしていき，そのぶん成功確率を$0$に近づけるような極限を適切に取ると指数分布になります。）

また、指数分布の形にはポイントが3つあります。

  • $\lambda$が小さいほど減少の仕方が緩やかになる
  • $\lambda$の値がどうであれ必ず単調減少する
  • $x$が$0$に近ければ近いほど確率密度は高くなる

ポイントの1つ目は「$1$時間で平均$5$人お客さんが来る店よりも、平均$15$人来る店の方が次のお客さんが間隔空けずにすぐに来るであろう」みたいなイメージなので感覚的にもわかりやすいと思います。

ポイントの2つ目と3つ目は「あるお客さんが来店した時、次のお客さん来店するのはその直後である確率が最も高い」ということです。これは何となく違和感のある方もいらっしゃるのではないでしょうか。実は指数分布の無記憶性に由来するもので、お客さんが$1$人来たからといって、次のお客さんは暫く来ないということも、連続で来店しやすいということもなく、お客さんの来店は完全にランダムであることを想定しています。従って発生しない時間がずっと続くよりもすぐに発生する確率の方が高いという考え方です。(この無記憶性は幾何分布の時に出てきたものと同じです。)

指数分布の累積密度関数

私たちの日常の感覚としてより知りたいのは、次のお客さんがきっかり$10$分後に来る確率、よりも$10$分以内に来る確率の方だと思います。その場合は次にお客さんが来店する間隔が$0$秒である確率から$10$分である確率までを全て足しあげる必要があります。

この時、確率密度関数を積分した結果の下記累積密度関数を使用します。
単位時間当たりに平均$\lambda$回起こる事象が発生するのが$x$単位時間以内である確率を表しています。


{{\begin{eqnarray}
F(x)  &=& 1 - \mathrm{e}^{-\lambda x} \\
\end{eqnarray}}
}

それでは$\lambda$の値(単位時間あたりに発生する事象の平均回数)が変化するにつれて累積密度関数のグラフの形がどのように変化するかを描画していきます。

先ほどと同様にお客さんが$1$時間当たり平均$30$人来店する店から平均$1$人しか来店しない店まで、次お客さんの来店するまでの時間が何分以内であるか、その確率を表している分布であると思って見てみましょう。


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

fig = plt.figure()

def update(a):
    plt.cla()

    x =  np.arange(0, 1.2, 0.01)

    y = [expon.cdf(i, scale = 1/a) for i in x]

    plt.plot(x, y, label="expon λ= %d" % a)

    plt.legend()
    plt.ylim(0, 1.2)
    plt.xlim(0, 1.2)

ani = animation.FuncAnimation(fig,
                              update,
                              interval=500,
                              frames = np.arange(30, 1, -1),
                              blit=True)
plt.show()
ani.save('Exponential_cumulative_distribution.gif', writer='pillow')

Exponential_cumulative_distribution.gif

$\lambda$が大きくなるにつれて分布のカーブが緩やかになっていくのがわかるかと思います。1時間に平均$30$人来店する店の方が$1$人しか来店しない店よりも$10$分以内に次のお客さんが来店する確率は高いと思うので納得感はあると思います。

負の二項分布

統計検定2級で出ることは殆どないですが今まで扱ってきた分布と関連性が深いので負の二項分布についても扱います。起こる結果が2つしかない互いに独立した試行(ベルヌーイ試行)が$r$回成功するまでに必要な試行回数$X$が従う分布を負の二項分布と呼びます。


P(X = k) = {}_{k-1} C _{r-1}p^r(1-p)^{k-r}
  • サイコロを振り続けて1が$3$回出るまでの試行回数
  • コインを投げ続けて表が$5$回出るまでの試行回数

などが負の二項分布に従います。

負の二項分布はその名前の通り二項分布の拡張版で下記のような違いがあります。

二項分布：試行回数を固定，成功回数が確率変数
負の二項分布：成功回数を固定，試行回数が確率変数

また$r=1$と置くと幾何分布と式になります。(幾何分布は初めて成功するまでの試行回数の確率分布であるため。)

更にポアソン分布の$\lambda$がガンマ分布に従う時に負の二項分布となります。(ガンマ分布については統計検定2級の範囲を超えるので本記事では扱いません。)

また確率変数$X$が負の二項分布に従う時、の期待値$E(X)$と分散$V(X)$は以下のようになります。


E(X) = \frac{r}{p}



V(X) = \frac{k(1-p)}{p^2}

例えばサイコロで1が5回出るまでに必要な試行回数の期待値は$\frac{5}{\frac{1}{6}} = 30$という感じです。

それではある事象が$10$回成功するまでに必要な試行回数の確率分布(負の二項分布)について成功確率を$10$%から$90$%の間で動かしながらみていきます。


import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

def comb_(n, k):
    result = math.factorial(n) / (np.math.factorial(n - k) * np.math.factorial(k))
    return result

def negative_binomial_dist(p, k, r):
    result = comb_(k - 1, r - 1) * (p**r) * ((1 - p) ** (k - r))
    return result

fig = plt.figure()

def update(a):
    plt.cla()

    r = 10

    x =  np.arange(r, 70, 1)

    y = [negative_binomial_dist(a, i, r) for i in x]

    plt.bar(x, y, align="center", width=0.4, color="blue", 
                 alpha=0.5, label="Negative binomial p= " + "{:.1f}".format(a))

    plt.legend()
    plt.ylim(0, 0.4)
    plt.xlim(10, 70)


ani = animation.FuncAnimation(fig,
                              update,
                              interval=1000,
                              frames = np.arange(0.1, 1, 0.1),
                              blit=True)
plt.show()
ani.save('Negative_binomial_distribution.gif', writer='pillow')

Negative_binomial_distribution.gif

成功確率が高ければ高いほど$10$回成功するまでの試行回数は$10$回に近づいていく様子がわかります。

確率分布間の関係性

統計検定2級の出題範囲となっている様々な確率分布を扱ってきましたが、実は全て関連しています。下記が今まで扱ってきた確率分布の相互の関係を図にしたものです。

分布間の関係性.png

数式だけをみているとよくわからないですが、実際に分布を描画しながら相互の確率分布の関係性を考えると理解が深まると思います。

NEXT

次回は正規分布やt分布について取り上げようと思います。

7
11
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
7
11