我试图对许多数据点进行高斯拟合。E、 我有一个256 x 262144的数据阵列。256点需要拟合高斯分布,我需要262144点。
有时高斯分布的峰值在数据范围之外,因此得到一个准确的平均结果曲线拟合是最好的方法。即使峰值在范围内,曲线拟合也会给出更好的西格玛,因为其他数据不在范围内。
我用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_做了一些修改,使之适合工作。
我用来让它工作的代码是
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发布的解决方案非常好,而且解决得非常快。但是,除非高斯的有效区域在适当的区域内,否则它似乎不起作用。不过,我得测试平均值是否仍然准确,因为这是我用它来做的主要事情。
最简单的事情就是把问题线性化。你用的是非线性迭代法,比线性最小二乘法要慢。
基本上,你有:
y = height * exp(-(x - mu)^2 / (2 * sigma
^2)要使其成为线性方程,取两边的(自然)对数:
然后将其简化为多项式:
我们可以用更简单的形式重铸:
其中:
然而,有一个陷阱。在分布的“尾部”存在噪声时,这将变得不稳定。
因此,我们只需要使用分布“峰值”附近的数据。在拟合中只包含低于某个阈值的数据是很容易的。在这个例子中,我只包括大于我们拟合的给定高斯曲线最大观测值20%的数据。
不过,一旦我们做到了,就相当快了。求解262144条不同的高斯曲线只需要大约1分钟(如果在这么大的东西上运行,请确保删除代码的绘图部分…)。它也很容易并行化,如果你想。。。
对于并行版本,我们只需要更改主函数。(我们还需要一个伪函数,因为
multiprocessing.Pool.imap
不能为它的函数提供额外的参数…)它看起来像这样:编辑:如果简单的多项式拟合效果不好,请尝试使用@tslisten共享的y值来加权问题,as mentioned in the link/paper(和Stefan van der Walt实现的问题,尽管我的实现有点不同)。
如果这仍然给你带来麻烦,那么尝试迭代地重新加权最小二乘问题(link@tslisten中提到的最后一个“最佳”推荐方法)。不过,请记住,这将相当缓慢。
相关问题 更多 >
编程相关推荐