相关资料:《Pattern Recognition and Machine Learning》 by Christopher M. Bishop
以“贝叶斯原理在硬币实验中的应用”为背景,对其进行代码实现,以助于对相关理论部分的理解
原理
掷硬币问题的过拟合
掷硬币问题描述
非常经典的掷硬币情景:将一枚硬币连续掷N次,每次的结果都会是正面朝上或反面朝上,将这些结果全部记录下来。现在要做的,是仅根据这些结果,来计算硬币正面朝上的概率
最大似然估计法的缺陷
撇开生活经验不谈,假定硬币正面朝上的概率应该为 μ \mu μ,且记硬币朝上时,观测值 x = 1 x=1 x=1,反之,观测值 x = 0 x=0 x=0。则每次观测值的分布可以记录为如下形式(Bernoulli分布):
B e r n ( x ∣ μ ) = μ x ( 1 − μ ) 1 − x Bern(x|\mu)=\mu^x(1-\mu)^{1-x} Bern(x∣μ)=μx(1−μ)1−x
∣ | ∣符号右边的 μ \mu μ表示这是一个控制分布的参数
现在,假定对此硬币进行了N次观测(掷了N次),得到数据集(观测值的集合) D = { x 1 , … , x N } \mathcal D=\{x_1,\ \dots,\ x_N\} D={x1, …, xN},则根据各次观测之间相互独立的特性,可以计算出数据集处于各种状态(即对于各种可能的N个观测值的组合)的概率:
p ( D ∣ μ ) = Π n = 1 N p ( x n ∣ μ ) p(\mathcal D |\mu)=\mathop\Pi\limits_{n=1}^Np(x_n|\mu) p(D∣μ)=n=1ΠNp(xn∣μ)
两边取对数,得到:
ln p ( D ∣ μ ) = ∑ n = 1 N ln p ( x n ∣ μ ) = ∑ n = 1 N { x n ln μ + ( 1 − x n ) ln ( 1 − μ ) } \ln p(\mathcal D|\mu)=\sum\limits_{n=1}^N\ln p(x_n|\mu)=\sum\limits_{n=1}^N\{x_n\ln\mu + (1-x_n)\ln(1-\mu)\} lnp(D∣μ)=n=1∑Nlnp(xn∣μ)=n=1∑N{xnlnμ+(1−xn)ln(1−μ)}
试图找到上述 ln p ( D ∣ μ ) \ln p(\mathcal D|\mu) lnp(D∣μ)的最大值,即求出满足下列等式的 μ M L \mu_{ML} μML:
∂ ln p ( D ∣ μ ) ∂ μ = 0 \dfrac{\partial\ln p(\mathcal D|\mu)}{\partial \mu}=0 ∂μ∂lnp(D∣μ)=0
得到: μ M L = 1 N ∑ n = 1 N x n \mu_{ML}=\dfrac{1}{N}\sum\limits_{n=1}^N x_n μML=N1n=1∑Nxn. 对于上述提及的掷硬币问题,假定有m次是正面朝上,则最后的结果会是 μ M L = m N \mu_{ML}=\dfrac{m}{N} μML=Nm。这便是使用极大似然估计法,计算估测的硬币正面朝上的概率
这个估算出来的概率,在N无限大的时候是准确的,但在N非常小的时候就会出现问题。例如:N=3,连续三次硬币都是正面朝上,则 μ M L = 1 \mu_{ML}=1 μML=1,这显然是不合理的。这是一种过拟合的问题
使用贝叶斯思想进行改进
贝叶斯思想可以表示为:
p o s t e r i o r ∝ l i k e l i h o o d × p r i o r posterior\propto likelihood\times prior posterior∝likelihood×prior
即后验概率正比于似然度与先验概率的乘积。在这里,参数 μ \mu μ对应的先验概率分布为 p ( μ ) p(\mu) p(μ)
在掷硬币问题中,受Bernoulli分布的启发,将先验概率的形式设为:
B e t a ( μ ∣ a , b ) = Γ ( a + b ) Γ ( a ) Γ ( b ) μ a − 1 ( 1 − μ ) b − 1 Beta(\mu|a,b)=\dfrac{\Gamma (a+b)}{\Gamma(a)\Gamma(b)}\mu^{a-1}(1-\mu)^{b-1} Beta(μ∣a,b)=Γ(a)Γ(b)Γ(a+b)μa−1(1−μ)b−1
其中, Γ \Gamma Γ为gamma函数,其定义为:
Γ ( x ) = ∫ 0 ∞ u x − 1 e − u d u \Gamma(x)=\int_0^\infty u^{x-1}e^{-u}du Γ(x)=∫0∞ux−1e−udu
这使得Beta函数满足归一化条件:
∫ 0 1 B e t a ( μ ∣ a , b ) d μ = 1 \int_0^1Beta(\mu|a,b)\ d\mu=1 ∫01Beta(μ∣a,b) dμ=1
此外,似然度选用binomial分布:
B i m ( m ∣ N , μ ) = C N m μ m ( 1 − μ ) N − m Bim(m|N,\mu)=C_N^m\ \mu^m(1-\mu)^{N-m} Bim(m∣N,μ)=CNm μm(1−μ)N−m
则后验概率也将满足Beta函数的形式:
p ( μ ∣ m , l , a , b ) = Γ ( m + a + l + b ) Γ ( m + a ) Γ ( l + b ) μ m + a − 1 ( 1 − μ ) l + b − 1 p(\mu|m,l,a,b)=\dfrac{\Gamma (m+a+l+b)}{\Gamma(m+a)\Gamma(l+b)}\mu^{m+a-1}(1-\mu)^{l+b-1} p(μ∣m,l,a,b)=Γ(m+a)Γ(l+b)Γ(m+a+l+b)μm+a−1(1−μ)l+b−1
其中, l = N − m l=N-m l=N−m
现在,对于掷硬币过程中的任意时刻,都有着当前的先验概率。在此基础上再进行投掷,则可以根据投掷的情况确定似然度,与之相乘从而得到后验概率,随后再拿后验概率来替换先验概率,即可使概率分布不断地被迭代更新
注意到初始时刻,没有任何记录,因此在 B e t a ( μ ∣ a , b ) Beta(\mu|a,b) Beta(μ∣a,b)中, a = b = 1 a=b=1 a=b=1,对应的函数为常函数,其值为1
实现
Gamma函数的计算
由于Gamma函数在先验概率和后验概率表达式中频繁出现,因此计算Gamma函数的值是不可避免的。因为Gamma函数是一个广义积分,使用解析的方式求解会非常地困难。所以在这里,借助了python的scipy模块来进行数值积分计算
例如,对 Γ ( 1 ) \Gamma(1) Γ(1)的求值相关代码如下:
from scipy import integrate
import numpy as np
import math
def gammaFunc(u, x):
return u ** (x - 1) * math.exp(-u)
S, err = integrate.quad(gammaFunc, 0, np.inf, args = (1))
print("积分结果:%f" % S)
gammaFunc(u, x)中的第一个参数u是积分变量,而第二个参数x则是任意指定的变量。任意指定的变量是在调用integrate.quad()方法的时候,由args参数进行传递。例如,这里的args = (1),就代表着令gammaFunc的第二个参数x=1
运行后,得到结果:
积分结果:1.000000
现在,查看x取不同整数时的计算结果:
from scipy import integrate
import numpy as np
import math
import matplotlib.pyplot as plt
def gammaFunc(u, x):
return u ** (x - 1) * math.exp(-u)
numList = np.arange(1, 51, 1)
resultList = []
for eachNum in numList:
S, err = integrate.quad(gammaFunc, 0, np.inf, args = (eachNum))
resultList.append(S)
fig, ax = plt.subplots(1, 1, figsize = (16, 9))
ax.stem(numList, resultList)
ax.set_yscale("log")
ax.set_ylim(bottom = 0.5)
fig.savefig("gamma_int.png")
plt.show()
观察到 Γ ( x ) \Gamma(x) Γ(x)的值是随着x的增加而急剧增加的(图中纵轴使用了对数尺度)
beta函数的计算
在能够计算 Γ ( x ) \Gamma(x) Γ(x)后,就可以对不同a, b下的beta函数进行计算了:
from scipy import integrate
import numpy as np
import math
import matplotlib.pyplot as plt
def gammaFunc(u, x):
return u ** (x - 1) * math.exp(-u)
def gammaFuncV(x):
S, err = integrate.quad(gammaFunc, 0, np.inf, args = (x))
return S
def betaFunc(mu, a, b):
beta = gammaFuncV(a + b)/(gammaFuncV(a) * gammaFuncV(b)) * (mu ** (a - 1)) * ((1 - mu) ** (b - 1))
return beta
fig, ax = plt.subplots(6, 6, figsize = (36, 36))
i = 0
for a in range(1, 18, 3):
j = 0
for b in range(1, 18, 3):
mu = np.arange(0, 1, 0.01)
beta = betaFunc(mu, a, b)
ax[i][j].plot(mu, beta, 'b-')
ax[i][j].set_title("a = %d, b = %d" % (a, b))
ax[i][j].set_xlim([0, 1])
j += 1
i += 1
fig.savefig("beta.png")
由此大致作出图像,以观察不同的a、b值对beta函数曲线的影响:
注意到a和b值越相近时,beta函数曲线的极值越靠近中间
掷硬币实验的贝叶斯方法模拟
标题可能不准确,因为不知道该怎么起标题了
因为先验概率和后验概率的函数形式都是beta函数,且后验概率是在先验概率的基础上改变a、b值。所以在利用贝叶斯原理计算后验概率时,只需在先验概率的基础上更改a、b值即可,而不必专门去乘似然度
在硬币实验中,当硬币正面朝上时,似然度函数为 μ \mu μ,反之为 1 − μ 1-\mu 1−μ. 即正面朝上时,更新a=a+1,反之,更新b=b+1
由此即可进行模拟实验:
import random
from scipy import integrate
import numpy as np
import math
import matplotlib.pyplot as plt
observeTS = [0, 1, 2, 3, 5, 7, 12, 18, 30] # 在掷硬币多少次后对当前的beta函数进行观察
def gammaFunc(u, x):
return u ** (x - 1) * math.exp(-u)
def gammaFuncV(x):
S, err = integrate.quad(gammaFunc, 0, np.inf, args = (x))
return S
def betaFunc(mu, a, b):
beta = gammaFuncV(a + b)/(gammaFuncV(a) * gammaFuncV(b)) * (mu ** (a - 1)) * ((1 - mu) ** (b - 1))
return beta
# 掷硬币,返回0或1
def throwCoin():
return 0 if(random.random() * 2 < 1) else 1
# idx转二维下标
def idx2xy(idx):
x = idx // 3
y = idx % 3
return [x, y]
# 初始化beta函数的横坐标
mu = np.arange(0, 1.01, 0.01)
# 初始化画布及子图下标
fig, ax = plt.subplots(3, 3, figsize = (9, 9))
idx = 0
# 初始化beta函数的参数a、b
a = 1
b = 1
# 初始化掷硬币的次数
N = 0
# 初始化beta函数
beta = betaFunc(mu, a, b)
if(N in observeTS):
[x, y] = idx2xy(idx)
ax[x][y].plot(mu, beta, 'b-')
ax[x][y].set_title("N = %d" % N)
idx += 1
for dummyIterIdx in range(31):
observeValue = throwCoin()
if(observeValue == 1):
a += 1
beta = betaFunc(mu, a, b)
else:
b += 1
beta = betaFunc(mu, a, b)
N += 1
if(N in observeTS):
[x, y] = idx2xy(idx)
ax[x][y].plot(mu, beta, 'b-')
ax[x][y].set_title("N = %d" % N)
idx += 1
fig.savefig("coin_exp.png")
当N(掷硬币的次数)逐渐增大时,beta函数的最大值也会逐渐趋近于中心。这里计算到N=30的时候,已经快达到这个数值计算所能承受的极限了(再大一些就会有OverflowError: (34, ‘Result too large’)这个报错出现了——数值太大了),所以就没有再往后模拟了
不过按照这个情况来看,感觉这个beta函数在N趋近于无穷大的时候,会趋近于 δ \delta δ函数(即在某个点处无穷大,在其他地方全为0)