non existing variable

[emilijapur]

Hello, I was trying to run code /example/Birthdays/model6_test.R

Running part:

if(fit_stan_flag){
  Stan_draws <- matrix(draws6, nrow = dim(draws6)[1])
  colnames(Stan_draws) <- dimnames(draws6)$variable
  unconstrained_draws <- unconstrain_cmd_draws(Stan_draws, posterior)
  
  lp_INV <- quantile(Stan_draws[, "lp__"], c(0.005, 0.995))
  
  save(file = "/home/rstudio/user_projects/Pathfinder/example/Birthdays/reference_model6.RData",
       list = c("unconstrained_draws", "lp_INV"))
}else{
  load("/home/rstudio/user_projects/Pathfinder/example/Birthdays/reference_model6.RData")
}

Shows error that file reference_model6.RData does not exist. And if I change fit_stant_flag=TRUE, then other error appears, that object draws6 does not exist:
Error in matrix(draws6, nrow = dim(draws6)[1]) : object 'draws6' not found

And if try to create object draws6 from the begining with code:

fit_stan_flag <- TRUE
if(fit_stan_flag){
  
  # use optimization to set initials #
  opt6 <- model6$optimize(data=standata6, init=0, algorithm='lbfgs',
                          history=100, tol_obj=10)
  #' Check whether parameters have reasonable values
  odraws6 <- opt6$draws()
  subset(odraws6, variable=c('intercept','sigma_','lengthscale_','sigma'), regex=TRUE)
  subset(odraws6, variable=c('beta_f3'))
  
  init6 <- sapply(c('intercept0','lengthscale_f1','lengthscale_f2',
                    'sigma_f1','sigma_f2','sigma_f4','sigma',
                    'beta_f1','beta_f2','beta_f3','beta_f4'),
                  function(variable) {as.numeric(subset(odraws6, variable=variable))})
  
  
  fit6 <- model6$sample(data=standata6,
                        chains=4, parallel_chains=4,
                        init=function() { init3 },
                        seed = 1948458383,
                        adapt_delta = 0.9,
                        thin = 100,
                        iter_warmup = 50000,
                        iter_sampling = 250000,
                        show_messages = TRUE,
                        sig_figs = 16)
  fit6$print("lp__")
  fit6$save_object(file = "/home/rstudio/user_projects/Pathfinder/example/Birthdays/birthday6_ref.RDS")
  p1 <- mcmc_trace(fit6$draws("lp__"), iter1 = 1) #c("lp__", "phi[1]", "lambda[1]", "theta1[1]")
  print(p1)
  
  #' Check whether parameters have reasonable values
  draws6 <- fit6$draws()
  summarise_draws(subset(draws6, variable=c('intercept','sigma_','lengthscale_',
                                            'sigma'), regex=TRUE))
  summarise_draws(subset(draws6, variable=c('beta_f3')))
  
  #' Compare the model to the data
  draws6 <- as_draws_matrix(draws6)
  Ef <- exp(apply(subset(draws6, variable='f'), 2, median))
  Ef1 <- apply(subset(draws6, variable='f1'), 2, median)
  Ef1 <- exp(Ef1 - mean(Ef1) + mean(log(data$births_relative100)))
  Ef2 <- apply(subset(draws6, variable='f2'), 2, median)
  Ef2 <- exp(Ef2 - mean(Ef2) + mean(log(data$births_relative100)))
  Ef_day_of_week <- apply(subset(draws6, variable='f_day_of_week'), 2, median)
  Ef_day_of_week <- exp(Ef_day_of_week - mean(Ef_day_of_week) + mean(log(data$births_relative100)))
  pf <- data %>%
    mutate(Ef = Ef) %>%
    ggplot(aes(x=date, y=births_relative100)) + geom_point(color=set1[2], alpha=0.2) +
    geom_line(aes(y=Ef), color=set1[1], alpha=0.75) +
    labs(x="Date", y="Relative number of births")
  pf1 <- data %>%
    mutate(Ef1 = Ef1) %>%
    ggplot(aes(x=date, y=births_relative100)) + geom_point(color=set1[2], alpha=0.2) +
    geom_line(aes(y=Ef1), color=set1[1]) +
    geom_hline(yintercept=100, color='gray') +
    labs(x="Date", y="Relative number of births")
  pf2 <- data %>%
    mutate(Ef2 = Ef2) %>%
    group_by(day_of_year2) %>%
    summarise(meanbirths=mean(births_relative100), meanEf2=mean(Ef2)) %>%
    ggplot(aes(x=as.Date("1987-12-31")+day_of_year2, y=meanbirths)) + geom_point(color=set1[2], alpha=0.2) +
    scale_x_date(date_breaks = "1 month", date_labels = "%b") +
    geom_line(aes(y=meanEf2), color=set1[1]) +
    geom_hline(yintercept=100, color='gray') +
    labs(x="Date", y="Relative number of births of year")
  pf3 <- ggplot(data=data, aes(x=day_of_week, y=births_relative100)) + geom_point(color=set1[2], alpha=0.2) +
    scale_x_continuous(breaks = 1:7, labels=c('Mon','Tue','Wed','Thu','Fri','Sat','Sun')) +
    geom_line(data=data.frame(x=1:7,y=Ef_day_of_week), aes(x=x, y=Ef_day_of_week), color=set1[1]) +
    geom_hline(yintercept=100, color='gray') +
    labs(x="Date", y="Relative number of births of week")
  (pf + pf1) / (pf2 + pf3)
} 

then I get error that object init3 does not exist:

Error in init() : object 'init3' not found

How shall I get init3 object?

[LuZhangstat]

Hello. Thank you for reporting the error. If fit_stan_flag = true then we fit the model by Stan to generate the reference posterior samples for the birthday example. This process takes several days. So I ran that part on my old desktop, and I edited the code with my laptop. I think that's why there is a typo in the code. The init3 should be init6. I have fixed that bug in the code. Please let me know if it still doesn't work.