Null included in results

[hyunjimoon]

Due to ADVI returning NULL, some results are NULL. NULL returns should be dropped and instead filled with the new attempt's nonNULL result.

> result$fits
[[1]]
NULL

[[2]]
    variable          mean      median           sd        mad            q5       q95
 lp__        -116164858.10 -2025040.00 1.120269e+09 2800980.55 -203134450.00 -51300.67
 lp_approx__         -7.76       -7.42 2.540000e+00       2.58        -12.28     -4.51
 a                    2.15        2.08 8.700000e-01       0.88          0.84      3.60
 b[1]                 0.99        0.91 8.400000e-01       0.90         -0.20      2.41
 b[2]                 1.71        1.76 9.200000e-01       0.93          0.18      3.13
 b[3]                 1.19        1.21 8.700000e-01       0.86         -0.16      2.64
 b[4]                -0.88       -0.89 8.700000e-01       0.92         -2.20      0.51
 b[5]                 0.26        0.33 8.400000e-01       0.75         -1.38      1.57
 b[6]                -1.49       -1.49 9.700000e-01       0.98         -3.04      0.13
 b[7]                -0.48       -0.53 8.500000e-01       0.86         -1.80      0.92

This is reproducible by:

set.seed(12)
generator <- function(clamp_val,clamp_dist, param, predictor = NULL){
  # fixed value across simultated datasets
  N <- clamp_val$N
  shape <- clamp_val$shape
  # fixed distribution across simultated datasets
  b <- clamp_dist$b
  # paramter - target, updated distribution 
  a_loc = mean(param$a)
  a_scale = sd(param$a)
  a <- rvar(rnorm(S, a_loc, a_scale))
  # predictor
  X = predictor$X
  # generate
  mu = exp(a + X %**% b)
  Y <- rfun(rgamma) (n = N, shape = shape, scale = mu/shape) 
  gen_rvars <- draws_rvars(N = N, Y = Y, K = K, X = X, 
                           a_loc = a_loc, a_scale = a_scale, shape = shape)
  # SDEdraws: rvar<S>[1] distributed to each simulation as rvar<1>[1]
  SBC_datasets(
    parameters = as_draws_matrix(param), 
    generated = draws_rvars_to_standata(gen_rvars)
  )
}
# dataset 
S = 10
N = 100
M = 300
K = 15
clamp_val <- list(N = N, shape = 1)
clamp_dist <- draws_rvars(b = rvar_rng(rnorm, ndraws = S, n = K, 0, 1)) #rgamma(S, shape = 0.01, scale = 0.01)
predictor <- draws_rvars(X = array(rnorm(S * N * K, mean = 1, sd = 1), dim = c(S, N, K)))
param_init_25 <- draws_rvars(a = rnorm(S, 2, 5)) # prior(normal(2, 2), class = "b", coef = "Intercept")
datasets_25 <- generator(
  clamp_val = clamp_val,
  clamp_dist = clamp_dist,
  param = param_init_25,
  predictor =  predictor
)
# backend 
mod_gr <- cmdstanr::cmdstan_model(stan_file = "gamma_reg")
backend_vi <- SBC_backend_cmdstan_variational(mod_gr, output_samples = M, algorithm = "fullrank")
result_25_vi <- compute_results(datasets_25, backend_vi, thin_ranks = 1)
>> result_25_vi[[3]]
# NULL

gamma_reg.stan file is here: https://github.com/hyunjimoon/SBC/blob/advi_self-calib/vignettes/models/gamma-reg.stan