diff --git a/Exercise10_revised.pdf b/Exercise10_revised.pdf index e282db8..8d4bebe 100644 Binary files a/Exercise10_revised.pdf and b/Exercise10_revised.pdf differ diff --git a/exercise_10.py b/exercise_10.py new file mode 100644 index 0000000..d6c39bd --- /dev/null +++ b/exercise_10.py @@ -0,0 +1,78 @@ +## exercise 10 ## + +# 1 +#Load packages +import pandas +import scipy +import scipy.integrate as si +from plotnine import * + + +#Define custom function +def popGrowth(y,t0,r,K,N): + N=y[0] + dNdt=r*(1-N/K)*N + return [dNdt] +times=range(0,100) +NO=[.01] +#Set a pool of values for growth rate +growthRates=[-.1,.1,.4,.8,1.0] +#Dataframe for storing model output +store_growthRates=pandas.DataFrame({"time":times,"r1":0,"r2":0,"r3":0,"r4":0,"r5":0}) + +#Using a for loop to make my life easier +for i in range(0,len(growthRates)): + pars=(growthRates[i],100,10) + sim=si.odeint(func=popGrowth,y0=NO,t=times,args=pars) + store_growthRates.iloc[:,i]=sim[:,0] +print(store_growthRates) + +#Plot 1- pop size as function of time +g=(ggplot(data=store_growthRates) + +geom_line(store_growthRates,aes(x="time",y="r1"))+theme_classic() + +ylab("population size") + +geom_line(store_growthRates,aes(x="time",y="r2"),color="green") + +geom_line(store_growthRates,aes(x="time",y="r3"),color="red") + +geom_line(store_growthRates,aes(x="time",y="r4"),color="orange") + +geom_line(store_growthRates,aes(x="time",y="r5"),color="purple")) +print(g) + +#set a pool of values for carrying capacity +carryCap=[10,50,100] +#Dataframe for storing model output +store_carryCap=pandas.DataFrame({"time":times,"K1":0,"K2":0,"K3":0}) +NO = 1 +#for loop for Plot 2 +for i in range(0,len(carryCap)): + pars=(.2,carryCap[i],1) + sim=si.odeint(func=popGrowth,y0=NO,t=times,args=pars) + store_carryCap.iloc[:,i]=sim[:,0] +print(store_carryCap) + +#Plot 2- +k=(ggplot(data=store_carryCap) + +geom_line(store_carryCap,aes(x="time",y="K1"),color="blue")+theme_classic() + +ylab("population size") + +geom_line(store_carryCap,aes(x="time",y="K2"),color="green") + +geom_line(store_carryCap,aes(x="time",y="K3"),color="purple")) +print(k) + +#Set a pool of values for init pop size +initPop=[1,50,100] +#Dataframe for storing model output +store_initPop=pandas.DataFrame({"time":times,"N1":0,"N2":0,"N3":0}) + +#Using a for loop to make my life easier +for i in range(0,len(initPop)): + pars=(.1,50,initPop[i]) + sim=si.odeint(func=popGrowth,y0=initPop[i],t=times,args=pars) + store_initPop.iloc[:,i]=sim[:,0] +print(store_initPop) + +#Plot 3- effect of initial Pop size differences +p=(ggplot(data=store_initPop) + +geom_line(store_initPop,aes(x="time",y="N1"),color="orange")+theme_classic() + +ylab("population size") + +geom_line(store_initPop,aes(x="time",y="N2"),color="yellow") + +geom_line(store_initPop,aes(x="time",y="N3"),color="red")) +print(p) diff --git a/exercise_10_B.py b/exercise_10_B.py new file mode 100644 index 0000000..c6df0fc --- /dev/null +++ b/exercise_10_B.py @@ -0,0 +1,141 @@ +#Load packages +import pandas +import scipy +import scipy.integrate as si +from plotnine import * + +# function +def SIR (y,t0,beta,gamma): + S = y[0] + I = y[1] + R = y[2] + dS = -1*(beta*I*S) + dI = (beta*I*S)-(gamma*I) + dR = (gamma*I) + return dS, dI, dR + +# initial conditions +times = range(0,500) +NO = [999, 1, 0] + +# dataframe of gamma and beta values +data = [{'beta' : .0005, 'gamma' : .05}, + {'beta': .005, 'gamma': .5}, + {'beta': .0001, 'gamma': .1}, + {'beta': .00005, 'gamma': .1}, + {'beta': .0001, 'gamma': .05}, + {'beta': .0002, 'gamma': .05}, + {'beta': .0001, 'gamma': .06}, + {'beta': .0001, 'gamma': .9}, + {'beta': .0001, 'gamma': .5}, + {'beta': .0001, 'gamma': .25}, + {'beta': .0001, 'gamma': .1}, + {'beta': .0001, 'gamma': .05}, + {'beta': .0001, 'gamma': .01}, + {'beta': .0001, 'gamma': .001}, + {'beta': .0001, 'gamma': .0001}, + {'beta': .9, 'gamma': .0001}, + {'beta': .5, 'gamma': .0001}, + {'beta': .25, 'gamma': .0001}, + {'beta': .1, 'gamma': .0001}, + {'beta': .05, 'gamma': .0001}, + {'beta': .01, 'gamma': .0001}, + {'beta': .001, 'gamma': .0001}, + {'beta': .0001, 'gamma': .0001}, + {'beta': .00001, 'gamma': .0001}, + ] +my_data = pandas.DataFrame(data) + +# make lists to hold the results +mdi = [] +mdp = [] +pa = [] +ro = [] +b = [] +g = [] + +# start big loop here +for line in range(0,len(my_data),): + q = my_data.iloc[line]['beta'] + p = my_data.iloc[line]['gamma'] + params = (q, p) # make tuple + + b.append(params[0]) # append list + g.append(params[1]) + + # make dataframe + infection = pandas.DataFrame({"time":times,"S":0,"I":0,"R":0}) + + # sim shite + sim = si.odeint(func=SIR, y0=NO, t=times, args=params) + + # fill dataframe + infection.iloc[:,2]=sim[:,0] + infection.iloc[:,0]=sim[:,1] + infection.iloc[:,1]=sim[:,2] + + # calc max daily incidence + daily_incidence = [] + for i in range(0,len(infection),): + if infection.time[i]==0: + continue + else: + I = infection.iloc[i]['I'] + Iold = infection.iloc[i-1]['I'] + incidence = I-Iold + daily_incidence.append(incidence) + max_daily_incidence = max(daily_incidence) + mdi.append(max_daily_incidence) + + # calc max daily prevalence + daily_prev = [] + for i in range(0,len(infection),): + I = infection.iloc[i]['I'] + R = infection.iloc[i]['R'] + S = infection.iloc[i]['S'] + prev = I/(S+I+R) + daily_prev.append(prev) + max_daily_prev = max(daily_prev) + mdp.append(max_daily_prev) + + #calc percent affected over simulation- use last time step (499) + I= infection.iloc[499]['I'] + R= infection.iloc[499]['R'] + S= infection.iloc[499]['S'] + percent_affected = (I+R)/(S+I+R) + pa.append(percent_affected) + + # basic reproduction number initial SIR + beta = params[0] + gamma = params[1] + I= infection.iloc[0]['I'] + R= infection.iloc[0]['R'] + S= infection.iloc[0]['S'] + repo_number = (beta*(S+I+R))/gamma + ro.append(repo_number) + +# make a dataframe for results from all the lists +results = pandas.DataFrame( + {'beta' : b, + 'gamma' : g, + 'max_daily_incide' : mdi, + 'max_daily_prev' : mdp, + 'percent_affect' : pa, + 'repo_num' : ro}) + +print results +#need these to fill into a list or a dataframe +# need to put all this intoa bigger loop + +''' +* observations * +We noticed that if beta is held constant while gamma gets smaller the max incidence, daily prevalence, percent infection, +and reproductive number all rise. Thus a bigger gamma causes, things to get smaller If we hold gamma constant, as beta +gets bigger the max incidence, daily prevalence, percent affected, and reproductive number all rise. Thus a small beta, +causes things to get smaller To summarize, a high beta and a small gamma cause the disease to have a higher rate and +srpead of infection in the population. + +''' + + + diff --git a/simData.xlsx b/simData.xlsx new file mode 100755 index 0000000..045e86f Binary files /dev/null and b/simData.xlsx differ