<jupyter_start><jupyter_text>Reed-Frost Model in R*Author*: <NAME>*Date*: 2018-10-01 - `xcurr` is the number of suspectible individuals at step j.- `ycurr` is the number of infected individuals at step j.- `q` is the probability that a suspectible individual escapes infection.- `ynext` is the number of infected persons at step j + 1.Everyone who is infected at step j becomes immune (or dies) at step j + 1. Thus at step j + 1, the number of susceptibles is `xcurr-ynext`. The following function gives the probability that the number of infected persons at $j+1$ is `ynext`.<jupyter_code>reed_frost <- function(xcurr, ycurr, q, ynext) {
  pescape <- q^ycurr
  pynext <- choose(xcurr, ynext) *
      (1 - pescape)^ynext *
      pescape^(xcurr - ynext)
  pynext
}<jupyter_output><empty_output><jupyter_text>The number of infected people at a time step can be no greater than the number of suspectibles at the current time step.So a sanity check for our code would for this number to be 0.<jupyter_code>reed_frost(10, 5, 0.5, 11)<jupyter_output><empty_output><jupyter_text>If the probability of escaping an infection is 1, then the probability that the number of infected at next time is 0 is 1.Similarly at the other extreme, if probability of escaping aninfection is 0, then every susceptible becomes infected and the number of infected at next time is `xcurr`.<jupyter_code>## Probability that ynext = 0 should be 1.
## 0 everywhere else.
reed_frost(10, 5, 1, 0:10)

## Probability that ynext = 10 should be 1.
## 0 everywhere else.
reed_frost(10, 5, 0, 0:10)<jupyter_output><empty_output><jupyter_text>That works as expected. Now for various values of q and \(y_curr\), we can check the behavior of the model.First we will load the needed libraries.<jupyter_code>library(ggplot2)
library(dplyr)
pescape <- 0.9
ycurrs <- c(20, 40, 60, 80, 100)
probs <- purrr::map_dfc(ycurrs,
                             ~ reed_frost(xcurr = 100,
                                          ycurr = .x,
                                          q = pescape,
                                          ynext = 0:100))
colnames(probs) <- ycurrs
probs %>%
    tidyr::gather(ycurr, prob, factor_key = TRUE) %>%
    ggplot(aes(rep(0:100, length(ycurrs)),
               log(prob),
               col = ycurr)) +
    geom_point() +
    xlab("ynext") +
    ylab("log(prob)")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09873949579831932, 'fraction_numerical': 0.029831932773109245, 'mean_word_length': 3.14808362369338, 'pattern_counts': {'":': 0, '<': 25, '<?xml version=': 0, '>': 20, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '480895', 'license': 'MIT', 'n_tokens_mistral': 809, 'n_tokens_neox': 731, 'n_words': 277, 'path': 'notebooks/reedfrost/r.ipynb', 'repo_name': 'epirecipes/epiRecipes', 'size': '2389'}