神奇的Gamma函数（scipy）

$$\Gamma(x)=\int_0^\infty t^{x-1}e^{-t}dt$$
对应于scipy（python库）的：

from scipy.special import gamma

通过分布积分的方法，进行如下的推导：

$$\Gamma(x+1)=\int_0^\infty t^xe^{-t}dt=-\int_0^\infty t^xd(e^{-t})=-\left [ t^xe^{-t}|_0^\infty-x\int_0^\infty t^{x-1} e^{-t}dt\right ]=x\Gamma(x)$$

可得该函数如下的递归性质：

$$\Gamma(x+1)=x\Gamma(x)$$

>>> gamma(5+1)
120.0
>>> 5*gamma(5)
120.0

于是很容易证明（$\Gamma(1)=1,\;\Gamma(n)=(n-1)\Gamma(n-1)$），$\Gamma(x)$ 函数可以看成是阶乘在实数集上的延拓，具有如下的性质：

$$\Gamma(n)=(n-1)!$$

$\Gamma(n)$ 与 $B(m,n)$的关系

$$B(m,n)=\int_0^1x^{m-1}(1-x)^{n-1}dx\\ B(m,n)=\frac{\Gamma(m)\Gamma(n)}{\Gamma(m+n)}$$

从二项分布到 Gamma 分布

对Gamma 函数的定义稍作变形，可得如下等式：

$$\int_0^{\infty}\frac{x^{\alpha-1}e^{-x}}{\Gamma(\alpha)}dx=1$$
于是 取积分中的函数作为概率密度（Probability Density Function，PDF），就得到一个形式最为简单的Gamma 分布的密度函数（density function）：
$$Gamma(x|\alpha)=\frac{x^{\alpha-1}e^{-x}}{\Gamma(\alpha)}$$
如果再做一个变换 $x=\beta t$，就得Gamma 分布的更一般形式：
$$Gamma(t|\alpha, \beta)=\frac{\beta^\alpha t^{\alpha-1}e^{-\beta t}}{\Gamma(\alpha)}$$
其中 $\alpha$ 称为形状参数（shape parameter），主要决定了分布曲线的形状，而 $\beta$ 称为 rate parameter（或叫 inverse scale parameter， $\frac1\beta$ 称为scale parameter），主要决定曲线有多陡。

Gamma 分布与Possion 分布

Gamma 分布首先与Possion 分布（离散型）、Possion 过程发生密切的联系。我们容易发现Gamma分布的概率密度和Possion分布在数学形式具有高度的一致性。参数为 $\lambda$ 的Possion分布，其概率（probability mass function，pmf）为：

$$Possion(X=k|\lambda)=\frac{\lambda^ke^{-\lambda}}{k!}$$

而Gamma分布的密度（$\alpha=k+1,\;\beta=1$）得到：

$$Gamma(x|\alpha=k+1)=\frac{x^\alpha e^{-x}}{\Gamma(k+1)}=\frac{x^ke^{-x}}{k!}$$
所以这两个分布在数学形式上是一致的，只是Possion分布是离散的，Gamma分布是连续的， 可以直观地认为Gamma分布是Possion分布在正实数集上的连续化版本。

我们在概率论与数理统计的课程中都学过，$Poisson(\lambda)$ 分布可以看成是二项分布 $B(n,p)$ 在 $np=\lambda,n\to\infty$ 条件下的极限分布：

$$B(k;n,p)=\binom{n}{k}p^k(1-p)^{n-k}\quad\stackrel{np=\lambda,n\to \infty}{\longrightarrow} \quad Poisson(X=k|np)=\frac{\lambda^ke^{-\lambda}}{k!}$$
如果你对二项分布的关注足够多，可能会知道二项分布的随机变量 $X~B(n,p)$满足下面一个奇妙的 恒等式：
$$P(X\leq k)=\frac{n!}{k!(n-k-1)!}\int_p^1t^k(1-t)^{n-k-1}dt$$
我们在等式右侧做一个变换 $t=\frac xn$可得：
$$\begin{split} P(X\leq k)=&\frac {n!}{k!(n-k-1)!}\int_{np}^n(\frac xn)^k(1-\frac xn)^{n-k-1}d(\frac xn)\\ =&\int_{np}^n\frac{(n-1)!}{k!(n-k-1)!}(\frac xn)^k(1-\frac xn)^{n-k-1}dx\\ =&\int_{np}^n\binom{n-1}{k}(\frac xn)^k(1-\frac xn)^{n-k-1}dx\\ =&\int_{np}^nB(k;n-1,\frac xn)dx\\ =&\int_{np}^nB(Y=k|n-1,\frac xn)dx \end{split}$$

上式左侧是二项分布 $B(n,p)$ 的累积分布函数（cumulative density function，cdf），而右侧为无穷多个二项分布 $B(n-1,\frac xn)$ 的积分和，所以可以写为：

$$B(X\leq k|n,p)=\int_{np}^nB(Y=k|n-1,\frac xn)dx$$
对上式两边在条件 $np=\lambda,\;n\to \infty$ 下取极限，则得到：
$$Poisson(X\leq k|\lambda)=\int_\lambda^\infty Poisson(Y=k|x)dx=\int_\lambda^\infty\frac{x^ke^{-x}}{k!}dx$$

到这，不妨先暂停，我们使用scipy做一个简单验证：

import scipy.stats as st
from scipy.misc import factorial
from scipy import integrate

lmbda, k = 2, 6
X = st.poisson(2)
                    # X ~ Poisson(2)
print(X.cdf(k))
                    # P(X<=k)

def poisson_pdf(x, k):
    return x**k*np.exp(-x)/factorial(k)
print(integrate.quad(possion_pdf, lmbda, np.inf, args=(k))[0])
                # 0.995466194474
                # 0.9954661944737513
                # 两者达到完美的相等

书归正传，我们再来看上面的公式，

$$Poisson(X\leq k|\lambda)=\int_\lambda^\infty Poisson(Y=k|x)dx=\int_\lambda^\infty \frac{x^ke^{-x}}{k!}dx$$

该等式即为著名的 Poisson-Gamma duality，接下来我们来点好玩的，对上面的等式两边在 $\lambda \to 0$取极限，左侧Poisson分布表示的是至少发生 $k$事件的概率，$\lambda \to 0$的时候就不可能有事件再发生了，故 $P(X\leq k)=0$，于是：

$$1=\lim_{\lambda \to 0}\int_\lambda^\infty \frac{x^ke^{-x}}{k!}dx=\int_{0}^\infty \frac{x^ke^{-x}}{k!}dx$$
该积分式子说明 $f(x)=\frac{x^ke^{-x}}{k!}$在实数集上是一个概率分布函数（probability density function，pdf），而这个函数恰好就是 Gamma分布（ $Gamma(x|\alpha)=\frac{x^{\alpha-1}e^{-x}}{\Gamma(\alpha)},\alpha=k+1$）（ 也即我们通过二项分布，再根据泊松定理，推导出了最后的Gamma分布）。我们继续把上式右边中的 $k!$移到左边，于是得到：

$$k!=\int_0^\infty x^ke^{-x}dx$$
于是我们得到了将 $k!$表示为积分（也即从此我们可以计算3.5的阶乘了）的方法。

好吧，文章已经够长的了，公式已经够繁复的了，下文，我们将重点实现阶乘在正实数轴上的连续化（或者叫插值），其实已经很明显了。