在Python中分离高斯混合物

12

这里有一个物理实验的结果,可以用直方图[i, amount_of(i)]表示。我认为这个结果可以用4-6个高斯函数的混合来估计。

是否有Python包能够将直方图作为输入,并返回混合分布中每个高斯分布的均值和方差?

例如原始数据如下:

Sample data


3
顺带一提,这是高斯混合模型(mixture of gaussians)而不是高斯求和(sum of gaussians)(多个独立的高斯函数相加也会得到正态分布)。你可能想要使用PyMix库(虽然我个人没有使用过),可在此处了解更多:http://www.pymix.org/pymix/index.php?n=PyMix.Tutorial。 - David Robinson
由于实验的物理感受 - 这应该是真正的总和,而不是混合物。此外,数学的最终目标是找到每个“子人口”(高斯曲线下的区域)在整个“人口”(曲线下的区域)中的百分比 - 据我所知,混合模型无法回答这个问题。 - Belegnar
2
当然可以,这就是它们的用途(或者更确切地说,它们可以估计-当然,由于涉及到随机事件,所以没有一种明确回答问题的方式)。除非我弄错了,我认为你的意思是一个混合物(除了分布的混合物类似于直方图的“总和”,其中一个直方图放在另一个上面)。 - David Robinson
在“每个点”指的是“Y轴上每个+1值”。例如,如果X接近56,则有25K个“点”。大多数这些点属于第一个高斯分布,少量属于第二个高斯分布,而更少的属于第三和第四个高斯分布。 - Belegnar
1
每个点应该属于一个且仅属于一个高斯分布——这正是混合模型的定义(请参见下面我的回答——它对您有用吗?)。您可能在考虑混合成员资格模型,其中每个点可以同时属于多个类别。 - David Robinson
显示剩余3条评论
1个回答

18

这是一个高斯混合模型,可以使用最大期望算法来估计(基本上,它同时找到分布的中心和均值以及它们如何混合)。

这在PyMix包中实现。下面我生成一个正态分布的混合示例,并使用PyMix对其进行混合模型拟合,包括确定您感兴趣的子群体的大小:

# requires numpy and PyMix (matplotlib is just for making a histogram)
import random
import numpy as np
from matplotlib import pyplot as plt
import mixture

random.seed(010713)  # to make it reproducible

# create a mixture of normals:
#  1000 from N(0, 1)
#  2000 from N(6, 2)
mix = np.concatenate([np.random.normal(0, 1, [1000]),
                      np.random.normal(6, 2, [2000])])

# histogram:
plt.hist(mix, bins=20)
plt.savefig("mixture.pdf")

以上代码只是生成和绘制混合结果,结果如下图所示:

enter image description here

现在要使用PyMix来确定百分比:

data = mixture.DataSet()
data.fromArray(mix)

# start them off with something arbitrary (probably based on a guess from the figure)
n1 = mixture.NormalDistribution(-1,1)
n2 = mixture.NormalDistribution(1,1)
m = mixture.MixtureModel(2,[0.5,0.5], [n1,n2])

# perform expectation maximization
m.EM(data, 40, .1)
print m

这个的输出模型是:

G = 2
p = 1
pi =[ 0.33307859  0.66692141]
compFix = [0, 0]
Component 0:
  ProductDist: 
  Normal:  [0.0360178848449, 1.03018725918]

Component 1:
  ProductDist: 
  Normal:  [5.86848468319, 2.0158608802]

注意它正确地发现了两个正态分布 (一个 N(0, 1) 和一个 N(6, 2), 大约), 它还估计了pi,这是两个分布中的比例 (你在评论中提到这是你最感兴趣的)。我们在第一个分布中有1000个样本,在第二个分布中有2000个样本,它几乎完全正确地得出了比例: [0.33307859 0.66692141]。如果您想直接获得此值,请运行 m.pi

一些注意事项:

  • 该方法需要一个值向量,而不是直方图。将数据转换为1D向量应该很容易 (即,将 [(1.4, 2), (2.6, 3)] 转换成 [1.4, 1.4, 2.6, 2.6, 2.6])
  • 我们必须提前猜测高斯分布的数量(如果您要求混合物为2,则它不会找到混合物的4)。
  • 我们必须对分布进行一些初始估计。 如果您做出合理的猜测,它应该会收敛到正确的估计。

非常感谢!抱歉,我感觉自己像个傻瓜 - 因为我的显示器太小了,直到现在才看到你的回答。 - Belegnar

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接