使用Scipy进行高斯和曲线拟合

Question

我是一名生物信息学研究员，我们在 mRNA 上映射小 RNA。我们有每个 mRNA 上蛋白质的映射坐标，并计算蛋白质结合 mRNA 的位置与被小 RNA 结合的位点之间的相对距离。

我获得了以下数据集：

dist    eff
-69 3
-68 2
-67 1
-66 1
-60 1
-59 1
-58 1
-57 2
-56 1
-55 1
-54 1
-52 1
-50 2
-48 3
-47 1
-46 3
-45 1
-43 1
0   1
1   2
2   12
3   18
4   18
5   13
6   9
7   7
8   5
9   3
10  1
13  2
14  3
15  2
16  2
17  2
18  2
19  2
20  2
21  3
22  1
24  1
25  1
26  1
28  2
31  1
38  1
40  2

当我绘制数据时，我得到了三张图片：一张在3-4左右，另一张在20左右，最后一张在-50左右。

我尝试了三次样条插值，但是对我的数据效果不好。

我的想法是使用高斯曲线的总和进行曲线拟合。例如，在我的情况下，估计点5、20和-50处的三个高斯曲线。

如何做到这一点？

我看了一下scipy.optimize.curve_fit()，但是如何在精确间隔上拟合曲线呢？ 如何添加曲线以得到单一曲线？

- user1845261

---

K-means聚类对你是否适用？ - Geoff

虽然Jaime的好答案对你的数据不太适合，但是可以参考PyMix中的separate-mixture-of-gaussians-in-python。 - denis

Answer 1

import numpy as np
import matplotlib.pyplot as plt
import scipy.stats
import scipy.optimize

data = np.array([-69,3, -68, 2, -67, 1, -66, 1, -60, 1, -59, 1,
                 -58, 1, -57, 2, -56, 1, -55, 1, -54, 1, -52, 1,
                 -50, 2, -48, 3, -47, 1, -46, 3, -45, 1, -43, 1,
                 0, 1, 1, 2, 2, 12, 3, 18, 4, 18, 5, 13, 6, 9,
                 7, 7, 8, 5, 9, 3, 10, 1, 13, 2, 14, 3, 15, 2,
                 16, 2, 17, 2, 18, 2, 19, 2, 20, 2, 21, 3, 22, 1,
                 24, 1, 25, 1, 26, 1, 28, 2, 31, 1, 38, 1, 40, 2])
x, y = data.reshape(-1, 2).T

def tri_norm(x, *args):
    m1, m2, m3, s1, s2, s3, k1, k2, k3 = args
    ret = k1*scipy.stats.norm.pdf(x, loc=m1 ,scale=s1)
    ret += k2*scipy.stats.norm.pdf(x, loc=m2 ,scale=s2)
    ret += k3*scipy.stats.norm.pdf(x, loc=m3 ,scale=s3)
    return ret


params = [-50, 3, 20, 1, 1, 1, 1, 1, 1]

fitted_params,_ = scipy.optimize.curve_fit(tri_norm,x, y, p0=params)

plt.plot(x, y, 'o')
xx = np.linspace(np.min(x), np.max(x), 1000)
plt.plot(xx, tri_norm(xx, *fitted_params))
plt.show()

>>> fitted_params
array([ -60.46845528,    3.801281  ,   13.66342073,   28.26485602,
          1.63256981,   10.31905367,  110.51392765,   69.11867159,
         63.2545624 ])

因此，您可以看到您对三峰函数的想法与您的实际数据并不完全一致。

- Jaime

---

这实际上并不是一个很糟糕的匹配。 - dor00012