贝叶斯线性回归及最大后验估计

贝叶斯推断

贝叶斯定理：通过观察到的数据DD，把先验概率p(θ)p(θ)转化为后验概率p(θ∣D)p(θ∣D)

p(θ∣D)=p(D∣θ)p(θ)∫p(D∣θ)p(θ)dθ=p(D∣θ)p(θ)p(D)p(θ∣D)=p(D∣θ)p(θ)∫p(D∣θ)p(θ)dθ=p(D∣θ)p(θ)p(D)

显然，分母是一个归一化常数。故有p(θ∣D)∝p(D∣θ)p(θ)p(θ∣D)∝p(D∣θ)p(θ) 即后验∝似然×先验后验∝似然×先验。

贝叶斯线性回归

问题是这样的，不能够一次性接收到整个数据集，而是不断接收到小的数据集Di,i=1,2,...,nDi,i=1,2,...,n，同时由于存储的限制不能存储已经接收到的所有数据集,每次可以处理的数据仅为DiDi。这就导致不能对所有数据做线性回归，但是可以通过贝叶斯线性回归达到同样的效果。

第ii个数据集DiDi中有mm 个训练样本，构成 (X(i),y(i))(X(i),y(i))

p(y(i)∣X(i),θ)=N(y(i);X(i)θ,I)∝exp(−12(y(i)−X(i)θ)T(y(i)−X(i)θ))p(y(i)∣X(i),θ)=N(y(i);X(i)θ,I)∝exp(−12(y(i)−X(i)θ)T(y(i)−X(i)θ))

为了确定模型参数向量 θθ的后验分布

假设其先验分布

p(θ)=N(θ;μ0,Λ0)∝exp(−12(θ−μ0)TΛ−10(θ−μ0))p(θ)=N(θ;μ0,Λ0)∝exp(−12(θ−μ0)TΛ0−1(θ−μ0))

其中μ0μ0,Λ0Λ0分别是先验分布的均值向量和协方差矩阵。通过贝叶斯回归得到的目标为θθ的期望。

模型参数的后验分布：

p(θ∣X(i),y(i)) ∝p(y(i)∣X(i),θ)p(θ) ∝exp(−12(y(i)−X(i)θ)T(y(i)−X(i)θ))exp(−12(θ−μ0)TΛ−10(θ−μ0)) ∝exp(−12(−2y(i)TX(i)θ+θTX(i)TX(i)θ+θTΛ−10θ−2μT0Λ−10θ))p(θ∣X(i),y(i)) ∝p(y(i)∣X(i),θ)p(θ) ∝exp(−12(y(i)−X(i)θ)T(y(i)−X(i)θ))exp(−12(θ−μ0)TΛ0−1(θ−μ0)) ∝exp(−12(−2y(i)TX(i)θ+θTX(i)TX(i)θ+θTΛ0−1θ−2μ0TΛ0−1θ))

Λi=(X(i)TX(i)+Λ−10)−1,μi=Λi(X(i)Ty(i)+Λ−10μ0)Λi=(X(i)TX(i)+Λ0−1)−1,μi=Λi(X(i)Ty(i)+Λ0−1μ0)

p(θ∣X(i),y(i))∝exp(−12(θ−μi)TΛ−1i(θ−μi))p(θ∣X(i),y(i))∝exp(−12(θ−μi)TΛi−1(θ−μi))

缺点：

1 参数先验分布的不同假设形式，可能会带来计算上的不便。

2 参数先验分布的假设有偏，对于小数据会有较大的影响。

解决方法：

1 参数的先验分布假设为数据分布假设的共轭先验

共轭先验：对于一个给定的概率分布p(x∣∣w)p(x∣∣w),能够寻找一个先验 p(w)p(w) 能够与似然函数共轭，从而后验分布的函数形式与先验分布相同。

Bern(x∣μ)=μx(1−μ)1−xBeta(μ∣a,b)=Γ(a+b)Γ(a)Γ(b)μa−1(1−μ)b−1Bern(x∣μ)=μx(1−μ)1−xBeta(μ∣a,b)=Γ(a+b)Γ(a)Γ(b)μa−1(1−μ)b−1

N(x∣m,Λ)=12π|Λ|exp(−12(x−m)TΛ−1(x−m))N(m∣μ,Λ′)=12π|Λ′|exp(−12(m−μ)T(Λ′)−1(m−μ))W(Λ∣W,v)=B|Λ|v−D−12exp(−12Tr(W−1Λ))N(x∣m,Λ)=12π|Λ|exp(−12(x−m)TΛ−1(x−m))N(m∣μ,Λ′)=12π|Λ′|exp(−12(m−μ)T(Λ′)−1(m−μ))W(Λ∣W,v)=B|Λ|v−D−12exp(−12Tr(W−1Λ))

2 合理初始化，迭代求解

对于接收到的第1个数据集有：

Λ1=(X(1)TX(1)+Λ−10)−1,μ1=Λ1(X(1)Ty(1)+Λ−10μ0)Λ1=(X(1)TX(1)+Λ0−1)−1,μ1=Λ1(X(1)Ty(1)+Λ0−1μ0)

p(θ∣X(1),y(1))∝exp(−12(θ−μ1)TΛ−11(θ−μ1))p(θ∣X(1),y(1))∝exp(−12(θ−μ1)TΛ1−1(θ−μ1))

这里根据极大似然估计得到的解θ=(X(1)TX(1))−1X(1)Ty(1)θ=(X(1)TX(1))−1X(1)Ty(1)，所以假设Λ−10=OΛ0−1=O, 此时极大似然的解和贝叶斯回归的参数期望一致。

对于接收到的第ii个数据集Di(i>1)Di(i>1)，将第i−1i−1 个数据集计算得到的参数后验作为先验，不断迭代。

Λi=(X(i)TX(i)+Λ−1i−1)−1,μi=Λi(X(i)Ty(i)+Λ−1i−1μi−1)Λi=(X(i)TX(i)+Λi−1−1)−1,μi=Λi(X(i)Ty(i)+Λi−1−1μi−1)

p(θ∣X(i),y(i))∝exp(−12(θ−μi)TΛ−1i(θ−μi))p(θ∣X(i),y(i))∝exp(−12(θ−μi)TΛi−1(θ−μi))

具体算法

输入：D1,D2,D3,...,DnD1,D2,D3,...,Dn 其中 Di=(X(i),y(i))Di=(X(i),y(i))

输出：μnμn

初始化

Λ1=(X(1)TX(1))−1μ1=Λ1(X(1)Ty(1))i+=1Λ1=(X(1)TX(1))−1μ1=Λ1(X(1)Ty(1))i+=1

while i<=ni<=n

Λi=(X(i)TX(i)+Λ−1i−1)−1μi=Λi(X(i)Ty(i)+Λ−1i−1μi−1)i+=1Λi=(X(i)TX(i)+Λi−1−1)−1μi=Λi(X(i)Ty(i)+Λi−1−1μi−1)i+=1

代码：

def BayesLR(path):
la=10
mu=np.mat(np.zeros(3)).T
gama=np.mat(np.eye(3)*la)
for i in range(n):
fileName = path + "%d.csv" % i
x0,y0 = loadDataFromFile(fileName)#从文件中加载数据
X, y = data2Mat(x0,y0)#将数据转换成np.mat的格式
mu0 = mu
gama0 = gama
if i==1:
gama = (X.T*X).I
mu = gama*(X.T*y)

else:
gama = (X.T*X+gama0.I).I
mu = gama*(X.T*y+gama0.I*mu0)

return np.array(mu)

最大后验估计（MAP）

MLE求的是找出一组能够使似然函数最大的参数，即

。 现在问题稍微复杂一点点，假如这个参数μμ有一个先验概率呢？比如说，在上面抛硬币的例子，假如我们的经验告诉我们，硬币一般都是匀称的，也就是μ=0.5μ=0.5的可能性最大，μ=0.2μ=0.2的可能性比较小，那么参数该怎么估计呢？这就是MAP要考虑的问题。 MAP优化的是一个后验概率，即给定了观测值后使概率最大：



把上式根据贝叶斯公式展开：