一枚不公平的硬币–最大似然估计和贝叶斯方法

问题介绍

想象一下，我们正在进行一个预测连续抛硬币的过程。我们可以使用伯努利过程来建模。在伯努利过程中，一个随机变量有两种可能的结果: {0,1} (我们可以把0看作是 “正面” 而把1看成是 “反面”。如果我们抛出一枚硬币，我们可以将抛出结果的概率建模成为:

$$
𝑃（𝑥=1）= \mu
$$

同时，当𝑃(𝑥=0):

$$
𝑃(𝑥=0)=1-𝑃(𝑥=1)
$$

如果$\mu$=0.5，这说明抛硬币是公平的。当我们抛出$𝑁$=100次时，我们可能会得到一个包含正面和反面结果的序列。由于每次抛出的单个元素是不可预测的，但我们可以估计大约有50个正面和50个反面。如果我们看到25个正面和75个反面，我们就很难相信$\mu$=0.5而可能猜测它更接近于$\mu$=0.75。

具体问题

假设我们现在已有的信息是：$\mu$有五个可能的值：0, 0.25, 0.5, 0.75, 1.0。并且我们得到了有限次数下投掷硬币的观察结果，将其表示为$𝑁$。在这些观察的基础上，我们需要对未来继续抛掷硬币可能产生的序列进行建模。我们将使用两种方法：最大似然估计法和贝叶斯法。

最大似然估计

首先，我们记下我们连续抛硬币后观察到的一系列结果，并写下该特定投掷序列的概率作为$\mu$的函数。我们可以这样做的原因是，我们已经根据𝑃(𝑥=1∣$\mu$)=$\mu$；𝑃(𝑥=0∣$\mu$)=1-𝑃(𝑥=1∣$\mu$) 得到了单个结果的概率，并且每次单次抛硬币的结果是独立的。因此我们可以将单个抛掷结果的概率相乘来得出一系列的概率。

当连续抛掷硬币序列的$𝑁$=5时，如果我们假设，第一次抛掷的结果是1，其他四次是0：
$$
10000
$$
那么，这个特定结果的概率是:
$$
\mu(1-\mu)^4
$$
当结果序列为:
$$
01000
$$
概率同理可得也是:
$$
\mu(1-\mu)^4
$$
不难看出，对于观察到的$N$次抛掷的序列，其中如果观察到抛掷结果为”1”的数量为$N_1$，那么抛掷结果是”0”的数量就为
$N_0 = N - N_1$
，这些结果的概率是（无论1和0是以什么顺序出现的）:

$$
\lambda = \mu^{N_{1}}(1- \mu)^{N-N_1},
$$

在最大似然估计法中，我们将上述表达式视为一种似然。在这个函数中只有一个参数，那就是$\mu$。而我们的目标是找到$\mu$的值，使$\lambda$的值达到最大。为了方便进行计算，我们通常将其转化为使用log的形式:
$$
\ln N = N_1 \ln \mu + (N - N_1) \ln (1 - \mu)
$$
当$\lambda$的值取到最大时，即使得$\ln \lambda$的值取到最小。所以，我们将它求导后的值取到0时来求达到该最小值时的$\mu$的值。
$$
\frac{N_1}{\mu} - \frac{(N - N_1)}{1 - \mu} = 0
$$
继续推得:
$$
N_1(1 - \mu) = (N - N_1) \mu
$$
最终可得:
\begin{equation}
\mu = \frac{N_1}{N}
\end{equation}

将上式进行二次求导可得：
$$
-\frac{N_1}{\mu^2} - \frac{(N - N_1)}{(1 - \mu)^2}
$$
由于$N - N_1 > 0$ 并且 $N >0$，因此得到的这个结果为负。从而我们可以推出$\mu = \frac{N_1}{N}$，是这个对数似然估计的最小值，并且是唯一的最小值。因此，$\mu$取这个值时得到最大的似然估计值。

直观地说，这个对$\mu$取值的估计是有意义的。

贝叶斯方法

贝叶斯方法要求我们对$\mu$的先验值做出假设。在贝叶斯的观点中，概率与个人主观的相信程度相对应。首先我们先假设认为$\mu$取到五个值，我们必须对这些值定义一个概率分布，这个分布即为先验分布。例如，我们可以先假设:

$$
\begin{align}
P_{prior}(\mu = 0. ) & = 0.05 \\
P_{prior}(\mu = 0.25) & = 0.05 \\
P_{prior}(\mu = 0.50 ) & = 0.7 \\
P_{prior}(\mu = 0.75 ) & = 0.15 \\
P_{prior}(\mu = 1.0 ) & = 0.05
\end{align}
$$

我们的目标是根据我们所观察到的一连串抛掷硬币的情况来重新评估这些概率。通过贝叶斯法则，我们可以得到:
$$
P(\mu \mid D ) = \frac{P(D \mid \mu)}{P(D)} P_{prior}(\mu)
$$
这里$P(D \mid \mu)$是上面所假设的数据的可能性:
$$
P(D \mid \mu) = \mu^{N_1}(1 - \mu)^{N - N_1},
$$
鉴于$\mu$它已经给出了一个特定序列的概率，它可以被归纳为唯一的两个相关数字:抛掷结果为0的数量$N_0$,和结果为1的数量$N_1$。其总数为$N= N_0 + N_1$。
通过对一系列的连续投掷结果观察，我们可以有统计的数字结果来计算后验概率$P(\mu \mid D)$。

归一化的复杂性

通过归一化系数计算一个单一的$\mu$的后验概率是比较简单的:
$$
P_{posterior}(\mu = 0.5) \sim P( D \mid \mu)P_{prior}(\mu = 0.5)
$$
考虑到可能性的简单形式和我们知道先验概率的情况，这是一个轻量级的计算。但是为了正确地归一化，我们需要直接计算$\sum_\mu P(D \mid \mu)P(\mu)$或者估计所有可能结果的后验概率，并使用它来归一化后验分布。

在后文中会具体演示这两种方法。无论是哪种方式，要计算一个结果的后验概率，我们都必须计算所有的结果。即使在这样一个小表格上，仅仅计算一个整体的归一化系数就需要很大量的工作。特别是对于大量的结果和多变量似然分布，这个问题可能变得很严重。但是通过使用参数化分布，如高斯分布，有时候可以避免这个问题。后文将更详细地讨论这个问题。

例一：

已有大量的$\mu = 0.75$的观测值。计算最大似然估计和后验分布。

数据样本

我们使用二项分布$N =1$， $p = 0.75$来生成 $\mbox{Ber}( x \mid \mu = 0.75)$的100个观察值的样本。

最大似然估计

import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = [4, 4]
plt.rcParams['figure.dpi'] = 200
import numpy as np

from numpy.random import default_rng
rng = default_rng()

N_bin  = 1 # binomial for N=1 is Bernouilli
p    = 0.75
size = 100
sample = rng.binomial(N_bin, p, size)

print(sample)

n_ones = sample.sum() # the sum is adequate for a total count

mu_ml = float(n_ones)/size

print('mu_ml = ', mu_ml)

结果为：
[1 0 1 1 1 0 1 1 0 1 1 1 0 1 1 1 1 1 0 1 1 1 1 1 1 0 0 1 1 1 1 1 1 0 1 1 1
1 1 1 1 0 1 1 0 0 1 1 0 1 1 1 0 1 0 1 0 0 1 0 0 1 1 0 1 1 1 1 1 1 1 1 1 0
1 1 0 1 1 0 1 1 0 1 0 0 1 1 1 1 1 1 1 0 1 1 0 0 0 1]
mu_ml = 0.71

提示

这是对原来$\mu$的合理估计。重要的是建模者只看到数据，而不能看到如何生成数据。为了执行预测，建模者将使用类似的代码(他们并不知道示例是如何创建的)，他们使用$p = \mu_{ml}$
这是$\mu$的估计值。重复单元格将在结果中提供更多的信息。
用真实数据做的实验往往不能随意重复。在这种情况下，最大似然估计仍然是一个点估计。

后验分布

import matplotlib.pyplot as plt

mu_values = np.array([0., 0.25, 0.5, 0.75, 1.])
prior   = np.array([0.05, 0.05,   0.7, 0.15,  0.05])

print('If the prior is a proper probability distribution function, it should be normalised:', prior.sum())

# we use the same sample as generated above that was used for the ML estimate

ones  = n_ones
zeros = size - ones

def likelihood(mu, ones, zeros):
    # multiply performs element wise multiplication
    return np.multiply(np.power(mu,ones),np.power(1-mu,zeros))

llh = likelihood(mu_values,ones,zeros)
posterior = np.multiply(llh,prior)
normalised_posterior = posterior/posterior.sum()

plt.plot(mu_values, prior,'b*',label='prior')
plt.plot(mu_values,normalised_posterior,'r+',label='posterior')
plt.xlabel('$\mu$')
plt.legend()
plt.savefig('falsecoin.pdf')