Codes for implementing algorithms mentioned in this chapter are available at my Github repository. Each script runs the respective algorithm with example data.
1. Gamma, Beta, Factorials, Binomial
The gamma function is defined in relation to the factorial function for integer $z$s (see part 1 in my notes on distributions for a simple proof)
\[\Gamma(z) = \int^{\infty}_0 t^{z-1}e^{-t} \mathop{dt} = (z-1)!\]In order to estimate the gamma function outputs, the Lanczos approximation is used, which accounts for all $\mathbb{Re} (z)>0$ (see Appendix A for the proof)
\[\Gamma(z+1) = (z+g+\frac{1}{2})^{z+\frac{1}{2}} e^{-(z+g+\frac{1}{2})} \sqrt{2\pi} (c_0 + \frac{c_1}{z+1} + \frac{c_2}{z+2} + ... + + \frac{c_3}{z+3} + \epsilon)\]As the value of $\Gamma(z)$ can often become overflowingly large, but is only required in intermediate steps, it is better to solve for $\ln \Gamma(z)$ instead in practice. Relatedly, the factorial functions can also be estimated as $exp(\ln \Gamma(z+1))$.
Again, if we can compute $\ln (z!)$, we can continue to estimate the binomial coefficients, which is defined as
\[{n \choose k} = \frac{n!}{k!(n-k)!} \quad 0 \le k \le n\]Finally, computing $\ln \Gamma(z)$ also lets us compute the beta function, which is defined as
\[B(z, x) = \int^1_0 t^{z-1} (1-t)^{x-1} \mathop{dt}\]and can be computed by (see Appendix B for the proof)
\[B(z, x) = \frac{\Gamma(z) \Gamma(x)}{\Gamma(z+x)}\]2. Incomplete Gamma, Error, Chi-Square, Cumlative-Poisson
The incomplete gamma function is bound in the range $[0, 1]$ with different slopes depending on the value of $x$. It is defined by
\[P(z, x) = \frac{\gamma(z, x)}{\Gamma(z)} = \frac{1}{\Gamma(z)} \int^{x}_0 t^{z-1}e^{-t} \mathop{dt} \quad (a>0) \\\]where $\Gamma(z)$ is the gamma function from the previous section. We can approximate $\gamma(z, x)$ by
\[\gamma(z, x) = \int^x_0 t^{z-1}e^{-t} \mathop{dt} = e^{-x} x^{z} \sum^{\infty}_{t=0} \frac{\Gamma(z)}{\Gamma(z+1+t)} x^t\]The incomplete gamma function’s complement is also sometimes called an incomplete gamma function
\[Q(z, x) = 1 - P(z, x) = \frac{\Gamma(z, x)}{\Gamma(z)} = \frac{1}{\Gamma(z)} \int^{\infty}_x t^{z-1}e^{-t} \mathop{dt} \quad (a>0)\]3. Incomplete Beta, Student’s F, Cumulative Binomial
Appendix
A. The Lanczos approximation1
This approximation can be derived with a series of substitutions. First, let $t = (z+g+1)(1-\ln v)$, where $g$ is an arbitrary constant. Then we know the following for applying substitution:
\[v = e^{1 - \frac{t}{z+g+1}}\\ \frac{\mathop{dt}}{\mathop{dv}} = -\frac{z+g+1}{v} \\ \text{When } t = \infty \text{, } v = e^{- \infty} = 0 \\ \text{When } t = 0 \text{, } v = e^{1 - 0} = e \\\]Now we can proceed with the substitution:
\[\begin{aligned} \Gamma(z+1) = z! &= \int^{\infty}_0 t^{z}e^{-t} \mathop{dt} =\int^0_e [(z+g+1)(1-\ln v)]^z e^{-(z+g+1)(1-\ln v)} (-\frac{z+g+1}{v}) {\mathop{dv}} \\ &= \int^e_0 (z+g+1)^{z+1} e^{-(z+g+1)}(1-\ln v)^z e^{(z+g+1)(\ln v)} \frac{1}{v} {\mathop{dv}} \\ &= (z+g+1)^{z+1} e^{-(z+g+1)} \int^e_0 (1-\ln v)^z v^{z+g} {\mathop{dv}} \\ &= (z+g+1)^{z+1} e^{-(z+g+1)} \int^e_0 [v(1-\ln v)]^z v^{g} {\mathop{dv}} \end{aligned}\]To make the later steps easier, instead of $\Gamma(z+1)$, we will actually compute $\Gamma(z+\frac{1}{2})$:
\[\Gamma(z+\frac{1}{2}) = (z-\frac{1}{2})! = (z+g+\frac{1}{2})^{z+\frac{1}{2}} e^{-(z+g+\frac{1}{2})} \int^e_0 [v(1-\ln v)]^{(z-\frac{1}{2})} v^{g} {\mathop{dv}}\]The first two terms look in place now; so we focus on the last integral term. The second substitution would be $v(1-\ln v) = 1 - x^2$. The intuition is that the function $v(1-\ln v)$ is concave and evaluates to $0$ both when $v=0$ and when $v=e$. This shape is similar to the inverted quadratic curve of $1-x^2$, which evaluates to $0$ both when $x=-1$ and when $x=1$. Additionally, we can also correspond $v=1$ to $x=0$.
For applying substitution, we need to know $\frac{\mathop{dv}}{\mathop{dx}}$:
\[\begin{aligned} \frac{\mathop{d}}{\mathop{dx}} v(1-\ln v) &= \frac{\mathop{d}}{\mathop{dx}} 1 - x^2 \\ \frac{\mathop{dv}}{\mathop{dx}} \frac{\mathop{d}}{\mathop{dv}} v(1-\ln v) &= -2x \\ \frac{\mathop{dv}}{\mathop{dx}} &= \frac{-2x}{v(-\frac{1}{v}) + (1-\ln v)} \\ \frac{\mathop{dv}}{\mathop{dx}} &= \frac{2x}{\ln v} \\ \end{aligned}\]We can also rearrange the equation to get $\ln v = 1 - \frac{1-x^2}{v}$. Combining the two equations gives us:
\[\begin{aligned} \frac{\mathop{dv}}{\mathop{dx}} &= \frac{2x}{1 - \frac{1-x^2}{v}} \\ \frac{\mathop{dv}}{\mathop{dx}} - (\frac{\mathop{dv}}{\mathop{dx}}) (\frac{1-x^2}{v}) - 2x &= 0 \\ (\frac{\mathop{dv}}{\mathop{dx}}) v - (\frac{\mathop{dv}}{\mathop{dx}}) (1-x^2) - 2xv &= 0 \end{aligned}\]We can approximate $v(x)$ by computing its Talyor series around $x=0$:
\[v(x) = \sum^{\infty}_{n=0} \frac{v^{(n)}(0)}{n!} x^n = 1 + v'(0)x + \frac{v''(0)}{2!} x^2 + \frac{v'''(0)}{3!}x^3 + ... \\\]Putting this approximation back into the combined euqation:
\[\begin{aligned} (\frac{\mathop{d}}{\mathop{dx}} 1 + v'(0)x + \frac{v''(0)}{2!} x^2 + \frac{v'''(0)}{3!}x^3 + ...) (1 + v'(0)x + \frac{v''(0)}{2!} x^2 + \frac{v'''(0)}{3!}x^3 + ...) - (\frac{\mathop{d}}{\mathop{dx}} 1 + v'(0)x + \frac{v''(0)}{2!} x^2 + \frac{v'''(0)}{3!}x^3 + ...) (1-x^2) - 2x (1 + v'(0)x + \frac{v''(0)}{2!} x^2 + \frac{v'''(0)}{3!}x^3 + ...) &= 0 \\ (v'(0) + v''(0) x + \frac{v'''(0)}{2}x^2 + ...) (1 + v'(0)x + \frac{v''(0)}{2!} x^2 + \frac{v'''(0)}{3!}x^3 + ...) - (v'(0) + v''(0) x + \frac{v'''(0)}{2}x^2 + ...) (1-x^2) - 2x (1 + v'(0)x + \frac{v''(0)}{2!} x^2 + \frac{v'''(0)}{3!}x^3 + ...) &= 0 \\ v'(0) (1+x) + v''(0) (x + \frac{x^2}{2}) + v'(0)v''(0) \frac{3}{2} x^2 + ... &= 0 \end{aligned}\]END: \(\Gamma(z+1) = (z+g+\frac{1}{2})^{z+\frac{1}{2}} e^{-(z+g+\frac{1}{2})} \sqrt{2\pi} (c_0 + \frac{c_1}{z+1} + \frac{c_2}{z+2} + ... + + \frac{c_3}{z+3} + \epsilon)\)
B. Relationship between beta and gamma functions
We start with the numerator of the right-hand side
\[\begin{aligned} \Gamma(z) \Gamma(x) &= \int^{\infty}_0 t^{z-1}e^{-t} \mathop{dt} \int^{\infty}_0 s^{x-1}e^{-s} \mathop{ds} \\ &= \int^{\infty}_0 \int^{\infty}_0 t^{z-1}s^{x-1}e^{-t-s} \mathop{ds} \mathop{dt} \\ \end{aligned}\]Applying a change of variable, where $t=uv$ and $s=u(1-v)$, we obtain the new differentials \(\mathop{ds} \mathop{dt} = |\det (D\phi)(u,v)|\mathop{du} \mathop{dv}\)
where $(D\phi)(u,v)$ is the Jacobian matrix of partial derivatives of $(s,t) = \phi (u, v)$ at point $(u,v)$. Hence
\[\begin{aligned} \mathop{ds} \mathop{dt} &= |\det \begin{bmatrix} \frac{\partial s}{\partial u} & \frac{\partial s}{\partial v} \\ \frac{\partial t}{\partial u} & \frac{\partial t}{\partial v} \end{bmatrix} | \mathop{du} \mathop{dv} \\ &=|\det \begin{bmatrix} 1-v & -u \\ v & u \end{bmatrix} | \mathop{du} \mathop{dv} \\ &= |u(1-v)-(-u)v| \mathop{du} \mathop{dv} \\ &= u \mathop{du} \mathop{dv} \end{aligned}\]Now we have
\[\begin{aligned} \Gamma(z) \Gamma(x) &= \int^{\infty}_0 \int^1_0 (uv)^{z-1}[u(1-v)]^{x-1}e^{-uv-u(1-v)} u \mathop{dv} \mathop{du}\\ &= \int^{\infty}_0 u^{z+x-1} e^{-u} \mathop{du} \int^1_0 v^{z-1} (1-v)^{x-1} \mathop{dv} \\ &= \Gamma(z+x)B(z,x) \end{aligned}\]-
Lanczos, C. (1964). A Precision Approximation of the Gamma Function. Journal of the Society for Industrial and Applied Mathematics Series B Numerical Analysis, 1(1), 86–96. doi:10.1137/0701008 ↩