星期日, 七月 20, 2014

python和ggplot2

python有个非常强大的工具,那就是ipython notebook。用户可以在浏览器中直接编写python脚本,并立即得到输出结果。这类文档可以存为ipynb分享给其它人,也可以存为html直接放在网站上,非常有利于学习交流。

在R语言方面就缺乏这类工具,不过ipython有一种“魔法”,可以在ipython中运行其它语言。在数据分析时,可以将python和R代码混编,充分利用两种语言的优势。以可视化为例,R的ggplot2图形语法可谓是独步江湖,python中虽然已经有不少优秀的绘图库。但总不及ggplot2用得习惯。下面的小例子就是示范在ipython notebook中画ggplot2。

首先是用numpy库建立两个向量,再用%load_ext建立python和R的连接机制。之后在ipython notebook的一个cell里面就可以使用%R后面接R代码行,或者使用%%R使用代码块。如果要在R代码块中读入python的对象,需要使用-i参数。

其它例子可以参见这个,python中也有人复制了ggplot语法,可参见这里

星期日, 七月 06, 2014

理解MCMC

在贝叶斯统计中,经常需要计算后验概率,概率计算就涉及到积分问题。一种解决方法是用解析式得到后验概率直接计算,另一种是利用统计模拟来计算近似值。 考虑一个简单问题,我们对一个硬币反复投掷,对于出现正面的概率theta先主观设定为一个均匀分布,然后实际投掷14次,得到11次正面,要根据这个信息data来更新后验概率。下面用统计模拟来计算。 因为实际的数据表现是正面出现次数较多,所以后验概率会做出相应调整。这个统计模拟的例子非常简单,容易实施。但如果涉及的参数个数很多时,计算量就会非常大。此时就可以尝试另一种统计模拟方法,即MCMC(Markov chain Monte Carlo)。我们先考虑一个思想实验:

一个岛国下面有7个岛,各岛人数不一样。旅行者已经在其中一个岛上,他要去各岛游历,准备在人口多的岛上游玩的时间长一些。但他并不知道具体的人口数,可以询问本岛和相邻岛的市长。他的旅行计划如下:
首先每天会扔一枚硬币,如果正面向上,就计划去左岛,如果反面向上,就计划去右岛。
然后比较当前所在岛人口数a和计划去的目的岛人口数b,若b小于a,则实际去的概率为b/a,若b大于a,则实际去的概率为1。所以实际去的概率归纳为min(1,b/a)。
就这样,旅行者每天不断的在不同岛间移动,最终他在各岛上呆的时间之比,就会收敛为各岛人口之比,这样就通过一个随机模拟得到了一个概率分布。