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

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.

``````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",
control = list(
max_treedepth = 20,
)

stancode(m)
``````
1 Like

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

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

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,
iter_warmup = 300,
iter_sampling = 300,
refresh = 10,
init = myinits,
max_treedepth = 20,
)
})

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",