Simple change detection model (2 means) being eerily slow

Modeling
, ,
1

I am reimplementing the change detection model from Lee&Wagenmakers Bayesian Cognitive Modeling (based on Smira’s earlier recoding).
The data (CognitiveModelingRecoded/13_data.txt at main · fusaroli/CognitiveModelingRecoded · GitHub) is generated by having a certain mean til time t and a different mean after. We want to infer t.

I am however only getting very slow sampling, both when coding it directly in Stan, and when using an equivalent brms model, and I am not sure whether I am missing more effective ways of coding the model.

Reprex is here: CognitiveModelingRecoded/13_ChangeDetection at main · fusaroli/CognitiveModelingRecoded · GitHub

data { 
  int n;
  vector[n] t;
  vector[n] c;
}

parameters {
  vector[2] mu;
  real<lower=0> sigma;
  real<lower=0,upper=n> tau;
  vector[2] muprior;
  real<lower=0> sigmaprior;
} 

model { 
  // Group Means
  mu ~ normal(0, 10);
  muprior ~ normal(0, 10);
  // Standard deviation
  sigma ~ normal(0, 10);
  sigmaprior ~ normal(0, 10);
    
  // Which Side is Time of Change Point?
  // Data Come From A Gaussian
  for (i in 1:n) {
    if ((t[i] - tau) < 0.0)
      c[i] ~ normal(mu[1], sigma);
    else 
      c[i] ~ normal(mu[2], sigma);
  }
}

Brms recoding

bform <- bf(
  c ~ Intercept1 * step(change - t) +  # Section 1
    Intercept2 * step(t - change),  # Section 2,
  Intercept1 + Intercept2 ~ 1,  # Fixed intercept and slopes
  change ~ 1,  # Per-person changepoints around pop mean
  nl = TRUE
)

# Priors
bprior <- prior(normal(35, 10), nlpar = "Intercept1") +
  prior(normal(35, 10), nlpar = "Intercept2") +
  prior(uniform(0, 1178), nlpar = "change")  # Within observed range



# Fit it!

m <- brm(
  bform,
  df,
  family = gaussian,
  prior = bprior,
  sample_prior = T,
  seed = 123,
  chains = 1,
  cores = 1,
  iter = 500,
  refresh=10,
  backend = "cmdstanr",
  threads = threading(2),
  control = list(
    max_treedepth = 20,
    adapt_delta = 0.99)
)

stancode(m)
1 Like
2

Afternoon! A few things might help to track down the problem. Can you post your Stan model call? And a snippet of the data or better some simulated data with know parameters?

A few things right off:

How much data is going into the model and what do those plots look like?

The tree depth and adapt_delta in the brms call look really high to me. Is the model throwing errors with the default settings (12 and 0.85 I think).

I usually run 6 cores and 4 chains on my laptop unless I am debugging.

Are the priors reasonable given prior knowledge and the data? The uniform prior stands out as a bad idea and might slow things down.

2 Likes
3

Be sure to see the SUG chapter on change-point models, where I suspect the marginalization trick permits a smoother posterior gradient for HMC to traverse.

3 Likes
4

Sorry, I had attached the reprex as a link. Now I played with the example a bit more and I put all the needed code here for both cmdstanr and brms.

  1. The data is generated to be very clear and neither too little nor too much

    `c <- c(rnorm(200, 30, 5), rnorm(200, 50, 5))`

  2. Priors and inits are set to be in the area of the generated data (cheating, of course, but I wanted to start with an ideal case)

  3. tree depth and adapt delta are high because I was getting divergences even with high iterations.

  4. I’ll check the marginalization trick next :-)

pacman::p_load(
  cmdstanr,
  posterior,
  here,
  job,
  brms
)


c <- c(rnorm(200, 30, 5), rnorm(200, 50, 5))

n <- length(c)
t <- 1:n

data <- list(c=c, n=n, t=t) # to be passed on to Stan

myinits <- list(
  list(mu=c(40, 40), sigma = 5, tau = n / 2))

stan_file <- write_stan_file("
data { 
  int n;
  vector[n] t;
  vector[n] c;
}

// The parameters accepted by the model. Our model
// accepts two parameters 'mu' and 'sigma'.
parameters {
  vector[2] mu;
  real<lower=0> sigma;
  real<lower=0,upper=n> tau;
  vector[2] muprior;
  real<lower=0> sigmaprior;
} 

model { 
  // Group Means
  mu ~ normal(40, 10);
  muprior ~ normal(40, 10);
  // Standard deviation
  sigma ~ normal(0, 10);
  sigmaprior ~ normal(0, 10);
  tau ~ uniform(0, 400);
    
  // Which Side is Time of Change Point?
  // Data Come From A Gaussian
  for (i in 1:n) {
    if ((t[i] - tau) < 0.0)
      c[i] ~ normal(mu[1], sigma);
    else 
      c[i] ~ normal(mu[2], sigma);
  }
}
")

mod <- cmdstan_model(stan_file, cpp_options = list(stan_threads = TRUE), pedantic = TRUE)


job::job({
  samples <- mod$sample(
    data = data,
    seed = 123,
    chains = 1,
    parallel_chains = 1,
    threads_per_chain = 1,
    iter_warmup = 300,
    iter_sampling = 300,
    refresh = 10,
    init = myinits,
    max_treedepth = 20,
    adapt_delta = 0.8,
  )
})

df <- tibble(c=c, t=t)

bform <- bf(
  c ~ Intercept1 * step(change - t) +  # Section 1
    Intercept2 * step(t - change),  # Section 2,
  Intercept1 + Intercept2 ~ 1,  # Fixed intercept and slopes
  change ~ 1,  # Per-person changepoints around pop mean
  nl = TRUE
)

# Priors
bprior <- prior(normal(40, 10), class = "b", coef="Intercept", nlpar = "Intercept1") +
  prior(normal(40, 10), class = "b", coef="Intercept", nlpar = "Intercept2") +
  prior(normal(0, 10), class="sigma") +
  prior(uniform(0, 1178), class = "b", coef="Intercept", nlpar = "change") # Within observed range



# Fit it!
job::job({
m <- brm(
  bform,
  df,
  family = gaussian,
  prior = bprior,
  sample_prior = T,
  seed = 123,
  chains = 1,
  cores = 1,
  iter = 500,
  refresh = 10,
  backend = "cmdstanr",
  threads = threading(2),
  control = list(
    max_treedepth = 20,
    adapt_delta = 0.8)
)})

stancode(m)
3 Likes