概率与统计-chapter0-三个重要函数
- Author: zhouyongsdzh@foxmail.com
- Date: 2015-07-04
- Sina Weibo: @周永_52ML
多产的数学家－欧拉
他是公认的数学史上4位最伟大的数学家之一，他的一生没有戏剧性的故事，但给后人留下宝贵的科学财富。
18世纪，数学家辈出的年代，他仍属于佼佼者，被认为是18世纪数学界最杰出的人物之一。
同时，他在数学领域最多产的一位，不仅在数学界做出伟大贡献，而且把数学应用到了几乎整个物理领域。他就是18世纪数学界之星-欧拉。
莱昂哈德·欧拉)（Leonhard Euler ，1707年4月15日～1783年9月18日），瑞士数学家、自然科学家。1707年出生于瑞士－巴塞尔的牧师家庭，15岁从巴塞尔大学毕业，翌年获得硕士学位。
简述其成就：
- _微分方程_: 18世纪中叶，欧拉与其它数学家在解决物理问题过程中，创立了微分方程学科。偏微分方程的纯数学研究的第一篇论文是欧拉写的《方程的积分法研究》。此外，还提出了函数用三角级数表示的方法和解微分方程的级数法等。
- _微分几何学_：引入了空间曲线的参数方程，给出了空间曲线曲率半径的解析表达式；1766年出版了《关于曲面上曲线的研究》，建立了曲面理论，是微分几何发展史上的里程碑。
- _分析学_：1729年引入了Gamma函数和Beta函数，证明了椭圆积分的加法定理，最早引入了二重积分。
- _数论_：欧拉的一系列成果奠定了该数学分支。
注：数学史上公认的4名最伟大的数学家分别是：阿基米德、牛顿、欧拉和高斯。
阿基米德有“翘起地球”的豪言壮语，牛顿因为苹果闻名世界，高斯少年时就显露出计算天赋，唯独欧拉没有戏剧性的故事让人印象深刻。
这里仅总结欧拉在分析学领域中的一个成就－欧拉积分及其对应的分布。欧拉积分是一种含参变量的积分，主要包括
- 第一类欧拉积分：又称为Beta函数（简称B函数）
- 第二类欧拉积分：又称为Gamma函数
Gamma函数－第二类欧拉积分
因为Beta函数和Dirichlet函数可以表示用Gamma函数表示，所以先介绍Gamma函数。
历史由来
1728年，哥德巴赫（Goldbach C. 德国数学家，1690-1764）在考虑数列插值的问题，通俗的说就是把数列的通项公式定义从整数集合延伸至实数集合，例如数列1,4,9,16…..可以用通项公式\(n^2\)自然地表达，即便\(n\)为实数时，此通项公式也具有良好的定义。另一个自然的解释就是可以找到一条平滑的曲线\(f(x)=x^2\)通过所有的整数点\([x, x^2]\)，从而可以把定义在整数集上的公式延伸至实数集合。
有一天，哥德巴赫开始处理阶乘序列\(1,2,6,24,120,720,\cdots\)，我们可以计算整数的阶乘（如\(2!,3!\)），那么是否可以计算实数的阶乘（如\(2.5!\)）呢？
遗憾的是，哥德巴赫当时无法解决阶乘的计算从整数集合延伸至实数集合上的问题，于是写信请教尼古拉斯.伯努利和其弟丹尼尔.伯努利，由于当时欧拉与丹尼尔.伯努利在一起，因此也得知了该问题。而欧拉与1729年完美的解决了这个问题，由此导致Gamma函数的诞生，当时的欧拉只有22岁。
Gamma发现之旅
其实最早发现\(n!\)差值计算的是丹尼尔.伯努利。他发现：如果\(m,n\)都是正整数，且\(m\to\infty\)时，有下面公式成立：
$$
\frac {1 \cdot 2 \cdot 3 \cdots m} {(1+n)(2+n) \cdots (m-1+n)} \left(m+ \frac{n}{2}\right)^{n-1} \to n!
$$
用此无穷乘积的方式可以把\(n!\)定义延伸至实数集合。如\(n=2.5,m\)足够大时，上式基本可近似计算出\(2.5!\)的值。
欧拉偶然地发现\(n!\)可用如下无穷乘积表达：
$$
\left[\left(\frac{2}{1}\right)^n \frac{1}{n+1}\right] \left[\left(\frac{3}{2}\right)^n \frac{2}{n+2}\right]
\left[\left(\frac{4}{3}\right)^n \frac{3}{n+3}\right]
\cdots = n!
$$
于是他用一些简单的例子做一些计算，寻找其规律。当\(n=1/2\)时，整理上式可得：
$$
\left(\frac{1}{2}\right)!=\sqrt{\frac{2\cdot4}{3\cdot3} \cdot \frac{4\cdot6}{5\cdot5} \cdot \frac{6\cdot8}{7\cdot7} \cdot \frac{8\cdot10}{9\cdot9} \cdots}
$$
此式恰好与著名的Wallis公式有关。
关于Wallis公式：
\( \qquad \lim\limits_{k\to\infty} \left(\frac{2^{2k}(k!)^2}{2k!}\right)^2 \frac{1}{2k+1} = \frac{\pi}{2} \)
Wallis公式是关于圆周率的无穷乘积的公式，Wallis本人在1665年使用插值方法计算半圆曲线\(y=\sqrt{x(1-x)}\)下的面积（即直径为1的半圆面积）时，发现\(\pi/4\)等于下式：
\(
\qquad\frac{2\cdot4}{3\cdot3} \cdot \frac{4\cdot6}{5\cdot5} \cdot \frac{6\cdot8}{7\cdot7} \cdot \frac{8\cdot10}{9\cdot9} \cdots = \frac{\pi}{4}
\)
于是欧拉根据Wallis公式得到如下的结果：
$$
\left(\frac{1}{2}\right)! = \frac{\sqrt{\pi}}{2}
$$
欧拉发现\((\frac{1}{2})!\)中含有\(\pi\)，而\(\pi\)与圆相关的积分有关系。欧拉猜测\(n!\)一定可以用积分形式表示。
Wallis时代微积分还没问世，当时是使用插值方法做推导计算的，但Wallis公式的推导过程基本上都是在处理积分\(\int_{0}^{1} x^\frac{1}{2} (1-x)^\frac{1}{2} dx\)的问题。
关于微积分诞生：
牛顿与莱布尼茨谁先发明的微积分在数学界有很大的争论。不过本人更认可是莱布尼茨，因为微积分的主要思想是他提出的，微分相关概念和符号是他定义的，并沿用至今。
1684年莱布尼茨发表了关于微积分的重要文献，标志着微积分作为独立学科正式诞生。莱布尼茨是微积分主要思想的独立发明人。当然，不可否认的是牛顿在此前同样做了很多微积分相关工作。
受Wallis启发，欧拉开始考虑如下一般形式的积分：
$$
J(e,n)=\int_{0}^{1} x^e(1-x)^n dx
$$
这里，\(n\)为正整数，\(e\)为正实数。利用分部积分法，可得：
$$
J(e,n)=\frac{n}{e+1}J(e+1, n-1)
$$
关于分部积分：
设\(u=u(x)\)和\(v=v(x)\)均为可导函数，分部积分形式如下：\( \qquad \int u \mathop{v’}dx = \int u dv = [uv] - \int v du \)
重复迭代上式，可得：
$$
J(e,n) = \frac {1\cdot2\cdot3\cdots n}{(e+1)(e+2)\cdots (e+n+1)}
$$
于是欧拉得到如下重要的式子：
$$
n!=(e+1)(e+2)\cdots (e+n+1) \int_{0}^{1} x^e(1-x)^n dx
$$
接下来，欧拉使用一些计算技巧，即取\(e=f/g\) 且 \(f\to1, g\to0\)，然后对上式右边计算极限，得到如下简洁的结果：
$$
n! = \int_{0}^{1}(-\log{t})^n dt
$$
到此，欧拉成功地把\(n!\)表达为了积分形式！令\(t=e^{-u}\)，可得常见的Gamma函数形式：
$$
n! = \int_{0}^{\infty} u^n e^{-u} du
$$
注意，此时的\(n\)仍然为正整数，利用上式把阶乘延伸至实数集上，就得到Gamma函数的一般形式：
$$
\Gamma(x) = \int_{0}^{1} (-\log t)^{x-1}dt = \int_{0}^{\infty} t^{x-1} e^{-t} dt
$$
注：使用MathJax引擎在浏览器上解析时，需要对”_”做转义，应写成下面形式：
为什么\(\Gamma(n)=(n-1)!\)，而不是\(\Gamma(n)=n!\)
如果对上述的Gamma函数式做一个调整，将\(t^{x-1}\)替换为\(t^x\)，可得如下函数：
$$
\Gamma(x) = \int_{0}^{\infty} t^x e^{-t}dx
$$
此时\(\Gamma(n)=n!\)，这才是欧拉最早的Gamma函数定义。但后来欧拉修改了原有Gamma函数的定义，使得\(\Gamma(n)=(n-1)!\)。后来的数学家们在对Gamma函数的进一步研究中，认可了这个定义，沿用至今。
Gamma分布
- 概率密度函数
$$
f(x)=
\begin{cases}
\displaystyle\frac {x^{\alpha-1} e^{x/\beta}} {\beta^\alpha \Gamma(\alpha)}, & 0 \leq x \leq \infty; \alpha > 0; \beta > 0 \\\
\;\;\quad 0, & others
\end{cases}
$$
Beta函数－第一类欧拉积分
在@rickjin的《LDA-math》系列中详细地解释了Beta和Dirichlet分布的由来和推导过程，堪称经典。原文见下面链接：
@rickjin在文中从一个魔鬼的游戏开始引入。魔鬼撒旦抓走一人，撒旦说：“你们人类很聪明，而我是很仁慈的，和你玩一个游戏，赢了就可以走，否则把灵魂出卖给我。”
游戏规则：
我有一个魔盒，上面有一个按钮，你每按一下按钮，就均匀的输出一个[0,1]之间的随机数，我现在按10下，我手上有10个数，你猜第7大的数是什么，偏离不超过0.01就算对。
数学抽象与推导
上述游戏实际在说随机变量\(X_1, X_2, \cdots, X_{10} \sim U(0,1)\), 把这\(n\)个随机变量排序后得到顺序统计量\(X_{(1)}, X_{(2)}, \cdots, X_{(10)}\)，然后问\(X_{(k)}\)的分布是什么？
因为如果知道随机变量\(X_{(k)}\)分布的概率密度，用概率密度的极值点作为猜测值是最好的策略。
对于上述游戏而言，\(n=10,k=7\)，\(X_{(k)}\)的分布如何计算？@rickjin在文中列举了通过尝试计算\(X_{(k)}\)落在区间\([x,x+\Delta x]\)的概率，也就是求如下概率值
$$
P(x \leq X_{(k)} \leq x + \Delta x) = \;?
$$
如果把[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+\Delta x,1]\)中应该有\(n-k\)个数。那么，符合上述要求的事件\(E\)可表示为：
$$
\begin{align*}
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{align*}
$$
注：
- 符号
*在markdown中有特殊含义，若要当作latex环境中的语法用，需要转义，即\ *.- 在Mathjax下，
$$E = \{x, y\}$$无法正常输出{}, 因为Latex要想输出{}，需要使用\{, \}。而\在markdown中要想正常表法其自身意义也需要再次转义，代码为:$$E = \\{x,y\\}$$。$$E = \\{ x,y\\}$$
那么事件\(E\)发生的概率，有：
$$
\begin{align*}
P(E) & = \prod_{i=1}^nP(X_i) \\\
& = x^{k-1}(1-x-\Delta x)^{n-k}\Delta x \\\
& = x^{k-1}(1-x)^{n-k}\Delta x + o(\Delta x)
\end{align*}
$$
\(o(\Delta x)\)表示\(\Delta x\)的高阶无穷小。显然，由于不同的排列组合，即\(n\)个数中有一个落在\([x, x+ \Delta x]\)区间的有\(n\)中取法，余下\(n-1\)个数中有\(k-1\)个落在\([0, x)\)区间有\(\binom{n-1}{k-1}\)种组合，所以与事件E等概率的事件一共有\(n \binom{n-1}{k-1}\)个（当只有1个数落在\([x, x+ \Delta x]\)时）。
继续考虑复杂一些的情形，假设\(n\)个数中两个数落在了区间\([x, x+ \Delta x]\)，此时事件\(E^{‘}\)可表示为：
$$
\begin{align*}
E’ = \\{
& X_1,X_2\in [x, x+\Delta x], \\\
& X_i \in [0,x) \quad (i=3,\cdots,k), \\\
& X_j \in (x+\Delta x,1] \quad (j=k+1,\cdots,n)
\\}
\end{align*}
$$
此时，事件\(E^{‘}\)发生的概率为：
$$
P(E^{‘}) = x^{k-2}(1-x-\Delta x)^{n-k}(\Delta x)^2 = o(\Delta x)
$$
很容易看出，只要落在\([x, x+ \Delta x]\)内的数字大于1个，则对应事件的概率就是\(\Delta x\)。于是随机变量\(X_{(k)}\)落在\([x, x+ \Delta x]\)区间的概率：
$$
\begin{align*}
& P( x \le X_{(k)} \le x+\Delta x) \\\
& = n\binom{n-1}{k-1}P(E) + o(\Delta x) \\\
& = n\binom{n-1}{k-1}x^{k-1}(1-x)^{n-k}\Delta x + o(\Delta x)
\end{align*}
$$
进一步可得\(X_{(k)}\)的概率密度函数为
$$
\begin{align*}
f(x) & = \lim_{\Delta x \to 0} \frac {P(x \leq X_{(x)} \leq x + \Delta x)}{\Delta x} \\\
& = n \binom{n-1}{k-1}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{align*}
$$
利用上述介绍的Gamma函数，可以把\(f(x)\)表示为
$$
f(x) = \frac {\Gamma(n+1)}{\Gamma(k) \Gamma(n-k+1)} x^{k-1} (1-x)^{n-k}
$$
\(\Gamma\)函数起源于人们希望把数学计算从整数集合拓展至实数集合。这里另\(\alpha=k, \beta=n-k+1\)，于是可以得到
$$
Beta(x|\alpha,\beta)=f(x)=\frac {\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)} x^{\alpha-1}(1-x)^{\beta-1}
$$
上式即为一般意义上的Beta分布！
\(\beta\)型概率分布
- \(\beta\)型随机变量\(X\)的概率密度函数：
$$
f(x)=
\begin{cases}
\displaystyle\frac {\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)} x^{\alpha-1}(1-x)^{\beta-1}, & 0 \leq x \leq 1;\;\alpha > 0, \beta > 0 \\\
0, & others
\end{cases}
\quad (**)
$$
关于Beta函数详细推导
$$
\begin{align*}
f(x;\alpha,\beta) & = \frac{x^{\alpha-1}(1-x)^{\beta-1}}{\int_{0}^{1}t^{\alpha-1}(1-t)^{\beta-1}dt} \\\
& = \frac{x^{\alpha-1}(1-x)^{\beta-1}}{Beta(\alpha,\beta)}
\end{align*}
$$其中，
$$
Beta(\alpha,\beta)= \int_{0}^{1}t^{\alpha-1}(1-t)^{\beta-1}dt= \frac{\Gamma(\alpha)\Gamma(\beta)}{\Gamma(\alpha+\beta)}
$$$$
\Gamma(\alpha)=\int_{0}^{\infty}t^{\alpha-1}e^{-t}dt
$$
且当\(\alpha\)是正整数时，\(\Gamma(\alpha)=(\alpha-1)!\)
- \(\beta\)型随机变量的期望与方差分别是：
$$\mu=\frac{\alpha}{\alpha+\beta} \;\qquad\qquad\qquad (期望)$$
$$\sigma^2=\frac{\alpha\beta}{(\alpha+\beta)^2(\alpha+\beta+1)} \quad (方差)$$
- \(\beta\)型分布的期望公式详细推导
如果\(p \sim Beta(t|\alpha,\beta)\)，那么\(p\)的期望表示如下：
$$
\begin{align*}
E(p)
& =\int_{0}^{1} t \ast Beta(t|\alpha,\beta)\;dt \\\
& = \int_{0}^{1} t \ast \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}t^{\alpha-1}(1-t)^{\beta-1}\;dt \\\
& =\frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}\int_{0}^{1}t^{\alpha}(1-t)^{\beta-1}\;dt
\end{align*}
$$
上式右边的积分\(\int_{0}^{1}t^{\alpha}(1-t)^{\beta-1}\;dt\)恰好对应到概率分布\(Beta(t|\alpha+1,\beta)\)。对于此分布，根据Beta型随机变量的密度函数在\([0,1]\)之间的累积分布等于1，可得:
$$
\int_{0}^{1} \frac{\Gamma(\alpha+\beta+1)}{\Gamma(\alpha+1)\Gamma(\beta)} t^\alpha (1-t)^{\beta-1}\;dt=1
$$
把上式带入\(E(p)\)计算式，得到
$$
\begin{align*}
E(p) & = \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)} \cdot
\frac{\Gamma(\alpha+1)\Gamma(\beta)}{\Gamma(\alpha+\beta+1)} \\\
& = \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha+\beta+1)}\frac{\Gamma(\alpha+1)}{\Gamma(\alpha)} \notag \\\
& = \frac{\alpha}{\alpha+\beta}
\end{align*}
$$
这说明，对于Beta分布的随机变量，其均值可以用\(\frac{\alpha}{\alpha+\beta}\)来估计。
- \(\beta\)型概率分布曲线
_(此处暂略~)_
回到魔鬼的游戏中，\(n=10, k=7\)时，按照密度分布的峰值去猜测是最好的策略。
$$
f(x)=\frac{10!}{(6)!(3)!} x^6 (1-x)^3 \quad x \in [0,1]
$$
即便按照密度函数分布的峰值作为猜测结果，第一次猜中的概率也不高。
游戏继续\(\cdots\)
很遗憾，根据上述最好的策略算出来的值竟然有偏差，没猜中，魔鬼微笑着说：“我再仁慈一点，再给你一个机会，你按5下这个机器，你就得到了5个\([0,1]\)之间的随机数，然后我可以告诉你这5个数中的每一个，和我的第7个数相比，谁大谁小，然后你继续猜我手头上的第7大的数是多少。” 这时候该如何猜测呢？
Beta-Binomial共轭
魔鬼的两个问题，数学抽象一下，就是：
- \(X_1,X_2, \cdots, X_n \in Uniform(0,1)\)，对应的顺序统计量为\(X_{(1)},X_{(2)}, \cdots, X_{(n)}\)，我们要猜测\(p=X_{(k)}\)；（第1个问题）
- \(Y_1,Y_2, \cdots, Y_m \in Uniform(0,1)\)，\(Y_i\)中有\(m_1\)个比\(p\)小，\(m_2\)个比\(p\)大；（第2个问题）
最后的问题：\(P(p|Y_1,Y_2, \cdots, Y_m)\)的分布是什么？
问题分析
由于\(p=X_(k)\)在\(X_1,X_2, \cdots, X_n\)中是第\(k\)大的，利用\(Y_i\)的信息，可以很容易的推理得到\(p=X_{(k)}\)在\(X_1,X_2, \cdots, X_n,\) \(Y_1, Y_2, \cdots, Y_m \in Uniform(0,1)\)这\((m+n)\)个独立随机变量中是第\(k+m_1\)大的。那么按照Beta分布（贝塔分布）小节的推理，此时\(p=X_{(k)}\)的概率密度函数是\(Beta(p|\,k+m_1, n-k+1+m_2)\)。
根据贝叶斯推理的逻辑，整理上述过程如下：
- \(p=X_{(k)}\)是我们需要猜测的参数，并且推导出\(p\)的分布为\(f(p)=Beta(p|k,n-k+1)\)，称为\(p\)的先验分布；
- 数据\(Y_i\)中有\(m_1\)个比\(p\)小，\(m_2\)个比\(p\)大，\(Y_i\)相当于做了\(m\)次贝努利试验，所以\(m_1\)服从二项分布\(B(m,p)\)；
- 在给定了来自数据提供的\((m_1, m_2)\)的知识后，\(p\)的后验分布变为:
$$f(p|\,m_1,m_2)=Beta(p|\,k+m_1, n-k+1+m_2)$$
关于贝叶斯参数估计的基本过程：
$$先验分布 ＋ 数据知识 ＝ 后验分布$$
在这里，贝叶斯分析过程的简单直观的表述就是：
$$Beta(p|\,k, n-k+1) + Count(m_1, m_2) = Beta(p|\,k+m_1,n-k+1+m_2)$$
其中\(m_1,m_2\)对应的是二项分布\(B(m_1+m_2, p)\)的计数。更为一般地数学表述：对于非负实数\(\alpha,\beta\)，存在如下关系
$$
Beta(p|\,\alpha, \beta) + Count(m_1, m_2) = Beta(p|\,\alpha+m_1, \beta+m_2) \quad (**)
$$
该式描述的就是Beta-Binomial共轭。
关于Beta-Binomial共轭
数据符合二项分布的时候，参数的先验分布和后验分布都能保持Beta分布的形式，这种形式不变的好处是：我们能够在先验分布中赋予参数很明确的物理意义，这个物理意义可以延续至后验分布中进行解释，同时从先验变换到后验的过程中从数据中补充的知识也容易给出物理上的解释。
推导过程中可以看到，\(Beta\)分布中的参数\(\alpha,\beta\)都可以理解为物理计数，这两个参数经常被称为伪计数（pseudo-count）。\(Beta(p|\,\alpha, \beta)\)可以写成如下表达式：
$$Beta(p|\,1,1) + Count(\alpha-1, \beta-1) = Beta(p|\,\alpha,\beta)$$
其中的\(Beta(p|\,1,1)\)恰好就是均匀分布\(Uniform(0,1)\)。
关于Beta分布与均匀分布的关系：
$$Beta(p|\,\alpha=1,\beta=1)= Uniform(a=0,b=1)$$
对于\((**)\)共轭公式，其实可纯粹从贝叶斯的角度来进行推导。理解过程：
- 假设有一个不均匀的硬币抛出正面的概率为\(p\)，抛\(m\)次后出现正面和反面的次数分别是\(m_1,m_2\)，那么按照传统的频率学派观点，\(p\)的估计值应该是\(\hat p = \frac{m_1}{m}\)。
- 而从贝叶斯学派的观点来看，开始对硬币不均匀性一无所知，所以应该假设\(p\ \in Uniform(0,1)\)，于是二项分布的计数\(m_1,m_2\)之后，按照贝叶斯公式计算\(p\)的后验分布：
$$
\begin{align*}
P(p|\,m_1,m_2) = & \frac {P(p) \cdot P(m_1,m_2\,|p)}{P(m_1,m_2)} \\\
= & \frac {1 \cdot P(m_1,m_2 |p)}{\int_{0}^{1} P(m_1,m_2|t)dt} \\\
= & \frac {\binom{m}{m_1}p^{m_1}(1-p)^{m_2}}{\int_{0}^{1} \binom{m}{m_1}t^{m_1}(1-t)^{m_2}dt} \\\
= & \frac {p^{m_1}(1-p)^{m_2}}{\int_{0}^{1}t^{m_1}(1-t)^{m_2}dt}
\end{align*}
$$
得到的后验分布正好是\(Beta(p|\;m_1 + 1, m_2 + 1)\)。
回到魔鬼的游戏，如果按出的5个随机数字中，魔鬼告诉你有2个（即\(m_1=2\)）小于他手中第7大的数，那么应该按照如下概率分布的峰值做猜测是最好的(\(\alpha+2=9,\beta+3=7\))：
$$
Beta(x|9,7)=\frac{15!}{(8)!(6)!}x^8(1-x)^6 \quad x \in [0,1]
$$
很幸运，这次猜中了，但是魔鬼开始耍赖，游戏不得不继续 \(\cdots\)
游戏3新规则：
魔说道：“这个游戏对你来说太简单了，我要加大点难度，我们重新来一次，我按20下生成20个随机数，你要同时给我猜第7大和第13大的数是什么？”
此时又该如何猜测呢？
Dirichlet分布
对于魔鬼变本加厉的新的游戏规则，数学抽象如下：
- \(X_1,X_2, \cdots, X_n \sim Uniform(0,1)\),
- 排序后对应的顺序统计量为\(X_{(1)},X_{(2)}, \cdots, X_{(n)}\),
- 问: \((X_{(k_1)},X_{(k_1 + k_2)})\)的联合分布是什么？
游戏3完全类似游戏1的推导过程，可进行如下的概率计算：
注：为了公式的简洁性，这里取\(x_3=1-x_1-x_2\),但只有\(x_1,x_2\)是变量。说明符号：
- n: 总的样本数
- \(X_{(k_1)}\)和\(X_{(k_1 + k_2)}\)将n分割为3段，长度分别为\(k_1, k_2和k_3\)，即
$$
分段结果：
\begin{cases}
[X_{(1)}, X_{(k-1)}], & 有(k_1-1)个取值，统一记为x_1 \\\
[X_{(k_1+1)}, X_{(k_1+k_2)}], & 有(k_2-1)个取值，统一记为x_2 \\\
[X_{(k_1+k_2+1)}, X_{(n)}] & 有(n-k_1-k_2)个取值，统一记为x_3
\end{cases}
$$
完全类似于游戏1的推导过程，\(X_{(k_1)},X_{(k_1+k_2)}\)的联合概率计算如下：
$$
\begin{align*}
& P\Bigl(X_{(k_1)} \in (x_1, x_1+\Delta x), X_{(k_1+k_2)} \in (x_2, x_2+\Delta x)\Bigr) \\\
& \quad = n(n-1)\binom{n-2}{k_1-1,k_2-1}x_1^{k_1-1}x_2^{k_2-1}x_3^{n-k_1-k_2}(\Delta x)^2 \\\
& \quad = \frac{n!}{(k_1-1)!(k_2-1)!(n-k_1-k_2)!}x_1^{k_1-1}x_2^{k_2-1}x_3^{n-k_1-k_2}(\Delta x)^2
\end{align*}
$$
式中的\(n(n-1)\binom{n-2}{k_1-1,k_2-1}\)是一个排列组合问题，比较容易理解。
于是得到\(X_{(k_1)},X_{(k_1+k_2)}\)的联合分布是：
$$
\begin{align*}
f(x_1,x_2,x_3) & =\frac{n!}{(k_1-)!(k_2-1)!(n-k_1-k_2)!}x_1^{k_1-1}x_2^{k_2-1}x_3^{n-k_1-k_2} \\\
& = \frac{\Gamma(n+1)}{\Gamma(k_1)\Gamma(k_2)\Gamma(n-k_1-k_2+1)}x_1^{k_1-1}x_2^{k_2-1}x_3^{n-k_1-k_2}
\end{align*} \quad (**)
$$
而公式\((**)\)则是3维形式的Dirichlet分布，即\(Dir(x_1,x_2,x_3\;|\;k_1,k_2,n-k_1-k_2+1)\)。这里令\(\alpha_1=k_1,\alpha_2=k_2,\alpha_3=n-k_1-k_2+1\)，于是分布密度函数可以写为：
$$
f(x_1, x_2, x_3) = \frac{\Gamma(\alpha_1+\alpha_2+\alpha_3)}{\Gamma(\alpha_1)\Gamma(\alpha_2)\Gamma(\alpha_3)}x_1^{\alpha_1-1}x_2^{\alpha_2-1}x_3^{\alpha_3-1} \qquad (0)
$$
公式\((0)\)即为一般形式的3维Dirichlet分布。即使\(\vec{\alpha}=(\alpha_1,\alpha_2,\alpha_3)\)延伸至非负实数集合，上述概率分布也是成立的。
在游戏3的基础上还可以往更高的维度上推，譬如\(X_{(1)},X_{(2)},\cdots,X_{(n)}\)中的4、5、…等更多个数，可以得到更高维度的Dirichlet分布和Dirichlet-Multinomial共轭。
- Dirichlet分布
如果\(\vec{p} \sim Dir(\vec{t}|\vec{\alpha})\),那么随机变量\(\vec{p}\)的概率密度函数为
$$
\begin{equation}
\displaystyle Dir(\overrightarrow{p}|\overrightarrow{\alpha}) =
\displaystyle \frac{\Gamma(\sum_{k=1}^{K} \alpha_k)}
{\prod_{k=1}^K \Gamma(\alpha_k)} \prod_{k=1}^{K} p_k^{\alpha_k -1}
\end{equation} \qquad (1)
$$
对于给定的\(\overrightarrow{p}\)和\(N\)，多项式分布定义为：
$$
Multi(\overrightarrow{n}|\overrightarrow{p},N)=
\binom{N}{\overrightarrow{n}} \prod_{k=1}^K p_k^{n_k} \qquad (2)
$$
而\(Multi(\overrightarrow{n}|\overrightarrow{p},N)\)和\(Dir((\overrightarrow{p}|\overrightarrow{\alpha})\)这两个分布是共轭关系。
- Dirichlet分布均值
$$
E(\vec{p}) = \Bigl(\frac{\alpha_1}{\sum_{i=1}^{K}\alpha_i},\frac{\alpha_2}{\sum_{i=1}^{K}\alpha_i}, \cdots, \frac{\alpha_K}{\sum_{i=1}^{K}\alpha_i}\Bigr)
$$
关于Dirichlet分布与Beta分布：
Dirichlet分布是Beta分布在高维度上的推广。当Dirichlet分布维度趋向无限时，便成为Dirichlet过程。
Dirichlet-Multinomial共轭
类似于魔鬼的游戏2，调整一下游戏3，可得游戏4的规则如下：
游戏4:
从魔盒中生成m个随机数\(Y_1,Y_2, \cdots, Y_m \sim Uniform(0,1)\)，魔鬼告诉我们\(Y_i\)和\(X_{(k_1)},X_{(k_1+k_2)}\)相比谁大谁小，然后再次猜测第7大和第13大的数是多少？
同样，对游戏4进行数学抽象，表示如下：
- \(X_1,X_2, \cdots, X_n \sim Uniform(0,1)\)，排序后对应的顺序统计量为\(X_{(1)},X_{(2)},\cdots,X_{(n)}\)；
- 令\(p_1=X_{(k_1)},p_2=X_{(k_1+k_2)},p_3=1-p_1-p_2(p_3是为了简洁的数学表达)\)，我们要猜测\(\vec{p}=(p_1,p_2,p_3)\)；
- \(Y_1,Y_2, \cdots, Y_m \sim Uniform(0,1)\),\(Y_i\)中落到\([0,p_1), [p_1,p_2),[p_2,1]\)三个区间的个数分别是\(m_1,m_2,m_3\),其中\(m=m_1+m_2+m_3\)；
- 问后验分布\(P(\vec{p}|Y_1,Y_2, \cdots, Y_m)\)的分布是什么？
为了计算方便，记：
$$
\vec{m}=(m_1, m_2, m_3), \quad \vec{k}=(k_1,k_2,n-k_1-k_2+1)
$$
从游戏中的信息可以得知：\(p_1,p_2\)在\(X_1,X_2, \cdots, X_n,Y_1,Y_2, \cdots, Y_m \sim Uniform(0,1)\)这\(m+n\)个数中分别成为了第\(k_1+m_1, k_2+m_2\)大的数，于是后验分布\(P(\vec{p}|Y_1,Y_2,\cdots, Y_m)\)应该是\(Dir(\vec{p}|k_1+m_1,k_2+m_2, n-k_1-k_2+1+m_3)\)，即\(Dir(\vec{p}|\vec{k}+\vec{m})\)。按照贝叶斯推理的逻辑，同样可以把上述过程整理如下：
- 首先，根据游戏3需要猜测参数\(\vec{p}=(p_1,p_2,p_3)\)，其先验分布为\(Dir(\vec{p}|\vec{k})\)；
- 其次，数据\(Y_i\)落到\([0,p_1),[p_1,p_2),[p_2,1]\)三个区间的个数分别为\(m_1,m_2,m_3\)，所以\(\vec{m}=(m_1, m_2, m_3)\)服从多项式分布\(Multi(\vec{m}|\vec{p})\)（理解这一点很重要）；
- 在给定了来自数据提供的知识\(\vec{m}\)后，\(\vec{p}\)的后验分布变为\(Dir(\vec{p}|\vec{k}+\vec{m})\)。
以上贝叶斯分析过程的最简单直接的表述：
$$
Dir(\overrightarrow{p}|\overrightarrow{k})+Multi(\vec{m})=Dir(\vec{p}|\vec{k}+\vec{m})
$$
令\(\overrightarrow{\alpha}=\vec{k}\)，把\(\vec{\alpha}\)从整数集合延拓至实数集合，可以证明如下关系仍然成立：
$$
Dir(\vec{p}|\vec{\alpha})+Multi(\vec{m})=Dir(\vec{p}|\vec{\alpha}+\vec{m}) \quad (11)
$$
公式\((11)\)描述就是Dirichlet-Multinomial共轭。同时，我们可以把Dirichlet分布中的\(\alpha\)都可以理解为物理计数。那么，类似于Beta分布，\(Dir(\vec{p}|\vec{\alpha})\)可做如下分解：
$$
Dir(\vec{p}|\vec{1})+Multi(\vec{m}-\vec{1}) = Dir(\vec{p}|\vec{\alpha})
$$
这里\(\vec{1}=(1,1,\cdots,1)\),上式同样可以类似的用纯粹贝叶斯的观点推导和解释。