Conditional Block Sampling Parameters

[Colton Padilla]

Hi again everyone,

I am currently working on a model that I was trying to implement via Reversible Jump, but after some further thought, I have reparameterized the model differently. I believe this formulation is actually a more realistic representation of the system I am attempting to characterize. After this reparameterization, I have ran the model with all standard samplers and found a posterior that looks very similar to my many other attempts at this model (plot attached). Given this behavior that I have seen, I do think that the model that I have currently can accurately capture what I am looking to capture which is when the mean of a time-series trend stops changing. However, the traceplots and the posterior samples always seem to look like this where there is poor mixing and multimodal sticking points. 

After reexamining my problem more in depth, I think that the normal nimble samplers may not be optimal for this problem. My reason for thinking this is that when I attempt to sample the "stable" parameter, whatever 'lag' parameter is changed from one process to the other along with 'beta0' are not being updated to match the 'stable' parameter. Since those two parameters are not being updated, there are points in the posterior where the sampler has trouble jumping out of. Although my thought process may be wrong, I think that block sampling 'beta0' and 'stable' along with the 'lag' values that stable is proposing to shift from one process to another may provide better sampling of the posterior density than the default nimble samplers purely because of the way the model is structured. 

All that being said, I have taken a stab at creating my own sampler, but I am still seeing poor mixing in the stable parameter itself. I am figuring that there is likely a better way to sample these parameters but that I just don't quite know how to do it. If anyone has ran into a similar problem, or has any suggestions for how to change this sampler to improve the sampling of the posterior or even reparameterize the model, I would really appreciate it! 

 # Model prop
  code <- nimbleCode({
   
    #-----------#
    # Priors ####
    #-----------#
   
    # Intercept prior
    beta0 ~ dnorm(0, 0.1)
   
    # Observation error prior
    tau ~ dgamma(1, 1)
    sd <- pow(tau, -0.5)
   
    # breakpoint prior
    stable ~ dunif(1, maxlag)
    z[1:maxlag] <- rjProb(stable = stable, maxlag = maxlag)
   
    # Place the disjoint process prior
    sd_rw ~ dexp(1)
    tau_rw <- pow(sd_rw, -2)
    sd_process ~ dexp(1)
    tau_process <- pow(sd_process, -2)
    lag[1] ~ dnorm(0, 0.01)
    for(i in 2:maxlag){
      mean[i] <- lag[i - 1] * z[i] + beta0 * (1 - z[i])
      tau_step[i] <- tau_rw * z[i] + tau_process * (1 - z[i])
      lag[i] ~ dnorm(mean[i], tau_step[i])
    }
 
    #---------------#
    # Likelihood ####
    #---------------#
   
    # Loop over
    for(i in 1:maxlag){
      mu[start[i]:end[i]] <- lag[i]
      response[start[i]:end[i]] ~ dnorm_vec(mu[start[i]:end[i]],
                                            sd)
    }
   
  })