# Simple change detection model (2 means) being eerily slow

**URL:** <https://discourse.mc-stan.org/t/simple-change-detection-model-2-means-being-eerily-slow/24329>\
**Category:** Modeling\
**Tags:** fitting-issues, cognitive-science, brms\
**Created:** [September 9, 2021, 8:46pm UTC](https://discourse.mc-stan.org/t/simple-change-detection-model-2-means-being-eerily-slow/24329 "2021-09-09T20:46:08Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![fusaroli](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/fusaroli/32/477_2.png) [@fusaroli](https://discourse.mc-stan.org/u/fusaroli)\
**Post date:** [September 9, 2021, 8:46pm UTC](https://discourse.mc-stan.org/t/simple-change-detection-model-2-means-being-eerily-slow/24329/1 "2021-09-09T20:46:08Z")

</div>

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](https://github.com/fusaroli/CognitiveModelingRecoded/blob/main/13_ChangeDetection/13_data.txt)) 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](https://github.com/fusaroli/CognitiveModelingRecoded/tree/main/13_ChangeDetection)

```stan
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)

```

---

<div class="post-metadata">

**Author:** ![Ara\_Winter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ara_winter/32/4284_2.png) [@Ara\_Winter](https://discourse.mc-stan.org/u/Ara_Winter)\
**Post date:** [September 10, 2021, 11:07pm UTC](https://discourse.mc-stan.org/t/simple-change-detection-model-2-means-being-eerily-slow/24329/2 "2021-09-10T23:07:45Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![mike-lawrence](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/mike-lawrence/32/59_2.png) [@mike-lawrence](https://discourse.mc-stan.org/u/mike-lawrence)\
**Post date:** [September 10, 2021, 11:44pm UTC](https://discourse.mc-stan.org/t/simple-change-detection-model-2-means-being-eerily-slow/24329/3 "2021-09-10T23:44:44Z")

</div>

Be sure to see the [SUG chapter on change-point models](https://mc-stan.org/docs/2_27/stan-users-guide/change-point-section.html), where I suspect the marginalization trick permits a smoother posterior gradient for HMC to traverse.

---

<div class="post-metadata">

**Author:** ![fusaroli](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/fusaroli/32/477_2.png) [@fusaroli](https://discourse.mc-stan.org/u/fusaroli)\
**Post date:** [September 13, 2021, 7:41am UTC](https://discourse.mc-stan.org/t/simple-change-detection-model-2-means-being-eerily-slow/24329/4 "2021-09-13T07:41:09Z")

</div>

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

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)

```