¿Cómo puedo realizar un ajuste de mínimos cuadrados sobre múltiples conjuntos de datos rápidamente?

Question

Estoy tratando de hacer un ajuste gaussiano sobre muchos puntos de datos. P.ej. Tengo una matriz de datos de 256 x 262144. Donde los 256 puntos deben ajustarse a una distribución gaussiana, y necesito 262144 de ellos.

A veces, el pico de la distribución gaussiana está fuera del rango de datos, por lo que para obtener un resultado de curva de resultado medio exacto es el mejor enfoque. Incluso si el pico está dentro del rango, el ajuste de curva da una sigma mejor porque otros datos no están en el rango.

Tengo esto funcionando para un punto de datos, usando el código de http://www.scipy.org/Cookbook/FittingData.

He tratado de repetir este algoritmo, pero parece que va a tomar algo del orden de 43 minutos para resolverlo. ¿Existe una forma rápida ya escrita de hacerlo en paralelo o de manera más eficiente?

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                                    

Tenga en cuenta que los datos no es necesariamente 256x262144 y he hecho algunas vueltas en fudging nd_fit para hacer este trabajo.

El código que utilizo para conseguir que esto funcione es

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) 

Nota: La solución publicado por debajo @JoeKington es grande y resuelve muy rápido. Sin embargo, no parece funcionar a menos que el área significativa del gaussiano esté dentro del área apropiada. Tendré que probar si la media todavía es precisa, ya que es lo principal para lo que uso esto.

Answer 1

Lo más fácil de hacer es linearizar el problema. Está utilizando un método iterativo no lineal que será más lento que una solución lineal de mínimos cuadrados.

Básicamente, usted tiene:

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

Para que esto sea una ecuación lineal, tomar el logaritmo (natural) de ambos lados:

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

Esto entonces se simplifica a el polinomio:

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

Podemos refundir esto en una forma un poco más simple :

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

donde:

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

Sin embargo, hay un inconveniente. Esto se volverá inestable en presencia de ruido en las "colas" de la distribución.

Por lo tanto, necesitamos usar solo los datos cerca de los "picos" de la distribución. Es bastante fácil incluir solo datos que caigan por encima de algún umbral en el ajuste. En este ejemplo, solo incluyo datos que son más del 20% del valor máximo observado para una curva gaussiana determinada que estamos ajustando.

Una vez que hemos hecho esto, sin embargo, es bastante rápido. La resolución de 262144 curvas gaussianas diferentes toma solo ~ 1 minuto (asegúrese de eliminar la parte de trazado del código si lo ejecuta en algo tan grande ...). También es bastante fácil de paralelizar, si quieres ...

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() 

Lo único que tendríamos que cambiar de una versión paralela es la función principal. (Necesitamos también una función ficticia porque multiprocessing.Pool.imap no puede proporcionar argumentos adicionales a su función de ...) se vería algo como esto:

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) 

Editar: En los casos en que el simple encaje polinomio no es funcionando bien, intente ponderar el problema con los valores y, as mentioned in the link/paper que @tslisten compartió (y Stefan van der Walt implementó, aunque mi implementación es un poco diferente).

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() 

Si aún así te está dando problemas, y luego tratar de forma iterativa-reponderación el problema de mínimos cuadrados (El último método "mejor" recomendados en el enlace @tslisten mencionado). Sin embargo, tenga en cuenta que esto será considerablemente más lento.

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