Week 11: SIR Model and Covid Data

In [1]:
import numpy as np
import matplotlib.pyplot as plt
import Pythonic_files as pf
import seaborn as sns
%config InlineBackend.figure_format='retina'
In [4]:
plt.figure(figsize=(13,6))

beta , gamma = 0.5,1/16

def sfunc(s,i):
    return -beta*s*i
def ifunc(s,i):
    return beta*s*i-gamma*i

I0=0.001
S0=1-I0
dt=0.01

t,s,i = pf.Ode_2D_solve(sfunc,ifunc,S0,I0,dt,100)

plt.plot(t,s,label='Susceptible')
plt.plot(t,i,label='Infectious')
plt.plot(t,1-s-i,label='Recovered')

plt.xlabel('time in days/months/mucci')
plt.ylabel('Proportion of Population')

plt.legend()
plt.show()


Days=[]
for x in range(len(t)):
    if t[x]//1-t[x-1]//1>=1:
        Days.append(x)
Betas = []
for x in Days[1:]:
    b = -(s[x]-s[x-1])/dt/(s[x]*i[x])
    Betas.append(b)
    
sns.distplot(Betas)
plt.show()

How do we predict models, given our Betas

First lets create a "clipped" graph

In [75]:
plt.figure(figsize=(13,6))

beta , gamma = 0.5,1/16

def sfunc(s,i):
    return -beta*s*i
def ifunc(s,i):
    return beta*s*i-gamma*i

I0=0.001
S0=1-I0
dt=0.01

t,s,i = pf.Ode_2D_solve(sfunc,ifunc,S0,I0,dt,30)

plt.plot(t,s,label='Susceptible')
plt.plot(t,i,label='Infectious')
plt.plot(t,1-s-i,label='Recovered')

plt.xlabel('time in days/months/mucci')
plt.ylabel('Proportion of Population')

plt.xlim(0,100)
plt.legend()
plt.show()
In [76]:
Days=[]
for x in range(len(t)):
    if t[x]//1-t[x-1]//1>=1:
        Days.append(x)
In [77]:
Betas = []
for x in Days[1:]:
    b = -(s[x]-s[x-1])/dt/(s[x]*i[x])
    Betas.append(b)
In [78]:
sns.distplot(Betas)
plt.show()
In [79]:
estimate_beta=np.mean(Betas)
In [80]:
initial_s = s[-1]
initial_i = i[-1]
In [84]:
plt.figure(figsize=(13,6))

beta , gamma = 0.5,1/16

def sfunc(s,i):
    return -beta*s*i
def ifunc(s,i):
    return beta*s*i-gamma*i

I0=0.001
S0=1-I0
dt=0.01

t,s,i = pf.Ode_2D_solve(sfunc,ifunc,S0,I0,dt,35)

beta , gamma = estimate_beta,1/16
def sfunc(s,i):
    return -beta*s*i
def ifunc(s,i):
    return beta*s*i-gamma*i



t1,s1,i1 = pf.Ode_2D_solve(sfunc,ifunc,initial_s,initial_i,dt,5)


plt.plot(t,i,label='exact Infectious')

plt.plot(t1+30,i1,label = 'approximate Infectious')

plt.xlabel('time in days/months/mucci')
plt.ylabel('Proportion of Population')

plt.xlim(29,35)
plt.legend()
plt.show()

Obviously this is a great approximation, but its dis-honest, lets get some noisy data.

Lets introduce noise

In [91]:
plt.figure(figsize=(13,6))

beta , gamma = 0.5,1/16

def sfunc(s,i):
    return -beta*s*i
def ifunc(s,i):
    return beta*s*i-gamma*i

I0=0.001
S0=1-I0
dt=0.01
noise=1


t,s,i = pf.Ode_2D_solve_noise(sfunc,ifunc,S0,I0,dt,22, noise)

plt.plot(t,s,label='Susceptible')
plt.plot(t,i,label='Infectious')
plt.plot(t,1-s-i,label='Recovered')

plt.xlabel('time in days/months/mucci')
plt.ylabel('Proportion of Population')

plt.xlim(0,20)
plt.legend()
plt.show()


print(s[-1])

0.07959173903235509

Lets only look at the first 16 days

In [92]:
t
Out[92]:
array([0.000e+00, 1.000e-02, 2.000e-02, ..., 2.197e+01, 2.198e+01,
       2.199e+01])
In [93]:
Days=[]
for x in range(len(t)):
    if t[x]//1-t[x-1]//1>=1:
        if t[x-1]<=16:
            Days.append(x)
Betas = []
for x in Days[1:]:
    b = -(s[x]-s[x-1])/dt/(s[x]*i[x])
    Betas.append(b)
In [94]:
sns.distplot(Betas)
print(np.mean(Betas))

0.48836515898184174
In [95]:
t[1600]
Out[95]:
16.0
In [97]:
initial_s = s[1601]
initial_i = i[1601]
estimate_beta = np.mean(Betas)
In [100]:
plt.figure(figsize=(13,6))

beta , gamma = estimate_beta,1/16

def sfunc(s,i):
    return -beta*s*i
def ifunc(s,i):
    return beta*s*i-gamma*i

I0=0.001
S0=1-I0
dt=0.01


t1,s1,i1 = pf.Ode_2D_solve(sfunc,ifunc,initial_s,initial_i,dt,5)


plt.plot(t[:1600],i[:1600],label='exact Infectious')
plt.plot(t[1600:],i[1600:],c='r',label='exact Infectious')
plt.plot(t1+16,i1,label = 'approximate Infectious')

plt.xlabel('time in days/months/mucci')
plt.ylabel('Proportion of Population')

#plt.xlim(29,35)
plt.legend()
plt.show()
In [102]:
plt.figure(figsize=(13,6))

beta , gamma = estimate_beta,1/16

def sfunc(s,i):
    return -beta*s*i
def ifunc(s,i):
    return beta*s*i-gamma*i

I0=0.001
S0=1-I0
dt=0.01


t1,s1,i1 = pf.Ode_2D_solve(sfunc,ifunc,initial_s,initial_i,dt,5)


plt.plot(t[:1600],i[:1600],label='exact Infectious')
plt.plot(t[1600:],i[1600:],c='r',label='exact Infectious')
plt.plot(t1+16,i1,label = 'approximate Infectious')

plt.xlabel('time in days/months/mucci')
plt.ylabel('Proportion of Population')

#plt.xlim(29,35)
plt.legend()
plt.show()
Out[102]:
15
In [109]:
def Ode_2D_solve_Sample(sfunc,ifunc,x0,y0,dt,T,Betas):
    
    def sfunc(s,i):
        return -beta*s*i
    def ifunc(s,i):
        return beta*s*i-gamma*i


    
    
    Num = int(T//dt)
    
    
    x=np.ones(Num+1)*x0
    y=np.ones(Num+1)*y0
    
    t=np.linspace(0,Num*dt,len(x))
    
    for count in range(Num):
        beta = np.random.choice(Betas)
        x[count+1] = x[count]+sfunc(x[count],y[count])*dt
        y[count+1] = y[count]+ifunc(x[count],y[count])*dt
        
    return t,x,y
In [127]:
plt.figure(figsize=(13,6))

t1,s1,i1 = Ode_2D_solve_Sample(sfunc,ifunc,initial_s,initial_i,dt,5,Betas)

plt.plot(t[:1600],i[:1600],c='g',linewidth=0.6,label='exact Infectious')
plt.plot(t[1600:],i[1600:],c='r',linewidth=0.6,label='exact Infectious')
plt.plot(t1+16,i1,c='b',alpha = 0.3,linewidth=0.6,label = 'approximate Infectious')

plt.xlabel('time in days/months/mucci')
plt.ylabel('Proportion of Population')

#plt.xlim(29,35)
plt.legend()
plt.show()
In [140]:
plt.figure(figsize=(13,6))

t1,s1,i1 = Ode_2D_solve_Sample(sfunc,ifunc,initial_s,initial_i,dt,5,Betas)

plt.plot(t[:1600],i[:1600],c='g',linewidth=0.6,label='original data Infectious')
plt.plot(t[1600:],i[1600:],c='r',linewidth=0.6,label='original future data Infectious')

#plt.plot(t[:1600],s[:1600],c='g',linewidth=0.6,label='original data Susceptible')
#plt.plot(t[1600:],s[1600:],c='r',linewidth=0.6,label='original future data Suceptible')



for count in range(0,100):
    t1,s1,i1 = Ode_2D_solve_Sample(sfunc,ifunc,initial_s,initial_i,dt,5,Betas)

    plt.plot(t1+16,i1,c='b',alpha = 0.1,linewidth=0.3)
    #plt.plot(t1+16,s1,c='b',alpha = 0.1,linewidth=0.3)

plt.xlabel('time in days/months/mucci')
plt.ylabel('Proportion of Population')

plt.xlim(15,20)
plt.ylim(0.3,0.6)
plt.legend()
plt.show()

lets look at real life data

In [143]:
Cases, Deaths, Populations = pf.Covid_Data("Platte County", "MO")
Counties = pf.Covid_Data_Keys('MO')
In [149]:
plt.plot(Cases[1],label='Confirmed Case')
plt.show()
In [150]:
len(Cases[1])
Out[150]:
277
In [152]:
for C in Counties:
    Cases, Deaths, Populations = pf.Covid_Data(C, "MO")
    if C != "Statewide Unallocated":
    
        plt.plot(Cases[1]/Populations, c='b',alpha=0.2)
    
plt.show()

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-152-b1a8796878f7> in <module>
      1 for C in Counties:
----> 2     Cases, Deaths, Populations = pf.Covid_Data(C, "MO")
      3     if C != "Statewide Unallocated":
      4 
      5         plt.plot(Cases[1]/Populations, c='b',alpha=0.2)

~/University Work/MAT217/Pythonic_files.py in Covid_Data(county, state)
     82     FIPS = dataC[(dataC.county_name==County)&(dataC.state==State)].county_fips_code.values[0]
     83 
---> 84     PopulationEstimate = dataX[(dataX.STATE==int(FIPS//1000))&(dataX.COUNTY==int(FIPS - FIPS//1000*1000))].POPESTIMATE2019.values[0]
     85 
     86 

IndexError: index 0 is out of bounds for axis 0 with size 0
In [153]:
3194*277
Out[153]:
884738
In [ ]: