EM算法的R实现和高斯混合模型

EM（Expectatioin-Maximalization）算法即期望最大算法，被誉为是数据挖掘的十大算法之一。它是在概率模型中寻找参数最大似然估计的算法，其中概率模型依赖于无法观测到的隐变量。最大期望算法经过两个步骤交替进行计算，第一步是计算期望（E），也就是将隐藏变量象能够观测到的一样包含在内，从而计算最大似然的期望值；另外一步是最大化（M），也就是最大化在 E 步上找到的最大似然的期望值从而计算参数的最大似然估计。M 步上找到的参数然后用于另外一个 E 步计算，这个过程不断交替进行。对于信息缺失的数据来说，EM算法是一种极有效的工具，我们先用一个简单例子来理解EM算法，并将其应用到GMM（高斯混合模型）中去。

幼儿园里老师给a，b，c，d四个小朋友发糖吃，但老师有点偏心，不同小朋友得到糖的概率不同，p(a)=0.5, p(b)=miu, p(c)=2*miu, p(d)=0.5-3*miu 如果确定了参数miu，概率分布就知道了。我们可以通过观察样本数据来推测参数。设四人实际得到糖果数目为(a,b,c,d)，可以用ML（极大似然），估计出miu=(b+c)/6*(b+c+d)，假如某一天四个小朋友分别得到了(40,18,0,23)个糖。根据公式可求出miu为0.1。在信息完整条件下，ML方法是很容易估计参数的。

假设情况再复杂一点：知道c和d二人得到的糖果数，也知道a与b二人的糖果数之和为h，如何来估计出参数miu呢？前面我们知道了，如果观察到a，b，c，d就可以用ML估计出miu。反之，如果miu已知，根据概率期望 a/b=0.5/miu，又有a+b=h。由两个式子可得到 a=0.5*h/(0.5+miu)和b=miu*h/(0.5+miu)。此时我们面临一个两难的境地，a，b需要miu才能求出来，miu需要a，b才能求出来。有点类似岸上的八戒和河里的沙僧的对话：你先上来，你先下来，你先上来......

那么一种思路就是一方先让一步，暂且先抛出一个随机的初值，然后用对方算出的数值反复迭代计算。直到计算结果收敛为止。这就是EM算法的思路，我们来看看用R来实现这种思路：

笨办法学R编程（5）

随着教程推进，基本的语法都接触得差不多了。当要解决某个具体问题时，只需要考虑用什么样的算法来整合运用这些函数和表达式。今天来解决Project Euler的第五个问题，该问题可以用很笨的暴力搜索法子来作，但是更聪明的作法是采用质因子分解的思路。即任何一个合数都可以分解为质数的乘积。为了完成这个题目，还需要学习一点点矩阵，以及和sapply函数相似的另一个函数apply。

# 预备练习
mat <- matrix(1:12,ncol=4)
print(mat)
t(mat)
colnames(mat) <- c('one','two','three','four')
rownames(mat) <- c('a','b','c')
print(mat)
apply(mat,1,sum)
apply(mat,2,sum)
sum(apply(mat,2,sum))
prod(apply(mat,2,sum))

笨办法学R编程（4）

看到各位对“笨办法系列”的东西还比较感兴趣，我也很乐意继续写下去。今天的示例将会用到数据框（data.frame）这种数据类型，并学习如何组合计算两个向量，以及如何排序。我们将用所学的东西来解决Project Euler的第四个问题，就是找出一个集合中最大的回文数。回文数是指一个像1534351这样“对称”的数，如果将这个数的数字按相反的顺序重新排列后，所得到的数和原来的数一样。开始啦！

# 预备练习
 
x <- y <- 1:9
data <- expand.grid(x=x,y=y)
print(data)
z <- data$x * data$y
# 一个九九乘法表
z <- matrix(z,ncol=9)
 
set.seed(1)
x <- round(runif(10),2)
print(x)
order(x)
x[order(x)[1]]
which.min(x)
x[which.min(x)]
x[order(x)]

笨办法学R编程（3）

经历了前面两个小挑战，你应该对R有点理解了。我们继续推进，今天的问题有点点复杂，复杂的不是R，而是一个数学概念：质数和质因子。任何一个合数都可以被几个质数所分解，这个性质很重要，我们将用它来解决Project Euler的第三个问题。还是和之前一样的，你需要自己在R控制台中敲打下面这些命令，根据结果自行揣摩其用处。

# 预备练习，学习for循环、建立自定义函数和其它一些函数
 
for (n in 1:10) {
    print(sqrt(n))
}
 
x <- c('hello','world','I','love','R')
for (n in x) {
    print(n)
}

新书推荐：数据新闻手册

如果你和我一样是数据爱好者，那么一定会经常造访卫报的数据博客栏目，去看看他们是如何用不同来源的数据制作有趣的新闻。你一定想知道，他们是如何提出问题的？如何采集数据？又如何整理和分析数据，并最终用可视化方法来展现到读者面前的。答案就在《数据新闻手册》这本书中。

笨办法学R编程（2）

本例将介绍R语言中的while循环和if条件。最终用它来解决Project Euler的第二个问题。除了练习之外你还需要了解一些斐波纳契数列的知识。废话不多说了，打开R控制台，跟着输入下面的代码，自行琢磨吧。

# 预备练习，while循环和if判断
x <- 1:10
print(x)
print(x[10])
print(x[-10])
 
i <- 1
while (i <= 10) {
  print(x[i])
  i <- i + 1
}

笨办法学R编程（1）

在倚天屠龙记中，有一人唤作火工头陀。此人练功不靠心法，只靠模仿他人招式，由外而内，自成一家。练习编程也有如此的法门，不看文字描述，只观察和模仿别人的代码。这样也可以由外而内学会编程。《笨办法学python》的作者Zed Shaw 就说过这种笨办法入门其实更简单。阳志平在他的文章《如何学习一门新的编程语言》中也讲到，初学编程要在学习区刻意的大量练习，少看理论书。

TED上一位教育家同样谈到这么一个故事，他把一个计算机扔在一个偏远的印度小村子里不去管它，在那里没有上过学的小孩就能自己学会英语和计算机的用法。实际上人脑是非常善于自我探索和学习的。因此本系列教程的特点就是只有演示代码加少量注释。通过反复模仿和练习，揣摩代码的变化和结果，你就能自行领悟其含义，并打下坚实的编程基础。

本系列每篇文章的目的都是用R语言编程来解决一个Project Euler的问题。Project Euler是一系列由易到难的计算机编程挑战，它提供了一个平台来激发我们解决问题的灵感和思路。本人写这个教程的目的有三：一是为了好玩，二是提高编程水平，三是示范说明以提供给需要的R初学者。另外从R-Blogger上了解，已经有两位高人用R在计算Project Euler，各位也可以参照他们的文章（博客1、博客2）。

Let's Go

2012伦敦奥运的金牌项目分布

伦敦奥运会尘埃落定，又到了数金牌的时候了。上图还是用的R中的treemap包绘制了这个大game中各项目的金牌分布情况，矩形的大小反映了项目的重要性。比如说在所有302块金牌中，Athletics（田径）所占的金牌比例就很高。与之相当的则是Aquatics（水上项目）。而中国所获金牌数则用不同的颜色深度来表示，例如说在Badminton（羽毛球）上拿到了5枚金牌，颜色就非常深。可以观察到在很多项目中，中国仍是一片空白。

参考资料：
http://www.london2012.com/country/china/medals/index.html
http://en.wikipedia.org/wiki/2012_Summer_Olympics#Sports
整理后的数据文档

代码：

在R语言中绘制Treemap

Treemap是一种流行的可视化技术，它用不同尺寸的嵌套矩形来表现层次数据。最常见的层次数据包括了文件目录、文章结构、组织结构等等。Treemap技术是通过矩形的大小来显示节点的重要性，并通过嵌套层次来显示结构的层级。R语言中的treemap包就可以实现这种可视化技术，下图就是对苹果公司财务报表的可视化。

中国的环境状况在全球的位置

由美国耶鲁大学和哥伦比亚大学联合推出的“年度全球环境绩效指数”（EPI）排名在2012年年初放出。此排名旨在评估一个国家的环境政策，环境卫生与生态系统之平衡的状态，涵盖10项领域共22项环境指标，涉及领域包括气候变化、农业、渔业、森林、水源、空气污染及环境负担等。详细报告和数据都可以在网站主页找到。本文就使用这个数据集来观察中国的环境状况在全球的位置。

上面的小提琴图即是22项指标的分布，叠加的红色圆点表示了中国得分位置。这些指标已经被研究人员进行了转换，所有的指标均是得分越高越好。从这些得分点可以看出，中国在森林保护(FORGROINV)和生物多样性(PACOV)方面做得还不错，CO2排放(CO2GDP)和饮用水(WATSUPINV)等方面做得一般，在PM2.5和室内空气(INDOOR)等方面处于垫底情况。更多的指标说明请参看报告主页。其中的EPI是环境绩效评价的综合得分。