Python中高斯的单独混合答案

【问题标题】：Separate mixture of gaussians in PythonPython中高斯的单独混合
【发布时间】：2012-12-20 19:51:27
【问题描述】：

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

Python中是否有一个包，它以直方图作为输入，并返回混合分布中每个高斯分布的均值和方差？

原始数据，例如：

【问题讨论】：

顺便说一下，这是一个mixture of gaussians，而不是高斯的总和（多个独立高斯的总和也是正常的）。您可能想使用PyMix 库（尽管我个人没有使用过）。
由于实验的物理意义 - 这应该是真正的总和，而不是混合。此外，数学的最终目标是找出每个“子种群”（高斯下的面积）在整个“种群”（曲线下的面积）中的百分比 - 据我了解，混合模型无法回答这个问题。
当然可以，这就是他们的目的（或者更确切地说，他们可以估计——当然没有办法明确地回答这个问题，因为涉及到随机机会）。除非我弄错了，否则我认为您的意思是混合（除非混合分布就像它们的直方图的“总和”，一个叠加在另一个之上）。
除非-这些是相互依赖的高斯的总和吗？
each point should belong to one and only one gaussian- 这正是混合模型的含义（请参阅下面的答案-它对您有用吗？）。您正在考虑一个混合成员资格模型，其中每个点可以同时属于多个类别。

标签： python statistics normal-distribution

【解决方案1】：

这是一个mixture of gaussians，可以使用expectation maximization 方法进行估计（基本上，它在估计它们如何混合在一起的同时找到分布的中心和均值）。

这是在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")

以上代码所做的就是生成并绘制混合物。它看起来像这样：

现在实际使用 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，这是两个分布中的每一个中的分数（您在 cmets 中提到这是您最感兴趣的）。我们在第一个分布中有 1000 个，在第二个分布中有 2000 个，它几乎完全正确地进行了划分：[ 0.33307859 0.66692141]。如果您想直接获取此值，请执行m.pi。

几点说明：

此方法采用值向量，而不是直方图。将数据转换为一维向量应该很容易（即将[(1.4, 2), (2.6, 3)] 转换为[1.4, 1.4, 2.6, 2.6, 2.6]）
我们必须提前猜测高斯分布的数量（如果您要求混合 2，它不会计算出 4 的混合）。
我们必须对分布进行一些初步估计。如果您做出了稍微合理的猜测，它应该会收敛到正确的估计值。

【讨论】：

非常感谢！抱歉，我觉得自己像个傻瓜——因为显示器太小，我现在才看到你的答案。