Wie führe ich eine exponentielle und logarithmische Kurvenanpassung in Python durch? Ich fand nur eine Polynomanpassung

157

Ich habe einen Datensatz und möchte vergleichen, welche Zeile ihn am besten beschreibt (Polynome unterschiedlicher Ordnung, exponentiell oder logarithmisch).

Ich benutze Python und Numpy und für die Polynomanpassung gibt es eine Funktion polyfit(). Aber ich habe keine solchen Funktionen für die exponentielle und logarithmische Anpassung gefunden.

Sind da irgendwelche? Oder wie kann man es anders lösen?

Tomas Novotny
quelle

Antworten:

222

Zum Anpassen von y = A + B log x passen Sie einfach y gegen (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

Nehmen Sie für die Anpassung von y = Ae Bx den Logarithmus beider Seiten und geben Sie log y = log A + Bx an . Passen Sie also (log y ) gegen x an .

Beachten Sie, dass die Anpassung (log y ) als linear kleine Werte von y hervorhebt , was zu einer großen Abweichung für großes y führt . Dies liegt daran, dass polyfit(lineare Regression) durch Minimieren von ∑ iY ) 2 = ∑ i ( Y i - Ŷ i ) 2 funktioniert . Wenn Y i = log y i ist , sind die Reste Δ Y i = Δ (log y i ) ≈ Δ y i / | y i |. Also auch wennpolyfittrifft eine sehr schlechte Entscheidung für großes y , das "dividieren durch | y |" Faktor wird dies kompensieren und polyfitkleine Werte begünstigen.

Dies könnte gemildert werden, indem jedem Eintrag ein "Gewicht" proportional zu y gegeben wird . polyfitunterstützt gewichtete kleinste Quadrate über das wSchlüsselwortargument.

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

Beachten Sie, dass Excel, LibreOffice und die meisten wissenschaftlichen Taschenrechner normalerweise die ungewichtete (voreingenommene) Formel für die exponentiellen Regressions- / Trendlinien verwenden. Wenn Sie möchten, dass Ihre Ergebnisse mit diesen Plattformen kompatibel sind, geben Sie die Gewichte nicht an, auch wenn dadurch bessere Ergebnisse erzielt werden.


Wenn Sie jetzt scipy verwenden können, können Sie scipy.optimize.curve_fitjedes Modell ohne Transformationen anpassen .

Für y = A + B log x ist das Ergebnis dasselbe wie für die Transformationsmethode:

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

Für y = Ae Bx können wir jedoch eine bessere Anpassung erhalten, da es Δ (log y ) direkt berechnet . Wir müssen jedoch eine Initialisierungsschätzung abgeben, damit curve_fitdas gewünschte lokale Minimum erreicht werden kann.

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

Vergleich der exponentiellen Regression

kennytm
quelle
2
@Tomas: Richtig. Durch Ändern der log-Basis wird lediglich eine Konstante mit log x oder log y multipliziert, was sich nicht auf r ^ 2 auswirkt.
Kennytm
4
Dies gibt Werten bei kleinem y ein größeres Gewicht. Daher ist es besser, Beiträge zu den Chi-Quadrat-Werten von y_i
Rupert Nash
17
Diese Lösung ist im herkömmlichen Sinne der Kurvenanpassung falsch. Das summierte Quadrat der Residuen im linearen Raum wird nicht minimiert, sondern im logarithmischen Raum. Wie bereits erwähnt, ändert dies effektiv die Gewichtung der Punkte - Beobachtungen, bei denen yes sich um kleine Punkte handelt, werden künstlich übergewichtet. Es ist besser, die Funktion (linear, nicht die Protokolltransformation) zu definieren und einen Kurvenmonteur oder Minimierer zu verwenden.
Santon
3
@santon Hat die Verzerrung der exponentiellen Regression behoben.
Kennytm
2
Vielen Dank für das Hinzufügen des Gewichts! Viele / die meisten Leute wissen nicht, dass Sie komisch schlechte Ergebnisse erzielen können, wenn Sie versuchen, nur ein Protokoll (Daten) zu nehmen und eine Zeile durch dieses zu führen (wie Excel). Wie ich es seit Jahren getan hatte. Als mein Bayesianischer Lehrer mir dies zeigte, sagte ich: "Aber unterrichten sie nicht den [falschen] Weg in phys?" - "Ja, wir nennen das 'Babyphysik', es ist eine Vereinfachung. Dies ist der richtige Weg, dies zu tun."
DeusXMachina
102

Sie können auch einen Datensatz an eine beliebige Funktion anpassen, die Sie verwenden curve_fitmöchten scipy.optimize. Zum Beispiel, wenn Sie eine Exponentialfunktion anpassen möchten (aus der Dokumentation ):

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)

Und wenn Sie dann zeichnen möchten, können Sie Folgendes tun:

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

(Hinweis: Die *vor , poptwenn Sie die Begriffe in die expandieren aus plot a, bund , cDas func. Erwartet)

IanVS
quelle
2
Nett. Gibt es eine Möglichkeit zu überprüfen, wie gut wir passen? R-Quadrat-Wert? Gibt es verschiedene Optimierungsalgorithmusparameter, mit denen Sie versuchen können, eine bessere (oder schnellere) Lösung zu erhalten?
user391339
Für eine gute Passform können Sie die angepassten optimierten Parameter in das Chisquare der Scipy-Optimierungsfunktion einfügen. es gibt 2 Werte zurück, von denen der 2. der p-Wert ist.
Jede Idee, wie die Parameter zu wählen a, bund c?
random9
47

Ich hatte einige Probleme damit, also lass mich sehr explizit sein, damit Noobs wie ich verstehen können.

Nehmen wir an, wir haben eine Datendatei oder ähnliches

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

Das Ergebnis ist: a = 0,849195983017, b = -1,18101681765, c = 2,24061176543, d = 0,816643894816

Rohdaten und angepasste Funktion

Leandro
quelle
8
y = [np.exp(i) for i in x]ist sehr langsam; Ein Grund, warum Numpy erstellt wurde, war, dass Sie schreiben konnten y=np.exp(x). Mit diesem Ersatz können Sie auch Ihren brutalen Kraftabschnitt loswerden. In Ipython gibt es die %timeitMagie, aus der 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
Vielen Dank, Sie haben Recht, aber der brutale Kraftteil, den ich noch verwenden muss, wenn ich mit Daten aus einer CSV-, XLS- oder anderen Formaten arbeite, mit denen ich mit diesem Algorithmus konfrontiert bin. Ich denke, dass die Verwendung nur dann sinnvoll ist, wenn jemand versucht, eine Funktion aus experimentellen oder Simulationsdaten anzupassen, und meiner Erfahrung nach kommen diese Daten immer in seltsamen Formaten vor.
Leandro
3
x = np.array(x, dtype=float)sollte es Ihnen ermöglichen, das langsame Listenverständnis loszuwerden.
Ajasja
8

Nun, ich denke, Sie können immer verwenden:

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

Die Antwort von IanVS leicht modifizieren :

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

Dies führt zu der folgenden Grafik:

Geben Sie hier die Bildbeschreibung ein

murphy1310
quelle
Gibt es einen Sättigungswert, den die Anpassung annähert? Wenn ja, wie kann man darauf zugreifen?
Ben
7

Hier ist eine Linearisierungsoption für einfache Daten, die Tools von scikit learn verwendet .

Gegeben

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

Exponentialdaten anpassen

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

Geben Sie hier die Bildbeschreibung ein

Protokolldaten anpassen

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

Geben Sie hier die Bildbeschreibung ein


Einzelheiten

Allgemeine Schritte

  1. Tragen Sie eine Log - Operation , um Datenwerte ( x, yoder beides)
  2. Regressieren Sie die Daten in ein linearisiertes Modell
  3. Zeichnen Sie, indem Sie alle Protokollvorgänge (mit np.exp()) "umkehren" und an die Originaldaten anpassen

Unter der Annahme, dass unsere Daten einem exponentiellen Trend folgen, kann eine allgemeine Gleichung + sein:

Geben Sie hier die Bildbeschreibung ein

Wir können die letztere Gleichung (z. B. y = Achsenabschnitt + Steigung * x) linearisieren, indem wir das Protokoll nehmen :

Geben Sie hier die Bildbeschreibung ein

Mit einer linearisierten Gleichung ++ und den Regressionsparametern könnten wir berechnen:

  • Avia intercept ( ln(A))
  • Büber Steigung ( B)

Zusammenfassung der Linearisierungstechniken

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)

+ Hinweis: Die Linearisierung von Exponentialfunktionen funktioniert am besten, wenn das Rauschen klein und C = 0 ist. Mit Vorsicht verwenden.

++ Hinweis: Während das Ändern von x-Daten zur Linearisierung exponentieller Daten beiträgt , hilft das Ändern von y-Daten zur Linearisierung von Protokolldaten .

Pylang
quelle
0

Wir demonstrieren Merkmale bei der lmfitLösung beider Probleme.

Gegeben

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

Ansatz 1 - lmfitModell

Exponentialdaten anpassen

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

Geben Sie hier die Bildbeschreibung ein

Ansatz 2 - Benutzerdefiniertes Modell

Protokolldaten anpassen

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

Geben Sie hier die Bildbeschreibung ein


Einzelheiten

  1. Wählen Sie eine Regressionsklasse
  2. Angebot benannt, erste Vermutungen, die die Domäne der Funktion berücksichtigen

Sie können die abgeleiteten Parameter aus dem Regressorobjekt ermitteln. Beispiel:

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

Hinweis: Es ExponentialModel()folgt eine Abklingfunktion , die zwei Parameter akzeptiert, von denen einer negativ ist.

Geben Sie hier die Bildbeschreibung ein

Siehe auch ExponentialGaussianModel(), das weitere Parameter akzeptiert .

Installieren Sie die Bibliothek über > pip install lmfit.

Pylang
quelle
0

Wolfram hat eine geschlossene Lösung zum Anpassen eines Exponentials . Sie haben auch ähnliche Lösungen für die Anpassung eines Logarithmus- und Potenzgesetzes .

Ich fand, dass dies besser funktioniert als Scipys Curve_Fit. Hier ist ein Beispiel:

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

Geben Sie hier die Bildbeschreibung ein

Ben
quelle