Come posso eseguire un montaggio dei minimi quadrati su più set di dati velocemente?

Question

Sto provando a fare un adattamento gaussiano su molti punti dati. Per esempio. Ho una matrice di dati 256 x 262144. Dove i 256 punti devono essere adattati a una distribuzione gaussiana, e ho bisogno di 262144 di essi.

A volte il picco della distribuzione gaussiana si trova al di fuori dell'intervallo di dati, quindi per ottenere un risultato medio accurato l'adattamento della curva rappresenta l'approccio migliore. Anche se il picco si trova all'interno dell'intervallo, l'adattamento della curva fornisce un sigma migliore perché altri dati non sono compresi nell'intervallo.

Ho questo funzionamento per un punto di dati, utilizzando il codice da http://www.scipy.org/Cookbook/FittingData.

Ho provato a ripetere questo algoritmo, ma sembra che occorrerà qualcosa nell'ordine di 43 minuti per risolvere questo problema. C'è un modo veloce già scritto per farlo in parallelo o in modo più efficiente?

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                                    

Si noti che i dati non sono necessariamente 256x262144 e ho fatto qualche giro in fudging nd_fit per fare questo lavoro.

Il codice che uso per ottenere questo al lavoro è

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 soluzione postato qui sotto per @JoeKington è grande e risolve veramente veloce. Tuttavia, non sembra funzionare a meno che l'area significativa del gaussiano sia all'interno dell'area appropriata. Dovrò verificare se la media è ancora precisa, dato che questa è la cosa principale per cui la uso.

Answer 1

La cosa più semplice da fare è risolvere il problema. Stai utilizzando un metodo iterativo non lineare che sarà più lento di una soluzione lineare di minimi quadrati.

In sostanza, si ha:

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

per rendere questa un'equazione lineare, prendere la) log (naturale di entrambe le parti:

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

Questo semplifica poi a il polinomio:

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

Possiamo riordinarlo in un modo un po 'più semplice :

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

dove:

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

Tuttavia, c'è una cattura. Questo diventerà instabile in presenza di rumore nelle "code" della distribuzione.

Pertanto, è necessario utilizzare solo i dati vicino ai "picchi" della distribuzione. È abbastanza facile includere solo i dati che cadono al di sopra di qualche soglia nel fitting. In questo esempio, sto includendo solo i dati che sono maggiori del 20% del valore massimo osservato per una data curva gaussiana che stiamo adattando.

Una volta terminato, è piuttosto veloce. La risoluzione di 262144 curve gaussiane diverse richiede solo ~ 1 minuto (assicurarsi di rimuovere la parte di tracciamento del codice se la si esegue su qualcosa di così grande ...). È anche abbastanza facile parallelizzare, se vuoi ...

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

L'unica cosa che avremmo bisogno di cambiare per una versione parallela è la funzione principale. (Abbiamo anche bisogno di una funzione fittizia perché multiprocessing.Pool.imap non può fornire ulteriori argomenti alla sua funzione ...) Sarebbe sembrare qualcosa di simile:

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) 

Edit: Nei casi in cui il semplice montaggio polinomiale non è lavorando bene, prova a soppesare il problema con i valori y, as mentioned in the link/paper condivisi da @tslisten (e Stefan van der Walt implementato, sebbene la mia implementazione sia un po 'diversa).

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

Se sei tu ancora dando problemi, quindi provare in modo iterativo-riponderazione il problema dei minimi quadrati (il "migliore" metodo CONSIGLIATI finale nel collegamento @tslisten citato). Tieni presente che questo sarà considerevolmente più lento, comunque.

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