<jupyter_start><jupyter_text><jupyter_code>!pip install pymc3==3.8

%matplotlib inline
import numpy as np
from IPython.display import display, Markdown
import matplotlib.pyplot as plt
import matplotlib
import pandas as pd
import seaborn as sns

import pymc3 as pm
from pymc3.ode import DifferentialEquation
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
import arviz as az
import theano

plt.style.use('seaborn-darkgrid')<jupyter_output>Requirement already satisfied: pymc3==3.8 in /usr/local/lib/python3.6/dist-packages (3.8)
Requirement already satisfied: theano>=1.0.4 in /usr/local/lib/python3.6/dist-packages (from pymc3==3.8) (1.0.4)
Requirement already satisfied: tqdm>=4.8.4 in /usr/local/lib/python3.6/dist-packages (from pymc3==3.8) (4.38.0)
Requirement already satisfied: h5py>=2.7.0 in /usr/local/lib/python3.6/dist-packages (from pymc3==3.8) (2.10.0)
Requirement already satisfied: patsy>=0.4.0 in /usr/local/lib/python3.6/dist-packages (from pymc3==3.8) (0.5.1)
Requirement already satisfied: numpy>=1.13.0 in /usr/local/lib/python3.6/dist-packages (from pymc3==3.8) (1.18.2)
Requirement already satisfied: arviz>=0.4.1 in /usr/local/lib/python3.6/dist-packages (from pymc3==3.8) (0.7.0)
Requirement already satisfied: scipy>=0.18.1 in /usr/local/lib/python3.6/dist-packages (from pymc3==3.8) (1.4.1)
Requirement already satisfied: pandas>=0.18.0 in /usr/local/lib/python3.6/dist-packages (from pymc3==3.8) (0.25.3)
Require[...]<jupyter_text>Definimos la ecuación diferencial SEIR - utilzamos un modelo simplificado para fittear y obtener la Ro<jupyter_code>def SEIR(y, t, p):
    """
    y[0] = S
    y[1] = E
    y[2] = I
    y[3] = R

    
    p[0] = Alpha
    p[1] = Beta
    p[2] = Gamma


    St = S[-1] - (beta*S[-1]*I[-1])*dt
    Et = E[-1] + (beta*S[-1]*I[-1] - alpha*E[-1])*dt
    It = I[-1] + rho*(alpha*E[-1] - gamma*I[-1])*dt
    Rt = gamma*I[-1]
    """



    ds = -p[1]*y[0]*y[2]
    de = p[1]*y[0]*y[2] - p[0]*y[1]
    di = p[0]*y[1] - p[2]*y[2]
    dr = p[2]*y[2]
    return [ds, de, di, dr]<jupyter_output><empty_output><jupyter_text>FUNCIONES DE UTILIDAD<jupyter_code>def TimeSeriesLoader(url):

  """
  Este Loader se puede cambiar según la fuente, en ese caso, tambien habría que cambiar un poco la función
  GetObs_and_times(estado,Pop)
  """
  df = pd.read_csv(url)

  df['Reporte'] = pd.to_datetime(df['Reporte'])
  df = (df.drop(["Procedencia", "Llegada","Inicio","Caso","Sexo","Edad","Confirmación"], axis=1)
          .dropna()
          .rename(columns = {"Reporte":"Dia"})
          .sort_values(by = ["Dia"])
          .reset_index(drop = True)
        )
  df["Casos_Acumulados"] = np.nan
  Estados = pd.DataFrame({
    "Estado":[],
    "Dia":[],
    "Casos_Acumulados":[]
  })

  for estado in df.Estado.unique():
      df.loc[(df.Estado == estado),"Casos_Acumulados"] =  np.arange(1,len(df.loc[(df.Estado == estado),"Casos_Acumulados"])+1,1)
      tmp = df.loc[df.Estado == estado].drop_duplicates("Dia", keep = 'last').reset_index(drop = True)
      Estados = pd.concat([Estados,tmp])
  return Estados
def GetObs_and_times(estado,Pop,r,d):
  from scipy.interpolate import InterpolatedUnivariateSpline
  """
  Esta funcion genera las condiciones inciales para PyMC3 a partir de las observaciones de 
  Infectados en un estado y dada la población del estado (se debe ingresar el número)
  como los registros de recuperados en mexico son nulos se asumen 0.
  """
  yobs = Estados.loc[Estados.Estado == estado,"Casos_Acumulados"].values
  days = Estados.loc[Estados.Estado == estado,"Dia"].dt.day.values
  days = (days - days[0])
  times = np.arange(0,days[-1]+1,1)
  
  interp = InterpolatedUnivariateSpline(days,yobs,k=2)
  obs_interp = interp(times)


  def SEIR_Day(obs,r=0, d=0):
    I = obs / Pop
    E = (obs*4) / Pop # Ojo aca con el 4
    D = d
    R = r / Pop
    S = 1 - (E-I-R)
    return [S, E, I, R]

  yobs_arr = np.array([SEIR_Day(obs,r,d) for obs in obs_interp])

## Crear un arreglo de N x 4 para enviar como estados al sistema


  return obs_interp, yobs_arr, times<jupyter_output><empty_output><jupyter_text>Condiciones iniciales para la ODE, verificamos salida de ODEINT<jupyter_code>np.random.seed(666)

# JALISCO (2020-03-14)
Population = 8000000
I_o = 2 / Population  # Tenemos 32 casos
E_o = (2*4)/ Population # Asumimos 4 expuestos por caso
S_o = (1) - (E_o+I_o) # El resto somos suceptibles
R_o = 0 # NO hay ningun recuperado

dias_evaluacion = 120
dt=1
periodo_evaluacion = np.arange(0,dias_evaluacion+dt,dt)

R_o = 3.5 #Escenario Base
alpha = 0.2
gamma = 0.5
beta = R_o * gamma

y0 = np.array([S_o,E_o,I_o,R_o])

args = ((alpha, beta, gamma,),)
 
y = odeint(SEIR,t=periodo_evaluacion,y0=y0, args=args, rtol=1e-08)
#@title Salida de ODEINT { display-mode: "form" }
import plotly.express as px
E_a = y[:,1]*Population
I_a = y[:,2]*Population
Clase = np.array(["Expuestos"]*y[:,1].shape[0] + ["Infectados"]*y[:,1].shape[0])
Dias = np.concatenate([periodo_evaluacion,periodo_evaluacion])
SEIR_df = pd.DataFrame({
    "Casos": np.concatenate([E_a,I_a]),
    "Clase": Clase,
    "Dias" : Dias
})

fig = px.line(SEIR_df, x="Dias", y="Casos", color='Clase',color_discrete_sequence=["green", "red", "goldenrod", "blue"], template = "ggplot2")

fig.update_layout(
    title=f"Predicción del modelo SEIR de la evolución de COVID-19 en Jalisco, Escenario: BASE Ro = 3.5",
    
    xaxis_title="Días",
    yaxis_title="Casos Totales",
    )

for trace in fig.data:
    trace.name = trace.name.split('=')[1]

fig.show()<jupyter_output><empty_output><jupyter_text>Por lo tanto está bien definido el sistema de ecuaciones.<jupyter_code>Estados = TimeSeriesLoader("https://gist.githubusercontent.com/said3427/18f39eab460766a17ea4802607dd5522/raw")
Estados.head()
y_interpol, obs,times = GetObs_and_times("Jalisco",8000000,0,0)
seir_model = DifferentialEquation(
    func=SEIR,
    times=times,
    n_states=4,
    n_theta=3,
)
sigma = 10 / Population # numero de error (prior) de casos por dia aquí se asumen 10
with pm.Model() as model:

    sigma = pm.HalfCauchy('sigma', sigma, shape=4) 
    #p_alpha = pm.Bound(pm.Normal, lower=0, upper=1)("alpha",.5, 0.3 )  
    p_gamma = pm.Bound(pm.Normal, lower=0, upper=1)("gamma",.5, 0.37 )  
    R0 = pm.Bound(pm.Normal, lower=1, upper=5)('R0', 3.5, 1)
    p_beta = pm.Deterministic('beta', gamma*R0)
    


    seir_curves = seir_model(y0=obs[0], theta=[alpha, p_beta, p_gamma])
    Y = pm.Normal('Y', mu=seir_curves, sd=sigma, observed=obs)
    
    prior = pm.sample_prior_predictive()
    trace = pm.sample(1000,tune=500, target_accept=0.95, cores=2)
    posterior_predictive = pm.sample_posterior_predictive(trace)

    data = az.from_pymc3(trace=trace, prior = prior, posterior_predictive = posterior_predictive)
az.plot_posterior(data,round_to=2, credible_interval=0.95)
ppc_samples = posterior_predictive["Y"]
mean_ppc = ppc_samples.mean(axis=0)
CriL_ppc = np.percentile(ppc_samples,q=2.5,axis=0)
CriU_ppc = np.percentile(ppc_samples,q=97.5,axis=0)
plt.plot(times,y_interpol,'o', color='b', lw=1, label='Infected cases observed')
plt.plot(times,mean_ppc[:,3]*Population, color='g', lw=4, label=f'mean of $I(t)$ ppc')
plt.plot(times,CriL_ppc[:,3]*Population, '--',  color='g', lw=2, label=f'$I(t)$ credible intervals')
plt.plot(times,CriU_ppc[:,3]*Population, '--',  color='g', lw=2)
plt.legend()
plt.figure(figsize = (7,5))
plt.hist(prior["R0"],histtype="stepfilled",bins=30,alpha = .80, label = f"Prior of $R_0$",color = "#A60628", density = True)
plt.hist(trace["R0"],histtype="stepfilled",bins=30,alpha = .80, label = f"Predictive Prior of $R_0$",color = "#467821", density = True)
plt.legend()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.16456026058631923, 'fraction_numerical': 0.043387622149837136, 'mean_word_length': 4.394237526352776, 'pattern_counts': {'":': 6, '<': 20, '<?xml version=': 0, '>': 28, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '397983', 'license': 'MIT', 'n_tokens_mistral': 3309, 'n_tokens_neox': 3149, 'n_words': 637, 'path': 'notebooks/Bayesian_MCMC_fitter_State_Agnostic.ipynb', 'repo_name': 'ricalanis/COVID19_MEX_MASTER', 'size': '7675'}