감마 함수

1. 개요

감마 함수 Gamma Function
정의	$\Gamma(z) = \int_0^\infty t^{z-1}e^{-t}\,dt$ ($\mathrm{Re}\,z > 0$)
계승	$\Gamma(n+1) = n!$
점화식	$\Gamma(z+1) = z\,\Gamma(z)$
반사 공식	$\Gamma(z)\Gamma(1-z) = \pi/\sin(\pi z)$
극	$z = 0, -1, -2, \dots$ 에서 1차 극, 영점 없음
배정밀도 넘침	$x > 171.6244$ 에서 $\Gamma(x)$ 가 오버플로
계산법	스털링 급수 + 점화식, 란초스 근사, 스푸지 근사

계승을 정수 밖으로 늘린 함수. 그리고 확률분포·고차원 기하·분수 미분 어디에나 상수처럼 끼어드는 함수.

감마 함수(gamma function)는 계승 $n!$ 을 복소수로 확장한 특수함수로, 실수부가 양수인 $z$ 에서

$$\Gamma(z) = \int_0^\infty t^{z-1} e^{-t}\,dt$$

로 정의된다. 부분적분 한 번으로 $\Gamma(z+1) = z\Gamma(z)$ 가 나오고 $\Gamma(1) = 1$ 이므로 $\Gamma(n+1) = n!$ 이다. 오일러가 1729년 골트바흐에게 보낸 편지에서 계승의 보간 문제를 무한곱으로 풀었고, 위 적분형과 $\Gamma$ 라는 기호는 그 뒤 르장드르를 거쳐 정착했다. 「$n!$ 이 아니라 $\Gamma(n+1)$」 이라는 한 칸 밀린 정의 때문에 지금까지도 계산 실수가 끊이지 않는다¹.

시뮬레이션 실무에서 감마 함수를 직접 공부할 일은 드물지만, 계산해야 할 일은 매일 있다. 포아송·이항 우도의 $\ln k!$, 카이제곱·스튜던트 t 분포의 정규화 상수와 p-값, 감마 분포 난수, $n$ 차원 공의 부피, 분수계 미분의 계수가 모두 감마 함수다. 이 문서는 그래서 성질은 쓰임새에 필요한 만큼만 정리하고, 어떻게 정확하고 넘치지 않게 계산하는가에 무게를 둔다.

2. 꼭 알아야 할 성질

• 특수값. $\Gamma(1/2) = \sqrt\pi$ (가우스 적분 $\int e^{-x^2}dx = \sqrt\pi$ 와 같은 식), 따라서 $\Gamma(3/2) = \sqrt\pi/2$, $\Gamma(-1/2) = -2\sqrt\pi \approx -3.5449$. 반정수 감마는 모두 $\sqrt\pi$ 의 유리수배다.
• 해석적 연속. 점화식을 거꾸로 $\Gamma(z) = \Gamma(z+1)/z$ 로 쓰면 실수부가 음수인 곳까지 정의가 넓어진다. $z = 0, -1, -2, \dots$ 에서 1차 극을 갖고, 영점은 하나도 없다. 그래서 $1/\Gamma(z)$ 는 복소평면 전체에서 해석적인 정함수이고, 수치 라이브러리가 종종 $1/\Gamma$ 를 따로 제공하는 이유가 이것이다(극 근처에서 0 으로 매끄럽게 간다).
• 반사 공식 $\Gamma(z)\Gamma(1-z) = \pi/\sin(\pi z)$. 음수 인자의 값을 양수 인자로 바꿔 계산하는 표준 경로다.
• 배가 공식(르장드르) $\Gamma(z)\Gamma(z+\tfrac12) = 2^{1-2z}\sqrt\pi\,\Gamma(2z)$.
• 양의 실축 위 모양. $x > 0$ 에서 볼록하고, $x_{\min} = 1.461632\ldots$ 에서 최솟값 $0.885603\ldots$ 을 찍는다. 보어-몰레루프 정리에 따르면 $\Gamma$ 는 「$f(1) = 1$, $f(x+1) = xf(x)$, $\ln f$ 가 볼록」 세 조건을 만족하는 유일한 양의 함수다. 계승의 보간 방법이 무수히 많은데도 다들 이것을 쓰는 이유다.
• 베타 함수와의 관계 $B(a, b) = \int_0^1 t^{a-1}(1-t)^{b-1}dt = \Gamma(a)\Gamma(b)/\Gamma(a+b)$. 베타 분포와 이항계수의 연속판이다.

3. 계산 — 넘치지 않게, 정확하게

3.1. 로그 공간에서 놀아라

$\Gamma(171) = 170! \approx 7.257\times10^{306}$ 이고, 배정밀도 최댓값을 넘는 경계는 $x = 171.62437\ldots$ 다. $\Gamma(172)$ 부터는 inf 다. 그런데 포아송 우도에 $k = 500$ 번 사건이 들어가는 일은 흔하다. 그래서 실무 코드는 $\Gamma$ 대신 $\ln\Gamma$ (C 의 lgamma, SciPy gammaln) 를 쓰고 끝까지 로그 공간에서 계산한다. 예컨대 $\ln\Gamma(1000) = 5905.22$ 는 아무 문제 없이 나온다. 이항계수 $\binom{n}{k}$ 도

$$\ln\binom{n}{k} = \ln\Gamma(n+1) - \ln\Gamma(k+1) - \ln\Gamma(n-k+1)$$

로 계산하고 마지막에 필요할 때만 지수를 취한다. 다만 $n$ 이 크고 $k$ 가 작으면 큰 수끼리 빼는 소거 오차(부동소수점 연산)가 생기므로, 정밀도가 중요하면 $\ln\Gamma(x) -$ 스털링 주항 꼴의 보정 함수(「스털링 오차항」)를 쓰는 구현이 더 정확하다(로더의 이항 확률 알고리즘이 이 방식이다)².

3.2. 큰 인자 — 스털링 급수

오일러-매클로린 공식을 $\sum\ln k$ 에 쓰면

$$\ln\Gamma(x) \sim \Bigl(x - \tfrac12\Bigr)\ln x - x + \tfrac12\ln 2\pi + \frac{1}{12x} - \frac{1}{360x^3} + \frac{1}{1260x^5} - \cdots$$

가 나온다. 이 급수는 발산하는 점근 급수 라 항을 늘려서 정확도를 사지 못하고, 대신 $x$ 가 크면 몇 항만으로 배정밀도에 닿는다. $x = 10$ 에서 주항 $\sqrt{2\pi/x}(x/e)^x$ 만 쓰면 상대오차가 $-8.3\times10^{-3}$, $1 + 1/(12x)$ 를 곱하면 $-3.2\times10^{-5}$ 다. 작은 $x$ 는 점화식 $\Gamma(x) = \Gamma(x+m)/[x(x+1)\cdots(x+m-1)]$ 로 큰 쪽으로 밀어 놓고 스털링을 쓴다. 이것이 고전적인 「밀고 나서 스털링」 전략이다.

3.3. 란초스 근사와 스푸지 근사

란초스(1964)는 점근 급수가 아니라 수렴하는 근사를 냈다.

$$\Gamma(z+1) = \sqrt{2\pi}\,\Bigl(z + g + \tfrac12\Bigr)^{z+1/2} e^{-(z+g+1/2)} \Bigl[c_0 + \sum_{k=1}^{n}\frac{c_k}{z+k}\Bigr] + \varepsilon$$

자유 매개변수 $g$ 와 항 수 $n$ 을 고르면 계수 $c_k$ 가 정해지고, 오른쪽 반평면 전체에서 오차가 균일하게 작다. 많이 퍼진 $g = 7$, $n = 9$ 계수 세트로 $0.6 \le x \le 140$ 을 재 보면 최대 상대오차가 $1.5\times10^{-14}$ 다. 같은 범위에서 시스템 math.gamma 는 $6.7\times10^{-16}$ 이니, 란초스 9항은 「배정밀도에서 두 자리 정도 아쉬운」 실용 근사다. 스푸지(1994)는 비슷한 꼴에서 계수를 닫힌꼴로 주고 오차 한계를 명시적으로 보장하는 변형을 냈다. 원하는 정밀도에 맞춰 항 수를 정할 수 있어 임의 정밀도 구현에서 쓰인다. 음수 인자는 둘 다 반사 공식으로 넘긴다.

현대 수학 라이브러리(glibc, Boost.Math 등)는 구간별 유리함수 근사 + 점화식 + 스털링을 섞어 마지막 비트까지 맞추는 쪽으로 간다. 직접 구현할 이유는 거의 없고, 구현해야 한다면 GPU 커널 안처럼 라이브러리가 없는 곳일 것이다. 그때 란초스가 짧고 분기가 없어 국룰이다.

4. 시뮬레이션에서 만나는 자리

확률분포. 감마 분포 $f(x) = x^{k-1}e^{-x/\theta}/(\Gamma(k)\theta^k)$ 는 대기 시간·강수량·켤레사전분포(포아송 비율의 사전분포)에 쓰인다. 그 난수는 마르살리아·창(2000)의 기각 표본추출이 표준이다 — 시도마다 정규 난수 하나와 균등 난수 하나를 쓰고, 형상 모수 $k \ge 1$ 에서 수락률이 95% 를 넘는다($k < 1$ 은 $k+1$ 로 뽑고 $U^{1/k}$ 를 곱해 변환한다). 카이제곱 분포의 누적확률은 정칙화 불완전 감마 함수 $P(a, x) = \gamma(a, x)/\Gamma(a)$ 이고, 몬테카를로 방법의 적합도 검정이나 마르코프 연쇄 몬테카를로의 수렴 진단에서 p-값을 낼 때 이것을 부른다.

고차원 기하. 반지름 1 인 $n$ 차원 공의 부피는

$$V_n = \frac{\pi^{n/2}}{\Gamma(n/2 + 1)}$$

이다. $n = 1, 2, 3, \dots$ 에서 $2,\ 3.1416,\ 4.1888,\ 4.9348,\ \mathbf{5.2638},\ 5.1677,\ 4.7248, \dots$ 로 5차원에서 최대를 찍고 그 뒤 0 으로 간다. $n = 14$ 이면 $0.599$ 다. 같은 공을 감싸는 정육면체 $[-1,1]^n$ 의 부피는 $2^n$ 으로 커지니, 고차원에서 「공 안에 떨어지는 균등 난수」의 비율은 순식간에 0 이 된다. 기각 표본추출이 고차원에서 죽는 이유이자 차원의 저주의 가장 짧은 증명이다.

분수계 미분. 리만-리우빌·카푸토 분수 미분의 정의에 $1/\Gamma(1-\alpha)$ 가, 그륀발트-레트니코프 이산화의 가중치 $w_k = (-1)^k\binom{\alpha}{k}$ 에 일반화 이항계수 $\binom{\alpha}{k} = \Gamma(\alpha+1)/(\Gamma(k+1)\Gamma(\alpha-k+1))$ 이 들어간다. 이상 확산·점탄성 모형에서 이 가중치를 직접 감마로 계산하면 $k$ 가 커질 때 넘치거나 소거가 생기므로, 실무에서는 점화식 $w_0 = 1$, $w_k = w_{k-1}(1 - (\alpha+1)/k)$ 로 만든다. 감마 함수를 계산하지 않는 것이 감마 함수를 다루는 최선의 방법인 경우다³.

양자·통계역학. 가우스형 기저함수의 겹침 적분($\int r^{n}e^{-\alpha r^2}$ 꼴), $d$ 차원 이상기체의 상태 수, 보스-페르미 적분의 정규화가 모두 반정수·정수 감마로 떨어진다. 대부분 닫힌꼴이 있으니 표로 박아 두는 게 보통이다.

5. 관련 문서

• 오일러-매클로린 공식 · 점근 급수 · 복소해석 · 테일러 급수
• 부동소수점 연산 · 몬테카를로 방법 · 마르코프 연쇄 몬테카를로 · 켤레사전분포
• 차원의 저주 · 중심극한정리 · 라플라스 근사
• 베타 함수 · 디감마 함수 · 감마 분포 · 분수계 미적분

6. Footnotes

1. 르장드르가 왜 한 칸 밀어 정의했는지는 아직도 정설이 없다. 가우스는 $\Pi(n) = n!$ 이라는 밀리지 않은 기호를 썼는데 졌다. 표준은 품질이 아니라 점유율로 정해진다는 수학판 VHS 대 베타맥스. ↩

2. 우도를 확률 그대로 곱하다가 0.0 이 나와서 MCMC 가 첫 스텝부터 전부 기각되는 사고는 거의 통과의례다. 로그 공간은 선택이 아니라 생존이다. ↩

3. 그륀발트 가중치를 gamma(alpha+1)/(gamma(k+1)*gamma(alpha-k+1)) 로 그대로 짜면 $k = 171$ 근처에서 inf/inf = nan 이 반겨 준다. 점화식 한 줄이면 끝날 일을 킹받게 디버깅하게 된다. ↩