compare_simulations.R

library(tidyverse)
library(rethinking)

# khench implementation -------------------
# updating funtion for an individual year
progress_year <- function(data, year, max_age = 65, n_births = 20, aom = 18){
  # initialize new cohort (new row for the population table)
  new_cohort <- tibble(
    age = 1,
    married = as.integer(0),
    happiness = seq(from = -2, to = 2, length.out = n_births),
    year_of_birth = year)
  
  # update the population table
  data %>% 
    # age the population
    mutate(age = age + 1) %>%
    # introduce new cohort into the population
    bind_rows(., new_cohort) %>%
    # determine this years marriages
    mutate(married = if_else(age >= aom & married == 0,
                             rbern(n(), inv_logit(happiness - 4)),
                             married )) %>%
    # move people "to spain"
    filter(age <= max_age)
}

sim_tidy <- function(seed = 1977, n_years = 1000, max_age = 65, n_births = 20, aom = 18){
  # initialize the population
  empty_tibble <- tibble(age = double(),
                         married = integer(),
                         happiness = double())
  
  # simulate the progressive years
  1:n_years %>% 
    reduce(.f = progress_year,
           .init = empty_tibble,
           max_age = max_age, n_births = n_births, aom = aom)
}


# rpruim implementation ----------------------------------
new_borns <- function(n = 20) {
  tibble(
    A = 1,   # 1 year old
    M = 0,   # not married
    H = seq(-2, 2, length.out = n)   # range of happiness scores
  )
}

# Everyone ages;
# Some unmarried adults get married (happier people are more likely to)
# Old people die;
# New people are born.

update_population <- function(
  pop, N_births = 20, aom = 18, max_age = 65) {
  
  pop %>%
    mutate(
      A = A + 1, # everyone gets one year older
      # some people get married
      M = ifelse(M >= 1, 1, (A >= aom) * rbinom(nrow(pop), 1, inv_logit(H - 4)))
    ) %>%
    filter(A <= max_age) %>%    # old people die
    bind_rows(
      new_borns(N_births)            # new people are born
    )
}

# wrapping  aproaches into simgle function for comparison
k_fun <- function(){
  set.seed(42)
  Pop <- sim_tidy(n_years = 65, n_births = 21)
  Pop
}

r_fun <- function(){
  set.seed(1997)
  Pop <- new_borns(20)          # year 1
  for(i in 2:65) {            # years 2 through 65
    Pop <- update_population(Pop, aom = 18, max_age = 65, N_births = 21)
  } 
  Pop
}

# running the comparison
tictoc::tic(); k_fun(); tictoc::toc()
#> 0.555 sec elapsed
tictoc::tic(); r_fun(); tictoc::toc()
#> 0.382 sec elapsed