Я новичок в Python и программировании. Я сделал приведенный ниже код, чтобы получить оптимальные параметры модели (R0, t_inc, t_rec, ex, teta) путем минимизации ошибки между данными и моделью (несколько дифференциальных уравнений). Я застрял в том, как определить функцию ошибки, как показано в приведенном ниже коде.
import numpy as np
import pandas as pd
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import datetime
from lmfit import Parameters, fit_report, minimize
totaldays = 93 # as of today
# This code is only to adjust the R0, t_infective, t_incubation to match data until to date
n_to_start = 0 # start data to fit
n_to_fit = totaldays # end data to fit -1,
NumberofTest = 136.73*n_to_fit**2 - 17609*n_to_fit + 605936
# getting data till to date
DataMalaysia = pd.read_csv('DataMalaysia.csv')
Dates = DataMalaysia.iloc[:,0].values
TotalCase = DataMalaysia.iloc[:,3].values
TotalRecovered = DataMalaysia.iloc[:, 5].values
TotalDeath = DataMalaysia.iloc[:,4].values
Days = DataMalaysia.iloc[:, 1].values
ActiveCase = DataMalaysia.iloc[:, 7].values
Daystofit = Days[:n_to_fit]
Dates = Dates[:n_to_fit]
start_date = datetime.date(2020,1,22)
ActiveCasetofit = ActiveCase[:n_to_fit]
TotalRecoveredtofit = TotalRecovered[:n_to_fit]
TotalDeathtofit = TotalDeath[:n_to_fit]
# parameter values including death and immigration
N = NumberofTest # number of population
i_initial = 4 # 4 people is infected at the beginning 25th Jan 2020
# parameters around the Susceptible population (possible to get infected)
immigrating_s = 0 # fraction of population immigrating into the infected location
death_s = 0 # fraction of population died due to other diseases
# parameters around the Exposed or Infected people (but not yet Infecting)
immigrating_e = 0 # fraction of the infected people immigrating into the infected location
death_e = 0 # fraction of the infected/exposed people die due to other diseases
#parameters around the Infectious population
immigrating_i = 0 # fraction of the infectious people immigrating into the infected location
# death_i_MCO = 0.0157 # fraction of the infectious people die due to the virus
# mitigation effort
# variables to fit
R0 = 2.96 # reproduction number. This number is relatively high
t_inc = 11.93 # incubation period (5-6 is most reported one)
t_rec = 1.24 # infectious period, gamma = 1/t_infectious is the recovery rate, typical 3-4 days
ex = 0.016 #death fraction
teta = 0.1 # recovered fraction without getting ill
# using population
e0 = 0
i0 = i_initial
r0 = 0
d0 = 0
rprime0 = 0
s0 = N - e0 - i0 - r0 - d0 - rprime0
# SEIR model including MCO
def SEIR(x, t, R0, t_inc, t_rec, ex, teta):
# introduction of the variables to calculate
s, e, i, r, rprime, d = x
alpha = 1/t_inc
gamma = 1/t_rec
R0t = R0/N
beta = R0t*gamma
# the differential equations
dsdt = -(1-u)*beta * s * i + (immigrating_s - death_s)*s
dedt = (1-u)*beta * s * i - alpha*e - teta*e + (immigrating_e - death_e)*e
didt = alpha * e - gamma * i + (immigrating_i - ex)*i
drdt = gamma*i
drprimedt = teta*e
dddt = ex*i
return [dsdt, dedt, didt, drdt, drprimedt, dddt]
# integrating the SEIR model
def integrate_i(t, R0, t_inc, t_rec, ex, teta):
x0 = s0, e0, i0, r0, rprime0, d0
solution = odeint(SEIR, x0, t, args = (R0, t_inc, t_rec, ex, teta)).T
solutiona = solution.T
return solutiona[:, 2]
def integrate_r(t, R0, t_inc, t_rec, ex, teta):
x0 = s0, e0, i0, r0, rprime0, d0
solution = odeint(SEIR, x0, t, args = (R0, t_inc, t_rec, ex, teta)).T
solutiona = solution.T
return solutiona[:, 3]
def integrate_d(t, R0, t_inc, t_rec, ex, teta):
x0 = s0, e0, i0, r0, rprime0, d0
solution = odeint(SEIR, x0, t, args = (R0, t_inc, t_rec, ex, teta)).T
solutiona = solution.T
return solutiona[:, 5]
def integrate_total(t_total, R0, t_inc, t_rec, ex, teta):
#slicing the time frame to each integration
ti = t_total[:n_to_fit]
td = t_total[len(ti)+len(ti):]
tr = t_total[len(ti):len(t_total)-len(td)]
result_i = integrate_i(ti, R0, t_inc, t_rec, ex, teta)
result_r = integrate_r(tr, R0, t_inc, t_rec, ex, teta)
result_d = integrate_d(td, R0, t_inc, t_rec, ex, teta)
return np.concatenate([result_i, result_r, result_d])
def error(t_total, R0, t_inc, t_rec, ex, teta):
R0 = 2.96 # reproduction number. This number is relatively high
t_inc = 11.93 # incubation period (5-6 is most reported one)
t_rec = 1.24 # infectious period, gamma = 1/t_infectious is the recovery rate, typical 3-4 days
ex = 0.016 #death fraction
teta = 0.1
total_error = (np.sum((integrate_total(t_total, R0, t_inc, t_rec, ex, teta)-y_total)**2))
return total_error
# fitting predictions with data points
start = n_to_start
end = n_to_fit
t = np.linspace(start, end, end)
ta = np.array(t)
# t_total = np.append(ta, ta, ta)
t_total = np.concatenate([ta, ta, ta])
y_total = np.concatenate([ActiveCasetofit, TotalRecoveredtofit, TotalDeathtofit])
p0=[R0, t_inc, t_rec, ex, teta]
params, extras = minimize(error, p0 ,method='BFGS',
options={'disp':True})
# Getting the optimized variables to plot what happens after to date
# getting the optimum values of R0, t_inc, t_rec, ex, teta
R0 = params[0]
t_incubation = params[1]
t_recovery = params[2]
death_ratio = params[3]
recovery_ratio = params[4]
#generation of the fitting curve
Predicted_ActiveCase = integrate_i(t, *params)
Predicted_RecoveredCase = integrate_r(t, *params)
Predicted_Death = integrate_d(t, *params)
print("Optimum R0 is {0:.2f}".format(params[0]))
print("Optimum Incubation Period is {0:.2f} days".format(params[1]))
print("Optimum Recovery Period is {0:.2f} days".format(params[2]))
print("Optimum Death Ratio is {0:.4f} ".format(params[3]))
print("Optimum Recovery Ratio Without Getting Ill is {0:.4f} ".format(params[4]))
# plotting data and results
fig1 = plt.figure()
t_fit = np.array([start_date+datetime.timedelta(days=i) for i in range(n_to_fit)])
plt.scatter(t_fit, ActiveCasetofit, c='red', label ='Data To Date')
plt.plot(t_fit, Predicted_ActiveCase, "r", label ='Fitted Active Case To Date')
plt.scatter(t_fit, TotalRecoveredtofit, c='blue', label ='Data To Date')
plt.plot(t_fit, Predicted_RecoveredCase, "b", label ='Fitted Total Recovered Case To Date')
plt.scatter(t_fit, TotalDeathtofit, c='green', label ='Data To Date')
plt.plot(t_fit, Predicted_Death, "g", label ='Fitted Death To Date')
plt.title('Fitting Data to Date')
plt.xlabel('Time/days')
plt.ylabel('Population')
plt.legend(loc='best')
Это дает мне следующую ошибку, я думаю, это связано с тем, что я не знаю, как ввести аргументы в код. В этом вся ошибка:
runfile ('C: /.../ Python / SEIR_v8.py', wdir = 'C: /.../ Python') Трассировка (последний вызов последним):
Файл "C: ... \ Python \ SEIR_v8.py", строка 171, в params, extras = minim (error, p0, method = 'BFGS')
Файл "C: \ Users ... \ AppData \ Local \ Continuum \ anaconda3 \ envs \ PythonNew \ lib \ site-packages \ lmfit \ minimizer.py", строка 2393, в сворачивании return fitter.minimize (method = method)
Файл "C: \ Users ... \ AppData \ Local \ Continuum \ anaconda3 \ envs \ PythonNew \ lib \ site-packages \ lmfit \ minimizer.py", строка 2176, в функции минимального возврата (** kwargs)
Файл "C: \ Users ... \ AppData \ Local \ Continuum \ anaconda3 \ envs \ PythonNew \ lib \ site-packages \ lmfit \ minimizer.py", строка 931, в scalar_minimize ret = scipy_minimize (self.penalty, переменные, ** fmin_kws)
Файл "C: \ Users ... \ AppData \ Local \ Continuum \ anaconda3 \ envs \ PythonNew \ lib \ site-packages \ scipy \ optimize_minimize.py", строка 604, в минимальном размере return _minimize_bfgs (fun, x0, args, jac , обратный вызов, ** варианты)
Файл "C: \ Users ... \ AppData \ Local \ Continuum \ anaconda3 \ envs \ PythonNew \ lib \ site-packages \ scipy \ optimize \ optimize.py", строка 1003, в _minimize_bfgs old_fval = f (x0)
Файл "C: \ Users ... \ AppData \ Local \ Continuum \ anaconda3 \ envs \ PythonNew \ lib \ site-packages \ scipy \ optimize \ optimize.py", строка 327, в функции возврата function_wrapper (* (wrapper_args + args ))
Файл "C: \ Users ... \ AppData \ Local \ Continuum \ anaconda3 \ envs \ PythonNew \ lib \ site-packages \ lmfit \ minimizer.py", строка 598, в штрафной r = self .__ остаток (fvars, apply_bounds_transformation)
Файл "C: \ Users ... \ AppData \ Local \ Continuum \ anaconda3 \ envs \ PythonNew \ lib \ site-packages \ lmfit \ minimizer.py", строка 530, в __residual out = self.userfcn (params, * self .userargs, ** self.userkws)
TypeError: error () отсутствуют 5 обязательных позиционных аргументов: 'R0', 't_inc', 't_rec', 'ex' и 'teta'
Пожалуйста помоги.
С уважением,
Зульфан
R0
и т. Д. В качестве параметров, а затем устанавливаете для них фиксированные значения? - person Code-Apprentice   schedule 25.04.2020minimize
изlmfit
иscipy.optimize
. Они работают по-разному и не взаимозаменяемы. - person M Newville   schedule 26.04.2020