配件與Python leastsq一個ODE時的初始條件作爲參數

Question

過去，我有一組數據，我想適合使用SciPy的的leastsq功能的ODE模型給出鑄造錯誤。我的ODE有參數beta和gamma，所以它看起來例如像這樣：

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

的想法是找到beta和gamma讓我ODE最好的系統的數值積分近似數據。如果我知道我的初始條件y0中的所有點，我可以使用leastsq來做到這一點。

現在，我試圖做同樣的事情，但現在通過的y0的條目之一作爲一個額外的參數。這裏是Python和我停止溝通的地方... 我做了一個函數，現在我傳遞給leastsq的參數的第一項是我的變量R的初始條件。 我收到以下消息：

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

這裏是我的代碼。這是一個涉及多一點什麼它需要在這個例子中，因爲在現實中我要適應另一種頌歌7個參數，並希望以適應多個數據集一次。但我想在這裏發佈一些更簡單的...任何幫助將非常感激！非常感謝你！

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

問題是與該變量initguess。該功能optimize.leastsq具有以下調用簽名：

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

這是第二個參數，X0，必須是一個數組。您的列表

initguess = [v0,0.5,0.5] 

不會轉換爲數組，因爲v0是列表而不是int或float。所以當你嘗試從列表中轉換initgsss到leastsq函數中的數組時，你會得到一個錯誤。

我會從

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

調整可變PARAMS，以便它是一個1-d陣列。做出前幾個條目的v0的值，然後附加beta和gamma。