Comment faire un ajustement de courbe exponentielle et logarithmique en Python? Je n'ai trouvé que l'ajustement polynomial

J'ai un ensemble de données et je veux comparer quelle ligne le décrit le mieux (polynômes d'ordres différents, exponentiels ou logarithmiques).

J'utilise Python et Numpy et pour l'ajustement polynomial, il y a une fonction polyfit(). Mais je n'ai trouvé aucune fonction de ce type pour l'ajustement exponentiel et logarithmique.

Y a-t-il? Ou comment le résoudre autrement?

Pour ajuster y = A + B log x , ajustez simplement y contre (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

Pour ajuster y = Ae ^Bx , prendre le logarithme des deux côtés donne log y = log A + Bx . Alors ajustez (log y ) contre x .

Notez que l'ajustement (log y ) comme s'il était linéaire accentuera les petites valeurs de y , provoquant un grand écart pour un y grand . En effet, la polyfit(régression linéaire) fonctionne en minimisant ∑ _i (Δ Y ) ² = ∑ _i ( Y _i - Ŷ _i ) ² . Lorsque Y _i = log y _i , les résidus Δ Y _i = Δ (log y _i ) ≈ Δ y _i / | y _i |. Donc même sipolyfitprend une très mauvaise décision pour un grand y , le "divide-by- | y |" facteur compensera cela, ce qui polyfitfavorise les petites valeurs.

Cela pourrait être atténué en donnant à chaque entrée un "poids" proportionnel à y . polyfitprend en charge les moindres carrés pondérés via l' wargument mot - clé.

>>> 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)

Notez qu'Excel, LibreOffice et la plupart des calculatrices scientifiques utilisent généralement la formule non pondérée (biaisée) pour les lignes de régression / tendance exponentielles. Si vous souhaitez que vos résultats soient compatibles avec ces plates-formes, n'incluez pas les poids même si cela donne de meilleurs résultats.

Maintenant, si vous pouvez utiliser scipy, vous pouvez utiliser scipy.optimize.curve_fitpour adapter n'importe quel modèle sans transformations.

Pour y = A + B log x, le résultat est le même que la méthode de transformation:

>>> 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)

Pour y = Ae ^Bx , cependant, nous pouvons obtenir un meilleur ajustement car il calcule Δ (log y ) directement. Mais nous devons fournir une estimation d'initialisation pour curve_fitpouvoir atteindre le minimum local souhaité.

>>> 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.

Vous pouvez également adapter un ensemble de données que vous fonctionne comme l' utilisation curve_fitde scipy.optimize. Par exemple si vous souhaitez adapter une fonction exponentielle (à partir de la documentation ):

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)

Et puis si vous voulez tracer, vous pouvez faire:

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()

(Note: le *devant poptlorsque vous tracer étendra les termes dans le a, bet cque func. ATTEND)

J'avais des problèmes avec cela, alors laissez-moi être très explicite pour que les noobs comme moi puissent comprendre.

Disons que nous avons un fichier de données ou quelque chose comme ça

# -*- 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()

le résultat est: a = 0,849195983017, b = -1,18101681765, c = 2,24061176543, d = 0,816643894816

Eh bien, je suppose que vous pouvez toujours utiliser:

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

Modification légèrement de la réponse d'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()

Cela donne le graphique suivant:

Voici une option de linéarisation sur des données simples qui utilise des outils de scikit learn .

Donné

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)

Code

Ajuster les données exponentielles

# 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")

Ajuster les données du journal

# 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")

Détails

Étapes générales

1. Appliquer une opération de journal des valeurs de données ( x, you les deux)
2. Régresser les données dans un modèle linéarisé
3. Tracez en "inversant" toutes les opérations du journal (avec np.exp()) et ajustez les données d'origine

En supposant que nos données suivent une tendance exponentielle, une équation générale ⁺ peut être:

On peut linéariser cette dernière équation (par exemple y = intersection + pente * x) en prenant le log :

Étant donné une équation linéarisée ⁺⁺ et les paramètres de régression, nous pourrions calculer:

• Avia interception ( ln(A))
• Bpar pente ( B)

Résumé des techniques de linéarisation

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)

_{⁺ Remarque: la linéarisation des fonctions exponentielles fonctionne mieux lorsque le bruit est petit et C = 0. Utiliser avec précaution.}

_{⁺⁺ Remarque: alors que la modification des données x aide à linéariser les données exponentielles , la modification des données y aide à linéariser les données du journal .}

Nous démontrons les fonctionnalités de lmfittout en résolvant les deux problèmes.

Donné

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

Code

Approche 1 - lmfitModèle

Ajuster les données exponentielles

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()

Approche 2 - Modèle personnalisé

Ajuster les données du journal

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()

Détails

1. Choisissez une classe de régression
2. Fournir nommé, estimations initiales qui respectent le domaine de la fonction

Vous pouvez déterminer les paramètres déduits à partir de l'objet régresseur. Exemple:

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

Remarque: le ExponentialModel()suit une fonction de décroissance , qui accepte deux paramètres, dont l'un est négatif.

Voir aussi ExponentialGaussianModel(), qui accepte plus de paramètres .

Installez la bibliothèque via > pip install lmfit.

Wolfram a une solution de forme fermée pour ajuster une exponentielle . Ils ont également des solutions similaires pour ajuster une loi logarithmique et de puissance .

J'ai trouvé que cela fonctionnait mieux que le curve_fit de scipy. Voici un exemple:

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()