Beta

It’s a beautiful thing when free data meets free algorithm.

话说到,去年在上海财经大学讲 R与统计图形时,提到了Edward Muybridge (1830-1904)的赛马动画。在准备材料的时候,我也比较八卦的翻了翻关于赛马动画的历史,结果发现:这幅图型不但是统计动画的鼻祖,同样是现代电影的先驱。

从 Edward Muybride 拍摄赛马动画后,美国的电影产业开始高速的发展,从此加利福尼亚州顺理成章地成为人类电影发展上的重镇,加州的好莱坞产生了大量的电影技术的创新,好莱坞电影也成为美国文化的主要代表之一。

电影、动画的原理,我就不八卦了,一般理科生大概都有些了解。关于这个赛马动画的产生,很有意思:

  • 1872年,前美国加州州长 Leland Stanford(也是斯坦福大学的创立者)是一个狂热的赛马爱好者,为了证明马在奔跑的时候会有一刻所有的蹄子同时悬空,和人打赌,赌金非常高,达到了$25,000。
  • 而在那个年代很难用肉眼确定马在奔跑时的状态(可以想象一下为什么马踏飞燕是那个样子?)。于是 Stanford 找到并雇佣 Muybridge 这个摄影家帮他解决这个问题。
  • Muybridge 本来在 1872 年的时候已经接受了 Stanford 的邀请,为 Stanford 提供那旷世赌博的摄像证据。但这家伙怀疑自己老婆有个情人Larkyns,并且冲动的枪杀了 Larkyns。
  • 直到1877年,Muybridge 被判无罪(Stanford 提供的辩护资助),才又继续他的奔马实验,于是有了这个:

后来 Muybridge 根据这些赛马的图片,创作了人类历史上的第一个小电影。下面这个动画就是用最上面的几张图合并而成的(因为偷懒用ImageMagick自动切割,所以这个小电影有点晃~~)

当然,还有一个效果更好的:

哈,这便是统计动画的原理,也是现代电影的最初版本~~

最近一个月所有的业余时间基本上都是在准备会议材料,开会中度过的。且不说上周去杭州,被客户折磨了3天,居然临走也没有时间去西湖逛一逛,稀里糊涂的又从杭州赶回北京,参加大数据技术大会。

工作会议就暂且不提了,说说最近参加的两场学术性质的会议。

第一场是在11月12日-13日,上海召开的中国R语言会议,本次会议各种重量级嘉宾:数据挖掘领域鼻祖级人物——谢邦昌,中国最年轻的教授——周涛,来自于奥克兰大学VGAM包的作者THomas W. Yee,等等guru。本着R语言会议一场不能落的原则,于是乎贡献了一篇 Data Mining With RWeka,主要讲了一下关于R语言和Weka之间的算法调用,以及附带讲了一下分类算法的评估。不过很可惜,本来要和舰哥仔细聊聊的,结果时间没排开,只好再找机会了。

第二场是在北京11月26日,由CSDN召集的大数据技术大会,和大家分享了一下R语言主题的《R You Ready?》,没想到结束之后又多位同行对R语言很感兴趣。这次会议有些主题非常有意思,比如淘宝数据团队负责人-赵昆的汇报《淘宝海量数据技术》,演示了一张淘宝卖家卖家交易客户的地域可视化图,等有时间也实现一下。

明天是年终工作总结会,又是一天的会议~~

今天看到老同学@JulieJulieJulieJulie浪漫求婚,真的很浪漫、很唯美、很感动。正如评论说的,我们又相信爱情了!于是,小兴奋,睡不着,爬起来补一篇文章。


正文开始

最近在数据挖掘专业网站 KDnuggets 上刊出了2011年度关于数据挖掘/分析语言流行度的调查,不出意料R、SQL、Python果然排在了前三位。当然有看官说了,参与调查的样本数量太少,而且以登录KDnuggets网站的用户为主,样本的信息显然是有偏的。实话说,我也对KDnuggets网站的Poll持保留态度,但它的结果毕竟代表了某一类人群的使用偏好,尤其是在语言角度。

我们看排名前5位的语言:

  • R:世界范围内的标准统计语言,以快速更新的算法,灵活的编程,广泛的扩展,绚丽的图形著称,遵循GPL协议的开源软件
  • SQL:大部分企业使用的,数据仓库、集市的通用查询语言,在大型数据应用上有极大的优势,同时也是数据分析/挖掘的基础
  • Python:传说中的Google的三大开发语言,适用于粘合一些复杂应用,我这里工作暂时没有涉及过
  • Java:太多的应用都基于Java的,不然Oracle也不会花上74亿美元收购SUN了
  • SAS:曾经的数据分析领域老大,当然现在市场份额依旧非常高。但SAS昂贵的使用费用迫使更多的分析工作者转到了开源领域,比如R

后四种语言同R语言还都有一些关系,闲扯起来还真是没完没了,这里就不再赘述,各位可以在搜索引擎上搜索R+XXX。 如果我们将范围限制在数据挖掘这个主题,R同SQL的关系则变得非常非常紧密。

众所周知,R的强项在于灵活的算法,以及开发速度,但其所有的计算都是在内存中进行,一旦数据量达到了内存上限,基本上就是叫天天不灵, 叫地地不应了。所以在使用R做数据挖掘时,就必须考虑使用其他的数据工具弥补R在这方面的劣势。尤其是在商业应用上, 不能搭建R环境的条件下,SQL语言是提供挖掘结果的不二选择。

支持SQL的商用数据库比如Oracle、DB2性能优异,但对系统的占用非常厉害,假如本地装了Oracle,又开了点其他应用,2G的内存很快就会吃到1.5G甚至以上, 再想用R做分析那只能用“捉襟见肘”这个词来形容了。当然如果在办公条件下有相应的服务器环境最好, 在某些应用环境下,甚至可以通过本地多开R进程来达到并行计算的目的。

或者本地分析比较多,但数据量又时常上到百兆,虽然R也能够处理,但依然建议将数据移植到本地构建的轻量数据库环境,比如MySQL环境。 从我的经验上看,虽然MySQL对比Oracle、DB2来说小巧很多,但在同R语言配合的本地应用上,性能更加有保证。

有了支持SQL的数据库环境,就要聊一聊R语言到底和SQL有什么关系:

  • 各大数据库厂商已经开发了相关的支持R语言的数据挖掘套件,比如Oracle的RODM,Teradata的 teradataR等。
  • R本身就可以通过扩展包来对数据库执行SQL,这时你可以把R语言作为调度环境。R的计算过程结果可以直接作为参数传递到数据库中,并将相应的结果返回,供R环境使用。
  • 通过sqldf包,在R内部使用标准SQL对数据进行预处理,包括group by,order by,join,where等操作。
  • 当然R最重要的用途是将数据挖掘的结果转义为标准SQL语言,利用数据库来实现挖掘结果。当然有人说了,不是有pmml可以将模型嵌入到数据库么?!扯!到现在我也没见pmml成为应用标准,老老实实的将模型结果转义到SQL才是王道。比如用于概率预测的Logistic回归或者分类模型的Tree-based Models,这些模型的转义工作都不难,这样最终的工程实施都脱离了R环境,更具通用性,且利用了数据库的高速性能。

说句题外话:不知道哪位看官见过70万字符长度的庞大SQL语句——是的,你没看错,70w,R转义的,可以执行,对于数据库而言不过是半分钟的事情。

在2008年,第一届中国R语言会议上,来自于艾瑞咨询的张翔为大家展示了一组极具震撼力的泡泡图,而这段视频便是Hans Rosling 在2006年 TED 上的演讲,讲述的是1956年-2006年之间各国间的经济发展变化。虽然个人认为泡泡图的实现的技术并没有太多技术含量,但惊讶于Hans Rosling大智若愚的演讲能力,甚至再看完五六遍之后,仍然还会被其感染。

再后来, Hans Rosling将这款数据展示产品卖给了Google,而Google又将其整合到Visualisation API 里,于是我们可以调用Google的API来使用这些有趣的图形展示功能。

更令人鼓舞的是:前不久,Markus Gesmann, Castillo两位老大,写了googleVis这个R扩展包,将Hans Rosling的泡泡图同R语言完美的整合在了一起:)

下面就是在R语言里,通过googleVis包,对我国2006年至今的货币供应量(M0、M1、M2)Motion Chart(将横坐标设置为time以后,点击Play)

碎碎念:从泡泡的跳动上看,央行投放的货币供应量(亿元)在2009-2011年间,增加的最为剧烈,当然这个时间段也是房价近十年间增长最为迅速的两年。 而现在呢,房市开始限购,国内房地产市场已经撑不住超发的货币,因此广泛的通货膨胀开始蔓延。引用郎咸平的几段话:

  • 我们近几年发了76.5万亿的货币,是GDP的2.5倍,而美国的货币量除上GDP只是0.6,只有我们的四分之一
  • 唯一的解释只能是货币的购买力在下降,这就是经济学里的通货膨胀

最后,关于googleVis包的使用请移步至这里

三天前,在统计之都论坛上问到了如何做Matrix67 博客上的平滑马赛克图,我是好事之徒,颠颠地跑去瞧了一眼。

恩,蛮有意思的,而且非常黄,非常暴力!但比较悲剧的是我不会用 Mathematica,只好用 R 实现了一下。

本来标题党一些,叫做《一千二百个女人和我的故事》,想想还是算了吧,虽说是用了 1200 个漂亮女人组成了我的头像,但我一个也不认识,哈哈。

用的原图我就不贴了,实际上我是戴着眼镜的,马赛克平滑以后,不明显了。最后是代码,非常简单,不到 20 行:

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
setwd('D:/doc/image/me')
library(ReadImages)
library(sqldf)
me <- read.jpeg('fun.jpg')
meid <- data.frame(z = 1:1200, y = as.numeric(me))
meid <- sqldf('select * from meid order by y')

setwd('D:/doc/image/others')
tmp <- NULL
for(i in dir()) tmp[[i]] <- read.jpeg(i)
id <- sapply(tmp, mean)
id <- data.frame(n = names(id), m = id)
id <- sqldf('select * from id order by m')
idx <- cbind(id, meid)
idx <- sqldf('select * from idx order by z')

setwd('D:/doc/image')
png('me.png', height = 1000, width = 750)
par(mfcol = c(40,30), mar = rep(0,4), xpd = NA)
for(i in idx$n) plot(tmp[[i]])
dev.off()

## 处理图片 ##
setwd('D:/doc/image/others')
shell("convert *.jpg -crop 120x120+10+5 thumbnail%03d.png")
shell("del *.jpg")
shell("convert -type Grayscale *.png thumbnail%03d.png")

大概所需要的时间:构思写代码 1 个小时,下载和整理图片时间长点,3个多小时(当然你本地资源和 Matrix67 一样丰富的话另说,哈)。

重更新代码

以前的包不能用了,另外不太优雅

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
library(glue)
library(imager) # mac 需要 X11 的支持
library(jpeg)

## 处理图片 ##
setwd('/Users/liusizhe/bitbucket/code/imagemagick')
## system("convert source_gray.png -colorspace gray source_gray.jpg")
## system("convert source_gray.jpg -resize 50 source40.jpg")

size = 300
s <- 1:size
center <- sample(1:20, size = size, replace = TRUE)
## 基本思想是将照片等距加亮,这样就可以表达所有主体图片了
## -crop 250x250+5+10 的意思是(裁剪)
## 图片的中心店为准,向右 5 个像素点,向下 10 个像素点,取 250x250的图片出来
## 这里为了保证图片不会过于相似,所以做了随机“晃动”
for (i in s){
cmd <- glue('convert input1.jpg -set option:modulate:colorspace hsb -modulate {s[i]},100 -resize 320 -crop 250x250+{center[i]}+{center[i]} image/modulate{s[i]}.jpg')
system(cmd)
}
me <- readJPEG("source40.jpg")

## 将前面生成的图片读到一个 list
tmp <- list()
for(i in dir('/Users/liusizhe/bitbucket/code/imagemagick/image')){
tmp[[i]] <- load.image(paste('image/', i, sep = ''))
}

## 找到主图每个像素点同生成的等距加亮图片,距离最近的那个
id_mean <- sapply(tmp, mean)
image_order <- outer(as.numeric(me), id_mean, '-') |>
abs() |> apply(1, which.min)

## 绘制图形
## 这里要注意 height 和 width 必须和 me 的像素一致
jpeg('me.png', height = nrow(me)* 50, width = ncol(me)* 50)
par(mfcol = dim(me), mai = rep(0,4), xpd = FALSE)
for(i in image_order) plot(tmp[[i]], axes = FALSE)
dev.off()

## 删除 image 目录下的临时文件
system("rm image/*")

0%