누적 이변량 정규분포의 재귀적 수치 평가

초록

본 논문은 Marsaglia가 제안한 단변량 정규분포 누적함수(Φ) 계산 기법을 확장하여, 이변량 정규분포 누적함수(Φ₂)의 고정밀·고속 수치 평가 알고리즘을 제시한다. 재귀적 전개와 유리함수 근사, 그리고 조건부 정규분포의 성질을 이용해 수식적으로 투명한 구조를 유지하면서도 기존 대표 알고리즘(Genz, Drezner 등)과 견줄 만한 실행 속도와 오차 한계를 제공한다. 또한 MPFR 기반 구현을 통해 임의 정밀도 연산이 가능하도록 설계되었다.

상세 분석

이 논문은 누적 이변량 정규분포 Φ₂(x, y; ρ)의 계산을 “조건부 누적 정규분포”라는 관점에서 재구성한다. ρ가 상관계수일 때, Φ₂는 한 변수에 대해 적분한 형태, 즉 Φ₂(x, y; ρ)=∫_{-∞}^{x} φ(t) Φ((y-ρt)/√(1-ρ²)) dt 로 표현된다. 여기서 φ는 표준 정규밀도, Φ는 단변량 누적 정규분포이다. Marsaglia가 제시한 Φ의 고속 근사식은 “rational approximation + continued fraction” 구조를 갖는데, 이는 입력값의 절대값에 따라 서로 다른 다항식 계수를 사용해 15~16자리 정확도를 보장한다. 저자는 이 구조를 그대로 차용하면서, 적분 구간을 재귀적으로 분할하고 각 구간에서 Φ를 호출하는 방식으로 Φ₂를 계산한다.

핵심 아이디어는 다음과 같다.

1. 대칭성 및 변환: (x, y, ρ) → (|x|, |y|, |ρ|) 로 변환해 양의 영역만 다루고, 필요 시 회전 변환을 통해 ρ가 0에 가까운 경우와 1에 가까운 경우를 별도로 처리한다.
2. 재귀적 구간 분할: 적분 구간을 가우스-레전드르(또는 등간) 노드로 샘플링하고, 각 노드에서 조건부 Φ 값을 계산한다. 노드 수는 ρ의 크기에 따라 동적으로 조정되며, ρ가 작을수록 적은 노드로 충분히 수렴한다.
3. Marsaglia식 재활용: Φ 호출은 기존 Marsaglia 구현을 그대로 사용한다. 이때, 고정밀 연산이 필요하면 MPFR 라이브러리로 래핑해 임의 정밀도를 제공한다.
4. 오차 제어: 재귀 깊이와 노드 수를 오차 추정식 ε≈C·e^{-k·n} (C, k는 상수, n은 노드 수) 로 관리한다. 실험에서는 n=8~12이면 10^{-12} 수준의 절대 오차를 얻는다.

알고리즘의 복잡도는 O(N·M) 형태인데, N은 재귀 단계(보통 2~3), M은 적분 노드 수다. 이는 기존 Genz 알고리즘의 O(N³) 복잡도에 비해 현저히 낮다. 또한, 구현이 단순해 메모리 사용량이 적고, 벡터화가 용이해 현대 CPU/GPU 환경에서 추가적인 가속이 가능하다.

정밀도 확장 측면에서 저자는 MPFR 기반 구현을 제시한다. Marsaglia식의 유리함수 계수를 고정소수점이 아닌 다중 정밀 부동소수점으로 재계산하면, 100자리 이상의 정확도도 손쉽게 달성한다. 이는 금융·보험·통계 분야에서 극단적인 신뢰구간을 요구하는 상황에 직접 적용 가능하다.

실험 결과는 다음과 같다.

• 속도: C++ 구현 기준, ρ=0.5, (x,y)=(1.0,1.0)에서 평균 0.45 µs, Genz 0.68 µs, Drezner 0.52 µs.
• 정밀도: 절대 오차 10^{-14} 이하를 일관적으로 달성, 특히 ρ≈±1 근처에서도 안정적인 수렴을 보임.
• 확장성: MPFR 128비트(≈34자리)에서 0.78 µs, 256비트에서 1.12 µs로, 정밀도 증가에 비례적인 시간 상승을 보인다.

결론적으로, 이 논문은 “조건부 정규분포 + Marsaglia식”이라는 두 가지 강력한 수학적 도구를 결합해, 구현 난이도는 낮추면서도 성능·정밀도 면에서 기존 최고 수준을 능가하는 알고리즘을 제시한다. 향후 다변량 정규분포(3차원 이상) 확장이나, 베이지안 추정에서의 사전·사후 분포 계산에도 자연스럽게 적용될 여지가 있다.