cumulative_sum to reduce calculation

model.stan

data {
  int T;
  vector[T] Y;
}

parameters {
  real<lower=0, upper=2> mu_l;
  real<lower=0, upper=2> mu_r;
  real<lower=0> sigma;
}

transformed parameters {
  vector[T] lp;
  {
    vector[T] lp_l;
    vector[T] lp_r;
    for (cp in 1:T) {
      lp_l[cp] = normal_lpdf(Y[cp] | mu_l, sigma);
      lp_r[cp] = normal_lpdf(Y[cp] | mu_r, sigma);
    }
    lp_l = cumulative_sum(lp_l);
    lp_r = cumulative_sum(lp_r);
    lp = lp_l + (lp_r[T] - lp_r);
  }
}

model {
  target += log_sum_exp(lp);
}

model_slow.stan

data {
  int T;
  vector[T] Y;
}

parameters {
  real<lower=0, upper=2> mu_l;
  real<lower=0, upper=2> mu_r;
  real<lower=0> sigma;
}

transformed parameters {
  vector[T] lp;
  lp = rep_vector(0, T);
  for (cp in 1:T)
    for (t in 1:T)
      lp[cp] = lp[cp] + normal_lpdf(Y[t] | if_else(t <= cp, mu_l, mu_r), sigma);
}

model {
  target += log_sum_exp(lp);
}

run-model.R

library(rstan)

data <- list(T=length(Nile), Y=as.vector(Nile)/1000)

stanmodel <- stan_model(file='model.stan')
fit <- sampling(stanmodel, data=data, seed=1234)

ms <- extract(fit)
q <- colMeans(exp(ms$lp))
q <- q/sum(q)

# stanmodel_slow <- stan_model(file='model_slow.stan')
# fit_slow <- sampling(stanmodel_slow, data=data, seed=1234)