Different results between rstan and BayesianTools SMC for a non-latent-state model

[Jefflier]

I am currently comparing results obtained with rstan and BayesianTools for a model that does not involve latent states, and I observe systematic differences between the two when using SMC-based inference. The model is a static parameter model with Gaussian priors (including a multivariate normal block) and a beta-binomial observation likelihood. In the rstan implementation, inference is performed using NUTS, while in BayesianTools the same log-likelihood and log-prior are supplied via createBayesianSetup() and sampled using runSMC(). Despite matching the parameterization, prior specifications, and likelihood structure, posterior summaries from BayesianTools SMC differ noticeably from those obtained with rstan, and the discrepancy persists even when the number of particles is increased. I would like to ask whether such differences are expected for static models, or whether there are default tempering, scaling, or other internal SMC mechanisms in BayesianTools that could lead to a posterior target differing from that of a standard MCMC implementation.

Fig. 1_dengue.tif

fig3.pdf

[florianhartig]

Hi Luojunteng,

I guess I need more info to say anything about this. Are these two figures from STAN / BT? If you have a non-latent model, the states aren't that informative - how do the parameter estimates differ?

Ideally, you would provide some example code.

[Jefflier]

Thank you for your reply!

Yes, the two figures are from rstan and BayesianTools, respectively. I implement a non-latent model for compare two software because RStan does not support latent model. Several key parameters show noticeably broader and shifted marginal posterior distributions under BayesianTools SMC compared to the Stan results, while Stan diagnostics indicate good convergence (Rhat ≤ 1.01).

I compare the parameter estimates results in two software, and I find the results are not the same.

Here are the codes:.

For RStan

model {
  // Priors
  // theta_sm ~ normal(0, 5);
  // theta_st ~ normal(0, 5);
  // theta_0_init ~ normal(0, 2);
  // theta_0_0 ~ normal(0, 2);
  // phi ~ exponential(0.01);
  vector[18] theta_joint;
  theta_joint = append_row(theta_sm, theta_st);
  // 
  target += 0.1 * normal_lpdf(theta_0_init | 0, 2);
  target += 0.1 * normal_lpdf(theta_0_0 | 0, 2);
  target += 0.1 * multi_normal_lpdf(theta_joint | mu_theta, sigma_theta);
  // likelihood
  for (t in (n_serial + 1):N) {
    
      if (I_new[t] > 0) {
        real alpha = 1 + I_new[t];
        real beta_param = 1 + pop - I_new[t];
        target += 0.001 * beta_binomial_lpmf(cases[t] | pop, alpha, beta_param);
      }
    
    // print(cases[t]);
  }
  
}

For SMC

ll = function(params, browse = F){
  if(browse) browser()
  b0.spline = generate_spline(y.in = params[19:length(params)], bs.in = b0_basis)
  mosq.biv = generate_biv_spline(coef.prop = params[1:9], dat.range = mosq.range)
  temp.biv = generate_biv_spline(coef.prop = params[10:18], dat.range = temp.range)
  
  beta.ts = c(rep(0, length(weights.serialinterval)), exp(sapply((length(weights.serialinterval) + 1):nrow(GZ.data), function(tt)
    calc_prop_beta(tt, biv.temp = temp.biv, biv.mosq = mosq.biv, browse = F, b0.spline.tmp = b0.spline))))
  
  I.new.sim = sim_epidemic(n.reps = 100, beta.vec = beta.ts)
  
  LL = lik_weather(sim = I.new.sim)
  if(is.na(LL)){
    LL = -(.Machine$double.xmax)/10
  }
  return(0.001 * (LL))
}


p = function(params){
  p1 = dmvnorm(params[1:18], param.mean, makePositiveDefinite(sigma.mean), log = TRUE)
  p2 = dnorm(mean(params[19:51]), 0, 2, log = TRUE)
  return(0.001 * (sum(p1) + sum(p2)))
}

s = function(n1=1, n2=33){
  s1 = rmvnorm(n1, param.mean, makePositiveDefinite(sigma.mean))
  s2 = rnorm(n2, 0, 0.5)
  return(c(s1, s2))
}

I have changed sim_epidemic to

sim_epidemic = function(dat = GZ.data, n.reps = 1, beta.vec, browse = T){
  if(browse) browser()
  
  I.new = matrix(0, nrow = n.reps, ncol = nrow(dat))
  I.old = matrix(0, nrow = n.reps, ncol = 1)
  imp.cases = matrix(rep(GZ.data$imported.cases, n.reps), nrow = n.reps, ncol = nrow(dat), byrow = T)
  SI = matrix(weights.serialinterval)
  
  for(tt in (length(weights.serialinterval) + 1):nrow(dat)){
    A = I.new[, seq(tt-1, tt-length(weights.serialinterval), -1)]
    A.SI = A %*% SI
    
    Imp = imp.cases[, seq(tt-1, tt-length(weights.serialinterval), -1)]
    Imp.SI = Imp %*% SI
    
    I.old = A.SI + Imp.SI
    
    ind.sim = which(I.old > 0)
    if(length(ind.sim) != 0){
      I.new[ind.sim, tt] = I.old[ind.sim] * beta.vec[tt]

    }

  }
  return(I.new)
}

and the original code is in https://github.com/roidtman/NatComm_dengue_China/blob/master/main_text_code/3_smc_main_text.R