前往小程序,Get更优阅读体验!
立即前往
首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >时隔四年,无意中看到了双重高斯分布拟合

时隔四年,无意中看到了双重高斯分布拟合

作者头像
生信技能树
发布2019-09-18 17:03:03
2.3K0
发布2019-09-18 17:03:03
举报
文章被收录于专栏:生信技能树

高斯分布在自然界非常常见,中心极限定理很好的说明了它,但事情往往不是那么纯粹,很多时候我们得到的结果里面会混入两个截然不同的样本数据集,它们虽然各自都是高斯分布,但是它们的均值和方差都不一样,如果拿到的是它们的混合数据,就不能简单的使用一个高斯拟合来处理它了。

如果我们有比较强的背景知识,或者看了如下分布的条形图,会下意识的猜想出是两个高斯分布的混合,但是想从数据的角度来探索,两个独立的高斯分布各自独立的均值和方差该如何推测出来呢?

这个难题早在四年前我处理免疫组库数据就遇到了,那个时候功力很浅,全网搜索各自资料也没有解决,后来换工作了,这个项目也就不了了之,最近看文献比较多,其中一个文章的描述了一个R包可以做:

The bimodal distribution was approximated by the ‘normalmixEM’ function in the ‘mixtools’ package of R.

生成测试数据

生成两个高斯分布,它们有各自独立的均值和方差

代码语言:javascript
复制
a=rnorm(1000)
b=rnorm(1000,mean = 1,sd=1.4)
par(mfrow=c(1,2))
hist(a);hist(b)

可视化如下:

使用mixtools包的normalmixEM函数

首先加载安装好的包:

代码语言:javascript
复制
> 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.

先看帮助文档的代码,

代码语言:javascript
复制
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函数处理它,就应该是需要把两个向量合并.

代码语言:javascript
复制
d=c(a,b)
out<-normalmixEM(d, arbvar = FALSE, epsilon = 1e-03)
hist(d);plot(out,2)

可以看到,其实这个函数并不是把我们生成的两个高斯分布完成拆解开来,前面我们模拟的是平均值分别是0和1的两个分布,但是函数拟合后是0和2的两个高斯分布,如下:

这个时候,有两种解决方案,首先你可以继续花时间去攻克这个R包的详细文档,看看有没有什么参数可以调整,其次你可以去看看其它R包或者算法。

因为这个不是我的课题了,仅仅是想分享给生信技能树的粉丝,符合你们的需求,所以接下来靠你们自己花费时间去摸索哈,比如 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.

最后,跟大家互动一下,说到双重高斯分布,大家首先想到的是什么呢?

最后,题外话,唠叨一句话

写了万篇教程后,尤其是有些脍炙人口的,广为传播的神级推文,会感觉后续写新的,很难超越以前优秀的自己,大家怎么看这件事?

本文参与 腾讯云自媒体同步曝光计划,分享自微信公众号。
原始发表:2019-09-16,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 生信技能树 微信公众号,前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体同步曝光计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
目录
  • 生成测试数据
  • 使用mixtools包的normalmixEM函数
  • 在我们的数据上面是使用
  • 后记
领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档