Motorcycle data
This example follows the HSGP model sequence in Aki Vehtari's motorcycle case study. The response is head acceleration after impact. Its mean changes nonlinearly over time, and the size of the residual fluctuations changes as the impact is absorbed, so a straight line with constant noise is a poor description.
The small random arrays below are build fixtures rather than the original measurements. Every comparison evaluates the exact displayed Julia source at documentation-build time and places the complete generated Stan program beside it, with the same expandable side-by-side modal as the feature atlas.
HSGP building block (hsgp below)
A Hilbert-space Gaussian-process (HSGP) approximation replaces a dense Gaussian-process covariance matrix with a finite basis expansion. That makes a smooth latent function a weighted matrix-vector product while preserving two interpretable hyperparameters: a length scale and a marginal amplitude.
The source first rescales x to a fixed interval and constructs a sine basis matrix X with 20 columns. x_scale controls how quickly the function can vary, y_scale controls its marginal standard deviation, and unit_weight ~ std_normal(; n=n_functions) introduces one standard-normal coefficient per basis function. The deterministic scale vector turns those coefficients into the GP approximation returned by the submodel.
Several StanBlocks features are doing work here:
@slictraces the Julia linear algebra and elementwise expressions rather than requiring separate Stan declarations.The
n=n_functionssampling keyword determines the length ofunit_weight; the productX * (scale .* unit_weight)therefore has its observation-length result inferred automatically.hsgp(; x)binds the free inputxas data. The returned vector can then be embedded in another@slicmodel withlhs ~ hsgp(; x).Positive support is encoded through the chosen distributions (
uniform(0, 2)andlognormal(0, 1)), and the Stan pane shows where the resulting parameters and transformed quantities are emitted.
using StanBlocks, Markdown
x = randn(10)
obs = randn(10)
hsgp = @slic begin
"Transforms to [-2, 2]"
xi = 2 * (x - min(x)) / (max(x) - min(x)) - 1.
L = 1.5
n_functions = 20
X = sin(pi/(2L) * (xi+L) * range(1,n_functions)')/sqrt(L)
"The GP lengthscale"
x_scale ~ uniform(0, 2)
"The GP marginal standard deviation"
y_scale ~ lognormal(0, 1)
"The scales for the basis functions weights"
scale = y_scale * sqrt(sqrt(2pi) * x_scale) * exp(-0.25*(x_scale*pi/2L)^2 * range(1,n_functions)^2)
"The basis functions weights"
unit_weight ~ std_normal(;n=n_functions)
"The final GP values"
return (X * (scale .* unit_weight))
end
hsgp_posterior = hsgp(;x)// Transforms to [-2, 2]
functions {
vector std_normal_vector_rng(
int anontok__1
) {
int n = anontok__1;
return to_vector(normal_rng(rep_vector(0, n), 1));
}
}
data {
int x_n;
vector[x_n] x;
}
transformed data {
vector[x_n] xi = (((2 * (x - min(x))) / (max(x) - min(x))) - 1.0);
real L = 1.5;
int n_functions = 20;
matrix[x_n, n_functions] X = (
sin(((3.141592653589793 / (2 * L)) * (xi + L) * (linspaced_vector(n_functions, 1, n_functions)'))) /
sqrt(L)
);
}
parameters {
}
transformed parameters {
}
model {
}
generated quantities {
// The GP lengthscale
real x_scale = uniform_rng(0, 2);
// The GP marginal standard deviation
real y_scale = lognormal_rng(0, 1);
// The scales for the basis functions weights
vector[n_functions] scale = (
y_scale *
sqrt((sqrt((2 * 3.141592653589793)) * x_scale)) *
exp(
(
-0.25 *
(((x_scale * 3.141592653589793) / (2 * L)) ^ 2) *
(linspaced_vector(n_functions, 1, n_functions) ^ 2)
)
)
);
// The basis functions weights
vector[n_functions] unit_weight = std_normal_vector_rng(n_functions);
// The final GP values
vector[x_n] MODEL_RV = (X * (scale .* unit_weight));
}Homoskedastic model
The first observation model uses one HSGP for the mean acceleration,
y_intercept sets the global level, dy ~ hsgp(; x) embeds the smooth deviation, and the scalar sigma is shared by every observation. This is the homoskedastic assumption: the mean may be highly nonlinear, but residual spread is constant. Binding homo(; x, obs) supplies both free names as data; the HSGP's local parameters are inlined into the parent without exposing them as manual arguments.
homo = @slic begin
y_intercept ~ std_normal()
dy ~ hsgp(;x)
sigma ~ lognormal(-2, 1)
obs ~ normal(y_intercept + dy, sigma)
end
homo_posterior = homo(;x, obs)functions {
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 obs_n;
vector[obs_n] obs;
}
transformed data {
vector[x_n] dy_xi = (((2 * (x - min(x))) / (max(x) - min(x))) - 1.0);
real dy_L = 1.5;
int dy_n_functions = 20;
matrix[x_n, dy_n_functions] dy_X = (
sin(
(
(3.141592653589793 / (2 * dy_L)) *
(dy_xi + dy_L) *
(linspaced_vector(dy_n_functions, 1, dy_n_functions)')
)
) /
sqrt(dy_L)
);
}
parameters {
real y_intercept;
// The GP lengthscale
real<lower=0, upper=2> dy_x_scale;
// The GP marginal standard deviation
real<lower=0.0> dy_y_scale;
// The basis functions weights
vector[dy_n_functions] dy_unit_weight;
real<lower=0.0> sigma;
}
transformed parameters {
// The scales for the basis functions weights
vector[dy_n_functions] dy_scale = (
dy_y_scale *
sqrt((sqrt((2 * 3.141592653589793)) * dy_x_scale)) *
exp(
(
-0.25 *
(((dy_x_scale * 3.141592653589793) / (2 * dy_L)) ^ 2) *
(linspaced_vector(dy_n_functions, 1, dy_n_functions) ^ 2)
)
)
);
// The final GP values
vector[x_n] dy = (dy_X * (dy_scale .* dy_unit_weight));
}
model {
y_intercept ~ std_normal();
// The GP lengthscale
dy_x_scale ~ uniform(0, 2);
// The GP marginal standard deviation
dy_y_scale ~ lognormal(0, 1);
// The basis functions weights
dy_unit_weight ~ std_normal();
sigma ~ lognormal(-2, 1);
obs ~ normal((y_intercept + dy), sigma);
}
generated quantities {
vector[obs_n] obs_likelihood = normal_lpdfs(obs, (y_intercept + dy), sigma);
vector[obs_n] obs_gen = normal_vector_rng(obs_n, (y_intercept + dy), sigma);
}Heteroskedastic model
The residual amplitude visibly changes over time, so the next model gives the log standard deviation its own smooth function:
Calling hsgp twice creates two independent, hygienically renamed sets of GP hyperparameters and basis weights: one for the mean and one for the log scale. exp(log_sigma_intercept + dlog_sigma) maps the second process to positive standard deviations. Relative to the homoskedastic model, the likelihood is still normal and the mean component is unchanged; only the scalar sigma is replaced by an observation-length vector.
hetero = @slic begin
y_intercept ~ std_normal()
dy ~ hsgp(;x)
log_sigma_intercept ~ std_normal()
dlog_sigma ~ hsgp(;x)
obs ~ normal(y_intercept + dy, exp(log_sigma_intercept + dlog_sigma))
end
hetero_posterior = hetero(;x,obs)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 x_n;
vector[x_n] x;
int obs_n;
vector[obs_n] obs;
}
transformed data {
vector[x_n] dy_xi = (((2 * (x - min(x))) / (max(x) - min(x))) - 1.0);
real dy_L = 1.5;
int dy_n_functions = 20;
matrix[x_n, dy_n_functions] dy_X = (
sin(
(
(3.141592653589793 / (2 * dy_L)) *
(dy_xi + dy_L) *
(linspaced_vector(dy_n_functions, 1, dy_n_functions)')
)
) /
sqrt(dy_L)
);
vector[x_n] dlog_sigma_xi = (((2 * (x - min(x))) / (max(x) - min(x))) - 1.0);
real dlog_sigma_L = 1.5;
int dlog_sigma_n_functions = 20;
matrix[x_n, dlog_sigma_n_functions] dlog_sigma_X = (
sin(
(
(3.141592653589793 / (2 * dlog_sigma_L)) *
(dlog_sigma_xi + dlog_sigma_L) *
(linspaced_vector(dlog_sigma_n_functions, 1, dlog_sigma_n_functions)')
)
) /
sqrt(dlog_sigma_L)
);
}
parameters {
real y_intercept;
// The GP lengthscale
real<lower=0, upper=2> dy_x_scale;
// The GP marginal standard deviation
real<lower=0.0> dy_y_scale;
// The basis functions weights
vector[dy_n_functions] dy_unit_weight;
real log_sigma_intercept;
// The GP lengthscale
real<lower=0, upper=2> dlog_sigma_x_scale;
// The GP marginal standard deviation
real<lower=0.0> dlog_sigma_y_scale;
// The basis functions weights
vector[dlog_sigma_n_functions] dlog_sigma_unit_weight;
}
transformed parameters {
// The scales for the basis functions weights
vector[dy_n_functions] dy_scale = (
dy_y_scale *
sqrt((sqrt((2 * 3.141592653589793)) * dy_x_scale)) *
exp(
(
-0.25 *
(((dy_x_scale * 3.141592653589793) / (2 * dy_L)) ^ 2) *
(linspaced_vector(dy_n_functions, 1, dy_n_functions) ^ 2)
)
)
);
// The final GP values
vector[x_n] dy = (dy_X * (dy_scale .* dy_unit_weight));
// The scales for the basis functions weights
vector[dlog_sigma_n_functions] dlog_sigma_scale = (
dlog_sigma_y_scale *
sqrt((sqrt((2 * 3.141592653589793)) * dlog_sigma_x_scale)) *
exp(
(
-0.25 *
(((dlog_sigma_x_scale * 3.141592653589793) / (2 * dlog_sigma_L)) ^ 2) *
(linspaced_vector(dlog_sigma_n_functions, 1, dlog_sigma_n_functions) ^ 2)
)
)
);
// The final GP values
vector[x_n] dlog_sigma = (dlog_sigma_X * (dlog_sigma_scale .* dlog_sigma_unit_weight));
}
model {
y_intercept ~ std_normal();
// The GP lengthscale
dy_x_scale ~ uniform(0, 2);
// The GP marginal standard deviation
dy_y_scale ~ lognormal(0, 1);
// The basis functions weights
dy_unit_weight ~ std_normal();
log_sigma_intercept ~ std_normal();
// The GP lengthscale
dlog_sigma_x_scale ~ uniform(0, 2);
// The GP marginal standard deviation
dlog_sigma_y_scale ~ lognormal(0, 1);
// The basis functions weights
dlog_sigma_unit_weight ~ std_normal();
obs ~ normal((y_intercept + dy), exp((log_sigma_intercept + dlog_sigma)));
}
generated quantities {
vector[obs_n] obs_likelihood = normal_lpdfs(obs, (y_intercept + dy), exp((log_sigma_intercept + dlog_sigma)));
vector[obs_n] obs_gen = normal_vector_rng(obs_n, (y_intercept + dy), exp((log_sigma_intercept + dlog_sigma)));
}Alternative heteroskedastic model using subsubmodels
The following formulation has the same statistical structure but factors the repeated “intercept plus HSGP” pattern into another submodel. It demonstrates nested SLIC composition rather than adding a new assumption.
Submodel with submodel (intercept_hsgp below)
intercept_hsgp samples an intercept, embeds hsgp, and returns their sum. The nested call keeps x as a keyword-bound data dependency and keeps all fresh parameters local to this use. StanBlocks follows the composition through both levels when it generates the final Stan program.
intercept_hsgp = @slic begin
intercept ~ std_normal()
"Submodel uses `hsgp` as a submodel"
d ~ hsgp(;x)
return intercept + d
end
intercept_hsgp_posterior = intercept_hsgp(;x)functions {
vector std_normal_vector_rng(
int anontok__1
) {
int n = anontok__1;
return to_vector(normal_rng(rep_vector(0, n), 1));
}
}
data {
int x_n;
vector[x_n] x;
}
transformed data {
vector[x_n] d_xi = (((2 * (x - min(x))) / (max(x) - min(x))) - 1.0);
real d_L = 1.5;
int d_n_functions = 20;
matrix[x_n, d_n_functions] d_X = (
sin(
(
(3.141592653589793 / (2 * d_L)) *
(d_xi + d_L) *
(linspaced_vector(d_n_functions, 1, d_n_functions)')
)
) /
sqrt(d_L)
);
}
parameters {
}
transformed parameters {
}
model {
}
generated quantities {
real intercept = std_normal_rng();
// The GP lengthscale
real d_x_scale = uniform_rng(0, 2);
// The GP marginal standard deviation
real d_y_scale = lognormal_rng(0, 1);
// The scales for the basis functions weights
vector[d_n_functions] d_scale = (
d_y_scale *
sqrt((sqrt((2 * 3.141592653589793)) * d_x_scale)) *
exp(
(
-0.25 *
(((d_x_scale * 3.141592653589793) / (2 * d_L)) ^ 2) *
(linspaced_vector(d_n_functions, 1, d_n_functions) ^ 2)
)
)
);
// The basis functions weights
vector[d_n_functions] d_unit_weight = std_normal_vector_rng(d_n_functions);
// The final GP values
vector[x_n] d = (d_X * (d_scale .* d_unit_weight));
vector[x_n] MODEL_RV = (intercept + d);
}Final model
The outer model calls intercept_hsgp twice. The first returned vector is the mean y; the second is log_sigma, exponentiated in the likelihood. Because each submodel use is hygienic, the two processes do not accidentally share an intercept, length scale, amplitude, or basis weights. Compare its Stan pane with hetero: the model is statistically equivalent, but the Julia source expresses the reusable structure directly.
hetero2 = @slic begin
y ~ intercept_hsgp(;x)
log_sigma ~ intercept_hsgp(;x)
obs ~ normal(y, exp(log_sigma))
end
hetero2_posterior = hetero2(;x,obs)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 x_n;
vector[x_n] x;
int obs_n;
vector[obs_n] obs;
}
transformed data {
vector[x_n] y_d_xi = (((2 * (x - min(x))) / (max(x) - min(x))) - 1.0);
real y_d_L = 1.5;
int y_d_n_functions = 20;
matrix[x_n, y_d_n_functions] y_d_X = (
sin(
(
(3.141592653589793 / (2 * y_d_L)) *
(y_d_xi + y_d_L) *
(linspaced_vector(y_d_n_functions, 1, y_d_n_functions)')
)
) /
sqrt(y_d_L)
);
vector[x_n] log_sigma_d_xi = (((2 * (x - min(x))) / (max(x) - min(x))) - 1.0);
real log_sigma_d_L = 1.5;
int log_sigma_d_n_functions = 20;
matrix[x_n, log_sigma_d_n_functions] log_sigma_d_X = (
sin(
(
(3.141592653589793 / (2 * log_sigma_d_L)) *
(log_sigma_d_xi + log_sigma_d_L) *
(linspaced_vector(log_sigma_d_n_functions, 1, log_sigma_d_n_functions)')
)
) /
sqrt(log_sigma_d_L)
);
}
parameters {
real y_intercept;
// The GP lengthscale
real<lower=0, upper=2> y_d_x_scale;
// The GP marginal standard deviation
real<lower=0.0> y_d_y_scale;
// The basis functions weights
vector[y_d_n_functions] y_d_unit_weight;
real log_sigma_intercept;
// The GP lengthscale
real<lower=0, upper=2> log_sigma_d_x_scale;
// The GP marginal standard deviation
real<lower=0.0> log_sigma_d_y_scale;
// The basis functions weights
vector[log_sigma_d_n_functions] log_sigma_d_unit_weight;
}
transformed parameters {
// The scales for the basis functions weights
vector[y_d_n_functions] y_d_scale = (
y_d_y_scale *
sqrt((sqrt((2 * 3.141592653589793)) * y_d_x_scale)) *
exp(
(
-0.25 *
(((y_d_x_scale * 3.141592653589793) / (2 * y_d_L)) ^ 2) *
(linspaced_vector(y_d_n_functions, 1, y_d_n_functions) ^ 2)
)
)
);
// The final GP values
vector[x_n] y_d = (y_d_X * (y_d_scale .* y_d_unit_weight));
vector[x_n] y = (y_intercept + y_d);
// The scales for the basis functions weights
vector[log_sigma_d_n_functions] log_sigma_d_scale = (
log_sigma_d_y_scale *
sqrt((sqrt((2 * 3.141592653589793)) * log_sigma_d_x_scale)) *
exp(
(
-0.25 *
(((log_sigma_d_x_scale * 3.141592653589793) / (2 * log_sigma_d_L)) ^ 2) *
(linspaced_vector(log_sigma_d_n_functions, 1, log_sigma_d_n_functions) ^ 2)
)
)
);
// The final GP values
vector[x_n] log_sigma_d = (log_sigma_d_X * (log_sigma_d_scale .* log_sigma_d_unit_weight));
vector[x_n] log_sigma = (log_sigma_intercept + log_sigma_d);
}
model {
y_intercept ~ std_normal();
// The GP lengthscale
y_d_x_scale ~ uniform(0, 2);
// The GP marginal standard deviation
y_d_y_scale ~ lognormal(0, 1);
// The basis functions weights
y_d_unit_weight ~ std_normal();
log_sigma_intercept ~ std_normal();
// The GP lengthscale
log_sigma_d_x_scale ~ uniform(0, 2);
// The GP marginal standard deviation
log_sigma_d_y_scale ~ lognormal(0, 1);
// The basis functions weights
log_sigma_d_unit_weight ~ std_normal();
obs ~ normal(y, exp(log_sigma));
}
generated quantities {
vector[obs_n] obs_likelihood = normal_lpdfs(obs, y, exp(log_sigma));
vector[obs_n] obs_gen = normal_vector_rng(obs_n, y, exp(log_sigma));
}