2016-07-12 8 views
3

Erste Frage mit: I experimentelle Daten mit der Funktion der folgenden Form zu passen bin versucht:unbestimmte Anzahl von Parametern in scipy Funktion curve_fit

f(x) = m_o*(1-exp(-t_o*x)) + ... + m_j*(1-exp(-t_j*x)) 

Derzeit finde ich keine Möglichkeit, ein haben unbestimmte Anzahl von Parametern m_j, t_j, ich bin so etwas zu tun gezwungen:

def fitting_function(x, m_1, t_1, m_2, t_2): 
    return m_1*(1.-numpy.exp(-t_1*x)) + m_2*(1.-numpy.exp(-t_2*x)) 

parameters, covariance = curve_fit(fitting_function, xExp, yExp, maxfev = 100000) 

(XEXP und YEXP sind meine experimentellen Punkte)

Ist es eine Möglichkeit, meine Anpassungsfunktion so zu schreiben:

def fitting_function(x, li): 
    res = 0. 
    for el in range(len(li)/2): 
     res += li[2*idx]*(1-numpy.exp(-li[2*idx+1]*x)) 
    return res 

wo li die Liste der Anpassungsparameter ist und tun dann ein curve_fitting? Ich weiß nicht, wie man curve_fitting sagen kann, wie viele Parameter passen. Wenn ich diese Art von Formular für fitting_function versuche, habe ich Fehler wie "ValueError: Anzahl der Fit-Parameter nicht ermitteln."

Zweite Frage: Gibt es eine Möglichkeit, meine Fitting-Parameter positiv zu erzwingen?

Jede Hilfe dankbar :)

+0

Also möchten Sie curve_fit auch herausfinden, wie viele Begriffe hinzuzufügen? Oder wollen Sie, dass eine generische Summe von Exponentialfunktionen dann eine andere Anzahl von Parametern erzeugt und sieht, wann Sie eine gute Anpassung erhalten? –

+0

* "Zweite Frage: ..." * Das sollte eine separate Frage sein. Aber suchen Sie zuerst, bevor Sie es erstellen - ich denke, es wurde schon einmal gefragt. –

Antwort

1

meine Frage Sehen und here beantworten. Ich habe auch ein minimales Arbeitsbeispiel gezeigt, das zeigt, wie es für Ihre Anwendung gemacht werden kann. Ich behaupte nicht, dass dies der beste Weg ist - ich durchwühle alles selbst, also werden alle Kritiken oder Vereinfachungen geschätzt.

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

def wrapper(x, *args): #take a list of arguments and break it down into two lists for the fit function to understand 
    N = len(args)/2 
    amplitudes = list(args[0:N]) 
    timeconstants = list(args[N:2*N]) 
    return fit_func(x, amplitudes, timeconstants) 


def fit_func(x, amplitudes, timeconstants): #the actual fit function 
    fit = np.zeros(len(x)) 
    for m,t in zip(amplitudes, timeconstants): 
     fit += m*(1.0-np.exp(-t*x)) 
    return fit 

def gen_data(x, amplitudes, timeconstants, noise=0.1): #generate some fake data 
    y = np.zeros(len(x)) 
    for m,t in zip(amplitudes, timeconstants): 
     y += m*(1.0-np.exp(-t*x)) 
    if noise: 
     y += np.random.normal(0, noise, size=len(x)) 
    return y 


def main(): 
    x = np.arange(0,100) 
    amplitudes = [1, 2, 3] 
    timeconstants = [0.5, 0.2, 0.1] 
    y = gen_data(x, amplitudes, timeconstants, noise=0.01) 

    p0 = [1, 2, 3, 0.5, 0.2, 0.1] 
    popt, pcov = curve_fit(lambda x, *p0: wrapper(x, *p0), x, y, p0=p0) #call with lambda function 
    yfit = gen_data(x, popt[0:3], popt[3:6], noise=0) 
    pl.plot(x,y,x,yfit) 
    pl.show() 
    print popt 
    print pcov 

if __name__=="__main__": 
    main() 

Ein Wort der Warnung, obwohl. Eine lineare Summe von Exponentialen macht die Anpassung EXTREM empfindlich für jegliches Rauschen, insbesondere für eine große Anzahl von Parametern. Sie können dies testen, indem Sie den im Skript erzeugten Daten sogar eine kleine Menge Rauschen hinzufügen - selbst kleine Abweichungen führen dazu, dass sie die falsche Antwort vollständig erhalten, während die Anpassung immer noch mit dem Auge einwandfrei aussieht (Test mit Rauschen = 0, 0.01 und 0,1). Seien Sie sehr vorsichtig bei der Interpretation Ihrer Ergebnisse, auch wenn die Passform gut aussieht. Es ist auch eine Form, die Variablenwechsel ermöglicht: Die am besten geeignete Lösung ist dieselbe, auch wenn Sie Paare von (m_i, t_i) gegen (m_j, t_j) tauschen, was bedeutet, dass Ihr Chi-Quadrat mehrere identische lokale Minima hat Variablen werden während der Anpassung abhängig von Ihren Anfangsbedingungen ausgetauscht. Dies ist wahrscheinlich kein numerisch robuster Weg, um diese Parameter zu extrahieren.

Um Ihre zweite Frage, ja, können Sie durch die Definition Ihrer Exponentialgrößen wie so:

m_0**2*(1.0-np.exp(-t_0**2*x)+...

Grundsätzlich quadrieren sie alle in Ihrer aktuellen Anpassungsfunktion, passen sie, und dann die Ergebnisse Quadrat (was negativ oder positiv sein könnte), um Ihre tatsächlichen Parameter zu erhalten. Sie können auch Variablen definieren, die zwischen einem bestimmten Bereich liegen, indem Sie verschiedene Proxy-Formulare verwenden.

Verwandte Themen