StanBlocks.jl feature atlas
StanCon 2026 companion: author source → corresponding Stan emission
This is the long-form companion to the talk. It inventories the current author-visible language surface reviewed against StanBlocks commit e6f607d (2026-08-04). The presentation files themselves are newer; the compiler source covered here is unchanged by those documentation commits.
Every StanBlocks authoring example is followed by its corresponding Stan emission. Small qualitative examples isolate the exact lowering under discussion; the end-to-end comparison panels reproduce the complete stan_code(model) result, including helper functions for vectorised pointwise densities, sized RNGs, ragged carriers, lifted closures, and distribution combinators. Descriptor and execution snippets are Julia API calls rather than transpilation examples and are labelled separately.
This atlas covers language and workflow features, not every registered Stan builtin. Several hundred Stan functions and distribution signatures are registered; the source of truth for those is src/slic_stan/builtin.jl.
Overview
The StanCon 2026 "quality-of-life improvements" from the opening slide, in that order — each with the one detail worth remembering and a link to its worked example: 2. Activity analysis — A backward likelihood-reachability pass marks every binding that can affect an observed quantity, which is what lets the compiler choose transformed data vs transformed parameters placement and move prior-only computations into generated quantities. → Core model
Submodels — Reusable model fragments (anonymous
@slic (…)blocks and named@slic f(…)=…functions) compose into larger models with automatic parameter namespacing (e.g.eta_beta). → CompositionUser-defined types —
@usertypelets author-defined structured types (currently NamedTuples and ragged containers) flow through tracing, data binding, and Stan emission as first-class shapes. → Structured dataHigher-order user-defined functions — Functions (distributions included) can take and return other functions, so combinators like
weighted/truncated/censoredand custom HOFs are authored once and reused. → Function signaturesAutomatic but optional type & shape inference — The compiler infers scalar types, array shapes, and constraints from data and usage, but you can always pin them explicitly (typed LHS, sized signatures). → Core model
(Ragged) plates —
plate(…) dois a compiler-owned sampling loop that promotes fresh per-cell variables to outer storage and supports ragged (uneven-length) cells without hand-written index bookkeeping. → Structured dataMetaprogramming — Caller-side macros expand before tracing, and
@inlinehelpers, trailing-!mutation,@stan_assert, andreturn_type_ofqueries let you compute model structure at transpile time. → Expansion toolsJulia-style multiple dispatch — Named submodels and helpers dispatch on positional argument types (variadic and function-typed included), so one name resolves to different fragments by the shapes it is called with. → Typed dispatch
…and more — Generated observation outputs (pointwise log-likelihood, predictive draws, missing-outcome imputation), custom distribution triads, fused GLMs, post-hoc
Base.mergevariants and cross-validation taint, executable descriptors, and the scientific-computing surface (ODEs, Torsten, GP,reduce_sum) — plus the honest current boundaries, in the sections below the capability walkthrough.
Feature map
| Area | Current surface | Detailed MWE |
|---|---|---|
| Model declarations | @slic, ~, =, typed LHS, bare flat-prior parameters | Core model |
| Analysis | block placement, type/shape/constraint inference, likelihood activity | Core model |
| Composition | anonymous and named submodels, typed dispatch, Base.merge, data rebinding | Composition |
| Structured data/models | RaggedVector, ragged constraints, EachCol/EachRow, fancy indexing, plate | Structured models |
| Julia-like functions | closures, HOFs, varargs, positional defaults, required/optional kwargs | Function signatures |
| Distribution abstraction | custom _lpdf/_lpdfs/_rng, @lpxf, @lhs, weighted/truncated/censored/interval | Distributions |
| Metaprogramming | caller macros, @inline, trailing !, @stan_assert, return_type_of | Expansion tools |
| Deterministic helpers | @deffun, loops, branches, mutation, comprehensions, iteration | @deffun |
| Observation workflow | pointwise log likelihoods, predictive draws, missing-outcome imputation | Generated outputs |
| Model variants | post-hoc overrides and lower-level cross-validation taint | Variants |
| Reflection/execution | descriptors, definition closure, derived operations, BridgeStan | Descriptor |
| Scientific models | ODE solvers, Torsten signatures, GP covariance, reduce_sum, fused GLMs | Scientific surface |
One model: block placement, types, shapes, and constraints
StanBlocks source
using StanBlocks
m = @slic (; x = [-1.0, 0.0, 1.0], y = [0.2, -0.1, 0.7]) begin
mx = mean(x)
alpha ~ normal(0, 2)
beta ~ normal(0, 1)
sigma ~ exponential(1)
mu = alpha + beta * (x - mx)
y ~ normal(mu, sigma)
endfunctions {
vector normal_lpdfs(
vector obs,
vector loc,
real scale
) {
return jbroadcasted_normal_lpdfs(obs, loc, scale);
}
vector jbroadcasted_normal_lpdfs(
vector x1,
vector x2,
real 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), x3);
}
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,
real b
) {
int n = anontok__1;
return to_vector(normal_rng(a, b));
}
}
data {
int x_n;
vector[x_n] x;
int y_n;
vector[y_n] y;
}
transformed data {
real mx = mean(x);
}
parameters {
real alpha;
real beta;
real<lower=0.0> sigma;
}
transformed parameters {
vector[x_n] mu = (alpha + (beta * (x - mx)));
}
model {
alpha ~ normal(0, 2);
beta ~ normal(0, 1);
sigma ~ exponential(1);
y ~ normal(mu, sigma);
}
generated quantities {
vector[y_n] y_likelihood = normal_lpdfs(y, mu, sigma);
vector[y_n] y_gen = normal_vector_rng(y_n, mu, sigma);
}What happened:
Julia vectors
xandyestablished Stan vector types and generated their data dimensionsx_nandy_n.mxdepends only on data, so it moved totransformed data.The prior family inferred
sigma's lower bound.mudepends on parameters and feeds the likelihood, so it lives intransformed parameters.The observed
ygenerated pointwise log-likelihood and posterior-predictive twins automatically.
Typed LHS and a prior-free parameter
Use explicit types when the prior call does not determine the desired container, or when the type itself carries Stan semantics:
X = [1.0 0.0; 1.0 1.0; 1.0 2.0]
y = [0.1, -0.2, 0.4]
typed_model = @slic (; X, y, k = size(X, 2)) begin
beta::vector[k] # improper-flat sampler parameter
sigma::real ~ exponential(1)
y ~ normal(X * beta, sigma)
endfunctions {
vector normal_lpdfs(
vector obs,
vector loc,
real scale
) {
return jbroadcasted_normal_lpdfs(obs, loc, scale);
}
vector jbroadcasted_normal_lpdfs(
vector x1,
vector x2,
real 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), x3);
}
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,
real b
) {
int n = anontok__1;
return to_vector(normal_rng(a, b));
}
}
data {
int k;
int y_n;
vector[y_n] y;
int X_m;
int X_n;
matrix[X_m, X_n] X;
}
transformed data {
}
parameters {
vector[k] beta;
real<lower=0.0> sigma;
}
transformed parameters {
}
model {
sigma ~ exponential(1);
y ~ normal((X * beta), sigma);
}
generated quantities {
vector[y_n] y_likelihood = normal_lpdfs(y, (X * beta), sigma);
vector[y_n] y_gen = normal_vector_rng(y_n, (X * beta), sigma);
}A bare typed declaration is allowed only in a model/submodel and means a parameter with no density statement. Inside @deffun, the same syntax declares a local that the function must fill. Native constrained types use the same LHS surface, for example L::cholesky_factor_corr[k] ~ lkj_corr_cholesky(2) when L participates in a downstream likelihood.
Composition and model families
Anonymous submodels: inputs by name
X = [1.0 0.0; 1.0 1.0; 1.0 2.0]
y = [0.1, -0.2, 0.4]
population_effect = @slic begin
k = dims(X)[2]
beta ~ normal(0, 1; n = k)
return X * beta
end
m = @slic (; X, y) begin
eta ~ population_effect(; X)
y ~ normal(eta, 1)
endfunctions {
vector normal_lpdfs(
vector obs,
vector loc,
int scale
) {
return jbroadcasted_normal_lpdfs(obs, loc, scale);
}
vector jbroadcasted_normal_lpdfs(
vector x1,
vector x2,
int 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), x3);
}
return rv;
}
real normal_lpdfs(
real args1,
real args2,
int 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,
int b
) {
int n = anontok__1;
return to_vector(normal_rng(a, b));
}
}
data {
int X_n;
int X_m;
matrix[X_m, X_n] X;
int y_n;
vector[y_n] y;
}
transformed data {
int eta_k = X_n;
}
parameters {
vector[eta_k] eta_beta;
}
transformed parameters {
vector[X_m] eta = (X * eta_beta);
}
model {
eta_beta ~ normal(0, 1);
y ~ normal(eta, 1);
}
generated quantities {
vector[y_n] y_likelihood = normal_lpdfs(y, eta, 1);
vector[y_n] y_gen = normal_vector_rng(y_n, eta, 1);
}Submodel parameters are namespaced under the receiving LHS. Anonymous submodels receive inputs through kwargs or enclosing scope; they deliberately have no positional calling form.
Named submodels: typed positional dispatch
x = [-1.0, 0.0, 1.0]
y = [0.1, -0.2, 0.4]
@slic random_slope(x::vector[n]) = begin
beta ~ normal(0, 1)
return beta * x
end
@slic offset(x::vector[n]) = begin
alpha ~ normal(0, 2)
return alpha + x
end
m = @slic (; x, y) begin
eta1 ~ random_slope(x)
eta2 ~ offset(eta1)
y ~ normal(eta2, 1)
endfunctions {
vector normal_lpdfs(
vector obs,
vector loc,
int scale
) {
return jbroadcasted_normal_lpdfs(obs, loc, scale);
}
vector jbroadcasted_normal_lpdfs(
vector x1,
vector x2,
int 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), x3);
}
return rv;
}
real normal_lpdfs(
real args1,
real args2,
int 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,
int b
) {
int n = anontok__1;
return to_vector(normal_rng(a, b));
}
}
data {
int x_n;
vector[x_n] x;
int y_n;
vector[y_n] y;
}
transformed data {
}
parameters {
real eta1_beta;
real eta2_alpha;
}
transformed parameters {
vector[x_n] eta1 = (eta1_beta * x);
vector[x_n] eta2 = (eta2_alpha + eta1);
}
model {
eta1_beta ~ normal(0, 1);
eta2_alpha ~ normal(0, 2);
y ~ normal(eta2, 1);
}
generated quantities {
vector[y_n] y_likelihood = normal_lpdfs(y, eta2, 1);
vector[y_n] y_gen = normal_vector_rng(y_n, eta2, 1);
}Typed arguments select Julia methods on the traced Stan center type. Multiple definitions of @slic f(::real) and @slic f(::int) are distinct methods; untyped positional arguments match any traced value.
Post-hoc variants with Base.merge
base = @slic begin
beta ~ normal(0, 1)
y ~ normal(beta * x, 1)
end
wide = Base.merge(base, quote
beta ~ normal(0, 5)
end)
wide_model = wide(; x=[-1.0, 0.0, 1.0], y=[0.1, -0.2, 0.4])
fixed_model = Base.merge(base, (; beta=0.25))(
; x=[-1.0, 0.0, 1.0], y=[0.1, -0.2, 0.4],
)functions {
vector normal_lpdfs(
vector obs,
vector loc,
int scale
) {
return jbroadcasted_normal_lpdfs(obs, loc, scale);
}
vector jbroadcasted_normal_lpdfs(
vector x1,
vector x2,
int 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), x3);
}
return rv;
}
real normal_lpdfs(
real args1,
real args2,
int 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,
int b
) {
int n = anontok__1;
return to_vector(normal_rng(a, b));
}
}
data {
int y_n;
vector[y_n] y;
int x_n;
vector[x_n] x;
}
transformed data {
}
parameters {
real beta;
}
transformed parameters {
}
model {
beta ~ normal(0, 5);
y ~ normal((beta * x), 1);
}
generated quantities {
vector[y_n] y_likelihood = normal_lpdfs(y, (beta * x), 1);
vector[y_n] y_gen = normal_vector_rng(y_n, (beta * x), 1);
}The merged statement matches by the bare LHS name and replaces the original. Typed and untyped versions of the same LHS still match. Additional names append new statements. A merged NamedTuple, such as (; beta=0.25), instead removes the matching statement and stores the supplied value as data: this fixes the name without retaining its prior as a likelihood. Statement replacements and fixed bindings may be combined in one Base.merge call. Ordinary model kwargs only bind data and do not remove statements. A merged SlicModel value can also be interpolated into the callee position of generated AST.
Structured data and compiler-owned plates
plate: an independent-cell loop that may introduce parameters
StanBlocks source:
plate_model = @slic (; y = randn(4), mu0 = 0.5) begin
sigma ~ exponential(1)
theta ~ plate(y; outer = 4) do yi
t ~ normal(mu0, 1)
yi ~ normal(t, sigma)
t
end
endfunctions {
real normal_lpdfs(
real args1,
real args2,
int args3
) {
return normal_lpdf(args1 | args2, args3);
}
real normal_lpdfs(
real args1,
real args2,
real args3
) {
return normal_lpdf(args1 | args2, args3);
}
}
data {
int y_n;
real mu0;
vector[y_n] y;
}
transformed data {
}
parameters {
real<lower=0.0> sigma;
vector[4] theta_t;
}
transformed parameters {
vector[4] theta;
for(plate_i__pl_1 in 1:4) {
theta[plate_i__pl_1] = theta_t[plate_i__pl_1];
}
}
model {
sigma ~ exponential(1);
for(plate_i__pl_1 in 1:4) {
theta_t[plate_i__pl_1] ~ normal(mu0, 1);
y[plate_i__pl_1] ~ normal(theta_t[plate_i__pl_1], sigma);
}
}
generated quantities {
vector[y_n] y_gen;
for(plate_i__pl_1 in 1:4) {
y_gen[plate_i__pl_1] = normal_rng(theta_t[plate_i__pl_1], sigma);
}
}Positional inputs are sliced; lexical captures are shared; fresh cell bindings are promoted to outer storage; the trailing expression becomes the cell output. Fixed vector cells become columns of a matrix. Additional outer axes add Stan array prefixes. Cells must remain independent: plate is not a scan.
Vector-valued plate cells
vector_plate_model = @slic (; n_groups = 8, k = 3) begin
L::cholesky_factor_corr[k] ~ lkj_corr_cholesky(2)
tau::vector[k] ~ normal(0, 1; lower=0)
b::vector[k] ~ plate(; outer=n_groups) do g
z::vector[k] ~ std_normal()
diag_pre_multiply(tau, L) * z
end
endfunctions {
vector std_normal_vector_rng(
int anontok__1
) {
int n = anontok__1;
return to_vector(normal_rng(rep_vector(0, n), 1));
}
}
data {
int k;
int n_groups;
}
transformed data {
}
parameters {
cholesky_factor_corr[k] L;
vector<lower=0>[k] tau;
matrix[k, n_groups] b_z;
}
transformed parameters {
matrix[k, k] b__pl_inv1_1 = diag_pre_multiply(tau, L);
matrix[k, n_groups] b;
for(plate_i__pl_1 in 1:n_groups) {
b[:, plate_i__pl_1] = (b__pl_inv1_1 * b_z[:, plate_i__pl_1]);
}
}
model {
L ~ lkj_corr_cholesky(2);
tau ~ normal(0, 1);
for(plate_i__pl_1 in 1:n_groups) {
b_z[:, plate_i__pl_1] ~ std_normal();
}
}
generated quantities {
}Both z and b have logical shape matrix[k,n_groups]; the per-group vectors occupy columns.
Ragged data
groups = [[1.0, 2.0], [3.0], [4.0, 5.0, 6.0]]
ragged_model = @slic (; groups, g=2, y=0.0) begin
first = groups[1]
selected = groups[g]
ng = length(groups)
y ~ normal(sum(first) + sum(selected) + ng, 1)
endfunctions {
int ragged_length_RaggedVector(
tuple(vector, array[] int) x,
int i
) {
return ((ragged_end_RaggedVector(x, i) - ragged_start_RaggedVector(x, i)) + 1);
}
int ragged_end_RaggedVector(tuple(vector, array[] int) x, int i) {
return x.2[i];
}
int ragged_start_RaggedVector(
tuple(vector, array[] int) x,
int i
) {
if((i == 1)) {
return 1;
} else {
return (1 + x.2[(i - 1)]);
}
}
vector getindex_RaggedVector(
tuple(vector, array[] int) rv,
int i
) {
return rv.1[ragged_start_RaggedVector(rv, i):ragged_end_RaggedVector(rv, i)];
}
int num_elements_RaggedVector(tuple(vector, array[] int) rv) {
return size(rv.2);
}
real normal_lpdfs(
real args1,
real args2,
int args3
) {
return normal_lpdf(args1 | args2, args3);
}
}
data {
int groups_mem_n;
int groups_ends_n;
tuple(vector[groups_mem_n], array[groups_ends_n] int) groups;
int g;
real y;
}
transformed data {
vector[ragged_length_RaggedVector(groups, 1)] first = getindex_RaggedVector(groups, 1);
vector[ragged_length_RaggedVector(groups, g)] selected = getindex_RaggedVector(groups, g);
int ng = num_elements_RaggedVector(groups);
}
parameters {
}
transformed parameters {
}
model {
y ~ normal((sum(first) + sum(selected) + ng), 1);
}
generated quantities {
real y_likelihood = normal_lpdfs(y, (sum(first) + sum(selected) + ng), 1);
real y_gen = normal_rng((sum(first) + sum(selected) + ng), 1);
}The nested Julia vectors become a nominal carrier with flat mem and inclusive group ends. groups[g] lowers to a Stan helper that reconstructs the requested slice; length(groups) counts groups, not scalar elements.
Ragged constrained parameters
ragged_constraint_model = @slic (; Ks=[2, 3, 4], y=0.3) begin
p::simplex[Ks] ~ flat()
y ~ normal(sum(p[1]), 0.1)
endfunctions {
array[] int jbroadcasted_sub(
array[] int x1,
int x2
) {
int n = dims(x1)[1];
array[n] int rv;
for(i in 1:n) {
rv[i] = (broadcasted_getindex(x1, i) - x2);
}
return rv;
}
int broadcasted_getindex(array[] int x, int i) {
return x[i];
}
real normal_lpdfs(
real args1,
real args2,
real args3
) {
return normal_lpdf(args1 | args2, args3);
}
int ragged_end(array[] int ends, int i) {
return ends[i];
}
int ragged_start(
array[] int ends,
int i
) {
if((i == 1)) {
return 1;
} else {
return (1 + ends[(i - 1)]);
}
}
}
data {
int Ks_n;
array[Ks_n] int Ks;
real y;
}
transformed data {
array[Ks_n] int c_end__rc_1 = cumulative_sum(Ks);
array[Ks_n] int f_end__rc_1 = cumulative_sum(jbroadcasted_sub(Ks, 1));
}
parameters {
vector[sum(jbroadcasted_sub(Ks, 1))] p_free__rc_1;
}
transformed parameters {
vector[sum(Ks)] p_mem__rc_1;
for(g__rc_1 in 1:num_elements(Ks)) {
p_mem__rc_1[((c_end__rc_1[g__rc_1] - Ks[g__rc_1]) + 1):c_end__rc_1[g__rc_1]] = simplex_jacobian(
p_free__rc_1[((f_end__rc_1[g__rc_1] - (Ks[g__rc_1] - 1)) + 1):f_end__rc_1[g__rc_1]]
);
}
}
model {
y ~ normal(sum(p_mem__rc_1[ragged_start(c_end__rc_1, 1):ragged_end(c_end__rc_1, 1)]), 0.1);
}
generated quantities {
real y_likelihood = normal_lpdfs(y, sum(p_mem__rc_1[ragged_start(c_end__rc_1, 1):ragged_end(c_end__rc_1, 1)]), 0.1);
real y_gen = normal_rng(sum(p_mem__rc_1[ragged_start(c_end__rc_1, 1):ragged_end(c_end__rc_1, 1)]), 0.1);
}Stan cannot declare a runtime-ragged pile of simplices. StanBlocks emits a flat unconstrained parameter, then a compiler-owned per-group constrain/Jacobian loop and a logical ragged view. The same route covers ordered, positive_ordered, cholesky_factor_corr, and square cholesky_factor_cov groups. This feature requires BridgeStan 2.9 / Stan 2.39 or newer.
Dense views and indexing
@deffun col_sums(X::matrix[m,k]) = [sum(c) for c in EachCol(X)]
views_model = @slic (; X=[1.0 2.0; 3.0 4.0], y=[0.0, 0.0]) begin
sums = col_sums(X)
y ~ normal(sums, 1)
endfunctions {
vector col_sums(
matrix X
) {
vector[(1 + (cols(X) - 1))] comprehension_result__lc_4;
for(value_index__vi_3 in 1:cols(X)) {
comprehension_result__lc_4[value_index__vi_3] = sum(col(X, value_index__vi_3));
}
return comprehension_result__lc_4;
}
vector normal_lpdfs(
vector obs,
vector loc,
int scale
) {
return jbroadcasted_normal_lpdfs(obs, loc, scale);
}
vector jbroadcasted_normal_lpdfs(
vector x1,
vector x2,
int 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), x3);
}
return rv;
}
real normal_lpdfs(
real args1,
real args2,
int 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,
int b
) {
int n = anontok__1;
return to_vector(normal_rng(a, b));
}
}
data {
int X_m;
int X_n;
matrix[X_m, X_n] X;
int y_n;
vector[y_n] y;
}
transformed data {
vector[(1 + (cols(X) - 1))] sums = col_sums(X);
}
parameters {
}
transformed parameters {
}
model {
y ~ normal(sums, 1);
}
generated quantities {
vector[y_n] y_likelihood = normal_lpdfs(y, sums, 1);
vector[y_n] y_gen = normal_vector_rng(y_n, sums, 1);
}| Source | Stan lowering |
|---|---|
EachCol(X)[j] | col(X,j) |
EachRow(X)[i] | row(X,i) |
length(EachCol(X)) | cols(X) |
length(EachRow(X)) | rows(X) |
alpha[county_idx] | vectorised integer-array indexing |
X[i,:,:], X[:,:,k], L[i,:,:] | resolved multi-colon slice |
Function signatures, defaults, keywords, varargs, and higher-order functions
Positional defaults and keyword arguments
@deffun affine(x::real, a::real = 2.0)::real = a * x + 1.0
@deffun scale(x::vector[n]; factor, bias)::vector[n] = factor * x + bias
@deffun shift(x::vector[n]; a, b = 2.0)::vector[n] = a * x + b
signature_model = @slic (; x = [1.0, 2.0], y = [0.0, 0.0]) begin
a1 = affine(3.0) # default a=2
a2 = affine(3.0, 4.0)
s = scale(x; factor=0.5, bias=0.0) # required kwargs
z = shift(s; a=1.5) # optional b defaults to 2
y ~ normal(z + a1 + a2, 1)
endfunctions {
real affine(real x, real a) {
return ((a * x) + 1.0);
}
vector kwcall_scale(
tuple(real, real) kw,
vector x
) {
int n = dims(x)[1];
real factor = kw.1;
real bias = kw.2;
return ((factor * x) + bias);
}
vector kwcall_shift(
tuple(real, real) kw,
vector x
) {
int n = dims(x)[1];
real a = kw.1;
real b = kw.2;
return ((a * x) + b);
}
vector normal_lpdfs(
vector obs,
vector loc,
int scale
) {
return jbroadcasted_normal_lpdfs(obs, loc, scale);
}
vector jbroadcasted_normal_lpdfs(
vector x1,
vector x2,
int 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), x3);
}
return rv;
}
real normal_lpdfs(
real args1,
real args2,
int 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,
int b
) {
int n = anontok__1;
return to_vector(normal_rng(a, b));
}
}
data {
int x_n;
vector[x_n] x;
int y_n;
vector[y_n] y;
}
transformed data {
real a1 = affine(3.0, 2.0);
real a2 = affine(3.0, 4.0);
vector[x_n] s = kwcall_scale((0.5, 0.0), x);
vector[x_n] z = kwcall_shift((1.5, 2.0), s);
}
parameters {
}
transformed parameters {
}
model {
y ~ normal((z + a1 + a2), 1);
}
generated quantities {
vector[y_n] y_likelihood = normal_lpdfs(y, (z + a1 + a2), 1);
vector[y_n] y_gen = normal_vector_rng(y_n, (z + a1 + a2), 1);
}Defaults are resolved at the call site; Stan receives a specialised call with ordinary positional arguments. Omitting a required keyword errors during tracing rather than emitting an under-specified function call.
Variadic and function-typed dispatch
import StanBlocks.stan: logit
@deffun begin
my_bernoulli_lpmf(y, ::typeof(logit), eta) =
bernoulli_logit_lpmf(y, eta)
my_bernoulli_lpmf(y, ::typeof(log), eta) =
bernoulli_lpmf(y, exp(eta))
my_bernoulli_lpmfs(y::int, args...) =
my_bernoulli_lpmf(y, args...)
apply_twice(f, x::real)::real = f(f(x))
end
dispatch_model = @slic (; y=1, eta=0.2, obs=0.0) begin
ll = my_bernoulli_lpmfs(y, logit, eta)
z = apply_twice(exp, eta)
obs ~ normal(ll + z, 1)
endfunctions {
real my_bernoulli_logit_lpmfs(int y, real args1) {
return my_bernoulli_logit_lpmf(y | args1);
}
real my_bernoulli_logit_lpmf(int y, real eta) {
return bernoulli_logit_lpmf(y | eta);
}
real apply_twice_exp(real x) {
return exp(exp(x));
}
real normal_lpdfs(
real args1,
real args2,
int args3
) {
return normal_lpdf(args1 | args2, args3);
}
}
data {
int y;
real eta;
real obs;
}
transformed data {
real ll = my_bernoulli_logit_lpmfs(y, eta);
real z = apply_twice_exp(eta);
}
parameters {
}
transformed parameters {
}
model {
obs ~ normal((ll + z), 1);
}
generated quantities {
real obs_likelihood = normal_lpdfs(obs, (ll + z), 1);
real obs_gen = normal_rng((ll + z), 1);
}Function values are compile-time tokens. StanBlocks selects and emits the specialised method; no runtime function object reaches Stan.
Closures
ts = [0.5, 1.0]
yobs = 0.6
closure_model = @slic (; ts, yobs) begin
lambda ~ exponential(1)
y = ode_rk45(
(t, state) -> -lambda * state,
[1.0], 0.0, to_array_1d(ts),
)
yobs ~ normal(y[1][1], 0.05)
endfunctions {
// lifted closure (id 1)
vector closure_1(
real t,
vector state,
real lambda
) {
return ((-lambda) * state);
}
real normal_lpdfs(
real args1,
real args2,
real args3
) {
return normal_lpdf(args1 | args2, args3);
}
}
data {
int ts_n;
vector[ts_n] ts;
real yobs;
}
transformed data {
}
parameters {
real<lower=0.0> lambda;
}
transformed parameters {
array[ts_n] vector[1] y = ode_rk45(closure_1, [1.0]', 0.0, to_array_1d(ts), lambda);
}
model {
lambda ~ exponential(1);
yobs ~ normal(y[1][1], 0.05);
}
generated quantities {
real yobs_likelihood = normal_lpdfs(yobs, y[1][1], 0.05);
real yobs_gen = normal_rng(y[1][1], 0.05);
}The closure is lifted into a generated Stan function. Captured data and parameters become explicit trailing arguments, and likelihood activity follows through those captures so lambda stays a live parameter.
Custom and higher-order distributions
Custom distribution triad
@deffun begin
@lhs @lpxf robust_lpdf(y::real, mu::real, sigma::real) =
student_t_lpdf(y, 4, mu, sigma)
robust_lpdfs(y::real, mu::real, sigma::real)::real =
robust_lpdf(y, mu, sigma)
robust_rng(mu::real, sigma::real)::real =
student_t_rng(4, mu, sigma)
end
custom_model = @slic (; y = 0.2) begin
mu ~ normal(0, 2)
sigma ~ exponential(1)
y ~ robust(mu, sigma)
endfunctions {
real robust_lpdf(real y, real mu, real sigma) {
return student_t_lpdf(y | 4, mu, sigma);
}
real robust_lpdfs(real y, real mu, real sigma) {
return robust_lpdf(y | mu, sigma);
}
real robust_rng(real mu, real sigma) {
return student_t_rng(4, mu, sigma);
}
}
data {
real y;
}
transformed data {
}
parameters {
real mu;
real<lower=0.0> sigma;
}
transformed parameters {
}
model {
mu ~ normal(0, 2);
sigma ~ exponential(1);
y ~ robust(mu, sigma);
}
generated quantities {
real y_likelihood = robust_lpdfs(y, mu, sigma);
real y_gen = robust_rng(mu, sigma);
}| Companion | Used for |
|---|---|
robust_lpdf | joint sampling statement |
robust_lpdfs | pointwise generated log likelihood |
robust_rng | posterior-predictive draw |
@lpxf | density/pointwise/RNG name registration |
@lhs | base-level LHS type inference for y ~ robust(...) |
This MWE deliberately declares a scalar observation. A vector observation needs a matching robust_lpdf(y::vector[n], ...) and pointwise companion. A ragged observation additionally needs the sized-token RNG form robust_rng(vector[n], args...)::vector[n].
Distribution higher-order functions
Each combinator takes a base distribution token and specializes a complete Stan distribution family. Keeping the examples separate makes the generated density, pointwise companion, predictive RNG, and transitive helper closure visible for each semantic choice. At regular width, switch between the StanBlocks and generated-Stan tabs; use Compare side by side for the full-width view.
Weighted observations
Use a data weight to scale an observation's log-density contribution without changing the posterior-predictive distribution.
Qualitatively:
y ~ weighted(normal, weight, mu, sigma)using StanBlocks
weighted_model = @slic (; y=0.1, w=2.0) begin
mu ~ normal(0.0, 1.0)
y ~ weighted(normal, w, mu, 1.0)
endfunctions {
real weighted_normal_lpdf(
real y,
real weight,
real args1,
real args2
) {
return (weight * normal_lpdf(y | args1, args2));
}
real weighted_normal_lpdfs(
real y,
real weight,
real args1,
real args2
) {
return (weight * normal_lpdfs(y, args1, args2));
}
real normal_lpdfs(
real args1,
real args2,
real args3
) {
return normal_lpdf(args1 | args2, args3);
}
real weighted_normal_rng(
real weight,
real args1,
real args2
) {
return normal_rng(args1, args2);
}
}
data {
real y;
real w;
}
transformed data {
}
parameters {
real mu;
}
transformed parameters {
}
model {
mu ~ normal(0.0, 1.0);
y ~ weighted_normal(w, mu, 1.0);
}
generated quantities {
real y_likelihood = weighted_normal_lpdfs(y, w, mu, 1.0);
real y_gen = weighted_normal_rng(w, mu, 1.0);
}The compiler specializes the normal token into weighted_normal_lpdf, its pointwise weighted_normal_lpdfs companion, and weighted_normal_rng. The RNG deliberately ignores the weight: weighting changes evidence, not the data-generating process.
Truncated observations
Truncation conditions a latent draw to lie inside the supplied bounds. Its normalizing constant belongs in the density, and prediction must draw from the same conditioned distribution.
Qualitatively:
y ~ truncated(normal, mu, sigma; lower=lo, upper=hi)using StanBlocks
truncated_model = @slic (; y=0.2, lo=-1.0, hi=1.0) begin
mu ~ normal(0.0, 1.0)
y ~ truncated(normal, mu, 1.0; lower=lo, upper=hi)
endfunctions {
real conditioning_normal_lpdf(
real y,
real lo,
real hi,
real args1,
real args2
) {
array[1] real rv;
rv[1] = negative_infinity();
if((lo >= hi)) {
reject("truncated: lower bound must be less than upper bound");
} else {
if((y >= lo)) {
if((y <= hi)) {
rv[1] = (
normal_lpdf(y | args1, args2) -
log_diff_exp(normal_lcdf_stable(hi, args1, args2), normal_lcdf_stable(lo, args1, args2))
);
}
}
}
return rv[1];
}
real normal_lcdf_stable(
real x,
real loc,
real scale
) {
return (log(erfc(((-(x - loc)) / (scale * sqrt(2.0))))) - log(2.0));
}
real conditioning_normal_lpdfs(
real y,
real lo,
real hi,
real args1,
real args2
) {
return conditioning_normal_lpdf(y | lo, hi, args1, args2);
}
real conditioning_normal_rng(
real lo,
real hi,
real args1,
real args2
) {
vector[1] draw;
array[1] int attempts;
draw[1] = normal_rng(args1, args2);
attempts[1] = 1;
while((conditioning_outside(draw[1], lo, hi) == 1)) {
if((attempts[1] >= 100000)) {
reject("truncated: rejection sampler exceeded 100000 draws");
}
draw[1] = normal_rng(args1, args2);
attempts[1] = (attempts[1] + 1);
}
return draw[1];
}
int conditioning_outside(
real x,
real lo,
real hi
) {
array[1] int rv;
rv[1] = 0;
if((x < lo)) {
rv[1] = 1;
} else {
if((x > hi)) {
rv[1] = 1;
}
}
return rv[1];
}
}
data {
real y;
real lo;
real hi;
}
transformed data {
}
parameters {
real mu;
}
transformed parameters {
}
model {
mu ~ normal(0.0, 1.0);
y ~ conditioning_normal(lo, hi, mu, 1.0);
}
generated quantities {
real y_likelihood = conditioning_normal_lpdfs(y, lo, hi, mu, 1.0);
real y_gen = conditioning_normal_rng(lo, hi, mu, 1.0);
}The density rejects invalid bounds, assigns negative infinity outside them, and subtracts a stable two-sided normalizing term inside. The predictive companion uses bounded rejection sampling and fails explicitly after 100,000 attempts.
Censored observations
Censoring records a threshold atom when a latent draw falls outside the bounds. The density therefore uses tail mass at an endpoint and ordinary density in the interior; prediction clamps an unconstrained draw.
Qualitatively:
y ~ censored(normal, mu, sigma; lower=lo, upper=hi)using StanBlocks
censored_model = @slic (; y=0.3, lo=-1.0, hi=1.0) begin
mu ~ normal(0.0, 1.0)
y ~ censored(normal, mu, 1.0; lower=lo, upper=hi)
endfunctions {
real clamping_normal_lpdf(
real y,
real lo,
real hi,
real args1,
real args2
) {
array[1] real rv;
if((lo >= hi)) {
reject("censored: lower bound must be less than upper bound");
}
rv[1] = negative_infinity();
if((y == lo)) {
rv[1] = normal_lcdf_stable(lo, args1, args2);
} else {
if((y == hi)) {
rv[1] = normal_lccdf_stable(hi, args1, args2);
} else {
if((y > lo)) {
if((y < hi)) {
rv[1] = normal_lpdf(y | args1, args2);
}
}
}
}
return rv[1];
}
real normal_lcdf_stable(
real x,
real loc,
real scale
) {
return (log(erfc(((-(x - loc)) / (scale * sqrt(2.0))))) - log(2.0));
}
real normal_lccdf_stable(
real x,
real loc,
real scale
) {
return (log(erfc(((x - loc) / (scale * sqrt(2.0))))) - log(2.0));
}
real clamping_normal_lpdfs(
real y,
real lo,
real hi,
real args1,
real args2
) {
return clamping_normal_lpdf(y | lo, hi, args1, args2);
}
real clamping_normal_rng(
real lo,
real hi,
real args1,
real args2
) {
vector[1] draw;
draw[1] = normal_rng(args1, args2);
if((draw[1] < lo)) {
draw[1] = lo;
} else {
if((draw[1] > hi)) {
draw[1] = hi;
}
}
return draw[1];
}
}
data {
real y;
real lo;
real hi;
}
transformed data {
}
parameters {
real mu;
}
transformed parameters {
}
model {
mu ~ normal(0.0, 1.0);
y ~ clamping_normal(lo, hi, mu, 1.0);
}
generated quantities {
real y_likelihood = clamping_normal_lpdfs(y, lo, hi, mu, 1.0);
real y_gen = clamping_normal_rng(lo, hi, mu, 1.0);
}The emitted family exposes both stable normal tail helpers. Lower and upper threshold observations contribute the matching tail mass; interior observations use normal_lpdf; generated draws are clamped to the observed support.
Interval-censored observations
An interval-censored value contributes evidence that a latent draw fell in (lo, hi]. The stored scalar is only a carrier for the observation statement; the interval endpoints determine the likelihood.
Qualitatively:
y ~ interval_censored(normal, lo, hi, mu, sigma)using StanBlocks
interval_model = @slic (; y=0.4, lo=-1.0, hi=1.0) begin
mu ~ normal(0.0, 1.0)
y ~ interval_censored(normal, lo, hi, mu, 1.0)
endfunctions {
real interval_evidence_impl_normal_lpdf(
real y,
real lo,
real hi,
real args1,
real args2
) {
array[1] real rv;
if((lo >= hi)) {
reject("interval_censored: lower bound must be less than upper bound");
}
rv[1] = log_diff_exp(normal_lcdf_stable(hi, args1, args2), normal_lcdf_stable(lo, args1, args2));
return rv[1];
}
real normal_lcdf_stable(
real x,
real loc,
real scale
) {
return (log(erfc(((-(x - loc)) / (scale * sqrt(2.0))))) - log(2.0));
}
real interval_evidence_impl_normal_lpdfs(
real y,
real lo,
real hi,
real args1,
real args2
) {
return interval_evidence_impl_normal_lpdf(y | lo, hi, args1, args2);
}
real interval_evidence_impl_normal_rng(
real lo,
real hi,
real args1,
real args2
) {
return normal_rng(args1, args2);
}
}
data {
real y;
real lo;
real hi;
}
transformed data {
}
parameters {
real mu;
}
transformed parameters {
}
model {
mu ~ normal(0.0, 1.0);
y ~ interval_evidence_impl_normal(lo, hi, mu, 1.0);
}
generated quantities {
real y_likelihood = interval_evidence_impl_normal_lpdfs(y, lo, hi, mu, 1.0);
real y_gen = interval_evidence_impl_normal_rng(lo, hi, mu, 1.0);
}The density is the stable log difference between the endpoint CDFs. Prediction returns the uncoarsened latent draw; consumers can apply their own interval reporting convention rather than losing information in the compiler.
Fused GLMs
Each fused likelihood remains a native Stan GLM in the model block. StanBlocks generates pointwise-density and predictive-RNG companions around it.
Normal identity-link GLM
X = [1.0 0.0; 1.0 1.0; 1.0 2.0]
y = [0.1, -0.2, 0.4]
normal_glm_model = @slic (; X, y, k=size(X, 2)) begin
alpha ~ normal(0, 1)
beta ~ normal(0, 1; n=k)
sigma ~ exponential(1)
y ~ normal_id_glm(X, alpha, beta, sigma)
endfunctions {
vector normal_id_glm_lpdfs(
vector y,
matrix X,
real alpha,
vector beta,
real sigma
) {
int n = dims(y)[1];
if (dims(X)[1] != n) reject("normal_id_glm_lpdfs: dim mismatch — `X` dim 1 (= ", dims(X)[1], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `X` dim 1 (= ", dims(X)[1], ").");
vector[n] rv;
for(i in 1:n) {
rv[i] = normal_lpdf(y[i] | (alpha + (X[i, :] * beta)), sigma);
}
return rv;
}
vector normal_id_glm_vector_rng(
int anontok__1,
matrix X,
real alpha,
vector beta,
real sigma
) {
int m = anontok__1;
if (dims(X)[1] != m) reject("normal_id_glm_rng: dim mismatch — `X` dim 1 (= ", dims(X)[1], ") does not match `m` (= ", m, "), inferred from `anontok__1` dim 1. `m` sizes: `anontok__1` dim 1 (= ", anontok__1, "), `X` dim 1 (= ", dims(X)[1], ").");
return normal_id_glm_rng(X, alpha, beta, sigma);
}
vector normal_id_glm_rng(
matrix X,
real alpha,
vector beta,
real sigma
) {
int m = dims(X)[1];
return to_vector(normal_rng((rep_vector(alpha, m) + (X * beta)), sigma));
}
}
data {
int k;
int y_n;
vector[y_n] y;
int X_m;
int X_n;
matrix[X_m, X_n] X;
}
transformed data {
}
parameters {
real alpha;
vector[k] beta;
real<lower=0.0> sigma;
}
transformed parameters {
}
model {
alpha ~ normal(0, 1);
beta ~ normal(0, 1);
sigma ~ exponential(1);
y ~ normal_id_glm(X, alpha, beta, sigma);
}
generated quantities {
vector[X_m] y_likelihood = normal_id_glm_lpdfs(y, X, alpha, beta, sigma);
vector[X_m] y_gen = normal_id_glm_vector_rng(y_n, X, alpha, beta, sigma);
}Bernoulli logit-link GLM
X = [1.0 0.0; 1.0 1.0; 1.0 2.0]
y = [0, 1, 1]
bernoulli_glm_model = @slic (; X, y, k=size(X, 2)) begin
alpha ~ normal(0, 1)
beta ~ normal(0, 1; n=k)
y ~ bernoulli_logit_glm(X, alpha, beta)
endfunctions {
vector bernoulli_logit_glm_lpmfs(
array[] int y,
matrix X,
real alpha,
vector beta
) {
int n = dims(y)[1];
if (dims(X)[1] != n) reject("bernoulli_logit_glm_lpmfs: dim mismatch — `X` dim 1 (= ", dims(X)[1], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `X` dim 1 (= ", dims(X)[1], ").");
vector[n] rv;
for(i in 1:n) {
rv[i] = bernoulli_logit_lpmf(y[i] | (alpha + (X[i, :] * beta)));
}
return rv;
}
array[] int bernoulli_logit_glm_int_rng(
int anontok__1,
matrix X,
real alpha,
vector beta
) {
int m = anontok__1;
if (dims(X)[1] != m) reject("bernoulli_logit_glm_rng: dim mismatch — `X` dim 1 (= ", dims(X)[1], ") does not match `m` (= ", m, "), inferred from `anontok__1` dim 1. `m` sizes: `anontok__1` dim 1 (= ", anontok__1, "), `X` dim 1 (= ", dims(X)[1], ").");
return bernoulli_logit_glm_rng(X, alpha, beta);
}
array[] int bernoulli_logit_glm_rng(
matrix X,
real alpha,
vector beta
) {
int m = dims(X)[1];
return bernoulli_logit_glm_rng(X, rep_vector(alpha, m), beta);
}
}
data {
int k;
int y_n;
array[y_n] int y;
int X_m;
int X_n;
matrix[X_m, X_n] X;
}
transformed data {
}
parameters {
real alpha;
vector[k] beta;
}
transformed parameters {
}
model {
alpha ~ normal(0, 1);
beta ~ normal(0, 1);
y ~ bernoulli_logit_glm(X, alpha, beta);
}
generated quantities {
vector[X_m] y_likelihood = bernoulli_logit_glm_lpmfs(y, X, alpha, beta);
array[X_m] int y_gen = bernoulli_logit_glm_int_rng(y_n, X, alpha, beta);
}Poisson log-link GLM
X = [1.0 0.0; 1.0 1.0; 1.0 2.0]
y = [0, 1, 2]
poisson_glm_model = @slic (; X, y, k=size(X, 2)) begin
alpha ~ normal(0, 1)
beta ~ normal(0, 1; n=k)
y ~ poisson_log_glm(X, alpha, beta)
endfunctions {
vector poisson_log_glm_lpmfs(
array[] int y,
matrix X,
real alpha,
vector beta
) {
int n = dims(y)[1];
if (dims(X)[1] != n) reject("poisson_log_glm_lpmfs: dim mismatch — `X` dim 1 (= ", dims(X)[1], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `X` dim 1 (= ", dims(X)[1], ").");
vector[n] rv;
for(i in 1:n) {
rv[i] = poisson_log_lpmf(y[i] | (alpha + (X[i, :] * beta)));
}
return rv;
}
array[] int poisson_log_glm_int_rng(
int anontok__1,
matrix X,
real alpha,
vector beta
) {
int m = anontok__1;
if (dims(X)[1] != m) reject("poisson_log_glm_rng: dim mismatch — `X` dim 1 (= ", dims(X)[1], ") does not match `m` (= ", m, "), inferred from `anontok__1` dim 1. `m` sizes: `anontok__1` dim 1 (= ", anontok__1, "), `X` dim 1 (= ", dims(X)[1], ").");
return poisson_log_glm_rng(X, alpha, beta);
}
array[] int poisson_log_glm_rng(
matrix X,
real alpha,
vector beta
) {
int m = dims(X)[1];
return poisson_log_rng((rep_vector(alpha, m) + (X * beta)));
}
}
data {
int k;
int y_n;
array[y_n] int y;
int X_m;
int X_n;
matrix[X_m, X_n] X;
}
transformed data {
}
parameters {
real alpha;
vector[k] beta;
}
transformed parameters {
}
model {
alpha ~ normal(0, 1);
beta ~ normal(0, 1);
y ~ poisson_log_glm(X, alpha, beta);
}
generated quantities {
vector[X_m] y_likelihood = poisson_log_glm_lpmfs(y, X, alpha, beta);
array[X_m] int y_gen = poisson_log_glm_int_rng(y_n, X, alpha, beta);
}Negative-binomial log-link GLM
X = [1.0 0.0; 1.0 1.0; 1.0 2.0]
y = [0, 1, 2]
negbin_glm_model = @slic (; X, y, k=size(X, 2)) begin
alpha ~ normal(0, 1)
beta ~ normal(0, 1; n=k)
phi ~ exponential(1)
y ~ neg_binomial_2_log_glm(X, alpha, beta, phi)
endfunctions {
vector neg_binomial_2_log_glm_lpmfs(
array[] int y,
matrix X,
real alpha,
vector beta,
real phi
) {
int n = dims(y)[1];
if (dims(X)[1] != n) reject("neg_binomial_2_log_glm_lpmfs: dim mismatch — `X` dim 1 (= ", dims(X)[1], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `X` dim 1 (= ", dims(X)[1], ").");
vector[n] rv;
for(i in 1:n) {
rv[i] = neg_binomial_2_lpmf(y[i] | exp((alpha + (X[i, :] * beta))), phi);
}
return rv;
}
array[] int neg_binomial_2_log_glm_int_rng(
int anontok__1,
matrix X,
real alpha,
vector beta,
real phi
) {
int m = anontok__1;
if (dims(X)[1] != m) reject("neg_binomial_2_log_glm_rng: dim mismatch — `X` dim 1 (= ", dims(X)[1], ") does not match `m` (= ", m, "), inferred from `anontok__1` dim 1. `m` sizes: `anontok__1` dim 1 (= ", anontok__1, "), `X` dim 1 (= ", dims(X)[1], ").");
return neg_binomial_2_log_glm_rng(X, alpha, beta, phi);
}
array[] int neg_binomial_2_log_glm_rng(
matrix X,
real alpha,
vector beta,
real phi
) {
int m = dims(X)[1];
return neg_binomial_2_log_rng((rep_vector(alpha, m) + (X * beta)), phi);
}
}
data {
int k;
int y_n;
array[y_n] int y;
int X_m;
int X_n;
matrix[X_m, X_n] X;
}
transformed data {
}
parameters {
real alpha;
vector[k] beta;
real<lower=0.0> phi;
}
transformed parameters {
}
model {
alpha ~ normal(0, 1);
beta ~ normal(0, 1);
phi ~ exponential(1);
y ~ neg_binomial_2_log_glm(X, alpha, beta, phi);
}
generated quantities {
vector[X_m] y_likelihood = neg_binomial_2_log_glm_lpmfs(y, X, alpha, beta, phi);
array[X_m] int y_gen = neg_binomial_2_log_glm_int_rng(y_n, X, alpha, beta, phi);
}The model block retains the fused Stan primitive. Where Stan has no matching fused RNG, StanBlocks expands only the predictive draw to the base-family RNG at alpha + X * beta; the likelihood remains fused.
Expansion, inlining, assertions, and return-type queries
Caller macros expand before tracing
macro center(x)
:($x - mean($x))
end
x = [-1.0, 0.0, 1.0]
y = [0.1, -0.2, 0.4]
macro_model = @slic (; x, y) begin
alpha ~ normal(0, 1)
beta ~ normal(0, 1)
y ~ normal(alpha + beta * @center(x), 1)
endfunctions {
vector normal_lpdfs(
vector obs,
vector loc,
int scale
) {
return jbroadcasted_normal_lpdfs(obs, loc, scale);
}
vector jbroadcasted_normal_lpdfs(
vector x1,
vector x2,
int 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), x3);
}
return rv;
}
real normal_lpdfs(
real args1,
real args2,
int 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,
int b
) {
int n = anontok__1;
return to_vector(normal_rng(a, b));
}
}
data {
int y_n;
vector[y_n] y;
int x_n;
vector[x_n] x;
}
transformed data {
}
parameters {
real alpha;
real beta;
}
transformed parameters {
}
model {
alpha ~ normal(0, 1);
beta ~ normal(0, 1);
y ~ normal((alpha + (beta * (x - mean(x)))), 1);
}
generated quantities {
vector[y_n] y_likelihood = normal_lpdfs(y, (alpha + (beta * (x - mean(x)))), 1);
vector[y_n] y_gen = normal_vector_rng(y_n, (alpha + (beta * (x - mean(x)))), 1);
}@views, @., @inbounds, and user-defined Julia macros use the same route. They expand in the caller module; StanBlocks traces the expanded syntax.
Inline helpers and caller-scope mutation
@deffun @inline scale(x::vector[n], s::real)::vector[n] = x * s
@deffun set_first!(buf::vector[n])::vector[n] = begin
buf[1] = 42.0
buf
end
@deffun mutate_scaled(x::vector[n])::vector[n] = begin
buf::vector[n] = scale(x, 2.0)
set_first!(buf)
end
inline_model = @slic (; x=[1.0, 2.0], y=0.0) begin
changed = mutate_scaled(x)
y ~ normal(sum(changed), 1)
endfunctions {
vector mutate_scaled(
vector x
) {
int n = dims(x)[1];
vector[n] buf = (x * 2.0);
buf[1] = 42.0;
return buf;
}
real normal_lpdfs(
real args1,
real args2,
int args3
) {
return normal_lpdf(args1 | args2, args3);
}
}
data {
int x_n;
vector[x_n] x;
real y;
}
transformed data {
vector[x_n] changed = mutate_scaled(x);
}
parameters {
}
transformed parameters {
}
model {
y ~ normal(sum(changed), 1);
}
generated quantities {
real y_likelihood = normal_lpdfs(y, sum(changed), 1);
real y_gen = normal_rng(sum(changed), 1);
}Calls expand at the call site rather than producing a Stan functions entry. Locals receive hygienic per-call names. A trailing ! is the Julia-convention spelling for the same inline route and makes caller-buffer mutation expressible. At a compile_slic_bundle boundary, a macro-free UDF metadata entry with markers=(:stanonly, :inline) lowers to this same path.
Runtime assertions
@deffun safe_log(x::real)::real = begin
@stan_assert x > 0 "safe_log: x must be positive"
log(x)
end
assertion_model = @slic (; x=1.0, y=0.0) begin
logged = safe_log(x)
y ~ normal(logged, 1)
endfunctions {
real safe_log(
real x
) {
if(!((x > 0))) {
reject("safe_log: x must be positive");
}
return log(x);
}
real normal_lpdfs(
real args1,
real args2,
int args3
) {
return normal_lpdf(args1 | args2, args3);
}
}
data {
real x;
real y;
}
transformed data {
real logged = safe_log(x);
}
parameters {
}
transformed parameters {
}
model {
y ~ normal(logged, 1);
}
generated quantities {
real y_likelihood = normal_lpdfs(y, logged, 1);
real y_gen = normal_rng(logged, 1);
}For compile_slic_bundle, the macro-free equivalent is an assertions record such as ((; condition=:(x > 0), message="safe_log: x must be positive"),) on the UDF metadata entry. The compiler prepends the validated record to the bodyful definition using the same @stan_assert lowering.
Transpile-time return type queries
@deffun element_type(x::real)::real = x
@deffun @stanonly copy_vec(x::vector[n]) = begin
out::return_type_of(element_type, x[1])[n]
for i in 1:n
out[i] = x[i]
end
out
end
return_type_model = @slic (; x=[1.0, 2.0], y=[0.0, 0.0]) begin
copied = copy_vec(x)
y ~ normal(copied, 1)
endfunctions {
vector copy_vec(
vector x
) {
int n = dims(x)[1];
vector[n] out;
for(i in 1:n) {
out[i] = x[i];
}
return out;
}
vector normal_lpdfs(
vector obs,
vector loc,
int scale
) {
return jbroadcasted_normal_lpdfs(obs, loc, scale);
}
vector jbroadcasted_normal_lpdfs(
vector x1,
vector x2,
int 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), x3);
}
return rv;
}
real normal_lpdfs(
real args1,
real args2,
int 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,
int b
) {
int n = anontok__1;
return to_vector(normal_rng(a, b));
}
}
data {
int x_n;
vector[x_n] x;
int y_n;
vector[y_n] y;
}
transformed data {
vector[x_n] copied = copy_vec(x);
}
parameters {
}
transformed parameters {
}
model {
y ~ normal(copied, 1);
}
generated quantities {
vector[y_n] y_likelihood = normal_lpdfs(y, copied, 1);
vector[y_n] y_gen = normal_vector_rng(y_n, copied, 1);
}return_type_of exposes the same registered inference table used by the transpiler. Computed typeof(f(x[1]))[n] annotations provide a related higher-order form: a real-valued f produces a vector, while an integer-valued f produces array[] int.
Opt-in Julia and Stan emission
@deffun definitions are Stan-only by default. Eligible deterministic, bodyful definitions annotated with @juliacompat also install one Julia method. That supports ordinary unit tests of shared deterministic helpers:
@deffun @juliacompat affine(x::real, a::real = 2.0)::real = a * x + 1.0
affine(3.0) == 7.0
dual_model = @slic (; y=0.0) begin
y ~ normal(affine(3.0), 1)
endfunctions {
real normal_lpdfs(
real args1,
real args2,
int args3
) {
return normal_lpdf(args1 | args2, args3);
}
real affine(real x, real a) {
return ((a * x) + 1.0);
}
}
data {
real y;
}
transformed data {
}
parameters {
}
transformed parameters {
}
model {
y ~ normal(affine(3.0, 2.0), 1);
}
generated quantities {
real y_likelihood = normal_lpdfs(y, affine(3.0, 2.0), 1);
real y_gen = normal_rng(affine(3.0, 2.0), 1);
}Probability, RNG, ODE, and parallel builtins are outside the bounded Julia target. Their own _lpdf/_rng-family definitions remain Stan-only even when annotated. @stanonly can document an intentionally Stan-only helper or opt one member out of a surrounding @juliacompat group.
Data binding and mock shapes
@slic captures a model body before it needs the final dataset. Bind minimal mock values to establish the types and shapes, then rebind real data later:
base = @slic (; y = [1], x = [0.0]) begin
alpha ~ normal(0, 1)
y ~ bernoulli_logit(alpha + x)
end
posterior = base(; y = [1, 0, 1], x = [-1.0, 0.0, 1.0])functions {
vector bernoulli_logit_lpmfs(
array[] int obs,
vector args1
) {
return jbroadcasted_bernoulli_logit_lpmfs(obs, args1);
}
vector jbroadcasted_bernoulli_logit_lpmfs(
array[] int x1,
vector x2
) {
int n = dims(x1)[1];
vector[n] rv;
for(i in 1:n) {
rv[i] = bernoulli_logit_lpmfs(broadcasted_getindex(x1, i), broadcasted_getindex(x2, i));
}
return rv;
}
real bernoulli_logit_lpmfs(
int args1,
real args2
) {
return bernoulli_logit_lpmf(args1 | args2);
}
int broadcasted_getindex(array[] int x, int i) {
return x[i];
}
real broadcasted_getindex(vector x, int i) {
return x[i];
}
array[] int bernoulli_logit_int_rng(
int anontok__1,
vector p
) {
int n = anontok__1;
return bernoulli_logit_rng(p);
}
}
data {
int y_n;
array[y_n] int y;
int x_n;
vector[x_n] x;
}
transformed data {
}
parameters {
real alpha;
}
transformed parameters {
}
model {
alpha ~ normal(0, 1);
y ~ bernoulli_logit((alpha + x));
}
generated quantities {
vector[y_n] y_likelihood = bernoulli_logit_lpmfs(y, (alpha + x));
array[y_n] int y_gen = bernoulli_logit_int_rng(y_n, (alpha + x));
}| Julia value | Inferred Stan carrier |
|---|---|
1 / 1.0 | int / real |
[1, 2] | array[2] int |
[1.0, 2.0] | vector[2] |
Matrix{Float64} | matrix[m,n] (transposed during Stan data preparation) |
Vector{Vector{<:Real}} | RaggedVector carrier, not a dense 2-D array |
Calling a traced StanModel with new data reuses the trace and replaces only the values. Changing a structure that affects tracing—such as whether an outcome contains missing entries—requires a new trace.
Generated observation outputs
For an ordinary top-level dense observation, qualitatively:
y ~ normal(mu, sigma)The complete minimal model makes the generated helper definitions and both derived outputs visible:
dense_output_model = @slic (; y=[0.1, -0.2, 0.4]) begin
mu ~ normal(0, 1)
sigma ~ exponential(1)
y ~ normal(mu, sigma)
endfunctions {
vector normal_lpdfs(
vector obs,
real loc,
real scale
) {
return jbroadcasted_normal_lpdfs(obs, loc, scale);
}
vector jbroadcasted_normal_lpdfs(
vector x1,
real x2,
real x3
) {
int n = dims(x1)[1];
vector[n] rv;
for(i in 1:n) {
rv[i] = normal_lpdfs(broadcasted_getindex(x1, i), x2, x3);
}
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,
real a,
real b
) {
int n = anontok__1;
return to_vector(normal_rng(rep_vector(a, n), b));
}
}
data {
int y_n;
vector[y_n] y;
}
transformed data {
}
parameters {
real mu;
real<lower=0.0> sigma;
}
transformed parameters {
}
model {
mu ~ normal(0, 1);
sigma ~ exponential(1);
y ~ normal(mu, sigma);
}
generated quantities {
vector[y_n] y_likelihood = normal_lpdfs(y, mu, sigma);
vector[y_n] y_gen = normal_vector_rng(y_n, mu, sigma);
}The current shape contract is:
| Observation form | Predictive output | Log-likelihood output |
|---|---|---|
| Dense, top-level | observation-shaped draw | elementwise |
Dense, inside plate | filled per cell | not synthesized |
| Ragged, top-level | flat draw + descriptor segments | one aggregate per group |
Ragged compiler-owned slice inside plate | flat draw + segments | one aggregate per group |
Consumers should use descriptor fields generative and source, not parse the _gen or _likelihood suffixes.
Missing continuous outcomes
Author source:
y = Union{Missing,Float64}[1.0, missing, 3.0]
m = @slic (; y) begin
mu ~ normal(0, 2)
sigma ~ exponential(1)
y ~ normal(mu, sigma)
endfunctions {
vector normal_vector_rng(
int anontok__1,
real a,
real b
) {
int n = anontok__1;
return to_vector(normal_rng(rep_vector(a, n), b));
}
vector normal_lpdfs(
vector obs,
real loc,
real scale
) {
return jbroadcasted_normal_lpdfs(obs, loc, scale);
}
vector jbroadcasted_normal_lpdfs(
vector x1,
real x2,
real x3
) {
int n = dims(x1)[1];
vector[n] rv;
for(i in 1:n) {
rv[i] = normal_lpdfs(broadcasted_getindex(x1, i), x2, x3);
}
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 merge_missing(
vector y_obs,
vector y_mis,
array[] int ii_obs,
array[] int ii_mis
) {
int n_obs = dims(y_obs)[1];
int n_mis = dims(y_mis)[1];
if (dims(ii_obs)[1] != n_obs) reject("merge_missing: dim mismatch — `ii_obs` dim 1 (= ", dims(ii_obs)[1], ") does not match `n_obs` (= ", n_obs, "), inferred from `y_obs` dim 1. `n_obs` sizes: `y_obs` dim 1 (= ", dims(y_obs)[1], "), `ii_obs` dim 1 (= ", dims(ii_obs)[1], ").");
if (dims(ii_mis)[1] != n_mis) reject("merge_missing: dim mismatch — `ii_mis` dim 1 (= ", dims(ii_mis)[1], ") does not match `n_mis` (= ", n_mis, "), inferred from `y_mis` dim 1. `n_mis` sizes: `y_mis` dim 1 (= ", dims(y_mis)[1], "), `ii_mis` dim 1 (= ", dims(ii_mis)[1], ").");
vector[(n_obs + n_mis)] y;
for(i in 1:n_obs) {
y[ii_obs[i]] = y_obs[i];
}
for(i in 1:n_mis) {
y[ii_mis[i]] = y_mis[i];
}
return y;
}
}
data {
int y_ii_mis_n;
int y_obs_n;
vector[y_obs_n] y_obs;
int y_ii_obs_n;
array[y_ii_obs_n] int y_ii_obs;
array[y_ii_mis_n] int y_ii_mis;
}
transformed data {
}
parameters {
real mu;
real<lower=0.0> sigma;
}
transformed parameters {
}
model {
mu ~ normal(0, 2);
sigma ~ exponential(1);
y_obs ~ normal(mu, sigma);
}
generated quantities {
vector[y_ii_mis_n] y_mis = normal_vector_rng(y_ii_mis_n, mu, sigma);
vector[y_obs_n] y_obs_likelihood = normal_lpdfs(y_obs, mu, sigma);
vector[y_obs_n] y_obs_gen = normal_vector_rng(y_obs_n, mu, sigma);
vector[(y_ii_obs_n + y_ii_mis_n)] y = merge_missing(y_obs, y_mis, y_ii_obs, y_ii_mis);
}The usual imputation does not enlarge the HMC parameter vector: missing values are posterior-predictive draws unless the completed outcome feeds another live likelihood. Missing predictors, discrete missing outcomes, and inherently joint outcomes such as multi_normal require an explicit model.
Deterministic programs with @deffun
@slic is a flat declaration language. @deffun is where deterministic control flow, mutation, and local allocation live. Type and shape annotations on UDF arguments, returns, and ordinary assigned locals are optional; inference specialises them from the call site and RHS. Annotate only for dispatch, named dimension locals, fresh uninitialised storage, or an otherwise ambiguous result.
Value iteration and accumulation
StanBlocks source:
@deffun centered_sum(x::vector[n])::real = begin
xbar = mean(x)
acc = 0.0
for xi in x
acc += xi - xbar
end
acc
end
centered_model = @slic (; x=[1.0, 2.0, 4.0], y=0.0) begin
centered = centered_sum(x)
y ~ normal(centered, 1)
endfunctions {
real centered_sum(
vector x
) {
real xbar = mean(x);
real acc = 0.0;
for(value_index__vi_1 in 1:num_elements(x)) {
real xi = x[value_index__vi_1];
acc += (xi - xbar);
}
return acc;
}
real normal_lpdfs(
real args1,
real args2,
int args3
) {
return normal_lpdf(args1 | args2, args3);
}
}
data {
int x_n;
vector[x_n] x;
real y;
}
transformed data {
real centered = centered_sum(x);
}
parameters {
}
transformed parameters {
}
model {
y ~ normal(centered, 1);
}
generated quantities {
real y_likelihood = normal_lpdfs(y, centered, 1);
real y_gen = normal_rng(centered, 1);
}Supported iteration forms
@deffun squares(x::vector[n]) = [xi * xi for xi in x]
@deffun weighted_sum(x::vector[n])::real = begin
acc = 0.0
for (i, xi) in enumerate(x)
acc += i * xi
end
acc
end
@deffun products(a::vector[n], b::vector[n]) =
[ai * bi for (ai, bi) in zip(a, b)]
@deffun outer(x::vector[n], y::vector[m]) =
[x[i] * y[j] for i in 1:n, j in 1:m]
iteration_model = @slic (
; x=[1.0, 2.0], a=[2.0, 3.0], b=[4.0, 5.0], y=0.0,
) begin
sq = squares(x)
ws = weighted_sum(x)
ps = products(a, b)
op = outer(x, b)
y ~ normal(sum(sq) + ws + sum(ps) + sum(to_vector(op)), 1)
endfunctions {
vector squares(
vector x
) {
vector[(1 + (num_elements(x) - 1))] comprehension_result__lc_4;
for(value_index__vi_3 in 1:num_elements(x)) {
comprehension_result__lc_4[value_index__vi_3] = (x[value_index__vi_3] * x[value_index__vi_3]);
}
return comprehension_result__lc_4;
}
real weighted_sum(
vector x
) {
real acc = 0.0;
for(i in 1:num_elements(x)) {
real xi = x[i];
acc += (i * xi);
}
return acc;
}
vector products(
vector a,
vector b
) {
int n = dims(a)[1];
if (dims(b)[1] != n) reject("products: dim mismatch — `b` dim 1 (= ", dims(b)[1], ") does not match `n` (= ", n, "), inferred from `a` dim 1. `n` sizes: `a` dim 1 (= ", dims(a)[1], "), `b` dim 1 (= ", dims(b)[1], ").");
vector[(1 + (min(num_elements(a), num_elements(b)) - 1))] comprehension_result__lc_4;
for(value_index__vi_3 in 1:min(num_elements(a), num_elements(b))) {
comprehension_result__lc_4[value_index__vi_3] = (a[value_index__vi_3] * b[value_index__vi_3]);
}
return comprehension_result__lc_4;
}
matrix outer(
vector x,
vector y
) {
int n = dims(x)[1];
int m = dims(y)[1];
matrix[(1 + (n - 1)), (1 + (m - 1))] comprehension_result__lc_2;
for(i in 1:n) {
for(j in 1:m) {
comprehension_result__lc_2[i, j] = (x[i] * y[j]);
}
}
return comprehension_result__lc_2;
}
real normal_lpdfs(
real args1,
real args2,
int args3
) {
return normal_lpdf(args1 | args2, args3);
}
}
data {
int x_n;
vector[x_n] x;
int a_n;
vector[a_n] a;
int b_n;
vector[b_n] b;
real y;
}
transformed data {
vector[(1 + (num_elements(x) - 1))] sq = squares(x);
real ws = weighted_sum(x);
vector[(1 + (min(num_elements(a), num_elements(b)) - 1))] ps = products(a, b);
matrix[(1 + (x_n - 1)), (1 + (b_n - 1))] op = outer(x, b);
}
parameters {
}
transformed parameters {
}
model {
y ~ normal((sum(sq) + ws + sum(ps) + sum(to_vector(op))), 1);
}
generated quantities {
real y_likelihood = normal_lpdfs(y, (sum(sq) + ws + sum(ps) + sum(to_vector(op))), 1);
real y_gen = normal_rng((sum(sq) + ws + sum(ps) + sum(to_vector(op))), 1);
}Emission patterns:
| Source form | Stan lowering |
|---|---|
for xi in x | index loop plus xi = x[i] |
enumerate(x) | 1-based index and element bindings |
zip(a,b,...) | loop to the shortest container |
| one-axis real comprehension | allocated vector plus fill loop |
| two-axis real comprehension | allocated matrix plus nested fill loops |
| integer-valued comprehension | array[] int, preserving index usability |
Nested if/else, for, while, =, +=, -=, *=, .= and indexed assignment work in @deffun. elseif chains, filtered/stepped generators, flattened ragged comprehensions, and 3-D+ comprehensions reject explicitly.
Post-hoc variants and cross-validation taint
Base.merge handles structural variants. Data kwargs rebind values. A lower-level cross-validation marker can additionally taint a held-out input:
person = [1, 2, 1, 2]
model = @slic (; person, y=[0.1, -0.2, 0.3, 0.0]) begin
n_person = maximum(person)
alpha ~ normal(0, 1; n=n_person)
y ~ normal(alpha[person], 1)
end
held_out = model(;
person=StanBlocks.stan.maybecv(:person, [1, 2, 1, 2]),
)functions {
vector normal_vector_rng(
int anontok__1,
int a,
int b
) {
int n = anontok__1;
return to_vector(normal_rng(rep_vector(a, n), b));
}
vector normal_lpdfs(
vector obs,
vector loc,
int scale
) {
return jbroadcasted_normal_lpdfs(obs, loc, scale);
}
vector jbroadcasted_normal_lpdfs(
vector x1,
vector x2,
int 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), x3);
}
return rv;
}
real normal_lpdfs(
real args1,
real args2,
int 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,
int b
) {
int n = anontok__1;
return to_vector(normal_rng(a, b));
}
}
data {
int person_n;
array[person_n] int person;
int y_n;
vector[y_n] y;
}
transformed data {
int n_person = max(person);
}
parameters {
}
transformed parameters {
}
model {
}
generated quantities {
vector[n_person] alpha = normal_vector_rng(n_person, 0, 1);
vector[y_n] y_likelihood = normal_lpdfs(y, alpha[person], 1);
vector[y_n] y_gen = normal_vector_rng(y_n, alpha[person], 1);
}Activity propagates from the marked input. Likelihood terms reached by the mark move out of the fit; affected latent variables are redrawn in generated quantities; unaffected population parameters stay fitted. The descriptor reports the resulting held_out state and removes :fit if no live likelihood remains.
This is compiler machinery rather than a complete user-facing CV workflow. BRM and other consumers decide how units, folds, posterior draws, and result objects are presented.
Executable model descriptors
This is Julia descriptor API, not a transpilation example:
d = stan_descriptor(m; name=:demo)For the core regression MWE, the current descriptor reports:
inputs:
y_n derived=true
y observed=true
x_n derived=true
x observed=false
outputs:
alpha parameter / posterior
beta parameter / posterior
sigma parameter / posterior
y_likelihood generated_quantity / pointwise_loglik / source=y
y_gen generated_quantity / draw / source=y
operations:
transpile · instantiate · fit · predict · pointwise_loglikThe descriptor also publishes the ordered inventory of emitted Stan function definitions, exact signatures, source spans, and dependency links. Consumers can select an exact definition and its transitive closure without parsing Stan text.
The execution calls below consume that descriptor and likewise emit no new Stan program:
problem = stan_execute(d, :fit)
pred = stan_execute(d, :predict; problem, draws=theta, seed=123)
ll = stan_execute(d, :pointwise_loglik; problem, draws=theta, seed=123)Operations are derived and fail closed. For example, :fit exists only if at least one parameter and one live likelihood reach the model block.
Scientific computing surface
ODE solvers
@deffun pk_rhs(t::real, y::vector[ny], ke::real)::vector[ny] = begin
dy::vector[ny]
dy[1] = -ke * y[1]
dy
end
ts = [0.5, 1.0]
conc = [0.8, 0.6]
ode_model = @slic (; ts, conc) begin
ke ~ lognormal(0, 0.5)
sigma ~ exponential(1)
pred = to_vector(
ode_rk45(pk_rhs, [1.0], 0.0, to_array_1d(ts), ke)[:, 1]
)
conc ~ lognormal(log(pred), sigma)
endfunctions {
vector pk_rhs(
real t,
vector y,
real ke
) {
int ny = dims(y)[1];
vector[ny] dy;
dy[1] = ((-ke) * y[1]);
return dy;
}
vector lognormal_lpdfs(
vector obs,
vector loc,
real scale
) {
return jbroadcasted_lognormal_lpdfs(obs, loc, scale);
}
vector jbroadcasted_lognormal_lpdfs(
vector x1,
vector x2,
real x3
) {
int n = dims(x1)[1];
vector[n] rv;
for(i in 1:n) {
rv[i] = lognormal_lpdfs(broadcasted_getindex(x1, i), broadcasted_getindex(x2, i), x3);
}
return rv;
}
real lognormal_lpdfs(
real args1,
real args2,
real args3
) {
return lognormal_lpdf(args1 | args2, args3);
}
real broadcasted_getindex(vector x, int i) {
return x[i];
}
vector lognormal_vector_rng(
int anontok__1,
vector a,
real b
) {
int n = anontok__1;
return to_vector(lognormal_rng(a, b));
}
}
data {
int ts_n;
vector[ts_n] ts;
int conc_n;
vector[conc_n] conc;
}
transformed data {
}
parameters {
real<lower=0.0> ke;
real<lower=0.0> sigma;
}
transformed parameters {
vector[ts_n] pred = to_vector(ode_rk45(pk_rhs, [1.0]', 0.0, to_array_1d(ts), ke)[:, 1]);
}
model {
ke ~ lognormal(0, 0.5);
sigma ~ exponential(1);
conc ~ lognormal(log(pred), sigma);
}
generated quantities {
vector[conc_n] conc_likelihood = lognormal_lpdfs(conc, log(pred), sigma);
vector[conc_n] conc_gen = lognormal_vector_rng(conc_n, log(pred), sigma);
}ode_rk45 is the default and ode_bdf the stiff escalation. Solver results are Stan arrays of state vectors, so extract one state with [:,k] and convert it back to a vector for a vector observation.
Torsten
The builtin signature layer includes the Torsten analytical pmx_solve_onecpt/pmx_solve_twocpt family. With a Torsten-enabled BridgeStan build, these calls have been exercised end to end through compilation and a gradient-correct log density. They remain an environment-dependent extension of the ordinary BridgeStan route, not a bundled StanBlocks compiler.
Other registered scientific primitives
reduce_sum,reduce_sum_static, andsimple_reduce_sum;Gaussian-process covariance functions, including multidimensional inputs;
Stan ODE variants and tolerance forms;
matrix decompositions, solves, eigenvalues, and covariance helpers;
fused GLMs and hundreds of scalar/vector/matrix probability signatures.
The compute path
@slic model + @deffun helpers
│
├── stan_model / stan_code ──→ inspectable Stan source
│
├── stanc_check / compiles ──→ syntax + Stan semantic checks
│
└── stan_instantiate ────────→ BridgeStan-backed StanProblem
│
└─ LogDensityProblems value + gradienttranspiles(model) checks only tracing/code generation. A meaningful compiler change still needs stanc and, where semantics matter, a BridgeStan log-density/gradient comparison. The generated Stan may also be compiled and sampled through the usual CmdStan workflow.
Honest current boundaries
| Requested construct | Current answer |
|---|---|
Ordinary for/if in @slic | No: put deterministic control flow in @deffun |
| A loop that introduces independent parameters | Use compiler-owned plate |
A scan where cell i consumes cell i-1 | No: write a deterministic recurrence in @deffun |
| Matrix-valued plate cell | No: shared matrix outside, scalar/vector cell result |
| General 3-D+ Julia container | No |
| Filtered/stepped/ragged comprehension | No; write an explicit supported loop |
| Missing continuous outcome vector | Automatic |
| Missing predictor or discrete outcome | Explicit model required |
| Ragged continuous observation | Yes, with groupwise log likelihood |
| Ragged integer observation/predictive draw | No carrier yet |
| Arbitrary existing Julia function | No auto-transpilation; register through @deffun/builtins |
Direct target += in a UDF | No; express a distribution through _lpdf/@lpxf |
| Full Julia runtime parity for every Stan builtin | No; deterministic @deffun subset only |
What changed since the submitted abstract
| Abstract wording | StanCon 2026 status |
|---|---|
| closures were “potentially planned” | shipped, including captured ODE parameters/data |
| keyword/default arguments were “potentially planned” | shipped: required/optional kwargs and positional defaults |
| Julia-style metaprogramming was “potentially planned” | caller macros expand before tracing; inline helpers and assertions shipped |
| composable submodels | expanded with named positional/typed submodel functions and first-class merged callees |
| custom struct-like types / named tuples | tuple/named-tuple results and nominal compiler carriers support structured lowering |
| automated predictive/log-likelihood output | broadened to descriptors, missing outcomes, plates, and ragged-group semantics |
| PKPD motivation | now includes numerical ODE and Torsten signatures, with SbPMX as a higher-level PKPD consumer |
Primary references
Maria I. Gorinova, Andrew D. Gordon & Charles Sutton, Probabilistic Programming with Densities in SlicStan: Efficient, Flexible and Deterministic, POPL 2019.
SlicStan public repository, including its explicit description as a blockless Stan-like language and its research-code caveat.
StanBlocks.jl, current authoring support, and API reference.