Python: ajustement gaussien à deux courbes avec les moindres carrés non linéaires

ma connaissance des mathématiques est limitée, c'est pourquoi je suis probablement coincé. J'ai un spectre auquel j'essaie de coller deux pics gaussiens. Je peux m'adapter au plus grand pic, mais je ne peux pas m'adapter au plus petit pic. Je comprends que j'ai besoin de résumer la fonction gaussienne pour les deux pics, mais je ne sais pas où je me suis trompé. Une image de ma sortie actuelle est affichée:

Current Output

La ligne bleue est mes données et la ligne verte est mon ajustement. Il y a un Epaule à gauche du pic principal de mes données que j'essaie actuellement d'ajuster, en utilisant le code suivant:

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()
17
demandé sur Usagi 2012-04-13 19:36:50

3 réponses

ce code a fonctionné pour moi à condition que vous ne positionniez qu'une fonction qui est une combinaison de deux distributions gaussiennes.

je viens de faire une fonction de résidus qui ajoute deux fonctions gaussiennes et puis les soustrait des données réelles.

les paramètres (p) que j'ai passés à la fonction des moindres carrés de Numpy comprennent: la moyenne de la première fonction gaussienne (m), la différence de la moyenne de la première et de la deuxième fonctions gaussienne (dm, c.-à-d. l'horizontale maj), l'écart-type de la première (sd1), et l'écart-type de la deuxième (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()

Results of the code.

16
répondu Usagi 2013-07-02 19:58:05

vous pouvez utiliser des modèles de mélange gaussien 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])

enter image description here

vous pouvez aussi utiliser la fonction ci-dessous pour ajuster le nombre de gaussiens que vous voulez avec le paramètre 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
12
répondu bougui 2014-05-28 07:47:53

les coefficients 0 et 4 sont dégénérés - il n'y a absolument rien dans les données qui puisse décider entre eux. vous devez utiliser un seul paramètre de niveau zéro au lieu de deux (c'est-à-dire supprimer l'un d'eux de votre code). c'est probablement ce qui arrête votre ajustement (ignorez les commentaires ici disant que ce n'est pas possible - il y a clairement au moins deux pics dans ces données et vous devriez certainement être en mesure de s'adapter à cela).

(il peut ne pas être clair pourquoi je suggère, mais ce qui se passe est-ce que les coefficients 0 et 4 peuvent s'annuler mutuellement. ils peuvent tous les deux être zéro, ou l'un pourrait être 100 et l'autre -100 - dans tous les cas, l'ajustement est tout aussi bon. cela "confond" la routine appropriée, qui passe son temps à essayer de déterminer ce qu'ils devraient être, quand il n'y a pas de réponse unique, parce que quelle que soit la valeur de l'un, l'autre peut être juste le négatif de cela, et l'ajustement sera le même).

en fait, de l'intrigue, on dirait qu'il ya peut-être pas besoin d'un niveau zéro à tous. j'essaierais de laisser tomber les deux et de voir à quoi ressemble la forme.

4
répondu andrew cooke 2012-04-14 01:54:55