Как быстро выполнить аппроксимацию методом наименьших квадратов для нескольких наборов данных?
Я пытаюсь выполнить гауссовское соответствие по многим точкам данных. Например. У меня есть массив данных 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_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, великолепно и решает очень быстро. Однако это не работает, если значительная область гаусса не находится внутри соответствующей области. Мне нужно будет проверить, является ли среднее значение точным, поскольку это главное, для чего я это использую.
