考虑二元随机变量 $x\in \{0,1\}$(抛硬币,正面为 1,反面为 0),其概率分布由参数 $\mu$ 决定:
$$ p(x=1)=\mu $$其中 $(0 \leq\mu \leq 1)$,并且有 $p(x=0)=1-\mu$。这就是伯努利分布(Bernoulli distribution
),其概率分布可以写成:
均值和方差为:
$$ \begin{align} \mathbb E[x] & = \mu \\ \text{var}[x] & = \mu(1-\mu) \end{align} $$考虑一组 $x$ 的观测数据 $\mathcal D = \{x_1, \dots, x_N\}$,在独立同分布的假设下,其似然函数为
$$ p(\mathcal D|\mu) = \prod_{n=1}^N p(x_n|\mu) = \prod_{n=1}^N \mu^{x_n}(1-\mu)^{1-x_n} $$对数似然为
$$ \ln p(\mathcal D|\mu) = \sum_{n=1}^N \ln p(x_n|\mu) = \sum_{n=1}^N \left\{x_n\ln \mu + (1-x_n)\ln (1-\mu)\right\} $$对数似然值只依赖于 $\sum_{n=1}^N x_n$ 的取值,而事实上 $\sum_{n=1}^N x_i$ 就是伯努利分布的一个充分统计量,它可以提供参数 $\mu$ 的全部信息。
对 $\mu$ 最大化对数似然,我们很容易得到
$$ \mu_{ML} = \frac 1 N \sum_{n=1}^N {x_n} $$即最大似然估计值为样本均值(sample mean
),若样本中 $x=1$ 的数目为 $m$ 则:
考虑抛三次硬币出现了三次正面的情况,此时 $N=m=3, \mu_{ML} = 1$。在这种情况下,最大似然估计会得到每次都是正面的结果,这显然违背了我们的正常认知。事实上,这是一种过拟合的典型表现。
为了解决这个问题,我们可以考虑引入先验知识。
给定数据总数 $N$,$x=1$ 的总次数 $m$ 满足一定的分布,这个分布叫做二项分布(binomial distribution
)。
从伯努利分布的似然函数中可以看出它应该正比于 $\mu^{m}(1-\mu)^{N-m}$,事实上它可以写成:
$$ \text{Bin}(m~|~N,\mu) = \begin{pmatrix} N \\m \end{pmatrix} \mu^{m}(1-\mu)^{N-m} $$其中
$$ \begin{pmatrix} N \\m \end{pmatrix} \equiv \frac{N!}{(N-m)!m!} $$是组合数。
验证它是一个概率分布,二项式定理给出:
$$ \sum_{m=0}^N \begin{pmatrix} N \\m \end{pmatrix} \mu^{m}(1-\mu)^{N-m} = (\mu + 1 - \mu)^N = 1 $$下图是 $N = 10, \mu=0.25$ 的分布的一个直方图示意。
import matplotlib.pyplot as plt
import numpy as np
import scipy as sp
%matplotlib inline
from scipy.stats import binom
n = 10
mu = 0.25
yy = binom.rvs(n, mu, size=1000)
fig, ax = plt.subplots()
ax.hist(yy, bins=range(11), normed=True, rwidth=0.8)
ax.set_xlabel("$m$", fontsize='x-large')
ax.set_xlim(0, 11)
ax.set_yticks(np.arange(0, 0.31, 0.1))
ax.set_xticks(np.arange(0.5, 10.6, 1))
ax.set_xticklabels(range(11))
ax.set_title(r'$N = 10, \mu=0.25$', fontsize='x-large')
plt.show()
其均值和方差分别为:
$$ \begin{align} \mathbb E[m] & = N\mu \\ \text{var}[m] & = N\mu(1-\mu) \end{align} $$计算均值时考虑下式对 $\mu$ 的导数,方差考虑对 $\mu$ 的二阶导数:
$$ \sum_{m=0}^N \begin{pmatrix} N \\m \end{pmatrix} \mu^{m}(1-\mu)^{N-m} = 1 $$之前看到,当数据量较少时,最大似然的结果很可能会过拟合。为了减少这样的情况,从 Bayes 概率的观点出发,我们引入一个关于 $\mu$ 的先验分布 $p(\mu)$。
我们观察到似然函数是一系列 $\mu^x(1-\mu)^{1-x}$ 形式的乘积,如果我们选择一个正比于 $\mu$ 的某个幂次和 $1-\mu$ 的某个幂次的先验分布,那么对 $\mu$ 来说,后验分布应当满足同样的形式。这样的性质叫做共轭性(conjugacy
)。
在这里,我们引入的是 $0-1$ 间的 beta
分布:
其中:
$$ \Gamma(x) \equiv \int_0^{\infty} u^{x-1}e^{-u} du $$满足如下性质:
$$ \begin{align} \Gamma(x+1) & = \int_0^{\infty} u^{x}e^{-u} du = \left[-e^{-u}u^x\right]_0^{\infty} + x \int_0^{\infty} u^{x-1}e^{-u} du = 0 + x\Gamma(x) = x\Gamma(x) \\ \Gamma(1) & = \int_0^{\infty} e^{-u} du = \left[-e^{-u}\right]_0^{\infty} = 1 \end{align} $$验证它是一个概率分布:
由定义
$$ \Gamma(a)\Gamma(b) = \int_0^\infty x^{a-1}e^{-x} dx \int_0^\infty y^{b-1}e^{-y} dy $$令 $t = y + x, dt = dy$,则有:
$$ \Gamma(a)\Gamma(b) = \int_0^\infty x^{a-1} \left\{\int_x^\infty (t-x)^{b-1}e^{-t} dt \right\} dx $$交换积分次序,原来 $t$ 是从 $x$ 积分到 $\infty$,现在 $x$ 是从 $0$ 积分到 $t$:
$$ \Gamma(a)\Gamma(b) = \int_0^\infty \int_0^t x^{a-1} (t-x)^{b-1}e^{-t} dxdt $$令 $x = t\mu, dx = td\mu$,则有
$$ \begin{align} \Gamma(a)\Gamma(b) & = \int_0^\infty e^{-t} t^{a-1} t^{b-1} tdt \int_0^1 \mu^{a-1} (1-\mu)^{b-1} d\mu \\ & = \Gamma(a + b) \int_0^1 \mu^{a-1} (1-\mu)^{b-1} d\mu \end{align} $$于是我们有:
$$ \int_0^1 \text{Beta}(\mu~|~a,b) d\mu = 1 $$其均值和方差为
$$ \begin{align} \mathbb E[\mu] & = \frac{a}{a+b} \\ \text{var}[\mu] & = \frac{ab}{(a+b)^2(a+b+1)} \end{align} $$求均值:
从归一化我们知道:
$$ \int_0^1 \mu^{a-1} (1-\mu)^{b-1} d\mu = \frac{\Gamma(a)\Gamma(b)}{\Gamma(a + b)} $$从而,利用 $\Gamma(x+1) = x\Gamma(x)$:
$$ \mathbb E[\mu] = \frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)} \int \mu^{a+1-1} (1-\mu)^{b-1} d\mu = \frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)} \frac{\Gamma(a+1)\Gamma(b)}{\Gamma(a+b+1)} = \frac{a}{a+b} $$类似可以得到:
$$ \mathbb E[\mu^2] = \frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)} \int \mu^{a+2-1} (1-\mu)^{b-1} d\mu = \frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)} \frac{\Gamma(a+2)\Gamma(b)}{\Gamma(a+b+2)} = \frac{a(a+1)}{(a+b)(a+b+1)} $$从而可以计算出方差。
$a$ 和 $b$ 被叫做超参,因为它们控制分布的参数 $\mu$。
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(\mu~|~m,l,a,b) \propto \mu^{m+a+1} (1-\mu)^{l+b-1} $$其中 $l = N - m$。
由共轭性,我们知道后验分布还是一个 beta
分布,从而有
至此,我们看到,如果我们观测到了一组 $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, 'b')
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|\mathcal D) = \int_{0}^1 p(x=1|\mu) p(\mu|\mathcal D) d\mu \int_0^1 \mu p(\mu|\mathcal D) = \mathbb E[\mu|\mathcal D] $$而我们知道后验概率的分布,所以可以计算出其均值:
$$ p(x=1|\mathcal D) = \frac{m+a}{m+a+l+b} = \frac{m+a}{N+a+b} $$当 $m, l \to \infty$ 时,有
$$ p(x=1|\mathcal D) = \frac{m+a}{m+a+l+b} = \frac{m+a}{N+a+b} \to \frac{m}{N} $$即趋向于最大似然的解。当 $N$ 有限时,我们的解总在先验均值和最大似然解之间。
另外我们观察到,随着 $(a,b)$ 的增加,函数分布变得越来越尖,即方差越来越小。事实上,随着观测数据的增加,我们对于后验分布的不确定性也在不断减小。
另外考虑条件期望和方差的性质,我们有:
$$ \mathbb E_\mathbf\theta[\mathbf\theta] = \mathbb E_\mathcal{D}[\mathbb E_\mathbf\theta[\mathbf\theta|\mathcal D]] $$和
$$ \text{var}_\mathbf\theta[\mathbf\theta] = \mathbb E_\mathcal{D}[\text{var}_\mathbf\theta[\mathbf\theta|\mathcal D]] + \text{var}_\mathcal{D}[\mathbb E_\mathbf\theta[\mathbf\theta|\mathcal D]] $$从而我们知道,从平均意义上来说,后验分布的方差要比先验分布的方差要小。后验分布在数据分布上的均值等于先验分布的均值