Create a gist now

Instantly share code, notes, and snippets.

@gmodena /changepoint.stan Secret
Last active Jul 31, 2016

What would you like to do?
// Built with stan 2.11
data {
int<lower=1> N;
real D[N];
}
// stan operates on log scale
transformed data {
real log_unif;
log_unif = log(N);
}
parameters {
real mu1;
real mu2;
real<lower=0> sigma1;
real<lower=0> sigma2;
}
// Marginalize out tau and
// calculate log_p(D | mu1, sd1, mu2, sd2)
// TODO: we can make this linear via dynamic programming
transformed parameters {
vector[N] log_p;
real mu;
real sigma;
log_p = rep_vector(log_unif, N);
for (tau in 1:N)
for (i in 1:N) {
mu = i < tau ? mu1 : mu2;
sigma = i < tau ? sigma1 : sigma2;
log_p[tau] = log_p[tau] + normal_lpdf(D[i] | mu, sigma);
}
}
model {
mu1 ~ normal(0, 100);
mu2 ~ normal(0, 100);
// scale parameters need to be > 0;
// we constrained sigma1, sigma2 to be positive
// so that stan interprets the following as half-normal priors
sigma1 ~ normal(0, 100);
sigma2 ~ normal(0, 100);
target += log_sum_exp(log_p);
}
//Draw the discrete parameter tau. This is highly inefficient
generated quantities {
int<lower=1,upper=N> tau;
tau = categorical_rng(softmax(log_p));
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment