Python ajustement de la courbe de bibliothèque qui me permet d'assigner des limites pour les paramètres

J'aimerais être capable d'effectuer des ajustements qui me permet d'ajuster une courbe arbitraire fonction des données, et me permet de fixer arbitrairement les limites sur les paramètres, par exemple: je veux ajustement de la fonction:

f(x) = a1(x-a2)^a3\cdot\exp(-\a4*x^a5)

et de dire:

  • a2 est dans la plage suivante: (-1, 1)
  • a3 et a5 sont positifs

Il est agréable scipy curve_fit fonction, mais il ne permet pas de spécifier le paramètre de limite. Il est aussi agréable http://code.google.com/p/pyminuit/ bibliothèque qui ne générique minimalization, et qu'elle permet de fixer des limites sur les paramètres, mais dans mon cas, il n'a pas coverge.

Voici une question similaire, avec quelques réponses utiles: stackoverflow.com/questions/22895794/...

OriginalL'auteur jb. | 2013-05-26

  1. 10

    Remarque: la Nouvelle version de 0,17 SciPy

    Supposons que vous souhaitez pour s'adapter à un modèle de données qui ressemble à ceci:

    y=a*t**alpha+b

    et avec la contrainte sur alpha

    0<alpha<2

    tandis que d'autres paramètres a et b reste libre. Ensuite, il faut utiliser les bornes de l'option de curve_fit de la manière suivante:

    import numpy as np
from scipy.optimize import curve_fit
def func(t, a,alpha,b):
     return a*t**alpha+b
param_bounds=([-np.inf,0,-np.inf],[np.inf,2,np.inf])
popt, pcov = curve_fit(func, xdata, ydata,bounds=param_bounds)

    Source est ici.

    Je suppose que c'est maintenant la bonne réponse! Merci!
    Puis-je définir le paramètre alpha gamme que le nombre entier d'intervalles au lieu de flotter? Exemple α = { 0, 1, 2, 3, 4, 5, ... }

    OriginalL'auteur Mohit Pandey

  2. 3

    Solution de contournement: utiliser la variable transformations comme a2=tanh(a2), a3=exp(a3') ou a5=a5'^2.

    OriginalL'auteur Joachim W

  3. 3

    Comme déjà mentionné par Rob Falck, vous pouvez utiliser, par exemple, la scipy d'optimisation non linéaire des routines dans scipy.minimiser pour réduire l'arbitraire de la fonction d'erreur, par exemple, l'erreur quadratique moyenne.

    Notez que la fonction que vous avez donné n'est pas nécessairement vraies valeurs - c'était peut être la raison de votre minimisation dans pyminuit ne converge pas. Vous avez à traiter un peu plus explicite, voir l'exemple 2.

    Les exemples ci-dessous utilisent tous les deux la L-BFGS-B la minimisation de la méthode, qui prend en charge délimitée paramètre régions. J'ai divisé cette réponse en deux parties:

    1. Une fonction avec de vrais arrivée, ressemblant à celui donné par vous. J'ai ajouté absolus pour assurer la fonction que vous avez donné renvoie des nombres réels dans le domaine [-3,3)
    2. De la fonction réelle que vous avez donné, qui a un complexe d'arrivée

    1. Véritable arrivée

    L'exemple ci-dessous montre l'optimisation de cette version légèrement modifiée de votre fonction.

    import numpy as np
import pylab as pl
from scipy.optimize import minimize
points = 500
xlim = 3.
def f(x,*p):
a1,a2,a3,a4,a5 = p
return a1*np.abs(x-a2)**a3 * np.exp(-a4 * np.abs(x)**a5)
# generate noisy data with known coefficients
p0 = [1.4,-.8,1.1,1.2,2.2]
x = (np.random.rand(points) * 2. - 1.) * xlim
x.sort()
y = f(x,*p0)
y_noise = y + np.random.randn(points) * .05
# mean squared error wrt. noisy data as a function of the parameters
err = lambda p: np.mean((f(x,*p)-y_noise)**2)
# bounded optimization using scipy.minimize
p_init = [1.,-1.,.5,.5,2.]
p_opt = minimize(
err, # minimize wrt to the noisy data
p_init, 
bounds=[(None,None),(-1,1),(None,None),(0,None),(None,None)], # set the bounds
method="L-BFGS-B" # this method supports bounds
).x
# plot everything
pl.scatter(x, y_noise, alpha=.2, label="f + noise")
pl.plot(x, y, c='#000000', lw=2., label="f")
pl.plot(x, f(x,*p_opt) ,'--', c='r', lw=2., label="fitted f")
pl.xlabel("x")
pl.ylabel("f(x)")
pl.legend(loc="best")
pl.xlim([-xlim*1.01,xlim*1.01])
pl.show()

    Python ajustement de la courbe de bibliothèque qui me permet d'assigner des limites pour les paramètres

    2. Extension du complexe arrivée

    Extension de la au-dessus de la minimisation de la domaine complexe qui peut être fait explicitement la conversion des nombres complexes et l'adaptation de la fonction d'erreur:

    Tout d'abord, vous exprimées explicitement la valeur de x à valeurs complexes f pour assurer des rendements des valeurs complexes et peuvent calculer exposant fractionnaire des nombres négatifs. Deuxièmement, nous calculons l'erreur de la fonction à la fois réel et imaginaire - un simple candidat est la moyenne du carré de la complexe absolus.

    import numpy as np
import pylab as pl
from scipy.optimize import minimize
points = 500
xlim = 3.
def f(x,*p):
a1,a2,a3,a4,a5 = p
x = x.astype(complex) # cast x explicitly to complex, to ensure complex valued f
return a1*(x-a2)**a3 * np.exp(-a4 * x**a5)
# generate noisy data with known coefficients
p0 = [1.4,-.8,1.1,1.2,2.2]
x = (np.random.rand(points) * 2. - 1.) * xlim
x.sort()
y = f(x,*p0)
y_noise = y + np.random.randn(points) * .05 + np.random.randn(points) * 1j*.05
# error function chosen as mean of squared absolutes
err = lambda p: np.mean(np.abs(f(x,*p)-y_noise)**2)
# bounded optimization using scipy.minimize
p_init = [1.,-1.,.5,.5,2.]
p_opt = minimize(
err, # minimize wrt to the noisy data
p_init, 
bounds=[(None,None),(-1,1),(None,None),(0,None),(None,None)], # set the bounds
method="L-BFGS-B" # this method supports bounds
).x
# plot everything
pl.scatter(x, np.real(y_noise), c='b',alpha=.2, label="re(f) + noise")
pl.scatter(x, np.imag(y_noise), c='r',alpha=.2, label="im(f) + noise")
pl.plot(x, np.real(y), c='b', lw=1., label="re(f)")
pl.plot(x, np.imag(y), c='r', lw=1., label="im(f)")
pl.plot(x, np.real(f(x,*p_opt)) ,'--', c='b', lw=2.5, label="fitted re(f)")
pl.plot(x, np.imag(f(x,*p_opt)) ,'--', c='r', lw=2.5, label="fitted im(f)")
pl.xlabel("x")
pl.ylabel("f(x)")
pl.legend(loc="best")
pl.xlim([-xlim*1.01,xlim*1.01])
pl.show()

    Python ajustement de la courbe de bibliothèque qui me permet d'assigner des limites pour les paramètres

    Notes

    Il semble que le minimiseur peut-être un peu sensible aux valeurs initiales - j'ai donc passé ma première hypothèse (p_init) pas trop loin de l'optimum. Devriez-vous avoir à lutter avec cela, vous pouvez utiliser la même procédure de minimisation en plus d'une optimisation globale de la boucle, par exemple bassin de saut ou brute.

    OriginalL'auteur flinz

  4. 3

    Vous pouvez utiliser lmfit pour ce genre de problèmes. Donc, j'ajoute un exemple (avec une autre fonction que vous utiliser, mais il peut facilement l'adapter) sur la façon de l'utiliser au cas où quelqu'un est intéressé par ce sujet, trop.

    Disons que vous avez un jeu de données comme suit:

    xdata = np.array([177.,180.,183.,187.,189.,190.,196.,197.,201.,202.,203.,204.,206.,218.,225.,231.,234.,
252.,262.,266.,267.,268.,277.,286.,303.])
ydata = np.array([0.81,0.74,0.78,0.75,0.77,0.81,0.73,0.76,0.71,0.74,0.81,0.71,0.74,0.71,
0.72,0.69,0.75,0.59,0.61,0.63,0.64,0.63,0.35,0.27,0.26])

    et que vous souhaitez adapter un modèle de données qui ressemble à ceci:

    model = n1 + (n2 * x + n3) * 1./ (1. + np.exp(n4 * (n5 - x)))

    avec les contraintes que 

    0.2 < n1 < 0.8
-0.3 < n2 < 0

    À l'aide de lmfit (version 0.8.3) vous avez alors obtenir le résultat suivant:

    n1:   0.26564921 +/- 0.024765 (9.32%) (init= 0.2)
n2:  -0.00195398 +/- 0.000311 (15.93%) (init=-0.005)
n3:   0.87261892 +/- 0.068601 (7.86%) (init= 1.0766)
n4:  -1.43507072 +/- 1.223086 (85.23%) (init=-0.36379)
n5:   277.684530 +/- 3.768676 (1.36%) (init= 274)

    Python ajustement de la courbe de bibliothèque qui me permet d'assigner des limites pour les paramètres

    Comme vous pouvez le voir, l'ajustement reproduit les données très bien et les paramètres sont dans l'demandé plages.

    Ici est l'ensemble du code qui reproduit l'intrigue avec quelques commentaires supplémentaires:

    from lmfit import minimize, Parameters, Parameter, report_fit
import numpy as np
xdata = np.array([177.,180.,183.,187.,189.,190.,196.,197.,201.,202.,203.,204.,206.,218.,225.,231.,234.,
252.,262.,266.,267.,268.,277.,286.,303.])
ydata = np.array([0.81,0.74,0.78,0.75,0.77,0.81,0.73,0.76,0.71,0.74,0.81,0.71,0.74,0.71,
0.72,0.69,0.75,0.59,0.61,0.63,0.64,0.63,0.35,0.27,0.26])
def fit_fc(params, x, data):
n1 = params['n1'].value
n2 = params['n2'].value
n3 = params['n3'].value
n4 = params['n4'].value
n5 = params['n5'].value
model = n1 + (n2 * x + n3) * 1./ (1. + np.exp(n4 * (n5 - x)))
return model - data #that's what you want to minimize
# create a set of Parameters
# 'value' is the initial condition
# 'min' and 'max' define your boundaries
params = Parameters()
params.add('n1', value= 0.2, min=0.2, max=0.8)
params.add('n2', value= -0.005, min=-0.3, max=10**(-10))
params.add('n3', value= 1.0766, min=-1000., max=1000.)
params.add('n4', value= -0.36379, min=-1000., max=1000.)
params.add('n5', value= 274.0, min=0., max=1000.)
# do fit, here with leastsq model
result = minimize(fit_fc, params, args=(xdata, ydata))
# write error report
report_fit(params)
xplot = np.linspace(min(xdata), max(xdata), 1000)
yplot = result.values['n1'] + (result.values['n2'] * xplot + result.values['n3']) * \
1./ (1. + np.exp(result.values['n4'] * (result.values['n5'] - xplot)))
#plot results
try:
import pylab
pylab.plot(xdata, ydata, 'k+')
pylab.plot(xplot, yplot, 'r')
pylab.show()
except:
pass

    EDIT:

    Si vous utilisez la version 0.9.x vous avez besoin de modifier le code en conséquence; vérifier ici dont les changements ont été faits à partir de 0.8.3 à 0,9.x.

    OriginalL'auteur Cleb

  5. 1

    Avez-vous envisagé de le traiter comme un problème d'optimisation et d'aide de l'une des routines d'optimisation non linéaire dans scipy pour minimiser la méthode des moindres carrés de l'erreur en faisant varier les coefficients de votre fonction? De nombreuses routines d'en optimiser permettre lié à des contraintes sur les variables indépendantes.

    OriginalL'auteur Rob Falck