星期四, 八月 30, 2012

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来实现这种思路:

星期六, 八月 25, 2012

笨办法学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))

星期三, 八月 22, 2012

笨办法学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)]

星期一, 八月 20, 2012

笨办法学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)
}

星期日, 八月 19, 2012

新书推荐:数据新闻手册



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

星期六, 八月 18, 2012

笨办法学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
}

星期五, 八月 17, 2012

笨办法学R编程(1)


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

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

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

Let's Go

星期一, 八月 13, 2012

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
整理后的数据文档

代码

星期五, 八月 10, 2012

在R语言中绘制Treemap

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


星期六, 八月 04, 2012

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

由美国耶鲁大学和哥伦比亚大学联合推出的“年度全球环境绩效指数”(EPI)排名在2012年年初放出。此排名旨在评估一个国家的环境政策,环境卫生与生态系统之平衡的状态,涵盖10项领域共22项环境指标,涉及领域包括气候变化、农业、渔业、森林、水源、空气污染及环境负担等。详细报告和数据都可以在网站主页找到。本文就使用这个数据集来观察中国的环境状况在全球的位置。
上面的小提琴图即是22项指标的分布,叠加的红色圆点表示了中国得分位置。这些指标已经被研究人员进行了转换,所有的指标均是得分越高越好。从这些得分点可以看出,中国在森林保护(FORGROINV)和生物多样性(PACOV)方面做得还不错,CO2排放(CO2GDP)和饮用水(WATSUPINV)等方面做得一般,在PM2.5和室内空气(INDOOR)等方面处于垫底情况。更多的指标说明请参看报告主页。其中的EPI是环境绩效评价的综合得分。