dati statistiche con funzione integrale

Question

Quando si utilizza curve_fit da scipy.optimize per adattare un alcuni dati in pitone, un primo definisce la funzione raccordo (ad esempio un 2 ° ordine polinomiale) come segue:

1. def f(x, a, b): return a*x**2+b*x
2. e poi procede con il raccordo popt, pcov = curve_fit(f,x,y)

Ma la domanda è ora, come si fa a definire la funzione nel punto 1. se la funzione contiene un integrale (o una somma discreta), ad esempio:

I dati sperimentali è ancora dato per x e f (x), quindi punto 2. sarebbe simile immagino volta posso definire f (x) in pitone. Tra l'altro ho dimenticato di dire che si presume che g (t) abbia una forma ben conosciuta qui e contiene i parametri di adattamento, cioè i parametri come aeb dati nell'esempio polinomiale. Ogni aiuto è molto apprezzato. La domanda dovrebbe essere generica e le funzioni utilizzate nel post sono solo esempi casuali.

Answer 1

Ecco un esempio di adattamento di una curva definita in termini di integrale. La curva è l'integrale di sin(t*w)/t+p su t da 0 a Pi. I nostri punti dati x corrispondono a w e stiamo modificando il parametro p per ottenere i dati adatti.

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 

che ti stampare qualcosa abbastanza vicino a 1,0, che è quello che abbiamo usato come p quando abbiamo creato il trueydata.

Si noti che utilizziamo numpy.vectorize sulla funzione curva per produrre una versione vettoriale compatibile con scipy.optimize.curve_fit.

Answer 2

A volte si può essere fortunati e si è in grado di valutare l'integrale analiticamente. Nell'esempio seguente il prodotto di h(t)=exp(-(t-x)**2/2) e un polinomio di secondo grado g(t) sono integrati da 0 a infinito. Sympy viene utilizzato per valutare l'Integrale e generare una funzione utilizzabile per 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 

Si noti che, in generale, il problema è spesso mal posto dal momento che l'integrale non neccesarily convergono in un grande quartiere abbastanza della soluzione (che è assunto da curve_fit()).