Skip to content

Instantly share code, notes, and snippets.

@brshallo
Created October 21, 2021 06:14
Show Gist options
  • Save brshallo/c198c8c3bc4f84b97e8a97f41a8445da to your computer and use it in GitHub Desktop.
Save brshallo/c198c8c3bc4f84b97e8a97f41a8445da to your computer and use it in GitHub Desktop.
library(ggplot2)
library(dplyr)
library(purrr)

# Binomial distribution method
probs_survive <- 
  tibble(survivors = 1:16) %>% 
  mutate(prob = dbinom(16 - survivors, size = 18, prob = 0.5))

bind_rows(
  tibble(survivors = 0,
         prob = 1 - sum(probs_survive$prob)),
  probs_survive
) %>% 
  print() %>% 
  ggplot(aes(x = survivors, y = prob))+
  geom_col()
#> # A tibble: 17 x 2
#>    survivors       prob
#>        <dbl>      <dbl>
#>  1         0 0.000656  
#>  2         1 0.00311   
#>  3         2 0.0117    
#>  4         3 0.0327    
#>  5         4 0.0708    
#>  6         5 0.121     
#>  7         6 0.167     
#>  8         7 0.185     
#>  9         8 0.167     
#> 10         9 0.121     
#> 11        10 0.0708    
#> 12        11 0.0327    
#> 13        12 0.0117    
#> 14        13 0.00311   
#> 15        14 0.000584  
#> 16        15 0.0000687 
#> 17        16 0.00000381

# Simulation method
set.seed(1234)
tibble(survivors = map_dbl(1:5000000, ~max(16 - sum(runif(18) < 0.5), 0))) %>% 
  count(survivors) %>% 
  mutate(prob = n / sum(n))
#> # A tibble: 17 x 3
#>    survivors      n      prob
#>        <dbl>  <int>     <dbl>
#>  1         0   3160 0.000632 
#>  2         1  15545 0.00311  
#>  3         2  58477 0.0117   
#>  4         3 163756 0.0328   
#>  5         4 354090 0.0708   
#>  6         5 609096 0.122    
#>  7         6 833650 0.167    
#>  8         7 926929 0.185    
#>  9         8 835403 0.167    
#> 10         9 605720 0.121    
#> 11        10 354176 0.0708   
#> 12        11 163388 0.0327   
#> 13        12  57875 0.0116   
#> 14        13  15610 0.00312  
#> 15        14   2792 0.000558 
#> 16        15    314 0.0000628
#> 17        16     19 0.0000038

Created on 2021-10-20 by the reprex package (v2.0.0)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment