디감마 함수

1. 개요

디감마 함수 Digamma Function
정의	$\psi(x) = \Gamma'(x)/\Gamma(x) = \frac{d}{dx}\ln\Gamma(x)$
점화식	$\psi(x+1) = \psi(x) + 1/x$
특수값	$\psi(1) = -\gamma = -0.577216$, $\psi(\tfrac12) = -\gamma - 2\ln 2$
조화수	$\psi(n+1) = H_n - \gamma$
양의 근	$x_0 = 1.4616321\ldots$ (감마의 최솟점)
큰 인자	$\psi(x) \sim \ln x - 1/(2x) - 1/(12x^2) + \cdots$
통계에서	$E[\ln X] = \psi(k) + \ln\theta$ ($X \sim$ 감마)

로그를 씌운 확률변수의 기댓값을 물으면 거의 언제나 이 함수가 대답한다.

디감마 함수(digamma function) $\psi(x)$ 는 감마 함수의 로그 도함수

$$\psi(x) = \frac{d}{dx}\ln\Gamma(x) = \frac{\Gamma'(x)}{\Gamma(x)}$$

다. 이름은 그리스 문자 디감마(Ϝ, 감마 두 개를 겹친 모양)에서 왔고, 「감마의 다음」 이라는 말장난이기도 하다¹. 다시 미분하면 트라이감마 $\psi_1 = \psi'$, 일반적으로 폴리감마 $\psi_m = \psi^{(m)}$ 이 된다.

순수수학에서는 조화수의 연속판 정도로 소개되지만, 통계·기계학습 코드에서 디감마는 매일 호출되는 함수다. 감마·베타·디리클레 분포에서 $\ln X$ 의 기댓값이 디감마이기 때문이다. 그 결과 감마 분포의 최대우도추정, 잠재 디리클레 할당(LDA)의 변분 추론, k-최근접 이웃 기반 엔트로피·상호정보량 추정기에 전부 디감마가 들어간다. 이 문서는 성질을 짧게 정리하고 계산법과 그 쓰임새를 실측과 함께 다룬다.

2. 성질

• 점화식. $\Gamma(x+1) = x\Gamma(x)$ 에 로그를 취해 미분하면 $\psi(x+1) = \psi(x) + 1/x$.
• 조화수. $\psi(1) = -\gamma$ ($\gamma = 0.5772156649\ldots$ 는 오일러-마스케로니 상수)이고, 점화식을 반복하면 $\psi(n+1) = H_n - \gamma$, $H_n = 1 + \tfrac12 + \cdots + \tfrac1n$. 그래서 $\psi$ 는 조화수를 실수로 확장한 것이라고 보면 된다.
• 반사 공식. $\psi(1-x) - \psi(x) = \pi\cot(\pi x)$. 음수 인자를 양수로 옮기는 표준 경로다.
• 반정수. $\psi(\tfrac12) = -\gamma - 2\ln 2 = -1.9635100\ldots$
• 극과 근. $x = 0, -1, -2, \dots$ 에서 1차 극(유수 $-1$). 양의 실축에서는 단조증가하고 근이 하나뿐인데, 그 근 $x_0 = 1.4616321449\ldots$ 가 바로 감마 함수의 최솟점이다. 도함수가 0 인 곳이니 당연하다.
• 급수. $\psi(x) = -\gamma + \sum_{n=0}^{\infty}\Bigl(\frac{1}{n+1} - \frac{1}{n+x}\Bigr)$. 수렴이 $1/n$ 급으로 느려서 계산용은 아니다.
• 트라이감마. $\psi_1(x) = \sum_{n\ge0}(x+n)^{-2}$, $\psi_1(1) = \pi^2/6 = \zeta(2)$. 일반적으로 $\psi_m(1) = (-1)^{m+1}m!\,\zeta(m+1)$ 로 리만 제타 함수와 이어진다.

3. 계산 — 밀고 나서 점근 전개

표준 알고리즘은 감마 함수의 「밀고 나서 스털링」 과 같은 구조다. 작은 $x$ 는 점화식으로 큰 쪽으로 밀어 두고

$$\psi(x) = \psi(x+m) - \sum_{j=0}^{m-1}\frac{1}{x+j},$$

큰 인자에서는 베르누이 수가 계수인 점근 급수

$$\psi(x) \sim \ln x - \frac{1}{2x} - \sum_{k\ge1}\frac{B_{2k}}{2k\,x^{2k}} = \ln x - \frac{1}{2x} - \frac{1}{12x^2} + \frac{1}{120x^4} - \frac{1}{252x^6} + \frac{1}{240x^8} - \frac{1}{132x^{10}} + \cdots$$

를 쓴다. 위의 다섯 보정항까지 쓰고 실측한 최대 상대오차(인자 $0.1 \sim 7$):

밀어 올리는 기준	최대 상대오차
$x \ge 6$ 까지	$1.5\times10^{-11}$
$x \ge 10$ 까지	$3.4\times10^{-14}$

다음 항이 $691/(32760x^{12})$ 이므로 $x = 6$ 에서 $10^{-11}$, $x = 10$ 에서 $2\times10^{-14}$ 로 잘린다는 예측과 맞는다. 배정밀도 끝자리까지 원하면 기준을 더 올리거나 항을 하나 더 넣는다. 대신 점화식으로 많이 밀수록 $1/(x+j)$ 를 빼는 횟수가 늘어 반올림이 쌓이니, 기준을 무한정 올리는 것도 답은 아니다.

함정 둘. 첫째, 양의 근 $x_0$ 근처에서는 $\psi$ 자체가 0 에 가까워 상대오차가 폭발한다. 이 구간까지 상대 정밀도를 원하면 $x_0$ 에서의 테일러 전개를 따로 둔다(Boost·Cephes 계열 구현이 그렇게 한다). 둘째, 음수 인자에서 반사 공식의 $\pi\cot(\pi x)$ 는 정수 근처에서 $\sin(\pi x)$ 를 정확히 계산해야 하므로, $x$ 에서 가장 가까운 정수를 먼저 빼고 삼각함수를 부른다.

4. 통계에서의 쓰임

로그 기댓값. 감마 분포 $X \sim \mathrm{Gamma}(k, \theta)$ 이면 $E[\ln X] = \psi(k) + \ln\theta$, $\mathrm{Var}[\ln X] = \psi_1(k)$ 다. 베타 분포 $\mathrm{Beta}(a, b)$ 에서는 $E[\ln X] = \psi(a) - \psi(a+b)$. 디리클레 분포 $\mathrm{Dir}(\boldsymbol\alpha)$ 에서는

$$E[\ln\theta_k] = \psi(\alpha_k) - \psi\Bigl(\sum_j\alpha_j\Bigr)$$

이다. 지수족의 충분통계량이 $\ln x$ 인 분포에서 로그 정규화 상수의 도함수가 감마 함수의 로그 도함수가 되기 때문이다.

감마 분포 최대우도추정. 형상 모수 $k$ 의 MLE 는

$$\ln k - \psi(k) = s,\qquad s = \ln\bar x - \overline{\ln x}$$

의 해다. 닫힌꼴이 없으니 뉴턴-랩슨법으로 풀고, 도함수에 트라이감마 $1/k - \psi_1(k)$ 가 들어간다. 민카(2000)의 초기값

$$k_0 = \frac{3 - s + \sqrt{(s-3)^2 + 24s}}{12s}$$

에서 출발하면 거의 즉시 끝난다. $k = 2.5$, $\theta = 2$ 에서 뽑은 표본 1000개($s = 0.20442$)로 실측하면 $k_0 = 2.59366$, 뉴턴 1회 2.6005400, 2회 2.6005593, 3회에 기계 정밀도로 수렴한다(SciPy gamma.fit 과 일치). 적률법 추정 $\bar x^2/s^2 = 2.608$ 과는 셋째 자리에서 갈린다. 자세한 것은 감마 분포와 최대우도추정에.

LDA 와 변분 추론. 잠재 디리클레 할당(블라이·응·조던 2003)의 변분 추론 갱신식은 $\phi_{nk} \propto \beta_{k,w_n}\exp\bigl(\psi(\gamma_k) - \psi(\sum_j\gamma_j)\bigr)$ 이다. 문서 수 × 단어 수 × 토픽 수만큼 디감마를 부르니, 대규모 토픽 모델에서는 디감마가 프로파일러 상위에 올라온다. 그래서 $x \ge 6$ 기준 + 짧은 점근식 같은 「적당히 정확하고 빠른」 구현이 이 분야에서 흔하다².

디리클레 모수 추정. 디리클레 MLE 고정점 반복은 $\psi(\alpha_k^{\text{new}}) = \psi(\sum_j\alpha_j) + \overline{\ln p_k}$ 형태라 역디감마 $\psi^{-1}(y)$ 가 필요하다. 민카의 처방은 초기값을 $y \ge -2.22$ 이면 $e^y + \tfrac12$, 아니면 $-1/(y + \gamma)$ 로 잡고 뉴턴 5회를 돌리는 것이다. $\psi$ 가 양의 축에서 단조라 해가 유일하다.

kNN 엔트로피·상호정보량 추정기. 코자첸코-레오넨코 엔트로피 추정기와 그 확장인 크라스코프-스퇴그바우어-그라스베르거(KSG) 상호정보량 추정기는

$$\hat I(X;Y) = \psi(k) + \psi(N) - \bigl\langle \psi(n_x + 1) + \psi(n_y + 1)\bigr\rangle$$

꼴이다($k$ 번째 이웃 거리 안에 든 주변 점 개수 $n_x$, $n_y$). 여기서 $\ln$ 대신 $\psi$ 가 들어가는 이유는 $k$ 번째 이웃 거리의 확률질량에 대한 $\ln$ 의 기댓값이 정확히 $\psi(k) - \psi(N)$ 이기 때문이다 — 작은 정수 개수에서 $\ln n$ 으로 바꾸면 편향이 생긴다. 시계열 인과 분석·전달 엔트로피 추정에서 이 꼴을 그대로 쓴다³.

5. 관련 문서

• 감마 함수 · 베타 함수 · 감마 분포 · 베르누이 수
• 점근 급수 · 리만 제타 함수 · 뉴턴-랩슨법
• 최대우도추정 · 변분 추론 · 상호정보량 · 피셔 정보
• 오일러-마스케로니 상수 · 디리클레 분포 · 잠재 디리클레 할당 · 전달 엔트로피

6. Footnotes

1. 디감마(Ϝ)는 고대 그리스 문자에서 사라진 글자다. 수학에서는 살아남았으니 문자 입장에서는 이직 성공. ↩

2. 토픽 모델 학습이 느리다고 GPU 를 질렀는데 프로파일러를 열어 보니 디감마가 1등. 수학 라이브러리 호출 한 줄이 클러스터 비용을 좌우하는 현생. ↩

3. KSG 추정기 논문의 식을 「어차피 $\psi(n) \approx \ln n$ 이니까」 하고 로그로 바꿔 짜면, 표본이 적은 구간에서 상호정보량이 음수로 나오는 킹받는 경험을 하게 된다. ↩