Hello everyone! Is anyone here familiar with pymc/sunode and sympy? Now I have a problem.
My code is an MCMC simulation of the number of new infections with the flu, and I use pymc and sunode. The data I have is the indoor temperature and humidity per day, and the number of new cases of flu per day. I wanted to use the formula of temperature and humidity to express its effect on the influenza transmission rate, bring the transmission rate into the equations of the SIRS model, and finally use the mcmc to fit and observe the parameters.
But I’m having some problems right now. In my code, the input t (time) of the sunode model is a symbol type, but I want to get the t at a certain step at the time of sampling, so that I can get the values of RH and Temp at that time t, and correspond. RH and Temp are a list, but t is a symbol.The core of the problem is that the symbolic variable t conflicts with the type of the numeric index, and I don’t know how to pass this expression into the β (infection rate) of the model.
All suggestions would be appreciated!
Here is my code:
import matplotlib.pyplot as plt
import numpy as np
import pymc as pm
import arviz as az
import sunode
import sunode.wrappers.as_pytensor
import sympy
import csv
class DiseaseData:
def __init__(self, file_path):
self.data = self.csv_to_city_dict(file_path)
def csv_to_city_dict(self, file_path):
#get data from csv file
city_data = {}
with open(file_path, mode='r', encoding='utf-8') as csvfile:
reader = csv.DictReader(csvfile)
for row in reader:
city = row['CITY']
if city not in city_data:
city_data[city] = {
'time': [],
'TEMP-Indoor': [],
'RH-Indoor': [],
'TEMP-Outdoor': [],
'RH-Outdoor': [],
'FLU-Incidence': [],
'FLU-ILI': [],
}
try:
city_data[city]['TEMP-Indoor'].append(float(row['TEMP-Indoor']))
city_data[city]['RH-Indoor'].append(float(row['RH-Indoor']))
city_data[city]['TEMP-Outdoor'].append(float(row['TEMP-Outdoor']))
city_data[city]['RH-Outdoor'].append(float(row['RH-Outdoor']))
city_data[city]['time'].append(int(row['TIME']))
except ValueError:
print(f"Skipping invalid value in row {row}")
city_data[city]['FLU-Incidence'].append(row['FLU-Incidence'])
city_data[city]['FLU-ILI'].append(row['FLU-ILI'])
'''city_time_data = {}
for i in range(len(city_data[city]['time'])):
city_time_data[city_data[city]['time'][i]] = {
'TEMP-Indoor': city_data[city]['TEMP-Indoor'][i],
'RH-Indoor': city_data[city]['RH-Indoor'][i],
'TEMP-Outdoor': city_data[city]['TEMP-Outdoor'][i],
'RH-Outdoor': city_data[city]['RH-Outdoor'][i],
'FLU-Incidence': city_data[city]['FLU-Incidence'][i],
'FLU-ILI': city_data[city]['FLU-ILI'][i],
}'''
return city_data
class DiseaseModel:
def __init__(self, data, city):
#initialize the data
self.data = data
self.city = city
self.TEMP = data[city]['TEMP-Indoor']
self.RH = data[city]['RH-Indoor']
self.FLU = data[city]['FLU-Incidence']
self.ILI = data[city]['FLU-ILI']
self.times = self.times = np.arange(0, len(self.ILI))
def get_RH(self,t):
return self.RH[t]
def get_TEMP(self,t):
return self.TEMP[t]
def SIRSmodel(self, t, y, p):
# define β = [aRH(t)2 + bRH(t) + c][Tc/T(t)]**Texp
# beta is the infection rate
# here is the crux, t presents a symbol type, but the function get_RH/get_TEMP needs a number type
beta = (p.a * self.get_RH(t)**2 + p.b * self.get_RH(t) + p.c) * (p.Tc / self.get_TEMP(t))**p.Texp
return {
"S": (p.n - y.S - y.I) / p.L - (beta * y.S * y.I) / p.n,
"I": (beta * y.S * y.I) / p.n - y.I / p.D,
"New_I": beta * y.S * y.I / p.n,
}
def fit(self):
with pm.Model() as self.model:
# Priors for unknown model parameters
a = pm.Normal('a', mu=0, sigma=1)
b = pm.Normal('b', mu=0, sigma=1)
c = pm.Normal('c', mu=0, sigma=1)
Tc = pm.Uniform('Tc', lower=20, upper=24)
Texp = pm.Uniform('Texp', lower=0.5, upper=2)
L = pm.Gamma('L', alpha=1, beta=1)
D = pm.Gamma('D', alpha=1, beta=1)
n = 1e6
res, _, problem, solver, _, _ = sunode.wrappers.as_pytensor.solve_ivp(
y0={
"S": (n - 1,()),
"I": (1,()),
"New_I": (0,()),
},
params={
"a": (a, ()),
"b": (b, ()),
"c": (c, ()),
"Tc": (Tc, ()),
"Texp": (Texp, ()),
"L": (L, ()),
"D": (D, ()),
"n": (n, ()),
'_dummy': (np.array(1.), ()),
},
tvals=self.times[1:],
rhs=self.SIRSmodel,
t0=self.times[0],
)
New_I = res['New_I']
New_I_obs = pm.Poisson('New_I_obs', mu=New_I, observed=self.ILI[1:])
step = pm.Metropolis()
trace = pm.sample(1000, step=step)
return trace
def analysis(self, trace):
# analysis
axes = az.plot_trace(trace)
fig = axes.ravel()[0].figure
fig.savefig('REplot.png')
res = str(az.summary(trace))
az.plot_forest(trace, r_hat=True)
az.plot_posterior(trace)
with open('rate.txt', 'a', encoding='utf-8') as file:
file.write('\r\n' + self.country + res)
with self.model:
ppc_samples = pm.sample_posterior_predictive(trace)
New_I_mean = ppc_samples.posterior_predictive['New_I_obs'].mean(axis=0)
plt.plot(self.times[1:] ,self.ILI, "o", color="r", lw=0.1, ms=1.5, label="Observed")
for i in range(New_I_mean.shape[0]):
plt.plot(self.times[1:], New_I_mean[i, :], color="b", lw=0.3)
plt.legend()
plt.savefig('ILI.png')
def evaluate(self):
pass
def plot(self):
pass
if __name__ == '__main__':
data = DiseaseData('Tianjing V2.csv').data
city = 'Tianjin'
model = DiseaseModel(data, city)
trace = model.fit()
model.analysis(trace)
model.evaluate()
model.plot()
File is like
CITY,TIME,TEMP-Indoor,RH-Indoor,TEMP-Outdoor,RH-Outdoor,FLU-Incidence,FLU-ILI
Tianjin,1,22.5,45.0,22.5,45.0,0.0,0.0
Tianjin,2,22.3,44.8,22.3,44.8,1.2,0.5
...