[확률과정] 푸아송 과정 실전편 - 발생 시간 간격의 관점에서

 이전 푸아송 과정을 설명하는 글에 이은 실전편. 푸아송 과정을 사건 발생 시간 간격의 관점에서 분석해보자.

시간 간격의 관점에서 분석하는 이유

---

 이전 푸아송 과정 글에서 최종적으로 유도한 식은 다음과 같다

 한 step(t=1)에서 λ의 기대 값을 가지는 사건이 발생할 경우, 시간 간격 t 내에 k번의 사건이 발생할 확률을 나타내는 식이다. 이처럼 사건의 발생 원인을 확실히 알고 있는 경우, 즉, 사건이 step별로 독립적이며 λ라는 기대값을 가진다는 것을 알고 있는 경우라면, 위의 식이 기대값을 도출하는데 도움이 될 수 있다.

 허나, 아무것도 모르는 상태에서 이 사건이 푸아송 과정인지 살펴보기 위해서는 위의 식으로는 불편한 점이 있다. 그래서 사건이 발생한 시간 간격의 관점에서 분석하는 것이 편리하다.

 이 관점이 편리한 이유는 다음과 같다.

 예를 들면, 사건이 다음과 같이 발생했다는 Data가 있다.

 t = 1, 3, 7, 13, 16, 18, ....

 그럼 시간 간격은 다음과 같다.

 T = 1, 2, 4, 6, 3, 2, ... 

 이제, T를 분석하면 푸아송 과정 여부를 유추할 수 있게 되는 것이다.

---

수식

---

 

 먼저 첫번째 식은 발생 시간 간격 T_i 가 특정 시간간격 t 보다 작을 확률을 나타낸다. 이 때, T_i 는 모든 i에 대해서 독립적이고 같은 분포를 따른다. 사건이 매 step 마다 같은 발생 확률을 지니므로 서로 다른 T_i가 얽힐 이유가 없고, 또한 T_i들의 발생 원인이 같기 때문이다. 

 첫째 식이 둘째 식과 같은 이유는 시간간격 t 내에 T_i가 존재할 확률은 시간간격 t에서 사건이 1번 이상 발생할 확률과 같은 의미 이기 때문이다.

 둘째 식은 시간간격 t 내에 사건이 발생하지 않을 확률을 1에서 뺀 값과 같고, 이는 네번째, 다섯번째 식으로 전개할 수 있다. (이전글 참조)

---

예제

---

 이제 우리가 해야할 일은 시간간격 별 빈도를 정규화된 histogram으로 나타내고(0~1의 확률), 이를 축적해 나가는 Cumulative distribution 으로 나타낸다. Cumulative distribution를 간단히 설명하자면,

  P(t=1)=0.1, P(t=2)=0.2, P(t=3)=0.3 ....

와 같은 이산 확률 분포는 Cumulative distribution에서

  P(t<=1)=0.1, P(t<=2)=0.1+0.2, P(t<=3)=0.1+0.2+0.3 ...

와 같은 값을 지니게 된다. 고로, P(T_i = t)를 P(T_i <= t)로 바꾸는 작업이라고 할 수 있다.

 이제 이 분포를 1 - e^(-λt) 식에 핏팅시켜본다. 핏팅이 잘 이루어진다면 우리는 이 분포가 푸아송 분푸라고 추측할 수 있고, 사건 발생 기대값에 대해 논해볼 수 있다. 이것이 사건 발생 Data로 부터 발생 원인을 추정해나가는 inverse statistics methodology이다. 

 시뮬레이션 결과와 수식으로 계산한 확률분포의 비교 결과는 다음과 같다.

λ = 0.1 

step = 100,000

 확률분포의 가시화와 fitting은 관습적으로 Cumulative를 반대방향, 즉, P(T<=t) 를 P(T>=t)로 표현하고, 1 - e^(-λt) 를 e^(-λt)로 표현한다. 또 분포가 너무 넓게 퍼져있어 형태가 불분명하다면 log스케일을 취해 형태를 알기 쉽게 표현하기도 한다. 여기선 exponential 함수이므로 y축을 log스케일로 표시해 보았다. Power law distribution과 같은 두터운 꼬리를 지닌 분포는 x축도 log스케일을 취해 log-log 스케일로 표현하기도 한다.

 시뮬레이션 값과 수식으로 계산한 값이 잘 일치하고 있는 것이 보인다. 

 이번엔 e^(-λt) 함수에 fitting 시켜 λ를 구해보자.

Simulation에서 설정한 0.1값이 비슷한 값을 구할 수 있다. 이렇게 사건이 내재한 원인을 모름에도 불구하고, 푸아송 과정이라고 추측함으로써, 이 사건의 발생 기대값을 역으로 유추할 수 있게 된다. 사실, 이 분포에는 많은 정보가 있지만 또한 많은 정보가 누락되기도 한다. 그래서 무조건 이 분포를 신뢰해서 이 사건은 독립적인 발생확률을 지닌 푸아송 과정이라고 추측하는 것은 옳지 않다.

 하지만, 사건을 들여다 볼 수 있는 그 첫번째 단서로써, 또한 다른 사건들과 비교할 수 있는 도구로써 중요한 의미를 지니고 있으며, 푸아송 분포 이외에도 다양한 분포를 통한 함수 fitting은 우리가 추측할 수 없던 복잡한 사건들에 대해 단서를 얻을 수 있는 주요한 방법이다.

Code

import matplotlib.pyplot as plt
import numpy as np
import scipy.optimize as op


ti = 0
c = 0.1   # 사건 발생 확률 λ
ti_log = []     # T_i 로그
for step in range(100000):
    p = c
    status = np.random.choice([1, 0], p=[p, 1-p])
    if status == 1:     # 사건 발생
        ti_log += [ti]
        ti = 0
    else:   # 사건 발생안함
        ti += 1


""" Simulation vs Expected """
fig = plt.figure()
plt.hist(ti_log, bins=max(ti_log)-min(ti_log), normed=True, cumulative=-1, label="simulated")
x = np.arange(0, 100, 1)
y = np.exp(-c*x)
plt.plot(x, y, label='expected')
plt.legend()
plt.yscale('log')
plt.xlabel('t')
plt.ylabel('P(T_i>=t)')
plt.show(block=False)


""" Simulation vs Fitting """
def func_exp(T, alpha):
    y = np.exp(-alpha * T)
    return y


fig = plt.figure()
q = plt.hist(ti_log, bins=max(ti_log)-min(ti_log), normed=True, cumulative=-1, label="simulated")
x = q[1][:-1]
y = q[0]
popt, pcov = op.curve_fit(func_exp, x, y)
plt.plot(x, func_exp(x, popt[0]), label='fitted lamda: '+str(popt[0]))
plt.legend()
plt.yscale('log')
plt.xlabel('t')
plt.ylabel('P(T_i>=t)')
plt.show(block=False)

---