Python-Fitting 二维高斯到数据集

Question

我在 .txt 文件中有数据点(分隔符 = 空格)，第一列是 x 轴，第二列是 y 轴。我想使用 Python 将 2D 高斯拟合到这些数据点。事实是，我不了解高斯拟合(一维或二维)背后的理论。我在 stackoverflow 上阅读过类似的帖子并获得了代码，但它不太合适。请有人帮忙。谢谢

下面是我在 .txt 文件中的内容:

3.369016418457e+02 3.761813938618e-01
3.369006652832e+02 4.078308343887e-01
3.368996887207e+02 4.220226705074e-01
3.368987121582e+02 4.200653433800e-01
3.368977355957e+02 4.454285204411e-01
3.368967590332e+02 4.156131148338e-01
3.368957824707e+02 3.989491164684e-01
3.368948059082e+02 4.512043893337e-01
3.368938293457e+02 4.565380811691e-01
3.368928527832e+02 4.095999598503e-01
3.368918762207e+02 4.196371734142e-01
3.368908996582e+02 4.002234041691e-01
3.368899230957e+02 4.133881926537e-01
3.368889465332e+02 4.394644796848e-01
3.368879699707e+02 4.504477381706e-01
3.368869934082e+02 3.946847021580e-01
3.368860168457e+02 4.214486181736e-01
3.368850402832e+02 3.753573596478e-01
3.368840637207e+02 3.673824667931e-01
3.368830871582e+02 4.088735878468e-01
3.368821105957e+02 4.351278841496e-01
3.368811340332e+02 4.393630325794e-01
3.368801574707e+02 4.210205972195e-01
3.368791809082e+02 4.322172403336e-01
3.368782043457e+02 4.652716219425e-01
3.368772277832e+02 5.251595377922e-01
3.368762512207e+02 5.873318314552e-01
3.368752746582e+02 6.823697686195e-01
3.368742980957e+02 8.375824093819e-01
3.368733215332e+02 9.335057735443e-01
3.368723449707e+02 1.083636641502e+00
3.368713684082e+02 1.170072913170e+00
3.368703918457e+02 1.224770784378e+00
3.368694152832e+02 1.158735513687e+00
3.368684387207e+02 1.131350398064e+00
3.368674621582e+02 1.073648810387e+00
3.368664855957e+02 9.659162163734e-01
3.368655090332e+02 8.495713472366e-01
3.368645324707e+02 7.637447714806e-01
3.368635559082e+02 6.956064105034e-01
3.368625793457e+02 6.713893413544e-01
3.368616027832e+02 5.285132527351e-01
3.368606262207e+02 4.968771338463e-01
3.368596496582e+02 5.077748298645e-01
3.368586730957e+02 4.686309695244e-01
3.368576965332e+02 4.693206846714e-01
3.368567199707e+02 4.462305009365e-01
3.368557434082e+02 3.872672021389e-01
3.368547668457e+02 4.243377447128e-01
3.368537902832e+02 3.918920457363e-01
3.368528137207e+02 3.848327994347e-01
3.368518371582e+02 4.093343317509e-01
3.368508605957e+02 4.321203231812e-01

下面是我试过的代码:

%pylab inline
import matplotlib.pyplot as plt
import numpy as np
import astropy
import scipy.optimize as opt
import pylab as plb
from scipy.optimize import curve_fit
from scipy import asarray as ar,exp
x,y=np.loadtxt('taper2reduced.txt', unpack= True, delimiter=' ')

mean = sum(x * y) / sum(y)
sigma = np.sqrt(sum(y * (x - mean)**2) / sum(y))

def Gauss(x, a, x0, sigma):
    <pre>return a * np.exp(-(x - x0)\**2 / (2 * sigma**2))<code>

popt,pcov = curve_fit(Gauss, x, y, p0=[max(y), mean, sigma])

plt.plot(x, y, 'b+:', label='data')
plt.plot(x, Gauss(x, *popt), 'r-', label='fit')
plt.legend()
plt.title('Fig. 1 - Fit for Frequency')
plt.xlabel('Frequecy (GHz)')
plt.ylabel('Flux Density (mJy)')
plt.show()

Answer 1

您的问题是您的函数不能很好地反射(reflect)您的数据集。您定义了 0 和 max_y 之间的分布，而实际上您的数据在 min_y 和 max_y 之间。像这样更改您的功能:

import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit
#function declaration with additional offset parameter
def Gauss(x, a, x0, sigma, offset):
    return a * np.exp(-(x - x0)**2 / (2 * sigma**2)) + offset
#loading x, y dataset
x, y = np.loadtxt('test.txt', unpack = True, delimiter=' ')
#calculate parameter of fit function with scipy's curve fitting algorithm
popt, pcov = curve_fit(Gauss, x, y, p0=[np.max(y), np.median(x), np.std(x), np.min(y)])
#plot original data
plt.plot(x, y, 'b+:', label='data')
#create different x value array for smooth fit function curve
x_fit = np.linspace(np.min(x), np.max(x), 1000)
#plot fit function
plt.plot(x_fit, Gauss(x_fit, *popt), 'r-', label='fit')
#beautify graph
plt.legend()
plt.title('Fig. 1 - Fit for Frequency')
plt.xlabel('Frequecy (GHz)')
plt.ylabel('Flux Density (mJy)')
plt.show()

输出:

你可能已经注意到我又改变了两件事。
我整理了进口商品。 It is not a good idea to load a lot of different unused functions and modules into your name space.
我改变了开始参数估计。我们不必在这里是正确的，通常可以使用近似值。不需要太多代码而且速度很快的东西。