programar - Python: ajuste gaussiano de dos curvas con mínimos cuadrados no lineales

regresion lineal python numpy

Question

Mi conocimiento de las matemáticas es limitado y es por eso que probablemente estoy atascado. Tengo un espectro al que intento ajustar dos picos Gaussianos. Puedo caber en el pico más grande, pero no puedo encajar en el pico más pequeño. Entiendo que necesito sumar la función gaussiana para los dos picos, pero no sé dónde me he equivocado. Se muestra una imagen de mi salida actual:

La línea azul es mi información y la línea verde es mi ajuste actual. Hay un hombro a la izquierda del pico principal en mis datos que estoy tratando de encajar, utilizando el siguiente código:

import matplotlib.pyplot as pt
import numpy as np
from scipy.optimize import leastsq
from pylab import *

time = []
counts = []


for i in open(''/some/folder/to/file.txt'', ''r''):
    segs = i.split()
    time.append(float(segs[0]))
    counts.append(segs[1])

time_array = arange(len(time), dtype=float)
counts_array = arange(len(counts))
time_array[0:] = time
counts_array[0:] = counts


def model(time_array0, coeffs0):
    a = coeffs0[0] + coeffs0[1] * np.exp( - ((time_array0-coeffs0[2])/coeffs0[3])**2 )
    b = coeffs0[4] + coeffs0[5] * np.exp( - ((time_array0-coeffs0[6])/coeffs0[7])**2 ) 
    c = a+b
    return c


def residuals(coeffs, counts_array, time_array):
    return counts_array - model(time_array, coeffs)

# 0 = baseline, 1 = amplitude, 2 = centre, 3 = width
peak1 = np.array([0,6337,16.2,4.47,0,2300,13.5,2], dtype=float)
#peak2 = np.array([0,2300,13.5,2], dtype=float)

x, flag = leastsq(residuals, peak1, args=(counts_array, time_array))
#z, flag = leastsq(residuals, peak2, args=(counts_array, time_array))

plt.plot(time_array, counts_array)
plt.plot(time_array, model(time_array, x), color = ''g'') 
#plt.plot(time_array, model(time_array, z), color = ''r'')
plt.show()

Answer 1

los coeficientes 0 y 4 son degenerados: no hay absolutamente nada en los datos que pueda decidir entre ellos. debe usar un solo parámetro de nivel cero en lugar de dos (es decir, eliminar uno de ellos de su código). esto es probablemente lo que está deteniendo su ajuste (ignore los comentarios aquí diciendo que esto no es posible - claramente hay al menos dos picos en esa información y ciertamente debería poder encajar en eso).

(Puede no estar claro por qué estoy sugiriendo esto, pero lo que está sucediendo es que los coeficientes 0 y 4 pueden cancelarse mutuamente. Ambos pueden ser cero, o uno podría ser 100 y el otro -100 - en cualquier caso, el ajuste es tan bueno. Esto "confunde" la rutina de ajuste, que pasa su tiempo tratando de resolver lo que deberían ser, cuando no hay una sola respuesta correcta, porque cualquiera que sea el valor uno, el otro puede ser el negativo de eso, y el ajuste será el mismo).

de hecho, de la trama, parece que no puede haber necesidad de un nivel cero en absoluto. Trataría de eliminar ambos y ver cómo se ve el ajuste.

Además, no es necesario ajustar los coeficientes 1 y 5 (o el punto cero) en los mínimos cuadrados. en cambio, debido a que el modelo es lineal en los que puede calcular sus valores en cada ciclo. esto hará las cosas más rápidas, pero no es crítico. Me acabo de dar cuenta de que dices que tus cálculos no son tan buenos, así que probablemente ignores este.

Answer 2

Este código funcionó para mí siempre que solo esté ajustando una función que es una combinación de dos distribuciones gaussianas.

Acabo de crear una función residual que agrega dos funciones gaussianas y luego las resta de los datos reales.

Los parámetros (p) que pasé a la función de mínimos cuadrados de Numpy incluyen: la media de la primera función gaussiana (m), la diferencia en la media de la primera y segunda función gaussiana (dm, es decir, el desplazamiento horizontal), la desviación estándar de la primera (sd1), y la desviación estándar de la segunda (sd2).

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

######################################
# Setting up test data
def norm(x, mean, sd):
  norm = []
  for i in range(x.size):
    norm += [1.0/(sd*np.sqrt(2*np.pi))*np.exp(-(x[i] - mean)**2/(2*sd**2))]
  return np.array(norm)

mean1, mean2 = 0, -2
std1, std2 = 0.5, 1 

x = np.linspace(-20, 20, 500)
y_real = norm(x, mean1, std1) + norm(x, mean2, std2)

######################################
# Solving
m, dm, sd1, sd2 = [5, 10, 1, 1]
p = [m, dm, sd1, sd2] # Initial guesses for leastsq
y_init = norm(x, m, sd1) + norm(x, m + dm, sd2) # For final comparison plot

def res(p, y, x):
  m, dm, sd1, sd2 = p
  m1 = m
  m2 = m1 + dm
  y_fit = norm(x, m1, sd1) + norm(x, m2, sd2)
  err = y - y_fit
  return err

plsq = leastsq(res, p, args = (y_real, x))

y_est = norm(x, plsq[0][0], plsq[0][2]) + norm(x, plsq[0][0] + plsq[0][1], plsq[0][3])

plt.plot(x, y_real, label=''Real Data'')
plt.plot(x, y_init, ''r.'', label=''Starting Guess'')
plt.plot(x, y_est, ''g.'', label=''Fitted'')
plt.legend()
plt.show()

Answer 3

Puede usar modelos de mezclas gaussianas de scikit-learn :

from sklearn import mixture
import matplotlib.pyplot
import matplotlib.mlab
import numpy as np
clf = mixture.GMM(n_components=2, covariance_type=''full'')
clf.fit(yourdata)
m1, m2 = clf.means_
w1, w2 = clf.weights_
c1, c2 = clf.covars_
histdist = matplotlib.pyplot.hist(yourdata, 100, normed=True)
plotgauss1 = lambda x: plot(x,w1*matplotlib.mlab.normpdf(x,m1,np.sqrt(c1))[0], linewidth=3)
plotgauss2 = lambda x: plot(x,w2*matplotlib.mlab.normpdf(x,m2,np.sqrt(c2))[0], linewidth=3)
plotgauss1(histdist[1])
plotgauss2(histdist[1])

También puede usar la función siguiente para ajustar el número de gaussianos que desea con el parámetro ncomp:

from sklearn import mixture
%pylab

def fit_mixture(data, ncomp=2, doplot=False):
    clf = mixture.GMM(n_components=ncomp, covariance_type=''full'')
    clf.fit(data)
    ml = clf.means_
    wl = clf.weights_
    cl = clf.covars_
    ms = [m[0] for m in ml]
    cs = [numpy.sqrt(c[0][0]) for c in cl]
    ws = [w for w in wl]
    if doplot == True:
        histo = hist(data, 200, normed=True)
        for w, m, c in zip(ws, ms, cs):
            plot(histo[1],w*matplotlib.mlab.normpdf(histo[1],m,np.sqrt(c)), linewidth=3)
    return ms, cs, ws