Skip to content

Instantly share code, notes, and snippets.

@gmodena
Last active July 31, 2016 19:58

Revisions

  1. gmodena renamed this gist Jul 31, 2016. 1 changed file with 0 additions and 0 deletions.
    File renamed without changes.
  2. gmodena created this gist Jul 31, 2016.
    55 changes: 55 additions & 0 deletions change point.stan
    Original file line number Diff line number Diff line change
    @@ -0,0 +1,55 @@
    // 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));
    }