从一个例子看频率学派与贝叶斯学派的不同（Python）

考虑如下的一个游戏场景：

setup：

Alice and Bob enter a room. Behind a curtain there is a billiard table, which they cannot see, but their friend Carol can. Carol rolls a ball down the table, and marks where it lands. Once this mark is in place, Carol begins rolling new balls down the table. If the ball lands to the left of the mark, Alice gets a point; if it lands to the right of the mark, Bob gets a point. We can assume for the sake of example that Carol’s rolls are unbiased: that is, the balls have an equal chance of ending up anywhere on the table. The first person to reach six points wins the game.

question：

In a particular game, after eight rolls, Alice has five points and Bob has three points. What is the probability that Bob will go on to win the game?

一个朴素（naïve）的频率学派的方法

确定最终的结果之前，我们需要对marker（标记）的位置有一个估计。我们将标记的位置量化为有利于Alice的概率$p$。已知8次中5次为Alice得分，使用频率的方法，我们很容易计算出$p$的最大似然估计（MLE：Maximum likelihood estimate）：

$$\hat{p}=5/8$$
紧接着，我们也很容易计算出Bob获胜的概率，此时需要Bob连赢三场（in a row）：
$$P(B)=(1-\hat{p})^3$$

p_hat = 5/8
        # 整数相除默认为float
        # 这是python3的语法支持
        # 5//8(python3) == 5/8(python2)
        # 5/8(python3) == 5./8(python2)
p_B = (1-p_hat)**3
print('Naïve frequentist probability of Bob wining: {0:.2f}'.format(p_B))

输出为：

Naïve frequentist probability of Bob wining: 0.05

根据几率公式：$o(p)=\frac{1-p}p$，我们计算Bob获胜的几率$o(P(B))$，

print('Odds against Bob wining: {0:.0f} to 1'.format((1-p_B)/p_B))

输出为：

Odds against Bob wining: 18 to 1

所谓几率18:1，也即Alice赢18次，而Bob只赢一次。

贝叶斯的方法

首先需要定义一些记号，考虑如下的随机变量（R.V.，Random Variable）：

• $B$ = Bob wins

• $D$ = observed data, i.e. $D=(n_A, n_B)=(5, 3)$

  贝叶斯方法的核心观点，将$D$视为已知数据，然后构造条件概率，而非仅仅得到一个简单的概率。

• $p$ = 未知的概率（球落在Alice一方的概率）

所以，贝叶斯的观点下问题变为求解：$P(B|D)$，也即给定当前观测数据（Alice赢五球，Bob赢3球）下Bob最终获胜的概率。

对待一些无关参数（nuisance parameters，如$P(B,p|D)$中的$p$），贝叶斯的思路是将之边际化（marginalization），或者对无关参数在其整个定义域上对联合分布进行积分。这里，我们需要首先计算联合分布：$P(B, p|D)$，然后使用如下的等式对$p$进行边际化：

$$P(B|D)\equiv \int_{-\infty}^{\infty}P(B,p|D)dp$$

又根据条件概率的推导，即$P(B,p|D)=P(B|p,D)P(p|D)$，得：

$$\begin{split} P(B|D)=&\int_{-\infty}^{\infty}P(B|p,D)P(p|D)dp\\ =&\int_{-\infty}^{\infty}P(B|p,D)\frac{P(p)P(D|p)}{P(D)}dp\\ =&\frac{\int P(B|p,D)P(p)P(D|p)dp}{p(D)}\\ =&\frac{\int P(B|p,D)P(p)P(D|p)dp}{\int P(D|p)P(p)dp} \end{split}$$

其中第一步推导用的是：条件概率的关系式，第二步是贝叶斯公式，第三步是提出与积分无关的项，第四步是对分母应用全概率公式。

这一系列转换、边缘化的最终目的是将不易计算的概率转换为容易计算的概率关系：

• $P(B|p,D)$：给定$p$，且在给定$D$（$(n_A,n_B)=(5, 3)$）条件下，Bob获胜的概率，$P(B|p,D)=(1-p)^3$。

• $P(D|p)$，也是一个容易计算的项，$P(D|p)\propto p^5(1-p)^3$

• $P(p)$，概率$p$的先验概率。$P(p)\propto 1$

所以，最终的计算结果为：

$$P(B|D)=\frac{\int_0^1(1-p)^6p^5dp}{\int_0^1(1-p)^3p^5dp}$$

积分并不好计算，好在有现成的$\beta$ function：

$$\beta(n, m)=\int_0^1(1-p)^{n-1}p^{m-1}dp$$ 。

>>> from scipy.special import beta
>>> bayes_p_B = beta(6+1, 5+1)/beta(3+1, 5+1)
>>> print('P(B|D) = {0:.2f}'.format(bayes_p_B))

输出为：

P(B|D) = 0.09

其相关的几率为：

print('Bayesian odds against Bob wining: {0:.0f}'.format((1-bayes_p_B)/bayes_p_B))

输出为：

Bayesian odds against Bob wining: 10 to 1

仿真验证（Monte Carlo）

p = np.random.rand(100000)
rolls = np.random.rand(11, len(p))

Alice_wins = rolls < p
Bob_wins = rolls >= p
total_wins = Alice_wins + Bob_wins
assert np.all(total_wins == 1)
print('sanity check passed')

good_games = np.sum(Alice_wins[0:8, :], 0) == 5
print('the respective probs:', p[good_games])
print('# of suitable games: {0}'.format(good_name.sum()))

# truncate our results to consider only these good games
# 表达一种条件概率的概念
given_Alice_wins = Alice_wins[:, good_games]
given_Bob_wins = Bob_wins[:, good_games]

# Monte Carlo Prob
P_B_mc = (np.sum(given_Bob_wins, 0) == 6).sum() / good_games.sum()
print('Monte Carlo probability of Bob wining: {:.2f}'.format(P_B_mc))
print('MC odds against Bob wining :{:.0f}'.format((1-P_B_mc)/P_B_mc))

输出发现，结果更接近于贝叶斯的方案。