变分推断学习笔记(3)——三硬币问题的变分推断解法

变分推断学习笔记系列:

  1. 变分推断学习笔记(1)——概念介绍
  2. 变分推断学习笔记(2)——一维高斯模型的例子
  3. 变分推断学习笔记(3)——三硬币问题的变分推断解法

其实三硬币的例子不写,前面的介绍也够了,写这个纯粹是吃撑了。这次我们采取更加普遍的假设,将原先假设的3枚硬币拓展开来。现在假设有\(K+1\)个骰子,第一个骰子有\(K\)个面,其余的骰子有\(T\)个面。进行如下实验:先掷第一个骰子,根据投出的结果\(Z_k\),选择第\(Z_k\)个骰子再投,观测到投出的\(N\)个结果,每个结果\(w_n\)可能是 \[ 1,3,7,8,3,2,6,9,... \]

可以看到现在第1个骰子投出的标签服从多项分布: \[Z_k \sim Multinomial(\pi)\] 然后剩余骰子投出的面也服从多项分布 \[W_{Z_{kt}} \sim Multinomial(\theta_{Z_k})\] 我们假设,随机变量\(\pi\)\(\theta\)的先验分布为狄利克雷分布,超参分别为\(\alpha\)\(\beta\)

Read More

概率图模型简介

介绍

定义

概率图模型(Graphical Model)是概率论和图论之间的桥梁,概率图模型的基本想法来自模块的概念,即一个复杂系统是由简单的部分所联合而成的。概率论部分告诉我们哪些部分是耦合在一起的,然后提供推断模型的方法,而图论部分给我们一个非常直观的认识,把拥有相互关系的变量看做数据结构,从而导出一般化的方法来解决这个问题。很多经典的多元概率系统,比如混合模型,因子分析,隐含马尔科夫模型,Kalman filter 和Ising model等,从概率图模型的角度来看,都可以看做普遍隐含形成过程的一种实例表现。这意味着,一旦在某个系统上有什么特别的方法发现,就很容易推广到一系列的模型中去。除此之外,概率图模型还非常自然地提供了设计新系统的方法。

在概率图模型中,点代表随机变量,点与点之间边的存在与否代表了点与点之间是存在条件依赖。点与边的组合描绘了联合概率分布的特征结构。假设有\(N\)个二元随机变量,在没有任何信息帮助的情况下,联合分布\(P(X_1,\ldots,X_N)\),需要\(O(2^N)\)个参数。而通过概率图描绘点与点之间的条件关系之后,表示联合分布,所需要的参数会减少很多,这对后面的模型的推断和学习是有帮助的。

概率图模型主要分为两种,无向图(也叫Markov random fields)和有向图(也叫Bayesian networks),前者多用于物理和图像领域,后者多用于AI和机器学习,具体的基本就不多介绍了。下面给出一个简单贝叶斯网络的例子。

QQ截图20140224133640

这里Cloudy指是否多云,Sprinkler指洒水器是否打开,Rain指是否有雨, WetGrass指草地是否湿了。

Read More

变分推断学习笔记(2)——一维高斯模型的例子

变分推断学习笔记系列:

  1. 变分推断学习笔记(1)——概念介绍
  2. 变分推断学习笔记(2)——一维高斯模型的例子
  3. 变分推断学习笔记(3)——三硬币问题的变分推断解法

举一个一元高斯模型的例子。假设我们有数据\(\mathbf{X}=\{x_1,\ldots,x_M\}\),要推断平均值\(\mu\)和精度\(\tau(1/\sigma)\)的后验概率分布。 写出似然 \[\begin{equation} p(\mathbf{X}|\mu,\tau)=(\frac{\tau}{2\pi})^{N/2}\exp\{-\frac{\tau}{2}\sum^N_{n=1}(x_n-\mu)^2\} \end{equation}\] 其中\(\mu,\tau\)各自服从先验分布 \[\begin{equation}p(\mu|\tau)=N(\mu|\mu,(\lambda_0\tau)^{-1})\end{equation}\] \[\begin{equation}p(\tau)=Gam(\tau|a_0,b_0)\end{equation}\] 其中Gam为Gamma分布(见备注1)。

通用的估计方法

好,我们现在假设\(q\)之间的分布都独立。 \[\begin{equation}q(\mu,\tau)=q_u(\mu)q_r(\tau)\end{equation}\]

Read More

变分推断学习笔记(1)——概念介绍

变分推断学习笔记系列:

  1. 变分推断学习笔记(1)——概念介绍
  2. 变分推断学习笔记(2)——一维高斯模型的例子
  3. 变分推断学习笔记(3)——三硬币问题的变分推断解法

问题描述

变分推断是一类用于贝叶斯估计和机器学习领域中近似计算复杂(intractable)积分的技术,它广泛应用于各种复杂模型的推断。本文是学习PRML第10章的一篇笔记,错误或不足的地方敬请指出。

先给出问题描述。记得在上一篇EM的文章中,我们有一个观察变量\(\mathbf{X}=\{x^{\{1\}},\ldots,x^{\{m\}}\}\)和隐藏变量\(\mathbf{Z}=\{z^{\{1\}},\ldots,z^{\{m\}}\}\), 整个模型\(p(\mathbf{X},\mathbf{Z})\)是个关于变量\(\mathbf{X},\mathbf{Z}\)的联合分布,我们的目标是得到后验分布\(P(\mathbf{Z}|\mathbf{X})\)的一个近似分布。

在之前介绍过Gibbs Sampling这一类Monte Carlo算法,它们的做法就是通过抽取大量的样本估计真实的后验分布。而变分推断不同,与此不同的是,变分推断限制近似分布的类型,从而得到一种局部最优,但具有确定解的近似后验分布。

之前在EM算法的介绍中我们有似然的式子如下: \[\begin{equation}\ln p(\mathbf{X})=L(q)+KL(q||p)\end{equation}\] 其中
\[\begin{equation}L(q)=\int q(\mathbf{Z})\ln{\frac{p(\mathbf{X},\mathbf{Z})}{q(\mathbf{Z})}}d\mathbf{Z}\end{equation}\]
\[\begin{equation}KL(q||p)=-\int q(\mathbf{Z}) \ln{\frac{p(\mathbf{Z}|\mathbf{X})}{q(\mathbf{Z})}}d\mathbf{Z}\end{equation}\]

这里公式中不再出现参数\(\theta\),因为参数不再是固定的值,而变成了随机变量,所以也是隐藏变量,包括在\(\mathbf{Z}\)之内了。

Read More

三硬币问题-一个EM算法和Gibbs Sampling的例子

讲一个EM算法和Gibbs 抽样的小例子,用于加深理解(变分推断版本请见变分推断学习笔记(3)——三硬币问题的变分推断解法)。

题目(引用自参考1):假设有3枚硬币,分别记做A,B,C。这些硬币正面出现的概率分别是\(\pi\),\(p\)\(q\)。进行如下掷硬币实验:先掷硬币A,根据其结果选出硬币B或C,正面选B,反面选硬币C;然后投掷选重中的硬币,出现正面记作1,反面记作0;独立地重复\(n\)次(n=10),结果为 \[1111110000\] 我们只能观察投掷硬币的结果,而不知其过程,估计这三个参数\(\pi\),\(p\)\(q\)

EM算法

可以看到投掷硬币时到底选择了B或者C是未知的。我们设隐藏变量Z 来指示来自于哪个硬币,\(Z=\{z_1,z_2,\ldots,z_n \}\),令\(\theta=\{\pi,p,q\}\),观察数据\(X=\{x_1,x_2,\ldots,x_n \}\)

写出生成一个硬币时的概率: \[\begin{split}P(x|\theta) & =\sum_z P(x,z|\theta)=\sum_z P(z|\pi)P(x|z,\theta) \\& =\pi p^x (1-p)^{1-x}+(1-\pi)q^x(1-q)^{1-x} \\\end{split}\] 有了一个硬币的概率,我们就可以写出所有观察数据的log似然函数: \[L(\theta|X)=\log P(X|\theta)=\sum^n_{j=1}\log[\pi p^{x_j} (1-p)^{1-{x_j}}+(1-\pi)q^{x_j}(1-q)^{1-{x_j}}]\] 然后求极大似然 \[\hat{\theta}=\arg \max L(\theta|X)\] 其中\(L(\theta|X)=\log P(X|\theta)=\log \sum_Z P(X,Z|\theta)\)。因为log里面带着加和所以这个极大似然是求不出解析解的。

Read More

《Gibbs Sampling for the UniniTiated》阅读笔记(下)---连续型参数求积分的思考

《Gibbs Sampling for the UniniTiated》阅读笔记结构:

  1.  参数估计方法及Gibbs Sampling简介
  2. 一个朴素贝叶斯文档模型例子
  3. 连续型参数求积分的思考

这篇是下篇,讨论中篇联合分布中对参数求积分来简化的问题。

之前存在的一个问题就是为啥我们可以对连续参数\(\pi\)求积分消去它,而不能对词分布\(\theta_0\)\(\theta_1\)求积分。这个主意看上去很美,但是实际做的时候,你会碰到一大把无法约掉的伽马函数。让我们看看具体的过程。

Read More

《Gibbs Sampling for the UniniTiated》阅读笔记(中)---一个朴素贝叶斯文档模型例子

《Gibbs Sampling for the UniniTiated》阅读笔记结构:

  1.  参数估计方法及Gibbs Sampling简介
  2. 一个朴素贝叶斯文档模型例子
  3. 连续型参数求积分的思考

这篇是中篇,介绍一个非常简单的朴素贝叶斯文档模型生成的例子,用来说明Gibbs Sampler具体是如何构造的。

文档生成的建模过程

首先我们有一批文档,文档里面有很多单词,这些单词都是无顺序可交换的(词袋模型),这些文档分成两类,类标签为0或者1。给予一篇未标记的文档\(W_j\),我们要做的工作就是预测文档的类标签是\(L_j=0\)还是\(L_j=1\)。为了方便起见,我们定了类标签所表示的类\(\mathbb{C}_0={W_j|L_j=0}\)\(\mathbb{C}_1={W_j|L_j=1}\)。一般来说预测这种事都是选择最有可能发生的,即找到\(W_j\)的后验概率\(P(L_j|W_j)\)最大的标签\(L_j\)。使用贝叶斯公式 \[\begin{equation} \begin{split} L_j=\arg \max \limits_{L}P(L|W_j)& =\arg \max \limits_{L}\frac{P(W_j|L)P(L)}{P(W_j)}\\& =\arg \max \limits_{L} P(W_j|L)P(L) \\\end{split} \end{equation}\] 因为分母\(P(W_j)\)\(L\)无关所以删去了。 通过贝叶斯公式的转换,我们可以想象这些文档的生成过程。首先,我们选择文档的类标签\(L_j\);假设这个过程是通过投硬币完成的(正面概率为\(\pi=P(L_j=1)\) ),正式地来说,就是服从贝努利分布 \[\begin{equation}L_j \sim Bernoulli(\pi)\end{equation}\] 然后,对于文档上\(R_j\)个“词位”中的每一个,我们根据一个概率分布\(\theta\),随机独立地抽样一个词\(w_i\)。因为每个类生成词的\(\theta\)分布都不同,所以应该有\(\theta_1\)\(\theta_2\),具体地生成词的时候,我们根据文档的标签\(L_j\)来决定由哪个类来生成 \[\begin{equation} W_j \sim Multinomial(R_j,\theta_{L_j}) \end{equation}\]

Read More

《Gibbs Sampling for the UniniTiated》阅读笔记(上)---参数估计方法及Gibbs Sampling简介

前一阵子折腾的事儿太多,写了点东西都没有传上来,是我偷懒了- -,下不为例。

这篇文章基本上是来自于《Gibbs Sampling for the UniniTiated》,说是笔记其实和翻译也差不多了。

整个结构分为上中下三部分:

  1.  参数估计方法及Gibbs Sampling简介
  2. 一个朴素贝叶斯文档模型例子
  3. 连续型参数求积分的思考

这篇是上部分,介绍基础参数估计和Gibbs Sampling概念。

为什么求积分—参数估计方法

很多概率模型的算法并不需要使用积分,只要对概率求和就行了(比如隐马尔科夫链的Baum-Welch算法),那么什么时候用到求积分呢?—— 当为了获得概率密度估计的时候,比如说根据一句话前面部分的文本估计下一个词的概率,根据email的内容估计它是否是垃圾邮件的概率等等。为了估计概率密度,一般有MLE(最大似然估计),MAP(最大后验估计),bayesian estimation(贝叶斯估计)三种方法。

最大似然估计

这里举一个例子来讲最大似然估计。假设我们有一个硬币,它扔出正面的概率\(\pi\)不确定,我们扔了10次,结果为HHHHTTTTTT(H为正面,T为反面)。利用最大似然估计的话,很容易得到下一次为正面的概率为0.4,因为它估计的是使观察数据产生的概率最大的参数。 first

\(\chi=\{HHHHTTTTTT\}\)代表观察到的数据,\(y\)为下一次抛硬币可能的结果,估计公式如下: \[\begin{equation}\begin{split}\tilde{\pi}_{MLE} &=\arg \max \limits_{\pi}P(\chi|\pi) \\P(y|\chi) & \approx \int_{\pi} p(y|\tilde{\pi}_{MLE})P(\pi|\chi) d\pi = p(y|\tilde{\pi}_{MLE})\end{split}\end{equation}\]

Read More

高斯混合模型的matlab实现(转)

高斯混合函数实现部分是基本上是转载的的pluskid大神文章里的里的代码,加了一点注释,并根据他给的方法二解决 covariance 矩阵 singular 的问题。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
function varargout = gmm(X, K_or_centroids)
% ============================================================
%转载自http://blog.pluskid.org/?p=39
% Expectation-Maximization iteration implementation of
% Gaussian Mixture Model.
%
% PX = GMM(X, K_OR_CENTROIDS)
% [PX MODEL] = GMM(X, K_OR_CENTROIDS)
%
% - X: N-by-D data matrix.%需要注意的是这里的X包括了全部
% - K_OR_CENTROIDS: either K indicating the number of
% components or a K-by-D matrix indicating the
% choosing of the initial K centroids.
%
% - PX: N-by-K matrix indicating the probability of each
% component generating each point.
% - MODEL: a structure containing the parameters for a GMM:
% MODEL.Miu: a K-by-D matrix.
% MODEL.Sigma: a D-by-D-by-K matrix.
% MODEL.Pi: a 1-by-K vector.
% ============================================================
threshold = 1e-15;
[N, D] = size(X);
if isscalar(K_or_centroids)
K = K_or_centroids;
% randomly pick centroids
rndp = randperm(N);
centroids = X(rndp(1:K),:);
else
K = size(K_or_centroids, 1);
centroids = K_or_centroids;
end
% initial values
[pMiu pPi pSigma] = init_params();
Lprev = -inf;
while true
Px = calc_prob();%计算N(x|mu,sigma)
% new value for pGamma
pGamma = Px .* repmat(pPi, N, 1);%估计 gamma 是个N*K的矩阵
pGamma = pGamma ./ repmat(sum(pGamma, 2), 1, K);%对矩阵的理解真是出神入化,
% new value for parameters of each Component
Nk = sum(pGamma, 1);%N_K
pMiu = diag(1./Nk) * pGamma' * X; %数字 *( K-by-N * N-by-D)加个括号有助理解
pPi = Nk/N;
for kk = 1:K
Xshift = X-repmat(pMiu(kk, : ), N, 1);%x-u
pSigma(:, :, kk) = (Xshift' * ...
(diag(pGamma(:, kk)) * Xshift)) / Nk(kk);%更新sigma
end
% check for convergence
L = sum(log(Px*pPi'));
if L-Lprev < threshold
break;
end
Lprev = L;
end
if nargout == 1
varargout = {Px};
else
model = [];
model.Miu = pMiu;
model.Sigma = pSigma;
model.Pi = pPi;
varargout = {pGamma, model};%注意!!!!!这里和大神代码不同,他返回的是px,而我是 pGamma
end
function [pMiu pPi pSigma] = init_params()%初始化参数
pMiu = centroids;% K-by-D matrix
pPi = zeros(1, K);%1-by-K matrix
pSigma = zeros(D, D, K);%
% hard assign x to each centroids
distmat = repmat(sum(X.*X, 2), 1, K) + ... % X is a N-by-D data matrix.
repmat(sum(pMiu.*pMiu, 2)', N, 1) - ...% X->K列 U->N行 XU^T is N-by-K
2*X*pMiu';%计算每个点到K个中心的距离
[~, labels] = min(distmat, [], 2);%找到离X最近的pMiu,[C,I] labels代表这个最小值是从那列选出来的
for k=1:K
Xk = X(labels == k, : );% Xk是所有被归到K类的X向量构成的矩阵
pPi(k) = size(Xk, 1)/N;% 数一数几个归到K类的
pSigma(:, :, k) = cov(Xk); %计算协方差矩阵,D-by-D matrix,最小方差无偏估计
end
end
function Px = calc_prob()
Px = zeros(N, K);
for k = 1:K
Xshift = X-repmat(pMiu(k, : ), N, 1);%x-u
lemda=1e-5;
conv=pSigma(:, :, k)+lemda*diag(diag(ones(D)));%这里处理singular问题,为协方差矩阵加上一个很小lemda*I
inv_pSigma = inv(conv);%协方差的逆
tmp = sum((Xshift*inv_pSigma) .* Xshift, 2);%(X-U_k)sigma.*(X-U_k),tmp是个N*1的向量
coef = (2*pi)^(-D/2) * sqrt(det(inv_pSigma));%前面的参数
Px(:, k) = coef * exp(-0.5*tmp);%把数据点 x 带入到 Gaussian model 里得到的值
end
end
end
%repmat 通过拓展向量到矩阵
%inv 求逆
%min 求矩阵最小值,可以返回标签
%X(labels == k, : ) 对行做筛选
% size(Xk, 1) 求矩阵的长或宽
%scatter 对二维向量绘图

注意:

pluskid大神这里最后返回的是px,我觉得非常奇怪,因为PRML里对点做hard assignment时是根据后验概率来判别的。于是我在大神博客上问了一下,他的解释是最大似然和最大后验的区别,前者是挑x被各个模型产生的概率最大的那个,而后者加上了先验知识,各有道理。一句话就茅塞顿开,真大神也~

Read More

高斯混合模型参数估计详细推导过程

已知多元高斯分布的公式: \[N(x|\mu,\Sigma)=\frac{1}{(2\pi)^{D/2}|\Sigma|^{1/2}}\exp(-\frac{1}{2}(x-\mu)^T\Sigma^{-1}(x-\mu))\] 其中\(D\)为维度,\(x\)\(\mu\)均为\(D\)维向量,协方差\(\Sigma\)为D维矩阵。我们求得后验概率: \[w^{(i)}_j=Q_i(Z^{i}=j)=P(z^{(i)}=j|x^{(i)};\phi,\mu,\Sigma)\] 在E步,\(w^{(i)}_j\)是一个固定值,然后我们用它来估计似然函数\(L(X,Z;\theta)\)(这里\(\theta=(\phi,\mu,\Sigma)\))在分布\(Z\sim P(Z|X;\theta)\)上的期望\(E_{Z|X,\theta_t}[L(X,Z;\theta)]\)(式子1): \[\begin{split} & \sum^m_{i=1}\sum_{z^{(i)}} Q_i(z^{(i)})\log{\frac{p(x^{(i)},z^{(i)};\phi,\mu,\Sigma)}{Q_i(z^{(i)})}} \\& =\sum^m_{i=1}\sum^k_{j=1} Q_i(z^{(i)}=j)\log{\frac{p(x^{(i)}|z^{(i)}=j;\mu,\Sigma)p(z^{(i)}=j;\phi)}{Q_i(z^{(i)})}} \\& =\sum^m_{i=1}\sum^k_{j=1} w^{(i)}_j\log{\frac{\frac{1}{(2\pi)^{D/2}|\Sigma|^{1/2}}\exp(-\frac{1}{2}(x^{(i)}-\mu_j)^T\Sigma_j^{-1}(x^{(i)}-\mu_j))\cdot\phi_j}{ w^{(i)}_j}} \\\end{split}\] 由于分母\(w^{(i)}_j\)在取对数之后是常数,与参数无关,求导时自然会变成0,所以我们写公式的时候为了简便舍去分母。

Read More