Почему параметры не меняются, когда я пытаюсь ускорить оценку параметра в Python? - PullRequest
Новая
       7

Почему параметры не меняются, когда я пытаюсь ускорить оценку параметра в Python?

0 голосов
/ 06 февраля 2020

Я пытаюсь с помощью этого поста: 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()

Я знаю, что это огромен, но если кто-то является экспертом в оценке параметров системы дифференциальных уравнений, я был бы рад услышать любую информацию. Большое спасибо!

1 Ответ

0 голосов
/ 08 февраля 2020

Вы не видите никаких изменений в параметрах, потому что вы не используете их, функция счета постоянна. Когда вы передаете parms в eq, внутри eq вы не распаковываете параметр par.

Добавьте eq в качестве первой строки 

beta,gamma,sigma,mu,nu = par

и алгоритм минимизации делает что-то нетривиальное. Все еще возможно, что никакое решение не будет найдено в разумное время. Установите maxiter на что-то более разумное. Если метод потерпит неудачу, это также может быть связано с формулировкой проблемы. Возможно, начальные условия также должны быть переменными оптимизации, или система ODE не очень подходит.

...