PRML：二元变量分布

伯努利分布

考虑二元随机变量 x∈{0,1}（抛硬币，正面为 1，反面为 0），其概率分布由参数 μ 决定：

p(x=1)=μ

其中 (0≤μ≤1)，并且有 p(x=0)=1−μ。这就是伯努利分布（Bernoulli distribution），其概率分布可以写成：

Bern(x|μ)=μx(1−μ)1−x

均值和方差为：

E[x]var[x]=μ=μ(1−μ)

伯努利分布的最大似然估计

考虑一组 x 的观测数据 D={x1,…,xN}，在独立同分布的假设下，其似然函数为

p(D|μ)=∏n=1Np(xn|μ)=∏n=1Nμxn(1−μ)1−xn

对数似然为

lnp(D|μ)=∑n=1Nlnp(xn|μ)=∑n=1N{xnlnμ+(1−xn)ln(1−μ)}

对数似然值只依赖于 ∑Nn=1xn 的取值，而事实上 ∑Nn=1xi 就是伯努利分布的一个充分统计量，它可以提供参数 μ 的全部信息。

对 μ 最大化对数似然，我们很容易得到

μML=1N∑n=1Nxn

即最大似然估计值为样本均值（sample mean），若样本中 x=1 的数目为 m 则：

μML=mN

考虑抛三次硬币出现了三次正面的情况，此时 N=m=3,μML=1。在这种情况下，最大似然估计会得到每次都是正面的结果，这显然违背了我们的正常认知。事实上，这是一种过拟合的典型表现。

为了解决这个问题，我们可以考虑引入先验知识。

二项分布

给定数据总数 N，x=1 的总次数 m 满足一定的分布，这个分布叫做二项分布（binomial distribution）。

从伯努利分布的似然函数中可以看出它应该正比于 μm(1−μ)N−m，事实上它可以写成：

Bin(m | N,μ)=(Nm)μm(1−μ)N−m

其中

(Nm)≡N!(N−m)!m!

是组合数。

验证它是一个概率分布，二项式定理给出：

∑m=0N(Nm)μm(1−μ)N−m=(μ+1−μ)N=1

例如下面 N=20,μ=0.6 的一个分布直方图。

import matplotlib.pyplot as plt
import numpy as np
import scipy as sp
%matplotlib inline
from scipy.stats import binom

n = 20
mu = 0.6

X = binom.rvs(n, mu, size=10000)
fig, ax = plt.subplots()
ax.hist(X, bins=range(21), rwidth=0.7)
ax.set_xlabel("$m$", fontsize='x-large')
# ax.set_xlim(0, 21)
# ax.set_yticks(np.arange(0, 0.31, 0.1))
# ax.set_xticks(np.arange(0.5, 10.6, 1))
# ax.set_xticklabels(range(21))
ax.set_title(r'$N = 20, \mu=0.6$', fontsize='x-large')
plt.show()

计算均值时考虑下式对 μ 的导数，方差考虑对 μ 的二阶导数，其均值和方差分别为：

E[m]var[m]=Nμ=Nμ(1−μ)

beta 分布

之前看到，当数据量较少时，最大似然的结果很可能会过拟合。为了减少这样的情况，从 Bayes 概率的观点出发，我们引入一个关于 μ 的先验分布 p(μ)。

我们观察到似然函数是一系列 μx(1−μ)1−x 形式的乘积，如果我们选择一个正比于 μ 的某个幂次和 1−μ 的某个幂次的先验分布，那么对 μ 来说，后验分布应当满足同样的形式。这样的性质叫做共轭性（conjugacy）。

在这里，我们引入的是 0−1 间的 beta 分布：

Beta(μ | a,b)=Γ(a+b)Γ(a)Γ(b)μa−1(1−μ)b−1

其中：

Γ(x)≡∫∞0ux−1e−udu

满足如下性质：

Γ(x+1)Γ(1)=∫∞0uxe−udu=[−e−uux]∞0+x∫∞0ux−1e−udu=0+xΓ(x)=xΓ(x)=∫∞0e−udu=[−e−u]∞0=1

验证它是一个概率分布：

由定义

Γ(a)Γ(b)=∫∞0xa−1e−xdx∫∞0yb−1e−ydy

令 t=y+x,dt=dy，则有：

Γ(a)Γ(b)=∫∞0xa−1{∫∞x(t−x)b−1e−tdt}dx

交换积分次序，原来 t 是从 x 积分到 ∞，现在 x 是从 0 积分到 t：

Γ(a)Γ(b)=∫∞0∫t0xa−1(t−x)b−1e−tdxdt

令 x=tμ,dx=tdμ，则有

Γ(a)Γ(b)=∫∞0e−tta−1tb−1tdt∫10μa−1(1−μ)b−1dμ=Γ(a+b)∫10μa−1(1−μ)b−1dμ

于是我们有：

∫10Beta(μ | a,b)dμ=1

其均值和方差为

E[μ]var[μ]=aa+b=ab(a+b)2(a+b+1)

求均值：

从归一化我们知道：

∫10μa−1(1−μ)b−1dμ=Γ(a)Γ(b)Γ(a+b)

从而，利用 Γ(x+1)=xΓ(x)：

E[μ]=Γ(a+b)Γ(a)Γ(b)∫μa+1−1(1−μ)b−1dμ=Γ(a+b)Γ(a)Γ(b)Γ(a+1)Γ(b)Γ(a+b+1)=aa+b

类似可以得到：

E[μ2]=Γ(a+b)Γ(a)Γ(b)∫μa+2−1(1−μ)b−1dμ=Γ(a+b)Γ(a)Γ(b)Γ(a+2)Γ(b)Γ(a+b+2)=a(a+1)(a+b)(a+b+1)

从而可以计算出方差。

a 和 b 叫做超参数，因为它们控制分布的参数 μ。

from scipy.stats import beta

fig, axes = plt.subplots(2, 2,figsize=(10, 7))

axes = axes.flatten()

A = (0.1, 1, 2, 8)
B = (0.1, 1, 3, 4)

xx = np.linspace(0, 1, 100)

for a, b, ax in zip(A, B, axes):
    yy = beta.pdf(xx, a, b)
    ax.plot(xx, yy, 'r')
    ax.set_ylim(0, 3)

    ax.set_xticks([0, 0.5, 1])
    ax.set_xticklabels(["$0$", "$0.5$", "$1$"], fontsize="large")
    ax.set_yticks([0, 1, 2, 3])
    ax.set_yticklabels(["$0$", "$1$", "$2$", "$3$"], fontsize="large")
    ax.set_xlabel("$\mu$", fontsize="x-large")

    ax.text(0.1, 2.5, r"$a={}$".format(a), fontsize="x-large")
    ax.text(0.1, 2.2, r"$b={}$".format(b), fontsize="x-large")

有了先验分布，我们的后验分布为

p(μ | m,l,a,b)∝μm+a+1(1−μ)l+b−1

其中 l=N−m。

由共轭性，我们知道后验分布还是一个 beta 分布，从而有

p(μ | m,l,a,b)∼Beta(μ | a+m,b+l)

至此，我们看到，如果我们观测到了一组 m 个 x=1 和 l 个 x=0 的数据，那么超参由 a,b 变成 a+m,b+l。因此，超参 a,b 可以看成是 x=1 和 x=0 的有效观测次数（注意：a,b 可以不是整数）。

更进一步，当有新的数据到来时，后验概率可以看成新的先验概率，因此这个过程可以序列化进行。下图表示观测到一个新的 x=1 数据时，先验到后验的变化。

xx = np.linspace(0, 1, 100)

fig, axes = plt.subplots(1, 3, figsize=(10, 2))

axes = axes.flatten()

axes[0].plot(xx, beta.pdf(xx, 2, 2), 'r')
axes[0].set_ylim(0, 2)
axes[0].text(0.1, 1.6, "prior", fontsize="x-large")
axes[0].set_xlabel("$\mu$", fontsize="x-large")
axes[0].set_xticks([0, 0.5, 1])
axes[0].set_yticks([0, 1, 2])

axes[1].plot(xx, xx)
axes[1].set_ylim(0, 2)
axes[1].text(0.1, 1.6, "likelihood", fontsize="x-large")
axes[1].set_xlabel("$\mu$", fontsize="x-large")
axes[1].set_xticks([0, 0.5, 1])
axes[1].set_yticks([0, 1, 2])

axes[2].plot(xx, beta.pdf(xx, 3, 2), 'r')
axes[2].set_ylim(0, 2)
axes[2].text(0.1, 1.6, "posterior", fontsize="x-large")
axes[2].set_xlabel("$\mu$", fontsize="x-large")
axes[2].set_xticks([0, 0.5, 1])
axes[2].set_yticks([0, 1, 2])

plt.show()

如果我们的目的是预测下一次实验的结果，那么我们有

p(x=1|D)=∫10p(x=1|μ)p(μ|D)dμ∫10μp(μ|D)=E[μ|D]

而我们知道后验概率的分布，所以可以计算出其均值：

p(x=1|D)=m+am+a+l+b=m+aN+a+b

当 m,l→∞ 时，有

p(x=1|D)=m+am+a+l+b=m+aN+a+b→mN

即趋向于最大似然的解。当 N 有限时，我们的解总在先验均值和最大似然解之间。

另外我们观察到，随着 (a,b) 的增加，函数分布变得越来越尖，即方差越来越小。事实上，随着观测数据的增加，我们对于后验分布的不确定性也在不断减小。

另外考虑条件期望和方差的性质，我们有：

Eθ[θ]=ED[Eθ[θ|D]]

和

varθ[θ]=ED[varθ[θ|D]]+varD[Eθ[θ|D]]

从而我们知道，从平均意义上来说，后验分布的方差要比先验分布的方差要小。后验分布在数据分布上的均值等于先验分布的均值