포아송 분포

1. 개요

포아송 분포 Poisson Distribution
확률질량	$P(N=k) = e^{-\lambda}\lambda^k/k!$, $k = 0, 1, 2, \dots$
모수	비율 $\lambda > 0$ (구간당 평균 사건 수)
평균 · 분산	둘 다 $\lambda$
유래	푸아송 (1837), 이항 분포의 희귀 사건 극한
대기 시간	사건 간격은 지수분포, $k$ 번째까지는 감마(얼랑)
켤레 사전	감마 분포
시뮬레이션	타우 도약, 광자·입자 계수 잡음, 결함·고장 수

「드물게, 서로 상관없이, 평균 속도는 일정하게」 일어나는 사건을 세면 거의 언제나 이 분포가 나온다.

포아송 분포(Poisson distribution)는 정해진 시간·공간 구간에서 사건이 일어나는 횟수 $N$ 의 분포로,

$$P(N = k) = \frac{e^{-\lambda}\lambda^k}{k!},\qquad k = 0, 1, 2, \dots$$

이다. 시메옹 드니 푸아송이 1837년 판결 확률에 관한 책에서 이항 분포의 극한으로 유도했다. 모수는 평균 사건 수 $\lambda$ 하나뿐이고, 평균과 분산이 둘 다 $\lambda$ 라는 강한 제약이 이 분포의 정체성이자 약점이다.

시뮬레이션에서 포아송은 세 군데서 튀어나온다. 계수 잡음(검출기에 들어오는 광자·중성자·전자 수, 즉 샷 노이즈), 확률적 사건 생성(대기행렬 이론의 도착, 몬테카를로 수송의 충돌, 화학 반응 시뮬레이션의 타우 도약), 그리고 결함·고장 수의 통계 모형이다. 이 문서는 극한의 정확한 오차, 평균=분산 제약의 결과, 그리고 포아송 난수를 만드는 법과 거기서 실제로 나는 사고를 다룬다. 비율에 감마 사전을 두는 베이즈 쪽 이야기와 감마-포아송 혼합(음이항)은 감마 분포 문서가 다뤘다.

2. 포아송 과정과 기본 성질

비율 $\lambda$ 인 포아송 과정은 (1) 겹치지 않는 구간의 사건 수가 독립이고, (2) 길이 $t$ 구간의 사건 수가 $\mathrm{Poisson}(\lambda t)$ 인 계수 과정이다. 동치로, 사건 간격이 서로 독립인 평균 $1/\lambda$ 의 지수분포다. 지수분포의 무기억성 때문에 「지금까지 얼마나 기다렸든 다음 사건까지의 대기는 새로 시작」이고, 이게 포아송이 「서로 상관없는 사건」의 모형이 되는 이유다.

실무에서 매일 쓰는 성질 네 가지:

• 합. 독립인 $\mathrm{Poisson}(\lambda_1)$ 과 $\mathrm{Poisson}(\lambda_2)$ 의 합은 $\mathrm{Poisson}(\lambda_1 + \lambda_2)$ 다. 검출기 여러 개의 계수를 합쳐도 포아송이다.
• 솎아내기(thinning). 각 사건을 확률 $p$ 로 독립적으로 남기면 남은 것은 $\mathrm{Poisson}(p\lambda)$ 다. 검출 효율 $p$ 인 검출기를 지난 광자 수가 그대로 포아송인 이유다. 비균질 포아송 과정을 만드는 루이스-셰들러 알고리즘도 이것으로 작동한다 — 최대 비율로 사건을 뽑고 $\lambda(t)/\lambda_{\max}$ 확률로 남긴다.
• 총합 조건부는 다항. $N_1, \dots, N_m$ 이 독립 포아송일 때 $\sum N_i = n$ 을 조건으로 주면 $(N_1, \dots, N_m)$ 은 확률 $\lambda_i/\sum\lambda_j$ 의 다항 분포다. 히스토그램 빈의 계수를 다룰 때 「총 개수 고정」 여부에 따라 두 모형을 오가는 근거다.
• 얼랑 쌍대성. 「$t$ 까지 사건이 $k$ 번 이상」과 「$k$ 번째 사건이 $t$ 이전」은 같은 사건이다. 그래서 포아송 누적분포를 감마 분포의 불완전 감마 함수로 계산한다(감마 분포 참고).

3. 이항 분포의 극한 — 얼마나 정확한가

$n$ 번 독립 시행에서 각 성공 확률이 $p = \lambda/n$ 이면, 성공 횟수 $\mathrm{Bin}(n, \lambda/n)$ 은 $n \to \infty$ 에서 $\mathrm{Poisson}(\lambda)$ 로 간다. 「희귀 사건의 법칙」이다. 극한이 얼마나 빨리 오는지는 르캉(1960)의 부등식이 정량한다. 성공 확률이 $p_i$ 로 서로 달라도 되며, 두 분포의 전변동 거리가

$$d_{\mathrm{TV}}\Bigl(\textstyle\sum_i \mathrm{Bernoulli}(p_i),\ \mathrm{Poisson}\bigl(\sum_i p_i\bigr)\Bigr) \le \sum_i p_i^2$$

를 넘지 않는다. $\lambda = 2$ 로 실측하면

$n$	10	100	1000
실제 $d_{\mathrm{TV}}$	0.0522	0.00457	0.000452
르캉 상한 $\lambda^2/n$	0.4	0.04	0.004

로, 거리가 정확히 $1/n$ 속도로 줄고 상한은 대략 9 배 보수적이다. 요점은 상한이 $n$ 이 아니라 $\sum p_i^2$ 에만 매달린다는 것이다. 시행이 아무리 많아도 성공 확률이 큰 시행이 몇 개 섞여 있으면 포아송 근사가 깨진다. 「사건이 드물다」는 조건의 정확한 의미가 이것이다.

4. 평균 = 분산이라는 제약

포아송의 분산은 평균과 같으므로, 계수 데이터에서 표본분산/표본평균 비(산포 지수)가 1 에서 크게 벗어나면 포아송 모형이 틀린 것이다. 1 보다 크면 과대산포(사건이 뭉쳐 온다, 비율 자체가 흔들린다 — 음이항 분포로), 1 보다 작으면 과소산포(사건이 서로 밀어낸다 — 검출기 불감시간, 규칙적 도착)다. 실제 계수 데이터는 거의 언제나 과대산포다.

분산 안정화. 분산이 평균에 매달려 있으니 계수가 큰 빈과 작은 빈의 오차 막대가 다르다. 앤스컴(1948) 변환 $2\sqrt{N + 3/8}$ 은 분산을 대략 1 로 고정한다. 40만 표본으로 재면 $\lambda = 1, 3, 10, 50$ 에서 변환 후 분산이 0.717, 0.983, 1.000, 1.000 이다. $\lambda \ge 3$ 쯤부터 잘 듣고, $\lambda = 1$ 같은 저계수에서는 아직 부족하다. 저계수 영상 잡음 제거에서 「앤스컴 변환 → 가우스 잡음용 필터 → 역변환」 파이프라인이 쓰이는 이유이자, 저계수에서 그 파이프라인이 한계를 보이는 이유다.

0 개 관측의 신뢰구간. 계수가 $x$ 일 때 정확한(가우드 1936) 95% 양측 구간은 카이제곱 분위수로

$$\Bigl[\tfrac12\chi^2_{0.025}(2x),\ \tfrac12\chi^2_{0.975}(2x+2)\Bigr]$$

이다. $x = 0$ 이면 $[0, 3.69]$, $x = 1$ 이면 $[0.025, 5.57]$, $x = 5$ 면 $[1.62, 11.67]$ 이다. 단측 95% 상한은 $x = 0$ 에서 $-\ln 0.05 = 2.996$ — 「사건이 0 번이었으면 비율 상한은 3」이라는 3의 규칙이 여기서 나온다. 「$\sqrt{N}$ 오차 막대」는 $N = 0$ 에서 폭 0 을 주니 저계수에서 쓰면 안 된다.¹

5. 포아송 난수 — 무엇을 써야 하나

곱셈법(크누스). 균등 난수를 계속 곱해 $e^{-\lambda}$ 아래로 떨어질 때까지의 횟수에서 1 을 뺀다. 지수 간격을 누적해 시간 1 을 넘는 순간을 찾는 것과 같다. 평균 $\lambda + 1$ 개의 균등 난수가 들고(2만 회 실측: $\lambda = 3$ 에서 3.997 개, $\lambda = 30$ 에서 31.03 개) 코드가 세 줄이라 작은 $\lambda$ 에서는 국룰이다.

그런데 큰 $\lambda$ 에서 조용히 틀린다. $e^{-\lambda}$ 는 $\lambda > 745$ 에서 배정밀도 0 으로 언더플로한다. 그러면 종료 조건 「곱 ≤ 0」은 곱 자체가 언더플로할 때에야 참이 되고, 그 시점은 $\lambda$ 와 무관하게 약 745 번째 곱셈이다. 실제로 $\lambda = 1000$ 으로 2000 번 뽑으면 평균이 745.3 이 나온다. 에러도 경고도 없다.² 로그 공간에서 $-\ln U$ 를 더하는 판으로 바꾸면 언더플로는 사라지지만 비용이 여전히 $O(\lambda)$ 다.

역변환. 역변환 표본추출로 누적확률을 $k = 0$ 부터 더해 가며 균등 난수 $U$ 를 넘는 첫 $k$ 를 찾는다. 균등 난수가 하나뿐이라 준몬테카를로·공통 난수 기법과 잘 맞지만, 순차 탐색 비용이 역시 $O(\lambda)$ 다. 최빈값 $\lfloor\lambda\rfloor$ 에서 시작해 양쪽으로 탐색하면 $O(\sqrt\lambda)$ 로 줄어든다.

변환 기각(PTRS). 큰 $\lambda$ 에서는 회르만(1993)의 변환 기각법이 표준이다. 포아송 질량을 감싸는 모자 함수를 역변환으로 뽑고 기각표본추출로 보정하며, 시도 한 번에 균등 난수 2 개를 쓰고 수락률이 높아 기대 비용이 $\lambda$ 와 무관한 상수다. NumPy 도 $\lambda \ge 10$ 에서 이 방법을 쓴다.

확률질량 자체의 계산. $\lambda^k/k!$ 을 그대로 계산하면 $\lambda = k = 1000$ 에서 OverflowError 가 난다. 로그로 $k\ln\lambda - \lambda - \ln\Gamma(k+1)$ 을 계산해 지수를 취하면 0.0126146113487 로, 라이브러리 값과 유효숫자 11 자리가 맞는다. $\ln\Gamma$ 는 감마 함수 문서의 스털링 급수로 계산된다. 꼬리 확률이 필요하면 질량을 더하지 말고 정칙화된 불완전 감마 함수를 쓴다.

정규 근사의 함정. $\lambda$ 가 크면 $\mathrm{Poisson}(\lambda) \approx \mathcal N(\lambda, \lambda)$ 이지만 오른쪽으로 치우친 꼬리가 오래 남는다. 연속성 보정을 해도 $P(N \le \lambda)$ 가 $\lambda = 5$ 에서 정확값 0.616 대 근사 0.588, $\lambda = 50$ 에서 0.5375 대 0.5282 로 차이가 $1/\sqrt\lambda$ 속도로만 줄어든다. 꼬리 확률(희귀 사건 경보 문턱)을 정규 근사로 정하면 체계적으로 틀린다.

6. 시뮬레이션에서 — 타우 도약

화학 반응계를 분자 수준에서 정확히 시뮬레이션하는 길레스피 알고리즘은 반응 하나하나를 지수 대기 시간으로 뽑는다. 분자 수가 많으면 반응이 너무 잦아 느리다. 길레스피(2001)의 타우 도약(τ-leaping)은 시간 간격 $\tau$ 동안 반응 성향 $a_j$ 가 거의 일정하다고 보고, 각 반응 채널의 발생 횟수를 $\mathrm{Poisson}(a_j\tau)$ 로 한 번에 뽑는다. 바로 위의 「포아송 과정 = 일정 비율 사건 계수」를 그대로 쓴 것이다.

함정은 포아송이 상한이 없다는 점이다. 분자가 3 개 남은 종을 소비하는 반응에 $\mathrm{Poisson}(2)$ 를 뽑으면 확률 14% 로 4 회 이상이 나와 분자 수가 음수가 된다. 그래서 실용 구현은 $\tau$ 를 분자 수 변화가 작도록 고르고, 분자가 적은 반응은 정확한 길레스피 단계로 처리하거나 이항 도약으로 바꾼다. $a_j\tau$ 가 크면 포아송을 정규로 바꿀 수 있는데, 그 극한이 화학 랑주뱅 방정식, 곧 확률미분방정식이다.³

7. 관련 문서

• 감마 분포 · 음이항 분포 · 이항 분포 · 지수족 · 켤레사전분포
• 대기행렬 이론 · 브라운 운동 · 확률미분방정식 · 몬테카를로 방법
• 난수 생성기 · 역변환 표본추출 · 기각표본추출 · 감마 함수
• 중심극한정리 · 최대우도추정 · 길레스피 알고리즘

8. Footnotes

1. 「0 건이니 위험 0」 보고서가 반려되는 이유. 30 번 시험해 고장이 0 번이면 95% 신뢰로 말할 수 있는 건 「고장률이 시험당 0.1 이하」까지다. 0 은 증거가 아니라 표본이 작다는 고백이다. ↩

2. 「평균이 λ 인지」 한 줄짜리 단위 테스트만 있었어도 잡히는 버그인데, 다들 λ=3 으로만 테스트한다. 테스트 입력에 극단값 하나씩은 꼭 넣어라. ↩

3. 분자 3 개짜리 반응에 τ 를 크게 잡았다가 분자 수 −2 를 찍고 NaN 으로 터지는 건 타우 도약 입문자의 통과의례다. 음수 분자는 화학자에게 보여 주지 마라. ↩