如何使用Python分离两个高斯曲线？

Question

我测量了成千上万个粒子的荧光强度并制作了直方图，其中显示了两个相邻的高斯曲线。如何使用Python或其包将它们分成两个高斯曲线并制作两个新图？谢谢。

- Steve Xu

Answer 1

基本上，您需要推断高斯混合的参数。我将为演示生成类似的数据集。

使用已知参数生成混合物

from itertools import starmap

import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
from matplotlib import mlab
sns.set(color_codes=True)
# inline plots in jupyter notebook
%matplotlib inline


# generate synthetic data from a mixture of two Gaussians with equal weights
# the solution below readily generalises to more components 
nsamples = 10000
means = [30, 120]
sds = [10, 50]
weights = [0.5, 0.5]
draws = np.random.multinomial(nsamples, weights)
samples = np.concatenate(
    list(starmap(np.random.normal, zip(means, sds, draws)))
)

绘制分布

sns.distplot(samples)

推断参数

from sklearn.mixture import GaussianMixture

mixture = GaussianMixture(n_components=2).fit(samples.reshape(-1, 1))
means_hat = mixture.means_.flatten()
weights_hat = mixture.weights_.flatten()
sds_hat = np.sqrt(mixture.covariances_).flatten()

print(mixture.converged_)
print(means_hat)
print(sds_hat)
print(weights_hat)

我们得到：

True
[ 122.57524745   29.97741112]
[ 48.18013893  10.44561398]
[ 0.48559771  0.51440229]

你可以调整GaussianMixture的超参数来改善拟合，但现在看起来已经足够好了。现在我们可以绘制每个组件（我只绘制了第一个）：

mu1_h, sd1_h = means_hat[0], sds_hat[0]
x_axis = np.linspace(mu1_h-3*sd1_h, mu1_h+3*sd1_h, 1000)
plt.plot(x_axis, mlab.normpdf(x_axis, mu1_h, sd1_h))

P.S.

顺便说一句，看起来你正在处理受限数据，并且你的观察结果非常接近于左约束（零）。虽然高斯分布可能足以近似你的数据，但你应该小心谨慎，因为高斯分布假定无约束几何。

- Eli Korvigo

---

@Eli 是否有一种方法可以找到数据比例确实遵循正态分布的概率？ - Riley

@Riley 这取决于你对“比例”的定义。你能进一步阐述吗？通常，您会使用Shapiro-Wilk测试来估计从单变量高斯分布中观察到随机经验偏差的概率。对于混合和多元高斯分布，情况要复杂得多。 - Eli Korvigo

@EliKorvigo，假设您构建了2个F分布并像一开始那样将它们混合在一起。然后，您可以应用相同的方法并获得高斯混合，但是分布本身不是正态的。 我的问题是：是否有一种方法可以查看正态分布/结果是否适用于分布的“相关”部分？ 我知道您可以使用Shapiro进行正态性测试，但是然后我们需要分布的“截止点”并单独测试它们，因此我想知道这个软件包是否可以提供更简便和/或更短的方式，类似于某种p值。 - Riley

@Riley 只是想澄清一下，我在这里使用的软件包没有提供此功能。如果您在处理此问题时遇到困难，我猜您应该在交叉验证上提出问题。 - Eli Korvigo

对于2021年的读者来说，可能有点晚了，但是需要注意的是，函数normpdf已经从mlab中删除。可以使用scipy.stats中的norm.pdf。 - eidal