StanBlocks.jl

Appearance

Golf putting ​

This guided example follows the model sequence in Andrew Gelman's golf putting case study. The question is not merely whether longer putts are harder, but which physical mechanisms make them harder. We therefore move from a descriptive regression to models of aiming error and distance control.

The small arrays below are build fixtures, not the original golf data. Each comparison evaluates the displayed Julia source while the documentation is built and shows the complete generated Stan program beside it. The comparison can also be opened as the same side-by-side modal used by the feature atlas; readers do not have to execute anything to see the Stan code.

Logistic regression ​

The baseline treats the success probability as a logistic curve of distance,

It is useful as a descriptive benchmark, but a and b do not say why a putt misses. In the SLIC source, the first use of the fresh names a and b on the left of ~ introduces scalar parameters. The call logistic(; y, n, x) binds the remaining free names as data. StanBlocks infers the scalar and vector shapes from those values and emits the data, parameter, and model declarations visible in the Stan pane.

julia
using StanBlocks
x = distance_tolerance = overshot = randn(10)
y = fill(1, 10)
n = fill(2, 10)
R = r = 1.

logistic = @slic begin 
    a ~ flat()
    b ~ flat() 
    y ~ binomial_logit(n, a + b * x)
end
logistic_posterior = logistic(;y,n,x)
stan
functions {
vector binomial_logit_lpmfs(
    array[] int y,
    array[] int args1,
    vector args2
) {
    return jbroadcasted_binomial_logit_lpmfs(y, args1, args2);
}
vector jbroadcasted_binomial_logit_lpmfs(
    array[] int x1,
    array[] int x2,
    vector x3
) {
    int n = dims(x1)[1];
    vector[n] rv;
    for(i in 1:n) {
        rv[i] = binomial_logit_lpmfs(
            broadcasted_getindex(x1, i),
            broadcasted_getindex(x2, i),
            broadcasted_getindex(x3, i)
        );
    }
    return rv;
}
real binomial_logit_lpmfs(
    int args1,
    int args2,
    real args3
) {
    return binomial_logit_lpmf(args1 | args2, args3);
}
int broadcasted_getindex(array[] int x, int i) {
    return x[i];
}
real broadcasted_getindex(vector x, int i) {
    return x[i];
}
array[] int binomial_logit_int_rng(
    int anontok__1,
    array[] int N,
    vector eta
) {
    int n = anontok__1;
    if (dims(N)[1] != n) reject("binomial_logit_rng: dim mismatch — `N` dim 1 (= ", dims(N)[1], ") does not match `n` (= ", n, "), inferred from `anontok__1` dim 1. `n` sizes: `anontok__1` dim 1 (= ", anontok__1, "), `N` dim 1 (= ", dims(N)[1], ").");
    return binomial_rng(N, inv_logit(eta));
}
}
data {
    int y_n;
    array[y_n] int y;
    int n_n;
    array[n_n] int n;
    int x_n;
    vector[x_n] x;
}
transformed data {
}
parameters {
    real a;
    real b;
}
transformed parameters {
}
model {
    y ~ binomial_logit(n, (a + (b * x)));
}
generated quantities {
    vector[y_n] y_likelihood = binomial_logit_lpmfs(y, n, (a + (b * x)));
    array[n_n] int y_gen = binomial_logit_int_rng(y_n, n, (a + (b * x)));
}

Modelling based on first principles ​

A successful putt must satisfy two conditions: it must start within the angular window subtended by the hole, and its speed must leave the ball within a tolerable stopping-distance range. Splitting those mechanisms into SLIC submodels makes the assumptions explicit and lets later models reuse them.

Submodels ​

Both definitions below are anonymous @slic values. Their free names become data when supplied as keyword arguments, their local sampled variables remain local to each use, and their final return value is the quantity embedded by the parent model. A parent writes p ~ submodel(; inputs...); this is SLIC composition, not a probability distribution named submodel. The compiler inlines the submodel hygienically and carries its parameters and transformed values into the generated Stan blocks.

Angle submodel (angle_submodel below) ​

For a hole of radius R, a ball of radius r, and putt distance x, the allowable aiming angle is asin((R-r)/x). The positive parameter sigma represents angular error. Assuming centered Gaussian aiming error, the returned vector is the probability of landing inside the angular window. The elementwise ./ and the vector-valued Phi call are enough for StanBlocks to infer a vector result; no Stan declaration is written by hand.

julia
angle_submodel = @slic begin 
    threshold_angle = asin((R - r) ./ x) 
    sigma ~ flat(;lower=0.)
    sigma_degrees = sigma * 180 / pi
    return 2 * Phi(threshold_angle / sigma) - 1
end

Distance submodel (distance_submodel below) ​

The second component models speed control. overshot locates the target stopping point beyond the hole, while distance_tolerance defines the accepted interval around it. sigma_distance is constrained positive by the lower=0. distribution keyword. The difference of two normal CDF values is the probability that the stopping distance lies inside that interval.

julia
distance_submodel = @slic begin 
    sigma_distance ~ std_normal(;lower=0.)
    return Phi(
        (distance_tolerance - overshot) ./ ((x + overshot) * sigma_distance)
    ) - Phi(
        (-overshot)./ ((x + overshot) * sigma_distance)
    )
end

Angle model ​

This is the first mechanistic replacement for logistic regression. The parent model obtains the vector p by embedding angle_submodel, then uses that probability in the same binomial observation model. Compare the generated Stan with the baseline: the regression coefficient disappears, while the angular error parameter and its deterministic probability calculation appear.

julia
angle = @slic begin 
    p ~ angle_submodel(;R,r,x)
    y ~ binomial(n, p)
end
angle_posterior = angle(;R,r,x,y,n)
stan
functions {
vector binomial_lpmfs(
    array[] int y,
    array[] int args1,
    vector args2
) {
    return jbroadcasted_binomial_lpmfs(y, args1, args2);
}
vector jbroadcasted_binomial_lpmfs(
    array[] int x1,
    array[] int x2,
    vector x3
) {
    int n = dims(x1)[1];
    vector[n] rv;
    for(i in 1:n) {
        rv[i] = binomial_lpmfs(
            broadcasted_getindex(x1, i),
            broadcasted_getindex(x2, i),
            broadcasted_getindex(x3, i)
        );
    }
    return rv;
}
real binomial_lpmfs(
    int args1,
    int args2,
    real args3
) {
    return binomial_lpmf(args1 | args2, args3);
}
int broadcasted_getindex(array[] int x, int i) {
    return x[i];
}
real broadcasted_getindex(vector x, int i) {
    return x[i];
}
array[] int binomial_int_rng(
    int anontok__1,
    array[] int N,
    vector p
) {
    int n = anontok__1;
    if (dims(N)[1] != n) reject("binomial_rng: dim mismatch — `N` dim 1 (= ", dims(N)[1], ") does not match `n` (= ", n, "), inferred from `anontok__1` dim 1. `n` sizes: `anontok__1` dim 1 (= ", anontok__1, "), `N` dim 1 (= ", dims(N)[1], ").");
    return binomial_rng(N, p);
}
}
data {
    int x_n;
    real R;
    real r;
    vector[x_n] x;
    int y_n;
    array[y_n] int y;
    int n_n;
    array[n_n] int n;
}
transformed data {
    vector[x_n] p_threshold_angle = asin(((R - r) ./ x));
}
parameters {
    real<lower=0.0> p_sigma;
}
transformed parameters {
    vector[x_n] p = ((2 * Phi((p_threshold_angle / p_sigma))) - 1);
}
model {
    y ~ binomial(n, p);
}
generated quantities {
    real p_sigma_degrees = ((p_sigma * 180) / 3.141592653589793);
    vector[y_n] y_likelihood = binomial_lpmfs(y, n, p);
    array[n_n] int y_gen = binomial_int_rng(y_n, n, p);
}

Angle + distance model ​

The next model assumes the aiming and distance-control events are independent, so their probabilities multiply elementwise. The two submodel calls introduce distinct, hygienically scoped parameters (sigma and sigma_distance) and return vectors of the same inferred length. The parent only states how those components combine and how successes are observed.

julia
second_principles = @slic begin 
    p_angle ~ angle_submodel(;R,r,x)
    p_distance ~ distance_submodel(;distance_tolerance, overshot, x)
    p = p_angle .* p_distance
    y ~ binomial(n, p)
end
second_principles_posterior = second_principles(;R,r,x,distance_tolerance, overshot,y,n,)
stan
functions {
vector binomial_lpmfs(
    array[] int y,
    array[] int args1,
    vector args2
) {
    return jbroadcasted_binomial_lpmfs(y, args1, args2);
}
vector jbroadcasted_binomial_lpmfs(
    array[] int x1,
    array[] int x2,
    vector x3
) {
    int n = dims(x1)[1];
    vector[n] rv;
    for(i in 1:n) {
        rv[i] = binomial_lpmfs(
            broadcasted_getindex(x1, i),
            broadcasted_getindex(x2, i),
            broadcasted_getindex(x3, i)
        );
    }
    return rv;
}
real binomial_lpmfs(
    int args1,
    int args2,
    real args3
) {
    return binomial_lpmf(args1 | args2, args3);
}
int broadcasted_getindex(array[] int x, int i) {
    return x[i];
}
real broadcasted_getindex(vector x, int i) {
    return x[i];
}
array[] int binomial_int_rng(
    int anontok__1,
    array[] int N,
    vector p
) {
    int n = anontok__1;
    if (dims(N)[1] != n) reject("binomial_rng: dim mismatch — `N` dim 1 (= ", dims(N)[1], ") does not match `n` (= ", n, "), inferred from `anontok__1` dim 1. `n` sizes: `anontok__1` dim 1 (= ", anontok__1, "), `N` dim 1 (= ", dims(N)[1], ").");
    return binomial_rng(N, p);
}
}
data {
    int x_n;
    real R;
    real r;
    vector[x_n] x;
    int distance_tolerance_n;
    vector[distance_tolerance_n] distance_tolerance;
    int overshot_n;
    vector[overshot_n] overshot;
    int y_n;
    array[y_n] int y;
    int n_n;
    array[n_n] int n;
}
transformed data {
    vector[x_n] p_angle_threshold_angle = asin(((R - r) ./ x));
}
parameters {
    real<lower=0.0> p_angle_sigma;
    real<lower=0.0> p_distance_sigma_distance;
}
transformed parameters {
    vector[x_n] p_angle = ((2 * Phi((p_angle_threshold_angle / p_angle_sigma))) - 1);
    vector[x_n] p_distance = (
        Phi(((distance_tolerance - overshot) ./ ((x + overshot) * p_distance_sigma_distance))) -
        Phi(((-overshot) ./ ((x + overshot) * p_distance_sigma_distance)))
    );
    vector[x_n] p = (p_angle .* p_distance);
}
model {
    p_distance_sigma_distance ~ std_normal();
    y ~ binomial(n, p);
}
generated quantities {
    real p_angle_sigma_degrees = ((p_angle_sigma * 180) / 3.141592653589793);
    vector[y_n] y_likelihood = binomial_lpmfs(y, n, p);
    array[n_n] int y_gen = binomial_int_rng(y_n, n, p);
}

Allowing extra variation ​

The final step changes the observation layer rather than the putting physics. It models the observed proportions y ./ n with a normal approximation whose variance contains the usual binomial term p .* (1-p) ./ n plus an extra scale sigma_y^2. That additional scale can absorb variation not explained by the two physical mechanisms.

to_vector converts the integer count arrays to Stan vectors before division; the dotted arithmetic then remains elementwise. StanBlocks infers raw_proportions, both component probabilities, and the combined p as vector-valued transformed quantities. The generated Stan pane makes the change from a binomial likelihood to this continuous approximation explicit.

julia
third_principles = @slic begin 
    raw_proportions = to_vector(y) ./ to_vector(n)

    p_angle ~ angle_submodel(;R,r,x)
    p_distance ~ distance_submodel(;distance_tolerance, overshot, x)
    sigma_y ~ std_normal(;lower=0.)

    p = p_angle .* p_distance
    raw_proportions ~ normal(p, sqrt(p .* (1 - p) ./ to_vector(n) + sigma_y ^ 2))
end
third_principles_posterior = third_principles(;R,r,x,distance_tolerance,overshot,y,n)
stan
functions {
vector normal_lpdfs(
    vector obs,
    vector loc,
    vector scale
) {
    return jbroadcasted_normal_lpdfs(obs, loc, scale);
}
vector jbroadcasted_normal_lpdfs(
    vector x1,
    vector x2,
    vector x3
) {
    int n = dims(x1)[1];
    vector[n] rv;
    for(i in 1:n) {
        rv[i] = normal_lpdfs(broadcasted_getindex(x1, i), broadcasted_getindex(x2, i), broadcasted_getindex(x3, i));
    }
    return rv;
}
real normal_lpdfs(
    real args1,
    real args2,
    real args3
) {
    return normal_lpdf(args1 | args2, args3);
}
real broadcasted_getindex(vector x, int i) {
    return x[i];
}
vector normal_vector_rng(
    int anontok__1,
    vector a,
    vector b
) {
    int n = anontok__1;
    return to_vector(normal_rng(a, b));
}
}
data {
    int n_n;
    int y_n;
    array[y_n] int y;
    array[n_n] int n;
    int x_n;
    real R;
    real r;
    vector[x_n] x;
    int distance_tolerance_n;
    vector[distance_tolerance_n] distance_tolerance;
    int overshot_n;
    vector[overshot_n] overshot;
}
transformed data {
    vector[n_n] raw_proportions = (to_vector(y) ./ to_vector(n));
    vector[x_n] p_angle_threshold_angle = asin(((R - r) ./ x));
}
parameters {
    real<lower=0.0> p_angle_sigma;
    real<lower=0.0> p_distance_sigma_distance;
    real<lower=0.0> sigma_y;
}
transformed parameters {
    vector[x_n] p_angle = ((2 * Phi((p_angle_threshold_angle / p_angle_sigma))) - 1);
    vector[x_n] p_distance = (
        Phi(((distance_tolerance - overshot) ./ ((x + overshot) * p_distance_sigma_distance))) -
        Phi(((-overshot) ./ ((x + overshot) * p_distance_sigma_distance)))
    );
    vector[x_n] p = (p_angle .* p_distance);
}
model {
    p_distance_sigma_distance ~ std_normal();
    sigma_y ~ std_normal();
    raw_proportions ~ normal(p, sqrt((((p .* (1 - p)) ./ to_vector(n)) + (sigma_y ^ 2))));
}
generated quantities {
    real p_angle_sigma_degrees = ((p_angle_sigma * 180) / 3.141592653589793);
    vector[n_n] raw_proportions_likelihood = normal_lpdfs(raw_proportions, p, sqrt((((p .* (1 - p)) ./ to_vector(n)) + (sigma_y ^ 2))));
    vector[n_n] raw_proportions_gen = normal_vector_rng(n_n, p, sqrt((((p .* (1 - p)) ./ to_vector(n)) + (sigma_y ^ 2))));
}
You are viewing the dev branch. This branch may include code written with Claude Code with less human supervision. Only human-approved code is merged into main.