Я новичок в Python и программировании. Я сделал приведенный ниже код, чтобы получить оптимальные параметры модели (R0, t_in c, t_re c, 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') Traceback ( последний вызов был последним):
Файл "C: ... \ Python \ SEIR_v8.py", строка 171, в параметрах, дополнительные = минимизировать (ошибка, p0, метод = 'BFGS')
Файл "C: \ Users ... \ AppData \ Local \ Continuum \ anaconda3 \ envs \ PythonNew \ lib \ site-packages \ lmfit \ minimizer.py", строка 2393, в установщике минимизации возврата. свернуть (метод = метод)
Файл "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, variable, ** fmin_kws)
Файл" C: \ Users ... \ AppData \ Local \ Continuum \ anaconda3 \ envs \ PythonNew \ lib \ site-packages \ scipy \ optimize_minimize.py ", строка 604, в минимизации возврата _minimize_bfgs (весело, x0, args, ja c, обратный вызов, ** опции)
Файл "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)
Ошибка типа: error () отсутствует 5 обязательных позиционных аргументов: 'R0', 't_in c', 't_re c', 'ex' и 'teta'
Пожалуйста, помогите.
С уважением,
Zulfan