[data analysis]지수분포(Exponential Distribution)

지수분포(Exponential Distribution)

지수분포는 가장 많이 사용되는 연속분포 중의 하나로 어떤 사건들 사이에 시간의 경과를 모형화 하는데 많이 사용됩니다. 이 분포의 확률밀도함수(PDF)는 식 1과 같이 정의되며 모수 λ는 단위 시간동안 평균 이벤트 발생 횟수를 의미합니다.

$$\begin{array}{rl}& X\sim \text{Exponential}(\lambda )\\ & f(x)=\{\begin{array}{ll}\lambda e^{-\lambda x}& x>0\\ 0& \text{otherwise}\end{array}\\ & \lambda >0\end{array}$$

(식 1)

sympy.integrate()함수를 적용하여 지수분포의 누적분포함수를 계산합니다.

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

x, l=symbols("x lambda", positive=True)
f=l*exp(-l*x)
f

$\lambda e^{-\lambda x}$

F=integrate(f, (x, 0, x))
F

$1-e^{-\lambda x}$

지수분포의 평균과 분산은 각각의 정의에 따라 식 2와 같이 정의됩니다.

$$\begin{array}{rl}& \begin{array}{rl}E(x)& =\mu \\ & ={\int }_0^{\mathrm{\infty }}xf(x)dx\\ & ={\int }_0^{\mathrm{\infty }}x\lambda e^{-\lambda x}dx\\ & =\frac{1}{\lambda }(-e^{-\lambda x}-x{e}^{-\lambda x}){|}_0^{\mathrm{\infty }}\\ & =\frac{1}{\lambda }\end{array}\\ & \begin{array}{rl}E(x^2)& ={\int }_0^{\mathrm{\infty }}{x}^{2}f(x)dx\\ & ={\int }_0^{\mathrm{\infty }}{x}^2\lambda e^{-\lambda x}dx\\ & =\frac{2}{{\lambda }^2}\end{array}\\ & \begin{array}{rl}\text{Var(X)}& =E(x^2)-(E(x){)}^2\\ & =\frac{2}{{\lambda }^2}-\frac{1}{{\lambda }^2}\\ & =\frac{1}{{\lambda }^2}\end{array}\end{array}$$

(식 2)

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

$\frac{1}{\lambda }$

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

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

var=E2-E**2
var

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

지수분포에 부합하는 랜덤변수의 통계량은 scipy.stats.expon() 클래스의 여러 메소드를 사용하여 계산할 수 있습니다. 이 클래스이 메소드의 인수 중 scale은 지수분포의 표준편차인 $\frac{1}{\lambda }$입니다. 즉, scale은 평균 빈도인 λ를 기반으로 하는 데이터의 퍼짐을 조절하는 모수(paramenter)라고 할 수 있습니다.

예 1)

어떤 기계를 수리하는데 걸리는 시간은 모수 λ = 3인 지수분포를 따른다고 합니다. 평균과 표준편차를 결정합니다.

그 통계량들은 지수분포의 확률밀도함수(식 1)를 사용하여 계산할 수 있습니다.

x, l=symbols("x lambda", positive=True)
f=l*exp(-l*x);f
E=integrate(x*f,(x, 0, oo))
E.subs(l, 3)

$\frac{1}{3}$

E2=integrate(x**2*f,(x, 0, oo)); E2
E2.subs(l, 3)

$\frac{2}{9}$

var=E2-E**2
var.subs(l, 3)

$\frac{1}{9}$

위 결과들을 stats.expon 클래스로 확인해 봅니다.

mu, var=stats.expon.stats(scale=1/3, moments="mv")
print(f"평균: {np.around(mu, 3)}, 분산: {np.around(var, 3)}")

평균: 0.333, 분산: 0.111

클래스의 moment() 메서드를 사용하여 평균과 분산을 계산해 봅니다.

E=stats.expon.moment(1,scale=1/3); round(E, 3)

0.333

E2=stats.expon.moment(2, scale=1/3); round(E2, 3)

0.222

var=E2-E**2; round(var, 3)

0.111

그림 1에서 나타낸 것과 같이 지수분포의 매개변수 λ는 분포의 모양을 결정합니다.

그림 1. λ에 따른 지수분포의 변화.

x=np.linspace(0.01, 6, 1000)
fig, ax=plt.subplots(figsize=(4,3))
col=["g", "b", "r"]
for i, j in enumerate([0.5, 1, 2]):
    y=stats.expon.pdf(x, j)
    idx=np.where(y>0)[0]
    ax.plot(x[idx], y[idx],  color=col[i], label=r"$\lambda$="+str(j))
ax.set_xlabel("x", loc="right")
ax.set_ylabel("pdf", loc="top")
ax.legend(loc="best", frameon=False)
plt.show()

위와 같은 지수분포의 형태는 이산변수의 기하분포와 유사합니다.

기하분포는 성공확률이 p인 베르누이 시행을 반복하여 처음으로 성공이 나올 확률을 나타내는 분포입니다. 이 분포가 연속변수의 분포인 지수분포와 유사하다는 것은 기하분포의 모수인 p와 지수분포의 모수인 λ사이에 연관성을 생각할 수 있습니다. 연속분포의 λ의 경우는 확률의 소멸 속도와 관계있는 모수로서 작을수록 확률의 변화가 작아집니다. 역시 기하분포의 p역시 작을수록 확률의 변화가 작습니다. 그러므로 지수분포는 어떠한 시행에서 처음으로 성공할 확률이 매우 작은 경우로 해석이 가능할 것입니다. 예를 들어 매우 짧은 시간 간격으로 고객이 상점을 방문하는 데이터에 대한 확률분포를 지수분포로 생각할 수 있을 것입니다.

이와같이 지수분포의 해석은 기하분포와의 특성들로 이해될 수 있습니다. 이러한 특성들 중에 중요한 점은 이전 조건에 어떠한 제약을 받지 않는다는 것입니다. 만약에 동전던지기를 반복 시행 할 경우 현재까지 앞면(성공)을 관찰하지 못한 경우 그 시점부터의 확률은 이전의 결과에 의해 제약을 받지 않습니다. 이러한 특성을 무기억성(memoryless)이라 합니다. 즉, λ > 0인 모수를 가진 지수분포를 따르는 연속변수 X는 무기억성 확률변수로서 식 3으로 증명됩니다.

$$\begin{array}{rl}P(X>x+a|X>a)& =P(X>x)\\ F(x)& ={\int }_0^x\lambda e^{-\lambda t}dt\\ & =1-e^{-\lambda x}\\ P(X>x+a|X>a)& =\frac{P(X>x+a)\cap P(X>a)}{P(X>a)}\\ & =\frac{P(X>x+a)}{P(X>a)}\\ & =\frac{1-F(x+a)}{1-F(a)}\\ & =\frac{e^{-\lambda (x+a)}}{e^{-\lambda a}}\\ & =e^{-\lambda x}\\ & =1-F(x)\\ & =P(X>x)\end{array}$$	(식 3)

예 2)

한 회사 직원의 분당 평균 통화 빈도는 0.6이라고 합니다. 통화 사이의 경과된 시간(분)을 확률변수 t로 할 경우 이 확률변수는 평균빈도 λ = 0.6인 지수분포에 적합시킬 수 있습니다. 이 분포의 확률밀도 함수는 다음과 같습니다.

f(t) = 0.6 exp(-0.6t),   t > 0

다음을 결정합니다.

1. P(t ≤ 5)
2. P(t ≥ 2)

a) 통화 후 다음 통화까지 5분 이내일 확률?

p=stats.expon.cdf(5, scale=1/0.6)
round(p, 3)

0.95

b) 통화 후 다음 통화가 2분 이상일 확률?

p2=stats.expon.sf(2, scale=1/0.6)
round(p2, 3)

0.301

예 3)

다음 일정기간의 kospi 주가의 고가와 저가의 일일변동율 기준으로 1% 이상이 되는 기간을 확률변수 x라고 할 때 그 변수를 지수분포에 적용하여 봅시다.

st=pd.Timestamp(2023,1, 30)
et=pd.Timestamp(2024, 5, 30)
da=fdr.DataReader("KS11", st, et)
change=(da['High']-da['Low'])/da['Low']*100
change.head(3)

Date
2023-01-30    1.820485
2023-01-31    1.407538
2023-02-01    0.818615
dtype: float64

위 자료로부터 x를 계산하기 위해 위 자료에서 1%이상인 경우를 1, 미만인 경우를 0으로 구분합니다. 변환을 위해 pd.cut() 함수를 적용합니다.

cat=pd.cut(change, bins=[-100, 1, 100], labels=[0,1])
cat.head(3)

Date
2023-01-30    1
2023-01-31    1
2023-02-01    0
dtype: category
Categories (2, int64): [0 < 1]

위 결과로부터 값이 1인 부분의 인덱스를 확인하고 인접한 인덱스 사이의 차이를 계산합니다. 이 차이는 1% 이상 상승이 일어나는 기간을 나타냅니다. 1%이상인 인덱스를 반환하기 위해 np.where() 함수, 그리고 인접한 두 값사이의 차이는 반복문을 적용합니다.

IncPos=np.where(cat==1)[0]
print(IncPos[:10])

[ 0  1  3  4  5 12 13 14 15 18]

x=np.array([1])
for i in range(1, len(IncPos)):
    indNum = IncPos[i]-IncPos[i-1]
    x=np.append(x,indNum)
print(x[:10])

[1 1 2 1 1 7 1 1 1 3]

위 결과의 고유값들과 각각에 대응하는 빈도수를 확인하기 위해 np.unique() 함수를 적용합니다. 각 고유값은 1% 상승을 나타내기 위한 기간이 되며 이 기간에 대한 확률은 빈도수를 사용하여 계산할 수 있습니다.

num, fre=np.unique(x, return_counts=True)
print(num)
print(fre)

[ 1  2  3  4  5  6  7 11 12 13 14]
[66 34  9  9  7  4  1  1  1  1  2]

prop=fre/np.sum(fre)
print(prop.round(3))

[0.489 0.252 0.067 0.067 0.052 0.03  0.007 0.007 0.007 0.007 0.015]

λ는 식 4와 같이 객체 x의 평균(기대값)으로부터 결정합니다.

$$E(X)\hspace{0.5em}=\hspace{0.5em}\frac{1}{\lambda }\to \lambda =\frac{1}{E(X)}$$

(식 4)

mu=np.mean(x)
print(mu.round(3))
lam=1/mu
print(lam.round(3))

2.422
0.413

그림 2는 위 결과들을 각 변수에 따른 확률을 나타내기 위해 막대그래프와 평균빈도(λ ≈ 0.413)를 모수로 하는 지수분포를 함께 나타냈습니다. 약간의 차이가 관찰되지만 전체적으로 지수분포에 적합됨을 알 수 있습니다.

그림 2. 예제 3의 자료에 대한 히스토그램과 지수분포(λ ≈ 0.413).

x=np.linspace(0.01, 14, 1000)
fig, ax=plt.subplots(figsize=(4,3))
y=stats.expon.pdf(x, lam, 1/lam)
idx=np.where(y>0)[0]
ax.bar(num, prop, color="g", alpha=0.3)
ax.plot(x[idx], y[idx],  color='g', label=r"$\lambda$="+str(lam.round(3)))
ax.set_xlabel("x", loc="right")
ax.set_ylabel("pdf", loc="top")
ax.legend(loc="best", frameon=False)
plt.show()

그림 3.2.7의 자료는 표준화하지 않은 것이므로 확률밀도함수는 stats.expon.pdf(x, loc, scale)를 적용하기 위해 loc과 scale을 명시하여야 합니다. 두 매개변수들은 각각 평균과 표준편차를 의미하는 것으로 λ와 $\frac{1}{\lambda }$를 사용합니다.

위 결과에 의하면 저가와 고가의 변화율이 1% 이상되는 사건간의 기간이 3일이하가 될 확률은 다음과 같습니다.

pMore3=stats.expon.cdf(3, lam, scale=1/lam)
pMore3.round(3)

0.656