星期日, 七月 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)。
就这样,旅行者每天不断的在不同岛间移动,最终他在各岛上呆的时间之比,就会收敛为各岛人口之比,这样就通过一个随机模拟得到了一个概率分布。

星期日, 六月 22, 2014

一个简单排队论问题的python实现

一个诊所只有一个医生,病人到来的时间是随机的,从早上九点开始,服从一个时间参数为10min的泊松过程,即每个人到来的时间服从独立同分布的指数分布,其期望为10min,每个病人到来之后,下一个病人到来的时间服从独立同分布的指数分布,期望为10min。当一个病人到来以后,将等待直到医生有空。每个医生在每个病人上花费的时间是一个随机变量,在5min到10min之间均匀分布。诊所从下午4点不再接受新病人,最后一个病人走后,诊所关门。
(1)模拟一天的病人到来和医生接诊情况,有多少个病人来诊所看病?其中有多少病人需要等待医生?平均等待时间是多少?
(2)模拟100天的情况,给出上面各个数值的分布

星期日, 二月 16, 2014

从模拟角度理解混合模型第三:广义和加性混合模型(完)

之前讨论的混合效应模型均是最常见的线性模型基础上的扩展。但在实际使用中会发现用到其它情况,例如要处理分类问题,需要logistic回归模型。再比如说线性关系不存在,需要用加性模型来处理复杂的非线性关系。这两种模型分别归于广义线性模型和加性模型。那么在这两种模型基础上再考虑解释变量的混合效应,要分别衍生出广义混合效应模型和加性混合效应模型。在R中lme4包可以处理广义混合效应模型,而mgcv包可以处理加性混合效应模型。在之前文章一样,我们先模拟这两种数据,再分别进行建模。代码来源于《Regress by simulation》

星期三, 一月 29, 2014

用模拟来理解混合效应模型之二:Random Intercept and slope model

在之前的这篇文章中,混合效应模型的意义已经说的比较清楚了,简言之,样本中不能穷尽总体level的变量都是随机效应。也可以这么认为,会影响目标变量,但我们不关心的解释变量都是随机效应。之前文章的随机效应只影响模型的intercept,那么也会有影响slope的随机效应。我们先来看一下这种混合效应模型的假设,再用假设来生成数据,并建模和绘图。

Yij = b0 + (b1+si)*Xij +  bi + eij 

* b0: fixed intercept
* b1: fixed slope
* X: fixed effect
* bi: random effect(influence intercept)
* eij: noise
* si: random effect(influence slope)