¿Cómo hacer un ajuste de curva exponencial y logarítmica en Python? Encontré solo ajuste polinómico

157

Tengo un conjunto de datos y quiero comparar qué línea lo describe mejor (polinomios de diferentes órdenes, exponenciales o logarítmicos).

Utilizo Python y Numpy y para el ajuste polinómico hay una función polyfit(). Pero no encontré tales funciones para el ajuste exponencial y logarítmico.

¿Hay alguna? ¿O cómo resolverlo de otra manera?

Tomás Novotny
fuente

Respuestas:

222

Para ajustar y = A + B log x , simplemente ajuste y contra (log x ).

>>> x = numpy.array([1, 7, 20, 50, 79])
>>> y = numpy.array([10, 19, 30, 35, 51])
>>> numpy.polyfit(numpy.log(x), y, 1)
array([ 8.46295607,  6.61867463])
# y ≈ 8.46 log(x) + 6.62

Para ajustar y = Ae Bx , tomar el logaritmo de ambos lados da log y = log A + Bx . Entonces ajuste (log y ) contra x .

Tenga en cuenta que el ajuste (log y ) como si fuera lineal enfatizará los valores pequeños de y , causando una gran desviación para y grande . Esto se debe a que polyfit(regresión lineal) funciona minimizando ∑ iY ) 2 = ∑ i ( Y i - Ŷ i ) 2 . Cuando Y i = log y i , los residuos Δ Y i = Δ (log y i ) ≈ Δ y i / | y yo |. Entonces, incluso sipolyfittoma una muy mala decisión para la gran y , la "división por- | y |" factor lo compensará, lo que polyfitfavorecerá los valores pequeños.

Esto podría aliviarse dando a cada entrada un "peso" proporcional a y . polyfitadmite mínimos cuadrados ponderados a través del wargumento de palabra clave.

>>> x = numpy.array([10, 19, 30, 35, 51])
>>> y = numpy.array([1, 7, 20, 50, 79])
>>> numpy.polyfit(x, numpy.log(y), 1)
array([ 0.10502711, -0.40116352])
#    y ≈ exp(-0.401) * exp(0.105 * x) = 0.670 * exp(0.105 * x)
# (^ biased towards small values)
>>> numpy.polyfit(x, numpy.log(y), 1, w=numpy.sqrt(y))
array([ 0.06009446,  1.41648096])
#    y ≈ exp(1.42) * exp(0.0601 * x) = 4.12 * exp(0.0601 * x)
# (^ not so biased)

Tenga en cuenta que Excel, LibreOffice y la mayoría de las calculadoras científicas suelen utilizar la fórmula no ponderada (sesgada) para las líneas de tendencia / regresión exponencial. Si desea que sus resultados sean compatibles con estas plataformas, no incluya los pesos incluso si proporciona mejores resultados.


Ahora, si puede usar scipy, podría usarlo scipy.optimize.curve_fitpara ajustar cualquier modelo sin transformaciones.

Para y = A + B log x el resultado es el mismo que el método de transformación:

>>> x = numpy.array([1, 7, 20, 50, 79])
>>> y = numpy.array([10, 19, 30, 35, 51])
>>> scipy.optimize.curve_fit(lambda t,a,b: a+b*numpy.log(t),  x,  y)
(array([ 6.61867467,  8.46295606]), 
 array([[ 28.15948002,  -7.89609542],
        [ -7.89609542,   2.9857172 ]]))
# y ≈ 6.62 + 8.46 log(x)

Sin embargo, para y = Ae Bx , podemos obtener un mejor ajuste ya que calcula Δ (log y ) directamente. Pero necesitamos proporcionar una conjetura de inicialización para que curve_fitpodamos alcanzar el mínimo local deseado.

>>> x = numpy.array([10, 19, 30, 35, 51])
>>> y = numpy.array([1, 7, 20, 50, 79])
>>> scipy.optimize.curve_fit(lambda t,a,b: a*numpy.exp(b*t),  x,  y)
(array([  5.60728326e-21,   9.99993501e-01]),
 array([[  4.14809412e-27,  -1.45078961e-08],
        [ -1.45078961e-08,   5.07411462e+10]]))
# oops, definitely wrong.
>>> scipy.optimize.curve_fit(lambda t,a,b: a*numpy.exp(b*t),  x,  y,  p0=(4, 0.1))
(array([ 4.88003249,  0.05531256]),
 array([[  1.01261314e+01,  -4.31940132e-02],
        [ -4.31940132e-02,   1.91188656e-04]]))
# y ≈ 4.88 exp(0.0553 x). much better.

comparación de regresión exponencial

kennytm
fuente
2
@Tomas: Correcto. Cambiar la base del registro simplemente multiplica una constante para registrar x o log y, lo que no afecta a r ^ 2.
kennytm
44
Esto dará mayor peso a los valores en y pequeño. Por lo tanto, es mejor ponderar las contribuciones a los valores de chi-cuadrado por y_i
Rupert Nash
17
Esta solución es incorrecta en el sentido tradicional de ajuste de curvas. No minimizará el cuadrado sumado de los residuos en el espacio lineal, sino en el espacio logarítmico. Como se mencionó anteriormente, esto cambia efectivamente la ponderación de los puntos: las observaciones donde ysea ​​pequeño se sobreponderarán artificialmente. Es mejor definir la función (lineal, no la transformación logarítmica) y usar un ajustador o minimizador de curvas.
santon
3
@santon Abordó el sesgo en la regresión exponencial.
kennytm
2
Gracias por agregar el peso! Mucha / la mayoría de la gente no sabe que puede obtener resultados cómicamente malos si intenta tomar un registro (datos) y ejecutar una línea a través de él (como Excel). Como lo había estado haciendo durante años. Cuando mi maestra bayesiana me mostró esto, pensé "¿Pero no enseñan la forma [incorrecta] en física?" - "Sí, lo llamamos 'física del bebé', es una simplificación. Esta es la forma correcta de hacerlo".
DeusXMachina
102

También puede adaptarse a un conjunto de datos a cualquier función te gusta usar curve_fita partir scipy.optimize. Por ejemplo, si desea ajustar una función exponencial (de la documentación ):

import numpy as np
import matplotlib.pyplot as plt
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)

Y luego, si quieres trazar, puedes hacer:

plt.figure()
plt.plot(x, yn, 'ko', label="Original Noised Data")
plt.plot(x, func(x, *popt), 'r-', label="Fitted Curve")
plt.legend()
plt.show()

(Nota: la *frente poptal trazar ampliará los términos en el a, by cque func. Espera)

IanVS
fuente
2
Agradable. ¿Hay alguna manera de verificar qué tan bien nos quedamos? Valor R cuadrado? ¿Existen diferentes parámetros del algoritmo de optimización que puede intentar para obtener una solución mejor (o más rápida)?
user391339
Para un mejor ajuste, puede lanzar los parámetros optimizados ajustados a la función de optimización de scipy chisquare; devuelve 2 valores, el segundo de los cuales es el valor p.
Alguna idea sobre cómo seleccionar los parámetros a, by c?
Te dije que el
47

Estaba teniendo algunos problemas con esto, así que déjame ser muy explícito para que novatos como yo puedan entender.

Digamos que tenemos un archivo de datos o algo así

# -*- coding: utf-8 -*-

import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
import numpy as np
import sympy as sym

"""
Generate some data, let's imagine that you already have this. 
"""
x = np.linspace(0, 3, 50)
y = np.exp(x)

"""
Plot your data
"""
plt.plot(x, y, 'ro',label="Original Data")

"""
brutal force to avoid errors
"""    
x = np.array(x, dtype=float) #transform your data in a numpy array of floats 
y = np.array(y, dtype=float) #so the curve_fit can work

"""
create a function to fit with your data. a, b, c and d are the coefficients
that curve_fit will calculate for you. 
In this part you need to guess and/or use mathematical knowledge to find
a function that resembles your data
"""
def func(x, a, b, c, d):
    return a*x**3 + b*x**2 +c*x + d

"""
make the curve_fit
"""
popt, pcov = curve_fit(func, x, y)

"""
The result is:
popt[0] = a , popt[1] = b, popt[2] = c and popt[3] = d of the function,
so f(x) = popt[0]*x**3 + popt[1]*x**2 + popt[2]*x + popt[3].
"""
print "a = %s , b = %s, c = %s, d = %s" % (popt[0], popt[1], popt[2], popt[3])

"""
Use sympy to generate the latex sintax of the function
"""
xs = sym.Symbol('\lambda')    
tex = sym.latex(func(xs,*popt)).replace('$', '')
plt.title(r'$f(\lambda)= %s$' %(tex),fontsize=16)

"""
Print the coefficients and plot the funcion.
"""

plt.plot(x, func(x, *popt), label="Fitted Curve") #same as line above \/
#plt.plot(x, popt[0]*x**3 + popt[1]*x**2 + popt[2]*x + popt[3], label="Fitted Curve") 

plt.legend(loc='upper left')
plt.show()

el resultado es: a = 0.849195983017, b = -1.18101681765, c = 2.24061176543, d = 0.816643894816

Datos sin procesar y función ajustada

Leandro
fuente
8
y = [np.exp(i) for i in x]es muy lento; Una de las razones por las que se creó Numpy fue para poder escribir y=np.exp(x). Además, con ese reemplazo, puedes deshacerte de tu sección de fuerza brutal. En ipython, existe la %timeitmagia de la cual In [27]: %timeit ylist=[exp(i) for i in x] 10000 loops, best of 3: 172 us per loop In [28]: %timeit yarr=exp(x) 100000 loops, best of 3: 2.85 us per loop
esmit
1
Gracias Esmit, tienes razón, pero la parte de fuerza brutal que todavía necesito usar cuando estoy tratando con datos de un csv, xls u otros formatos que he enfrentado usando este algoritmo. Creo que su uso solo tiene sentido cuando alguien intenta ajustar una función a partir de datos experimentales o de simulación, y en mi experiencia, estos datos siempre vienen en formatos extraños.
Leandro
3
x = np.array(x, dtype=float)debería permitirle deshacerse de la comprensión lenta de la lista.
Ajasja
8

Bueno, supongo que siempre puedes usar:

np.log   -->  natural log
np.log10 -->  base 10
np.log2  -->  base 2

Modificando ligeramente la respuesta de IanVS :

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

def func(x, a, b, c):
  #return a * np.exp(-b * x) + c
  return a * np.log(b * x) + c

x = np.linspace(1,5,50)   # changed boundary conditions to avoid division by 0
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)

plt.figure()
plt.plot(x, yn, 'ko', label="Original Noised Data")
plt.plot(x, func(x, *popt), 'r-', label="Fitted Curve")
plt.legend()
plt.show()

Esto da como resultado el siguiente gráfico:

ingrese la descripción de la imagen aquí

murphy1310
fuente
¿Hay un valor de saturación que el ajuste se aproxime? Si es así, ¿cómo puede acceder?
Ben
7

Aquí hay una opción de linealización en datos simples que utiliza herramientas de scikit learn .

Dado

import numpy as np

import matplotlib.pyplot as plt

from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import FunctionTransformer


np.random.seed(123)

# General Functions
def func_exp(x, a, b, c):
    """Return values from a general exponential function."""
    return a * np.exp(b * x) + c


def func_log(x, a, b, c):
    """Return values from a general log function."""
    return a * np.log(b * x) + c


# Helper
def generate_data(func, *args, jitter=0):
    """Return a tuple of arrays with random data along a general function."""
    xs = np.linspace(1, 5, 50)
    ys = func(xs, *args)
    noise = jitter * np.random.normal(size=len(xs)) + jitter
    xs = xs.reshape(-1, 1)                                  # xs[:, np.newaxis]
    ys = (ys + noise).reshape(-1, 1)
    return xs, ys
transformer = FunctionTransformer(np.log, validate=True)

Código

Ajustar datos exponenciales

# Data
x_samp, y_samp = generate_data(func_exp, 2.5, 1.2, 0.7, jitter=3)
y_trans = transformer.fit_transform(y_samp)             # 1

# Regression
regressor = LinearRegression()
results = regressor.fit(x_samp, y_trans)                # 2
model = results.predict
y_fit = model(x_samp)

# Visualization
plt.scatter(x_samp, y_samp)
plt.plot(x_samp, np.exp(y_fit), "k--", label="Fit")     # 3
plt.title("Exponential Fit")

ingrese la descripción de la imagen aquí

Ajustar datos de registro

# Data
x_samp, y_samp = generate_data(func_log, 2.5, 1.2, 0.7, jitter=0.15)
x_trans = transformer.fit_transform(x_samp)             # 1

# Regression
regressor = LinearRegression()
results = regressor.fit(x_trans, y_samp)                # 2
model = results.predict
y_fit = model(x_trans)

# Visualization
plt.scatter(x_samp, y_samp)
plt.plot(x_samp, y_fit, "k--", label="Fit")             # 3
plt.title("Logarithmic Fit")

ingrese la descripción de la imagen aquí


Detalles

Pasos generales

  1. Aplicar una operación de registro de valores de datos ( x, yo ambos)
  2. Regrese los datos a un modelo linealizado
  3. Trace "invirtiendo" cualquier operación de registro (con np.exp()) y ajústela a los datos originales

Suponiendo que nuestros datos siguen una tendencia exponencial, una ecuación general + puede ser:

ingrese la descripción de la imagen aquí

Podemos linealizar la última ecuación (por ejemplo, y = intercepción + pendiente * x) tomando el registro :

ingrese la descripción de la imagen aquí

Dada una ecuación linealizada ++ y los parámetros de regresión, podríamos calcular:

  • Aa través de intercept ( ln(A))
  • Bvía pendiente ( B)

Resumen de técnicas de linealización

Relationship |  Example   |     General Eqn.     |  Altered Var.  |        Linearized Eqn.  
-------------|------------|----------------------|----------------|------------------------------------------
Linear       | x          | y =     B * x    + C | -              |        y =   C    + B * x
Logarithmic  | log(x)     | y = A * log(B*x) + C | log(x)         |        y =   C    + A * (log(B) + log(x))
Exponential  | 2**x, e**x | y = A * exp(B*x) + C | log(y)         | log(y-C) = log(A) + B * x
Power        | x**2       | y =     B * x**N + C | log(x), log(y) | log(y-C) = log(B) + N * log(x)

+ Nota: la linealización de las funciones exponenciales funciona mejor cuando el ruido es pequeño y C = 0. Usar con precaución.

++ Nota: mientras que la modificación de datos x ayuda a linealizar datos exponenciales , la modificación de datos y ayuda a linealizar datos de registro .

pylang
fuente
0

Demostramos características de lmfital resolver ambos problemas.

Dado

import lmfit

import numpy as np

import matplotlib.pyplot as plt


%matplotlib inline
np.random.seed(123)

# General Functions
def func_log(x, a, b, c):
    """Return values from a general log function."""
    return a * np.log(b * x) + c


# Data
x_samp = np.linspace(1, 5, 50)
_noise = np.random.normal(size=len(x_samp), scale=0.06)
y_samp = 2.5 * np.exp(1.2 * x_samp) + 0.7 + _noise
y_samp2 = 2.5 * np.log(1.2 * x_samp) + 0.7 + _noise

Código

Enfoque 1 - lmfitModelo

Ajustar datos exponenciales

regressor = lmfit.models.ExponentialModel()                # 1    
initial_guess = dict(amplitude=1, decay=-1)                # 2
results = regressor.fit(y_samp, x=x_samp, **initial_guess)
y_fit = results.best_fit    

plt.plot(x_samp, y_samp, "o", label="Data")
plt.plot(x_samp, y_fit, "k--", label="Fit")
plt.legend()

ingrese la descripción de la imagen aquí

Enfoque 2 - Modelo personalizado

Ajustar datos de registro

regressor = lmfit.Model(func_log)                          # 1
initial_guess = dict(a=1, b=.1, c=.1)                      # 2
results = regressor.fit(y_samp2, x=x_samp, **initial_guess)
y_fit = results.best_fit

plt.plot(x_samp, y_samp2, "o", label="Data")
plt.plot(x_samp, y_fit, "k--", label="Fit")
plt.legend()

ingrese la descripción de la imagen aquí


Detalles

  1. Elige una clase de regresión
  2. Suministre conjeturas iniciales con nombre que respeten el dominio de la función

Puede determinar los parámetros inferidos del objeto regresor. Ejemplo:

regressor.param_names
# ['decay', 'amplitude']

Nota: ExponentialModel()sigue una función de disminución , que acepta dos parámetros, uno de los cuales es negativo.

ingrese la descripción de la imagen aquí

Ver también ExponentialGaussianModel(), que acepta más parámetros. .

Instalar la biblioteca a través de > pip install lmfit.

pylang
fuente
0

Wolfram tiene una solución de forma cerrada para ajustar una exponencial . También tienen soluciones similares para ajustar una ley logarítmica y de poder .

Encontré que esto funciona mejor que scipy's curve_fit. Aquí hay un ejemplo:

import numpy as np
import matplotlib.pyplot as plt

# Fit the function y = A * exp(B * x) to the data
# returns (A, B)
# From: https://mathworld.wolfram.com/LeastSquaresFittingExponential.html
def fit_exp(xs, ys):
    S_x2_y = 0.0
    S_y_lny = 0.0
    S_x_y = 0.0
    S_x_y_lny = 0.0
    S_y = 0.0
    for (x,y) in zip(xs, ys):
        S_x2_y += x * x * y
        S_y_lny += y * np.log(y)
        S_x_y += x * y
        S_x_y_lny += x * y * np.log(y)
        S_y += y
    #end
    a = (S_x2_y * S_y_lny - S_x_y * S_x_y_lny) / (S_y * S_x2_y - S_x_y * S_x_y)
    b = (S_y * S_x_y_lny - S_x_y * S_y_lny) / (S_y * S_x2_y - S_x_y * S_x_y)
    return (np.exp(a), b)


xs = [33, 34, 35, 36, 37, 38, 39, 40, 41, 42]
ys = [3187, 3545, 4045, 4447, 4872, 5660, 5983, 6254, 6681, 7206]

(A, B) = fit_exp(xs, ys)

plt.figure()
plt.plot(xs, ys, 'o-', label='Raw Data')
plt.plot(xs, [A * np.exp(B *x) for x in xs], 'o-', label='Fit')

plt.title('Exponential Fit Test')
plt.xlabel('X')
plt.ylabel('Y')
plt.legend(loc='best')
plt.tight_layout()
plt.show()

ingrese la descripción de la imagen aquí

Ben
fuente