Game of Life glider, animated

Open Operational Research

Poking data, mostly Open

James Riley

3-Minute Read

Sometimes I see a Riddler Classic that makes me just want to go attack it with Monte Carlo analysis.

Let’s make a die and roll it 6 times, recording the outcome:

die <- 1:6

faces <- sample(seq_along(die), length(die), replace=T) # referring to the length of the die vector will help us if we do the extra credit

die <- die[faces]

die
## [1] 3 2 1 5 1 3

That’s 1 iteration. Loops will be necessary since I need the 2nd die to produce the 3rd die, so R’s clever vectorisation won’t help.

I’d like to track a ‘score’ - how far we are from the target of all faces equal. Let’s say the score is the number of distinct numbers still on the die, and the target is a score of 1.


update_die <- function(die){
  
  faces <- sample(seq_along(die), length(die), replace=T) # referring to the length of the die vector will help us if we do the extra credit

  die[faces]
}
update_score <- function(die){
  unique(die) %>% 
  length()
}

die <- 1:6
score <- update_score(die)

while(last(score) >1){
  die <- update_die(die)
  
  score <- c(score, update_score(die))
}

Anything we can do 1ce, we can do many times.

simulate_1_game <- function(index, size_of_die){
  die <- seq_len(size_of_die)
  score <- update_score(die)
  
  while(last(score) >1){
    die <- update_die(die)
    
    score <- c(score, update_score(die))
  }
  
  games <- seq_along(score)
  tibble(index=index, score=score, game=games)
}

monte_carlo_d6 <- map_dfr(seq_len(1e4), ~simulate_1_game(.x,6)) 
monte_carlo_d6 %>% 
    ggplot(aes(x=game, after_stat(density))) + geom_freqpoly(binwidth=1) + facet_wrap(~score, ncol=1, scales = "free_y") +
  theme + scale_y_continuous(labels=scales::percent)

monte_carlo_d6 %>% 
  tabyl(game, score) %>% 
  pivot_longer(-game, names_to = "score") %>% 

  ggplot(aes(x=game,y=value,fill=score)) + geom_area() + theme + ggthemes::scale_fill_few()

monte_carlo_d6 %>% 
  group_by(index) %>% 
  summarise(games_until_win = max(game)) %>% 
  ggplot(aes(x=games_until_win, after_stat(density)))  + geom_freqpoly(binwidth = 1) + 
    theme + scale_y_continuous(labels=scales::percent)

monte_carlo_d6 %>% 
  group_by(index) %>% 
  summarise(games_until_win = max(game)) %>% 
  tabyl(games_until_win) %>% 
  mutate(percent = scales::percent(percent)) %>% 
  knitr::kable()
games_until_win n percent
3 124 1.2%
4 499 5.0%
5 826 8.3%
6 979 9.8%
7 1018 10.2%
8 956 9.6%
9 823 8.2%
10 757 7.6%
11 660 6.6%
12 555 5.6%
13 439 4.4%
14 407 4.1%
15 330 3.3%
16 297 3.0%
17 229 2.3%
18 179 1.8%
19 137 1.4%
20 125 1.2%
21 100 1.0%
22 85 0.8%
23 77 0.8%
24 70 0.7%
25 56 0.6%
26 43 0.4%
27 41 0.4%
28 38 0.4%
29 23 0.2%
30 23 0.2%
31 17 0.2%
32 7 0.1%
33 13 0.1%
34 15 0.2%
35 15 0.2%
36 8 0.1%
37 2 0.0%
38 4 0.0%
39 3 0.0%
40 4 0.0%
41 2 0.0%
42 3 0.0%
43 2 0.0%
44 2 0.0%
45 1 0.0%
46 2 0.0%
47 2 0.0%
48 1 0.0%
50 1 0.0%

Other dice

We’d like to keep track of the size of our dice, tweak the simulate function. Just going to track average games until win to save on memory.

rm(monte_carlo_d6)


estimate_victory <- function(size_of_die){

  original_die <- seq_len(size_of_die)
  die <- original_die
  
  dice_pool <- sample(original_die, 1000*size_of_die, replace=T)
  
  play_1_game <- function(...){
    n=0L
    score <- 10 
   
    while(last(score) > 1){
      if(length(dice_pool) < size_of_die){
          dice_pool <<- sample(original_die, 100*size_of_die, replace=T)
      }
      
      n <- n+1L
      die <- die[dice_pool[original_die]] 
      dice_pool <<- dice_pool[-original_die] #Pop 1:n from the dice pool.
      score <- c(score, update_score(die))
    }   
    return(n)
  }
  
  map_int(1:1000, play_1_game) %>% 
    mean()
    
}

results <- future_map_dbl(1:100, estimate_victory)

Looks pretty much linear, and might end up more linear if I gave it more simulation time.

tibble(size_of_die = 1:100, result=results) %>% 
  ggplot(aes(x=size_of_die,y=result)) + theme + geom_line()

Recent Posts

categories