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语法，可参见这里。

	import numpy as np
	x = np.random.randn(100)
	y = 2*x + np.random.randn(100)
	%load_ext rpy2.ipython
	%%R -i x,y -w 500 -h 300
	df <- data.frame(x,y)
	m <- lm(y~x)
	inter <- m$coef[1]
	slop <- m$coef[2]
	library(ggplot2)
	p <- ggplot(df, aes(x,y))
	p <- p +
	geom_point() +
	geom_abline(intercept = inter, slope = slop)
	print(p)

view raw rmagic.py hosted with ❤ by GitHub

理解MCMC

在贝叶斯统计中，经常需要计算后验概率，概率计算就涉及到积分问题。一种解决方法是用解析式得到后验概率直接计算，另一种是利用统计模拟来计算近似值。 考虑一个简单问题，我们对一个硬币反复投掷，对于出现正面的概率theta先主观设定为一个均匀分布，然后实际投掷14次，得到11次正面，要根据这个信息data来更新后验概率。下面用统计模拟来计算。

	myData=c(1,1,1,1,1,1,1,1,1,1,1,0,0,0)
	# 尝试1万个不同的参数
	tryn = 1e4
	Theta = sort(runif(tryn))
	pTheta = 1/tryn
	z = sum( myData==1 )
	N = length( myData )
	# 似然函数
	pDataGivenTheta = Theta^z * (1-Theta)^(N-z)
	pData = sum( pDataGivenTheta * pTheta )
	# 后验概率
	pThetaGivenData = pDataGivenTheta * pTheta / pData
	plot(x=Theta,y=pThetaGivenData,type='l')

view raw simulation1.R hosted with ❤ by GitHub

因为实际的数据表现是正面出现次数较多，所以后验概率会做出相应调整。这个统计模拟的例子非常简单，容易实施。但如果涉及的参数个数很多时，计算量就会非常大。此时就可以尝试另一种统计模拟方法，即MCMC(Markov chain Monte Carlo)。我们先考虑一个思想实验：

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