Sunday, April 28, 2013

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.

15 comments :

  1. Fascinating. However, wouldn't a well-armed population increased the likelihood of human-human killings? For having watched far too many zombie movies, this danger can often prevail...

    ReplyDelete
    Replies
    1. Excellent point. I had assumed that humans and zombies collaborated among themselves. I had resisted modelling human/human conflicts because under any kind of randomness of murder scenario nearly the entire human population will be wiped out no matter what happens. I guess I would have to model also replacement humans (child-birth) in order to hope to change the overall outcome.

      Delete
  2. This post has earned my subscription. :)

    ReplyDelete
  3. Hi Francis
    what is this line doing:
    x <- order="" p="" runif="" x="" xleng="">
    it does not work in R
    (Error: unexpected symbol in:
    " # Randomize x
    x <- order="" p")

    thanks,
    Aga

    ReplyDelete
    Replies
    1. Hi Aga, thanks for writing. Sometimes blogger mangles my code in an attempt to interpret html tags which are not there. Try downloading the attached r script file with the link at the beginning of the document. Are you still having problems? Here is the link as well:
      R Script

      Delete
  4. I had done some simple compartmental epidemic models of the zombie apocalypse for a workshop I taught - amusingly, they have remarkably similar findings. Zombie defense classes are the only way :)

    ReplyDelete
    Replies
    1. Awesome :) Looks like a developing interdisciplinary field. At Michigan State we now have a course on How to survive a Zombie Apocalypse

      Delete
  5. This comment has been removed by the author.

    ReplyDelete
  6. One of my friends sent me a link to this web-comic. Perhaps I should do a follow up simulation: Dinosaur vs Zombies vs Humans

    Dinosaur Comics

    ReplyDelete

    ReplyDelete
  7. I have a question that relates to all zombie apocalypse scenarios. Why does an encounter with a zombie lead to either human death or zombie creation? Surely as the number of zombies increased, they would just totally consume any human victims thus creating no new zombies. Eventually, the zombies would rot away and people could come out of whatever caves they've been hiding in. Am I missing something?

    ReplyDelete
    Replies
    1. A very good point. There seems to be some yet unidentified law of biology that allows zombies to only partially rot such that they seem to be capable of lingering effectively forever until killed by human means. This is the struggle that any zombie movie or show seems to face. If you are smart, you build a bunker. But then you are safe from zombies. Now what? Generate plot twists by some means.

      Delete
    2. Right, but that still doesn't answer the total consumption of new victims problem...

      Delete
  8. This comment has been removed by a blog administrator.

    ReplyDelete