library python numpy scipy

library - scipy python



ajuste de la caída exponencial sin adivinanzas iniciales (8)

¿Alguien sabe un módulo scipy / numpy que permita ajustar la caída exponencial a los datos?

La búsqueda de Google devolvió algunas publicaciones de blog, por ejemplo, http://exnumerus.blogspot.com/2010/04/how-to-fit-exponential-decay-example-in.html , pero esa solución requiere que y-offset sea pre-especificado, que no siempre es posible

EDITAR:

curve_fit funciona, pero puede fallar de manera muy miserable sin una estimación inicial de los parámetros, y eso a veces es necesario. El código con el que estoy trabajando es

#!/usr/bin/env python import numpy as np import scipy as sp import pylab as pl from scipy.optimize.minpack import curve_fit x = np.array([ 50., 110., 170., 230., 290., 350., 410., 470., 530., 590.]) y = np.array([ 3173., 2391., 1726., 1388., 1057., 786., 598., 443., 339., 263.]) smoothx = np.linspace(x[0], x[-1], 20) guess_a, guess_b, guess_c = 4000, -0.005, 100 guess = [guess_a, guess_b, guess_c] exp_decay = lambda x, A, t, y0: A * np.exp(x * t) + y0 params, cov = curve_fit(exp_decay, x, y, p0=guess) A, t, y0 = params print "A = %s/nt = %s/ny0 = %s/n" % (A, t, y0) pl.clf() best_fit = lambda x: A * np.exp(t * x) + y0 pl.plot(x, y, ''b.'') pl.plot(smoothx, best_fit(smoothx), ''r-'') pl.show()

que funciona, pero si eliminamos "p0 = guess", falla miserablemente.

Implementación en Python de la solución de @JJacquelin. Necesitaba una solución aproximada no basada en la resolución sin adivinanzas iniciales, por lo que la respuesta de @JJacquelin fue realmente útil. La pregunta original se planteó como una solicitud de python numpy / scipy. Tomé el código R limpio y agradable de @ johanvdw y lo reformulé como python / numpy. Esperemos que sea útil para alguien: https://gist.github.com/friendtogeoff/00b89fa8d9acc1b2bdf3bdb675178a29

import numpy as np """ compute an exponential decay fit to two vectors of x and y data result is in form y = a + b * exp(c*x). ref. https://gist.github.com/johanvdw/443a820a7f4ffa7e9f8997481d7ca8b3 """ def exp_est(x,y): n = np.size(x) # sort the data into ascending x order y = y[np.argsort(x)] x = x[np.argsort(x)] Sk = np.zeros(n) for n in range(1,n): Sk[n] = Sk[n-1] + (y[n] + y[n-1])*(x[n]-x[n-1])/2 dx = x - x[0] dy = y - y[0] m1 = np.matrix([[np.sum(dx**2), np.sum(dx*Sk)], [np.sum(dx*Sk), np.sum(Sk**2)]]) m2 = np.matrix([np.sum(dx*dy), np.sum(dy*Sk)]) [d, c] = (m1.I * m2.T).flat m3 = np.matrix([[n, np.sum(np.exp( c*x))], [np.sum(np.exp(c*x)),np.sum(np.exp(2*c*x))]]) m4 = np.matrix([np.sum(y), np.sum(y*np.exp(c*x).T)]) [a, b] = (m3.I * m4.T).flat return [a,b,c]

La forma correcta de hacerlo es hacer una estimación de Prony y usar el resultado como la conjetura inicial para el ajuste de mínimos cuadrados (o alguna otra rutina de ajuste más robusta). La estimación de Prony no necesita una estimación inicial, pero sí muchos puntos para obtener una buena estimación.

Aquí hay un resumen

http://www.statsci.org/other/prony.html

En Octave, esto se implementa como expfit , por lo que puede escribir su propia rutina basada en la función de biblioteca Octave.

La estimación de Prony necesita que se conozca el desplazamiento, pero si vas "lo suficientemente lejos" en tu decadencia, tienes una estimación razonable del desplazamiento, por lo que simplemente puedes cambiar los datos para colocar el desplazamiento en 0. En cualquier caso, Prony la estimación es solo una forma de obtener una estimación inicial razonable para otras rutinas de ajuste.

No sé python, pero sí conozco una manera simple de estimar de manera no iterativa los coeficientes de decaimiento exponencial con un desplazamiento, dados tres puntos de datos con una diferencia fija en su coordenada independiente. Sus puntos de datos tienen una diferencia fija en sus coordenadas independientes (sus valores de x están espaciados en un intervalo de 60), por lo que mi método se puede aplicar a ellos. Seguramente puedes traducir las matemáticas a python.

Asumir

y = A + B*exp(-c*x) = A + B*C^x

donde C = exp(-c)

Dado y_0, y_1, y_2, para x = 0, 1, 2, resolvemos

y_0 = A + B y_1 = A + B*C y_2 = A + B*C^2

para encontrar A, B, C como sigue:

A = (y_0*y_2 - y_1^2)/(y_0 + y_2 - 2*y_1) B = (y_1 - y_0)^2/(y_0 + y_2 - 2*y_1) C = (y_2 - y_1)/(y_1 - y_0)

El exponencial correspondiente pasa exactamente a través de los tres puntos (0, y_0), (1, y_1) y (2, y_2). Si sus puntos de datos no están en las coordenadas x 0, 1, 2 sino en k, k + s, y k + 2 * s, entonces

y = A′ + B′*C′^(k + s*x) = A′ + B′*C′^k*(C′^s)^x = A + B*C^x

por lo que puede utilizar las fórmulas anteriores para encontrar A, B, C y luego calcular

A′ = A C′ = C^(1/s) B′ = B/(C′^k)

Los coeficientes resultantes son muy sensibles a los errores en las coordenadas y, lo que puede llevar a grandes errores si extrapola más allá del rango definido por los tres puntos de datos utilizados, por lo que es mejor calcular A, B, C a partir de tres puntos de datos que son tan lejos como sea posible (sin dejar de tener una distancia fija entre ellos).

Su conjunto de datos tiene 10 puntos de datos equidistantes. Seleccionemos los tres puntos de datos (110, 2391), (350, 786), (590, 263) para su uso, estos tienen la mayor distancia fija posible (240) en la coordenada independiente. Entonces, y_0 = 2391, y_1 = 786, y_2 = 263, k = 110, s = 240. Entonces A = 10.20055, B = 2380.799, C = 0.3258567, A ′ = 10.20055, B ′ = 3980.329, C ′ = 0.9953388. El exponencial es

y = 10.20055 + 3980.329*0.9953388^x = 10.20055 + 3980.329*exp(-0.004672073*x)

Puede usar esta exponencial como la conjetura inicial en un algoritmo de ajuste no lineal.

La fórmula para calcular A es la misma que la utilizada por la transformación de Shanks ( http://en.wikipedia.org/wiki/Shanks_transformation ).

Nunca conseguí que curve_fit funcione correctamente, como dices, no quiero adivinar nada. Estaba tratando de simplificar el ejemplo de Joe Kington y esto es lo que conseguí trabajando. La idea es traducir los datos ''ruidosos'' en el registro y luego volver a traducirlos y usar polyfit y polyval para determinar los parámetros:

model = np.polyfit(xVals, np.log(yVals) , 1); splineYs = np.exp(np.polyval(model,xVals[0])); pyplot.plot(xVals,yVals,'',''); #show scatter plot of original data pyplot.plot(xVals,splineYs(''b-''); #show fitted line pyplot.show()

donde xVals y yVals son solo listas.

Procedimiento para ajustar exponencial sin adivinación inicial ni proceso iterativo:

Esto viene del artículo (pp.16-17): https://fr.scribd.com/doc/14674814/Regressions-et-equations-integrales

Si es necesario, se puede utilizar para iniciar un cálculo de regresión no lineal para elegir un criterio específico de optimización.

EJEMPLO:

El ejemplo dado por Joe Kington es interesante. Desafortunadamente no se muestran los datos, solo el gráfico. Por lo tanto, los datos (x, y) a continuación provienen de una exploración gráfica de la gráfica y, como consecuencia, los valores numéricos probablemente no son exactamente los que usa Joe Kington. Sin embargo, las respectivas ecuaciones de las curvas "ajustadas" están muy cerca una de la otra, considerando la amplia dispersión de los puntos.

La figura superior es la copia de la gráfica de Kington.

La figura inferior muestra los resultados obtenidos con el procedimiento presentado anteriormente.

Si tu decaimiento no comienza desde 0 usa:

popt, pcov = curve_fit(self.func, x-x0, y)

donde x0 es el inicio de la descomposición (donde desea iniciar el ajuste). Y luego otra vez use x0 para trazar:

plt.plot(x, self.func(x-x0, *popt),''--r'', label=''Fit'')

donde la función es:

def func(self, x, a, tau, c): return a * np.exp(-x/tau) + c

Tienes dos opciones:

  1. Linealice el sistema y ajuste una línea al registro de los datos.
  2. Utilice un solucionador no lineal (por ejemplo, scipy.optimize.curve_fit

La primera opción es, con mucho, la más rápida y robusta. Sin embargo, requiere que se conozca el desplazamiento y a priori, de lo contrario es imposible linealizar la ecuación. (es decir, y = A * exp(K * t) puede ser linealizado ajustando y = log(A * exp(K * t)) = K * t + log(A) , pero y = A*exp(K*t) + C solo se puede linealizar ajustando y - C = K*t + log(A) , y como y es su variable independiente, C debe ser conocido de antemano para que sea un sistema lineal.

Si usa un método no lineal, a) no se garantiza que converja y produzca una solución, b) será mucho más lento, c) proporcione una estimación mucho menor de la incertidumbre en sus parámetros, yd) a menudo es mucho menos preciso . Sin embargo, un método no lineal tiene una gran ventaja sobre una inversión lineal: puede resolver un sistema no lineal de ecuaciones. En su caso, esto significa que no tiene que conocer C antemano.

Solo para dar un ejemplo, resolvamos para y = A * exp (K * t) con algunos datos ruidosos usando métodos tanto lineales como no lineales:

import numpy as np import matplotlib.pyplot as plt import scipy as sp import scipy.optimize def main(): # Actual parameters A0, K0, C0 = 2.5, -4.0, 2.0 # Generate some data based on these tmin, tmax = 0, 0.5 num = 20 t = np.linspace(tmin, tmax, num) y = model_func(t, A0, K0, C0) # Add noise noisy_y = y + 0.5 * (np.random.random(num) - 0.5) fig = plt.figure() ax1 = fig.add_subplot(2,1,1) ax2 = fig.add_subplot(2,1,2) # Non-linear Fit A, K, C = fit_exp_nonlinear(t, noisy_y) fit_y = model_func(t, A, K, C) plot(ax1, t, y, noisy_y, fit_y, (A0, K0, C0), (A, K, C0)) ax1.set_title(''Non-linear Fit'') # Linear Fit (Note that we have to provide the y-offset ("C") value!! A, K = fit_exp_linear(t, y, C0) fit_y = model_func(t, A, K, C0) plot(ax2, t, y, noisy_y, fit_y, (A0, K0, C0), (A, K, 0)) ax2.set_title(''Linear Fit'') plt.show() def model_func(t, A, K, C): return A * np.exp(K * t) + C def fit_exp_linear(t, y, C=0): y = y - C y = np.log(y) K, A_log = np.polyfit(t, y, 1) A = np.exp(A_log) return A, K def fit_exp_nonlinear(t, y): opt_parms, parm_cov = sp.optimize.curve_fit(model_func, t, y, maxfev=1000) A, K, C = opt_parms return A, K, C def plot(ax, t, y, noisy_y, fit_y, orig_parms, fit_parms): A0, K0, C0 = orig_parms A, K, C = fit_parms ax.plot(t, y, ''k--'', label=''Actual Function:/n $y = %0.2f e^{%0.2f t} + %0.2f$'' % (A0, K0, C0)) ax.plot(t, fit_y, ''b-'', label=''Fitted Function:/n $y = %0.2f e^{%0.2f t} + %0.2f$'' % (A, K, C)) ax.plot(t, noisy_y, ''ro'') ax.legend(bbox_to_anchor=(1.05, 1.1), fancybox=True, shadow=True) if __name__ == ''__main__'': main()

Tenga en cuenta que la solución lineal proporciona un resultado mucho más cercano a los valores reales. Sin embargo, tenemos que proporcionar el valor de desplazamiento de y para utilizar una solución lineal. La solución no lineal no requiere este conocimiento a priori.

Yo usaría la función scipy.optimize.curve_fit . La cadena de documentación para él incluso tiene un ejemplo de ajuste de un decaimiento exponencial que copiaré aquí:

>>> import numpy as np >>> from scipy.optimize import curve_fit >>> def func(x, a, b, c): ... return a*np.exp(-b*x) + c >>> x = np.linspace(0,4,50) >>> y = func(x, 2.5, 1.3, 0.5) >>> yn = y + 0.2*np.random.normal(size=len(x)) >>> popt, pcov = curve_fit(func, x, yn)

Los parámetros ajustados variarán debido al ruido aleatorio agregado, pero obtuve 2.47990495, 1.40709306, 0.53753635 como a, b y c, por lo que no es tan malo con el ruido allí. Si me adapto a y en lugar de yn, obtengo los valores exactos a, b y c.