我正在尝试对许多数据点进行高斯拟合。例如,我有一个 256 x 262144 的数据数组。需要将 256 个点拟合到高斯分布,我需要其中的 262144 个。
有时高斯分布的峰值超出数据范围,因此要获得准确的平均结果曲线拟合是最好的方法。即使峰值在范围内,曲线拟合也会给出更好的 sigma,因为其他数据不在范围内。
我使用http://www.scipy.org/Cookbook/FittingData中的代码为一个数据点工作。
我试图重复这个算法,但看起来它需要大约 43 分钟才能解决这个问题。是否有一种已经编写好的快速方法可以并行或更有效地执行此操作?
from scipy import optimize
from numpy import *
import numpy
# Fitting code taken from: http://www.scipy.org/Cookbook/FittingData
class Parameter:
def __init__(self, value):
self.value = value
def set(self, value):
self.value = value
def __call__(self):
return self.value
def fit(function, parameters, y, x = None):
def f(params):
i = 0
for p in parameters:
p.set(params[i])
i += 1
return y - function(x)
if x is None: x = arange(y.shape[0])
p = [param() for param in parameters]
optimize.leastsq(f, p)
def nd_fit(function, parameters, y, x = None, axis=0):
"""
Tries to an n-dimensional array to the data as though each point is a new dataset valid across the appropriate axis.
"""
y = y.swapaxes(0, axis)
shape = y.shape
axis_of_interest_len = shape[0]
prod = numpy.array(shape[1:]).prod()
y = y.reshape(axis_of_interest_len, prod)
params = numpy.zeros([len(parameters), prod])
for i in range(prod):
print "at %d of %d"%(i, prod)
fit(function, parameters, y[:,i], x)
for p in range(len(parameters)):
params[p, i] = parameters[p]()
shape[0] = len(parameters)
params = params.reshape(shape)
return params
请注意,数据不一定是 256x262144,我已经在 nd_fit 中做了一些捏造来完成这项工作。
我用来让它工作的代码是
from curve_fitting import *
import numpy
frames = numpy.load("data.npy")
y = frames[:,0,0,20,40]
x = range(0, 512, 2)
mu = Parameter(x[argmax(y)])
height = Parameter(max(y))
sigma = Parameter(50)
def f(x): return height() * exp (-((x - mu()) / sigma()) ** 2)
ls_data = nd_fit(f, [mu, sigma, height], frames, x, 0)
注意:@JoeKington 下面发布的解决方案很棒,而且解决得非常快。然而,除非高斯的重要区域在适当的区域内,否则它似乎不起作用。我将不得不测试平均值是否仍然准确,因为这是我使用它的主要目的。