webentwicklung-frage-antwort-db.com.de

exponentieller Zerfall ohne anfängliche Vermutung

Kennt jemand ein Scipy/Numpy-Modul, das einen exponentiellen Zerfall der Daten ermöglicht? 

Die Google-Suche hat ein paar Blogeinträge zurückgegeben, z. B. - http://exnumerus.blogspot.com/2010/04/how-to-fit-fit-exponential-decay-example-in.html , diese Lösung erfordert jedoch y -offset ist vorzugeben, was nicht immer möglich ist 

BEARBEITEN:

curve_fit funktioniert, aber es kann ziemlich unglücklich ausfallen, ohne dass eine anfängliche Schätzung der Parameter vorliegt, was manchmal erforderlich ist. Der Code, mit dem ich arbeite, ist

#!/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()

was funktioniert, aber wenn wir "p0 = rate" entfernen, scheitert es kläglich.

20

Sie haben zwei Möglichkeiten:

  1. Linearisieren Sie das System und passen Sie eine Linie an das Protokoll der Daten an.
  2. Verwenden Sie einen nichtlinearen Solver (z. B. scipy.optimize.curve_fit

Die erste Option ist bei weitem die schnellste und robusteste. Es ist jedoch erforderlich, dass Sie den y-Versatz a-priori kennen, andernfalls kann die Gleichung nicht linearisiert werden. (Das heißt y = A * exp(K * t) kann durch Anpassen von y = log(A * exp(K * t)) = K * t + log(A) linearisiert werden, aber y = A*exp(K*t) + C kann nur durch Anpassen von y - C = K*t + log(A) linearisiert werden. Da y Ihre unabhängige Variable ist, muss C vorher bekannt sein, damit dies ein lineares System ist.

Wenn Sie eine nichtlineare Methode verwenden, ist es nicht garantiert, dass sie konvergiert und eine Lösung liefert, b) viel langsamer ist, c) die Unsicherheit in Ihren Parametern viel schlechter einschätzt und d) die Genauigkeit viel geringer ist . Ein nichtlineares Verfahren hat jedoch einen großen Vorteil gegenüber einer linearen Inversion: Es kann ein nichtlineares Gleichungssystem lösen. In Ihrem Fall bedeutet dies, dass Sie C nicht vorher kennen müssen.

Um ein Beispiel zu geben, lösen wir y = A * exp (K * t) mit einigen verrauschten Daten, wobei sowohl lineare als auch nichtlineare Methoden verwendet werden:

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 __== '__main__':
    main()

Fitting exp

Beachten Sie, dass die lineare Lösung ein Ergebnis liefert, das den tatsächlichen Werten viel näher kommt. Wir müssen jedoch den y-Versatzwert angeben, um eine lineare Lösung verwenden zu können. Die nichtlineare Lösung erfordert dieses a-priori-Wissen nicht. 

42
Joe Kington

Ich würde die scipy.optimize.curve_fit-Funktion verwenden. Der doc-String dafür hat sogar ein Beispiel für das Anpassen eines exponentiellen Zerfalls, das ich hier kopieren werde:

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

Die angepassten Parameter variieren aufgrund des zufälligen Rauschens, aber ich habe 2.47990495, 1.40709306, 0.53753635 als a, b und c erhalten, das ist also nicht so schlimm. Wenn ich zu y und nicht zu yn passe, bekomme ich die exakten Werte für a, b und c.

7
Justin Peel

Vorgehensweise, um exponentiell ohne anfängliche Vermutung und nicht iterativen Prozess zu passen:

 enter image description here

Dies kommt aus dem Papier (S.16-17): https://fr.scribd.com/doc/14674814/Regressions-et-equations-integrales

Bei Bedarf kann hiermit ein nichtlinearer Regressionskalkül initialisiert werden, um ein bestimmtes Optimierungskriterium zu wählen.

BEISPIEL

Das Beispiel von Joe Kington ist interessant. Leider werden die Daten nicht angezeigt, nur die Grafik. Die untenstehenden Daten (x, y) stammen also von einer grafischen Untersuchung des Graphen. Daher sind die numerischen Werte wahrscheinlich nicht genau die von Joe Kington verwendeten. Dennoch liegen die jeweiligen Gleichungen der "angepassten" Kurven in Anbetracht der breiten Streuung der Punkte sehr nahe beieinander.

 enter image description here

Die obere Abbildung ist die Kopie des Kington-Diagramms.

Die untere Abbildung zeigt die Ergebnisse, die mit dem oben dargestellten Verfahren erhalten wurden.

6
JJacquelin

Der richtige Weg, dies zu tun, ist, eine Prony-Schätzung durchzuführen und das Ergebnis als anfängliche Schätzung für die Anpassung der kleinsten Quadrate (oder eine andere, robustere Anpassungsroutine) zu verwenden. Bei der Schätzung von Prony ist keine erste Schätzung erforderlich, aber es sind viele Punkte erforderlich, um eine gute Schätzung zu erhalten.

Hier ist eine Übersicht

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

In Octave ist dies als expfit implementiert, sodass Sie Ihre eigene Routine basierend auf der Octave-Bibliotheksfunktion schreiben können.

Bei der Prony-Schätzung muss der Versatz bekannt sein. Wenn Sie jedoch "weit genug" in Ihren Zerfall gehen, haben Sie eine vernünftige Schätzung des Versatzes, sodass Sie die Daten einfach verschieben können, um den Versatz auf 0 zu setzen. Jedenfalls Prony Schätzung ist nur ein Weg, um eine vernünftige anfängliche Schätzung für andere Anpassungsroutinen zu erhalten. 

3
Marcus P S

Ich habe curve_fit nie dazu gebracht, richtig zu funktionieren, wie Sie sagen, ich möchte nichts erraten. Ich habe versucht, das Beispiel von Joe Kington zu vereinfachen. Die Idee ist, die "verrauschten" Daten in ein Protokoll zu übersetzen und sie dann zurückzurücken und polyfit und polyval zu verwenden, um die Parameter herauszufinden:

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

dabei sind xVals und yVals nur Listen. 

1
Elendurwen

Ich kenne Python nicht, aber ich kenne einen einfachen Weg, um die Koeffizienten des exponentiellen Abfalls nicht-iterativ mit einem Offset zu schätzen, wenn drei Datenpunkte mit einer festen Differenz in ihrer unabhängigen Koordinate angegeben werden. Ihre Datenpunkte haben eine feste Differenz in ihrer unabhängigen Koordinate (Ihre X-Werte sind in einem Abstand von 60 angeordnet), sodass meine Methode auf sie angewendet werden kann. Sie können die Mathematik sicher in Python übersetzen.

Annehmen

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

wo C = exp(-c)

Mit y_0, y_1, y_2 für x = 0, 1, 2 lösen wir 

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

a, B, C wie folgt finden:

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)

Das entsprechende Exponential verläuft genau durch die drei Punkte (0, y_0), (1, y_1) und (2, y_2). Wenn Ihre Datenpunkte nicht an den x-Koordinaten 0, 1, 2 liegen, sondern eher bei k, k + s und k + 2 * s, dann 

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

sie können also die obigen Formeln verwenden, um A, B, C zu finden und dann zu berechnen

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

Die resultierenden Koeffizienten sind sehr empfindlich für Fehler in den y-Koordinaten, was zu großen Fehlern führen kann, wenn Sie über den durch die drei verwendeten Datenpunkte definierten Bereich hinaus extrapolieren. Daher ist es am besten, A, B, C aus drei Datenpunkten zu berechnen so weit wie möglich voneinander entfernt (während immer noch ein fester Abstand zwischen ihnen besteht).

Ihr Datensatz hat 10 äquidistante Datenpunkte. Wählen wir die drei Datenpunkte (110, 2391), (350, 786), (590, 263) zur Verwendung aus - diese haben den größtmöglichen festen Abstand (240) in der unabhängigen Koordinate. Y_0 = 2391, y_1 = 786, y_2 = 263, k = 110, s = 240. Dann ist A = 10.20055, B = 2380.799, C = 0,3258567, A '= 10.20055, B' = 3980.329, C '= 0,9953388. Das Exponential ist

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

Sie können dieses Exponential als anfänglichen Schätzwert in einem nichtlinearen Anpassungsalgorithmus verwenden.

Die Formel für die Berechnung von A ist die gleiche, die von der Shanks-Transformation verwendet wird ( http://en.wikipedia.org/wiki/Shanks_transformation ).

1
Louis Strous

Python-Implementierung der @ JJacquelin-Lösung. Ich brauchte eine ungefähre nicht lösungsbasierte Lösung ohne anfängliche Vermutungen, daher war die Antwort von @JJacquelin wirklich hilfreich. Die ursprüngliche Frage wurde als Python-Anfrage von Numpy/Scipy gestellt. Ich habe @ johanvdws Nice-Clean-R-Code genommen und als Python/Numpy überarbeitet. Hoffentlich nützlich für jemanden: 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]
0
FriendToGeoff

Wenn Ihr Verfall nicht bei 0 beginnt, verwenden Sie:

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

wobei x0 der Beginn des Zerfalls ist (wo der Fit beginnen soll). Und dann erneut x0 zum Plotten:

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

wo ist die Funktion:

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