多くのデータポイントにガウスフィットを適用しようとしています。たとえば、256 x 262144 のデータ配列があります。256 ポイントをガウス分布に適合させる必要があり、そのうち 262144 個が必要です。
ガウス分布のピークがデータ範囲外にある場合があるため、正確な平均結果を得るには、曲線近似が最適な方法です。ピークが範囲内にある場合でも、他のデータが範囲内にないため、カーブ フィッティングによりシグマが改善されます。
http://www.scipy.org/Cookbook/FittingDataのコードを使用して、1 つのデータ ポイントに対してこれを機能させています。
このアルゴリズムを繰り返してみましたが、これを解決するには 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 であるとは限らないことに注意してください。
これを機能させるために使用するコードは
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 によって以下に投稿されたソリューションは素晴らしく、非常に高速に解決します。ただし、ガウスの重要な領域が適切な領域内にない限り、機能しないようです。ただし、これを使用する主な目的であるため、平均がまだ正確であるかどうかをテストする必要があります。
