베타 함수

1. 개요

베타 함수 Beta Function
정의	$B(a,b) = \int_0^1 t^{a-1}(1-t)^{b-1}\,dt$
감마 표현	$B(a,b) = \Gamma(a)\Gamma(b)/\Gamma(a+b)$
대칭	$B(a,b) = B(b,a)$
특수값	$B(\tfrac12,\tfrac12) = \pi$, $B(1,b) = 1/b$
불완전판	$I_x(a,b)$ = 베타 분포의 누적분포함수
통계 연결	이항·스튜던트 t·F 분포의 누적확률
표준 구현	TOMS 708 (DiDonato-Morris 1992), 연분수 + 대칭 전환

이항계수의 연속판. 통계 소프트웨어의 p-값 절반은 결국 이 함수의 불완전판을 부른다.

베타 함수(beta function) $B(a, b)$ 는 양의 실수부를 갖는 $a$, $b$ 에 대해

$$B(a, b) = \int_0^1 t^{a-1}(1-t)^{b-1}\,dt = \frac{\Gamma(a)\,\Gamma(b)}{\Gamma(a+b)}$$

로 정의되는 특수함수다. 정수 인자에서는 $B(m, n) = (m-1)!(n-1)!/(m+n-1)!$ 로, 이항계수의 역수와 거의 같다. 오일러가 1730년대에 감마 함수와 함께 다뤘고, 「제1종 오일러 적분」 이라고도 부른다(감마 함수가 제2종)¹.

감마 함수 문서가 계산법의 큰 틀을 다뤘으니, 이 문서는 불완전 베타 함수 $I_x(a, b)$ 와 그 계산에 집중한다. 실무에서 베타 함수 자체보다 불완전판을 훨씬 자주 부르기 때문이다. 이항검정의 p-값, 스튜던트 t 검정, F 검정, 베이지안 A/B 테스트의 사후확률이 모두 $I_x$ 한 함수로 귀결된다.

2. 성질

• 삼각함수 꼴. $t = \sin^2\theta$ 로 치환하면 $B(a, b) = 2\int_0^{\pi/2}\sin^{2a-1}\theta\,\cos^{2b-1}\theta\,d\theta$. 그래서 $B(\tfrac12, \tfrac12) = \pi$ 이고, 이것이 $\Gamma(\tfrac12)^2 = \pi$ 의 다른 증명이다.
• 점화식. $B(a+1, b) = B(a, b)\cdot a/(a+b)$. 정수 쪽으로 밀거나 당길 때 쓴다.
• 이항계수. $\binom{n}{k} = 1/\bigl((n+1)B(n-k+1, k+1)\bigr)$. 정수가 아닌 $n$, $k$ 로 이항계수를 늘리는 표준 경로다.
• 다변수판. $B(\boldsymbol\alpha) = \prod_i\Gamma(\alpha_i)/\Gamma(\sum_i\alpha_i)$ 는 디리클레 분포의 정규화 상수다. 심플렉스 위에서 $\prod x_i^{\alpha_i - 1}$ 을 적분한 값이다.

3. 불완전 베타 함수와 통계 분포

정칙화 불완전 베타 함수는

$$I_x(a, b) = \frac{1}{B(a,b)}\int_0^x t^{a-1}(1-t)^{b-1}\,dt,\qquad 0 \le x \le 1$$

이다. $I_0 = 0$, $I_1 = 1$ 이고 대칭식 $I_x(a, b) = 1 - I_{1-x}(b, a)$ 가 성립한다. 이것이 곧 베타 분포 $\mathrm{Beta}(a, b)$ 의 누적분포함수다. 베타 분포의 평균은 $a/(a+b)$, 분산은 $ab/((a+b)^2(a+b+1))$ 이다.

여러 분포의 누적확률이 $I_x$ 로 바뀐다. 모두 수치로 확인했다.

분포	누적확률	확인 (예)
이항 $\mathrm{Bin}(n, p)$	$P(X \le k) = I_{1-p}(n-k,\ k+1)$	$n=20, k=5, p=0.3$ → 0.416371
스튜던트 t ($\nu$ 자유도, $t > 0$)	$F(t) = 1 - \tfrac12 I_{\nu/(\nu+t^2)}(\nu/2,\ 1/2)$	$\nu=5, t=2$ → 0.949030
F 분포 $(d_1, d_2)$	$F(x) = I_{d_1x/(d_1x+d_2)}(d_1/2,\ d_2/2)$	—
음이항	$P(X \le k) = I_p(r,\ k+1)$	—

이항 누적확률을 $\sum_{j\le k}\binom{n}{j}p^j(1-p)^{n-j}$ 로 직접 더하면 $n$ 이 클 때 항이 넘치거나 언더플로하고 $O(k)$ 번의 연산이 든다. $I_x$ 한 번이면 $n = 10^9$ 이어도 상수 시간이다.

베이지안 A/B 테스트. 베타 분포는 이항 우도의 켤레사전분포라, 균등 사전 $\mathrm{Beta}(1,1)$ 에서 1000명 중 120명 전환을 관측하면 사후분포가 $\mathrm{Beta}(121, 881)$ 이다. 사후 평균 0.1208, 95% 신용구간 $[0.1013,\ 0.1416]$ 은 $I_x$ 의 역함수로 구한다. B 안이 1000명 중 140명이면 $P(p_B > p_A) = 0.908$ 이다. 이 확률은 베타 함수 비의 유한합으로 닫힌꼴 계산이 가능하고(백만 쌍 몬테카를로 0.9075 와 일치), 표본이 작을 때 정규 근사보다 믿을 만하다².

4. 어떻게 계산하나

4.1. 완전 베타 — 로그 감마의 소거

$\ln B(a, b) = \ln\Gamma(a) + \ln\Gamma(b) - \ln\Gamma(a+b)$ 는 공식으로는 끝이지만, $a \gg b$ 이면 거의 같은 두 큰 수를 뺀다. 실측하면

$a$	$b$	순진한 lgamma 조합의 절대오차	정답 $\ln B$
$10^6$	2.5	$1.9\times10^{-9}$	−34.254095
$10^{10}$	2.5	$2.4\times10^{-5}$	−57.279944

$\ln\Gamma(10^{10}) \approx 2.2\times10^{11}$ 이라 배정밀도 반올림 자체가 $10^{-5}$ 급이고, 그 오차가 결과에 그대로 남는다(부동소수점 연산). $a \gg b$ 에서는 $\ln B(a, b) \approx \ln\Gamma(b) - b\ln a$ 같은 점근식을 쓰거나, 스털링 보정항끼리 따로 빼는 구현(SciPy betaln, Boost beta)을 쓴다. 이 구간에서 자작 lgamma 조합은 금지다.

4.2. 불완전 베타 — 연분수와 대칭 전환

$I_x(a, b)$ 의 표준 수치 경로는 연분수다(DLMF 8.17.22).

$$I_x(a,b) = \frac{x^a(1-x)^b}{a\,B(a,b)}\cdot\cfrac{1}{1+\cfrac{d_1}{1+\cfrac{d_2}{1+\cdots}}}$$

계수 $d_{2m+1} = -\frac{(a+m)(a+b+m)x}{(a+2m)(a+2m+1)}$, $d_{2m} = \frac{m(b-m)x}{(a+2m-1)(a+2m)}$ 이다. 이 연분수는 $x < (a+1)/(a+b+2)$ 일 때 빨리 수렴하고, 그 반대면 대칭식으로 $I_{1-x}(b, a)$ 를 계산해 1 에서 뺀다. 평가는 렌츠(1976)의 앞에서-뒤로 알고리즘(톰슨-바넷 수정판)으로 하면 미리 깊이를 정하지 않고 수렴할 때 멈출 수 있다.

실측한 반복 횟수(상대 증분 $10^{-15}$ 에서 정지):

$(a, b, x)$	경계 $(a+1)/(a+b+2)$	반복	상대오차
(2, 3, 0.3)	0.43	3	$10^{-16}$
(50, 50, 0.45)	0.50	16	$2\times10^{-16}$
(50, 50, 0.55) — 경계 넘음	0.50	24	$1.5\times10^{-14}$
(500, 500, 0.49)	0.50	37	$5.5\times10^{-13}$

경계를 넘겨도 수렴은 하지만 반복이 늘고 정밀도가 깎인다. 큰 $a, b$ 에서 정밀도가 떨어지는 원인은 연분수가 아니라 앞의 전치인자 $x^a(1-x)^b/B(a,b)$ 다 — 지수가 수백인 로그들을 더하고 빼는 과정에서 자릿수가 빠진다. $a$, $b$ 가 모두 크면 $I_x$ 는 정규분포 쪽으로 다가가므로, 전용 균등 점근 전개(템메의 전개)를 쓰는 것이 낫다³.

TOMS 708(DiDonato-Morris 1992)은 이런 구간을 $(a, b, x)$ 공간에서 십여 개로 나눠 거듭제곱 급수·연분수·점근 전개를 골라 쓰는 알고리즘이고, R 의 pbeta 가 지금도 이것을 쓴다. 꼬리 확률이 $10^{-300}$ 근처여도 상대 정밀도를 지키도록 설계돼 있다. 자작 연분수 구현이 「대부분 맞는데 가끔 틀리는」 구간을 다 막아 둔 결과물이라, 직접 구현할 이유는 GPU 커널처럼 라이브러리가 없는 곳뿐이다.

역함수. 신용구간·분위수에는 $I_x(a,b) = p$ 를 $x$ 에 대해 풀어야 한다. 도함수가 베타 밀도 그 자체라 뉴턴-랩슨법이 자연스럽지만, 밀도가 0 이나 무한대로 가는 양 끝에서 튀므로 초기값을 정규·윌슨-힐퍼티 근사로 잡고 $[0, 1]$ 밖으로 나가면 이분법으로 되돌리는 안전장치를 둔다.

5. 관련 문서

• 감마 함수 · 디감마 함수 · 감마 분포 · 연분수
• 켤레사전분포 · 최대우도추정 · 몬테카를로 방법
• 부동소수점 연산 · 뉴턴-랩슨법 · 점근 급수
• 베타 분포 · 디리클레 분포 · 스튜던트 t 분포 · 이항분포

6. Footnotes

1. 「제1종」「제2종」 이름은 르장드르가 붙였다. 감마 함수의 기호 Γ 도 르장드르 작품이다. 오일러가 만들고 르장드르가 이름 붙이는 게 18~19세기 해석학의 국룰. ↩

2. 「p-값 0.05 넘었으니 효과 없음」 대신 「B 가 더 나을 확률 91%」 라고 말할 수 있어서 기획자들이 베이지안 A/B 를 좋아한다. 사전분포를 뭘로 잡았는지 물어보면 대화가 길어진다. ↩

3. 1968년 베네치아노가 강입자 산란 진폭을 맞추려고 쓴 식이 그냥 $B(-\alpha(s), -\alpha(t))$ 였고, 이걸 해석하다가 끈 이론이 태어났다. 특수함수 하나 잘못 건드리면 물리학 한 분야가 생긴다. ↩