2015-05-19 20 views
10

Przy zastosowaniu curve_fit z scipy.optimize dopasowania do niektórych danych pythonowa najpierw określa funkcję dopasowania (na przykład 2-go rzędu wielomian) w następujący sposób:dopasowanie danych z wbudowaną funkcję

  1. def f(x, a, b): return a*x**2+b*x
  2. a następnie przechodzi z okucia popt, pcov = curve_fit(f,x,y)

ale pytanie jest teraz, w jaki sposób można przejść o definiowaniu funkcji w punkcie 1. jeżeli funkcja zawiera integralną (lub dyskretnych sumę), np:

enter image description here

Dane doświadczalne są jeszcze podane dla X i F (x), a więc punkt 2. byłaby podobna sobie wyobrazić po można zdefiniować f (x) w pytona. Przy okazji zapomniałem powiedzieć, że zakłada się, że g (t) ma tutaj dobrze znaną postać i zawiera parametry dopasowania, tj. Parametry takie jak a i b podane w przykładzie wielomianowym. Każda pomoc jest doceniana. Pytanie naprawdę ma być ogólne, a funkcje używane w poście są po prostu przypadkowymi przykładami.

+0

Oczywistą odpowiedzią jest: potrzebujesz sposobu na ocenę tej całki, albo poprzez znalezienie rozwiązania zamkniętego, albo za pomocą kwadratury liczbowej. Nie ma ogólnego rozwiązania tego problemu. – cfh

+0

@ cfh oh Widzę, to prawda, ale jeśli nie ma żadnego rozwiązania zamkniętego, co dokładnie oznacza kwadratura liczbowa? czy nie zakłada, że ​​wszystkie parametry powinny być wtedy znane? –

+0

Tak, ale w momencie wywołania 'f' znasz wszystkie parametry, ponieważ są one przekazywane jako argumenty. – cfh

Odpowiedz

6

Oto przykład dopasowania krzywej zdefiniowanej jako integralna. Krzywa jest całką sin(t*w)/t+p przez t od 0 do Pi. Nasze x punkty danych odpowiadają w, a my dostosowujemy parametr p, aby dopasować dane.

import math, numpy, scipy.optimize, scipy.integrate 

def integrand(t, args): 
    w, p = args 
    return math.sin(t * w)/t + p 

def curve(w, p): 
    res = scipy.integrate.quad(integrand, 0.0, math.pi, [w, p]) 
    return res[0] 

vcurve = numpy.vectorize(curve, excluded=set([1])) 

truexdata = numpy.asarray([0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0]) 
trueydata = vcurve(truexdata, 1.0) 

xdata = truexdata + 0.1 * numpy.random.randn(8) 
ydata = trueydata + 0.1 * numpy.random.randn(8) 

popt, pcov = scipy.optimize.curve_fit(vcurve, 
             xdata, ydata, 
             p0=[2.0]) 
print popt 

To będzie wydrukować coś z dość blisko do 1,0, czyli co użyliśmy jako p kiedy stworzyliśmy trueydata.

Należy pamiętać, że używamy funkcji numpy.vectorize w funkcji krzywej w celu utworzenia wektorowej wersji zgodnej z scipy.optimize.curve_fit.

+2

to jest interesujące, użyteczne i myślę, że generyczne wystarczy, ale czy możemy zrobić trochę lepiej wiedząc, że jest to integralne? W twoim rozwiązaniu "krzywa" może być wszystkim i nie korzysta z faktu, że jest integralna. Może pytanie powinno zostać zmienione na "dopasowanie danych z nie-wektoryzowanymi funkcjami"? – dashesy

3

Czasami możesz mieć szczęście i możesz ocenić całkę analitycznie. W poniższym przykładzie produkt h(t)=exp(-(t-x)**2/2) i wielomian drugiego stopnia g(t) jest zintegrowany od 0 do nieskończoności. Sympy służy do oceny Integral i wygenerować funkcję użytkową dla curve_fit():

import sympy as sy 
sy.init_printing() # LaTeX-like pretty printing of IPython 


t, x = sy.symbols("t, x", real=True) 

h = sy.exp(-(t-x)**2/2) 

a0, a1, a2 = sy.symbols('a:3', real=True) # unknown coefficients 
g = a0 + a1*t + a2*t**2 

gh = (g*h).simplify() # the intgrand 
G = sy.integrate(gh, (t, 0, sy.oo)).simplify() # integrate from 0 to infinty 

# Generate numeric function to be usable by curve_fit() 
G_opt = sy.lambdify((x, t, a0, a1, a2), G) 

print(G_opt(1, 2, 3, 4, 5)) # example usage 

Zauważ, że w ogóle problem jest często źle postawione, ponieważ całka nie neccesarily zbiegają się w wystarczająco dużej dzielnicy roztworu (co jest przyjmowane przez curve_fit()).