Anpassen einer ODE mit Python lostsq gibt einen Umwandlungsfehler, wenn die Anfangsbedingungen als Parameter übergeben werden

Question

Ich habe eine Reihe von Daten, die ich an ein ODE-Modell mit scipy's Lostsq-Funktion anpassen möchte. Meine ODE hat Parameter beta und gamma, so dass es zum Beispiel wie folgt aussieht:

# dS/dt = -betaSI 
# dI/dt = betaSI - gammaI 
# dR/dt = gammaI 
# with y0 = y(t=0) = (S(0),I(0),R(0)) 

Die Idee ist beta und gamma so dass die numerische Integration meines Systems von ODE ist zu finden am besten die Daten annähert. Ich bin in der Lage, dies gut zu machen, indem ich lostsq benutze, wenn ich alle Punkte in meinem Anfangszustand kenne y0.

Nun versuche ich das gleiche zu tun, aber jetzt einen der Einträge von y0 als einen zusätzlichen Parameter zu übergeben. Hier ist, wo der Python und ich die Kommunikation stoppen ... habe ich eine Funktion, so dass nun der erste Eintrag des Parameters, die ich zu leastsq passiere die Anfangsbedingung meiner Variable R. ist ich die folgende Meldung:

*Traceback (most recent call last): 
    File "/Users/Laura/Dropbox/SHIV/shivmodels/test.py", line 73, in <module> 
    p1,success = optimize.leastsq(errfunc, initguess, args=(simpleSIR,[y0[0]],[Tx],[mydata])) 
    File "/Library/Frameworks/Python.framework/Versions/7.2/lib/python2.7/site-packages/scipy/optimize/minpack.py", line 283, in leastsq 
    gtol, maxfev, epsfcn, factor, diag) 
TypeError: array cannot be safely cast to required type* 

Hier ist mein Code. Es ist ein bisschen komplizierter, was es für dieses Beispiel sein muss, weil ich in Wirklichkeit eine andere Ode mit 7 Parametern anpassen möchte und gleichzeitig mehrere Datensätze bearbeiten möchte. Aber ich wollte hier etwas einfacheres posten ... Jede Hilfe wird sehr sehr geschätzt! Vielen Dank!

import numpy as np 
from matplotlib import pyplot as plt 
from scipy import optimize 
from scipy.integrate import odeint 

#define the time span for the ODE integration: 
Tx = np.arange(0,50,1) 
num_points = len(Tx) 

#define a simple ODE to fit: 
def simpleSIR(y,t,params): 
    dydt0 = -params[0]*y[0]*y[1] 
    dydt1 = params[0]*y[0]*y[1] - params[1]*y[1] 
    dydt2 = params[1]*y[1] 
    dydt = [dydt0,dydt1,dydt2] 
    return dydt 


#generate noisy data: 
y0 = [1000.,1.,0.] 
beta = 12*0.06/1000.0 
gamma = 0.25 
myparam = [beta,gamma] 
sir = odeint(simpleSIR, y0, Tx, (myparam,)) 

mydata0 = sir[:,0] + 0.05*(-1)**(np.random.randint(num_points,size=num_points))*sir[:,0] 
mydata1 = sir[:,1] + 0.05*(-1)**(np.random.randint(num_points,size=num_points))*sir[:,1] 
mydata2 = sir[:,2] + 0.05*(-1)**(np.random.randint(num_points,size=num_points))*sir[:,2] 
mydata = np.array([mydata0,mydata1,mydata2]).transpose() 


#define a function that will run the ode and fit it, the reason I am doing this 
#is because I will use several ODE's to see which one fits the data the best. 
def fitfunc(myfun,y0,Tx,params): 
    myfit = odeint(myfun, y0, Tx, args=(params,)) 
    return myfit 

#define a function that will measure the error between the fit and the real data: 
def errfunc(params,myfun,y0,Tx,y): 
    """ 
    INPUTS: 
    params are the parameters for the ODE 
    myfun is the function to be integrated by odeint 
    y0 vector of initial conditions, so that y(t0) = y0 
    Tx is the vector over which integration occurs, since I have several data sets and each 
    one has its own vector of time points, Tx is a list of arrays. 
    y is the data, it is a list of arrays since I want to fit to multiple data sets at once 
    """ 
    res = [] 
    for i in range(len(y)): 
     V0 = params[0][i] 
     myparams = params[1:] 
     initCond = np.zeros([3,]) 
     initCond[:2] = y0[i] 
     initCond[2] = V0 
     myfit = fitfunc(myfun,initCond,Tx[i],myparams) 
     res.append(myfit[:,0] - y[i][:,0]) 
     res.append(myfit[:,1] - y[i][:,1]) 
     res.append(myfit[1:,2] - y[i][1:,2]) 
    #end for 
    all_residuals = np.hstack(res).ravel() 
    return all_residuals 
#end errfunc 


#example of the problem: 

V0 = [0] 
params = [V0,beta,gamma] 
y0 = [1000,1] 

#this is just to test that my errfunc does work well. 
errfunc(params,simpleSIR,[y0],[Tx],[mydata]) 
initguess = [V0,0.5,0.5] 

p1,success = optimize.leastsq(errfunc, initguess, args=(simpleSIR,[y0[0]],[Tx],[mydata])) 

Answer 1

Das Problem ist mit den Variablen initigess. Die Funktion hat den folgenden optimize.leastsq Anruf Signatur:

http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.leastsq.html

Es ist zweites Argument, X0, weist ein Array sein. Ihre Liste

initguess = [v0,0.5,0.5] 

wird nicht in ein Array konvertiert, da v0 eine Liste anstelle von int oder float ist. Sie erhalten also eine Fehlermeldung, wenn Sie versuchen, initguess von einer Liste in ein Array in der Funktion "lostsq" zu konvertieren.

würde ich die Variable params von

def errfunc(params,myfun,y0,Tx,y): 

einzustellen, so dass es sich um ein 1-D-Array ist. Machen Sie die ersten paar Einträge mit den Werten von v0 und fügen Sie dann Beta und Gamma hinzu.