The @slic macro

Features

• Stan blocks, types, dimensions and constraints are inferred for each parameter from the combination of model and data.
• Models can define default priors, which are not used when the corresponding parameter is provided.
• Models can return expressions, with the return value assigned to the lhs of a sampling statement involving that model.

Missing features

• Prettyfication of generated code
• Control flow (loops, if-else-blocks)
• User defined functions
• Broadcasting
• More

Case studies

The below case studies only use dummy data to generate the stan models. Some case studies are not functional yet, as in they use currently unimplemented features.

Multilevel regression modeling (Radon data)

Qualitatively reproduces Mitzi Morris’ radon case study.

Complete pooling (radon_cp below)

SlicModel

Specification

begin
    alpha ~ normal(0, 10)
    beta ~ normal(0, 10)
    sigma ~ normal(0, 10; lower = 0.0)
    y ~ normal(alpha + beta * x, sigma)
end

Stan

functions {
vector normal_lpdfs(vector obs, vector loc, real scale){
    int n = dims(loc)[1];
        vector[n] rv;
    for(i in 1:n){
        rv[i] = normal_lpdf(obs[i] | loc[i], scale);
    }
    return rv;

}

}

data {
int y_n;
vector[y_n] y;
int x_n;
vector[x_n] x;
}

parameters {
real alpha;
real beta;
real<lower=0.0> sigma;
}

model {
alpha ~ normal(0, 10);
beta ~ normal(0, 10);
sigma ~ normal(0, 10);
y ~ normal((alpha + (beta * x)), sigma);
}

generated quantities {
vector[x_n] y_likelihood = normal_lpdfs(y, (alpha + (beta * x)), sigma);
vector[y_n] y_gen = to_vector(normal_rng((alpha + (beta * x)), sigma));
}

No pooling

SlicModel

Specification

begin
    n_counties = max(county)
    alpha ~ normal(0, 10; n = n_counties)
    beta ~ normal(0, 10)
    sigma ~ normal(0, 10; lower = 0.0)
    y ~ normal(alpha[county] + beta * x, sigma)
end

Stan

functions {
vector normal_lpdfs(vector obs, vector loc, real scale){
    int n = dims(loc)[1];
        vector[n] rv;
    for(i in 1:n){
        rv[i] = normal_lpdf(obs[i] | loc[i], scale);
    }
    return rv;

}

}

data {
int county_n;
array [county_n] int county;
int y_n;
vector[y_n] y;
int x_n;
vector[x_n] x;
}

transformed data {
int n_counties = max(county);
}

parameters {
vector[n_counties] alpha;
real beta;
real<lower=0.0> sigma;
}

model {
alpha ~ normal(0, 10);
beta ~ normal(0, 10);
sigma ~ normal(0, 10);
y ~ normal((alpha[county] + (beta * x)), sigma);
}

generated quantities {
vector[county_n] y_likelihood = normal_lpdfs(y, (alpha[county] + (beta * x)), sigma);
vector[y_n] y_gen = to_vector(normal_rng((alpha[county] + (beta * x)), sigma));
}

Partial pooling

SlicModel

Specification

begin
    n_counties = max(county)
    mu_alpha ~ normal(0, 10)
    sigma_alpha ~ normal(0, 10; lower = 0)
    alpha ~ normal(mu_alpha, sigma_alpha; n = n_counties)
    beta ~ normal(0, 10)
    sigma ~ normal(0, 10; lower = 0.0)
    y ~ normal(alpha[county] + beta * x, sigma)
end

Stan

functions {
vector normal_lpdfs(vector obs, vector loc, real scale){
    int n = dims(loc)[1];
        vector[n] rv;
    for(i in 1:n){
        rv[i] = normal_lpdf(obs[i] | loc[i], scale);
    }
    return rv;

}

}

data {
int county_n;
array [county_n] int county;
int y_n;
vector[y_n] y;
int x_n;
vector[x_n] x;
}

transformed data {
int n_counties = max(county);
}

parameters {
real mu_alpha;
real<lower=0> sigma_alpha;
vector[n_counties] alpha;
real beta;
real<lower=0.0> sigma;
}

model {
mu_alpha ~ normal(0, 10);
sigma_alpha ~ normal(0, 10);
alpha ~ normal(mu_alpha, sigma_alpha);
beta ~ normal(0, 10);
sigma ~ normal(0, 10);
y ~ normal((alpha[county] + (beta * x)), sigma);
}

generated quantities {
vector[county_n] y_likelihood = normal_lpdfs(y, (alpha[county] + (beta * x)), sigma);
vector[y_n] y_gen = to_vector(normal_rng((alpha[county] + (beta * x)), sigma));
}

Cross validation

SlicModel

Specification

begin
    n_counties = max(county)
    mu_alpha ~ normal(0, 10)
    sigma_alpha ~ normal(0, 10; lower = 0)
    alpha ~ normal(mu_alpha, sigma_alpha; n = n_counties)
    beta ~ normal(0, 10)
    sigma ~ normal(0, 10; lower = 0.0)
    y ~ normal(alpha[county] + beta * x, sigma)
end

Stan

functions {
vector normal_lpdfs(vector obs, vector loc, real scale){
    int n = dims(loc)[1];
        vector[n] rv;
    for(i in 1:n){
        rv[i] = normal_lpdf(obs[i] | loc[i], scale);
    }
    return rv;

}

}

data {
int county_n;
array [county_n] int county;
int y_n;
vector[y_n] y;
int x_n;
vector[x_n] x;
}

transformed data {
int n_counties = max(county);
}

parameters {
real mu_alpha;
real<lower=0> sigma_alpha;
real beta;
real<lower=0.0> sigma;
}

model {
mu_alpha ~ normal(0, 10);
sigma_alpha ~ normal(0, 10);
beta ~ normal(0, 10);
sigma ~ normal(0, 10);
}

generated quantities {
vector[n_counties] alpha = to_vector(normal_rng(rep_vector(mu_alpha, n_counties), sigma_alpha));
vector[county_n] y_likelihood = normal_lpdfs(y, (alpha[county] + (beta * x)), sigma);
vector[y_n] y_gen = to_vector(normal_rng((alpha[county] + (beta * x)), sigma));
}