讲一个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》阅读笔记结构：

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

   ---

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

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

Read More

《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》，说是笔记其实和翻译也差不多了。

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

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

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

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

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

最大似然估计

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

令\(\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

SyntaxHighlighter是个WordPress（也支持别的博客）上非常好用的语法高亮插件，但是在该插件的设置里没给出默认支持的高亮语言列表（让你自己去官网上找）。查了一下插件的文件，发现除去官网上所列的以外，还支持latex和matlab：

Read More

高斯混合函数实现部分是基本上是转载的的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

1. 社区发现简介

社区，从直观上来看，是指网络中的一些密集群体，每个社区内部的结点间的联系相对紧密，但是各个社区之间的连接相对来说却比较稀疏（图1，当然社区的定义不止有这一种）。这样的社区现象被研究已经很多年了，最早期的记录甚至来自于80年前。

比较经典的社区研究案例包括对空手道俱乐部(karate club),科学家合作网络(Collaboration network) 和斑马群体(zebras) 的社交行为研究等（见图2），其中著名的空手道俱乐部社区已经成为通常检验社区发现算法效果的标准(benchmark)之一。

Read More

2013年10月10日更新。

LDA的概率图如下图1所示：

参数的意思如图2所示：

根据模型，文章m的第n个词t是这样生成的：先从文章m的doc-topic分布中生成一个topic编号\(z_{m,n}\)，在根据编号第\(z_{m,n}\)个的topic-word分布中生成这个词，总够有\(K\)个topic，所以总的概率为： \[ p(w_{m,n}=t|\vec{\theta}_m,\underline{\Phi})=\sum^K_{k=1}p(w_{m,n}=t|\vec{\phi}_k)p(z_{m,n}=k|\vec{\theta}_m)\] 如果我们写出这篇文章的complete-data的联合分布（意思就是所以变量都已知的情况下），那么式子就是这样的：

通过对\(\vec{\vartheta_m}\)（doc-topic分布）和\(\underline{\Phi}\)（topic-word分布）积分以及\(z_{m,n}\)求和，我们可以求得\(\vec{w_m}\)的边缘分布：

(实际上这个边缘分布是求不出来的，因为\(z_{m,n}\)是隐藏变量，从而导致\(\underline{\vartheta}\)与\(\underline{\Phi}\)存在耦合现象，无法积分得到。要注意联合分布和边缘分布对Z乘积与加和的区别）

因为一个语料库有很多篇文章，而且文章之间都是相互独立的，所以整个语料库的似然为 \[ p(\mathcal{W}|\vec{\alpha},\vec{\beta})=\prod^{M}_{m=1}p(\vec{w_m}|\vec{\alpha},\vec{\beta})\]

虽然LDA（latent Dirichlet allocation)是个相对简单的模型，对它直接推断一般也是不可行的，所以我们要采用近似推断的方法，比如Gibbs sampling。

Read More

三门问题（Monty Hall problem），是一个源自博弈论的数学游戏问题，大致出自美国的电视游戏节目Let’s Make a Deal。问题的名字来自该节目的主持人蒙提·霍尔（Monty Hall）。问题非常的有意思^_^，给出叙述如下：

在三扇门中的某扇门以后有一个奖品，选中这扇门就能拿到门后的奖品。你选定了一扇门，具体地说，假设你选择了1号门。这时候主持人蒙提·霍尔会打开剩下两扇门的其中一扇，你看到门后没有奖品。这时候他给你一个机会选择要不要换另外一扇没有打开的门。你是选择换还是不换呢？

答：因为我之前就选了，换或者不换机会都是均等的，所以换不换无关紧要╮(╯▽╰)╭。 ↓ ↓ ↓ ↓ ↓ ↓ ↓ ↓ ↓ ↓ 真的是这样么？仔细分析一下，游戏过程中你做了2个操作： 第一，你选择了一扇门。第二，你选择了换门或者不换。 定义事件\(A\)为你第一次就选中奖品。 定义事件\(B\)为你换门选中奖品。 那么\(A^c\)为集合A的余集，即第一次没有选中奖品。同理，\(B^c\)为换门没有选中奖品。整个游戏过程中，\(A\)或\(A^c\)先发生，\(B\)或\(B^c\)再发生。 显而易见的是 \[P(A)=1/3,P(A^c)=2/3\] 要注意的是，在第一次操作之后，还有一件事——主持人打开了一扇没有奖品的门。值得注意的是，这里主持人的动作是跟你第一次选择有关系的：

1. 如果你一开始就选中了奖品，即事件\(A\)发生了，那么他就在剩下的两扇没有奖品的门之间任选一门打开。接下来，如果你选择换门，那么抽中的概率为0，不换抽中的概率为1。 即 \[P(B|A)=0,P(B^c|A)=1\] 因为事件\(B\)或\(B^c\)是在事件A发生之后发生的，所以这里的概率是基于事件\(A\)的条件概率。

2. 如果你开始没有选中奖品，即事件\(A^c\)发生了，那么他只能打开另一扇没有奖品的门。这时候，如果你选择换门，那么抽中的概率为1，不换抽中的概率为0。 即 \[P(B|A^c)=1,P(B^c|A^c)=0\] 这里的概率是基于事件\(A^c\)的条件概率。

由上我们可以发现一点，事件\(A\)会影响事件\(B\)的概率，即事件\(A\)和事件\(B\)并不是相互独立的。\(A^c\)和\(B^c\)也是同理。 于是，利用全概率公式，我们可以求得换门选中奖品的概率为 \[P(B)=P(B|A)P(A)+P(B|A^c)P(A^c)=0*1/3+1*2/3=2/3\] \[P(B^c)=P(B^c|A)P(A)+P(B^c|A^c)P(A^c)=1*1/3+0*2/3=1/3\] 所以事实上你换门能得到奖品的概率为2/3，是不换门的2倍（懂点数学真好啊）。是不是和直觉不太一样？个人认为，直观上觉得\(P(B)=1/2\)的原因是忽略了第一次选择时，通过主持人的动作改变了换门的事件概率这一客观过程。 2013年7月3日更新
从信息论的角度来说，B事件的熵为H(B)，在A事件发生之后B事件的条件熵为H(B|A)，可以证明 \[H(B) \geq H(B|A) \] 也就是说，在给予了A的信息之后，B的不确定性下降了。