拟合部分高斯函数

Question

我正在尝试使用scikit-learn拟合高斯函数的总和，因为与使用curve_fit相比，scikit-learn GaussianMixture似乎更加健壮。 问题：它不能很好地拟合甚至是单个高斯峰截断部分的情况：

from sklearn import mixture
import matplotlib.pyplot
import matplotlib.mlab
import numpy as np

clf = mixture.GaussianMixture(n_components=1, covariance_type='full')
data = np.random.randn(10000)
data = [[x] for x in data]
clf.fit(data)
data = [item for sublist in data for item in sublist]
rangeMin = int(np.floor(np.min(data)))
rangeMax = int(np.ceil(np.max(data)))
h = matplotlib.pyplot.hist(data, range=(rangeMin, rangeMax), normed=True);
plt.plot(np.linspace(rangeMin, rangeMax),
         mlab.normpdf(np.linspace(rangeMin, rangeMax),
                      clf.means_, np.sqrt(clf.covariances_[0]))[0])

给出，现在将data = [[x] for x in data]更改为data = [[x] for x in data if x <0]以截断分布返回。有什么想法可以使截断适当地拟合吗？

注意：该分布不一定在中间被截断，可能留下完整分布的50％至100％之间的任何内容。

如果有人能指向替代软件包，我也会很高兴。我只尝试过curve_fit，但在涉及两个以上峰值时无法得到有用的结果。

- lhcgeneva

---

你期望什么？高斯分布的一半并不是高斯分布，因此一个简单的高斯拟合无法产生良好的结果。如果你正在尝试高斯混合模型，可以使用混合物。但对于多个高斯分布的混合物，你需要超过n_components=1。这就是整个想法。 - sascha

不确定这是正确的。如果你使用curve_fit对单个高斯函数进行拟合，它会非常好地工作。只是对于多个峰值，这种方法并不稳健。另外，为了简单起见，n_components=1... - lhcgeneva

你对“curve_fit”方法有什么具体的困扰？你尝试过scipy.optimize中的其他最小化方法吗？它是一个非常强大的工具箱，我认为不应该轻易地被忽视。 - Vlas Sokolov

这个有进展了吗？ - O.rka

很遗憾，目前还没有适用于多个峰值的强大解决方案。 - lhcgeneva

Answer 1

一个有点粗糙但简单的解决方案是将曲线分成两半(data = [[x] for x in data if x < 0])，将左半部分镜像(data.append([-data[d][0]]))，然后进行常规高斯拟合。

import numpy as np
from sklearn import mixture
import matplotlib.pyplot as plt
import matplotlib.mlab as mlab

np.random.seed(seed=42)
n = 10000

clf = mixture.GaussianMixture(n_components=1, covariance_type='full')

#split the data and mirror it
data = np.random.randn(n)
data = [[x] for x in data if x < 0]
n = len(data)
for d in range(n):
    data.append([-data[d][0]])

clf.fit(data)
data = [item for sublist in data for item in sublist]
rangeMin = int(np.floor(np.min(data)))
rangeMax = int(np.ceil(np.max(data)))
h = plt.hist(data[0:n], bins=20, range=(rangeMin, rangeMax), normed=True);
plt.plot(np.linspace(rangeMin, rangeMax),
         mlab.normpdf(np.linspace(rangeMin, rangeMax),
                      clf.means_, np.sqrt(clf.covariances_[0]))[0] * 2)

plt.show()

- Maximilian Peters

---

这看起来很不错，但我应该提到的是它不一定是分布的一半，而可以是任意截断（通常留下比分布多一点的部分以适应）。对此感到抱歉，我已经修改了问题。 - lhcgeneva

Answer 2

lhcgeneva，问题在于一旦你有了不包括曲线最大值的数据，越来越多的高斯函数就可以拟合：

黑点代表数据，红点代表拟合的高斯函数

在这张图中，黑点代表需要拟合曲线的数据，红点代表拟合结果。使用一个简单的高斯函数拟合算法得到了这个结果。

- Ed Breen