Building up a model of golf putting

This example follows the model-building sequence from the Stan case study Model building and expansion for golf putting. The data record the distance x, number of attempts n, and number of successful putts y at each distance. Rather than presenting five unrelated programs, the StanBlocks implementation starts with a small model and reuses it as the physical assumptions become richer.

The first two models use the shorter dataset1. The angle-and-distance models use dataset12, the row-wise concatenation of both source datasets, because the longer-distance observations are where the additional mechanism matters.

The five-model progression

Display label	Statistical model and what changes
Logistic regression	A descriptive baseline: the success log-odds are a + b * x. The flat priors make a and b free coefficients, but the model does not encode the geometry of a ball and cup.
Angle model	Replaces the linear predictor with a geometric tolerance angle, asin((R-r)/x). A putt succeeds when angular error falls inside that interval, with normally distributed error scale sigma. 2Phi(threshold/sigma)-1 is the corresponding success probability. sigma_degrees is a generated quantity for interpretation.
Angle and distance	Keeps the angle calculation and adds an independent distance-success term. The two positive components of sigma control angular and distance error. In this variant, overshoot and distance tolerance are supplied constants. The final probability is p_angle .* p_distance.
Angle, distance, and residual variation	Keeps the two-mechanism success probability but uses a normal approximation for the observed count. sigma_y adds extra-binomial variation to the count-scale standard deviation, allowing more dispersion than the binomial sampling term alone.
Estimated overshoot and tolerance	Returns to the binomial angle-and-distance likelihood, but promotes overshot and distance_tolerance from fixed constants to positive parameters with weak normal priors.

The labels in the generated-code selector are therefore descriptions of the increment added at each stage. In the Julia source, the longer historical identifiers (golf_angle_distance_2, golf_angle_distance_3_with_resids, and golf_angle_distance_4) preserve the numbering used while the family was developed.

What StanBlocks is doing

Each @slic block is a reusable model fragment. Names such as x, n, y, r, and R remain inputs until the fragment is instantiated with keyword arguments. Splatting dataset1... or dataset12... binds the data by name; the same mechanism supplies physical constants and the fixed overshoot assumptions.

Several StanBlocks features make the progression concise:

• a ~ flat() and b ~ flat() declare unconstrained parameters without adding a prior density. Keyword constraints such as lower=0 become Stan parameter constraints.

• Ordinary Julia scalar, vector, and element-wise expressions are traced into Stan. Deterministic values are placed in the earliest legal Stan block, while sigma_degrees is emitted after sampling as a generated quantity.

• Base.merge derives a model by replacing or adding statements in a quoted fragment. The angle calculation is written once and reused by all three descendants; only the assumptions that differ appear in each merge.

• The final named tuple is a family of independently instantiated models. The build evaluates this exact displayed source, calls stan_code on every member, and places each complete Stan program behind its descriptive label.

Generated Stan models

Use the StanBlocks and Generated Stan tabs to switch views. Compare side by side opens the same pair in a wide modal; the generated pane keeps all five labelled programs together.

using StanBlocks

"Downloads and preprocesses the golf datasets"
golf_data(url) = begin 
    x, n, y = eachcol(mapreduce(row->(parse.(Float64, split(row))), hcat, filter(startswith(r"[0-9]"), readlines(download(url))))')
    (;x, n=Int.(n), y=Int.(y))
end

# Fetch the first dataset
dataset1 = golf_data("https://gist.githubusercontent.com/nsiccha/553d6ce6a784142e87fbfd7aaf4c5e99/raw/668b86ddd84b45b707e594dd313f717bef1a553d/dataset1.txt")

# Fetch the second dataset
dataset2 = golf_data("https://gist.githubusercontent.com/nsiccha/13fd14e1cc5e1520aa688c1f4df2f9a7/raw/726764870b46c72911e89668eef48a22a6c13035/dataset2.txt")

# Combine the two datasets
dataset12 = map(vcat, dataset1, dataset2)


golf_logistic = @slic begin 
    a ~ flat()
    b ~ flat()
    y ~ binomial_logit(n, a + b * x)
end
golf_angle = @slic begin 
    threshold_angle = asin((R - r) ./ x)
    sigma ~ flat(;lower=0)
    p = 2 * Phi(threshold_angle / sigma) - 1
    y ~ binomial(n, p)
    sigma_degrees = sigma * 180 / pi
end
r = (1.68/2)/12
R = (4.25/2)/12
golf_angle_distance_2 = Base.merge(golf_angle, quote 
    sigma ~ normal(0, 1; n=2, lower=0.)
    sigma_angle = sigma[1]
    sigma_distance = sigma[2]
    p_angle = 2 * Phi(threshold_angle / sigma_angle) - 1
    p_distance = (
        Phi((distance_tolerance - overshot) ./ ((x + overshot) * sigma_distance))
         - Phi((-overshot) ./ ((x + overshot) * sigma_distance))
    )
    p = p_angle .* p_distance
    sigma_degrees = sigma_angle * 180 / pi
end)

overshot = 1.
distance_tolerance = 3.
golf_angle_distance_3_with_resids = Base.merge(golf_angle_distance_2, quote
    vec_n = to_vector(n)
    sigma_y ~ normal(0, 1)
    y ~ normal(vec_n .* p, vec_n .* sqrt(p .* (1 - p) ./ to_vector(n) + sigma_y ^ 2))
end)
golf_angle_distance_4 = Base.merge(golf_angle_distance_2, quote 
    overshot ~ normal(1, 5; lower=0)
    distance_tolerance ~ normal(3, 5; lower=0)
end)

golf_models = (;
    golf_logistic=golf_logistic(;dataset1...),
    golf_angle=golf_angle(;r, R, dataset1...),
    golf_angle_distance_2=golf_angle_distance_2(;r, R, overshot, distance_tolerance, dataset12...),
    golf_angle_distance_3_with_resids=golf_angle_distance_3_with_resids(;r, R, overshot, distance_tolerance, dataset12...),
    golf_angle_distance_4=golf_angle_distance_4(;r, R, overshot, distance_tolerance, dataset12...),
)

Logistic regression

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)));
}

Angle model

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] threshold_angle = asin(((R - r) ./ x));
}
parameters {
    real<lower=0> sigma;
}
transformed parameters {
    vector[x_n] p = ((2 * Phi((threshold_angle / sigma))) - 1);
}
model {
    y ~ binomial(n, p);
}
generated quantities {
    vector[y_n] y_likelihood = binomial_lpmfs(y, n, p);
    array[n_n] int y_gen = binomial_int_rng(y_n, n, p);
    real sigma_degrees = ((sigma * 180) / 3.141592653589793);
}

Angle and distance

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;
    real distance_tolerance;
    real overshot;
    int y_n;
    array[y_n] int y;
    int n_n;
    array[n_n] int n;
}
transformed data {
    vector[x_n] threshold_angle = asin(((R - r) ./ x));
}
parameters {
    vector<lower=0.0>[2] sigma;
}
transformed parameters {
    real sigma_angle = sigma[1];
    vector[x_n] p_angle = ((2 * Phi((threshold_angle / sigma_angle))) - 1);
    real sigma_distance = sigma[2];
    vector[x_n] p_distance = (
        Phi(((distance_tolerance - overshot) ./ ((x + overshot) * sigma_distance))) -
        Phi(((-overshot) ./ ((x + overshot) * sigma_distance)))
    );
    vector[x_n] p = (p_angle .* p_distance);
}
model {
    sigma ~ normal(0, 1);
    y ~ binomial(n, p);
}
generated quantities {
    vector[y_n] y_likelihood = binomial_lpmfs(y, n, p);
    array[n_n] int y_gen = binomial_int_rng(y_n, n, p);
    real sigma_degrees = ((sigma_angle * 180) / 3.141592653589793);
}

Angle, distance, and residual variation

functions {
vector normal_lpdfs(
    array[] int obs,
    vector loc,
    vector scale
) {
    return jbroadcasted_normal_lpdfs(obs, loc, scale);
}
vector jbroadcasted_normal_lpdfs(
    array[] int 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(
    int args1,
    real args2,
    real args3
) {
    return normal_lpdf(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[] real normal_int_rng(
    int anontok__1,
    vector a,
    vector b
) {
    int n = anontok__1;
    return normal_rng(a, b);
}
}
data {
    int x_n;
    real R;
    real r;
    vector[x_n] x;
    real distance_tolerance;
    real overshot;
    int n_n;
    array[n_n] int n;
    int y_n;
    array[y_n] int y;
}
transformed data {
    vector[x_n] threshold_angle = asin(((R - r) ./ x));
    vector[n_n] vec_n = to_vector(n);
}
parameters {
    vector<lower=0.0>[2] sigma;
    real sigma_y;
}
transformed parameters {
    real sigma_angle = sigma[1];
    vector[x_n] p_angle = ((2 * Phi((threshold_angle / sigma_angle))) - 1);
    real sigma_distance = sigma[2];
    vector[x_n] p_distance = (
        Phi(((distance_tolerance - overshot) ./ ((x + overshot) * sigma_distance))) -
        Phi(((-overshot) ./ ((x + overshot) * sigma_distance)))
    );
    vector[x_n] p = (p_angle .* p_distance);
}
model {
    sigma ~ normal(0, 1);
    sigma_y ~ normal(0, 1);
    y ~ normal((vec_n .* p), (vec_n .* sqrt((((p .* (1 - p)) ./ to_vector(n)) + (sigma_y ^ 2)))));
}
generated quantities {
    vector[y_n] y_likelihood = normal_lpdfs(y, (vec_n .* p), (vec_n .* sqrt((((p .* (1 - p)) ./ to_vector(n)) + (sigma_y ^ 2)))));
    array[y_n] real y_gen = normal_int_rng(
        y_n,
        (vec_n .* p),
        (vec_n .* sqrt((((p .* (1 - p)) ./ to_vector(n)) + (sigma_y ^ 2))))
    );
    real sigma_degrees = ((sigma_angle * 180) / 3.141592653589793);
}

Estimated overshoot and tolerance

functions {
real normal_lpdfs(
    real args1,
    int args2,
    int args3
) {
    return normal_lpdf(args1 | args2, args3);
}
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;
    real overshot;
    real distance_tolerance;
    int y_n;
    array[y_n] int y;
    int n_n;
    array[n_n] int n;
}
transformed data {
    vector[x_n] threshold_angle = asin(((R - r) ./ x));
}
parameters {
    vector<lower=0.0>[2] sigma;
}
transformed parameters {
    real sigma_angle = sigma[1];
    vector[x_n] p_angle = ((2 * Phi((threshold_angle / sigma_angle))) - 1);
    real sigma_distance = sigma[2];
    vector[x_n] p_distance = (
        Phi(((distance_tolerance - overshot) ./ ((x + overshot) * sigma_distance))) -
        Phi(((-overshot) ./ ((x + overshot) * sigma_distance)))
    );
    vector[x_n] p = (p_angle .* p_distance);
}
model {
    sigma ~ normal(0, 1);
    overshot ~ normal(1, 5);
    distance_tolerance ~ normal(3, 5);
    y ~ binomial(n, p);
}
generated quantities {
    real overshot_likelihood = normal_lpdfs(overshot, 1, 5);
    real overshot_gen = normal_rng(1, 5);
    real distance_tolerance_likelihood = normal_lpdfs(distance_tolerance, 3, 5);
    real distance_tolerance_gen = normal_rng(3, 5);
    vector[y_n] y_likelihood = binomial_lpmfs(y, n, p);
    array[n_n] int y_gen = binomial_int_rng(y_n, n, p);
    real sigma_degrees = ((sigma_angle * 180) / 3.141592653589793);
}