[R] Poisson Process

Definition of Poisson Process

임의의 시점에서 사건이 발생한다고 가정하고, $N(t)$를 구간 [0,t]에서 발생한 사건의 수라고 하자.
이 사건들이 비율(Rate)이 $\lambda(>0)$인 푸아송 과정을 형성한다고 할 때, 다음의 조건들을 만족한다.

• $N(0)=0.$
• 서로 다른 구간에서 발생한 사건의 수는 독립.
• 주어진 구간에서 발생한 사건 수의 분포는 그 구간의 위치와 상관없이 오직 구간의 길이에만 의존
• $\lim_{h \rightarrow0}\frac{P\{N(h)=1\}}{h}=\lambda$
• $\lim_{h \rightarrow0}\frac{P\{N(h)\geq 2\}}{h}=0$

첫 번째 조건은 이 과정이 0에서 시작함을 의미한다.
두 번째 조건은 $N(t)$와 $N(t+s)-N(t)$가 독립임을 의미한다.
세 번째 조건은 $N(t+s)-N(t)$가 $t$값에 관계 없이 오직 구간의 길이 $s$에만 의존함을 의미한다.
네 번째와 다섯 번째 조건은 매우 작은 길이의 $h$에서 사건이 1번 발생할 확률이 $\lambda h$에 수렴하는 반면,
2번 발생할 확률은 0에 수렴함을 의미한다.

이러한 가정을 만족하는 길이 $t$의 구간에서 발생한 사건 수 $N(t)$가 평균이 $\lambda t$인 포아송 분포를 따른다는 사실이 알려져 있으므로, 우리는 이 사실을 받아들이고 포아송 과정의 성질을 알아보자.

포아송 과정의 성질

• 모든 $t$에 대하여 $N(t)\sim Poission(\lambda t)$
• 도착간격시간(Inter-arrival times)의 분포는 $iid \ exp(\lambda)$
• $N(t) = n$이 주어졌을 때, 도착시각(arrival times) $S_1,...,S_n$의 분포는 구간 $(0,t)$에서 추출한 $n$개의 균일확률변수의 순서통계량의 분포와 같다.
• 마지막 도착시각 $S_n$의 값이 $S_n=t$로 주어졌을 때, 나머지 $(n-1)$개의 도착시각의 분포는 구간$(0,t)$에서 추출한 $(n-1)$개의 균일확률변수의 순서통계량의 분포와 같다.

Generating a Poisson Process

[방법 1] Poission Process의 두 번째 성질을 이용한 방법이다.
$\lambda=1$, 구간은 $(0,T=10)$이다.

> lambda = 1
> T = 10
> t = 0; I = 0
> S = numeric()
> repeat{
+   t = t + rexp(1, lambda)
+   if (t > T) break
+   I = I + 1; S[I] = t
+ }
> cat("Total number of events is", I, "\nArrival times are \n", signif(S,3), "\n")
Total number of events is 10 
Arrival times are 
 0.566 0.672 0.731 1.31 5.27 6.44 7.44 8.87 8.91 9.24 

S에는 도착시각을 저장한다.
도착간격시간(구간의 길이)가 지수분포를 따른다는 점을 고려하여 다음 $t$를 생성한다.

[방법 2] 다음은 Poission Process의 첫 번째 성질과 세 번째 성질을 이용한 방법이다.
구간 $(0,T)$에서의 사건의 수가 포아송 분포를 따르는 것을 이용해 $n$을 생성한 뒤, 그 개수만큼 균일확률난수를 만들어 정렬하면 된다.

> lambda = 1
> T = 10
> n = rpois(1, lambda*T)
> S = sort(T*runif(n))
> cat("Total number of events is", n, "\nArrival times are \n", signif(S,3), "\n")
Total number of events is 9 
Arrival times are 
 1.29 1.61 4.32 4.4 5.53 6.09 6.56 8.14 9.92 

어떤 방법이 더 효율적일까? [방법 2]는 repeat를 사용하지도 않고 로그 연산을 할 필요가 없어 더 효율적이다. 하지만 도착시간을 순서대로 알고 싶다면 sorting을 해야하는 점을 기억하자.

이제 아래 예제를 풀어보자.

버스가 시간당 5대의 비율로 포아송 과정에 따라 도착하는 사건이 있다고 하자. 각 버스는 20명, 21명, …, 40명의 사람을 태우고 올 확률이 동일하며, 서로 다른 버스에 있는 사람의 수는 독립적이다. 시간 $t=1$까지 도착하는 사람들의 수를 시뮬레이션하라.

1. $\lambda=5$, 전체 구간은 $(0,1)$이다.
2. 포아송 과정의 첫 번째 성질에 따라 $n$은 평균이 5인 포아송 분포를 따를 것이고,
   $S$는 $n$개의 Uniform random variable가 정렬된 벡터이다.
3. 버스에 탑승할 수 있는 사람의 수는 20~40명이 가능한데, 각 경우의 발생 확률은 동일하다.
4. 따라서 20에서 40 사이의 수에서 $n$개 만큼 복원추출한다.
5. 각 시점(도착시간)마다의 누적합을 구한다.

위 절차를 R 코드로 옮기면 다음과 같다.

> T = 1; lamb = 5
> n.bus = rpois(1, T*lamb)
> S = sort(T*runif(n.bus))
> X = sample(20:40, n.bus, replace = TRUE)
> Y = cumsum(X)
> cat("도착시간은", signif(S,3),"이고, \n각 시각까지 도착한 사람들의 누적합은",Y,"이다.")
도착시간은 0.865 이고, 
각 시각까지 도착한 사람들의 누적합은 23 이다.

Definition of Nonhomogeneous Poisson Process

포아송 과정의 단점은 모든 동일한 길이의 구간에서 사건 발생 확률이 동일하다는 가정이 필요하다는 것인데, 이 가정을 일반화한 형태를 비동질적 포아송 과정(Nonhomogeneous Poisson Process) 라고 한다.

만약 사건들이 구간에서 무작위로 발생하고 $N(t)$가 시간 $t$까지 발생한 사건의 수라고 한다면, $\{N(t),\ t\geq 0\}$가 강도 함수(intensity function) $\lambda(t)$를 갖는 비동질적 포아송 과정을 형성한다고 한다.

고정된 비율($\lambda$)를 가정하지 않아 덜 제한적인 탓에 응용범위가 넓으나, 해석적(Analytical) 결과를 얻기가 어렵다. 따라서 시뮬레이션을 이용한 연구가 필요하다.

함수 $m(t) = \int_0^t\lambda(s)ds,\ t\geq0$을 평균값 함수(mean value function)이라 부르는데,
다음과 같은 결과가 성립하는 것이 알려져 있다.

Generating a Nonhomogeneous Poisson Process

[방법 1] Thinning Algorithm
$\lambda(t)\leq\lambda$ 일 때 비율이 $\lambda$인 동질적 포아송 과정에서 생성한 사건을 확률 $\lambda(t)/\lambda$로 받아들이는 방법.

다음 예제를 R 코드로 구현해보자.

Intensity function이 주어진 Nonhomogeneous Poisson Process의 첫 10개의 Arrival times를 생성하기 위해 thinning 알고리즘을 사용하는 프로그램을 작성하시오.

> lambda = 7; T = 10
> lt = function(t) 3 + 4/(t+1)
> t = 0; I = 0
> S = numeric()
> repeat{
+   t = t + (-log(runif(1))/lambda)
+   if (t > T) break
+   if (runif(1) < lt(t)/lambda) {I = I + 1; S[I] = t}
+ }
> cat("Total number of events is", I, "\nFirst 10 arrival times are \n", signif(S[1:10],3), "\n")
Total number of events is 51 
First 10 arrival times are 
 0.0267 0.176 0.226 0.344 0.507 0.661 0.674 0.899 1 1.67 

$\lambda(t)\leq \lambda$를 만족해야 하므로 구간 $(0,t)$에서 $\lambda(t)$의 최댓값을 $\lambda$로 선택한다.

[방법 2] Superposition of two Poisson Processes
$\{N(t),\ t\geq 0\}$와 $\{M(t),\ t\geq 0\}$가 각각 강도함수 $\lambda(t),\mu(t)$인 비동질적 포아송 과정이고 독립이면, $\{N(t)+M(t),\ t\geq 0\}$는 강도함수가 $\lambda(t)+\mu(t)$인 비동질적 포아송 과정이 된다.

이 방법으로 [방법 1]에서 풀었던 예제를 다시 구현해보자.

hpp = function(lambda, t0 = 0, T) {
  t = t0; I = 0L # 64비트 정수형
  S = numeric()
  repeat {
    t = t + (-log(runif(1)))/lambda # Inverse Transform Algorithm of exp
    if (t > T) break
    I = I + 1; S[I] = t
  }
  list(n_events = I, arrival_times = S)
}

thin = function(max_lambda, t0 = 0, T,
                 lt = function(t) {}) {
  t = t0; I = 0L
  S = numeric()
  repeat {
    t = t + (-log(runif(1))/max_lambda)
    if (t > T) break
    if (runif(1) < lt(t)/max_lambda ) {I = I + 1; S[I] = t}
  }
  list(n_events = I, arrival_times = S)
}

> out1 = hpp(lambda = 3, T = 10)
> out2 = thin(max_lambda = 4, T = 10, lt = function(t) 4/(t+1))
> (n_events = out1$n_events + out2$n_events)
[1] 47
> (S = sort(c(out1$arrival_times, out2$arrival_times)))
 [1] 0.06743639 0.13243070 0.16844648 0.17034790 0.21404090 0.28172300 0.47478777 0.80299329 0.94983608 0.99192496
[11] 1.11286585 1.29620768 1.32263685 1.67114594 1.79144230 1.90126647 2.04186626 2.11612510 2.47939565 2.85760740
[21] 3.13653807 3.19606082 4.15696837 4.17959084 5.04824600 5.11351029 5.23141325 5.23296518 6.00666030 6.05085542
[31] 6.41131756 6.45372540 6.47980505 6.69270699 6.70126079 6.96807585 7.32528133 7.60655091 7.87403296 8.23904285
[41] 8.27249489 8.66432102 8.97968172 9.37242501 9.41551035 9.71788683 9.96074969

$\lambda(t)=3+\frac{4}{t+1}, \ 0\leq t \leq 10$을 $\lambda(t)=3$, $\mu(t)=\frac{4}{t+1}$인 두 비동질적 포아송 과정으로 각각 생성한 후 병합하면 된다.

Arrival times를 병합만 하면 사건 발생 순서대로 정렬되지 않으니 주의하자.

[방법 3] 구간을 쪼개서 각 구간에 Thinning Algorithm을 적용하는 방법
$i$번째 구간에서 $\lambda(t)\leq \lambda_i$를 만족할 때 적용할 수 있으며, [방법 1]보다 빠르다. 그 이유는 지수분포의 무기억성을 이용해 생성하는 지수난수의 개수를 줄일 수 있기 때문이다.

위 방법을 이용하여 다음 예제를 구현해보자.

다음과 같은 Intensity function을 갖는 Nonhomogeneous Poisson Process의 첫 10개의 Arrival times를 효율적으로 생성하는 알고리즘을 제시하라.

각 구간$(t_0,T)$에서 강도함수의 최댓값을 $\lambda_i$라고 한 뒤, [방법 2]처럼 각각의 비동질적 포아송 과정을 생성한 뒤 병합한다.

> out1 = thin(max_lambda = 1, t0 = 0, T = 5,
+             lt = function(t) t/5)
> out2 = thin(max_lambda = 26, t0 = 5, T = 10,
+             lt = function(t) 1 + 5*(t - 5))
> (n_events = out1$n_events + out2$n_events)
[1] 74
> (S = c(out1$arrival_times, out2$arrival_times))
 [1] 2.274273 4.369036 5.145196 5.904052 6.019631 6.027422 6.040848 6.170788 6.441836 6.455370 6.510889 6.568357
[13] 6.589603 6.637936 6.711909 6.730458 6.840531 6.900103 6.981688 7.410973 7.421870 7.450231 7.507743 7.755879
[25] 7.776126 7.787563 7.799275 7.885053 7.960399 7.966983 8.027413 8.076848 8.147335 8.151058 8.201475 8.239681
[37] 8.608568 8.614472 8.621317 8.655787 8.656616 8.685017 8.757989 8.854468 8.947832 8.973323 9.052801 9.164173
[49] 9.206218 9.255868 9.271815 9.341236 9.345004 9.351239 9.402605 9.405267 9.406409 9.455420 9.497963 9.521220
[61] 9.522033 9.580863 9.633107 9.647237 9.688557 9.707904 9.721669 9.746601 9.805789 9.847475 9.891562 9.901327
[73] 9.953214 9.970617