bolivia go brrr

bolivia-go-brrr.py

from scipy.integrate import solve_ivp
import pandas as pd
from datetime import datetime, timedelta
from numpy import linspace

# Estima el R0 que mejor se ajusta al número de casos confirmados reportados ayer y produce el número de casos para una fecha en el futuro

def modelo_seir(current_timepoint, state_values, parameters):
  s,e,i,r = state_values
  dS = -parameters['beta'] * s * i
  dE = (parameters['beta'] * s * i) - (parameters['delta'] * e)
  dI = (parameters['delta'] * e) - (parameters['gamma'] * i)
  dR = parameters['gamma'] * i
  return [dS, dE, dI, dR]

def ro_ayer(casos_confirmados, start_date, target_date):
  days = (datetime.today() - start_date).days
  target = (target_date - start_date).days
  x = 27 # infectados
  y = 0 # removidos
  z = x * 10 # expuestos
  w = 11595000 - x - z
  n = w + x + y + z
  valores_iniciales = [w/n, x/n, y/n, z/n]
  duracion = list(range(1,target + 2))
  probabilidad_transmision = 0.04
  periodo_infeccion = 7.5
  periodo_latencia = 6
  valor_gamma = 1 / periodo_infeccion
  valor_delta = 1 / periodo_latencia
  hoy = pd.DataFrame(columns = ['contacto', 'ro', 'ayer', 'final'])
  for ratio_contacto in linspace(3,13,201):
    valor_beta = ratio_contacto * probabilidad_transmision
    parameters = {'beta': valor_beta, 'gamma': valor_gamma, 'delta': valor_delta}
    ro = valor_beta / valor_gamma
    resultado = solve_ivp(modelo_seir, t_span=(0, target + 2), t_eval=duracion, y0=valores_iniciales, args=[parameters], method='LSODA', rtol=1e-6)
    df = pd.DataFrame(resultado.y).T
    df.columns = ['s','e','i','r']
    casos = pd.DataFrame((df['i'] + df['r']).apply(lambda x: int(x * 11595000)), columns=['casos'], index=duracion)
    hoy_casos = [ratio_contacto, ro]
    hoy_casos.extend(casos.iloc[days - 1].to_list())
    hoy_casos.extend(casos.iloc[target - 1].to_list())
    hoy = hoy.append(pd.DataFrame([hoy_casos], columns = ['contacto', 'ro', 'ayer', 'final']))
  print(hoy[(abs(hoy['ayer']-casos_confirmados) < 25)].to_markdown(tablefmt='orgtbl', showindex=False))

ro_ayer(1110, datetime(year=2020,month=3,day=22), datetime(year=2020,month=5,day=30))

Sample result with data from April 30

| contacto |    ro | ayer | final |
|----------+-------+------+-------|
|      8.5 |  2.55 | 1092 | 12514 |
|     8.55 | 2.565 | 1110 | 13032 |
|      8.6 |  2.58 | 1128 | 13570 |

Sample result with data from May 3

|   contacto |    ro |   ayer |   final |
|------------+-------+--------+---------|
|       8.7  | 2.61  |   1579 |   14709 |
|       8.75 | 2.625 |   1610 |   15314 |

Sample result with data from May 10

|   contacto |    ro |   ayer |   final |
|------------+-------+--------+---------|
|       8.6  | 2.58  |   2430 |   13570 |
|       8.65 | 2.595 |   2491 |   14128 |

Sample result with data from May 21

|   contacto |    ro |   ayer |   final |
|------------+-------+--------+---------|
|       8.25 | 2.475 |   4865 |   10209 |
|       8.3  | 2.49  |   5032 |   10632 |