概率论在解决模式识别问题时具有重要作用,它是构成更复杂模型的基石。
概率分布的一个作用是在给定有限次观测x1, . . . , xN的前提下,对随机变量x的概率分布p(x)建模。这个问题被称为密度估计(density estimation)。选择一个合适的分布与模型选择的问题相关,这是模式识别领域的一个中心问题。
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
考虑一个二元随机变量$x\in\{0,1\}$,$x=1$的概率被记作参数$\mu$,因此$p(x=1|\mu)=\mu$,$p(x=0|\mu)=1-\mu$,其中$\mu\leq1$。
x的概率分布可以写成$Bern(x|\mu)=\mu^x(1-\mu)^{1-x}$,称为伯努利分布($Bernoulli$ $distribution$)。
伯努利分布的均值和方差分别为$\mathbb{E}[x]=\mu$,$var[x]=\mu(1-\mu)$。
假设有x的观测值数据集$\mathbb{D}={x_1,...,x_N}$。假设每次观测都是独立地从$p(x|\mu)$中抽取的,因此关于$\mu$的似然函数是$p(\mathbb{D}|\mu)=\prod\limits_{n=1}^N p(x_n|\mu)=\prod\limits_{n=1}^N \mu^{x_n}(1-\mu)^{1-x_n}$。
其最大似然的估计值是$\mu_{ML}=\frac{1}{N}\sum\limits_{n=1}^N x_n=\frac{m}{N}$,如果把数据集里$x=1$的观测的数量记为m。
因此,在最大似然的框架中,x=1的概率是数据集中x=1的观测所占的比例。假设扔一枚硬币5次,碰巧5次硬币都正面朝上,那么$\mu_{ML}=1$。最大似然的结果会预测所有未来的观测值都是正面向上。这是最大似然中过拟合现象的一个极端例子。我们通过引入$\mu$的先验分布,来得到更加合理的解释。
from scipy.stats import bernoulli
伯努利分布直方图,u=0.6
u = 0.6
x = np.arange(2)
bern = bernoulli.pmf(x, u)
plt.bar(x, bern, facecolor='green', alpha=0.5)
plt.xlabel('x')
plt.xlim(-0.2, 2)
plt.title('Histogram of the nernoulli distribution')
plt.show()
给定数据集规模N的条件下,x = 1的观测出现的数量m的概率分布,称为二项分布(binomial distribution)。
$Bin(m|N, \mu)=\binom{N}{m}\mu^m(1-\mu)^{N-m}$,其中$\binom{N}{m}=\frac{N!}{(N-m)!m!}$。
其均值和方差分别是$\mathbb{E}[m]=\sum\limits_{m=0}^{N}mBin(m|N,\mu)=N\mu$,$var[m]=\sum\limits_{m=0}^{N}(m-\mathbb{E}[x])^2Bin(m|N,\mu)=N\mu(1-\mu)$。
from scipy.stats import binom
N = 10
u = 0.25
m = np.arange(N+1)
binData = binom.pmf(m, N, u)
binData
array([ 5.63135147e-02, 1.87711716e-01, 2.81567574e-01, 2.50282288e-01, 1.45998001e-01, 5.83992004e-02, 1.62220001e-02, 3.08990479e-03, 3.86238098e-04, 2.86102295e-05, 9.53674316e-07])
关于m的函数的二项分布直方图,其中N=10且u=0.25
plt.bar(m, binData, facecolor='green', alpha=0.5)
plt.xlabel('m')
plt.xlim(-0.2, 10.5)
plt.title('Histogram of the binomial distribution')
plt.show()
伯努利分布的似然函数是某个因子与$\mu^x(1-\mu)^{1-x}$的乘积的形式,如果我们选择一个正比于$\mu$和$(1-\mu)$的幂指数的先验概率分布,那么后验概率分布(正比于先验和似然函数的乘积)就会有着与先验分布相同的函数形式。这种性质叫做共轭性(conjugacy)。
Beta分布定义如下:
$$Beta(\mu|a, b) = \frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)}\mu^{a-1}(1-\mu)^{b-1}$$$\Gamma(x)$是Gamma函数,保证Beta分布归一化。
Beta分布的均值和方差为$$\mathbb{E}[\mu]=\frac{a}{a+b},var[\mu]=\frac{ab}{(a+b)^2(a+b+1)}$$
参数a和b被称为超参数(hyperparameter),因为它们控制了参数$\mu$的概率分布。
from scipy.stats import beta
fig = plt.figure(figsize=(12,8))
ax = fig.add_subplot(221)
x = np.linspace(0, 1, 100)
pbeta = beta.pdf(x, 0.1, 0.1)
ax.plot(x, pbeta, label='a=0.1\nb=0.1',lw=2)
plt.ylim(0,3)
plt.xlabel('$\mu$')
plt.grid(True)
plt.legend(loc="best")
ax = fig.add_subplot(222)
x = np.linspace(0, 1, 100)
pbeta = beta.pdf(x, 1, 1)
ax.plot(x, pbeta, 'g', label='a=1\nb=1',lw=2)
plt.xlabel('$\mu$')
plt.ylim(0,3)
plt.grid(True)
plt.legend(loc="best")
ax = fig.add_subplot(223)
x = np.linspace(0, 1, 100)
pbeta = beta.pdf(x, 2, 3)
ax.plot(x, pbeta, 'r', label='a=2\nb=3',lw=2)
plt.xlabel('$\mu$')
plt.ylim(0,3)
plt.grid(True)
plt.legend(loc="best")
ax = fig.add_subplot(224)
x = np.linspace(0, 1, 100)
pbeta = beta.pdf(x, 8, 4)
ax.plot(x, pbeta, 'm', label='a=8\nb=4',lw=2)
plt.xlabel('$\mu$')
plt.ylim(0,3)
plt.grid(True)
plt.legend(loc="best")
plt.show()
D:\Python27_32\lib\site-packages\matplotlib\transforms.py:656: RuntimeWarning: invalid value encountered in absolute inside = ((abs(dx0 + dx1) + abs(dy0 + dy1)) == 0)
将Beta先验和二项似然函数相乘,得到后验概率分布的形式为:
$$p(\mu|m,l,a,b)\propto\mu^{m+a-1}(1-\mu)^{l+b-1}$$其中$l=N-m$对应硬币反面朝上的样本数量。
具体得到归一化系数,公式如下:
$$p(\mu|m,l,a,b)=\frac{\Gamma(m+a+l+b)}{\Gamma(m+a)\Gamma(l+b)}\mu^{m+a-1}(1-\mu)^{l+b-1}$$我们可以发现,如果数据集里有m次观测为x=1,有l次观测为x=0,那么从先验概率到后验概率,a的值变大了m,b的值变大了l。
from scipy.misc import comb
pbinom = lambda m, N, p: comb(N, m) * p**m * (1-p)**(N-m)
上面代码是二项式分布的分布函数$Bin(m|N, p)=\binom{N}{m}p^m(1-p)^{N-m}$,其似然函数和该分布函数具有相同的形式,只是其参数p是未知的。
fig = plt.figure(figsize=(15,3))
ax = fig.add_subplot(131)
x = np.linspace(0, 1, 100)
pbeta = beta.pdf(x, 2, 2)
ax.plot(x, pbeta,lw=2)
plt.xlim(0,1)
plt.ylim(0,2)
plt.xlabel('$\mu$')
plt.grid(True)
plt.text(0.05, 1.7, "prior", fontsize=15)
ax = fig.add_subplot(132)
x = np.linspace(0, 1, 100)
pbin = pbinom(1, 1, x)
ax.plot(x, pbin,lw=2)
plt.xlim(0,1)
plt.ylim(0,2)
plt.xlabel('$\mu$')
plt.grid(True)
plt.text(0.05, 1.7, "likelihood function", fontsize=15)
ax = fig.add_subplot(133)
x = np.linspace(0, 1, 100)
post = pbeta * x
ax.plot(x, post,lw=2)
ax.plot(x, beta.pdf(x, 3, 2), lw=2)
plt.xlim(0,1)
plt.ylim(0,2)
plt.xlabel('$\mu$')
plt.grid(True)
plt.text(0.05, 1.7, "posterior", fontsize=15)
plt.show()
上面第三幅图中,绿色的线是实际的后验概率的曲线,蓝色曲线是直接将先验概率和似然函数相乘后,没有进行归一化的结果。
图中,beta先验分布的参数是a=2,b=2,似然函数N=m=1,对应获得一个x=1的样本数据,修正的beta后验分布参数是a=3,b=2