luzhangstat/pathfinder

pathfinder's Issues

gpbf6.stan does not compile

In code Birthdays/model6_test.R file "/home/luzhang/Downloads/casestudies-master/Birthdays/gpbf6.stan", is read. Is it the same file as in working_directory/example/Birthdays/gpbf6.stan ? If so, while executing code:

## compile model and load data ##
model6 <- cmdstan_model(stan_file = "/Pathfinder/example/Birthdays/gpbf6.stan",
                        include_paths = "/Pathfinder/example/Birthdays")

Error appears:

Compiling Stan program...
Semantic error in '/tmp/RtmpoztcEC/model-12919d333b6.stan', line 24, column 26 to column 51:
-------------------------------------------------
22: // Basis functions for f1
23: real L_f1 = c_f1*max(xn);
24: matrix[N,M_f1] PHI_f1 = PHI_EQ(N, M_f1, L_f1, xn);
^
25: // Basis functions for f2
26: real period_year = 365.25/xsd;
-------------------------------------------------

A returning function was expected but an undeclared identifier 'PHI_EQ' was supplied.

make: *** [make/program:50: /tmp/RtmpoztcEC/model-12919d333b6.hpp] Error 1

Error: An error occured during compilation! See the message above for more information.

System information:

sessionInfo()
R version 4.1.0 (2021-05-18)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.3 LTS
Matrix products: default
BLAS/LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.8.so
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=C LC_PAPER=en_US.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] parallel stats graphics grDevices utils datasets methods base
other attached packages:
[1] optimx_2021-10.12 matrixStats_0.61.0 Matrix_1.3-3 patchwork_1.1.1 bayesplot_1.8.1 loo_2.4.1
[7] posterior_1.1.0 cmdstanr_0.4.0.9001 forcats_0.5.1 stringr_1.4.0 dplyr_1.0.7 purrr_0.3.4
[13] readr_2.1.0 tidyr_1.1.3 tibble_3.1.6 ggplot2_3.3.5 tidyverse_1.3.1 foreach_1.5.1
[19] rstan_2.26.4 StanHeaders_2.26.4

non existing variable

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?