как согласовать данные с уравнениями с помощью минимизации в Python для получения параметров модели

Я новичок в 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'

Пожалуйста помоги.

С уважением,

Зульфан


person Zulfan    schedule 24.04.2020    source источник
comment
Какая строка вызывает ошибку?   -  person Code-Apprentice    schedule 24.04.2020
comment
Это один? Файл C: \ Users \ ... \ AppData \ Local \ Continuum \ anaconda3 \ envs \ PythonNew \ lib \ site-packages \ scipy \ optimize \ optimize.py, строка 327, в функции возврата function_wrapper (* (wrapper_args + args) ) TypeError: error () отсутствуют 6 обязательных позиционных аргументов: 't_inc', 't_rec', 'ex', 'teta', 't_total' и 'y_total'   -  person Zulfan    schedule 24.04.2020
comment
Пожалуйста, отредактируйте свой вопрос, чтобы показать полную трассировку стека.   -  person Code-Apprentice    schedule 24.04.2020
comment
вопрос обновлен   -  person Zulfan    schedule 25.04.2020
comment
Почему вы берете R0 и т. Д. В качестве параметров, а затем устанавливаете для них фиксированные значения?   -  person Code-Apprentice    schedule 25.04.2020
comment
ваш пример смешивает сигнатуру вызова с функцией minimize из lmfit и scipy.optimize. Они работают по-разному и не взаимозаменяемы.   -  person M Newville    schedule 26.04.2020


Ответы (2)


Вы звоните minimize() неправильно. Я слишком устал, чтобы разбираться в деталях. Предлагаю вам внимательно прочитать документацию scipi. Вам нужно передать вектор x0 в качестве начального предположения, а также args, которые являются фиксированными параметрами функции, которую вы минимизируете. В существующем виде вы передаете только x0, но ваша функция error ожидает дополнительных параметров.

person Code-Apprentice    schedule 25.04.2020

Я изменил функцию ошибки, как показано ниже, и она работает.

def error(params, t_total, y_total):
    R0 = params['R0'].value
    t_inc = params['t_inc'].value
    t_rec = params['t_rec'].value
    ex = params['ex'].value
    teta = params['teta'].value

    y_model = integrate_total(t_total, R0, t_inc, t_rec, ex, teta)
    return y_model - y_total


params = lmfit.Parameters()
params.add('R0', 3.65, min=0, max=5.0)
params.add('t_inc', 15, min=0, max=20.0)
params.add('t_rec', 1, min=0, max=2.0)
params.add('ex', 0.02, min=0, max=1.0)
params.add('teta', 0.1, min=0, max=1.0)

# 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.concatenate([ta, ta, ta])
y_total = np.concatenate([ActiveCasetofit, TotalRecoveredtofit, TotalDeathtofit])

o1 = lmfit.minimize(error, params, args=(t_total, y_total), method='powell')
print("# Fit using leastsq:")
lmfit.report_fit(o1)

Затем я могу изменить решающую программу и начальные точки, чтобы добиться лучшего соответствия.

person Zulfan    schedule 25.04.2020