当我尝试在Python中执行参数估计时为什么参数不更改？

Question

我正在本文的帮助下尝试：http://adventuresinpython.blogspot.com/2012/08/fitting-differential-equation-system-to.html，但是参数保持不变，与选择何种初始条件无关。

#Zombie Display
# zombie apocalypse modeling
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
#The Model
#=======================================================
def eq(par,initial_cond,start_t,end_t,incr):
     #-time-grid-----------------------------------
     t  = np.linspace(start_t, end_t,incr)
     #differential-eq-system----------------------
     def funct(phi,t):
         S=phi[0]
         E=phi[1]
         I=phi[2]
         R=phi[3]
         dS_dt = mu*(N-S)-beta*S*I/N-nu*S
         dE_dt = beta*S*I/N-(mu+sigma)*E
         dI_dt = sigma*E-(mu+gamma)*I
         dR_dt=gamma*I-mu*R+nu*S
         dphi_dt = [dS_dt,dE_dt,dI_dt,dR_dt]
         return dphi_dt
     #integrate------------------------------------
     ds = integrate.odeint(funct,initial_cond,t)
     return (ds[:,0],ds[:,1],ds[:,2],ds[:,3],t)
#=======================================================

#=====================================================

#1.Get Data
#====================================================
Td=np.array([1,2,3,4,5,6,7,8,9,10,11,12,13,14,15])#time
Zd=np.array([274,326,547,639,2000,2700,4400,6000,7700,9700,11200,14300,17200,19200,20000])#zombie pop
#====================================================

#2.Set up Info for Model System
#===================================================
# model parameters
#----------------------------------------------------
beta=1.4
gamma=0.06
sigma=0.15
mu=0
nu=0
rates=(beta,gamma,sigma,mu,nu)

# model initial conditions
#---------------------------------------------------
N=11*pow(10,6)
S0 = N-1               # initial population
E0 = 105.1                  # initial zombie population
I0 = 27.679                  # initial death population
R0= 2.0
y0 = [S0, E0, I0, R0]      # initial condition vector

# model steps
#---------------------------------------------------
start_time=0.0
end_time=140
intervals=5.0
mt=np.linspace(start_time,end_time,intervals)

# model index to compare to data
#----------------------------------------------------
findindex=lambda x:np.where(mt>=x)[0][0]
mindex=list(map(findindex,Td))
#=======================================================



#3.Score Fit of System
#=========================================================
def score(parms):
    #a.Get Solution to system
    F0,F1,F2,F3,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]
beta, gamma, sigma, mu, nu=answ[0]
newrates=(beta,gamma,sigma,mu,nu)
#=======================================================

#5.Generate Solution to System
#=======================================================
F0,F1,F2,F3,T=eq(newrates,y0,start_time,end_time,intervals)
Zm=F1[mindex]
Tm=T[mindex]
#======================================================




#6. Plot Solution to System
#=========================================================
plt.figure()
plt.plot(T,F1,'b-',Tm,Zm,'ro',Td,Zd,'go')
plt.xlabel('days')
plt.ylabel('population')
title='Zombie Apocalypse  Fit Score: '+str(bestscore)
plt.title(title)
plt.show()

我知道这是巨大的，但是如果有人是微分方程系统的参数估计专家，那么我将很高兴听到任何信息。非常感谢！

Answer 1

您不会看到参数的任何变化，因为未使用它们，分数函数是恒定的。当您将parms传递给eq时，在eq内部，您无需解压缩par参数的包装。

在eq中添加为第一行

beta,gamma,sigma,mu,nu = par

并且最小化算法做了不平凡的事情。仍然有可能在合理的时间内找不到解决方案。将maxiter设置为更合理的值。如果该方法失败了，那也可能是问题制定的问题。也许初始条件也需要是优化变量，否则ODE系统并不十分适合。