Error in example

[mkamenet3]

Hello there!

I was working through your example code in example(sampler) and when I run out <- sampler(simu_dat,model_options0,mcmc_options0) for the first example, I get the following error: Error in log_p[t + 1] <- log_v[t + 1] - log_v[t] + log(b) + log_marginal(dat[i, : replacement has length zero. I was wondering if you have any insights on how to resolve this? Below is all the code I've run and sessionInfo

library(rewind)
example(sampler)

library(matrixStats)
library(ars) #adaptive rejection sampling


# color palette for heatmaps:
YlGnBu5   <- c('#ffffd9','#c7e9b4','#41b6c4','#225ea8','#081d58',"#092d94")
hmcols    <- colorRampPalette(YlGnBu5)(256)

# simulate data:
L0 <- 100
options_sim0  <- list(N = 200,  # sample size.
                      M = 3,   
                      L = L0,   
                      K = 8,    
                      theta = rep(0.9,L0), # true positive rates.
                      psi   = rep(0.1,L0), # false positive rates.
                      alpha1 = 1, # half of the people have the first machine.
                      frac = 0.2, # fraction of positive dimensions (L-2M) in Q.
                      #pop_frac = rep(1/K0,K0) # population prevalences.
                      pop_frac = c(rep(2,4),rep(1,4)) # population prevalences.
)
str(options_sim0)


image(simulate_data(options_sim0,SETSEED = TRUE)$datmat,col=hmcols)
simu     <- simulate_data(options_sim0, SETSEED=TRUE)
str(simu)
simu_dat <- simu$datmat
str(simu_dat)

rle(simu$Z) 

#
# visualize truth:
#
plot_truth(simu,options_sim0)

#
# specifying options:
#

# model options:
m_max0 <- 5
model_options0 <- list(
  n   = nrow(simu_dat),
  t_max  = 40,
  m_max  = m_max0,
  b  = 1, # Dirichlet hyperparameter; in the functions above,
  # we used "b" - also can be called "gamma"!.
  #Q  = simu$Q,
  a_theta = replicate(L0, c(9,1)),
  a_psi   = replicate(L0, c(1,9)),
  #theta = options_sim0$theta,
  #psi   = options_sim0$psi,
  alpha   = 0.1,
  frac = 0.2,
  #p_both      = rep(0.5,3),#,c(0.5,0.5^2,0.5^3,0.5^4,0.5^5)
  #p0 = rep(0.5,m_max0), # <--- this seems to make a difference in convergence.
  log_pk = "function(k) {log(0.1) + (k-1)*log(0.9)}"# Geometric(0.1).
  #Prior for the number of components.
)

# pre-compute the log of coefficients in MFM:
model_options0$log_v<-mfm_coefficients(eval(parse(text=model_options0$log_pk)),#response probability/Bernoulli
                                       model_options0$b, #Dirichelet hyperparameter
                                       model_options0$n,#200, number of subjects
                                       model_options0$t_max+1) #upper limit to compute the coefficient
# mcmc options:
mcmc_options0 <- list(
  n_total = 200,
  n_keep  = 50,
  n_split = 5, #??? number of splits
  print_mod = 10,
  constrained = TRUE, # <-- need to write a manual about when these options are okay.
  block_update_H = TRUE,
  block_update_Q = !TRUE,
  ALL_IN_ONE_START =!TRUE,  # <--- TRUE for putting all subjects in one cluster,
  # FALSE by starting from a hierechical clustering
  # (complete linkage) and cut
  # to produce floor(t_max/4). Consider this as a warm start.
  MORE_SPLIT = TRUE,
  print_Q    = TRUE,
  init_frac  = c(0.51,0.99)
)

# run posterior algorithm for simulated data:
out <- sampler(simu_dat,model_options0,mcmc_options0)

 sessionInfo()
R version 4.2.1 (2022-06-23)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS Monterey 12.6.1

Matrix products: default
LAPACK: /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] ars_0.6            matrixStats_0.62.0 rewind_0.7.0      

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.8.3      compiler_4.2.1    later_1.3.0       urlchecker_1.0.1  prettyunits_1.1.1 profvis_0.3.7    
 [7] remotes_2.4.2     tools_4.2.1       digest_0.6.29     pkgbuild_1.3.1    pkgload_1.3.0     evaluate_0.15    
[13] memoise_2.0.1     lifecycle_1.0.1   rlang_1.0.5       shiny_1.7.1       cli_3.3.0         rstudioapi_0.13  
[19] curl_4.3.2        yaml_2.3.5        xfun_0.31         fastmap_1.1.0     withr_2.5.0       stringr_1.4.0    
[25] knitr_1.39        fs_1.5.2          htmlwidgets_1.5.4 devtools_2.4.4    rprojroot_2.0.3   glue_1.6.2       
[31] R6_2.5.1          processx_3.6.1    rmarkdown_2.14    sessioninfo_1.2.2 purrr_0.3.4       callr_3.7.0      
[37] magrittr_2.0.3    usethis_2.1.6     promises_1.2.0.1  ps_1.7.1          ellipsis_0.3.2    htmltools_0.5.2  
[43] mime_0.12         xtable_1.8-4      httpuv_1.6.5      stringi_1.7.6     miniUI_0.1.1.1    cachem_1.0.6     
[49] crayon_1.5.1