Week 11: SIR Model and Covid Data

import numpy as np
import matplotlib.pyplot as plt
import Pythonic_files as pf
import seaborn as sns
%config InlineBackend.figure_format='retina'

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

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()

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()

estimate_beta=np.mean(Betas)

initial_s = s[-1]
initial_i = i[-1]

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

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

t

array([0.000e+00, 1.000e-02, 2.000e-02, ..., 2.197e+01, 2.198e+01,
       2.199e+01])

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)

sns.distplot(Betas)
print(np.mean(Betas))

0.48836515898184174

t[1600]

16.0

initial_s = s[1601]
initial_i = i[1601]
estimate_beta = np.mean(Betas)

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()

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()

15

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
        
    
    

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()

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

Cases, Deaths, Populations = pf.Covid_Data("Jackson County (including other portions of Kansas City)", "MO")
Counties = pf.Covid_Data_Keys('MO')

plt.figure(figsize=(13,6))
plt.plot(Cases[0],label='Confirmed Case exact')
plt.plot(Cases[1],label='Confirmed Case rolling')
plt.legend()


plt.show()

plt.figure(figsize=(13,6))
plt.plot(Deaths[0],label='Confirmed Death exact')
plt.plot(Deaths[1],label='Confirmed Death rolling')
plt.legend()


plt.show()

plt.figure(figsize=(13,6))
plt.plot(Cases[0]/Populations,label='Confirmed Case exact')
plt.plot(Cases[1]/Populations,label='Confirmed Case rolling')
plt.legend()


plt.show()

plt.figure(figsize=(13,6))

County="Jackson County (including other portions of Kansas City)"
State = "MO"
Cases, Deaths, Populations = pf.Covid_Data(County, State)

plt.plot(Cases[1]/Populations,label=County)

County="Buchanan County"
Cases, Deaths, Populations = pf.Covid_Data(County, State)
plt.plot(Cases[1]/Populations,label=County)

County="Taney County"
Cases, Deaths, Populations = pf.Covid_Data(County, State)
plt.plot(Cases[1]/Populations,label=County)



plt.legend()


plt.show()

plt.figure(figsize=(13,6))

County="Jackson County (including other portions of Kansas City)"
State = "MO"
Cases, Deaths, Populations = pf.Covid_Data(County, State)

plt.plot(Deaths[1]/Populations,label=County)

County="Buchanan County"
Cases, Deaths, Populations = pf.Covid_Data(County, State)
plt.plot(Deaths[1]/Populations,label=County)

County="Taney County"
Cases, Deaths, Populations = pf.Covid_Data(County, State)
plt.plot(Deaths[1]/Populations,label=County)



plt.legend()


plt.show()

Lets take a look at all of the counties in MO, as see if there are outliers, or trends

plt.figure(figsize=(13,6))

State = "MO"
Counties = pf.Covid_Data_Keys(State)

for C in Counties:
    if C != 'Statewide Unallocated':
        Cases, Deaths, Populations = pf.Covid_Data(C, State)
        plt.plot(Cases[0]/Populations*100, c='b',linewidth=0.4, alpha =0.2)
        if np.max(Cases[0]/Populations*100)>=4:
            print(C)


plt.title('Percentage of Confirmed Cases for each county in '+State)
plt.ylabel('Percentage')
plt.xlabel('Days since Jan 22')
plt.show()

Perry County
Saline County
Sullivan County
Pemiscot County
Jasper County
Madison County
Nodaway County
New Madrid County
McDonald County
St. Francois County

plt.figure(figsize=(13,6))

State = "ND"
Counties = pf.Covid_Data_Keys(State)

for C in Counties:
    if C != 'Statewide Unallocated':
        Cases, Deaths, Populations = pf.Covid_Data(C, State)
        plt.plot(Cases[0]/Populations*100, c='b',linewidth=0.4, alpha =0.2)
        if np.max(Cases[0]/Populations*100)>=4:
            print(C)


plt.title('Percentage of Confirmed Cases for each county in '+State)
plt.ylabel('Percentage')
plt.xlabel('Days since Jan 22')
plt.show()

Burleigh County
Cass County
Stark County
Morton County
Grand Forks County
Mountrail County
Golden Valley County
Hettinger County
Logan County
Towner County
Dickey County
LaMoure County
Nelson County
Foster County
McIntosh County
Benson County
Sioux County
Renville County
Emmons County
Eddy County
Mercer County
McLean County
Ramsey County
Dunn County
Williams County

Can we look at the daily increases, how freaky its getting

plt.figure(figsize=(13,6))

State = "MO"
Counties = pf.Covid_Data_Keys(State)
for C in Counties:
    if C != 'Statewide Unallocated':
        Cases, Deaths, Populations = pf.Covid_Data(C, State)
        plt.plot(Cases[0][1:]/Populations*100-Cases[0][:-1]/Populations*100,'.',c='b', alpha =0.1,label='Confirmed Case exact')



plt.title('Percentage daily change of Confirmed Cases for each county in '+State)
plt.ylabel('Percentage')
plt.xlabel('Days since Jan 22')


plt.show()

Lets try forcast, lets only look at Buchanan County and stop it at 260

Cases, Deaths, Populations = pf.Covid_Data("Buchanan County", "MO")
Counties = pf.Covid_Data_Keys('MO')


plt.plot(Cases[1][0:250]/Populations*100)

[<matplotlib.lines.Line2D at 0x7fd988ba1d30>]

plt.plot(Cases[1][14:250]/Populations*100 - Cases[1][0:236]/Populations*100)

[<matplotlib.lines.Line2D at 0x7fd9d9740f40>]

Infectious=(Cases[1][9:250]-Cases[1][0:241])/Populations*100
Susceptible = (np.ones(len(Cases[1][9:250]))*Populations - Cases[1][9:250])/Populations*100

plt.plot(Infectious[60:])

#plt.plot(Cases[1][14:250])

[<matplotlib.lines.Line2D at 0x7fd9c8bc8730>]

Betas = []
for x in range(60,len(Infectious)):
    b = -(Susceptible[x]-Susceptible[x-1])/(Susceptible[x]*Infectious[x])
    Betas.append(b)

sns.distplot(Betas)

/opt/anaconda3/lib/python3.8/site-packages/seaborn/distributions.py:2551: FutureWarning: `distplot` is a deprecated function and will be removed in a future version. Please adapt your code to use either `displot` (a figure-level function with similar flexibility) or `histplot` (an axes-level function for histograms).
  warnings.warn(msg, FutureWarning)

<AxesSubplot:ylabel='Density'>

Infectious2=(Cases[1][9:277]-Cases[1][0:268])/Populations*100

plt.figure(figsize=(13,6))

Times = range(0,277)

plt.plot(Times[60:241],Infectious[60:241],c='g',linewidth=0.6,label='original data Infectious')
plt.plot(Times[240:268],Infectious2[240:268],c='r',linewidth=0.6,label='original future data Infectious')

[<matplotlib.lines.Line2D at 0x7fd990b8d610>]

def Ode_2D_solve_Sample(x0,y0,dt,T,Betas):
    
    def sfunc(s,i):
        return -beta*s*i
    def ifunc(s,i):
        return beta*s*i-gamma*i

    gamma = 1/9

    
    
    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
        
    
    

plt.figure(figsize=(13,6))

Times = range(0,277)

plt.plot(Times[60:241],Infectious[60:241]*Populations,c='g',linewidth=0.6,label='original data Infectious')
plt.plot(Times[240:268],Infectious2[240:268]*Populations,c='r',linewidth=0.6,label='original future data Infectious')

dt=1

initial_s=Susceptible[-1]
initial_i=Infectious[-1]

for count in range(0,200):
    t1,s1,i1 = Ode_2D_solve_Sample(initial_s,initial_i,dt,28,Betas)

    plt.plot(t1+240,i1*Populations,c='b',alpha = 0.3,linewidth=0.3)
    

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

plt.xlim(235,277)
plt.legend()
plt.show()