SIR Model - The Flue Season - Dynamic Programming

# The SIR Model (susceptible, infected, and recovered) model is a common and useful tool in epidemiological modelling.

# In this post and in future posts I hope to explore how this basic model can be enriched by including different population groups or disease vectors.

# Simulation Population Parameters:
  # Proportion Susceptible
  Sp = .9

  # Proportion Infected
  Ip = .1

  # Population
  N = 1000

  # Number of periods
  r = 200

  # Number of pieces in each time period.
  # A dynamic model can be simulated by dividing each dynamic period into a sufficient number of discrete pieces.
  # As the number of pieces approaches infinity then the differences between the simulated outcome and the outcome achieved by solving the dynamic equations approaches zero.
  np = 1

# Model - Dynamic Change
  DS = function() -B*C*S*I/N
  DI = function() (B*C*S*I/N) - v*I
  DZ = function() v*I
  # I is the number of people infected, N the number of people in total, S is the number of people susceptible for infection, and Z is the number of people immune to the infection (from already recovering from the infection).

# Model Parameters:
  # Transmition rate from contact with an infected individual.
  B = .2
  # Contact rate.  The number of people that someone becomes in contact with sufficiently to recieve transmition.
  C = .5
  # Recovery rate. Meaning the average person will recover in 20 days (3 weeks).
  # This would have to be a particularly virolent form of the flu (not impossible at all).
  v = .05

# Initial populations:

  # Sesceptible population, Sv is a vector while S is the population values as the current period
  Sv = S = Sp*N

  # Infected, Iv is a vector while I is the population values as the current period
  Iv = I = Ip*N

  # Initial immunity.
  Zv = Z = 0

# Now let's how the model works.
  # Loop through periods
  for (p in 1:r) {
    # Loop through parts of periods
    for (pp in 1:np) {
   
      # Calculate the change values
      ds = DS()/np
      di = DI()/np
      dz = DZ()/np
 
      # Change the total populations
      S = S + ds
      I = I + di
      Z = Z + dz
   
      # Save the changes in vector form
      Sv = c(Sv, S)
      Iv = c(Iv, I)
      Zv = c(Zv, Z)
    }
  }

# ggplot2 generates easily high quality graphics
require(ggplot2)

# Save the data to a data frame for easy manipulation with ggplot
mydata = data.frame(Period=rep((1:length(Sv))/np,3), Population = c(Sv, Iv, Zv), Indicator=rep(c("Uninfected", "Infected", "Recovered"), each=length(Sv)))

# This sets up the plot but does not actually graph anything yet.

p = ggplot(mydata, aes(x=Period, y=Population, group=Indicator))  

 # This graphs the first plot just by the use of the p command.
 # Adding the geom_line plots the lines changing the color or the plot for each indicator (population group)
 p + geom_line(aes(colour = Indicator)) + ggtitle("Flu Season")

 # Save initial graph:
 ggsave(file="2013-05-14flu.png")

 # Let's do some back of the envelope cost calculations.
 # Let's say the cost of being infected with the flu is about $10 a day (a low estimate) in terms of lost productivity as well as expenses on treatment.
 # This amounts to:
 sum(Iv/np)*10
 # Which is a cost of $165,663.40 over an entire flu season for the thousand people in our simulated sample.
 # Or about $165 per person.

 # Imagine if we could now do a public service intervention.
 # Telling people to wash their hands, practice social distancing, and avoid touching their noses and eyes, and staying at home when ill.
 # Let's say people take up these practices and it reduces the number of potential exposure periods per contact by half.
 C = .25

 # ....

 p + geom_line(aes(colour = Indicator)) + ggtitle("Flu Season with Prevention")

  # Save initial graph:
 ggsave(file="2013-05-14flu2.png")

 # ....

  sum(Iv/np)*10
 # Which is a cost of $76,331.58 over an entire flu season  for the thousand people in our simulated sample or about 76 dollars per person.

 # The difference in costs is about 89 thousand dollars for the whole population or on average 89 per person.  The argument is therefore, so long as a public service intervention that reduces personal contact costs less than 89 thousand dollars for those 1000 people, then it is an efficient intervention (at least by the made up parameters I have here).

The Power of the Evaluate, Parse, Paste Combination

# See comments below.

# One of the powerful features of Stata which I have missed the most when working with R is the absence of the Stata Macros that allow the user to construct bits of Stata code from anything and combine them into commands or variable names.

# I know many a programmer who has modest to little use of Stata might sneer at this ability since it seems to imply some kind of laziness or lack of precision.  However, I have found myself on many an occation forced to use inefficient structures in coding in R that could have easily been simplified in Stata.

# At last I have stumbled upon a solution!

# By combining the command eval(parse(text=paste("String Command"))) I am able to do exactly that I want.

# For example, let's say I wanted to construct a dataframe A with elements a through z which are populated by 100 random normal variables.

# First I need to construct a data frame.
A = data.frame(id=1:10)
# Now to populate it.
for (i in letters) eval(parse(text=paste("A$",i,"=rnorm(10)",sep="")))

# This could be simplified a bit.
teval = function(...) eval(parse(text=paste(...,sep="")))

A = data.frame(id=1:10)
# Now to populate it using almost an identical command to that above.
# Except now I need to use the dreaded double arrow assign command because we are attempting to assign from within a function.
for (i in letters) teval("A$",i,"<<-rnorm(10)")

for (i in letters[seq(1,26,2)]) teval("A$",i,"<<-runif(10)")

# Unfortunately, this command is still not as powerful as using Stata's macros.
# However, it is a lot closer.

Spatial Critter Swarming Simulation

# I am interested in how small bits of individualized instructions can create collective action.

# In this simulation I will give a single instruction to each individual in the swarm.

# Choose another individual who is not too close, then accelerate towards that individual.

# I also control momentum causing the previous movement and direction to only decay at a small rate.

# TO SEE Original Script


# Critters are initially distributed randomly on a 1 x 1 grid.

ncritters = 40

xypos = matrix(runif(ncritters*2),ncol=2)
plot(xypos, main="Critters are Initially Distributed Randomly"
          , xlab="X", ylab="Y")

# Now let's imagine that each critter has an ideal safe distance from each other critter.

safe.dist = .3

critter.speed = .001

# If another critter is not at that safe distance than the critter will move towards the closest nearby critter.

# Let's see how this works.

# First let's check how close each critter is to each other critter.
# We will accomplish this by going through each critter and checking how far away each other critter is.
distances = NULL
for (i in 1:ncritters) distances = rbind(distances, apply((xypos[i,]-t(xypos))^2,2,sum))

# In order to prevent critters from always chasing whatever is closest to them (and themselves) we drop anything which is closer than the safe.distance.
distances[abs(distances)closest =  matrix(1:ncritters, ncol=ncritters, nrow=ncritters)[apply(abs(distances), 1, order)[1,]]
  # The apply command will apply the order command to each row whiel the [1,] selects only the critter that is closes.

# Plot the
plot(xypos, xlab = "X", ylab = "Y")
for (i in 1:ncritters) arrows(x0=xypos[i,1], y0=xypos[i,2],
                              x1=xypos[closest,][i,1],
                              y1=xypos[closest,][i,2],
                              length=.1)

# This calculates the difference between the current position of each critter and that of the closest critter.
ab = xypos-xypos[closest,]

# To see how this is calculated, see my previous post simulating a werewolf attack.

# Now calculate the difference in the horizontal and vertical axes that the critters will move as a projection into the direction of the closest critter outside of the safe zone.
a.prime = critter.speed/(1 + (ab[,2]^2)/(ab[,1]^2))^.5
b.prime = (critter.speed^2-a.prime^2)^.5

# This corrects the movement to ensure that the critters are flying at each other rather than away from each other.
movement = cbind(a.prime * sign(ab[,2]), b.prime * sign(ab[,1]))
between = function(xy1,xy2,point) (point>xy1&pointxy2&pointmovement = movement*(-1)^between(xypos,xypos[closest,], xypos-movement)

# Set the new xypos
xypos1 = xypos+movement

points(xypos1, col="red")

# ------------------------------------------------------
# Let's turn this into an animation.

library(animation)

# loopnum = 100; ncritters=40; inertia = .5; show.grid=T; ani.pause=F; plot.fixed=F; plot.centered=F; brownian = F; arrow = T
flocking <- ani.pause="F," arrow="T)" brownian="F," function="" inertia=".5," loopnum="100," ncritters="40," p="" plot.centered="F," plot.fixed="F," show.grid="T,">
  # Generate xy initial positions.
  # xypos will hold the current critter position while
  # xypos0 will hold the position of the critters the previous time.
  xypos = xypos0 = matrix(runif(ncritters*2),ncol=2)-.5

  movement0 = 0

  # Loop though all of the loops.
  for (i in 1:loopnum) {

  # This specifies the range to be graphed.
  if (plot.fixed) rangex=rangey = -.5:.5
  if (!plot.fixed) {
    rangex = c(min(xypos[,1]), max(xypos[,1]))
    rangey = c(min(xypos[,2]), max(xypos[,2]))
  }

  #  This handles the grid size when
  if (plot.centered&!plot.fixed) {
    rangex=c(-max(abs(xypos[,1])), max(abs(xypos[,1])))
    rangey=c(-max(abs(xypos[,2])), max(abs(xypos[,2])))
  }

  # This centers the plot at the middle (0,0) if the plot width is also set to be fixed.
  if (plot.centered&plot.fixed) rangex=rangex-mean(rangex)
  if (plot.centered&plot.fixed) rangey=rangey-mean(rangey)

  # Draw critters
  plot(xypos, main="Swarming Animation", xlab="X", ylab="Y", axes=F, ylim=rangey, xlim=rangex, type="p")
  # Draw arrows.  The start of the arrows is the previous periods location.
  if (arrow&i>1) arrows(x0=xypos0[,1],y0=xypos0[,2],x1=xypos[,1],y1=xypos[,2], length = .1)

  # Show the grid in the background.
  if (show.grid) {
    abline(v=seq(-10,10,.1))
    abline(h=seq(-10,10,.1))
  }

  # Show the origin
  text(0,0, "(0,0)")

  # Calculate each critters distance from each other
  distances = NULL
  for (i in 1:ncritters) distances = rbind(distances, apply((xypos[i,]-t(xypos))^2,2,sum))

  # Drop those within the safe zone.
  distances[abs(distances)
  # This selects the critter closest to the selected critter.
  closest =  matrix(1:ncritters, ncol=ncritters, nrow=ncritters)[apply(abs(distances), 1, order)[1,]]

#   distances[as.apply(!apply(distances, 1, is.na),1,sum)==0,]=0

  # As done above
  ab = xypos-xypos[closest,]

  a.prime = critter.speed/(1 + (ab[,2]^2)/(ab[,1]^2))^.5
  b.prime = (critter.speed^2-a.prime^2)^.5

  movement = cbind(a.prime * sign(ab[,2]), b.prime * sign(ab[,1]))

  between = function(xy1,xy2,point) (point>xy1&pointxy2&point
  movement = movement*(-1)^between(xypos,xypos[closest,], xypos-movement)

  movement[is.na(movement)]=0

  movement0 = movement0*inertia + movement

  # This fancy dodad allows half of the change in movement to be due to random variation.
  if (brownian) movement0=movement0+matrix(rnorm(ncritters*2),ncol=2)*critter.speed/2

  # Set the previous round's xy position to be equal to the current round's.
  xypos0 = xypos

  # Update the current round's.
  xypos = xypos+movement0

  # This is only used in the event that the animate package is in use.
  if (ani.pause) ani.pause()
  }

}

# This generates a GIF animation demonstrating smoothly how these GIFs can be incorper
  ani.options(ani.width=400, ani.height=400, interval=.1)

# You must have imagemagick installed for this to work.
  saveGIF(flocking(300,100,.999, ani.pause=T), movie.name = "Swarming.gif", replace=T)

# Here are two different graphs generated by the previous command(though the one on the bottom uses 200 frames while the one on the top uses 300)

# Let's see how this works.
flocking()
flocking(400,100,.99)
flocking(400,100,.99, plot.fixed=T)

Quandl Package - 5,000,000 free datasets at the tip of your fingers!

# Yes, you read that correctly and no Quandl (http://www.quandl.com/) did not pay me anything.

# Quandl is a new database management tool which seeks to become the place to find datasets.  They boast of having over 5x10^6 data sets available though after examining them, I have decided that they are not entirely what everybody might think of as data sets.  That is, each unique indicator is considered an independent data set.  This helps them to seem to have a ginormous quantity of data sets.

# That said, they are not wrong in calling each indicator its own data set since much of their data, like financial data or government data is collected by disjoint teams.  The scope of their ambition is fantastic yet it is doable and frankly someone needed to do it.

# Currently, data seekers can access the Inter-University Consortium for Political and Social Research (IPCSR).  This great resource is composed mostly of cross section and panel data sets which are great for much analysis but IPCSR resricts access to data to member universities.  In addition, the kind of data that Quandl is indexing is a lot of data that would not show up on IPCSR database.  In addition, Quandl is integrating an automated structure that will be self-updating.

# For an example of how Quandl is a good step ahead of the game take a look at this search quiery:

http://www.quandl.com/search/lansing,%20michigan

# In this search, I searched out Lansing, Michigan where I live and returned results of data for the last decade or earlier up to today from sources such as the Federal Reserve and the US Energy Information Administration.

http://www.icpsr.umich.edu/icpsrweb/ICPSR/studies?q=Lansing%2C+Michigan&permit%5B0%5D=AVAILABLE

# In constrast when queirying ICPSR, I found a few databases listed but they were historical databases that spanned back generally between 30 and 70 years.  That said both sources could provide valuable information depending upon what I am interested in modeling.

# Quandl is very clever for a number of reasons.  One of these reasons is that they have simultaneously released 8 software packages that can be used in a number of statistical packages such as R, Stata, and Excel.

# In order to demonstrate the use of Quandl I will grab a few data sets from the Lansing quiery drawn from the Federal Reserve.

install.packages("Quandl")
library(Quandl)

# Employment numbers (thousands of people") for Lansing, Michigan
NonFarm = Quandl("FRED/LANS626NAN")
CivLaborForce = Quandl("FRED/LANS626LFN")
PerCapitaIncome = Quandl("FRED/LANS626PCPI")

# Now let's combine the data so that we can related data values.
Labor = merge(NonFarm, CivLaborForce, by="Date")
Combined = merge(Labor, PerCapitaIncome, by="Date")
colnames(Combined) = c("Date", "NonFarm", "CivLaborForce", "PerCapitaIncome")
  # Notice that though our data had many more data points, the default option of merge only keeps data that exists in both data sets.  In this case, it is per capital income that has the least number of data points.

# Let's see if we can predict income as a function of employment:
summary(lm(PerCapitaIncome~NonFarm+CivLaborForce, data=Combined))

# Our naive prediction as a result of this is that as the Civilian Labor Force increases, wages rise.  This is of course a naive example ignoring completely issues of causation and endogeneity not to mention probable random walks and other challenging features of this kind of data.

# The overall take away though, should be "cool", I think.  Maybe this data bank does not provide information currently on many issues of interest to those looking for data.  But it does make things easier and self-updating, which are great features.

Concerto Simple Demonstration Test v4.b

# This is my first attempt to post a example of a test using Concerto 4 beta.

# This is an extremely simple test.  Too simple to be useful to most people yet I believe it is didactically helpful.

---

The Following is the HTML code for Template 9. It is the form that the user first inputs answers to the three questions.

Template Description/Name: FixedSimpleQuestions

<p>Please Answer:</p>

<!--This displays each question.  The important thing to notice is the name of each text box, txt_1, txt_2, and txt_3-->
<p>1. 1+1=<input name="txt_1" type="text" /></p>

<p>2. 2+2=<input name="txt_2" type="text" /></p>

<p>3. 3+3=<input name="txt_3" type="text" /></p>

<p><input name="done" type="button" value="done" />​</p>

---

The Following is the HTML code for Template 10. It is the form that the user receives feedback from.

Template Description/Name: FixedSimpleResults

<p>You got questions:</p>

<!-- Notice the values in {{}}.  These are the parameters passed to this template -->

<p>1. {{score1}}</p>

<p>2. {{score2}}</p>

<p>3. {{score3}}</p>

---

# The Following is the R code that the test actually uses to run.

# We can see how R controls the flow of the user interface.

# Concerto Simple Test 1

# We call template 9, FixedSimpleQuestions, to be shown.  It can either be referenced by the template number of the name of the template.
items = concerto.template.show("FixedSimpleQuestions")

# Check if the items results are correct.  We can see that the responses to the items from the first template are saved within the list object item with the names txt_1, txt_2, and txt_3.
TF1 = items$txt_1==2
TF2 = items$txt_2==4
TF3 = items$txt_3==6

# Show the results to the user by feeding our item values as parameters back into the second template.
concerto.template.show("FixedSimpleResults", params=list(score1=TF1, score2=TF2, score3=TF3))

# I hope this gives you a feel for how easy it is to use Concerto to handle web based inputs.

# A powerful feature of Concerto is that it has functions that and query and modify tables using My SQL.

# More in future posts.

# I used an HTML encoder to translate my html code to post into blogger.

A Command for Randomly Creating Sets of Elements from a Vector

# This command pairs (of arbitrary length) together a list of items randomly.
# I think it is useful for simulations matching individuals together.

pairup = function(x, ncol=2, unmatched.self=T) {

  # Calculate the length of the x vector.
  xleng = length(x)

  # This checks if the input vector is a scalar.
  # If it is then it assumes that the number of that scalar is to be the last number is the seqence 1:x
  if (xleng==1) x = 1:(xleng=x)
  # This command is pretty weird.  I combined two commands to make a single line command.
  # It sets the xleng equal to x and x equal to the sequence 1:x

  # Calculate the length of each column that will be created.
  hleng = floor(xleng/ncol)

  # Randomize x
  x = x[order(runif(xleng))]

  pair = x[1:hleng]
  for (i in 2:ncol) pair = cbind(pair,x[((i-1)*hleng+1):((i)*hleng)])

  # If there is a odd number of xs then this will match the remaining unmatched x with itself if unmatched.self is T.
  if ((unmatched.self)&(floor(xleng/ncol)!=xleng/ncol)) pair=rbind(pair, c(x[ncol*hleng+1]))
  if ((!unmatched.self)&(floor(xleng/ncol)!=xleng/ncol)) print(x[!(x %in% pair)])

  return (pair)
}

pairup(1:10)
pairup(99)
pairup(letters)
pairup(letters, ncol=3)
pairup(letters, ncol=3, unmatched.self=F)

A Dynamic Simulation of a Zombie Apocalypse

# Zombies vs Humans Agent Based Simulation (the r script file in case blogger mangled my code)

# I also wrote a Spatial Simulation of a zombie infection in Stata previously.

# This simulation is a simple repeated matching simulation in which each agent is matched with a random different agent.

# I believe this is an appropriate way of modelling a human zombie exchange as portrayed in the movies.  Generally speaking, each encounter can be thought of as a probabilistic draw in which either the human becomes zombified or the zombie is permanently killed.

# Agents all start out as humans (except a random percent which are initially zombified).  Each human is

defined on a scale of 0 to 90 percentile with increments of 10 which reflect the chance of that human vanquishing a zombie that the human encountered or being in turn vanquished.

# If a human is vanquished a zombie is added to the zombie population.

# number of humans (approximated due to rounding issues)
start.humans = 100000

# matrix of human types
htypes = seq(0,90,10)

# Frequency of each time of human from 0 to 90 percentile.
# This number is only relative to the scale of the other frequencies so long as it is positive.
# Thus if only type of human had a 5 then it would be 5 times more likely than a type of human with a frequency of 1.
freq = rep(1,length(htypes))

# frequency is the most important parameter choice in the model as will be seen below.

# Bind the information into a single matrix.
human.types = cbind(htypes,freq)

# Now we calculate what percentage of our start.humans are each type.
perc = round((freq/sum(freq))*start.humans)

# Finally we generate our moving things data.
# Initially the only thing moving is humans.
walking.things = rep(htypes,perc)

# Looking good
walking.things

# Some percentage of the humans are initially infected.
infected.per = .025

# Calculate the  number initially that become infected.
nselected = round(start.humans*infected.per)

# Now we randomly select the humans infected initially.
initial.zombies = sample(1:length(walking.things), nselected)

walking.things=c(walking.things[-initial.zombies], rep(-77,nselected))
           
# -77 Is the number for a zombie.
   
# After the intial infection phase peopole get their guns out and start acting defensively.
walking.things

# Percent zombies
perc.zombies = mean(walking.things==-77)

# Total population (living and dead)
nthings.vector = nthings = length(walking.things)

# Count the number of zombies
nzombies = sum(walking.things==-77)

# Count the number of humans
nhumans = sum(walking.things!=-77)

  nhumans0  = sum(walking.things==00)
  nhumans10 = sum(walking.things==10)
  nhumans20 = sum(walking.things==20)
  nhumans30 = sum(walking.things==30)
  nhumans40 = sum(walking.things==40)
  nhumans50 = sum(walking.things==50)
  nhumans60 = sum(walking.things==60)
  nhumans70 = sum(walking.things==70)
  nhumans80 = sum(walking.things==80)
  nhumans90 = sum(walking.things==90)


# This command pairs up a vector.  It is used to match humans with zombies.
pairup = function(x, unmatched.self=T) {

  # Calculate the length of the x vector.
  xleng = length(x)

  # This checks if the input vector is a scalar.
  if (xleng==1) x = 1:(xleng=x)

  # Half the length of x rounded down.
  hleng = floor(xleng/2)

  # Randomize x
  x = x[order(runif(xleng))]
  pairs = cbind(x[1:hleng],x[(hleng+1):(2*hleng)])

  # If there is a odd number of xs then this will match the remaining unmatched x with itself if unmatched.self is T.
  if ((unmatched.self)&(xleng/2!=hleng)) pairs=rbind(pairs, c(x[2*hleng+1]))

  return (pairs)
}

#
max.rounds = 45

# Let's start the simulation:
n = 1
while  (nzombies[n]>0 & nhumans[n]>0 & n  n = n+1

  # This calls the previously defined function pairup to match two different individuals together.
  # This matches them by position in the walking.things vector.
  encounter=pairup(nthings)

  # This assigns to the matrix the values
  types = cbind(walking.things[encounter[,1]],walking.things[encounter[,2]])

  # Create a vector of terminated or zombified things
  conflict = types*0
     # 0 Unresolved
     # 1 Zombified
     # 2 Permenent Death
     # 3 No conflict
     # 4 win conflict
 
   # This code will check if a zombie is in the right column and human in the left.
   hvz = (types[,2]==-77)&(types[,1]>=0)
     # If so, the human and zombie places will be switched.
     types.temp = types
     types[hvz,1]=types.temp[hvz,2]
     types[hvz,2]=types.temp[hvz,1]
   
     encounter.temp = encounter
     encounter[hvz,1]=encounter.temp[hvz,2]
     encounter[hvz,2]=encounter.temp[hvz,1]
   
   # Zombie encounters human
   zvh = (types[,1]==-77)&(types[,2]>=0)
     # Calculate the win count of the conflict
     win.zvh = (runif(sum(zvh))>types[zvh,2]/100)
   
     # Translate a zombie win onto the conflict map
     conflict[zvh,1][win.zvh]=4
     conflict[zvh,2][win.zvh]=1  
   
     # Translate a human win onto the conflict map
     conflict[zvh,1][!win.zvh]=2
     conflict[zvh,2][!win.zvh]=4

     # Resolve non-conflict. Zombies don't fight zombies and humans don't fight humans.
     conflict[types[,1]==types[,2],] = 3
     conflict[(types[,1]>=0)&(types[,2]>=0),] = 3
   
   # Finally, adjust the walking.things vector to adjust for the changes.
   # Zombify some
   walking.things[encounter[conflict==1]] = -77
   
   # Remove others
   walking.things=walking.things[-encounter[conflict==2]]

  # Store stats
  # Percent zombies
  perc.zombies = c(perc.zombies, mean(walking.things==-77))

  # Total population (living and dead)
  nthings = length(walking.things)
  nthings.vector = c(nthings.vector, nthings)

  # Count the number of zombies
  nzombies = c(nzombies, sum(walking.things==-77))

  # Count the number of humans and save them to vectors
  nhumans = c(nhumans,  sum(walking.things!=-77))

  nhumans0  = c(nhumans0,   sum(walking.things==0))
  nhumans10 = c(nhumans10,  sum(walking.things==10))
  nhumans20 = c(nhumans20,  sum(walking.things==20))
  nhumans30 = c(nhumans30,  sum(walking.things==30))
  nhumans40 = c(nhumans40,  sum(walking.things==40))
  nhumans50 = c(nhumans50,  sum(walking.things==50))
  nhumans60 = c(nhumans60,  sum(walking.things==60))
  nhumans70 = c(nhumans70,  sum(walking.things==70))
  nhumans80 = c(nhumans80,  sum(walking.things==80))
  nhumans90 = c(nhumans90,  sum(walking.things==90))

}

# Count the number of rounds completed.
nrounds.completed = length(nhumans0)

plot(c(1,nrounds.completed), c(0,max(nhumans0,nzombies)), type="n",
     ylab="Population", xlab="Round\n*Note: The number at the end of each line is the probability that an individual\nin this population group will kill a zombie when encountering one. Z is the zombie population"
     , main="Population During a Zombie Attack\nWith a well armed human population a zombie apocalypse is easily prevented")
for (i in seq(0,90,10)) {
  lines(get(paste("nhumans",i,sep="")))
  text(nrounds.completed+1, get(paste("nhumans",i,sep=""))[nrounds.completed], i)
}
lines(nzombies, lwd=3)
  text(nrounds.completed+1, nzombies[nrounds.completed], "Z")

# However, this equal proportion of highely effective zombie killers to very ineffective zombie killers is nonrepresentational of a typical zombie movies or games.

# I will instead rerun the simulation with a larger percentage of low ability humans.
htypes = seq(0,90,10)
freq = 10:1

# .... using same code as above but with new human population proportions

nrounds.completed = length(nhumans0)

plot(c(1,nrounds.completed), c(0,max(nhumans0,nzombies)), type="n",
     ylab="Population", xlab="Round\n*Note: Even the best trained individual will be overwhelmed eventually\nif there are too many easy zombie victims. Z is the zombie population"
     , main="Population During a Zombie Attack\nA poorly armed population is ill-equiped to survive a zombie attack")
for (i in seq(0,90,10)) {
  lines(get(paste("nhumans",i,sep="")))
  text(nrounds.completed+1, get(paste("nhumans",i,sep=""))[nrounds.completed], i)
}
lines(nzombies, lwd=3)
  text(nrounds.completed+1, nzombies[nrounds.completed], "Z")

# Let's try one more variant with still more weak humans but a little better ratios.

htypes = seq(0,90,10)
freq = seq(5,1, length.out=10)

# .... using same code as above but with new human population proportions

nrounds.completed = length(nhumans0)

plot(c(1,nrounds.completed), c(0,max(nhumans0,nzombies)), type="n",
     ylab="Population", xlab="Round\n*Note: In this scenario only the top 10 to 20% most effective zombie killers survive."
     , main="Population During a Zombie Attack-When the population of dangerous humans\n is sufficiently large, it is possible humanity survives, just barely")
   
for (i in seq(0,90,10)) {
  lines(get(paste("nhumans",i,sep="")))
  text(nrounds.completed+1, get(paste("nhumans",i,sep=""))[nrounds.completed], i)
}
lines(nzombies, lwd=3)
  text(nrounds.completed+1, nzombies[nrounds.completed], "Z")

# Overall conclusion? Mandating zombie defense classes is the only way to be certain humanity will survive.