Link to Bayesian Lilacs Lab
Simple linear regression
data {
int<lower=0> n; //
vector[n] y; // Y Vector
vector[n] x; // X Vector
}
parameters {
real beta0; // intercept
real beta1; // slope
real<lower=0> sigma; // standard deviation
}
model {
//Declare variable
vector[n] yhat;
//Priors
0.001, 0.001);
sigma ~ inv_gamma(0,1e2);
beta0 ~ normal(0,1e2);
beta1 ~ normal(
//model
yhat = beta0 + beta1 * (x - mean(x));
y ~ normal(yhat, sigma); }
Hockeystick (segmented regression) model
data {
int<lower=0> n; //
vector[n] y; // Y Vector
vector[n] x; // X Vector
}
parameters {
real beta0; // intercept
real beta1; // slope
real<lower=min(x), upper=max(x)> tau; // changepoint
real<lower=0> sigma; // variance
}
model {
vector[n] yhat;
// Priors
0.001, 0.001);
sigma ~ inv_gamma(112,30);
beta0 ~ normal(0, 2);
beta1 ~ normal(
tau ~ uniform(min(x), max(x));// Define the mean
for(i in 1:n)
yhat[i] = beta0 + int_step(x[i]-tau)*beta1*(x[i] - tau);
// likelihood / probability model
y ~ normal(yhat, sigma); }
Hierarchical hockey stick
data {
int<lower=0> n; // number of total observations
int<lower=1> nsites; // number of sites
int<lower=1,upper=nsites> site[n]; // vector of site IDs
real y[n]; // flowering date
real x[n]; // year
}
parameters {
real mu_beta0; // mean intercept
real<lower=0> sd_beta0; // s.d. of intercept
real mu_beta1; // mean slope
real<lower=0> sd_beta1; // s.d. of intercept
// overall changepoint and s.d.
real<lower = min(x), upper = max(x)> tau;
real<lower=0> sigma;
vector[nsites] beta0; //
vector[nsites] beta1; // unique intercept/slope site
}
model {
vector[n] yhat;
// Priors
0.001, 0.001);
sigma ~ inv_gamma(
130,20);
mu_beta0 ~ normal(0.001, 0.001);
sd_beta0 ~ inv_gamma(
5, .5);
mu_beta1 ~ normal(-.0.001, 0.001);
sd_beta1 ~ inv_gamma(
tau ~ uniform(min(x), max(x));
for (i in 1:nsites){
beta0[i] ~ normal(mu_beta0, sd_beta0);
beta1[i] ~ normal(mu_beta1, sd_beta1);
}
// The model!
for (i in 1:n)
yhat[i] = beta0[site[i]] + int_step(x[i]-tau) * beta1[site[i]] * (x[i] - tau);
// likelihood / probability model
y ~ normal(yhat, sigma); }
Hierarchical hockey stick with covariates
data {
int<lower=1> ncovars; // 3 = intercept + elevation + latitude
int<lower=0> n; // number of total observations
int<lower=1> nsites; // number of sites
int<lower=1,upper=nsites> site[n]; // vector of site ID's
real y[n]; // flowering date
real x[n]; // year
row_vector[ncovars] covars[n]; // matrix of: cbind(1, elevation, latitude)
}
parameters {
real mu_beta0[ncovars]; // regression coefficients of beta0
real<lower=0> sd_beta0; // s.d. of beta0
real mu_beta1; // mean of beta1
real<lower=0> sd_beta1; // s.d. of beta1
// overall change point and s.d. of process
real<lower=min(x), upper = max(x)> tau;
real<lower=0> sigma;
vector[ncovars] beta0[nsites];
vector[nsites] beta1;
}
model {
vector[n] yhat;
vector[n] x_beta0;
// Priors
8, 4);
sigma ~ normal(
// Establish (fairly well-informed) hierarchical priors
1] ~ normal(130,20);
mu_beta0[2] ~ normal(0,20); // per degree latitude
mu_beta0[3] ~ normal(0,20); // per 1km elevation
mu_beta0[4,4);
sd_beta0 ~ normal(
0, 10);
mu_beta1 ~ normal(0.2, 0.2);
sd_beta1 ~ normal(
tau ~ uniform(min(x), max(x));
for (i in 1:nsites){
beta0[i] ~ normal(mu_beta0, sd_beta0);
beta1[i] ~ normal(mu_beta1, sd_beta1);
}
// The model!
for (i in 1:n){
x_beta0[n] = covars[i] * beta0[site[i]];
yhat[i] = x_beta0[n] + int_step(x[i]-tau) * beta1[site[i]] * (x[i] - tau);
}
// likelihood / probability model
y ~ normal(yhat, sigma); }