引言
概率一直是个微妙的话题,概率究竟是什么,不同的学派和价值观给出了不同的解释,以至于数学上干脆只给出了满足条件的函数作为定义,而没有任何对概率本身意义的阐述。关于概率两种最重要的方法论学派就是 frequentist 和 Bayenisan。前者认为概率是确定的值,虽然可能无法直接得到,但可以通过大量重复实验估计,根据大数定理,频率将会趋近于概率。而后者只把概率当作一种确定性的描述,根据获取的信息和知识不同,这一概率也可以不断的调整变动。对于同一陈述,不同的人由于知识经验的不同,给出的概率可能也不同。
本文讨论一个著名的概率问题,并通过该问题对不同的统计和概率方法的哲学进行分析和对比。本文主要参考了一些维基百科的相关词条,主要包括 prior probability,Beta distribution,sunrise problem,rule of succession。
太阳升起的概率
明天太阳升起的概率是多少?对于 frequentist 的世界观来说,这种问题是莫名其妙的。一个由已知物理定律描述的确定性事件何谈概率?不过对于 Bayenisan 来说,这就是一个 well defined 的问题。任何事情包括客观陈述都可以计算概率。就像引言中所说,概率只是反应某种情况下某人对该陈述的确定性的大小。比如你早饭是否吃了米饭,这本来是已经发生的确定的事件,但贝叶斯学派依旧可以根据你今天的一些迹象来推断出这一事件的“概率”。很明显,本质上两种学派下概率的意义完全不同。
数学描述
抽象一下太阳升起问题,对于一个试验(最多有成功失败两种可能),独立重复进行了 \(N\) 次,\(S\) 次成功,\(F=N-S\) 次失败,问该试验每次成功的概率,也即下一次试验成功的概率 \(P\) 是多少。太阳升起问题就是以上问题 \(S=N\) 的情形。
如果是 frequentist 回答这个问题,就特别简洁明了,频率趋向概率,那么能够给出的对概率最好的估计就是\(P\approx \frac{S}{N}.\) 对于太阳升起问题,概率就是 1, 和物理定律的确定性相符。更严格一点,就是所谓的 MLE(Maximum Likelihood Estimation) 方法来求的概率,这一方法里将概率 \(P\) 视为固定的未知参数,那么上述问题的可能性为(这里我们将每次试验记为\(X_i\),对应结果记为\(x_i=1,0\),试验成功次数由二项分布描述。)
\[L_P(\mathbf{X})=P(\mathbf{X}\vert P=\theta)=C_N^S \theta^S(1-\theta)^{(N-S)}.\]于是我们通过最大化这一可能性来得到概率的值(对上式关于 \(\theta\) 求导取零点即可),得到的结果 \(P=\theta=\frac{S}{N}\),与直接按频率估计相同。
这一处理框架有很多不足。首先我们无法对估算的概率的不确定度有任何定量的描述。虽然我们通过大量试验估计出了一个概率值,但无法估算这一概率值的误差。其次,对于一个成功 \(N\) 次的试验武断地认为成功概率为 1,不太符合这门科学的哲学。
贝叶斯推断
我们现在由贝叶斯学派的视角,来分析这一概率问题。这一过程的核心就是贝叶斯概率公式,哲学是将概率 \(P\) 不再视为固定参数,而是视为一个随机变量。那么我们通过分析 \((P\vert \mathbf{X})\) 这一随机变量的概率分布,就得到了当我们知道试验数据 \(\mathbf{X}\) 后,概率参数 \(P\) 的分布。一旦我们得到的是 \(P\) 的整个分布,而不是一个估值,获取的信息量自然大大增加。我们既可以和上文方案一样,取这一分布的概率极值 (离散分布时的众数)处的 \(P\) 值作为概率参数的估计,这就是所谓的 MAP (Maximum a Posterior)方法。也可以取这一分布的均值作为对概率的估计。同时我们还可以计算该分布的方差作为对我们求的概率 \(P\) 的误差的估计。
总之我们要计算的核心就是
\[P(P\vert \mathbf{X})\propto P(P)P(\mathbf{X}\vert P).\]分母上的 \(P(X)\) 并不需要额外计算,对于分子最后做一个归一化即可,或者如果只是比较两种可能性的概率的比值的话,直接用分子相比即可,因此我们这里的贝叶斯公式省略了分母。对于\(P(\mathbf{X}\vert P)\) 我们称为可能性,这一值已在前文用二项式分布给出。对于 \(P(P)\) 称为先验,表达了我们没做这一系列试验时,我们对于概率参数 \(P\) 的分布的信念。最终计算出的结果则是后验。
如果对于问题的模型结构有比较清楚的认识,可能性部分通常较容易求得。贝叶斯推断最玄学的地方,就在于先验的选取。对于概率类参数的分布,我们往往选取的描述参数分布的先验是贝塔分布的形式,也即
\[P(P)=Beta(\alpha,\beta)=\frac{p^{\alpha-1}(1-p)^{\beta-1}}{B(\alpha,\beta)}\]其中 B 是 Beta 函数,我默认读者对于 Beta 函数和 Gamma 函数的联系,相关的积分表达式等简单内容都很熟悉,这里不再赘述。
几种做法
最后我们把求太阳升起的概率 \(P\) 这一问题,归化到了其先验分布的估计上,而在这一过程中我们又额外引入了 \( \alpha,\beta\) 两个超参数(参数 \(P\) 的参数),问题不仅没有减少反而增加了。对于如何处理超参数,在贝叶斯学派里又分化出了几种不同的看法和观念。
最正宗的就是 Hierarchical Bayes Inference。这种方案将超参数视为新的随机变量,再一次应用贝叶斯公式进行转化,从而可能引入超超参数,这一过程可以理论上一直做下去。不过这种方式,计算复杂,应用性差,大多数时候也没有必要。
还有一种是 Empirical Bayes Method。这种方案里,有很多不同的实现方式和算法。但基本精神就是,通过观测的试验数据来估计超参数和先验,这里引入试验数据来估计,就有了某种迭代的意味。经验贝叶斯因此有点像 frequentist 和 Bayenian 的结合。
当然我们也可以靠一些 argue 来直接把最适合该参数分布的超参数从开始就固定下来,也可以说成从一开始就选取一个合适的没有超参数的先验分布。这种做法又有两种信念,一种是先验分布(即超参数)完全取决于每个人的观察,知识和想法,没有一个统一的标准,从而计算出的后验概率完全取决于不同人的不同世界观,这种信念被称为主观贝叶斯。而另一种信念,被称为客观贝叶斯,认为至少在对于先验一无所知(uninformative)的时候,可以找到一种确定唯一的对于一无所知的先验分布描述(uninformative prior)。
无信息先验
在简单了解贝叶斯推断复杂的方法论和变种后,我们现在采取客观贝叶斯的方法论,鉴于我们对该试验的情况一无所知,我们现在要寻找一个代表了这种无知的不含超参数的先验分布。不过这件事情理论上还没能有足够完美和决定性的结论。既然是最无知的先验,自然的想法就是可能取值区间的均匀分布。但这件事情比较微妙,因为很多时候,坐标系一改变,原来的所谓均匀分布就不均匀了(比如运动速度的均匀分布不是到达时间的均匀分布)。另一方面,注意到这里的“一无所知”都有两种含义。一种是确定的知道该试验可以发生成功和失败两种结果,只不过对成功概率 \(P\) 的值完全没有先验知识。另一种含义是,对这一试验是否会发生成功和失败两种结果都一无所知(比如也许试验只是会一直失败)。这两种含义下的一无所知是完全不同的。对于前者,可以认为观察者至少有成功失败两次试验结果来 refine 先验了。确信试验既可以成功也可以失败就已经是很强的信息了。
Jaynes 说明了对于后者的无知先验(完全不知道该试验是否既可以成功也可以失败)公式为 \(\frac{1}{p(1-p)}\)。这一式子被称为 Haldane prior。Haldane prior 看起来很奇怪,注意到该先验在0到1区间的积分不收敛,因此该概率表达式并未归一化。但如果将这种积分不收敛的先验带入贝叶斯公式,最后可以得出有效的后验的话,这种发散的先验被称为 improper prior,由于我们的目的是计算合理的后验,这种发散先验也可以被用于推导而不会有什么本质的困难。
对于 Haldane prior 的进一步分析,显示对于完全无知的情况,先验应该集中在 \(P=0,1\) 两点,这表明完全未知的试验极大可能结果是确定的。另一方面考虑到成功失败两者的定义可以交换,完全无信息的先验自然应该关于 \(P=\frac{1}{2}\) 对称,Haldane prior 满足该对称性。 Haldane prior 可以视为 \(Beta(0,0)\) 分布,虽然该分布在 \(\alpha,\beta =0\) 这种情况系数也是发散的,但比例关系还在,类似 improper prior 还可以形式上用来推导。
如果想计算前者情形的无信息先验(已知该试验可以成功和失败),那么该先验就等价于第二种情形的先验补充上两次试验一次成功一次失败结果对该先验修正最终得出的后验。也即
\[P_1(P)=P(P\vert X_1=1,X_2=0)\propto P_2(P) P(X_1=1,X_2=0\vert P)=\frac{1}{P(1-P)}P(1-P)=1.\]也就是说,这样得出的第一种情形的先验恰好是在\((0,1)\)的均匀分布,这也和直觉相符,侧面验证了完全无知时使用 Haldane prior 的合理性。这一均匀分布也可形式上表示为 \(Beta(1,1)\)。
和概率问题各处细节的不同流派一样,无信息先验这点上,大家依旧无法达成统一。 Jefferys 通过对于完全无信息应该在重参数化下不改变这一性质(类似于我们根据对称性得出的 prior 应关于 \(P=1/2\) 对称的分析,不过稍复杂些,具体推导参见wiki),推导出完全没有已知信息的先验应该为 Jefferys Prior: \(Beta(\frac{1}{2},\frac{1}{2})\propto \frac{1}{\sqrt{p(1-p)}}\)。究竟哪种先验更无知的争论可能还没完全解决,甚至是否有所谓的完全无信息的先验都充满争议(可以参考这个回答)。这时也许就显示出了主观贝叶斯哲学的好处:先验就是个人的观点,你有多无知,得出什么形式的先验,不分对错,和我无关,无需争论,没必要一致,算出的后验靠不靠谱自己负责。
后验计算
既然知道了各种不同的先验,我们现在就可以来计算太阳升起的概率参数的后验分布了。为了方便,我们直接用共轭分布的性质来计算。一个事实:如果可能性是二项分布,那么先验是Beta分布则后验也是Beta分布,两个Beta分布被称为共轭分布。公式为(推导很简单,这里省略)
\[Beta(x+a,n-x+b)\propto Binomial(n,\theta)*Beta(a,b).\]注意到 \(Beta(\alpha,\beta)\) 分布的极值位置(众数)是 \(P_{max}=\frac{a-1}{(a-1)+(b-1)}\)。 这对于类似多次独立试验模型,先验就等价的相当于在数据集中添加了 \(a-1\) 次成功试验和 \(b-1\) 次失败试验。另一方面 Beta 分布的期望值为 \(E(P)=\frac{a}{a+b}\),所以如果用期望估计代替点估计,那么先验相当于在数据集添加了 \(a\) 次成功试验和 \(b\) 次失败试验。在这个意义上,上小节那些可以约化到 Beta 分布的先验,可以解释为在这一组试验前我们观察到的一系列试验情况(这也是先验是我们以前已知的信息相符)。比如 \(Beta(1,1)\) 这一先验,就表示了之前我们已经有了一次成功一次失败的经验,这恰好和我们上文的推导相一致。Beta 分布的参数 \(\alpha,\beta\) 越大,就代表我们之前的试验经验越多,我们的先验信息就越强。从这个意义上,Haldane prior \(Beta(0,0)\) 代表了无信息的先验也是合理的。
根据以上事实和公式,现在可以轻松计算出各种无信息先验形式的情况下,概率的分布以及根据点估计或均值估计给出的概率。
均匀分布:\(Beta(S+1,F+1)\),均值 \(\frac{S+1}{N+2}\),众数 \(\frac{S}{N}\)。
Jefferys prior: \(Beta(S+1/2,F+1/2)\),均值 \(\frac{S+1/2}{N+1}\),众数 \(\frac{S-1/2}{N-1}\)。
Haldane prior: \(Beta(S,F)\),均值 \(\frac{S}{N}\), 众数 \(\frac{S-1}{N-2}\)。
结果各不相同,不过好在当试验次数和成功次数 \(N,S\) 足够大时,以上不同先验不同方法给出的估计非常接近。你也可以计算每种先验对应后验的方差,来给出概率参数估值的准确程度。
一些反思
具体到太阳升起问题,\(S=0\),这时最经典的回答是拉格朗日的答案,对应我们这里现代版本中的均匀分布先验加期望值估计参数。也就是拉格朗日的出发点(先验)就默认了太阳升起和不升起都一定会发生,这一假设是否合理,我们只能按照主观贝叶斯的佛系哲学置之不顾了。即使不考虑先验,拉格朗日的日出概率问题的解也招致很多争议。如果我们认为明天日出的概率是 \(\frac{N+1}{N+2}\),那么如果明天日出后天日出的概率则为 \(\frac{N+2}{N+3}\),以此类推,接下来N天总是太阳升起的概率恰好为 \(\frac{1}{2}\)。比如观察者已经确定过去一万年太阳总是会升起,由此推断出未来一万年太阳总是会升起的概率是0.5,这么小概率显然是有悖常识的。这也是初始先验(表达了存在失败试验的信息)选择不好的结果。
另外,对于 Haldane 先验,在日出问题里,后验发散,先验无法作为 improper prior 计算。此时正确的概率估计,也许可以通过使用 \(Beta(\epsilon,\epsilon)\) 的分布作为先验来逼近,也可参考下文不放回情形的极限。
无放回取样版本
问题模型
刚才的独立试验,也可以看作有放回的取小球(二项分布),这节我们看一下无放回取小球的情况(超几何分布),两者在球总数足够多,成功球比例固定的极限下,是相同的。
这节的问题情境如下:黑箱里有14亿个小球,已知小球可能存在红球或白球(具体是否有白球不知道,也就是 Haldane prior)。假设严格按照分层随机抽样,选取了足够有代表性的2970个球,并且发现2970个球都是红色,问,黑箱里有多少白球?
简单推导
由于是 Haldane prior,我们取先验为 \(\frac{1}{S(N-S)}\),那么后验为(N 为总球数,n为取出的球数)
\[P(S|N,n,s)\propto {1 \over S(N-S)}{S \choose s}{N-S \choose n-s}\propto {S!(N-S)! \over S(N-S)(S-s)!(N-S-[n-s])!}.\]注意到 \(S\) 和 \(S!\) 的首项正好相消,这也是超几何分布可以考虑 \(S=0\) 情形的原因。归一化后成功球数的期望为
\[E\left({S \over N}|n,s=0,N\right)={1 \over N}\sum _{S=1}^{N-n}SP(S|N,n=1,s=0)={1 \over N}{\sum _{S=1}^{N-n}\prod _{j=1}^{n-1}(N-S-j) \over \sum _{R=1}^{N-n}{\prod_{j=1}^{n-1}(N-R-j) \over R}}.\]超几何分布的解析性质自然比二项分布差了许多。因此为了得到结果我们需要对阶乘(乘方估计)和求和(积分估计)做一些近似(参考wiki),N等于14亿这么大的数字,使得做了很多近似误差也不大。
最后我们有
\[E\left({S \over N}\vert n,s=0,N\right) \approx \frac{0.434}{n \log_{10}N}.\]对于我们的小红球问题,\(N\approx 1.4\times 10^9, n=2970\),代入上式可得,黑箱里白球数目的期望为
\[E(S)=N E({S \over N})\approx 22000.\]也就是说如果抓球真的是随机有代表性的,那么黑箱中非红色的球只有大概22000个左右。考虑2970个球中出现一个白球时,按频率主义,对应的黑箱中白球则是47万(其他方案给出估计数量级相同)。
当然如果这个试验每次取球没那么随机,就不是数学能处理的问题了,也不在我的理解力和想象力之内了。
聪明的你,读完这篇文字,也许看懂点了概率背后的数学和哲学,但看得懂为什么白球那么少吗?看不懂也没什么关系。因为,拉普拉斯告诉了你,明天的太阳,大概率还是照常升起。
EOF