[data analysis] 감마분포(Gamma Distribution)

감마분포(Gamma Distribution)

확률은 전체 경우의 수 중에서 특정 시건의 비율이므로 확률을 계산하는데 있어 전체 경우의 수를 산정하는 것이 중요합니다. 이산변수의 경우 전체 경우의 수는 계승(factorial)을 사용하여 계산합니다. 연속변수의 경우에서 랜덤변수는 셀수 있는 수가 아니므로 계승을 직접적으로 사용할 수 없습니다. 대신에 계승에 대응하는 적분식을 사용하는데 이 부분을 감마함수(식 1)로 대체할 수 있습니다. 감마함수를 기본으로 하는 감마분포는 지수 분포와 정규분포와 관련된 분포입니다.

감마함수

Γ(x)로 나타내는 감마함수는 양의 정수 영역에서 factorial 함수의 형태를 가지며 식 1과 같이 정의합니다. 이 식의 첫번째는 이산변수, 두번째는 연속변수를 위한 함수입니다.

\begin{align}\tag{1} \Gamma(n)&=(n-1)! , \quad\quad n \in \{1,2,3, \cdots \}\\ &=(n-1)\Gamma(n-1)\\ &=\int^n_0 x^{n-1}e^{-x} ,\quad n > 0 \end{align}

(식 1)

예 1)

n = 10인 이산변수와 연속변수의 factorial과 감마함수를 계산합니다.

factorial은 numpy 또는 scipy의 math 모듈의 factorial()함수로 계산할 수 있으며 위의 감마함수는 scipy.special의 gamma() 함수로 계산할 수 있습니다. 또한 이 계산은 sympy 패키지의 symbol()함수와 intergrate()함수를 사용하여 수행할 수 있습니다.

import numpy as np
import pandas as pd 
import matplotlib.pyplot as plt
from scipy import stats
from sympy import *

np.math.factorial(10-1)

362880

from scipy import special
special.gamma(10)

362880.0

x=symbols('x')
g=x**9*exp(-x)
gam=g.integrate((x, 0, oo))
gam

362880

감마함수의 형태는 그림 1과 같습니다.

그림 1. 감마함수.

x=np.linspace(0.1, 4.1, 100)
y=scipy.special.gamma(x)
plt.figure(figsize=(4,3))
plt.plot(x, y)
plt.xlabel("x")
plt.ylabel(r"$\Gamma(x)$")
plt.show()

식 2와 같이 연속변수의 전체 표본공간을 감마함수(Γ(α))로 계산하는 확률밀도함수를 따르는 분포를 감마분포라고 합니다. 결과적으로 감마분포는 감마함수에 역비례인 분포로서 모수 α (전체발생수)와 λ(단위 시간당 평균발생수)에 의해 특정됩니다.

\begin{align}f(x) &=\frac{\lambda^\alpha x^{\alpha-1} \exp(-\lambda x)}{\Gamma(\alpha)}\\& \alpha,\, \lambda\, x \gt 0 \end{align}

(식 2)

이 분포를 따르는 연속랜덤변수는 식 3과 같이 나타냅니다.

X ~ Gamma(α, λ)	(식 3)

만약 α = 1인 경우 감마분포의 PDF는 식 4와 같이 지수 분포의 확률밀도함수와 같게 됩니다.

\begin{align}f(x)_{\text{Gamma} (1, λ)} &= λ \exp(-λx)\\&=f(x)_{\text{Exponential}(\lambda)}\end{align}

(식 4)

감마분포는 scipy.stats.gamma() 클래스를 사용하여 다양한 통계량을 계산할 수 있습니다. 예로서 감마분포의 확률밀도함수는 stats.gamma.pdf(x, a, loc=0, scale=1)를 사용하여 계산합니다. 이 메서드의 'a'는 α를 의미하고 분포의 이동과 규모는 매개변수 "loc"과 "scale"로 지정됩니다. 표준화된 데이터의 경우 "loc"과 "scale"은 0과 1이며 위 메소드의 기본값입니다. 분포의 규모를 나타내는 "scale"은 자료의 표준편차와 관련된 인수로서 표준화되지 않은 자료의 경우 $\text{scale} = \frac{1}{\lambda}$로 지정합니다. 즉, 감마함수의 모수인 λ는 분포의 scale에 영향을 주는 인자 입니다. 그림 2는 감마분포에서의 모수 α와 λ의 영향을 나타낸 것입니다.

그림 2. 감마분포의 모수 (a)α, (b)λ의 영향.

x=np.linspace(0, 20, 1000)
ab=[1, 3,5,7]
col=['g','b','r','k']
fig, (ax1, ax2)=plt.subplots(nrows=1, ncols=2, figsize=(8, 3))
for i, j in zip(col, ab):
    y=stats.gamma.pdf(x, j, 1)
    ax1.plot(x, y, color=i, label=f"Gamma({j}, 1)")
    y2=stats.gamma.pdf(x, 10, 0, 1/j)
    ax2.plot(x, y2, color=i, label=f"Gamma(10, {j})")
ax1.set_xlabel("x\n(a)", loc="center" , size="12")
ax1.set_ylabel("pdf")
ax1.legend(loc="best")
ax2.set_xlabel("x\n(b)", loc="center", size="12" )
ax2.legend(loc="best")
plt.show()

그림 2에서 나타낸 것과 같이 α가 1인 경우는 지수분포가 되며 α가 증가할수록 정규분포에 근사함을 알 수 있습니다. 이것은 샘플수가 증가하면 정규분포를 근접한다는 중심극한정리에 부합합니다. 또한 λ의 증가로 데이터의 퍼짐이 넓어지는 것으로 분포의 형태를 결정하는 인자가 됩니다.

감마분포의 확률밀도함수를 사용하여 계산한 평균과 분산은 각각 식 5와 같이 정의합니다.

\begin{align}E(X)&=\int^\infty_0 x\frac{\lambda^\alpha x^{\alpha -1}\exp(-\lambda x)}{\Gamma(\alpha)}\, dx\\ &=\int^\infty_0 \frac{\lambda^\alpha}{\Gamma(\alpha)} x^\alpha \exp(-\lambda x)\, dx\\ &= \int^\infty_0 \frac{\alpha \lambda^{\alpha+1}}{\Gamma(\alpha+1)\lambda}x^\alpha \exp(-\lambda x)\\ &=\frac{\alpha}{\lambda}\int^\infty_0 f(x)\, dx\\&=\frac{\alpha}{\lambda}\\\\ E(X^2)&=\int^\infty_0 x^2\frac{\lambda^\alpha x^{\alpha -1}\exp(-\lambda x)}{\Gamma(\alpha)}\, dx\\ &=\int^\infty_0 \frac{\lambda^\alpha}{\Gamma(\alpha)} x^{\alpha+1} \exp(-\lambda x)\, dx\\ &= \int^\infty_0 \frac{\alpha \lambda^{\alpha+2}}{\Gamma(\alpha+1)\lambda^2}x^{\alpha+1} \exp(-\lambda x)\\ &=\frac{\alpha(\alpha+1)}{\lambda^2}\int^\infty_0 f(x)\, dx\\&=\frac{\alpha(\alpha+1)}{\lambda^2}\\\\ \text{Var}(x)&=E(X^2) -(E(X))^2\\&=\frac{\alpha(\alpha+1)}{\lambda^2}-\left(\frac{\alpha}{\lambda}\right)^2\\&=\frac{\alpha}{\lambda^2}\end{align}

(식 5)

a, l, x=symbols("alpha, lambda, x", positive=True)
f=l**a*x**(a-1)*exp(-l*x)/gamma(a); f

$\frac{\lambda^{\alpha} x^{\alpha - 1} e^{- \lambda x}}{\Gamma\left(\alpha\right)}$

E=simplify(integrate(x*f,(x, 0, oo)));E

$\frac{\alpha}{\lambda}$

E2=simplify(integrate(x**2*f,(x, 0, oo)));E2

$\frac{\alpha \left(\alpha + 1\right)}{\lambda^{2}}$

var=E2-E**2
simplify(var)

$\frac{\alpha}{\lambda^{2}}$

예 2)

시료 분석기기에 예비용을 포함하여 2개의 퓨즈를 포함하고 있습니다. 일반적으로 이 퓨즈는 최대 50시간이 유지되어야 합니다. 그러나 퓨즈 제조사로부터의 사양은 각 100시간에 1번 끈어진다고 보고하고 있습니다. 두 퓨즈가 50시간 동안 유지 되지 않을 확률?

정해진 기간동안 사건의 평균 발생수는 λ = 0.01 hour^-1 입니다. 또한 두 번째 휴즈가 고장 날 때까지 경과 시간을 랜덤변수를 X라고하면 전체 발생수인 α = 2와 λ를 모수로하는 감마분포(X ~ Gamma (2, 0.01))로 나타낼 수 있습니다.

이 예는 P(X ≤ 50)을 계산하는 것과 같습니다.

round(stats.gamma.cdf(50, 2, scale=1/0.01),4)

0.0902