Dati di caricamento Python e adattamento multi gaussiano

Question

Ho cercato un modo per eseguire più adattamenti gaussiani ai miei dati. La maggior parte degli esempi che ho trovato finora usano una distribuzione normale per creare numeri casuali. Ma sono interessato a guardare la trama dei miei dati e verificare se ci sono 1-3 picchi.

Posso farlo per un picco, ma non so come farlo per altro.

Ad esempio, ho questi dati: http://www.filedropper.com/data_11

Ho provato a usando lmfit, e, naturalmente SciPy, ma senza risultati piacevoli.

Grazie per qualsiasi aiuto!

Answer 1

È sufficiente eseguire funzioni di modello parametrizzate della somma di singoli gaussiani. Scegli un buon valore per la tua ipotesi iniziale (questo è un passaggio davvero cruciale) e poi hai scipy.optimize di modificare quei numeri un po '.

Ecco come si potrebbe farlo:

import numpy as np 
import matplotlib.pyplot as plt 
from scipy import optimize 

data = np.genfromtxt('data.txt') 
def gaussian(x, height, center, width, offset): 
    return height*np.exp(-(x - center)**2/(2*width**2)) + offset 
def three_gaussians(x, h1, c1, w1, h2, c2, w2, h3, c3, w3, offset): 
    return (gaussian(x, h1, c1, w1, offset=0) + 
     gaussian(x, h2, c2, w2, offset=0) + 
     gaussian(x, h3, c3, w3, offset=0) + offset) 

def two_gaussians(x, h1, c1, w1, h2, c2, w2, offset): 
    return three_gaussians(x, h1, c1, w1, h2, c2, w2, 0,0,1, offset) 

errfunc3 = lambda p, x, y: (three_gaussians(x, *p) - y)**2 
errfunc2 = lambda p, x, y: (two_gaussians(x, *p) - y)**2 

guess3 = [0.49, 0.55, 0.01, 0.6, 0.61, 0.01, 1, 0.64, 0.01, 0] # I guess there are 3 peaks, 2 are clear, but between them there seems to be another one, based on the change in slope smoothness there 
guess2 = [0.49, 0.55, 0.01, 1, 0.64, 0.01, 0] # I removed the peak I'm not too sure about 
optim3, success = optimize.leastsq(errfunc3, guess3[:], args=(data[:,0], data[:,1])) 
optim2, success = optimize.leastsq(errfunc2, guess2[:], args=(data[:,0], data[:,1])) 
optim3 

plt.plot(data[:,0], data[:,1], lw=5, c='g', label='measurement') 
plt.plot(data[:,0], three_gaussians(data[:,0], *optim3), 
    lw=3, c='b', label='fit of 3 Gaussians') 
plt.plot(data[:,0], two_gaussians(data[:,0], *optim2), 
    lw=1, c='r', ls='--', label='fit of 2 Gaussians') 
plt.legend(loc='best') 
plt.savefig('result.png') 

Come si può vedere, non v'è quasi alcuna differenza tra i due attacchi (visivamente). Così non si può sapere con certezza se c'erano 3 gaussiane presenti nella sorgente o solo 2. Tuttavia, se si dovesse fare una congettura, quindi controllare per il più piccolo residuo:

err3 = np.sqrt(errfunc3(optim3, data[:,0], data[:,1])).sum() 
err2 = np.sqrt(errfunc2(optim2, data[:,0], data[:,1])).sum() 
print('Residual error when fitting 3 Gaussians: {}\n' 
    'Residual error when fitting 2 Gaussians: {}'.format(err3, err2)) 
# Residual error when fitting 3 Gaussians: 3.52000910965 
# Residual error when fitting 2 Gaussians: 3.82054499044 

In questo caso, 3 I gaussiani danno un risultato migliore, ma ho anche fatto la mia ipotesi iniziale abbastanza accurata.