7月 17

Error BackPropagation 实践 & PRML 5.3

本文主要记录BP矩阵实现。鉴于GPU服务器还没买来,现在暂时用C++编码,反正矩阵上的实现都是一样的。若公式显示有问题需要多刷新几次。

Back Propagation虽然很简单,但是每次要实现就又不太清楚了,这次可以记下来,便于温习。主要是矩阵形式有点费脑子。正好要弄个Autoencoder,就以它为例。网络如下。弄好了可以用来搞Semi-Supervised Recursive Autoencoders for Predicting Sentiment Distributions (Socher)。

正好在这个例子中,推导跟PRML5.3是一模一样的,直接抄PRML。把上面那个图一层Autoencoder抠出来,配上PRML Figure 5.7的符号,构成下图。

主要公式:$\frac{\partial E_n}{\partial w_{ji}}=\frac{\partial E_n}{\partial a_j}\frac{\partial a_j}{\partial w_{ji}}$,其中右侧$\frac{\partial a_j}{\partial w_{ji}}=z_i=x_i$

​左侧:$\delta_j=\frac{\partial E_n}{\partial a_j}=\sum\limits_k \frac{\partial E_n}{\partial a_k}\frac{\partial a_k}{\partial a_{j}}=f'(a_j)\sum\limits_k w_{kj}\delta_k$

现利用上面公式计算 PRML 5.3.2 中的例子,准备工作如下:

语言:C++,VS2012,矩阵库:MTL 4,依赖:boost,其中initializer_list出错,要不安装Microsoft Visual C++ Compiler Nov 2012 CTP,要不直接下载本文中的那个文件到include即可。

符号表示:
​WS = WORD_EMBEDDING_SIZE​
$w^{(1)}$ = mtr_encode
$w^{(2)}$ = mtr_decode
$\delta_k$ = vec_delta_out
$\delta_j$ = vec_delta_hidden
$x_i$ = vec_input
$z_i$ = vec_rep
$y_k$ = vec_output
 

图例如下,所有向量均为列向量:

首先应正向传播

mtl_col_vector vec_input(2*WORD_EMBEDDING_SIZE, 0.0);
mtl_col_vector vec_output(2*WORD_EMBEDDING_SIZE, 0.0);
mtl_col_vector vec_a(WORD_EMBEDDING_SIZE, 0.0);

//the size of row is the size of result vector
//the size of column is the size of input vector
mtl_row_matrix mtr_encode(WORD_EMBEDDING_SIZE, 2*WORD_EMBEDDING_SIZE, 0.0);
mtl_row_matrix mtr_decode(2*WORD_EMBEDDING_SIZE, WORD_EMBEDDING_SIZE, 0.0);
vec_a = mtr_encode * vec_input;

for(int j=0; j<WORD_EMBEDDING_SIZE; j++)
	    vec_rep[j] = (exp(vec_a[j])-exp(-1.0*vec_a[j]))/(exp(vec_a[j])+exp(-1.0*vec_a[j]));

vec_output = mtr_decode * vec_rep;

开始Back Propagation,第一步:

$\delta_k=y_k-t_k$  (PRML 5.65)

vec_delta_out = vec_output - vec_input

第二步,最关键的一步,计算$\delta_j$,

$\delta_j=(1-z_j^2)\sum\limits_{k=1}^{K}w_{kj}\delta_k$  (PRML 5.66)

//计算左半部分 1-z_j^2
mtl_col_vector vec_rep_square(WORD_EMBEDDING_SIZE, 0.0);
for(int j=0; j<word_embedding_size; j++)
    vec_rep_square[j] = 1.0 - vec_rep[j]*vec_rep[j];

//计算右半部分,如下图
mtl_col_vector vec_sigma_kj;
vec_sigma_kj = mtl::matrix::trans(mtr_decode)*vec_delta_out;

//不想查接口了,我自己直接替他把两个向量乘一下
mtl_col_vector vec_delta_hidden;
for(int j=0; j<WORD_EMBEDDING_SIZE; j++)
    vec_delta_hidden[j] = vec_rep_square[j] * vec_sigma_kj[j];

第三步,调整权重,

$\frac{\partial E_n}{\partial w_{ji}^{(1)}} = \delta_j x_i$  $\frac{\partial E_n}{\partial w_{kj}^{(2)}} = \delta_k z_j$(PRML 5.67)

mtl_row_matrix mtr_adjust_encode;
mtl_row_matrix mtr_adjust_decode;

mtr_adjust_encode = (LEARN_RATE*vec_delta_hidden) * mtl::vector::trans(vec_input);
mtr_adjust_decode = (LEARN_RATE*vec_delta_out) * mtl::vector::trans(vec_rep);

mtr_encode -= mtr_adjust_encode;
mtr_decode -= mtr_adjust_decode;

 

最后应该Check Gradient,见Standford UFLDL课程

// gradient checking
mtl_row_matrix w1, w2;
w1 = mtr_adjust_encode/LEARN_RATE;
w2 = mtr_adjust_decode/LEARN_RATE;
double epsilon = 1e-5;

for(int i=0; i<WORD_EMBEDDING_SIZE; i++)
{
	for(int j=0; j<WORD_EMBEDDING_SIZE*2; j++)
    {
		mtl_row_matrix w_offset = mtr_encode;
		w_offset[i][j] -= epsilon;

		//calc forword error
		mtl_col_vector vec_tmp_a(WORD_EMBEDDING_SIZE, 0.0);
		mtl_col_vector vec_tmp_rep(WORD_EMBEDDING_SIZE, 0.0);
		vec_tmp_a = w_offset * vec_input;

		for(int k=0; k<WORD_EMBEDDING_SIZE; k++)
			vec_tmp_rep[k] = (exp(vec_tmp_a[k])-exp(-1.0*vec_tmp_a[k]))/(exp(vec_tmp_a[k])+exp(-1.0*vec_tmp_a[k]));

		mtl_col_vector vec_tmp_out(2*WORD_EMBEDDING_SIZE, 0.0);
		vec_tmp_out = mtr_decode * vec_tmp_rep;

		double tmp_error = 0.0;
		for(int k=0; k<2*WORD_EMBEDDING_SIZE; k++) {
			tmp_error += (vec_input[k]-vec_tmp_out[k])*(vec_input[k]-vec_tmp_out[k])/2.0;
		}

		double gi= (derror - tmp_error) / (epsilon);

		cout<<gi<<" "<<w1[i][j]<<endl;

	}
}

TestDemo:

12月 04

PRML Chapter 2.3 The Gaussian Distribution

部分内容参考《Methods of Multivariate Analysis》简称MMA

1. preliminary (Matrix Algebra — Methods of Multivariate Analysis C H A P T E R 2)

加入一些其他内容,调整至http://www.xperseverance.net/blogs/2012/12/1579/

2.多元情况下的协方差

我觉得PRML上2.3中的多元高斯分布为啥让人觉得虎,就是因为相对于单元高斯分布,多元情况下的方差令人迷惑和费解,所以本节只记录《MMA》中第三章讲述的多元协方差矩阵。

bivariate case:现在假设只在二维情况下讨论问题,则每个随机变量表示为$(x_i,y_i)^T$,则两者的协方差:

$\sigma_{xy}=E[(x-\mu_x)(y-\mu_y)]$ 在任何统计或概率书上都有定义。

同时相关系数:$\rho_{xy}=corr(x,y)=\frac{\sigma_{xy}}{\sigma_x\sigma_y}=\frac{E[(x-\mu_x)(y-\mu_y)]}{\sqrt{E(x-\mu_x)^2}\sqrt{E(y-\mu_y)^2}}$

对于两个变量的相关性的理解可以用下面一个例子来描述:

设有二维随机变量(x,y),x表示升高,y表示体重,则凭经验就可以想到,身高是和体重相关的,所以讲这个随机变量画成二维点图应该如下(点集中在两个象限,展现出很高的相关性):

而如果把x换成智力,y换成身高,那么就会变成下面这个样子(所有点四个象限都有,展现出无关性):

3. 多元高斯分布

啊啊待写

4月 09

PRML Chapter 11.2 Markov Chain Monte Carlo

本篇标题虽然是PRML,但是内容来自三本书,另外两本是《Introduction to Probability Models》(以下简称IPM)和《Markov Chain Monte Carlo In Practice》(以下简称MCMCP),都是经典教材,其中前一本内容讲的比较基础,例子很丰富,但是翻译版《应用随机过程-概率模型导论》个人不推荐,翻译较粗糙,只能用来对比中英文数学概念,呵呵。

引用清华大学教授@刘洋THU大牛在微博里的话:“PRML不适合初学者”。我觉得很在理,PRML虽然内容看起来诱人,但真的不太适合偶等初学的,书中有的地方吧,解释的很细腻,但是大部分内容都解释的不够基础和详细。所以光看一本书是不行的……

1. Markov chains

马尔科夫链的基本定义和性质就不废话了,只说一下,一个马尔科夫链可以用如下两项定义:

(1)初始状态概率$p(z^{(0)})$:其中z的上标(i)表示马尔科夫链跳转的第i步,也就是时间i,$z^{(0)}$表示在时间0这个马尔科夫链所在的初始状态,那么$p(z^{(0)})$表示什么意思呢?我的理解大概是比如一个马尔科夫链有3个状态1、2、3,那么$p(z^{(0)})$就是这个马尔科夫链一开始所在的初始状态的概率分布,比如{0.1,0.6,0.3}

(2)转移概率$T_m(z^{m}, z^{(m+1)})\equiv p(z^{(m+1)}|z^{(m)})$:这个PRML上叫做transition probability,MCMCP上叫transition kernel,就是转移矩阵啦。注意PRML上$T_m(z^{m}, z^{(m+1)})$这样的写法可不是表示联合概率,他是表示一个条件概率,表示从状态$z^{(m)}$转移到$z^{(m+1)}$的概率,这个感觉挺容易误会。这里只要记住是从前一个转移到后一个就行。一般随机过程都用$P_{ij}$书,觉得PRML上用$z^{(m)}$和$z^{(m+1)}$更科学,应为i,j是状态,而这里的z是随机变量。

对于上述一个特定的随机变量$z^{(m+1)}$,它的概率等于它所有上一步状态的条件概率取边缘:

$p(z^{(m+1)})=\sum\limits_{z^{(m)}}p(z^{(m+1)}|z^{(m)})p(z^{(m)})$

以前看随机过程都一直没有深入理解好这里,这个式子就是求第m+1步时的随机变量z在某一个状态时的概率。当$m\to\infty$时,它就成了极限概率,也就是通常随机过程书中写的$\pi$,事实上$\pi_j$本来也就是这样求的:$\pi_j=\lim\limits_{n\to\infty}P_{ij}^n=\sum\limits_{i=0}^{\infty}\pi_iPij$,这回要记好了这些个理解以免以后又忘记。

Continue reading

3月 30

PRML Chapter 11.1.1 Basic Sampling Algorithms (for Standard distributions)

        每次看新的部分一开始总会看不懂,然后就开始加快速度看,然后就越看越糊涂。像这部分,现在觉得还算不难,但是第一次看的时候愣是根本没看懂。在重复看了3遍,翻了3本其他的书之后,现在终于把它看懂了,哈哈。现在有战胜它的喜悦感,下次再碰到看不懂的东西有挫折感的时候,不用烦躁,总会看明白的~耐心点就是了。

———————————— 进入正题 ——————————————

        首先,sample算法干的事情就是要产生一些概率分布的具体事件的例子。比如有一个给定参数$(\mu, \sigma^2)$的高斯分布,产生100个具体数据点实例近似符合这个高斯分布,这就是sample算法。最简单的sample算法就是随机数产生算法,C++的函数rand()就干这事,产生[0,1)均匀分布,搜了一下这个,貌似有很多种什么“同余法”、“Mersenne twister”等干这事。Bishop大牛说了,咱不关心这个,假设[0,1)均匀分布已经搞定了,我们产生了一组真正的[0,1)均匀分布的实例z,我们用这组实例z来生成其他概率分布实例。

         要完成上述目标,可以通过一个函数$y=f(z)$来将z转换到一个y的分布上,要做的事情就是要挑选好这函数f(z)以使得转换后得到的y符合我们期待的分布,也就是p(y)符合我们需要sample的概率密度。The distribution of y will be governed by:

$p(y)=p(z)\Big | \frac{dz}{dy} \Big |$      式11.5

把上面这个式子取积分(p(z)=1)得到(这里积分的过程想了半天,还问人问了半天,我决定把微积分重新看看):

$z=h(y)\equiv\int_{-\infty}^{y}p(\widehat{y})d\widehat{y}$    式11.6

从11.6可知$y=h^{-1}(z)$,而函数h是y的分布的CDF函数,其中p(y)是我们想要sample的概率分布密度(the desired distribution),从11.6我们就可以将z的值通过$y=h^{-1}(z)$映射到y的值,并且y的值的分布符合the desired distribution。

下面给出一个例子,假设要我们要sample的分布为exponential distribution,其密度符合:

$p(y)=\lambda exp(-\lambda y)$       $0\le y < \infty$   式11.7

在这个例子里,11.7的积分下界是0,所以$h(y)=1-exp(-\lambda y)$,$h(y)=z$,$y=h^{-1}=-\lambda^{-1}ln(1-z)$,这样只要把每个z的值通过这个映射transform成y,则y的分布就会符合exponential distribution

Continue reading

3月 21

The Dirichlet Distribution 狄利克雷分布 (PRML 2.2.1)

        Dirichlet分布可以看做是分布之上的分布。如何理解这句话,我们可以先举个例子:假设我们有一个骰子,其有六面,分别为{1,2,3,4,5,6}。现在我们做了10000次投掷的实验,得到的实验结果是六面分别出现了{2000,2000,2000,2000,1000,1000}次,如果用每一面出现的次数与试验总数的比值估计这个面出现的概率,则我们得到六面出现的概率,分别为{0.2,0.2,0.2,0.2,0.1,0.1}。现在,我们还不满足,我们想要做10000次试验,每次试验中我们都投掷骰子10000次。我们想知道,出现这样的情况使得我们认为,骰子六面出现概率为{0.2,0.2,0.2,0.2,0.1,0.1}的概率是多少(说不定下次试验统计得到的概率为{0.1, 0.1, 0.2, 0.2, 0.2, 0.2}这样了)。这样我们就在思考骰子六面出现概率分布这样的分布之上的分布。而这样一个分布就是Dirichlet分布。

首先用上面这一段来点直观印象,然后列一些资料:

          维基里面对于狄利克雷分布貌似介绍的挺复杂,不够基础。我找到了一个CMU的PPT:Dirichlet Distribution, Dirichlet Process and Dirichlet Process Mixture,找到一篇华盛顿大学的《Introduction to the Dirichlet Distribution and Related Processes》介绍。

       发现CMU那个ppt里面讲到,Beta is the conjugate prior of Binomial,有一种原来如此的感觉。嗯,原来贝塔分布是二项分布的共轭先验分布,那么狄利克雷分布就是多项分布的共轭先验分布。所以要看狄利克雷分布,就要先了解多项分布,然后呢,想要了解狄利克雷之于多元的关系,就要先看贝塔分布和伯努利分布的关系。所以,二项分布、beta分布、以及共轭这三点是理解狄利克雷分布的关键基础知识,这个基础知识记录在这里(PRML2.1整小章介绍了这个)。

       下面正式进入狄利克雷分布介绍,首先说一下这个多项分布的参数μ。在伯努利分布里,参数μ就是抛硬币取某一面的概率,因为伯努利分布的状态空间只有{0,1}。但是在多项分布里,因为状态空间有K个取值,因此μ变成了向量$\vec{\mu}=(\mu_1,~…,~\mu_k)^T$。多项分布的likelihood函数形式是$\prod\limits_{k=1}^{K}\mu_k^{m_k}$,因此就像选择伯努利分布的共轭先验贝塔函数时那样,狄利克雷分布的函数形式应该如下:

$p(\mu|\alpha)\propto\prod\limits_{k=1}^{K}\mu_k^{\alpha_k-1}$  式2.37

上式中,$\sum_k\mu_k=1$,$\vec{\alpha}=(\alpha_1,~…,~\alpha_k)^T$是狄利克雷分布的参数。最后把2.37归一化成为真正的狄利克雷分布:

$Dir(\mu|\alpha)=\frac{\Gamma(\alpha_0)}{\Gamma(\alpha_1)…\Gamma(\alpha_k)}\prod\limits_{k=1}^{K}\mu_k^{\alpha_k-1}$

其中$\alpha_0=\sum\limits_{k=1}^{K}\alpha_k$。这个函数跟贝塔分布有点像(取K=2时就是Beta分布)。跟多项分布也有点像。就像Beta分布那样,狄利克雷分布就是它所对应的后验多项分布的参数$\vec{\mu}$的分布,只不过μ是一个向量,下图是当$\vec{\mu}=(\mu_1,\mu_2,\mu_3)$时,即只有三个值时狄利克雷概率密度函数的例子。其中中间那个图的三角形表示一个平放的Simplex,三角形三个顶点分别表示$\vec{\mu}=(1,0,0)$,$\vec{\mu}=(0,1,0)$和$\vec{\mu}=(0,0,1)$,因此三角形中间部分的任意一个点就是$\vec{\mu}$的一个取值,纵轴就是这个$\vec{\mu}$的Simplex上的概率密度值(PDF)。

对于参数$\vec{\mu}$的估计时,可知 后验=似然*先验 的函数形式如下:

$p(\mathbf{\mu}|D,\mathbf{\alpha})\propto(D|\mathbf{\mu})p(\mathbf{\mu}|\mathbf{\alpha})\propto\prod\limits_{k=1}^{K}\mu_k^{\alpha_k+m_k-1}$

从这个形式可以看出,后验也是狄利克雷分布。类似于贝塔分布归一化后验的方法,我们把这个后验归一化一下,得到:

$\begin{equation}\begin{split}p(\mathbf{\mu}|D,\mathbf{\alpha}) & =Dir(\mathbf{\mu}|\mathbf{\alpha}+\mathbf{m}) & =\frac{\Gamma(\alpha_0+N)}{\Gamma(\alpha_1+m_1)…\Gamma(\alpha_K+m_K)}\prod\limits_{k=1}^{K}\mu_k^{\alpha_k+m_k-1}\end{split}\end{equation}$

3月 17

PRML Chapter 9.2.1 Maximum likelihood (for GMM – Gaussian mixture model)

设观察到$\{x_1,…,x_N\}$共N个数据点独立产生自混合高斯分布,则这些独立同分布数据(i.i.d 意思就是后一个数据$x_{i+1}$的产生与前一个数据$x_i$独立)的概率图模型可以用图9.6表示

这些数据的似然函数可以表示为:

$ln~p(X|\pi,\mu,\Sigma)=\sum\limits_{n=1}^{N}ln\{\sum\limits_{k=1}^{K}\pi_kN(x_n|\mu_k,\Sigma_k)\}$

关于似然函数这玩意,我打算下次忘记的时候再翻翻统计学,然后记下来,我现在还记得是怎么回事,暂时先不多写了。

P433提到用最大似然估计GMM(Gaussian mixture model)的参数会有一个严重的问题,就是奇异点的问题(the presence of singularities)。设当第j个component的均值$\mu_j$正好出现在观察数据中时,即$x_n=\mu_j$,那么这个$x_n$就会成为一个奇异点,此时这第j个component就会“崩塌”到这个点上($x_n=\mu_j$均值相等,而方差$\sigma_j\to 0$)。书上说这也是一个最大似然估计方法严重over-fitting的例子。(3@21补充如下)今天回看第二章的内容时,在P70页提到了另一个over-fitting的例子,就是当观察数据量(训练数据)很少时,最大似然估计会严重过拟合,比如估计硬币正面的概率,当抛硬币很少次而正好全部正面时,会出现概率为1的情况。

按照我的理解,不是因为第j个component的均值正好出现在数据中导致崩塌,情况应该是反过来的,只要存在奇异点,那么在最大似然估计方法估计参数的时候,其中某一个component就会因为存在这样一个奇异点而在估计参数时崩塌到这个点上。

值得一提的是单高斯分布就不会出现这样的问题,因为其他数据点会阻止崩塌(这里看得不是很明白,所以只能这样理解),而一旦有多个高斯分布混合的时候,假设就两个高斯分布a和b,那么其他数据点就会都分布到a上面去,而只留下一个奇异点分给b,而导致崩塌。图9.7正是说明了这一点。

避免崩塌的方法是当发现一个component的方差变得太小的时候给它重置一个随机的均值,并且把方差改大点。

3月 16

PRML Chapter 9.2 Mixtures of Gaussians 混合高斯分布简介

单元高斯分布(The univariate Gaussian),我们高中时就知道了,其表达式如下:

$N(x|\mu , \sigma^2)=\frac{1}{{(2\pi \sigma^2)}^{1/2}}exp\{-\frac{1}{2\sigma^2}(x-\mu)^2\}$

而多元高斯分布(Multivariate Gaussian Distribution)就是有多个参数控制的高斯分布,其均值是一个均值向量$\mu$,设均值向量维度为D,而方差则是方差矩阵$\Sigma$,因此其表达式如下:

$N(x|\mu, \Sigma)=\frac{1}{(2\pi)^{D/2}}\frac{1}{|\Sigma|^{1/2}}exp\{-\frac{1}{2}(x-\mu)^T\Sigma^{-1}(x-\mu)\}$

书中P84,P111对于单个的高斯分布对数据建模的缺点作了描述。由于单个高斯模型是一个unimodal,意即模型只有一个最大值的(看到某页面上称之为“单模态”)(P84:has a single maximum),所以如P111所说的,单个高斯模型表示能力很有限,当数据有两个峰时就没法很好表示啦。因此需要混合高斯分布。

设有K个高斯分布混合在一起,每一个高斯分布称为一个Component,那么混合高斯分布表达式为(这个表达式来自P111):

$p(x)=\sum\limits_{k=1}^{K}{\pi}_kN(x|\mu_k,\Sigma_k)$

其中$\pi_k$称为mixing coefficients,并且$\sum\limits_{k=1}^{K}\pi_k=1$

P430引入了一个二元随机变量z,它的表示方法跟多项分布的表示方法一样,都是1-of-K representation,即$z_k\in\{0,1\}$并且$\sum_{k}z_k=1$ (意即对于K维的向量z只有一个值是1,其他所有值都是0),这样就可以用z来表示$\pi_k$,它们之间的关系是$p(z_k=1)=\pi_k$。所以z的概率分布类似于多项分布的表示,即

$p(z)=\prod\limits_{k=1}^{K}{\pi_k}^{z_k}$          式9.10

$(=\pi_k 当对应k的那个z_k是1时,因为其他指数都是0)$

已知z之后,x对于z的条件分布为

$p(x|z_k=1)=N(x|\mu_k,\Sigma_k)$     式9.11

根据P431页图模型(z->x),x与z的联合概率为p(x,z)=p(z)p(x|z),而将这个联合概率对x取边缘概率,就得到了如下(结合9.10和9.11):

$p(x)=\sum_{z}p(z)p(x|z)=\sum\limits_{k=1}^{K}{\pi}_kN(x|\mu_k,\Sigma_k)$

这就是混合高斯分布表达式的推导过程。之所以要引入这个隐含变量z是为了引入EM算法,这个之后再说。

另一个重要变量是$\gamma(z_k)$,表示当观察到x时它由第k个Component产生的概率,即$\gamma(z_k)=p(z_k=1|x)$,更具体的

$\gamma(z_k)=p(z_k=1|x)=\frac{p(z_k=1)p(x|z_k=1)}{\sum\limits_{j=1}^{K}p(z_j=1)p(x|z_j=1)}$

$=\frac{\pi_kN(x|\mu_k,\Sigma_k)}{\sum\limits_{j=1}^{K}\pi_jN(x|\mu_j,\Sigma_j)}$    式9.13

这个式9.13后面介绍EM算法的时候要用到,这里先写一下。

P432讲的是从混合高斯分布Sample的过程,使用了图模型那章讲的ancestral sampling(P365)(其实就是根据图模型的拓扑结构依次sample随机变量值)。混合高斯分布的sample过程很简单(z->x),现根据p(z)选择z,即选择第$z_k$个Component,此时已知z,然后再sample p(x|z),即式9.11啦,所以此时从第k个Component Sample 出来的x值已经和z无关了(z已知)

呵呵,整理了一下这篇博客之后,关于这些个内容确实清晰了很多,要坚持写下去!

3月 07

PRML Chapter 9. Mixture Models and EM

今天从网上搜的 EM 算法的原始论文叫做 Maximum likelihood from incomplete data via the EM algorithm,下次仔细学习的时候可以看看,最近比较功利,就不弄得那么明白了。(2012@3@21)

9.1 K-means Clustering

主要介绍了 K-means 和 EM 算法之间的关系,第一次听说原来 K-means 就是 EM 算法,不知道的东东还真是多。

接下来第二小结是 Mixtures of Gaussians,看来不得不回2.3去看了。

9.2 Mixtures of Gaussians

这章内容太多,所以单独写成分日志。

第一部分链接:Mixtures of Gaussians

第二部分链接:9.2.1 Maximum likelihood

9.2.2 的推导看不明白,直接上结论

1. 对于GMM的EM算法,由于EM计算量大而且收敛慢,可以先用K-means找到几个中心,然后均值、方差和混合参数都可以根据K-means的结果设定初始值。

2. 需要防止模型在奇异点崩塌

3. EM不能保证找到全局最优值

下面不打字了,直接把书上截图了:

 

3月 05

PRML Chapter 8.2.1 Three example graphs

 

8.2.1 Three example graphs

书中说明了三种基本条件独立的结构

(1) tail-to-tail   (从a,b的角度看,a和b都是c的tail)

如图 Figure 8.15

整个图表示abc的联合概率 p(a,b,c),此时可以根据图8.15的结构推出 p(a,b,c) = p(a|c)p(b|c)p(c)

如果同时两边对c取边缘概率(marginalizing both sides with respect to c),则可以得到式8.24

$p(a,b)=\sum\limits_cp(a|c)p(b|c)p(c)$

通常情况下,式8.24右边那部分是不会等于p(a)*p(b)的,所以a,b是不相互独立的。

但是如果c是作为输入而被已观察到了,情况就变成了图Figure 8.16

此时,整个图如果要关心a,b之间的联合概率的话,就要表示成 p(a,b|c)了,因为c已被观察到,所以a,b之间的联合概率是基于c的条件概率。

所以 $p(a,b|c)=\frac{p(a,b,c)}{p(c)}=p(a|c)p(b|c)$ 因此由这个式子可以看出a,b是关于c条件独立的

 

(2) head-to-tail  (从a,b的角度看,a是c的head,b是c的tail)

如图 Figure 8.17  变化至  8.18

   ->        

对于图8.17,$p(a,b,c)=p(a)p(c|a)p(b|c)$

两边对c取边缘概率得  $p(a,b)=p(a)\sum\limits_{c}p(c|a)p(b|c)=p(a)p(b|a)$  所以a,b不独立。

而对于图8.18,可知:

$p(a,b,|c)=\frac{p(a,b,c)}{p(c)}=\frac{p(a)p(c|a)p(b|c)}{p(c)}=p(a|c)p(b|c)$,因此图8.18中a与b相互独立。

 

(3) head-to-head (从a,b的角度看,a和b都是c的head)

如图8.19变化至8.20

 -> 

        对于图8.19,p(a,b,c)=p(a)p(b)p(c|a,b),两边对c取边缘直接可以得到 p(a,b)=p(a)p(b),因此这个时候没有观察到c反而倒是a,b相互独立,而如果观察到c呢?就变成了图8.20中的p(a,b|c)=p(a,b,c)/p(c)=p(a)p(b)p(a,b|c)/p(c),不等于p(a|c)p(b|c),因此a和b反而不独立。

        由以上三个基本图可以看出,三个随机变量a,b,c的联合概率的分解方法实际上是和它们之间的依存关系有关的,不能随便分解,如果分解成p(a,b,c)=p(a,b|c)p(c),那么就默认了tail-to-tail这种结构,如果分解成其他,则有其他相应的结构。同样关于a,b,c三者之间的概率,不同结构下,观察到c的情况下,a与b可能独立,也可能不独立!

P377

提供了一个例子,说明head-to-head的情况下,两个图8.20中a,b是如何在已知c的情况下相互影响的。

 

3月 02

PRML Chapter 2. Probability Distributions

注:mathjax公式显示常会有问题,若遇到问题,请刷新或换个浏览器查看。

P68

conjugate priors

      In Bayesian probability theory, if the posterior distributions p(θ|x) are in the same family as the prior probability distributionp(θ), the prior and posterior are then called conjugate distributions, and the prior is called a conjugate prior for the likelihood. For example, the Gaussian family is conjugate to itself (or self-conjugate) with respect to a Gaussian likelihood function: if the likelihood function is Gaussian, choosing a Gaussian prior over the mean will ensure that the posterior distribution is also Gaussian.

exponential family

      The exponential families include many of the most common distributions, including the normalexponentialgammachi-squaredbetaDirichletBernoullibinomial,multinomialPoissonWishartInverse Wishart and many others.

2012@3@21补充:起先第二章果真没有仔细看好,现在看到狄利克雷分布了,回来再了解一下贝塔分布,发现里面这么关键的内容愣是都没看出来,汗。不过今天22点半了,明天在写咯。

2.1.1 The beta distribution

    如果忘记伯努利分布和二项分布是怎么回事了,看这里

    书中引出贝塔分布的理由:P70提到,由于最大似然估计在观察数据很少时,会出现严重over-fitting(比如估计抛硬币正反面概率,只有3次抛硬币观察数据,且结果正好都是正面,则模型预测以后所有抛硬币都将是正面)。为了解决这个问题,可以考虑贝叶斯方法,即引入一个先验知识(先验分布p(μ))来控制参数μ,那么如何挑选这个分布呢?

   考虑到伯努利分布的似然函数的形式是$\mu^x(1-\mu)^{1-x}$错!!原先这里看了个似懂非懂,完全写错了,囧死了,得到一个教训,写日志还是要多来回看看,看懂了再写,否则留下笑柄!现在重写如下:应该是,二项分布的似然函数是:$\mu^m(1-\mu)^n$ (就是二项分布除归一化参数之外的后面那部分,似然函数之所以不是pdf,是因为它不需要归一化),这个函数的形式是$\mu$的m次方乘以$1-\mu$的n次方,记住这个形式,下面要用到。

根据:

$posterior=likelihood * prior$

已经知道似然函数的形式,如果选择的先验分布也与 μ 和 (1-μ) 两者的乘方成比例,那么后验分布的函数形式就会跟它的先验函数形式一样了。具体来说,选择prior的形式是$w_1*\mu^a(1-\mu)^b$,那么posterior就会变成$w_2*\mu^{m+a}(1-\mu)^{n+b}$这个样子了($w_1,w_2$为pdf的归一化参数),所以posterior和prior具有相同的函数形式(都是μ和(1-μ)的次方的乘积),这就是所谓的conjugacy。

    最终这里的先验和后验就都是贝塔分布了,其中先验的形式如下:

$Beta(\mu|a,b)=\frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)}\mu^{a-1}(1-\mu)^{b-1}$  式2.13

其中$\frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)}$这玩意就是w1,是为了把整个分布概率归一化,从而使:

$\int_{0}^{1}Beta(\mu|a,b)d\mu=1$    式2.14

在维基里面,有这么一个式子:

$B(\alpha,\beta)=\frac{\Gamma(\alpha)\Gamma(\beta)}{\Gamma(\alpha+\beta)}\simeq {\alpha-1 \choose \alpha+\beta-2}$

瞬间觉得世界清晰了,因为$\Gamma(n)=(n-1)!$,所以其实当上式中$\alpha,\beta$为整数时,就是$C^{\alpha-1}_{\alpha+\beta-2}$。因此,其实beta分布就是二项分布推广成实数域上的情况而已!注意,这里曾经把Beta函数写反过,Beta function 是指$B(x,y)=\frac{\Gamma(x)\Gamma(y)}{\Gamma(x+y)}$,而Beta distribution的pdf公式为$Beta(\mu|a,b)=\frac{1}{B(a,b)}\mu^{a-1}(1-\mu)^{b-1}$

从式2.14看出,Beta分布就是一个μ的PDF(概率密度函数)(这个昨天(3@21)刚仔细看过哈),μ本身是二项分布的参数,而a,b由于2.14的归一化过程可以被视作μ的控制参数,因此贝塔分布的a和b就被称作hyperparameters。下面的图是Beta分布的几个例子,其中横轴是μ取值范围,纵轴是PDF取值,PDF的值可以大于1哦。

最后得到的posterior如下:

$p(\mu|m,l,a,b)\propto\mu^{m+a-1}(1-\mu)^{l+b-1}$  式2.17,其中l=N-m

要把这个posterior归一化其实可以参照式2.13,式2.17中的m+a等同于2.13中那个a,而l+b就是2.13中那个b,所以:

$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}$

最后,如果我们已经有观察数据D,要估计μ,即$p(\mu|D)$,我们可以得到:

$p(x=1|D)=\frac{m+a}{m+a+l+b}$  式2.20

可以发现这个式子比最大似然估计的结果m/(m+l)多了a和b,也就是先验知识的影响。

2012@4@17补充:事实上这种给参数加先验的方法叫做Maximum a posterior(MAP),而PRML里的这章没有明显的提到,刚从《Parameter estimation of text analysis》里看到,在这里提一下

2013@4@22 在评论中有同学指正,2.20的结果不应该是MAP,输出期望的应是Bayesian inference。MAP结果应为(m+a-1)/(m+a+l+b-2)。当时并没有搞清楚这个两者的区别。在此纠正。非常感谢评论中的同学!!鉴于经常被发现错误,我打算近期写完论文检查一下所有博客。

2.2 Multinomial Variables

        Multinomial Variables说白了就是多种选择选其一。比如随机变量X有三种取值x1,x2,x3,那么用一个三维向量表示Multinomial 的取值就是{1,0,0},{0,1,0},{0,0,1}分别代表选中x1,x2,x3,即必须选中一个,同时只能选一个这样的意思。

如果用μk表示xk=1时的概率,那么对于随机变量x的取值的概率分布可以表示为:

$p(x|\mu)=\prod\limits_{k=1}^{K}\mu_{k}^{x_k}$

        其实这个式子的意思就是当K取值k的时候,只有xk是1,其他都是0,所以这个p(x|μ)的值就是μk的值而已,因为一个数的0次方是1,所以对于其他xi(i≠k)的那部分μi全部都乘以了一个1而已。搞了这么一个玄乎的式子,应该是为了数学表示全面点,事实上直接理解就是p(x|μ) = μk。

       上面所讲的这些其实只是多项分布的一次事件(或一次观察),如果有N多次观察,那么就需要用多项分布来描述了。就像伯努利分布只是描述一次抛硬币,而二项分布是描述N次抛硬币的一样。

        对于Multinomial 的极大似然估计其实可想而知,就是数数xk的个数然后取占整个集合的比例作为概率了。式(2.31)给了数学上的likelihood的式子,但是那个什么拉格朗日乘子λ我已经没啥概念了,只知道是用来求函数极值的,这里记着点以后到高数里去看2012@4@4补充,大致看了一下拉格朗日乘数法,没有想象中的复杂,就是用来求一个条件极值,在这里

2.3 The Gaussian Distribution

        这章的公式太唬人了,西方人的数学功底感觉很好,虽然这些数学都挺简单,但是我没有全局了解的那种观念,所以有些东西没法一下子理解了。多元高斯分布里面的矩阵的东西不太好懂。

       P80讲的是多元高斯分布通过正交向量$u_i$ 线性变换为以 $y=U(x-\mu )$   为基的向量空间,大概也只能明白这点了。在这样一个基上,高斯分布是一个椭圆(当所有特征值$\lambda _i$都是正数的时候)

       P82有讲到关于 moment 这个就是矩,k阶中心矩(The kth central moment)为 $E((X-\mu )^k)$

       P84讲高斯分布的缺点,,然后说到可以引入隐含随机变量来解决高斯分布模型表示能力不强的问题。可以用一系列离散的隐含变量混合高斯分布,可以用连续隐含变量,用来干嘛不太理解。

       再后面已经完全不知所云了。所以只能Skip了。

2012@11@30

时隔8月,偶终于回来把这些东西都看懂了!其实并没有那么虎,只不过PRML讲得不是那么细而已。只要翻《Methods of Multivariate Analysis》,这些东西其实相对来说还是easy的!链接

2.5 Nonparametric Methods

      这章主要介绍两种无参方法:核方法和近邻法。

       P122 开始介绍核密度估计时,从公式2.242到2.246都是为了推导未知概率密度p(x)的估计。

      最后推导得到式2.246如下:

$p(x)=\frac{K}{NV}$

      其中V是x附近区域R的体积,K则是落入x附近区域R中的数据点个数,由此导出了两种不同的密度估计方法:

      (1)如果固定K而估计V的大小,那么就是kNN算法(k固定而根据选定最近的k个数据来评估R的体积)

      (2)如果固定V而估计K的大小,那么就是核密度估计,用一个核函数表示一个固定的体积V,然后数数这个体积里面数据点K的个数

     这就是对这两种无参数方法比较深入的理解,很好。

      关于核密度估计,虽然还不是很清楚,但是可以知道其实它的道理跟P121的histogram approach是一样的,只不过核密度估计是高维的而已

 

           2012年3月15日,大致模糊的看完第二章