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

Mes connaissances en maths est limitée, c'est pourquoi je suis probablement bloqué. J'ai un spectres pour laquelle je suis en train de monter deux pics Gaussiens. Je peux adapter pour le plus grand pic, mais je ne peux pas adapter à la plus petite pointe. Je comprends que j'ai besoin de la somme de la fonction Gaussienne pour les deux pics, mais je ne sais pas où je suis allé mal. Une image de ma sortie de courant est indiqué:

La ligne bleue est mes données et la ligne verte est mon ajustement. Il y a une épaule à gauche du pic principal dans mes données, je suis en train d'essayer de s'adapter, 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()

source d'informationauteur Harpal

Ce code a fonctionné pour moi à condition que vous soyez seul côté d'une fonction qui est une combinaison de deux distributions Gaussiennes.

Je viens de faire un résidus de la fonction qui ajoute deux fonctions Gaussiennes, puis soustraire les données réelles.

Les paramètres (p) que j'ai passé à Numpy la méthode des moindres carrés: la moyenne de la première fonction Gaussienne (m), la différence dans la moyenne de la première et de la deuxième Gaussienne fonctions (dm, c'est à dire le décalage horizontal), 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()

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

Vous pouvez utiliser des modèles de mélanges Gaussiens 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])

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

Vous pouvez également utiliser la fonction ci-dessous pour ajuster le nombre de Gaussiennes vous voulez avec ncomp paramètre:

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

4

coeffs 0 et 4 sont dégénérés - il n'y a absolument rien dans les données qui peuvent décider entre eux de la. vous devez utiliser un seul niveau zéro du paramètre au lieu de deux (c'est à dire enlever l'un d'entre eux à partir de votre code). c'est probablement ce qui est de l'arrêt de votre ajustement (ignorer les commentaires ici en disant que ce n'est pas possible - il y a clairement des au moins deux pics dans les données et vous devez certainement être en mesure d'ajustement).

(il peut ne pas être clair pourquoi je suggère, mais ce qui se passe est que les coeffs 0 et 4 peuvent annuler les uns les autres. ils peuvent à la fois être de zéro, ou l'on pouvait être à 100, et les autres -100 - de toute façon, l'ajustement est tout aussi bon. cette "confond" l'aménagement de la routine, qui passe son temps à essayer de travailler sur ce qu'ils devraient être, quand il n'y a pas de bonne réponse, parce que quelle que soit la valeur de l'un est, l'autre peut-être juste le négatif, et l'ajustement sera le même).

en fait, à partir de l'intrigue, on dirait qu'il ya peut-être pas besoin d'un niveau zéro à tous. je voudrais essayer de faire tomber les deux de ceux-ci et de voir comment l'ajustement semble.

aussi, il n'est pas nécessaire pour s'adapter à coeffs 1 et 5 (ou le point zéro) dans la méthode des moindres carrés. au lieu de cela, parce que le modèle est linéaire dans ceux que vous pourriez calculer leurs valeurs de chaque boucle. cela va rendre les choses plus vite, mais n'est pas critique. je viens de remarquer que vous dites en mathématiques n'est pas très bon, donc probablement ignorer celui-ci.

Vous devez vous connecter pour publier un commentaire.