时隔四年,无意中看到了双重高斯分布拟合
高斯分布在自然界非常常见,中心极限定理很好的说明了它,但事情往往不是那么纯粹,很多时候我们得到的结果里面会混入两个截然不同的样本数据集,它们虽然各自都是高斯分布,但是它们的均值和方差都不一样,如果拿到的是它们的混合数据,就不能简单的使用一个高斯拟合来处理它了。
如果我们有比较强的背景知识,或者看了如下分布的条形图,会下意识的猜想出是两个高斯分布的混合,但是想从数据的角度来探索,两个独立的高斯分布各自独立的均值和方差该如何推测出来呢?
The bimodal distribution was approximated by the 'normalmixEM’ function in the 'mixtools’ package of R.
生成测试数据
生成两个高斯分布,它们有各自独立的均值和方差
a=rnorm(1000)
b=rnorm(1000,mean = 1,sd=1.4)
par(mfrow=c(1,2))
hist(a);hist(b)
可视化如下:
使用mixtools包的normalmixEM函数
首先加载安装好的包:
> library(mixtools)
mixtools package, version 1.1.0, Released 2017-03-10
This package is based upon work supported by the National Science Foundation under Grant No. SES-0518772.
先看帮助文档的代码,
data(faithful)
attach(faithful)
set.seed(100)
system.time(out<-normalmixEM(waiting, arbvar = FALSE, epsilon = 1e-03))
out
hist(waiting)
plot(out,2)
可以看到,很简单一个函数,就可以把faithful这个数据框里面的waiting列的数据进行双重高斯分布拟合!
在我们的数据上面是使用
前面我们根据R包说明书进行了示例数据分析,那么理所当然就学会了它,然后就应该是应用于自己的数据,我们测试数据是两个高斯分布的向量,如果要使用mixtools包的normalmixEM函数处理它,就应该是需要把两个向量合并.
d=c(a,b)
out<-normalmixEM(d, arbvar = FALSE, epsilon = 1e-03)
hist(d);plot(out,2)
可以看到,其实这个函数并不是把我们生成的两个高斯分布完成拆解开来,前面我们模拟的是平均值分别是0和1的两个分布,但是函数拟合后是0和2的两个高斯分布,如下:
因为这个不是我的课题了,仅仅是想分享给生信技能树的粉丝,符合你们的需求,所以接下来靠你们自己花费时间去摸索哈,比如 http://exploringdatablog.blogspot.com/2011/08/fitting-mixture-distributions-with-r.html
后记
其实这个R包早在2009就发表了,不知道为什么我四年前没有搜索到它,那个时候还没有生信技能树公众号呢
Benaglia, T., Chauveau, D., Hunter, D. R., and Young, D. mixtools: An R package for analyzing finite mixture models. Journal of Statistical Software, 32(6):1-29, 2009.
最后,跟大家互动一下,说到双重高斯分布,大家首先想到的是什么呢?
写了万篇教程后,尤其是有些脍炙人口的,广为传播的神级推文,会感觉后续写新的,很难超越以前优秀的自己,大家怎么看这件事?