11

我正在尝试对许多数据点进行高斯拟合。例如,我有一个 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 下面发布的解决方案很棒,而且解决得非常快。然而,除非高斯的重要区域在适当的区域内,否则它似乎不起作用。我将不得不测试平均值是否仍然准确,因为这是我使用它的主要目的。 高斯分布估计分析

4

1 回答 1

18

最简单的方法是将问题线性化。您正在使用非线性迭代方法,该方法将比线性最小二乘解决方案慢。

基本上,你有:

y = height * exp(-(x - mu)^2 / (2 * sigma^2))

为了使它成为一个线性方程,取两边的(自然)对数:

ln(y) = ln(height) - (x - mu)^2 / (2 * sigma^2)

然后将其简化为多项式:

ln(y) = -x^2 / (2 * sigma^2) + x * mu / sigma^2 - mu^2 / sigma^2 + ln(height)

我们可以用更简单的形式重新定义它:

ln(y) = A * x^2 + B * x + C

在哪里:

A = 1 / (2 * sigma^2)
B = mu / (2 * sigma^2)
C = mu^2 / sigma^2 + ln(height)

但是,有一个问题。如果分布的“尾部”中存在噪声,这将变得不稳定。

因此,我们只需要使用分布“峰值”附近的数据。在拟合中仅包含高于某个阈值的数据很容易。在这个例子中,我只包括大于我们正在拟合的给定高斯曲线的最大观察值的 20% 的数据。

但是,一旦我们这样做了,它就会相当快。求解 262144 条不同的高斯曲线只需要大约 1 分钟(如果你在这么大的东西上运行它,请务必删除代码的绘图部分......)。如果您愿意,并行化也很容易...

import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import itertools

def main():
    x, data = generate_data(256, 6)
    model = [invert(x, y) for y in data.T]
    sigma, mu, height = [np.array(item) for item in zip(*model)]
    prediction = gaussian(x, sigma, mu, height)

    plot(x, data, linestyle='none', marker='o')
    plot(x, prediction, linestyle='-')
    plt.show()

def invert(x, y):
    # Use only data within the "peak" (20% of the max value...)
    key_points = y > (0.2 * y.max())
    x = x[key_points]
    y = y[key_points]

    # Fit a 2nd order polynomial to the log of the observed values
    A, B, C = np.polyfit(x, np.log(y), 2)

    # Solve for the desired parameters...
    sigma = np.sqrt(-1 / (2.0 * A))
    mu = B * sigma**2
    height = np.exp(C + 0.5 * mu**2 / sigma**2)
    return sigma, mu, height

def generate_data(numpoints, numcurves):
    np.random.seed(3)
    x = np.linspace(0, 500, numpoints)

    height = 100 * np.random.random(numcurves)
    mu = 200 * np.random.random(numcurves) + 200
    sigma = 100 * np.random.random(numcurves) + 0.1
    data = gaussian(x, sigma, mu, height)

    noise = 5 * (np.random.random(data.shape) - 0.5)
    return x, data + noise

def gaussian(x, sigma, mu, height):
    data = -np.subtract.outer(x, mu)**2 / (2 * sigma**2)
    return height * np.exp(data)

def plot(x, ydata, ax=None, **kwargs):
    if ax is None:
        ax = plt.gca()
    colorcycle = itertools.cycle(mpl.rcParams['axes.color_cycle'])
    for y, color in zip(ydata.T, colorcycle):
        ax.plot(x, y, color=color, **kwargs)

main()

在此处输入图像描述

对于并行版本,我们唯一需要更改的是 main 函数。(我们还需要一个虚拟函数,因为multiprocessing.Pool.imap无法为其函数提供额外的参数......)它看起来像这样:

def parallel_main():
    import multiprocessing
    p = multiprocessing.Pool()
    x, data = generate_data(256, 262144)
    args = itertools.izip(itertools.repeat(x), data.T)
    model = p.imap(parallel_func, args, chunksize=500)
    sigma, mu, height = [np.array(item) for item in zip(*model)]
    prediction = gaussian(x, sigma, mu, height)

def parallel_func(args):
    return invert(*args)

编辑:如果简单的多项式拟合效果不佳,请尝试通过 y 值对问题进行加权,如@tslisten 共享的链接/论文中所述(虽然我的实现有点不同的)。

import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import itertools

def main():
    def run(x, data, func, threshold=0):
        model = [func(x, y, threshold=threshold) for y in data.T]
        sigma, mu, height = [np.array(item) for item in zip(*model)]
        prediction = gaussian(x, sigma, mu, height)

        plt.figure()
        plot(x, data, linestyle='none', marker='o', markersize=4)
        plot(x, prediction, linestyle='-', lw=2)

    x, data = generate_data(256, 6, noise=100)
    threshold = 50

    run(x, data, weighted_invert, threshold=threshold)
    plt.title('Weighted by Y-Value')

    run(x, data, invert, threshold=threshold)
    plt.title('Un-weighted Linear Inverse'

    plt.show()

def invert(x, y, threshold=0):
    mask = y > threshold
    x, y = x[mask], y[mask]

    # Fit a 2nd order polynomial to the log of the observed values
    A, B, C = np.polyfit(x, np.log(y), 2)

    # Solve for the desired parameters...
    sigma, mu, height = poly_to_gauss(A,B,C)
    return sigma, mu, height

def poly_to_gauss(A,B,C):
    sigma = np.sqrt(-1 / (2.0 * A))
    mu = B * sigma**2
    height = np.exp(C + 0.5 * mu**2 / sigma**2)
    return sigma, mu, height

def weighted_invert(x, y, weights=None, threshold=0):
    mask = y > threshold
    x,y = x[mask], y[mask]
    if weights is None:
        weights = y
    else:
        weights = weights[mask]

    d = np.log(y)
    G = np.ones((x.size, 3), dtype=np.float)
    G[:,0] = x**2
    G[:,1] = x

    model,_,_,_ = np.linalg.lstsq((G.T*weights**2).T, d*weights**2)
    return poly_to_gauss(*model)

def generate_data(numpoints, numcurves, noise=None):
    np.random.seed(3)
    x = np.linspace(0, 500, numpoints)

    height = 7000 * np.random.random(numcurves)
    mu = 1100 * np.random.random(numcurves) 
    sigma = 100 * np.random.random(numcurves) + 0.1
    data = gaussian(x, sigma, mu, height)

    if noise is None:
        noise = 0.1 * height.max()
    noise = noise * (np.random.random(data.shape) - 0.5)
    return x, data + noise

def gaussian(x, sigma, mu, height):
    data = -np.subtract.outer(x, mu)**2 / (2 * sigma**2)
    return height * np.exp(data)

def plot(x, ydata, ax=None, **kwargs):
    if ax is None:
        ax = plt.gca()
    colorcycle = itertools.cycle(mpl.rcParams['axes.color_cycle'])
    for y, color in zip(ydata.T, colorcycle):
        #kwargs['color'] = kwargs.get('color', color)
        ax.plot(x, y, color=color, **kwargs)

main()

在此处输入图像描述 在此处输入图像描述

如果这仍然给您带来麻烦,请尝试迭代地重新加权最小二乘问题(@tslisten 提到的链接中的最终“最佳”推荐方法)。但是请记住,这会慢得多。

def iterative_weighted_invert(x, y, threshold=None, numiter=5):
    last_y = y
    for _ in range(numiter):
        model = weighted_invert(x, y, weights=last_y, threshold=threshold)
        last_y = gaussian(x, *model)
    return model
于 2012-01-09T03:31:13.530 回答