​对伽玛函数的定义做一个变形,就可以得到如下等式

\int_{0}^{\infty} \frac{x^{\alpha-1} e^{-x}}{\Gamma(\alpha)} dx=1 \\

​于是,取积分中的函数作为概率密度,就得到一个形式最简单的伽玛分布的密度函数

\operatorname{Gamma}(x \mid \alpha)=\frac{x^{\alpha-1} e^{-x}}{\Gamma(\alpha)} \\

​如果做一个变换 x=\beta t , 就得到伽玛分布的更一般的形式

\operatorname{Gamma}(t \mid \alpha, \beta)=\frac{\beta^{\alpha} t^{\alpha-1} e^{-\beta t}}{\Gamma(\alpha)} . \\

​其中 \alpha 被称为形状参数(shape parameter),主要决定了分布曲线的形状;\beta 称为率参数/逆尺度参数( rate parameter / inverse scale parameter)(\frac{1}{\beta}被称为尺度参数)。

​这里尺度参数的定义和人大那本《数据科学统计基础》不一样,在该书的第52页中,伽马分布的定义如下:若X\sim Ga(\alpha,\lambda),密度函数为 p(x)=\frac{\lambda^{\alpha} }{\Gamma(\alpha)}x^{\alpha-1} e^{-\lambda x},x\geq 0,则\alpha为形状参数,\lambda为尺度参数。

Gamma分布不同\alpha值下的图形比较(\lambda=1)

形状参数\alpha决定了分布曲线的形状吗,随着\alpha的变化,伽马分布的图像的单调性,凹凸性发生了变化。

Gamma分布不同\lambda值下图形比较(\alpha=10)

尺度参数\lambda决定了曲线有多陡,\lambda取值越小,伽马分布越陡。

这里我们采用《数据科学统计基础》\Gamma分布密度函数的定义:若X\sim Ga(\alpha,\lambda),密度函数为 p(x)=\frac{\lambda^{\alpha} }{\Gamma(\alpha)}x^{\alpha-1} e^{-\lambda x},x\geq 0,分布函数为F(x)=\int_0^{x}\frac{\lambda^{\alpha} }{\Gamma(\alpha)}t^{\alpha-1} e^{-\lambda t},x\geq 0.

1.期望与方差

E(X)=\frac{\alpha}{\lambda},\quad Var(X)=\frac{\alpha}{\lambda^2} \\

证明:由密度函数积分等于1可以得出

\int_0^{\infty}x^{\alpha-1} e^{-\lambda x}=\frac{\Gamma(\alpha)}{\lambda^{\alpha}}. \\

因此,

\begin{align} E(X)=\int_{0}^{\infty}xp(x)\mathrm{~d}x &=\int_{0}^{\infty}\frac{\lambda^{\alpha}}{\Gamma(\alpha)}x^{\alpha} e^{-\lambda x} \mathrm{~d}x \\ &=\frac{\lambda^{\alpha}}{\Gamma(\alpha)}\frac{\Gamma(\alpha+1)}{\lambda^{\alpha+1}}\int_{0}^{\infty}\frac{\lambda^{\alpha+1}}{\Gamma(\alpha+1)}x^{(\alpha+1)-1} e^{-\lambda x} \mathrm{~d}x \\ &=\frac{\lambda^{\alpha}}{\Gamma(\alpha)}\frac{\alpha\Gamma(\alpha)}{\lambda^{\alpha+1}}\times1\\ &=\frac{\alpha}{\lambda}. \end{align} \\\begin{align} Var(X)=E(X^2)-[E(X)]^2 &=\int_{0}^{\infty}x^2p(x)\mathrm{~d}x-(\frac{\alpha}{\lambda})^2 \mathrm{~d}x\\ &=\int_{0}^{\infty}\frac{\lambda^{\alpha}}{\Gamma(\alpha)}x^{\alpha+1}e^{-\lambda x} \mathrm{~d}x-(\frac{\alpha}{\lambda})^2 \mathrm{~d}x\\ &=\frac{\lambda^{\alpha}}{\Gamma(\alpha)}\frac{\Gamma(\alpha+2)}{\lambda^{\alpha+2}}{\lambda^{\alpha+1}}\int_{0}^{\infty}\frac{\lambda^{\alpha+2}}{\Gamma(\alpha+2)}x^{(\alpha+2)-1} e^{-\lambda x} \mathrm{~d}x -(\frac{\alpha}{\lambda})^2 \\ &=\frac{\lambda^{\alpha}}{\Gamma(\alpha)}\frac{(\alpha+1)\alpha\Gamma(\alpha)}{\lambda^{\alpha+2}}-\frac{\alpha^2}{\lambda^2}\\ &=\frac{\alpha}{\lambda^2}. \end{align} \\

\begin{align} E(X^k)=\int_{0}^{\infty}x^kp(x)\mathrm{~d}x &=\int_{0}^{\infty}\frac{\lambda^{\alpha}}{\Gamma(\alpha)}x^{\alpha+k-1} e^{-\lambda x} \mathrm{~d}x\\ &=\frac{\lambda^{\alpha}}{\Gamma(\alpha)}\frac{\Gamma(\alpha+k)}{\lambda^{\alpha+k}}\\ &=\frac{\Gamma(\alpha+k)}{\lambda^{k}\Gamma(\alpha)} \end{align} \\

3.偏度\beta_s,峰度\beta_k

\beta_s=\frac{E(X-E(X))^3}{[E(X-E(X))^2]^{3/2}}=\frac{2}{\sqrt{\alpha}}, \quad \beta_k=\frac{E(X-E(X))^4}{[E(X-E(X))^2]^{2}}-3=\frac{6}{\alpha} \\

证明:由于E(X^4)=\frac{\Gamma(\alpha+4)}{\lambda^{4}\Gamma(\alpha)},E(X^3)=\frac{\Gamma(\alpha+3)}{\lambda^{3}\Gamma(\alpha)},E(X^2)=\frac{\Gamma(\alpha+2)}{\lambda^{2}\Gamma(\alpha)},

记\mu=E(X),\sigma^2=Var(X),因此,

\begin{align} \beta_s=\frac{E(X-E(X))^3}{[E(X-E(X))^2]^{3/2}}= \frac{E(X^{3})-3 \mu \sigma^{2}-\mu^{3}}{\sigma^{3}}=& \frac{\frac{\alpha(\alpha+1)(\alpha+2)}{\lambda^{3}}-3\left(\frac{\alpha}{\lambda}\right)\left(\frac{\alpha}{\lambda^{2}}\right)-\left(\frac{\alpha}{\lambda}\right)^{3}}{\left(\frac{\alpha \sqrt{\alpha}}{\lambda^{3}}\right)} \\ =&\frac{\alpha^{3}+3 \alpha^{2}+2 \alpha-3 \alpha^{2}-\alpha^{3}}{(\alpha \sqrt{\alpha})} \\ =& \frac{2}{\sqrt{\alpha}} \end{align} \\\begin{aligned} \beta_k=\frac{E(X-E(X))^4}{[E(X-E(X))^2]^{2}}-3&= \frac{\mathrm{E}\left(X^{4}\right)-4 \mu \mathrm{E}\left(X^{3}\right)+6 \mu^{2} \mathrm{E}\left(X^{2}\right)-3 \mu^{4}}{\sigma^{4}}-3\\ &=\frac{\frac{\alpha(\alpha+1)(\alpha+2)(\alpha+3)}{\lambda^{4}}-4\left(\frac{\alpha}{\lambda}\right)\left(\frac{\alpha(\alpha+1)(\alpha+2)}{\lambda^{3}}\right)+6\left(\frac{\alpha}{\lambda}\right)^{2}\left(\frac{\alpha(\alpha+1)}{\lambda^{2}}\right)-3\left(\frac{\alpha}{\lambda}\right)^{4}}{\left(\frac{\alpha^{2}}{\lambda^{4}}\right)}-3 \\ &=\frac{\left(\alpha^{4}+6 \alpha^{3}+11 \alpha^{2}+6 \alpha\right)-4\left(\alpha^{4}+3 \alpha^{3}+2 \alpha^{2}\right)+6\left(\alpha^{4}+\alpha^{3}\right)-3 \alpha^{4}}{\left(\alpha^{2}\right)}-3 \\ &=\frac{(1-4+6-3) \alpha^{4}+(6-12+6) \alpha^{3}+(11-8) \alpha^{2}+6 \alpha}{\left(\alpha^{2}\right)}-3 \\ &=\frac{6}{\alpha} \end{aligned} \\

偏度是描述分布偏离对称性程度的一个特征数。偏度的衡量是相对于正态分布来说,正态分布的偏度为0,即若数据分布是对称的,偏度为0。若偏度大于0,则分布右偏,即分布有一条长尾在右;若偏度小于0,则分布为左偏,即分布有一条长尾在左;同时偏度的绝对值越大,说明分布的偏移程度越严重。

峰度是描述分布尖峭程度和(或)尾部粗细的一个特征数。峰度的取值范围为[1,+∞),完全服从正态分布的数据的峰度值为 3,峰度值越大,概率分布图越高尖,峰度值越小,越矮胖。

由上述计算结果可以得知,伽马分布的偏度和峰度只与\alpha有关,与\lambda无关,这也是\alpha称为形状参数的原因,随着\alpha的增大,\beta_s与\beta_k越来越小,最后趋于正态分布的状态:\beta_s=0,\beta_k=0.

4.矩母函数(mgf)

M_x(t)=E(e^{xt})=\frac{\lambda^{\alpha}}{\Gamma(\alpha)} \int_{0}^{\infty} e^{x t} x^{\alpha-1} e^{-\lambda x} d x=\left(\frac{\lambda}{\lambda-t}\right)^{\alpha}=(1-\frac{t}{\lambda})^{-\alpha} \\

5.特征函数

\varphi_{X}(t)=E(e^{i t X})=(1-\frac{it}{\lambda})^{-\alpha} \\

6.可加性:设X_1\sim Ga(\alpha_1,\lambda),X_2\sim Ga(\alpha_2,\lambda),且X_1与X_2​独立,则

X_1+X_2\sim Ga(\alpha_1+\alpha_2,\lambda) \\

证明:伽马分布Ga(\alpha_i,\lambda)的特征函数为:

\varphi_{k}(t)=(1-\frac{it}{\lambda})^{-\alpha},\quad k=1,2 \\

从而\lambda相同的伽马变量和的特征函数为:

\varphi(t)=\varphi_{1}(t)\varphi_{2}(t)=\varphi_{k}(t)=(1-\frac{it}{\lambda})^{-(\alpha_1+\alpha_2)} \\

这个结果可以推广到n个独立、具有相同尺度参数的伽玛随机变量:

\forall k \in\{1, \ldots, n\}: X_{k} \sim Ga\left(\alpha_{k}, \lambda\right) \quad \Rightarrow \quad \sum_{k=1}^{n} X_{k} \sim Ga\left(\sum_{i=1}^{n} \alpha_{i}, \lambda\right) \\

7.设X\sim Ga(\alpha,\lambda),则

Y=kX\sim Ga(\alpha,\frac{\lambda}{k}),k\neq0 \\

8.伽马分布的两个特例:指数分布和卡方分布。

​(1)\alpha=1时的伽马分布就是指数分布,即Ga(1,\lambda)=\exp(\lambda),指数分布\exp(\lambda)的密度函数为

p(x)=\lambda e^{-\lambda x}=\frac{\lambda^{1}x^{1-1}e^{-\lambda x}}{\Gamma(1)},\quad x\geq0 \\

​(2)称\alpha=n/2,\lambda=1/2时的伽马分布是自由度为n的\chi^2分布,记为\chi^2(n),即Ga(\frac{n}{2},\frac{1}{2})=\chi^2(n),其密度函数为:

p(x)=\frac{(1/2)^{n/2}x^{n/2-1}e^{-x/2}}{\Gamma(n/2)},\quad x\geq0 \\

9.与其他分布的关系

指数分布,伽马分布,泊松分布之间的关系。

从意义上看:

指数分布解决的问题是“要等到下一次随机事件发生,需要经历多久时间”;

泊松分布解决的是 “在特定时间里发生k个事件的概率”。

(1)泊松(Poisson)分布是预测在一个固定时间间隔内,随机事件发生n次的概率。而指数分布是预测下 一次时间发生需要等待多长的时间。比如,下面几种情况: 直到客户下次购买商品为止 (成功) 的时间、机器下次发生故障的时间、下次公交车到达需要等待的时间。

推导过程:

在泊松分布中有一个参数 \lambda ,其表示在单位时间区间内事件发生次数的平均值,其是一个单位时间 的比例(rate)值。那么 \lambda 倒数 1 / \lambda 是什么含义呢? 1 / \lambda 就表示 事件发生一次需要的时间的平均值。比如,当 \lambda= 时,表示在单位时间内平均发生了次,倒过来 (1 / \lambda=1 / ) 就是发生一次需要4个单位时间。 指数概率分布表示,在一个泊松过程中,两次事件发生的间隔时间的概率分布。用泊松分布的表达 就是在等待的之间范围内一次事件都没发生,这意味着有 \operatorname{Poisson}(x=0) 。

\begin{array}{c} p(x=k)=\frac{\lambda^{k} e^{-\lambda}}{k !} \\ p(x=0)=\frac{\lambda^{0} e^{-\lambda}}{0 !}=e^{-\lambda} \end{array} \\

需要注意的是,在泊松分布的概率函数中时间间隔仅仅是1个单位时间(unit time)。如果想要建立一 个在任意时间区间 (而不是一个单位时间)无事件发生的概率分布,我们需要做些什么呢? 泊松分布假设事件之间是相互独立的,计算t个单位时间内 0 次发生的概率可以是每个单位时间内0次发生概率的连乘。

\begin{aligned} \underbrace{p(T>t)}_{\text {在 } \mathrm{t} \text { 个单位时间内没有发生 }} &=\underbrace{p(x=0)}_{\text {第 } 1 \text { 个单位时间 }} \times \underbrace{p(x=0)}_{\text {第 } 2 \text { 个单位时间 }} \cdots \underbrace{p(x=0)}_{\text {第 } \mathrm{t} \text { 个单位时间 }} \\ &=e^{-\lambda t} \end{aligned} \\

T 表示下次事件发生的时间, p(T>t) 就表示下次事件发生的时间晚xxx的概率,换句话说,就是 t 时间内没有发生的概率。那么 t 时间xxx生的概率就是:

p(T \leq t)=1-p(T>t)=1-e^{-\lambda t} \\

上式是累计分布函数,通过对CDF进行微分(求导)就得到了分布的概率密度函数。

f(t)=\frac{d p(T \leq t)}{d t}=\lambda e^{-\lambda t} \\

如果单位时间内事件发生次数符合泊松分布,那么事件发生的间隔时间就服从指数分布。

(2)指数分布是预测下一次时间发生等待的时间,更进一步,如果需要预测第 k 次事件发生需要等待的时间呢? 这就是 Gamma 分布。

我们用符号 k 表示事件发生的次数,符号 T 表示直到 k 次事件发生等待的时间,也就是目标随机变量。 \lambda 表示泊松过程中事件发生的频率(单位时间xxx生的次数)。 p(T>t ; k) 表示直到第 k 次事件发生等待的时间 T 大于 t 的概率。用符号 p(k ; T=t) 表示在泊松过程中 t 个时间单元内事件发生 k 次的概率。

泊松分布表示单位时间窗口内事件发生次数的概率分布,泊松分布的概率分布函数为:

p(k)=\frac{\lambda^{k}}{k !} e^{-\lambda} \\

其表示在一个单位时间内事件发生 k 次的概率,唯一的参数 \lambda 表示单位时间内事件发生次数的平均 值,也就是 k 的期望值, \mathbb{E}[x]=\lambda 。现在要想计算 t 个单位时间内事件发生的次数,只需要用 \lambda t 替换上式中的 \lambda 即可, t 个时间单元内事件发生 k 次的概率为:

\begin{array}{c} p(k ; t)=\frac{(\lambda t)^{k}}{k !} e^{-\lambda t} \end{array} \\

泊松分布表示的是 一个单位时间片段内随机事件发生 次数 的概率分布, 指数分布表示的是随机事件 第一次 发生的需要的 时间 的概率分布。 泊松分布描述的是 次数,是离散变量的分布; 指数分布描述的是 时间 ,是连续值变量的分布。 然而指数分布仅仅描述了事件首次发生,不能表示更多次事件发生,不具备一般化, Gamma 分布就是指数分布的扩展,能够表达事件发生任意次数所需 时间 的概率分布。

我们令符号 p(T>t ; k) 表示随机事件发生 k 所需时间大于 t 的概率。如果随机事件第 k 发生的时 间大于 t ,意味着在 t 个时间单元内,事件发生次数一定是小于等于 k-1 次。换句话说, t 个时间单元内最多只发生 k-1 次,第 k 次发生的时间一定是大于 t 。因此 p(T>t ; k) 可以看成是 在 t 个时间单元内,事件发生 0,1,2, \ldots, k-1 次的概率之和。

p(T>t ; k)=\sum_{i=0}^{k-1} p(i ; t)=\sum_{i=0}^{k-1} \frac{(\lambda t)^{i}}{i !} e^{-\lambda t} \\

现在我们得到了第 k 次事件发生时间大于 t 的概率 p(T>t ; k) ,反过来,第 k 次事件发生时间小 于等于 t 的概率为:

\begin{aligned} p(T \leq t ; k) &=1-p(T>t ; k) \\ &=1-\sum_{i=0}^{k-1} \frac{(\lambda t)^{i}}{i !} e^{-\lambda t} \end{aligned} \\

显然是一个累积分布函数,是概率密度函数的积分,对其进行微分可以得到概率密度函数,这里省略微分的过程,直接给出结果

\begin{aligned} f(t) &=\frac{\lambda e^{-\lambda t}(\lambda t)^{k-1}}{(k-1) !} \\ &=\frac{\lambda^{k} e^{-\lambda t} t^{k-1}}{\Gamma(k)} \end{aligned} \\

上式就是 Gamma 分布的概率密度函数,其表示随机事件发生 k 次所需时间 t 的概率分布。 注意, f(t) 不是具体的时间值,而是时间 t 的概率分布。参数 k 表示事件发生的次数,又被称为 形状参数(shape parameter)。参数 \lambda xxx泊松分布,表示随机事件发生的频次,即单位时间发生的平均次数,又被称为速率参数(rate parameter)。

10.伽马分布是指数型分布族,拥有指数型分布族的全部性质。

我们将伽马分布的密度形式表示为:

\begin{align} p_{\alpha,\lambda}(x)=\frac{\lambda^{\alpha} }{\Gamma(\alpha)}x^{\alpha-1} e^{-\lambda x} &=\frac{\lambda^{\alpha} }{\Gamma(\alpha)}\exp \{(\alpha-1)\ln x-\lambda x\}\\&=c(\alpha,\lambda)\exp\{c_1(\alpha,\lambda)T_1(x)+c_2(\alpha,\lambda)T_2(x)\}h(x) \end{align} \\

其支撑集为\{x>0\},与参数\alpha,\lambda无关,且

c(\alpha,\lambda)=\frac{\lambda^{\alpha} }{\Gamma(\alpha)},h(x)=1 \\c_1(\alpha,\lambda)=\alpha-1,T_1(x)=\ln x \\c_2(\alpha,\lambda)=-\lambda,T_2(x)=x \\

并且参数空间\Theta=\{(\alpha,\lambda):\alpha>0,\lambda>0\}含有内点,根据指数型分布族的结论,伽马分布族参数(\alpha,\lambda)的充分完备统计量为(\sum^n_{i=1}\ln x_i,\sum^n_{i=1}x_i).

1.设X_1,X_2,\cdots,X_n为来自\Gamma(\alpha,\lambda)的独立随机样本,则\bar{X}=\frac{\sum^n_{i=1}X_i}{n}是\frac{\alpha}{\lambda}的UMVUE,并且当\alpha已知时,求\lambda的UMVUE.

证明:注意到T=(\sum^n_{i=1}\ln X_i,\sum^n_{i=1}X_i)是参数(\alpha,\lambda)的充分完备统计量,记g(x,y)=\frac{y}{n},则E(g(T))=E(\bar{X})=\frac{\alpha}{\lambda},由Lehmann-Scheffe定理(充分完备统计量的函数如果是无偏估计,那么也是UMVUE),\bar{X}是\frac{\alpha}{\lambda}的UMVUE.

当\alpha已知时,由伽马函数的可加性,S=\sum^n_{i=1}X_i\sim\Gamma(n\alpha,\lambda),因此E(\frac{1}{S})=\int^{\infty}_{0}\frac{1}{s}\frac{\lambda^{n\alpha}}{\Gamma(n\alpha)}s^{n\alpha-1}e^{-\lambda s}\mathrm{~d}s=\frac{\lambda}{n\alpha-1},则知E(\frac{n\alpha-1}{S})=\lambda,令h(x,y)=\frac{n\alpha-1}{y},则E(h(T))=E(\frac{n\alpha-1}{S})=\lambda,由Lehmann-Scheffe定理,\frac{n\alpha-1}{S}是\lambda的UMVUE。

2.(《数据科学统计基础》P54 例)

电子产品的失效常常是由于外界的 “冲击引起”. 若在 (0, t) xxx生冲击的次数 N(t) 服从参数为 \lambda t 的泊松分布, 试证第 n 次冲击来到的时间 S_{n} 服从伽马分布 \operatorname{Ga}(n, \lambda) . 证明:因为事件“第 n 次冲击来到的时间 S_{n} 小于等于 t ” 等价于事件 “ (0, t) xxx 生冲击的次数 N(t) 大于等于 n _, 即

\left\{S_{n} \leqslant t\right\}=\{N(t) \geqslant n\} . \\

于是, S_{n} 的分布函数为

F(t)=P\left(S_{n} \leqslant t\right)=P(N(t) \geqslant n)=\sum_{k=n}^{\infty} \frac{(\lambda t)^{k}}{k !} \mathrm{e}^{-\lambda t} . \\

由分部积分

\begin{aligned} \frac{\lambda^{n}}{\Gamma(n)} \int_{t}^{\infty} x^{n-1} \mathrm{e}^{-\lambda x} \mathrm{~d} x &=\frac{\lambda^{n}}{(n-1)!} \int_{t}^{\infty}(-\frac{1}{\lambda}) x^{n-1} d\left(e^{-\lambda x}\right) \\ &=\left.\frac{\lambda^{n}}{(n-1) !}(-\frac{1}{\lambda})x^{n-1} e^{-\lambda x}\right|_{t} ^{\infty}+\frac{ \lambda^{n-1}}{(n-1) !} \int_{t}^{\infty}(n-1) x^{n-2} e^{-\lambda x} d x \\ &=\frac{(\lambda t)^{n-1}}{(n-1) !} \mathrm{e}^{-\lambda t} +\frac{ \lambda^{n-1}}{(n-2) !} \int_{t}^{\infty} x^{n-2} e^{-\lambda x} d x\\ &=\cdots \\ &=\sum_{i=1}^{n-1} \frac{(\lambda t)^{k}}{k !} e^{-\lambda t}+\frac{\lambda}{1 !} \int_{t}^{\infty} e^{-\lambda x} dx\\ &=\sum_{k=0}^{n-1} \frac{(\lambda t)^{k}}{k !} \mathrm{e}^{-\lambda t} \end{aligned} \\

\sum_{k=0}^{n-1} \frac{(\lambda t)^{k}}{k !} \mathrm{e}^{-\lambda t}=\frac{\lambda^{n}}{\Gamma(n)} \int_{t}^{\infty} x^{n-1} \mathrm{e}^{-\lambda x} \mathrm{~d} x . \\

F(t)=1-\sum_{k=0}^{n-1} \frac{(\lambda t)^{k}}{k !} \mathrm{e}^{-\lambda t}=\frac{\lambda^{n}}{\Gamma(n)} \int_{0}^{t} x^{n-1} \mathrm{e}^{-\lambda x} \mathrm{~d} x, \\

这就表明 S_{n} \sim G a(n, \lambda) . 证毕.

—END—

Recap:

gamma分布的参数及作用

dgamma/qgamma/pgamma/rgamma函数的用法

gamma分布的应用

若随机变量X的密度函数为

p(x)=\frac {\Gamma (a+b)}{\Gamma (a)\Gamma (b)}x^{a-1}(1-x)^{b-1},0

则称X服从贝塔分布,记作X\sim Be(a,b),其中a>0,b>0都是形状参数,故贝塔分布族可以表示为\{Be(a,b):a>0,b>0\}.

Beta 分布的概率密度我们把它画成图,会发现它是个百变星君,它可以是凹的、凸的、单调上升的、单调下降的; 可以是曲线也可以是直线,而均匀分布也是特殊的 Beta 分布。由于Beta 分布能够拟合如此之多的形状, 因此它在统计数据拟合和xxx分析中被广泛使用。

又因为服从贝塔分布Be(a,b)的随机变量是仅在区间(0,1)取值的,所以不合格率、机器的维修率、市场占有率、射击的命中率等各种比率选用贝塔分布作为他们的概率分布是可能的,只要选择适合的参数a和b即可。

若随机变量X满足Beta分布Be(a,b),则它有以下性质

1.期望

\mu=EX=\frac{a}{a+b} \\

证明:

E(X)=\frac{\Gamma(a+b)}{\Gamma(a) \Gamma(b)} \int_{0}^{1} x^{a}(1-x)^{b-1} \mathrm{~d} x=\frac{\Gamma(a+b)}{\Gamma(a) \Gamma(b)} \cdot \frac{\Gamma(a+1) \Gamma(b)}{\Gamma(a+b+1)}=\frac{a}{a+b} \\

2.方差

\operatorname{Var}(X)=\frac{a b}{(a+b)^{2}(a+b+1)} \\

证明:先计算E(X^2)

\begin{aligned} E\left(X^{2}\right) &=\frac{\Gamma(a+b)}{\Gamma(a) \Gamma(b)} \int_{0}^{1} x^{a+1}(1-x)^{b-1} \mathrm{~d} x\\ &=\frac{\Gamma(a+b)}{\Gamma(a) \Gamma(b)} \cdot \frac{\Gamma(a+2) \Gamma(b)}{\Gamma(a+b+2)} \\ &=\frac{a(a+1)}{(a+b)(a+b+1)} . \end{aligned}d \\

则可得随机变量X​的方差为

\operatorname{Var}(X)=\frac{a(a+1)}{(a+b)(a+b+1)}-\left(\frac{a}{a+b}\right)^{2}=\frac{a b}{(a+b)^{2}(a+b+1)} \\

3.偏度\beta_s,峰度\beta_k

\beta_s=\frac{\mathrm{E}(X-\mu)^{3}}{\left[\mathrm{E}(X-\mu)^{2}\right]^{3 / 2}}=\frac{2(b-a) \sqrt{a+b+1}}{(a+b+2) \sqrt{a b}} \\\beta_k=\frac{\mathrm{E}(X-\mu)^{4}}{\left[\mathrm{E}(X-\mu)^{2}\right]^{2}}-3=\frac{6\left[a^{3}-a^{2}(2 b-1)+b^{2}(b+1)-2 a b(b+2)\right]}{a b(a+b+2)(a+b+3)} \\

\frac{6\left[(a-b)^{2}(a+b+1)-a b(a+b+2)\right]}{a b(a+b+2)(a+b+3)} \\

\mathrm{E}\left(X^{k}\right)=\frac{\mathrm{B}(a+k, b)}{\mathrm{B}(a, b)}=\frac{\Gamma(a+b)\Gamma(a+k)}{\Gamma(a+b+k) \Gamma(a)} \\

k阶矩还可以递归地表示为:

\mathrm{E}\left(X^{k}\right)=\frac{a+k-1}{a+b+k-1} \mathrm{E}\left(X^{k-1}\right) \\

共轭:

假定二项分布 b(n, p) 的参数 p 服从 B e(\alpha, \beta) 先验分布, 然后又做 了 n_{1}+n_{2} 次伯努利实验 (记为 \\right) , 成功 n_{1} 次, 失败 n_{2} 次, 于是后验分布为

\begin{aligned} P(p \mid W) &=\frac{P(p, W)}{P(W)}=\frac{P(W \mid p) P(p)}{\int_{0}^{1} P(W \mid p) P(p) d p} \\ &=\frac{C_{n_{1}+n_{2}}^{n_{1}} p^{n_{1}}(1-p)^{n_{2}} \frac{1}{B(\alpha, \beta)} p^{\alpha-1}(1-p)^{\beta-1}}{\int_{0}^{1} C_{n_{1}+n_{2}}^{n_{1}} p^{n_{1}}(1-p)^{n_{2}} \frac{1}{B(\alpha, \beta)} p^{\alpha-1}(1-p)^{\beta-1} d p} \\ &=\frac{p^{n_{1}+\alpha-1}(1-p)^{n_{2}+\beta-1}}{\int_{0}^{1} p^{n_{1}+\alpha-1}(1-p)^{n_{2}+\beta-1} d p}=\frac{p^{n_{1}+\alpha-1}(1-p)^{n_{2}+\beta-1}}{B\left(n_{1}+\alpha, n_{2}+\beta\right)} \end{aligned} \\

即此时 p 服从 B e\left(n_{1}+\alpha, n_{2}+\beta\right) 分布,我们知道xxx参数估计的基本过程是

先验分布+数据的知识=后验分布 \\

以上xxx分析过程的简单直观的表述就是

\operatorname{Be}(\alpha, \beta)+\operatorname{BinomCount}\left(n_{1}, n_{2}\right)=\operatorname{Be}\left(n_{1}+\alpha, n_{2}+\beta\right) \\

对于上式, 其实也可以纯粹从xxx的角度来进行推导和理解。 假设有一个不均匀的硬币抛出正面的概率为 p , 抛 m 次后出现正面和反面的次数分别是 m_{1}, m_{2} , 那么按传统的频率学派观点, p 的估计值应该为 \hat{p}=\frac{m_{1}}{m} . 而从xxx学派的观点来看, 开始对硬币不均匀性一无所知, 所以应该假设 p \sim U(0,1) , 于是有了二项分布的计数 \left(m_{1}, m_{2}\right) 之后, 按照xxx公式如下计算 p 的后验分布

\begin{aligned} P\left(p \mid m_{1}, m_{2}\right) &=\frac{P(p) \cdot P\left(m_{1}, m_{2} \mid p\right)}{P\left(m_{1}, m_{2}\right)}=\frac{1 \cdot P\left(m_{1}, m_{2} \mid p\right)}{\int_{0}^{1} P\left(m_{1}, m_{2} \mid t\right) d t}=\frac{C_{m}^{m_{1}} p^{m_{1}}(1-p)^{m_{2}}}{\int_{0}^{1} C_{m}^{m_{1}} t^{m_{1}}(1-t)^{m_{2}} d t} \\ &=\frac{p^{m_{1}}(1-p)^{m_{2}}}{\int_{0}^{1} t^{m_{1}}(1-t)^{m_{2}} d t}=\frac{\Gamma\left(m_{1}+m_{2}+2\right)}{\Gamma\left(m_{1}+1\right) \Gamma\left(m_{2}+1\right)} p^{m_{1}}(1-p)^{m_{2}} \end{aligned} \\

计算得到的后验分布正好是 p \sim B e\left(m_{1}+1, m_{2}+1\right) .(由于均匀分布是贝塔分布的一个特例,所以上述推导的成立是自然的)

分布高维形式的推广——Dirichlet分布

我们给出Dirichlet Distribution(xxx雷分布) 的定义。假设 X=(x_{1}, \ldots, x_{k}),\left(0 服从Dirichlet分布 \operatorname{Dir}\left(\alpha_{1}, \ldots, \alpha_{m}\right) ,则密度函数

p\left(x_{1}, \ldots, x_{m}\right)=\frac{\Gamma\left(\sum_{k=1}^m \alpha_{k}\right)}{\prod_{k=1}^m \Gamma\left(\alpha_{k}\right)} \prod_{k=1}^{m} x_{k}^{\alpha_{k}-1} \\

可以直观看出,Beta分布是Dirichlet分布的二元形式。类似于Beta-Binomial 共轭,Dirichlet分布是多项式分布(二项分布的高维推广)的共轭先验分布,Dirichlet-Multinomial共轭表示如下:

\operatorname{Dir}(\vec{\alpha})+M u l t \operatorname{Count}(\vec{m})=\operatorname{Dir}( \vec{\alpha}+\vec{m}) \\

1.顺序统计量建模

有一天你被魔鬼撒旦抓走了,撒旦说:“你们人类很聪明,而我是很仁慈的,和你玩一个游戏,赢了就可以走,否则把灵魂出卖给我。游戏的规则很简单,我有一个魔盒,上面有一个按钮,你每按一下按钮,就均匀的输出一个[0,1]之间的随机数,我现在按10下,我手上有10个数,你猜第7大的数是什么,偏离不超过就算对。”你应该怎么猜呢?

从数学的角度抽象一下,上面这个游戏描述如下:

1: X_{1}, X_{2}, \cdots, X_{n} \stackrel{\mathrm{iid}}{\sim} U n i \operatorname{form}(0,1)

2: 把这 n 个随机变量排序后得到顺序统计量 X_{(1)}, X_{(2)} \cdots, X_{(n)}

3: 问 X_{(k)} 的分布是什么?

对于上面的游戏而言 n=10, k=7 , 如果我们能求出 X_{(7)} 的分布的概率密度, 那么用概率密度的极值点去做猜测就是最好的策略。对于一般的情形, X_{(k)} 的分布是什么呢? 那我们尝试计算一下 X_{(k)} 落在一个区间 [x, x+\Delta x] 的概率, 也就是求如下概率值

P\left(x \leq X_{(k)} \leq x+\Delta x\right)=? \\

把 [0,1] 区间分成三段 [0, x),[x, x+\Delta x],(x+\Delta x, 1] , 我们先考虑简单的情形, 假设 n 个数中只有一个落在了区间 [x, x+\Delta x] 内, 则因为这个区间内的数 X_{(k)} 是第 k 大的, 则 [0, x) 中应该有 k-1 个数, (x, 1] 这个区间中应该有 n-k 个数。不失一般性, 我们先考虑如下一个符合上述要求的事件 E

\begin{aligned} E=\{& X_{1} \in[x, x+\Delta x], \\ X_{i} & \in[0, x) \quad(i=2, \cdots, k), \\ X_{j} &\in(x+\Delta x, 1] \quad(j=k+1, \cdots, n)\} \end{aligned} \\

\begin{aligned} P(E) &=\prod_{i=1}^{n} P\left(X_{i}\right) \\ &=x^{k-1}(1-x-\Delta x)^{n-k} \Delta x \\ &=x^{k-1}(1-x)^{n-k} \Delta x+o(\Delta x) \end{aligned} \\

o(\Delta x) 表示 \Delta x 的高阶无穷小。显然, 由于不同的排列组合, 即 n 个数中有一 个落在 [x, x+\Delta x] 区间的有 n 种取法, 余下 n-1 个数中有 k-1 个落在 [0, x) 的 有 \left(\begin{array}{c}n-1 \ k-1\end{array}\right) 种组合, 所以和事件 E 等价的事件一共有 n\left(\begin{array}{c}n-1 \ k-1\end{array}\right) 个。 继续考虑稍微复杂一点情形, 假设 n 个数中有两个数落在了区间 [x, x+ \Delta x] ,

\begin{aligned} E^{\prime}=\{& X_{1}, X_{2} \in[x, x+\Delta x] \\ & X_{i} \in[0, x) \quad(i=3, \cdots, k) \\ &\{j} \in(x+\Delta x, 1] \quad(j=k+1, \cdots, n)\right\} \end{aligned} \\

P\left(E^{\prime}\right)=x^{k-2}(1-x-\Delta x)^{n-k}(\Delta x)^{2}=o(\Delta x) \\

从以上分析我们很容易看出, 只要落在 [x, x+\Delta x] 内的数字超过一个, 则对应的事件的概率就是 o(\Delta x) 。于是

\begin{array}{l} P\left(x \leq X_{(k)} \leq x+\Delta x\right) \\ =n\left(\begin{array}{l} n-1 \\ k-1 \end{array}\right) P(E)+o(\Delta x) \\ =n\left(\begin{array}{l} n-1 \\ k-1 \end{array}\right) x^{k-1}(1-x)^{n-k} \Delta x+o(\Delta x) \end{array} \\

所以, 可以得到 X_{(k)} 的概率密度函数为

\begin{aligned} f(x) &=\lim _{\Delta x \rightarrow 0} \frac{P\left(x \leq X_{(k)} \leq x+\Delta x\right)}{\Delta x} \\ &=n\left(\begin{array}{l} n-1 \\ k-1 \end{array}\right) x^{k-1}(1-x)^{n-k} \\ &=\frac{n !}{(k-1) !(n-k) !} x^{k-1}(1-x)^{n-k} \quad x \in[0,1] \end{aligned} \\

这其实就是Beta分布Be(k,n-k+1).

好, 我们回到魔鬼的游戏, 在 n=10, k=7 这个具体的实例中, 我们按照如下密度分布的峰值去猜测才是最有把握的。

f(x)=\frac{10 !}{(6) !(3) !} x^{6}(1-x)^{3} \quad x \in[0,1] \\