What I am doing: I modified the code from the zombie invasion system to demonstrate how it should be written and tried to optimize the least square error (defined as score function) with the fmin function.
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from scipy import integrate
from scipy.optimize import fmin
#=====================================================
#Notice we must import the Model Definition
from zombiewithdata import eq
#=====================================================
#1.Get Data
#====================================================
Td=np.array([0.5,1,1.5,2,2.2,3,3.5,4,4.5,5])#time
Zd=np.array([0,2,2,5,2,10,15,50,250,400])#zombie pop
#====================================================
#2.Set up Info for Model System
#===================================================
# model parameters
#----------------------------------------------------
P = 0 # birth rate
d = 0.0001 # natural death percent (per day)
B = 0.0095 # transmission percent (per day)
G = 0.0001 # resurect percent (per day)
A = 0.0001 # destroy perecent (per day)
rates=(P,d,B,G,A)
# model initial conditions
#---------------------------------------------------
S0 = 500. # initial population
Z0 = 0 # initial zombie population
R0 = 0 # initial death population
y0 = [S0, Z0, R0] # initial condition vector
# model steps
#---------------------------------------------------
start_time=0.0
end_time=5.0
intervals=1000
mt=np.linspace(start_time,end_time,intervals)
# model index to compare to data
#----------------------------------------------------
findindex=lambda x:np.where(mt>=x)[0][0]
mindex=map(findindex,Td)
#=======================================================
#3.Score Fit of System
#=========================================================
def score(parms):
#a.Get Solution to system
F0,F1,F2,T=eq(parms,y0,start_time,end_time,intervals)
#b.Pick of Model Points to Compare
Zm=F1[mindex]
#c.Score Difference between model and data points
ss=lambda data,model:((data-model)**2).sum()
return ss(Zd,Zm)
#========================================================
#4.Optimize Fit
#=======================================================
fit_score=score(rates)
answ=fmin(score,(rates),full_output=1,maxiter=1000000)
bestrates=answ[0]
bestscore=answ[1]
P,d,B,G,A=answ[0]
newrates=(P,d,B,G,A)
#=======================================================
#5.Generate Solution to System
#=======================================================
F0,F1,F2,T=eq(newrates,y0,start_time,end_time,intervals)
Zm=F1[mindex]
Tm=T[mindex]
#======================================================
Now in the #optimize fit section, is there any way I can get best possible values of bestrates when I restrict the values of "rates" like lb <= P, d, B, G, A <= ub where lb=lower bound and ub=upper bound and manage to get minimum of score in that restricted region? It need not be the most optimized value. fmin uses Nelder-Mead (simplex) algorithm.
I am quite new to this, so any help in the right direction would be awesome. Feel free to ask any doubts regarding the code and I will answer to best of my knowledge. . Thank You.
I'm not sure why the original author of Adventures in Python : Fitting a Differential Equation System to Data jumps through hoops to get at the samples corresponding to the given data points, the procedure can be greatly simplified by passing to eq a time array instead of its construction parameters
#=======================================================
def eq(par,initial_cond,t):
#differential-eq-system----------------------
def funct(y,t):
Si, Zi, Ri=y
P,d,B,G,A=par
# the model equations (see Munz et al. 2009)
f0 = P - B*Si*Zi - d*Si
f1 = B*Si*Zi + G*Ri - A*Si*Zi
f2 = d*Si + A*Si*Zi - G*Ri
return [f0, f1, f2]
#integrate------------------------------------
ds = odeint(funct,initial_cond,t)
return ds.T
#=======================================================
This can then be called as
T = np.linspace(0, 5.0, 1000+1)
S,Z,R=eq(rates,y0,T)
but also in a way to produce only the values needed in the score function
Tm=np.append([0],Td)
Sm,Zm,Rm=eq(rates,y0,Tm)
This then simplifies the score function to
def score(parms):
#a.Get Solution to system
Sm,Zm,Rm=eq(parms,y0,Tm)
#c.Score Difference between model and data points
ss=lambda data,model:((data-model)**2).sum()
return ss(Zd,Zm[1:])
Now if you want for example strongly reject negative parameters, then you could change the return value to
return ss(Zd,Zm[1:]) + 1e6*sum(max(0,-x)**2 for x in parms)
which indeed renders all parameters positive (where previously in my notebook there was a negative first parameter).