Likelihoods
Complete built-in catalogue
This is the exhaustive public likelihood catalogue for the default SBBRMI Stan backend. A Distributions.jl subtype not named here is not accepted merely because it is a distribution: BRM rejects it until its Julia constructor has an explicit Stan mapping. The pure-Julia VBRMI backend is independent and does not implement every specialized family or response wrapper below.
Direct Distributions.jl families
| Outcome | Accepted constructors |
|---|---|
| Continuous | Normal, NormalCanon, Cauchy, TDist, Logistic, Gumbel, Chisq, Exponential, Gamma, Erlang, Beta, Uniform, LogNormal, Laplace, Frechet, Rayleigh, SkewNormal, Pareto, Weibull, InverseGamma, InverseGaussian, VonMises |
| Continuous, restricted parameterization | Arcsine() — the standard [0, 1] form only; SkewedExponentialPower(mu, sigma, 1, alpha) — only the literal shape 1 |
| Discrete | Bernoulli, BernoulliLogit, Binomial, BinomialLogit, BetaBinomial, Poisson, NegativeBinomial |
BRM normalizes constructor conventions where Julia and Stan differ. In particular, Exponential, Gamma, and Erlang use scale in Distributions.jl but rate in Stan; Pareto(shape, scale) is reordered to Stan's (minimum, shape) convention; NormalCanon(eta, lambda) becomes normal(eta / lambda, inv(sqrt(lambda))); NegativeBinomial(r, p) is translated to Stan's shape/inverse-scale form; TDist(nu) becomes student_t(nu, 0, 1); and Laplace becomes Stan's double_exponential. The standard Arcsine() is exactly beta(1/2, 1/2). Shifted/scaled Arcsine(a, b) constructors need a Jacobian-aware custom implementation and are rejected rather than silently treated as a standard beta likelihood.
BRM families and structured likelihoods
| Constructor | Meaning / boundary |
|---|---|
SkewDoubleExponential(mu, sigma, tau) | Stan-native asymmetric-Laplace parameterization |
LocationScale(mu, sigma, TDist(nu)) | location-scale Student-t regression |
ZeroInflatedPoisson(lambda, zi) | zero-inflated Poisson |
HurdlePoisson(lambda, p_zero) | hurdle Poisson with zero-truncated positive component |
NegativeBinomial2(mu, phi) | mean/precision negative binomial |
BetaBinomial2(n, mean, precision) | mean/precision beta-binomial |
CategoricalLogit(eta2, eta3, ...) or CategoricalLogit(@brm(...)) | reference-class categorical logit |
OrderedLogistic(eta) | legacy cumulative-logit ordinal model |
Ordinal(structure, link, eta; ...) | Cumulative() or StoppingRatio() crossed with LogitLink(), ProbitLink(), or CloglogLink() |
CircularVonMises(mu, kappa; interval=(-pi, pi)) | von Mises on a fixed principal interval |
TruncatedNormal(mu, sigma, lower, upper) | legacy censored-Normal marker; new models should use censored below |
[y1, y2, ...] ~ MvNormalCholesky([mu1, mu2, ...], L) | row-wise correlated Gaussian outcomes using a declared LKJCovarianceFactor |
MixtureModel([D, ...], weights) | finite mixture over K same-family scalar components |
Response compositions and modifiers
| Form | Exact supported surface |
|---|---|
truncated(d; lower, upper) | base family Normal, LogNormal, Exponential, Weibull, or Poisson |
censored(d; lower, upper) | the same five base families |
interval_censored(d; upper) | the same five base families; the response is the lower endpoint |
weighted(d, aweights(w)) | analytic/precision weights for Normal only |
weighted(d, fweights(w)) | frequency weights for direct mapped families in the first table |
weighted(d, weights(w)) | power-likelihood weights for direct mapped families in the first table |
mi(y) ~ d | partly-missing continuous response imputation |
Every specialized family is expected to supply the fitted density, pointwise log likelihood, and posterior-predictive RNG used by BRM's generated quantities. Truncation and censoring additionally require matching CDF/CCDF paths; ragged observations additionally require a sized RNG.
Hurdle counts with HurdlePoisson
Use HurdlePoisson(lambda, p_zero) when zeros arise from a separate hurdle and every count from the Poisson component is strictly positive:
brm-comparison
Hurdle-Poisson countsjulia
using LogExpFunctions: logit
hurdle_counts = (@brm begin
log(lambda) ~ 1 + x
logit(p_zero) ~ 1 + factor(group)
y ~ HurdlePoisson(lambda, p_zero)
end)((;
x=[0.0, 0.5, 1.0, 1.5, 2.0, 2.5],
group=[1, 1, 2, 2, 3, 3],
y=[0, 1, 0, 3, 2, 5],
))julia
BRMI:
x: data (eltype=Float64, n=6)
log(lambda) ~ 1 + x
group: data (eltype=Int64, n=6)
logit(p_zero) ~ 1 + factor(group)
y ~ HurdlePoisson(lambda, p_zero)julia
SBBRMI with data keys = [:group, :group_idx, :group_n_levels, :x, :y]
emitted @slic body:
begin
X_log_lambda = hcat(rep_vector(1.0, num_elements(x)), x)
pop_log_lambda ~ popefs(; X = X_log_lambda)
log_lambda = pop_log_lambda
lambda = exp(log_lambda)
X_logit_p_zero = hcat(rep_vector(1.0, num_elements(group)))
pop_logit_p_zero ~ popefs(; X = X_logit_p_zero)
cat_logit_p_zero_group ~ _sb_cat(; x = group_idx, n_levels = group_n_levels)
logit_p_zero = pop_logit_p_zero + cat_logit_p_zero_group
p_zero = inv_logit(logit_p_zero)
y ~ hurdle_poisson(lambda, p_zero)
endstan
functions {
matrix hcat(
vector x,
vector y
) {
int n = dims(x)[1];
if (dims(y)[1] != n) reject("hcat: dim mismatch — `y` dim 1 (= ", dims(y)[1], ") does not match `n` (= ", n, "), inferred from `x` dim 1. `n` sizes: `x` dim 1 (= ", dims(x)[1], "), `y` dim 1 (= ", dims(y)[1], ").");
return append_col(x, y);
}
matrix hcat(vector x) {
int n = dims(x)[1];
return to_matrix(x, n, 1);
}
real hurdle_poisson_lpmf(
array[] int y,
vector lambda,
vector p_zero
) {
int n = dims(y)[1];
if (dims(lambda)[1] != n) reject("hurdle_poisson_lpmf: dim mismatch — `lambda` dim 1 (= ", dims(lambda)[1], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `lambda` dim 1 (= ", dims(lambda)[1], "), `p_zero` dim 1 (= ", dims(p_zero)[1], ").");
if (dims(p_zero)[1] != n) reject("hurdle_poisson_lpmf: dim mismatch — `p_zero` dim 1 (= ", dims(p_zero)[1], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `lambda` dim 1 (= ", dims(lambda)[1], "), `p_zero` dim 1 (= ", dims(p_zero)[1], ").");
real rv = 0.0;
for(i in 1:n) {
rv += hurdle_poisson_lpmf(y[i] | lambda[i], p_zero[i]);
}
return rv;
}
real hurdle_poisson_lpmf(
int y,
real lambda,
real p_zero
) {
if((y == 0)) {
return log(p_zero);
} else {
return ((log1m(p_zero) + poisson_lpmf(y | lambda)) - poisson_lccdf(0 | lambda));
}
}
vector hurdle_poisson_lpmfs(
array[] int y,
vector lambda,
vector p_zero
) {
int n = dims(y)[1];
if (dims(lambda)[1] != n) reject("hurdle_poisson_lpmfs: dim mismatch — `lambda` dim 1 (= ", dims(lambda)[1], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `lambda` dim 1 (= ", dims(lambda)[1], "), `p_zero` dim 1 (= ", dims(p_zero)[1], ").");
if (dims(p_zero)[1] != n) reject("hurdle_poisson_lpmfs: dim mismatch — `p_zero` dim 1 (= ", dims(p_zero)[1], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `lambda` dim 1 (= ", dims(lambda)[1], "), `p_zero` dim 1 (= ", dims(p_zero)[1], ").");
vector[n] rv;
for(i in 1:n) {
rv[i] = hurdle_poisson_lpmf(y[i] | lambda[i], p_zero[i]);
}
return rv;
}
array[] int hurdle_poisson_int_rng(
int anontok__1,
vector lambda,
vector p_zero
) {
int n = anontok__1;
if (dims(lambda)[1] != n) reject("hurdle_poisson_rng: dim mismatch — `lambda` dim 1 (= ", dims(lambda)[1], ") does not match `n` (= ", n, "), inferred from `anontok__1` dim 1. `n` sizes: `anontok__1` dim 1 (= ", anontok__1, "), `lambda` dim 1 (= ", dims(lambda)[1], "), `p_zero` dim 1 (= ", dims(p_zero)[1], ").");
if (dims(p_zero)[1] != n) reject("hurdle_poisson_rng: dim mismatch — `p_zero` dim 1 (= ", dims(p_zero)[1], ") does not match `n` (= ", n, "), inferred from `anontok__1` dim 1. `n` sizes: `anontok__1` dim 1 (= ", anontok__1, "), `lambda` dim 1 (= ", dims(lambda)[1], "), `p_zero` dim 1 (= ", dims(p_zero)[1], ").");
array[n] int rv;
for(i in 1:n) {
rv[i] = hurdle_poisson_rng(lambda[i], p_zero[i]);
}
return rv;
}
int hurdle_poisson_rng(
real lambda,
real p_zero
) {
if((bernoulli_rng(p_zero) == 1)) {
return 0;
} else {
return hurdle_poisson_positive_rng(lambda);
}
}
int hurdle_poisson_positive_rng(
real lambda
) {
array[1] int draw;
draw[1] = poisson_rng(lambda);
while((draw[1] == 0)) {
draw[1] = poisson_rng(lambda);
}
return draw[1];
}
}
data {
int x_n;
vector[x_n] x;
int group_n;
array[group_n] int group;
int group_n_levels;
int group_idx_n;
array[group_idx_n] int group_idx;
int y_n;
array[y_n] int y;
}
transformed data {
matrix[x_n, 2] X_log_lambda = hcat(rep_vector(1.0, num_elements(x)), x);
int pop_log_lambda_n_covariates = 2;
matrix[num_elements(group), 1] X_logit_p_zero = hcat(rep_vector(1.0, num_elements(group)));
int pop_logit_p_zero_n_covariates = 1;
}
parameters {
vector[pop_log_lambda_n_covariates] pop_log_lambda_beta_pop;
vector[pop_logit_p_zero_n_covariates] pop_logit_p_zero_beta_pop;
vector[(group_n_levels - 1)] cat_logit_p_zero_group_beta;
}
transformed parameters {
vector[x_n] pop_log_lambda = (X_log_lambda * pop_log_lambda_beta_pop);
vector[x_n] log_lambda = pop_log_lambda;
vector[x_n] lambda = exp(log_lambda);
vector[num_elements(group)] pop_logit_p_zero = (X_logit_p_zero * pop_logit_p_zero_beta_pop);
vector[group_idx_n] cat_logit_p_zero_group = append_row(0.0, cat_logit_p_zero_group_beta)[group_idx];
vector[num_elements(group)] logit_p_zero = (pop_logit_p_zero + cat_logit_p_zero_group);
vector[num_elements(group)] p_zero = inv_logit(logit_p_zero);
}
model {
pop_log_lambda_beta_pop ~ std_normal();
pop_logit_p_zero_beta_pop ~ std_normal();
cat_logit_p_zero_group_beta ~ std_normal();
y ~ hurdle_poisson(lambda, p_zero);
}
generated quantities {
vector[num_elements(group)] y_likelihood = hurdle_poisson_lpmfs(y, lambda, p_zero);
array[num_elements(group)] int y_gen = hurdle_poisson_int_rng(y_n, lambda, p_zero);
}julia
#= line 0 =# Turing.@model(function brm_model(y, X_lambda, X_p_zero)
beta_pop ~ Distributions.product_distribution(Base.fill(Distributions.Normal(), 2))
eta_lambda = X_lambda * beta_pop
lambda = Base.exp.(eta_lambda)
beta_pop_p_zero ~ Distributions.product_distribution(Base.fill(Distributions.Normal(), 3))
eta_p_zero = X_p_zero * beta_pop_p_zero
p_zero = LogExpFunctions.logistic.(eta_p_zero)
begin
for i = Base.eachindex(y)
y[i] ~ BayesianRegressionModels.HurdlePoisson(lambda[i], p_zero[i])
end
end
(; lambda = lambda, p_zero = p_zero, response = y)
end)Here p_zero is exactly
This differs from ZeroInflatedPoisson(lambda, zi): zero inflation mixes a structural-zero component with an ordinary Poisson component, so both mixture components can produce zeros. HurdlePoisson is also an executable Distributions.jl distribution with matching params, logpdf, and rand semantics outside a formula. It requires finite lambda > 0 and 0 <= p_zero <= 1; BRM supplies no implicit link or prior.
Wald regression with InverseGaussian
Use Distributions.jl's InverseGaussian(mu, lambda) for positive continuous outcomes whose variance grows with the cube of the mean — the Wald GLM with a log link:
brm-comparison
Log-link Wald costsjulia
wald_costs = (@brm begin
log(lambda) ~ 1
eta ~ 1 + x
y ~ InverseGaussian(exp(eta), lambda)
end)((;
x=[-1.0, -0.5, 0.0, 0.5, 1.0],
y=[1.2, 0.8, 1.1, 2.0, 1.6],
))julia
BRMI:
log(lambda) ~ 1
x: data (eltype=Float64, n=5)
eta ~ 1 + x
y ~ InverseGaussian(exp(eta), lambda)julia
SBBRMI with data keys = [:x, :y]
emitted @slic body:
begin
X_log_lambda = hcat(rep_vector(1.0, num_elements(y)))
pop_log_lambda ~ popefs(; X = X_log_lambda)
log_lambda = pop_log_lambda
lambda = exp(log_lambda)
X_eta = hcat(rep_vector(1.0, num_elements(x)), x)
pop_eta ~ popefs(; X = X_eta)
eta = pop_eta
y ~ brm_inverse_gaussian((exp)(eta), lambda)
endstan
functions {
matrix hcat(vector x) {
int n = dims(x)[1];
return to_matrix(x, n, 1);
}
matrix hcat(
vector x,
vector y
) {
int n = dims(x)[1];
if (dims(y)[1] != n) reject("hcat: dim mismatch — `y` dim 1 (= ", dims(y)[1], ") does not match `n` (= ", n, "), inferred from `x` dim 1. `n` sizes: `x` dim 1 (= ", dims(x)[1], "), `y` dim 1 (= ", dims(y)[1], ").");
return append_col(x, y);
}
real brm_inverse_gaussian_lpdf(
vector y,
vector mu,
vector lambda
) {
int n = dims(y)[1];
if (dims(mu)[1] != n) reject("brm_inverse_gaussian_lpdf: dim mismatch — `mu` dim 1 (= ", dims(mu)[1], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `mu` dim 1 (= ", dims(mu)[1], "), `lambda` dim 1 (= ", dims(lambda)[1], ").");
if (dims(lambda)[1] != n) reject("brm_inverse_gaussian_lpdf: dim mismatch — `lambda` dim 1 (= ", dims(lambda)[1], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `mu` dim 1 (= ", dims(mu)[1], "), `lambda` dim 1 (= ", dims(lambda)[1], ").");
real rv = 0.0;
for(i in 1:n) {
rv += brm_inverse_gaussian_lpdf(y[i] | mu[i], lambda[i]);
}
return rv;
}
real brm_inverse_gaussian_lpdf(
real y,
real mu,
real lambda
) {
if((mu <= 0.0)) {
return negative_infinity();
} else {
if((lambda <= 0.0)) {
return negative_infinity();
} else {
if((y <= 0.0)) {
return negative_infinity();
} else {
return (
(
(log(lambda) - (1.8378770664093456 + (3.0 * log(y)))) -
((lambda * (y - mu) * (y - mu)) / (mu * mu * y))
) /
2.0
);
}
}
}
}
vector brm_inverse_gaussian_lpdfs(
vector y,
vector mu,
vector lambda
) {
int n = dims(y)[1];
if (dims(mu)[1] != n) reject("brm_inverse_gaussian_lpdfs: dim mismatch — `mu` dim 1 (= ", dims(mu)[1], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `mu` dim 1 (= ", dims(mu)[1], "), `lambda` dim 1 (= ", dims(lambda)[1], ").");
if (dims(lambda)[1] != n) reject("brm_inverse_gaussian_lpdfs: dim mismatch — `lambda` dim 1 (= ", dims(lambda)[1], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `mu` dim 1 (= ", dims(mu)[1], "), `lambda` dim 1 (= ", dims(lambda)[1], ").");
vector[n] rv;
for(i in 1:n) {
rv[i] = brm_inverse_gaussian_lpdf(y[i] | mu[i], lambda[i]);
}
return rv;
}
vector brm_inverse_gaussian_vector_rng(
int anontok__1,
vector mu,
vector lambda
) {
int n = anontok__1;
if (dims(mu)[1] != n) reject("brm_inverse_gaussian_rng: dim mismatch — `mu` dim 1 (= ", dims(mu)[1], ") does not match `n` (= ", n, "), inferred from `anontok__1` dim 1. `n` sizes: `anontok__1` dim 1 (= ", anontok__1, "), `mu` dim 1 (= ", dims(mu)[1], "), `lambda` dim 1 (= ", dims(lambda)[1], ").");
if (dims(lambda)[1] != n) reject("brm_inverse_gaussian_rng: dim mismatch — `lambda` dim 1 (= ", dims(lambda)[1], ") does not match `n` (= ", n, "), inferred from `anontok__1` dim 1. `n` sizes: `anontok__1` dim 1 (= ", anontok__1, "), `mu` dim 1 (= ", dims(mu)[1], "), `lambda` dim 1 (= ", dims(lambda)[1], ").");
vector[n] rv;
for(i in 1:n) {
rv[i] = brm_inverse_gaussian_rng(mu[i], lambda[i]);
}
return rv;
}
real brm_inverse_gaussian_rng(
real mu,
real lambda
) {
if((mu <= 0.0)) {
reject("brm_inverse_gaussian_rng: mu must be strictly positive");
return 0.0;
} else {
if((lambda <= 0.0)) {
reject("brm_inverse_gaussian_rng: lambda must be strictly positive");
return 0.0;
} else {
real z = normal_rng(0.0, 1.0);
real v = (z * z);
real w = (mu * v);
real x = (mu + ((mu / (2.0 * lambda)) * (w - sqrt((w * ((4.0 * lambda) + w))))));
real u = uniform_rng(0.0, 1.0);
if((u < (mu / (mu + x)))) {
return x;
} else {
return ((mu * mu) / x);
}
}
}
}
}
data {
int y_n;
vector[y_n] y;
int x_n;
vector[x_n] x;
}
transformed data {
matrix[num_elements(y), 1] X_log_lambda = hcat(rep_vector(1.0, num_elements(y)));
int pop_log_lambda_n_covariates = 1;
matrix[x_n, 2] X_eta = hcat(rep_vector(1.0, num_elements(x)), x);
int pop_eta_n_covariates = 2;
}
parameters {
vector[pop_log_lambda_n_covariates] pop_log_lambda_beta_pop;
vector[pop_eta_n_covariates] pop_eta_beta_pop;
}
transformed parameters {
vector[num_elements(y)] pop_log_lambda = (X_log_lambda * pop_log_lambda_beta_pop);
vector[num_elements(y)] log_lambda = pop_log_lambda;
vector[num_elements(y)] lambda = exp(log_lambda);
vector[x_n] pop_eta = (X_eta * pop_eta_beta_pop);
vector[x_n] eta = pop_eta;
}
model {
pop_log_lambda_beta_pop ~ std_normal();
pop_eta_beta_pop ~ std_normal();
y ~ brm_inverse_gaussian(exp(eta), lambda);
}
generated quantities {
vector[num_elements(y)] y_likelihood = brm_inverse_gaussian_lpdfs(y, exp(eta), lambda);
vector[num_elements(y)] y_gen = brm_inverse_gaussian_vector_rng(y_n, exp(eta), lambda);
}julia
#= line 0 =# Turing.@model(function brm_model(y, X_lambda, X_eta)
beta_pop ~ Distributions.product_distribution(Base.fill(Distributions.Normal(), 1))
eta_lambda = X_lambda * beta_pop
lambda = Base.exp.(eta_lambda)
beta_pop_eta ~ Distributions.product_distribution(Base.fill(Distributions.Normal(), 2))
eta_eta = X_eta * beta_pop_eta
eta = eta_eta
begin
for i = Base.eachindex(y)
y[i] ~ Distributions.InverseGaussian(Base.exp(eta[i]), lambda[i])
end
end
(; lambda = lambda, eta = eta, response = y)
end)This preserves the constructor order (mu, lambda) — already Stan's order, so no parameterization translation applies — the shorthand InverseGaussian(mu) == InverseGaussian(mu, 1), and the strict y > 0, mu > 0, lambda > 0 domain. Density, pointwise log likelihood, and predictive RNG all use the same closed form; BRM supplies no implicit link or prior.
Correlated Gaussian outcomes
Use one ordered vector likelihood when several measurements from the same row have experimental residual covariance that should be estimated rather than treated as independent:
brm-comparison
Estimated experimental covariancejulia
correlated_measurements = (@brm begin
L_res ~ LKJCovarianceFactor(
2; scale_prior=Exponential(1), shape=2,
)
shared_log_rate ~ Normal(0, 1)
concentration_mu = exp(-exp(shared_log_rate) * time)
response_mu = 1.0 - concentration_mu
[concentration, response] ~ MvNormalCholesky(
[concentration_mu, response_mu], L_res)
end)((;
time=[0.0, 1.0, 2.0, 4.0],
concentration=[1.0, 0.72, 0.51, 0.27],
response=[0.03, 0.18, 0.43, 0.79],
))julia
BRMI:
L_res ~ LKJCovarianceFactor(2; scale_prior=Exponential(1), shape=2)
shared_log_rate ~ Normal(0, 1)
time: data (eltype=Float64, n=4)
:concentration_mu = exp((-(exp(shared_log_rate)) * time))
:response_mu = -(1.0, concentration_mu)
[concentration, response] ~ MvNormalCholesky(NamedColumn{Symbol}[concentration_mu, response_mu], L_res)julia
SBBRMI with data keys = [:L_res_n, :brm_joint_concentration__response_n, :brm_joint_concentration__response_observed, :time]
emitted @slic body:
begin
L_res_scales ~ exponential(1.0; n = L_res_n)
L_res_L_corr::cholesky_factor_corr[L_res_n] ~ lkj_corr_cholesky(2.0)
L_res = diag_pre_multiply(L_res_scales, L_res_L_corr)
shared_log_rate ~ normal(0, 1)
concentration_mu = (exp)((-)((exp)(shared_log_rate)) .* time)
response_mu = (-)(1.0, concentration_mu)
brm_joint_concentration__response_means ~ plate(brm_joint_concentration__response_observed, concentration_mu, response_mu; outer = (brm_joint_concentration__response_n,)) do brm_joint_concentration__response_observed_cell, brm_joint_concentration__response_mean_cell_1, brm_joint_concentration__response_mean_cell_2
brm_joint_concentration__response_mean_vector = 0.0 .* brm_joint_concentration__response_observed_cell + [brm_joint_concentration__response_mean_cell_1, brm_joint_concentration__response_mean_cell_2]
brm_joint_concentration__response_mean_vector
end
brm_joint_concentration__response_observed ~ multi_normal_cholesky(brm_joint_concentration__response_means, L_res)
endstan
functions {
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)]);
}
}
int num_elements_RaggedVector(tuple(vector, array[] int) rv) {
return size(rv.2);
}
real multi_normal_cholesky_lpdfs(
vector args1,
vector args2,
matrix args3
) {
return multi_normal_cholesky_lpdf(args1 | args2, args3);
}
vector multi_normal_cholesky_vector_rng(
int anontok__1,
vector loc,
matrix scale
) {
int n = anontok__1;
if (dims(loc)[1] != n) reject("multi_normal_cholesky_rng: dim mismatch — `loc` dim 1 (= ", dims(loc)[1], ") does not match `n` (= ", n, "), inferred from `anontok__1` dim 1. `n` sizes: `anontok__1` dim 1 (= ", anontok__1, "), `loc` dim 1 (= ", dims(loc)[1], ").");
return multi_normal_cholesky_rng(loc, scale);
}
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)];
}
}
data {
int L_res_n;
int time_n;
vector[time_n] time;
int brm_joint_concentration__response_n;
int brm_joint_concentration__response_observed_ends_n;
int brm_joint_concentration__response_observed_mem_n;
tuple(
vector[brm_joint_concentration__response_observed_mem_n],
array[brm_joint_concentration__response_observed_ends_n] int
) brm_joint_concentration__response_observed;
}
transformed data {
array[brm_joint_concentration__response_n] int brm_joint_concentration__response_means__pl_len_1;
array[brm_joint_concentration__response_n] int brm_joint_concentration__response_means_brm_joint_concentration__response_mean_vector__pl_len_1;
for(plate_i__pl_1 in 1:brm_joint_concentration__response_n) {
brm_joint_concentration__response_means__pl_len_1[plate_i__pl_1] = (
1 +
(
ragged_end(brm_joint_concentration__response_observed.2, plate_i__pl_1) -
ragged_start(brm_joint_concentration__response_observed.2, plate_i__pl_1)
)
);
brm_joint_concentration__response_means_brm_joint_concentration__response_mean_vector__pl_len_1[
plate_i__pl_1
] = (
1 +
(
ragged_end(brm_joint_concentration__response_observed.2, plate_i__pl_1) -
ragged_start(brm_joint_concentration__response_observed.2, plate_i__pl_1)
)
);
}
array[brm_joint_concentration__response_n] int brm_joint_concentration__response_means__pl_end_1 = cumulative_sum(brm_joint_concentration__response_means__pl_len_1);
array[brm_joint_concentration__response_n] int brm_joint_concentration__response_means_brm_joint_concentration__response_mean_vector__pl_end_1 = cumulative_sum(
brm_joint_concentration__response_means_brm_joint_concentration__response_mean_vector__pl_len_1
);
}
parameters {
vector<lower=0.0>[L_res_n] L_res_scales;
cholesky_factor_corr[L_res_n] L_res_L_corr;
real shared_log_rate;
}
transformed parameters {
matrix[L_res_n, L_res_n] L_res = diag_pre_multiply(L_res_scales, L_res_L_corr);
vector[time_n] concentration_mu = exp(((-exp(shared_log_rate)) .* time));
vector[time_n] response_mu = (1.0 - concentration_mu);
vector[sum(brm_joint_concentration__response_means__pl_len_1)] brm_joint_concentration__response_means__pl_mem_1;
vector[
sum(brm_joint_concentration__response_means_brm_joint_concentration__response_mean_vector__pl_len_1)
] brm_joint_concentration__response_means_brm_joint_concentration__response_mean_vector__pl_mem_1;
for(plate_i__pl_1 in 1:brm_joint_concentration__response_n) {
brm_joint_concentration__response_means_brm_joint_concentration__response_mean_vector__pl_mem_1[
ragged_start(
brm_joint_concentration__response_means_brm_joint_concentration__response_mean_vector__pl_end_1,
plate_i__pl_1
):ragged_end(
brm_joint_concentration__response_means_brm_joint_concentration__response_mean_vector__pl_end_1,
plate_i__pl_1
)
] = (
(
0.0 .*
brm_joint_concentration__response_observed.1[
ragged_start(brm_joint_concentration__response_observed.2, plate_i__pl_1):ragged_end(brm_joint_concentration__response_observed.2, plate_i__pl_1)
]
) +
[concentration_mu[plate_i__pl_1], response_mu[plate_i__pl_1]]'
);
brm_joint_concentration__response_means__pl_mem_1[
ragged_start(brm_joint_concentration__response_means__pl_end_1, plate_i__pl_1):ragged_end(brm_joint_concentration__response_means__pl_end_1, plate_i__pl_1)
] = brm_joint_concentration__response_means_brm_joint_concentration__response_mean_vector__pl_mem_1[
ragged_start(
brm_joint_concentration__response_means_brm_joint_concentration__response_mean_vector__pl_end_1,
plate_i__pl_1
):ragged_end(
brm_joint_concentration__response_means_brm_joint_concentration__response_mean_vector__pl_end_1,
plate_i__pl_1
)
];
}
}
model {
L_res_scales ~ exponential(1.0);
L_res_L_corr ~ lkj_corr_cholesky(2.0);
shared_log_rate ~ normal(0, 1);
for(g__ro_2 in 1:num_elements_RaggedVector(brm_joint_concentration__response_observed)) {
getindex_RaggedVector(brm_joint_concentration__response_observed, g__ro_2) ~ multi_normal_cholesky(
brm_joint_concentration__response_means__pl_mem_1[
ragged_start(brm_joint_concentration__response_means__pl_end_1, g__ro_2):ragged_end(brm_joint_concentration__response_means__pl_end_1, g__ro_2)
],
L_res
);
}
}
generated quantities {
vector[num_elements(brm_joint_concentration__response_observed.1)] brm_joint_concentration__response_observed_gen;
vector[num_elements_RaggedVector(brm_joint_concentration__response_observed)] brm_joint_concentration__response_observed_likelihood;
for(g__rq_2 in 1:num_elements_RaggedVector(brm_joint_concentration__response_observed)) {
brm_joint_concentration__response_observed_gen[
ragged_start(brm_joint_concentration__response_observed.2, g__rq_2):ragged_end(brm_joint_concentration__response_observed.2, g__rq_2)
] = multi_normal_cholesky_vector_rng(
(1 + (ragged_end_RaggedVector(brm_joint_concentration__response_observed, g__rq_2) - ragged_start_RaggedVector(brm_joint_concentration__response_observed, g__rq_2))),
brm_joint_concentration__response_means__pl_mem_1[
ragged_start(brm_joint_concentration__response_means__pl_end_1, g__rq_2):ragged_end(brm_joint_concentration__response_means__pl_end_1, g__rq_2)
],
L_res
);
brm_joint_concentration__response_observed_likelihood[g__rq_2] = multi_normal_cholesky_lpdf(getindex_RaggedVector(brm_joint_concentration__response_observed, g__rq_2) |
brm_joint_concentration__response_means__pl_mem_1[
ragged_start(brm_joint_concentration__response_means__pl_end_1, g__rq_2):ragged_end(brm_joint_concentration__response_means__pl_end_1, g__rq_2)
],
L_res
);
}
}julia
#= line 0 =# Turing.@model(function brm_model(y, time)
L_res ~ DynamicPPL.to_submodel(BayesianRegressionModels._brm_turing_covariance_prior(2; scale_prior = Distributions.Exponential(1), shape = 2))
shared_log_rate ~ Distributions.Normal(0, 1)
concentration_mu = [Base.exp(-(Base.exp(shared_log_rate)) * time[i]) for i = Base.eachindex(y)]
response_mu = [1.0 - concentration_mu[i] for i = Base.eachindex(y)]
begin
for i = Base.eachindex(y)
y[i] ~ BayesianRegressionModelsTuringExt._brm_mvn_cholesky(Base.vect(concentration_mu[i], response_mu[i]), L_res)
end
end
(; shared_log_rate = shared_log_rate, L_res = L_res, response_mu = response_mu, concentration_mu = concentration_mu, response = y)
end)LKJCovarianceFactor(K; scale_prior=Exponential(1), shape=1) samples K positive marginal scales and an LKJ Cholesky correlation factor, then returns diag_pre_multiply(scales, L_corr). The likelihood lowers to Stan's native multi_normal_cholesky: each aligned data row contributes one joint scalar log likelihood and one ordered predictive vector. The example deliberately uses one sampled rate in both mechanistic means; shared parameters need no special syntax beyond ordinary formula-block assignments.
This first surface is deliberately strict. Every outcome and mean has the same row axis, the number of means and factor dimension must match the left-hand side, and outcome rows must be complete, finite, and nonempty. A missing value is rejected; BRM never silently drops the row or replaces the joint density with conditionally independent pieces. Correlated outcomes are currently an SBBRMI-only feature.
Finite mixtures with MixtureModel
Use Distributions.jl's MixtureModel when each observation comes from one of K latent subpopulations with its own parameters:
brm-comparison
Two-component Gaussian mixturejulia
mixture_gaussians = (@brm begin
mu1 ~ Normal(-2, 0.1)
mu2 ~ Normal(2, 0.1)
log(sigma) ~ 1
y ~ MixtureModel([
Normal(mu1, exp(log(sigma))),
Normal(mu2, exp(log(sigma))),
], [0.4, 0.6])
end)((;
y=[-2.0, -1.8, 1.9, 2.2],
))julia
BRMI:
mu1 ~ Normal(-2, 0.1)
mu2 ~ Normal(2, 0.1)
log(sigma) ~ 1
y ~ MixtureModel(ExprColumn{Type{Normal}, Tuple{NamedColumn{Symbol, ExprColumn{typeof(~), Tuple{NamedColumn{Symbol, MissingColumn}, ExprColumn{Type{Normal}, Tuple{Int64, Float64}, @NamedTuple{}}}, @NamedTuple{}}}, ExprColumn{typeof(exp), Tuple{ExprColumn{typeof(log), Tuple{NamedColumn{Symbol, ExprColumn{typeof(~), Tuple{ExprColumn{typeof(log), Tuple{NamedColumn{Symbol, MissingColumn}}, @NamedTuple{}}, Int64}, @NamedTuple{}}}}, @NamedTuple{}}}, @NamedTuple{}}}, @NamedTuple{}}[Normal(mu1, exp(log(sigma))), Normal(mu2, exp(log(sigma)))], [0.4, 0.6])julia
SBBRMI with data keys = [:y, :y_mixture_weights]
emitted @slic body:
begin
mu1 ~ normal(-2, 0.1)
mu2 ~ normal(2, 0.1)
X_log_sigma = hcat(rep_vector(1.0, num_elements(y)))
pop_log_sigma ~ popefs(; X = X_log_sigma)
log_sigma = pop_log_sigma
sigma = exp(log_sigma)
y ~ (ValueFamily(brm_mixture_c0badac3fba85005))(y_mixture_weights, brm_joint_mean_rows(mu1, num_elements(y)), brm_joint_mean_rows((exp)((log)(sigma)), num_elements(y)), brm_joint_mean_rows(mu2, num_elements(y)), brm_joint_mean_rows((exp)((log)(sigma)), num_elements(y)))
endstan
functions {
matrix hcat(vector x) {
int n = dims(x)[1];
return to_matrix(x, n, 1);
}
// value UDF brm_mixture_c0badac3fba85005_lpdf
real brm_mixture_c0badac3fba85005_lpdf(
vector y,
vector weights,
vector c1_1,
vector c1_2,
vector c2_1,
vector c2_2
) {
int n = dims(y)[1];
if (dims(c1_1)[1] != n) reject("brm_mixture_c0badac3fba85005_lpdf: dim mismatch — `c1_1` dim 1 (= ", dims(c1_1)[1], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `c1_1` dim 1 (= ", dims(c1_1)[1], "), `c1_2` dim 1 (= ", dims(c1_2)[1], "), `c2_1` dim 1 (= ", dims(c2_1)[1], "), `c2_2` dim 1 (= ", dims(c2_2)[1], ").");
if (dims(c1_2)[1] != n) reject("brm_mixture_c0badac3fba85005_lpdf: dim mismatch — `c1_2` dim 1 (= ", dims(c1_2)[1], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `c1_1` dim 1 (= ", dims(c1_1)[1], "), `c1_2` dim 1 (= ", dims(c1_2)[1], "), `c2_1` dim 1 (= ", dims(c2_1)[1], "), `c2_2` dim 1 (= ", dims(c2_2)[1], ").");
if (dims(c2_1)[1] != n) reject("brm_mixture_c0badac3fba85005_lpdf: dim mismatch — `c2_1` dim 1 (= ", dims(c2_1)[1], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `c1_1` dim 1 (= ", dims(c1_1)[1], "), `c1_2` dim 1 (= ", dims(c1_2)[1], "), `c2_1` dim 1 (= ", dims(c2_1)[1], "), `c2_2` dim 1 (= ", dims(c2_2)[1], ").");
if (dims(c2_2)[1] != n) reject("brm_mixture_c0badac3fba85005_lpdf: dim mismatch — `c2_2` dim 1 (= ", dims(c2_2)[1], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `c1_1` dim 1 (= ", dims(c1_1)[1], "), `c1_2` dim 1 (= ", dims(c1_2)[1], "), `c2_1` dim 1 (= ", dims(c2_1)[1], "), `c2_2` dim 1 (= ", dims(c2_2)[1], ").");
real rv = 0.0;
for(i in 1:n) {
rv = (rv + brm_mixture_c0badac3fba85005_lpdf(y[i] | weights, c1_1[i], c1_2[i], c2_1[i], c2_2[i]));
}
return rv;
}
// value UDF brm_mixture_c0badac3fba85005_lpdf
real brm_mixture_c0badac3fba85005_lpdf(
real y,
vector weights,
real c1_1,
real c1_2,
real c2_1,
real c2_2
) {
int K = dims(weights)[1];
vector[K] terms;
terms[1] = (log(weights[1]) + normal_lpdf(y | c1_1, c1_2));
terms[2] = (log(weights[2]) + normal_lpdf(y | c2_1, c2_2));
return log_sum_exp(terms);
}
// value UDF brm_mixture_c0badac3fba85005_lpdfs
vector brm_mixture_c0badac3fba85005_lpdfs(
vector y,
vector weights,
vector c1_1,
vector c1_2,
vector c2_1,
vector c2_2
) {
int n = dims(y)[1];
if (dims(c1_1)[1] != n) reject("brm_mixture_c0badac3fba85005_lpdfs: dim mismatch — `c1_1` dim 1 (= ", dims(c1_1)[1], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `c1_1` dim 1 (= ", dims(c1_1)[1], "), `c1_2` dim 1 (= ", dims(c1_2)[1], "), `c2_1` dim 1 (= ", dims(c2_1)[1], "), `c2_2` dim 1 (= ", dims(c2_2)[1], ").");
if (dims(c1_2)[1] != n) reject("brm_mixture_c0badac3fba85005_lpdfs: dim mismatch — `c1_2` dim 1 (= ", dims(c1_2)[1], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `c1_1` dim 1 (= ", dims(c1_1)[1], "), `c1_2` dim 1 (= ", dims(c1_2)[1], "), `c2_1` dim 1 (= ", dims(c2_1)[1], "), `c2_2` dim 1 (= ", dims(c2_2)[1], ").");
if (dims(c2_1)[1] != n) reject("brm_mixture_c0badac3fba85005_lpdfs: dim mismatch — `c2_1` dim 1 (= ", dims(c2_1)[1], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `c1_1` dim 1 (= ", dims(c1_1)[1], "), `c1_2` dim 1 (= ", dims(c1_2)[1], "), `c2_1` dim 1 (= ", dims(c2_1)[1], "), `c2_2` dim 1 (= ", dims(c2_2)[1], ").");
if (dims(c2_2)[1] != n) reject("brm_mixture_c0badac3fba85005_lpdfs: dim mismatch — `c2_2` dim 1 (= ", dims(c2_2)[1], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `c1_1` dim 1 (= ", dims(c1_1)[1], "), `c1_2` dim 1 (= ", dims(c1_2)[1], "), `c2_1` dim 1 (= ", dims(c2_1)[1], "), `c2_2` dim 1 (= ", dims(c2_2)[1], ").");
vector[n] rv;
for(i in 1:n) {
rv[i] = brm_mixture_c0badac3fba85005_lpdf(y[i] | weights, c1_1[i], c1_2[i], c2_1[i], c2_2[i]);
}
return rv;
}
// value UDF brm_mixture_c0badac3fba85005_rng
vector brm_mixture_c0badac3fba85005_vector_rng(
int anontok__1,
vector weights,
vector c1_1,
vector c1_2,
vector c2_1,
vector c2_2
) {
int n = anontok__1;
if (dims(c1_1)[1] != n) reject("brm_mixture_c0badac3fba85005_rng: dim mismatch — `c1_1` dim 1 (= ", dims(c1_1)[1], ") does not match `n` (= ", n, "), inferred from `anontok__1` dim 1. `n` sizes: `anontok__1` dim 1 (= ", anontok__1, "), `c1_1` dim 1 (= ", dims(c1_1)[1], "), `c1_2` dim 1 (= ", dims(c1_2)[1], "), `c2_1` dim 1 (= ", dims(c2_1)[1], "), `c2_2` dim 1 (= ", dims(c2_2)[1], ").");
if (dims(c1_2)[1] != n) reject("brm_mixture_c0badac3fba85005_rng: dim mismatch — `c1_2` dim 1 (= ", dims(c1_2)[1], ") does not match `n` (= ", n, "), inferred from `anontok__1` dim 1. `n` sizes: `anontok__1` dim 1 (= ", anontok__1, "), `c1_1` dim 1 (= ", dims(c1_1)[1], "), `c1_2` dim 1 (= ", dims(c1_2)[1], "), `c2_1` dim 1 (= ", dims(c2_1)[1], "), `c2_2` dim 1 (= ", dims(c2_2)[1], ").");
if (dims(c2_1)[1] != n) reject("brm_mixture_c0badac3fba85005_rng: dim mismatch — `c2_1` dim 1 (= ", dims(c2_1)[1], ") does not match `n` (= ", n, "), inferred from `anontok__1` dim 1. `n` sizes: `anontok__1` dim 1 (= ", anontok__1, "), `c1_1` dim 1 (= ", dims(c1_1)[1], "), `c1_2` dim 1 (= ", dims(c1_2)[1], "), `c2_1` dim 1 (= ", dims(c2_1)[1], "), `c2_2` dim 1 (= ", dims(c2_2)[1], ").");
if (dims(c2_2)[1] != n) reject("brm_mixture_c0badac3fba85005_rng: dim mismatch — `c2_2` dim 1 (= ", dims(c2_2)[1], ") does not match `n` (= ", n, "), inferred from `anontok__1` dim 1. `n` sizes: `anontok__1` dim 1 (= ", anontok__1, "), `c1_1` dim 1 (= ", dims(c1_1)[1], "), `c1_2` dim 1 (= ", dims(c1_2)[1], "), `c2_1` dim 1 (= ", dims(c2_1)[1], "), `c2_2` dim 1 (= ", dims(c2_2)[1], ").");
vector[n] rv;
for(i in 1:n) {
rv[i] = brm_mixture_c0badac3fba85005_rng(weights, c1_1[i], c1_2[i], c2_1[i], c2_2[i]);
}
return rv;
}
// value UDF brm_mixture_c0badac3fba85005_rng
real brm_mixture_c0badac3fba85005_rng(
vector weights,
real c1_1,
real c1_2,
real c2_1,
real c2_2
) {
int k = categorical_rng(weights);
if((k == 1)) {
return normal_rng(c1_1, c1_2);
} else {
return normal_rng(c2_1, c2_2);
}
}
vector brm_joint_mean_rows(real value, int rows) {
return rep_vector(value, rows);
}
vector brm_joint_mean_rows(
vector value,
int rows
) {
int n = dims(value)[1];
if(!((n == rows))) {
reject("assertion failed: n == rows");
}
return value;
}
}
data {
int y_n;
vector[y_n] y;
int y_mixture_weights_n;
vector[y_mixture_weights_n] y_mixture_weights;
}
transformed data {
matrix[num_elements(y), 1] X_log_sigma = hcat(rep_vector(1.0, num_elements(y)));
int pop_log_sigma_n_covariates = 1;
}
parameters {
real mu1;
real mu2;
vector[pop_log_sigma_n_covariates] pop_log_sigma_beta_pop;
}
transformed parameters {
vector[num_elements(y)] pop_log_sigma = (X_log_sigma * pop_log_sigma_beta_pop);
vector[num_elements(y)] log_sigma = pop_log_sigma;
vector[num_elements(y)] sigma = exp(log_sigma);
}
model {
mu1 ~ normal(-2, 0.1);
mu2 ~ normal(2, 0.1);
pop_log_sigma_beta_pop ~ std_normal();
y ~ brm_mixture_c0badac3fba85005(
y_mixture_weights,
brm_joint_mean_rows(mu1, num_elements(y)),
brm_joint_mean_rows(exp(log(sigma)), num_elements(y)),
brm_joint_mean_rows(mu2, num_elements(y)),
brm_joint_mean_rows(exp(log(sigma)), num_elements(y))
);
}
generated quantities {
vector[num_elements(y)] y_likelihood = brm_mixture_c0badac3fba85005_lpdfs(
y,
y_mixture_weights,
brm_joint_mean_rows(mu1, num_elements(y)),
brm_joint_mean_rows(exp(log(sigma)), num_elements(y)),
brm_joint_mean_rows(mu2, num_elements(y)),
brm_joint_mean_rows(exp(log(sigma)), num_elements(y))
);
vector[num_elements(y)] y_gen = brm_mixture_c0badac3fba85005_vector_rng(
y_n,
y_mixture_weights,
brm_joint_mean_rows(mu1, num_elements(y)),
brm_joint_mean_rows(exp(log(sigma)), num_elements(y)),
brm_joint_mean_rows(mu2, num_elements(y)),
brm_joint_mean_rows(exp(log(sigma)), num_elements(y))
);
}julia
#= line 0 =# Turing.@model(function brm_model(y, X_sigma)
mu1 ~ Distributions.Normal(-2, 0.1)
mu2 ~ Distributions.Normal(2, 0.1)
beta_pop ~ Distributions.product_distribution(Base.fill(Distributions.Normal(), 1))
eta_sigma = X_sigma * beta_pop
sigma = Base.exp.(eta_sigma)
begin
for i = Base.eachindex(y)
y[i] ~ Distributions.MixtureModel(Base.vect(Distributions.Normal(mu1, Base.exp(Base.log(sigma[i]))), Distributions.Normal(mu2, Base.exp(Base.log(sigma[i])))), Base.vect(0.4, 0.6))
end
end
(; sigma = sigma, mu1 = mu1, mu2 = mu2, response = y)
end)Each row contributes log_sum_exp(log(weights) + component_lpdf), with matching pointwise log likelihoods and posterior-predictive draws (select a component per row, then draw from it). The contract is deliberately narrow. Every component must be a scalar (Univariate) call of ONE Julia family: heterogeneous families are rejected because Stan's discrete densities throw on out-of-support values where Turing returns -Inf, so mixed-support mixtures would crash Stan where Turing stays finite. The family must be directly Stan-mapped, or NegativeBinomial2 / BetaBinomial2 (which are native translations); bespoke families, Uniform / Pareto (parameter-dependent continuous support), and Categorical (simplex parameters) are rejected. Binomial-family components must share one identical trial-count expression. Weights are a numeric vector summing to 1, a length-K numeric data column, or a Dirichlet-backed simplex parameter.
Adding another likelihood
Yes—when Stan already has the distribution, adding it is usually small. BRM needs an explicit Julia-type-to-Stan-name entry, plus an argument translation when the two libraries use different parameterizations. It then inherits the ordinary model, pointwise-log-likelihood, and predictive-RNG paths from StanBlocks.
The work becomes larger when Stan has no native family: the implementation must provide and test a density/mass function, a pointwise companion, and an RNG. Supporting truncated or censored also needs the relevant CDF/CCDF functions, and ragged responses need a vector-sized RNG. These are finite implementation tasks, not an architectural prohibition; the explicit gates prevent a family from appearing to fit while prediction or likelihood diagnostics silently mean something else.
Truncation and censoring
BRM preserves the standard Distributions.jl RHS composition for mathematical truncation and threshold censoring:
brm-comparison
Truncated, censored, and interval evidencejulia
bounded_evidence = (@brm begin
mu ~ 1 + x
log(sigma) ~ 1
y_truncated ~ truncated(Normal(mu, sigma); lower=0.0, upper=2.0)
y_clamped ~ censored(LogNormal(mu, sigma); lower=0.25, upper=1.8)
# Genuine interval evidence: y_lower stores the open lower endpoint and
# y_upper stores the closed upper endpoint.
y_lower ~ interval_censored(Normal(mu, sigma); upper=y_upper)
end)((;
x=[-1.0, 0.0, 1.0],
y_truncated=[0.2, 0.8, 1.4],
y_clamped=[0.25, 0.9, 1.8],
y_lower=[-0.4, 0.1, 0.8],
y_upper=[-0.1, 0.4, 1.2],
))julia
BRMI:
x: data (eltype=Float64, n=3)
mu ~ 1 + x
log(sigma) ~ 1
y_truncated ~ truncated(Normal(mu, sigma); lower=0.0, upper=2.0)
y_clamped ~ censored(LogNormal(mu, sigma); lower=0.25, upper=1.8)
y_upper: data (eltype=Float64, n=3)
y_lower ~ interval_censored(Normal(mu, sigma); upper=y_upper)julia
SBBRMI with data keys = [:x, :y_clamped, :y_lower, :y_truncated, :y_upper]
emitted @slic body:
begin
X_mu = hcat(rep_vector(1.0, num_elements(x)), x)
pop_mu ~ popefs(; X = X_mu)
mu = pop_mu
X_log_sigma = hcat(rep_vector(1.0, num_elements(y_truncated)))
pop_log_sigma ~ popefs(; X = X_log_sigma)
log_sigma = pop_log_sigma
sigma = exp(log_sigma)
y_truncated ~ truncated(normal, mu, sigma; lower = 0.0, upper = 2.0)
y_clamped ~ censored(lognormal, mu, sigma; lower = 0.25, upper = 1.8)
y_lower ~ interval_censored(normal, y_lower, y_upper, mu, sigma)
endstan
functions {
matrix hcat(
vector x,
vector y
) {
int n = dims(x)[1];
if (dims(y)[1] != n) reject("hcat: dim mismatch — `y` dim 1 (= ", dims(y)[1], ") does not match `n` (= ", n, "), inferred from `x` dim 1. `n` sizes: `x` dim 1 (= ", dims(x)[1], "), `y` dim 1 (= ", dims(y)[1], ").");
return append_col(x, y);
}
matrix hcat(vector x) {
int n = dims(x)[1];
return to_matrix(x, n, 1);
}
real conditioning_normal_lpdf(
vector y,
real lo,
real hi,
vector args1,
vector args2
) {
return sum(conditioning_normal_lpdfs(y, lo, hi, args1, args2));
}
vector conditioning_normal_lpdfs(
vector y,
real lo,
real hi,
vector args1,
vector args2
) {
int n = dims(y)[1];
return jbroadcasted_conditioning_lpdf_normal(y, lo, hi, args1, args2);
}
vector jbroadcasted_conditioning_lpdf_normal(
vector x1,
real x3,
real x4,
vector x5,
vector x6
) {
int n = dims(x1)[1];
vector[n] rv;
for(i in 1:n) {
rv[i] = conditioning_normal_lpdf(broadcasted_getindex(x1, i) |
x3,
x4,
broadcasted_getindex(x5, i),
broadcasted_getindex(x6, i)
);
}
return rv;
}
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 broadcasted_getindex(vector x, int i) {
return x[i];
}
vector conditioning_vector_normal_rng(
int anontok__1,
real lo,
real hi,
vector args1,
vector args2
) {
int n = anontok__1;
return jbroadcasted_conditioning_cell_rng_normal_rng(rep_vector(0.0, n), lo, hi, args1, args2);
}
vector jbroadcasted_conditioning_cell_rng_normal_rng(
vector x1,
real x3,
real x4,
vector x5,
vector x6
) {
int n = dims(x1)[1];
vector[n] rv;
for(i in 1:n) {
rv[i] = conditioning_cell_normal_rng(
broadcasted_getindex(x1, i),
x3,
x4,
broadcasted_getindex(x5, i),
broadcasted_getindex(x6, i)
);
}
return rv;
}
real conditioning_cell_normal_rng(
real dummy,
real lo,
real hi,
real args1,
real args2
) {
return conditioning_normal_rng(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];
}
real clamping_lognormal_lpdf(
vector y,
real lo,
real hi,
vector args1,
vector args2
) {
return sum(clamping_lognormal_lpdfs(y, lo, hi, args1, args2));
}
vector clamping_lognormal_lpdfs(
vector y,
real lo,
real hi,
vector args1,
vector args2
) {
int n = dims(y)[1];
return jbroadcasted_clamping_lpdf_lognormal(y, lo, hi, args1, args2);
}
vector jbroadcasted_clamping_lpdf_lognormal(
vector x1,
real x3,
real x4,
vector x5,
vector x6
) {
int n = dims(x1)[1];
vector[n] rv;
for(i in 1:n) {
rv[i] = clamping_lognormal_lpdf(broadcasted_getindex(x1, i) |
x3,
x4,
broadcasted_getindex(x5, i),
broadcasted_getindex(x6, i)
);
}
return rv;
}
real clamping_lognormal_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] = lognormal_lcdf(lo | args1, args2);
} else {
if((y == hi)) {
rv[1] = lognormal_lccdf(hi | args1, args2);
} else {
if((y > lo)) {
if((y < hi)) {
rv[1] = lognormal_lpdf(y | args1, args2);
}
}
}
}
return rv[1];
}
vector clamping_vector_lognormal_rng(
int anontok__1,
real lo,
real hi,
vector args1,
vector args2
) {
int n = anontok__1;
return jbroadcasted_clamping_cell_rng_lognormal_rng(rep_vector(0.0, n), lo, hi, args1, args2);
}
vector jbroadcasted_clamping_cell_rng_lognormal_rng(
vector x1,
real x3,
real x4,
vector x5,
vector x6
) {
int n = dims(x1)[1];
vector[n] rv;
for(i in 1:n) {
rv[i] = clamping_cell_lognormal_rng(
broadcasted_getindex(x1, i),
x3,
x4,
broadcasted_getindex(x5, i),
broadcasted_getindex(x6, i)
);
}
return rv;
}
real clamping_cell_lognormal_rng(
real dummy,
real lo,
real hi,
real args1,
real args2
) {
return clamping_lognormal_rng(lo, hi, args1, args2);
}
real clamping_lognormal_rng(
real lo,
real hi,
real args1,
real args2
) {
vector[1] draw;
draw[1] = lognormal_rng(args1, args2);
if((draw[1] < lo)) {
draw[1] = lo;
} else {
if((draw[1] > hi)) {
draw[1] = hi;
}
}
return draw[1];
}
real interval_evidence_impl_normal_lpdf(
vector y,
vector lo,
vector hi,
vector args1,
vector args2
) {
return sum(interval_evidence_impl_normal_lpdfs(y, lo, hi, args1, args2));
}
vector interval_evidence_impl_normal_lpdfs(
vector y,
vector lo,
vector hi,
vector args1,
vector args2
) {
int n = dims(y)[1];
return jbroadcasted_interval_evidence_impl_lpdf_normal(y, lo, hi, args1, args2);
}
vector jbroadcasted_interval_evidence_impl_lpdf_normal(
vector x1,
vector x3,
vector x4,
vector x5,
vector x6
) {
int n = dims(x1)[1];
vector[n] rv;
for(i in 1:n) {
rv[i] = interval_evidence_impl_normal_lpdf(broadcasted_getindex(x1, i) |
broadcasted_getindex(x3, i),
broadcasted_getindex(x4, i),
broadcasted_getindex(x5, i),
broadcasted_getindex(x6, i)
);
}
return rv;
}
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];
}
vector interval_evidence_impl_vector_normal_rng(
int anontok__1,
vector lo,
vector hi,
vector args1,
vector args2
) {
int n = anontok__1;
return normal_vector_rng(n, args1, args2);
}
vector normal_vector_rng(
int anontok__1,
vector a,
vector b
) {
int n = anontok__1;
if((n == 0)) {
vector[n] rv;
return rv;
} else {
return to_vector(normal_rng(a, b));
}
}
}
data {
int x_n;
vector[x_n] x;
int y_truncated_n;
vector[y_truncated_n] y_truncated;
int y_clamped_n;
vector[y_clamped_n] y_clamped;
int y_lower_n;
vector[y_lower_n] y_lower;
int y_upper_n;
vector[y_upper_n] y_upper;
}
transformed data {
matrix[x_n, 2] X_mu = hcat(rep_vector(1.0, num_elements(x)), x);
int pop_mu_n_covariates = 2;
matrix[num_elements(y_truncated), 1] X_log_sigma = hcat(rep_vector(1.0, num_elements(y_truncated)));
int pop_log_sigma_n_covariates = 1;
}
parameters {
vector[pop_mu_n_covariates] pop_mu_beta_pop;
vector[pop_log_sigma_n_covariates] pop_log_sigma_beta_pop;
}
transformed parameters {
vector[x_n] pop_mu = (X_mu * pop_mu_beta_pop);
vector[x_n] mu = pop_mu;
vector[num_elements(y_truncated)] pop_log_sigma = (X_log_sigma * pop_log_sigma_beta_pop);
vector[num_elements(y_truncated)] log_sigma = pop_log_sigma;
vector[num_elements(y_truncated)] sigma = exp(log_sigma);
}
model {
pop_mu_beta_pop ~ std_normal();
pop_log_sigma_beta_pop ~ std_normal();
y_truncated ~ conditioning_normal(0.0, 2.0, mu, sigma);
y_clamped ~ clamping_lognormal(0.25, 1.8, mu, sigma);
y_lower ~ interval_evidence_impl_normal(y_lower, y_upper, mu, sigma);
}
generated quantities {
vector[y_truncated_n] y_truncated_likelihood = conditioning_normal_lpdfs(y_truncated, 0.0, 2.0, mu, sigma);
vector[y_truncated_n] y_truncated_gen = conditioning_vector_normal_rng(y_truncated_n, 0.0, 2.0, mu, sigma);
vector[y_clamped_n] y_clamped_likelihood = clamping_lognormal_lpdfs(y_clamped, 0.25, 1.8, mu, sigma);
vector[y_clamped_n] y_clamped_gen = clamping_vector_lognormal_rng(y_clamped_n, 0.25, 1.8, mu, sigma);
vector[y_lower_n] y_lower_likelihood = interval_evidence_impl_normal_lpdfs(y_lower, y_lower, y_upper, mu, sigma);
vector[y_lower_n] y_lower_gen = interval_evidence_impl_vector_normal_rng(y_lower_n, y_lower, y_upper, mu, sigma);
}julia
#= line 0 =# Turing.@model(function brm_multi_model(y_1, y_2, y_3, X_mu, X_sigma, lower_y_truncated, upper_y_truncated, lower_y_clamped, upper_y_clamped, upper_y_lower)
beta_pop_mu ~ Distributions.product_distribution(Base.fill(Distributions.Normal(), 2))
eta_mu = X_mu * beta_pop_mu
mu = eta_mu
beta_pop_sigma ~ Distributions.product_distribution(Base.fill(Distributions.Normal(), 1))
eta_sigma = X_sigma * beta_pop_sigma
sigma = Base.exp.(eta_sigma)
begin
for i = Base.eachindex(y_1)
y_1[i] ~ Distributions.truncated(Distributions.Normal(mu[i], sigma[i]); lower = lower_y_truncated, upper = upper_y_truncated)
end
end
begin
for i = Base.eachindex(y_2)
y_2[i] ~ Distributions.censored(Distributions.LogNormal(mu[i], sigma[i]); lower = lower_y_clamped, upper = upper_y_clamped)
end
end
begin
for i = Base.eachindex(y_3)
y_3[i] ~ BayesianRegressionModelsTuringExt._brm_interval_evidence(Distributions.Normal(mu[i], sigma[i]), upper_y_lower[i])
end
end
(; responses = ((; mu = mu, sigma = sigma, response = y_1), (; mu = mu, sigma = sigma, response = y_2), (; mu = mu, sigma = sigma, response = y_3)))
end)These are three different likelihood contracts:
truncated(d; lower, upper)conditionsdon the inclusive bounds and predicts from that conditional distribution;censored(d; lower, upper)is the distribution ofclamp(X, lower, upper)and predicts clamped values;interval_censored(d; upper)contributeslog(CDF(upper) - CDF(response))for the genuine interval observation(response, upper], while prediction remains on the uncoarsened base scale.
The same marker has a separate predictor-side form, interval_censored(x; upper=lloq), for a quantified/BLOQ covariate. It allocates bounded latent predictor values on rows where x == lloq rather than changing a response likelihood; see Interval-censored predictor.
Bounds may be numeric literals or observed row-wise columns. The initial family-gated surface covers Normal, LogNormal, Exponential, Weibull, and Poisson; BRM rejects other base families until their aggregate density, pointwise likelihood, CDF/CCDF, generated prediction, and stanc paths are all tested. BRM's eager two-sided bound check accepts lower <= upper, while the StanBlocks producer requires a non-degenerate interval with lower < upper; equal bounds are therefore rejected during Stan lowering.
This composition is implemented only by the SBBRMI Stan backend. Do not use VBRMI for these formulas: it currently does not reject truncated or censored at construction and can return log densities for a misinterpreted model; interval_censored may fail only when the density is evaluated. A VBRMI result is therefore not a valid cross-check of an SBBRMI fit. The legacy TruncatedNormal marker remains a separate censored-Normal compatibility surface.
The backend compatibility floor for this surface is StanBlocks 0eaebfae904d3bffab150dfa2c59632ac783b992, where the public distribution-HOF tokens became truncated, censored, and interval_censored with no aliases. BRM revisions containing this lowering must be co-pinned with that StanBlocks commit or later; the preceding StanBlocks 9b879d5e expects the older internal token spellings and is intentionally incompatible.
Concise categorical regression
CategoricalLogit accepts an explicit nested @brm(...) predictor formula:
brm-comparison
Nested categorical-logit formulajulia
categorical_data = (;
x = [-1.0, -0.5, 0.0, 0.5, 1.0, 1.5],
y = ["b", "a", "c", "b", "c", "a"],
)
categorical_model = @brm categorical_data begin
y ~ CategoricalLogit(@brm(1 + x))
endjulia
BRMI:
x: data (eltype=Float64, n=6)
y_nested_arg1_class2 ~ 1 + x
y_nested_arg1_class3 ~ 1 + x
y ~ CategoricalLogit(y_nested_arg1_class2, y_nested_arg1_class3)julia
SBBRMI with data keys = [:x, :y]
emitted @slic body:
begin
X_y_nested_arg1_class2 = hcat(rep_vector(1.0, num_elements(x)), x)
pop_y_nested_arg1_class2 ~ popefs(; X = X_y_nested_arg1_class2)
y_nested_arg1_class2 = pop_y_nested_arg1_class2
X_y_nested_arg1_class3 = hcat(rep_vector(1.0, num_elements(x)), x)
pop_y_nested_arg1_class3 ~ popefs(; X = X_y_nested_arg1_class3)
y_nested_arg1_class3 = pop_y_nested_arg1_class3
y_categorical_logits = adjoint(hcat(rep_vector(0.0, num_elements(y)), y_nested_arg1_class2, y_nested_arg1_class3))
y ~ categorical_logit(y_categorical_logits)
endstan
functions {
matrix hcat(
vector x,
vector y
) {
int n = dims(x)[1];
if (dims(y)[1] != n) reject("hcat: dim mismatch — `y` dim 1 (= ", dims(y)[1], ") does not match `n` (= ", n, "), inferred from `x` dim 1. `n` sizes: `x` dim 1 (= ", dims(x)[1], "), `y` dim 1 (= ", dims(y)[1], ").");
return append_col(x, y);
}
matrix hcat(vector x, vector y, vector z) {
return hcat(hcat(x, y), z);
}
matrix hcat(
matrix x,
vector y
) {
int m = dims(x)[1];
int n = dims(x)[2];
if (dims(y)[1] != m) reject("hcat: dim mismatch — `y` dim 1 (= ", dims(y)[1], ") does not match `m` (= ", m, "), inferred from `x` dim 1. `m` sizes: `x` dim 1 (= ", dims(x)[1], "), `y` dim 1 (= ", dims(y)[1], ").");
return append_col(x, y);
}
real categorical_logit_lpmf(
array[] int y,
matrix eta
) {
int n = dims(y)[1];
if (dims(eta)[2] != n) reject("categorical_logit_lpmf: dim mismatch — `eta` dim 2 (= ", dims(eta)[2], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `eta` dim 2 (= ", dims(eta)[2], ").");
real rv = 0.0;
for(i in 1:n) {
rv += categorical_logit_lpmf(y[i] | eta[:, i]);
}
return rv;
}
vector categorical_logit_lpmfs(
array[] int y,
matrix eta
) {
int n = dims(y)[1];
if (dims(eta)[2] != n) reject("categorical_logit_lpmfs: dim mismatch — `eta` dim 2 (= ", dims(eta)[2], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `eta` dim 2 (= ", dims(eta)[2], ").");
vector[n] rv;
for(i in 1:n) {
rv[i] = categorical_logit_lpmf(y[i] | eta[:, i]);
}
return rv;
}
array[] int categorical_logit_int_rng(
int anontok__1,
matrix eta
) {
int n = anontok__1;
if (dims(eta)[2] != n) reject("categorical_logit_rng: dim mismatch — `eta` dim 2 (= ", dims(eta)[2], ") does not match `n` (= ", n, "), inferred from `anontok__1` dim 1. `n` sizes: `anontok__1` dim 1 (= ", anontok__1, "), `eta` dim 2 (= ", dims(eta)[2], ").");
array[n] int rv;
for(i in 1:n) {
rv[i] = categorical_logit_rng(eta[:, i]);
}
return rv;
}
}
data {
int x_n;
vector[x_n] x;
int y_n;
array[y_n] int y;
}
transformed data {
matrix[x_n, 2] X_y_nested_arg1_class2 = hcat(rep_vector(1.0, num_elements(x)), x);
int pop_y_nested_arg1_class2_n_covariates = 2;
matrix[x_n, 2] X_y_nested_arg1_class3 = hcat(rep_vector(1.0, num_elements(x)), x);
int pop_y_nested_arg1_class3_n_covariates = 2;
}
parameters {
vector[pop_y_nested_arg1_class2_n_covariates] pop_y_nested_arg1_class2_beta_pop;
vector[pop_y_nested_arg1_class3_n_covariates] pop_y_nested_arg1_class3_beta_pop;
}
transformed parameters {
vector[x_n] pop_y_nested_arg1_class2 = (X_y_nested_arg1_class2 * pop_y_nested_arg1_class2_beta_pop);
vector[x_n] y_nested_arg1_class2 = pop_y_nested_arg1_class2;
vector[x_n] pop_y_nested_arg1_class3 = (X_y_nested_arg1_class3 * pop_y_nested_arg1_class3_beta_pop);
vector[x_n] y_nested_arg1_class3 = pop_y_nested_arg1_class3;
matrix[(2 + 1), x_n] y_categorical_logits = (hcat(rep_vector(0.0, num_elements(y)), y_nested_arg1_class2, y_nested_arg1_class3)');
}
model {
pop_y_nested_arg1_class2_beta_pop ~ std_normal();
pop_y_nested_arg1_class3_beta_pop ~ std_normal();
y ~ categorical_logit(y_categorical_logits);
}
generated quantities {
vector[x_n] y_likelihood = categorical_logit_lpmfs(y, y_categorical_logits);
array[x_n] int y_gen = categorical_logit_int_rng(y_n, y_categorical_logits);
}julia
#= line 0 =# Turing.@model(function brm_model(y, X_y_nested_arg1_class2, X_y_nested_arg1_class3)
beta_pop ~ Distributions.product_distribution(Base.fill(Distributions.Normal(), 2))
eta_y_nested_arg1_class2 = X_y_nested_arg1_class2 * beta_pop
y_nested_arg1_class2 = eta_y_nested_arg1_class2
beta_pop_y_nested_arg1_class3 ~ Distributions.product_distribution(Base.fill(Distributions.Normal(), 2))
eta_y_nested_arg1_class3 = X_y_nested_arg1_class3 * beta_pop_y_nested_arg1_class3
y_nested_arg1_class3 = eta_y_nested_arg1_class3
begin
for i = Base.eachindex(y)
y[i] ~ BayesianRegressionModels.CategoricalLogit(y_nested_arg1_class2[i], y_nested_arg1_class3[i])
end
end
(; y_nested_arg1_class2 = y_nested_arg1_class2, y_nested_arg1_class3 = y_nested_arg1_class3, response = y)
end)For an outcome with sort(unique(y)) for the fitted order; a CategoricalVector uses its declared level order. The latter is the way to select a reference level deliberately.
For example, a three-level outcome above is equivalent in model structure to:
brm-comparison
Explicit categorical-logit predictorsjulia
explicit_data = (;
x=[-1.0, -0.5, 0.0, 0.5, 1.0, 1.5],
y=["b", "a", "c", "b", "c", "a"],
)
explicit_model = @brm explicit_data begin
y_nested_arg1_class2 ~ 1 + x
y_nested_arg1_class3 ~ 1 + x
y ~ CategoricalLogit(y_nested_arg1_class2, y_nested_arg1_class3)
endjulia
BRMI:
x: data (eltype=Float64, n=6)
y_nested_arg1_class2 ~ 1 + x
y_nested_arg1_class3 ~ 1 + x
y ~ CategoricalLogit(y_nested_arg1_class2, y_nested_arg1_class3)julia
SBBRMI with data keys = [:x, :y]
emitted @slic body:
begin
X_y_nested_arg1_class2 = hcat(rep_vector(1.0, num_elements(x)), x)
pop_y_nested_arg1_class2 ~ popefs(; X = X_y_nested_arg1_class2)
y_nested_arg1_class2 = pop_y_nested_arg1_class2
X_y_nested_arg1_class3 = hcat(rep_vector(1.0, num_elements(x)), x)
pop_y_nested_arg1_class3 ~ popefs(; X = X_y_nested_arg1_class3)
y_nested_arg1_class3 = pop_y_nested_arg1_class3
y_categorical_logits = adjoint(hcat(rep_vector(0.0, num_elements(y)), y_nested_arg1_class2, y_nested_arg1_class3))
y ~ categorical_logit(y_categorical_logits)
endstan
functions {
matrix hcat(
vector x,
vector y
) {
int n = dims(x)[1];
if (dims(y)[1] != n) reject("hcat: dim mismatch — `y` dim 1 (= ", dims(y)[1], ") does not match `n` (= ", n, "), inferred from `x` dim 1. `n` sizes: `x` dim 1 (= ", dims(x)[1], "), `y` dim 1 (= ", dims(y)[1], ").");
return append_col(x, y);
}
matrix hcat(vector x, vector y, vector z) {
return hcat(hcat(x, y), z);
}
matrix hcat(
matrix x,
vector y
) {
int m = dims(x)[1];
int n = dims(x)[2];
if (dims(y)[1] != m) reject("hcat: dim mismatch — `y` dim 1 (= ", dims(y)[1], ") does not match `m` (= ", m, "), inferred from `x` dim 1. `m` sizes: `x` dim 1 (= ", dims(x)[1], "), `y` dim 1 (= ", dims(y)[1], ").");
return append_col(x, y);
}
real categorical_logit_lpmf(
array[] int y,
matrix eta
) {
int n = dims(y)[1];
if (dims(eta)[2] != n) reject("categorical_logit_lpmf: dim mismatch — `eta` dim 2 (= ", dims(eta)[2], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `eta` dim 2 (= ", dims(eta)[2], ").");
real rv = 0.0;
for(i in 1:n) {
rv += categorical_logit_lpmf(y[i] | eta[:, i]);
}
return rv;
}
vector categorical_logit_lpmfs(
array[] int y,
matrix eta
) {
int n = dims(y)[1];
if (dims(eta)[2] != n) reject("categorical_logit_lpmfs: dim mismatch — `eta` dim 2 (= ", dims(eta)[2], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `eta` dim 2 (= ", dims(eta)[2], ").");
vector[n] rv;
for(i in 1:n) {
rv[i] = categorical_logit_lpmf(y[i] | eta[:, i]);
}
return rv;
}
array[] int categorical_logit_int_rng(
int anontok__1,
matrix eta
) {
int n = anontok__1;
if (dims(eta)[2] != n) reject("categorical_logit_rng: dim mismatch — `eta` dim 2 (= ", dims(eta)[2], ") does not match `n` (= ", n, "), inferred from `anontok__1` dim 1. `n` sizes: `anontok__1` dim 1 (= ", anontok__1, "), `eta` dim 2 (= ", dims(eta)[2], ").");
array[n] int rv;
for(i in 1:n) {
rv[i] = categorical_logit_rng(eta[:, i]);
}
return rv;
}
}
data {
int x_n;
vector[x_n] x;
int y_n;
array[y_n] int y;
}
transformed data {
matrix[x_n, 2] X_y_nested_arg1_class2 = hcat(rep_vector(1.0, num_elements(x)), x);
int pop_y_nested_arg1_class2_n_covariates = 2;
matrix[x_n, 2] X_y_nested_arg1_class3 = hcat(rep_vector(1.0, num_elements(x)), x);
int pop_y_nested_arg1_class3_n_covariates = 2;
}
parameters {
vector[pop_y_nested_arg1_class2_n_covariates] pop_y_nested_arg1_class2_beta_pop;
vector[pop_y_nested_arg1_class3_n_covariates] pop_y_nested_arg1_class3_beta_pop;
}
transformed parameters {
vector[x_n] pop_y_nested_arg1_class2 = (X_y_nested_arg1_class2 * pop_y_nested_arg1_class2_beta_pop);
vector[x_n] y_nested_arg1_class2 = pop_y_nested_arg1_class2;
vector[x_n] pop_y_nested_arg1_class3 = (X_y_nested_arg1_class3 * pop_y_nested_arg1_class3_beta_pop);
vector[x_n] y_nested_arg1_class3 = pop_y_nested_arg1_class3;
matrix[(2 + 1), x_n] y_categorical_logits = (hcat(rep_vector(0.0, num_elements(y)), y_nested_arg1_class2, y_nested_arg1_class3)');
}
model {
pop_y_nested_arg1_class2_beta_pop ~ std_normal();
pop_y_nested_arg1_class3_beta_pop ~ std_normal();
y ~ categorical_logit(y_categorical_logits);
}
generated quantities {
vector[x_n] y_likelihood = categorical_logit_lpmfs(y, y_categorical_logits);
array[x_n] int y_gen = categorical_logit_int_rng(y_n, y_categorical_logits);
}julia
#= line 0 =# Turing.@model(function brm_model(y, X_y_nested_arg1_class2, X_y_nested_arg1_class3)
beta_pop ~ Distributions.product_distribution(Base.fill(Distributions.Normal(), 2))
eta_y_nested_arg1_class2 = X_y_nested_arg1_class2 * beta_pop
y_nested_arg1_class2 = eta_y_nested_arg1_class2
beta_pop_y_nested_arg1_class3 ~ Distributions.product_distribution(Base.fill(Distributions.Normal(), 2))
eta_y_nested_arg1_class3 = X_y_nested_arg1_class3 * beta_pop_y_nested_arg1_class3
y_nested_arg1_class3 = eta_y_nested_arg1_class3
begin
for i = Base.eachindex(y)
y[i] ~ BayesianRegressionModels.CategoricalLogit(y_nested_arg1_class2[i], y_nested_arg1_class3[i])
end
end
(; y_nested_arg1_class2 = y_nested_arg1_class2, y_nested_arg1_class3 = y_nested_arg1_class3, response = y)
end)The generated names are deterministic implementation names; use the explicit form when those predictor names are part of another formula. Only nested @brm(...) opts into predictor-formula interpretation. Thus CategoricalLogit(1 + x) remains an ordinary expression and is rejected by the categorical backend, rather than silently acquiring coefficients.
The marker is not categorical-specific. It selects formula interpretation at one family-argument position while surrounding expressions retain their usual meaning. For example, a distributional Normal model can make both predictors concise while keeping the positive scale link explicit:
brm-comparison
Nested distributional predictorsjulia
distributional_data = (; x=[-1.0, 0.0, 1.0], y=[-0.2, 0.3, 1.1])
distributional_model = @brm distributional_data begin
y ~ Normal(@brm(1 + x), exp(@brm(1)))
endjulia
BRMI:
x: data (eltype=Float64, n=3)
y_nested_arg1 ~ 1 + x
y_nested_arg2_arg1 ~ 1
y ~ Normal(y_nested_arg1, exp(y_nested_arg2_arg1))julia
SBBRMI with data keys = [:x, :y]
emitted @slic body:
begin
X_y_nested_arg1 = hcat(rep_vector(1.0, num_elements(x)), x)
pop_y_nested_arg1 ~ _popefs_coefs(; X = X_y_nested_arg1)
X_y_nested_arg2_arg1 = hcat(rep_vector(1.0, num_elements(y)))
pop_y_nested_arg2_arg1 ~ popefs(; X = X_y_nested_arg2_arg1)
y_nested_arg2_arg1 = pop_y_nested_arg2_arg1
y ~ normal_id_glm(X_y_nested_arg1, 0.0, pop_y_nested_arg1, (exp)(y_nested_arg2_arg1))
y_nested_arg1 = X_y_nested_arg1 * pop_y_nested_arg1
endstan
functions {
matrix hcat(
vector x,
vector y
) {
int n = dims(x)[1];
if (dims(y)[1] != n) reject("hcat: dim mismatch — `y` dim 1 (= ", dims(y)[1], ") does not match `n` (= ", n, "), inferred from `x` dim 1. `n` sizes: `x` dim 1 (= ", dims(x)[1], "), `y` dim 1 (= ", dims(y)[1], ").");
return append_col(x, y);
}
matrix hcat(vector x) {
int n = dims(x)[1];
return to_matrix(x, n, 1);
}
vector normal_id_glm_lpdfs(
vector y,
matrix X,
real alpha,
vector beta,
vector 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,
vector 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], ").");
if((m == 0)) {
vector[m] rv;
return rv;
} else {
return normal_id_glm_rng(X, alpha, beta, sigma);
}
}
vector normal_id_glm_rng(
matrix X,
real alpha,
vector beta,
vector sigma
) {
int m = dims(X)[1];
if((m == 0)) {
vector[m] rv;
return rv;
} else {
return to_vector(normal_rng((rep_vector(alpha, m) + (X * beta)), sigma));
}
}
}
data {
int x_n;
vector[x_n] x;
int y_n;
vector[y_n] y;
}
transformed data {
matrix[x_n, 2] X_y_nested_arg1 = hcat(rep_vector(1.0, num_elements(x)), x);
int pop_y_nested_arg1_n_covariates = 2;
matrix[num_elements(y), 1] X_y_nested_arg2_arg1 = hcat(rep_vector(1.0, num_elements(y)));
int pop_y_nested_arg2_arg1_n_covariates = 1;
}
parameters {
vector[pop_y_nested_arg1_n_covariates] pop_y_nested_arg1_beta_pop;
vector[pop_y_nested_arg2_arg1_n_covariates] pop_y_nested_arg2_arg1_beta_pop;
}
transformed parameters {
vector[pop_y_nested_arg1_n_covariates] pop_y_nested_arg1 = pop_y_nested_arg1_beta_pop;
vector[num_elements(y)] pop_y_nested_arg2_arg1 = (X_y_nested_arg2_arg1 * pop_y_nested_arg2_arg1_beta_pop);
vector[num_elements(y)] y_nested_arg2_arg1 = pop_y_nested_arg2_arg1;
}
model {
pop_y_nested_arg1_beta_pop ~ std_normal();
pop_y_nested_arg2_arg1_beta_pop ~ std_normal();
y ~ normal_id_glm(X_y_nested_arg1, 0.0, pop_y_nested_arg1, exp(y_nested_arg2_arg1));
}
generated quantities {
vector[x_n] y_likelihood = normal_id_glm_lpdfs(y, X_y_nested_arg1, 0.0, pop_y_nested_arg1, exp(y_nested_arg2_arg1));
vector[x_n] y_gen = normal_id_glm_vector_rng(y_n, X_y_nested_arg1, 0.0, pop_y_nested_arg1, exp(y_nested_arg2_arg1));
vector[x_n] y_nested_arg1 = (X_y_nested_arg1 * pop_y_nested_arg1);
}julia
#= line 0 =# Turing.@model(function brm_model(y, X_y_nested_arg1, X_y_nested_arg2_arg1)
beta_pop ~ Distributions.product_distribution(Base.fill(Distributions.Normal(), 2))
eta_y_nested_arg1 = X_y_nested_arg1 * beta_pop
y_nested_arg1 = eta_y_nested_arg1
beta_pop_y_nested_arg2_arg1 ~ Distributions.product_distribution(Base.fill(Distributions.Normal(), 1))
eta_y_nested_arg2_arg1 = X_y_nested_arg2_arg1 * beta_pop_y_nested_arg2_arg1
y_nested_arg2_arg1 = eta_y_nested_arg2_arg1
begin
for i = Base.eachindex(y)
y[i] ~ Distributions.Normal(y_nested_arg1[i], Base.exp(y_nested_arg2_arg1[i]))
end
end
(; y_nested_arg1 = y_nested_arg1, y_nested_arg2_arg1 = y_nested_arg2_arg1, response = y)
end)This introduces distinct scalar predictors for location and log-scale, then passes exp(log_scale) to Normal; nested @brm never inserts a link. A standalone fragment such as @brm(1 + x) is not yet a first-class value and produces a targeted error outside an enclosing model.
This is the same broad model structure expressed by a top-level categorical formula in brms (y ~ 1 + x, family = categorical(link = "logit")) or Bambi ("y ~ 1 + x", family="categorical"). Defaults for priors, contrasts, and reference-level selection are package-specific; BRM does not import those defaults implicitly.
Ordinal structure and link composition
BRM treats the ordinal probability construction and inverse link as separate typed choices. A cumulative probit model is:
brm-comparison
Ordinal cumulative-probit modeljulia
ordinal_data = (; x=[-1.0, -0.5, 0.0, 0.5, 1.0], y=[1, 1, 2, 3, 3])
ordinal_model = @brm ordinal_data begin
eta ~ 0 + x
y ~ Ordinal(Cumulative(), ProbitLink(), eta)
endjulia
BRMI:
x: data (eltype=Float64, n=5)
eta ~ 0 + x
y ~ Ordinal(Cumulative(), ProbitLink(), eta)julia
SBBRMI with data keys = [:x, :y]
emitted @slic body:
begin
X_eta = hcat(x)
pop_eta ~ popefs(; X = X_eta)
eta = pop_eta
y_thresholds::ordered[2] ~ std_normal()
y_threshold_effect = rep_matrix(0.0, num_elements(y), 2)
y ~ brm_ordinal(eta, y_thresholds, 1.0, 1, 2, y_threshold_effect)
endstan
functions {
matrix hcat(vector x) {
int n = dims(x)[1];
return to_matrix(x, n, 1);
}
real brm_ordinal_lpmf(
array[] int y,
vector eta,
vector thresholds,
real discrimination,
int structure,
int link,
matrix threshold_effect
) {
int n = dims(y)[1];
int k = dims(thresholds)[1];
if (dims(eta)[1] != n) reject("brm_ordinal_lpmf: dim mismatch — `eta` dim 1 (= ", dims(eta)[1], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `eta` dim 1 (= ", dims(eta)[1], "), `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ").");
if (dims(threshold_effect)[1] != n) reject("brm_ordinal_lpmf: dim mismatch — `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `eta` dim 1 (= ", dims(eta)[1], "), `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ").");
if (dims(threshold_effect)[2] != k) reject("brm_ordinal_lpmf: dim mismatch — `threshold_effect` dim 2 (= ", dims(threshold_effect)[2], ") does not match `k` (= ", k, "), inferred from `thresholds` dim 1. `k` sizes: `thresholds` dim 1 (= ", dims(thresholds)[1], "), `threshold_effect` dim 2 (= ", dims(threshold_effect)[2], ").");
return brm_ordinal_lpmf(y |
eta,
thresholds,
rep_vector(discrimination, n),
structure,
link,
threshold_effect
);
}
real brm_ordinal_lpmf(
array[] int y,
vector eta,
vector thresholds,
vector discrimination,
int structure,
int link,
matrix threshold_effect
) {
int n = dims(y)[1];
int k = dims(thresholds)[1];
if (dims(eta)[1] != n) reject("brm_ordinal_lpmf: dim mismatch — `eta` dim 1 (= ", dims(eta)[1], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `eta` dim 1 (= ", dims(eta)[1], "), `discrimination` dim 1 (= ", dims(discrimination)[1], "), `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ").");
if (dims(discrimination)[1] != n) reject("brm_ordinal_lpmf: dim mismatch — `discrimination` dim 1 (= ", dims(discrimination)[1], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `eta` dim 1 (= ", dims(eta)[1], "), `discrimination` dim 1 (= ", dims(discrimination)[1], "), `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ").");
if (dims(threshold_effect)[1] != n) reject("brm_ordinal_lpmf: dim mismatch — `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `eta` dim 1 (= ", dims(eta)[1], "), `discrimination` dim 1 (= ", dims(discrimination)[1], "), `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ").");
if (dims(threshold_effect)[2] != k) reject("brm_ordinal_lpmf: dim mismatch — `threshold_effect` dim 2 (= ", dims(threshold_effect)[2], ") does not match `k` (= ", k, "), inferred from `thresholds` dim 1. `k` sizes: `thresholds` dim 1 (= ", dims(thresholds)[1], "), `threshold_effect` dim 2 (= ", dims(threshold_effect)[2], ").");
real rv = 0.0;
for(i in 1:n) {
rv += brm_ordinal_lpmf(y[i] |
eta[i],
thresholds,
discrimination[i],
structure,
link,
to_vector(threshold_effect[i, :])
);
}
return rv;
}
real brm_ordinal_lpmf(
int y,
real eta,
vector thresholds,
real discrimination,
int structure,
int link,
vector threshold_effect
) {
int k = dims(thresholds)[1];
if (dims(threshold_effect)[1] != k) reject("brm_ordinal_lpmf: dim mismatch — `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ") does not match `k` (= ", k, "), inferred from `thresholds` dim 1. `k` sizes: `thresholds` dim 1 (= ", dims(thresholds)[1], "), `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ").");
int K = (k + 1);
if((discrimination <= 0.0)) {
return negative_infinity();
} else {
if((y < 1)) {
return negative_infinity();
} else {
if((y > K)) {
return negative_infinity();
} else {
if((structure == 1)) {
if((link == 1)) {
return ordered_logistic_lpmf(y | (discrimination * eta), (discrimination .* thresholds));
} else {
if((y == 1)) {
real z_first = (discrimination * (thresholds[1] - eta));
return brm_ordinal_logcdf(z_first, link);
} else {
if((y == K)) {
real z_last = (discrimination * (thresholds[k] - eta));
return brm_ordinal_logccdf(z_last, link);
} else {
real z_hi = (discrimination * (thresholds[y] - eta));
real z_lo = (discrimination * (thresholds[(y - 1)] - eta));
return log_diff_exp(brm_ordinal_logcdf(z_hi, link), brm_ordinal_logcdf(z_lo, link));
}
}
}
} else {
real rv = 0.0;
for(j in 1:k) {
real z_stage = (discrimination * ((thresholds[j] - eta) - threshold_effect[j]));
if((j < y)) {
rv += brm_ordinal_logccdf(z_stage, link);
} else {
if((j == y)) {
rv += brm_ordinal_logcdf(z_stage, link);
}
}
}
return rv;
}
}
}
}
}
real brm_ordinal_logcdf(
real z,
int link
) {
if((link == 1)) {
return log_inv_logit(z);
} else {
if((link == 2)) {
return normal_lcdf(z | 0.0, 1.0);
} else {
return log1m_exp((-exp(z)));
}
}
}
real brm_ordinal_logccdf(
real z,
int link
) {
if((link == 1)) {
return log_inv_logit((-z));
} else {
if((link == 2)) {
return normal_lccdf(z | 0.0, 1.0);
} else {
return (-exp(z));
}
}
}
vector brm_ordinal_lpmfs(
array[] int y,
vector eta,
vector thresholds,
real discrimination,
int structure,
int link,
matrix threshold_effect
) {
int n = dims(y)[1];
int k = dims(thresholds)[1];
if (dims(eta)[1] != n) reject("brm_ordinal_lpmfs: dim mismatch — `eta` dim 1 (= ", dims(eta)[1], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `eta` dim 1 (= ", dims(eta)[1], "), `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ").");
if (dims(threshold_effect)[1] != n) reject("brm_ordinal_lpmfs: dim mismatch — `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `eta` dim 1 (= ", dims(eta)[1], "), `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ").");
if (dims(threshold_effect)[2] != k) reject("brm_ordinal_lpmfs: dim mismatch — `threshold_effect` dim 2 (= ", dims(threshold_effect)[2], ") does not match `k` (= ", k, "), inferred from `thresholds` dim 1. `k` sizes: `thresholds` dim 1 (= ", dims(thresholds)[1], "), `threshold_effect` dim 2 (= ", dims(threshold_effect)[2], ").");
return brm_ordinal_lpmfs(
y,
eta,
thresholds,
rep_vector(discrimination, n),
structure,
link,
threshold_effect
);
}
vector brm_ordinal_lpmfs(
array[] int y,
vector eta,
vector thresholds,
vector discrimination,
int structure,
int link,
matrix threshold_effect
) {
int n = dims(y)[1];
int k = dims(thresholds)[1];
if (dims(eta)[1] != n) reject("brm_ordinal_lpmfs: dim mismatch — `eta` dim 1 (= ", dims(eta)[1], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `eta` dim 1 (= ", dims(eta)[1], "), `discrimination` dim 1 (= ", dims(discrimination)[1], "), `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ").");
if (dims(discrimination)[1] != n) reject("brm_ordinal_lpmfs: dim mismatch — `discrimination` dim 1 (= ", dims(discrimination)[1], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `eta` dim 1 (= ", dims(eta)[1], "), `discrimination` dim 1 (= ", dims(discrimination)[1], "), `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ").");
if (dims(threshold_effect)[1] != n) reject("brm_ordinal_lpmfs: dim mismatch — `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `eta` dim 1 (= ", dims(eta)[1], "), `discrimination` dim 1 (= ", dims(discrimination)[1], "), `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ").");
if (dims(threshold_effect)[2] != k) reject("brm_ordinal_lpmfs: dim mismatch — `threshold_effect` dim 2 (= ", dims(threshold_effect)[2], ") does not match `k` (= ", k, "), inferred from `thresholds` dim 1. `k` sizes: `thresholds` dim 1 (= ", dims(thresholds)[1], "), `threshold_effect` dim 2 (= ", dims(threshold_effect)[2], ").");
vector[n] rv;
for(i in 1:n) {
rv[i] = brm_ordinal_lpmf(y[i] |
eta[i],
thresholds,
discrimination[i],
structure,
link,
to_vector(threshold_effect[i, :])
);
}
return rv;
}
array[] int brm_ordinal_int_rng(
int anontok__1,
vector eta,
vector thresholds,
real discrimination,
int structure,
int link,
matrix threshold_effect
) {
int n = anontok__1;
int k = dims(thresholds)[1];
if (dims(eta)[1] != n) reject("brm_ordinal_rng: dim mismatch — `eta` dim 1 (= ", dims(eta)[1], ") does not match `n` (= ", n, "), inferred from `anontok__1` dim 1. `n` sizes: `anontok__1` dim 1 (= ", anontok__1, "), `eta` dim 1 (= ", dims(eta)[1], "), `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ").");
if (dims(threshold_effect)[1] != n) reject("brm_ordinal_rng: dim mismatch — `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ") does not match `n` (= ", n, "), inferred from `anontok__1` dim 1. `n` sizes: `anontok__1` dim 1 (= ", anontok__1, "), `eta` dim 1 (= ", dims(eta)[1], "), `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ").");
if (dims(threshold_effect)[2] != k) reject("brm_ordinal_rng: dim mismatch — `threshold_effect` dim 2 (= ", dims(threshold_effect)[2], ") does not match `k` (= ", k, "), inferred from `thresholds` dim 1. `k` sizes: `thresholds` dim 1 (= ", dims(thresholds)[1], "), `threshold_effect` dim 2 (= ", dims(threshold_effect)[2], ").");
return brm_ordinal_int_rng(
n,
eta,
thresholds,
rep_vector(discrimination, n),
structure,
link,
threshold_effect
);
}
array[] int brm_ordinal_int_rng(
int anontok__1,
vector eta,
vector thresholds,
vector discrimination,
int structure,
int link,
matrix threshold_effect
) {
int n = anontok__1;
int k = dims(thresholds)[1];
if (dims(eta)[1] != n) reject("brm_ordinal_rng: dim mismatch — `eta` dim 1 (= ", dims(eta)[1], ") does not match `n` (= ", n, "), inferred from `anontok__1` dim 1. `n` sizes: `anontok__1` dim 1 (= ", anontok__1, "), `eta` dim 1 (= ", dims(eta)[1], "), `discrimination` dim 1 (= ", dims(discrimination)[1], "), `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ").");
if (dims(discrimination)[1] != n) reject("brm_ordinal_rng: dim mismatch — `discrimination` dim 1 (= ", dims(discrimination)[1], ") does not match `n` (= ", n, "), inferred from `anontok__1` dim 1. `n` sizes: `anontok__1` dim 1 (= ", anontok__1, "), `eta` dim 1 (= ", dims(eta)[1], "), `discrimination` dim 1 (= ", dims(discrimination)[1], "), `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ").");
if (dims(threshold_effect)[1] != n) reject("brm_ordinal_rng: dim mismatch — `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ") does not match `n` (= ", n, "), inferred from `anontok__1` dim 1. `n` sizes: `anontok__1` dim 1 (= ", anontok__1, "), `eta` dim 1 (= ", dims(eta)[1], "), `discrimination` dim 1 (= ", dims(discrimination)[1], "), `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ").");
if (dims(threshold_effect)[2] != k) reject("brm_ordinal_rng: dim mismatch — `threshold_effect` dim 2 (= ", dims(threshold_effect)[2], ") does not match `k` (= ", k, "), inferred from `thresholds` dim 1. `k` sizes: `thresholds` dim 1 (= ", dims(thresholds)[1], "), `threshold_effect` dim 2 (= ", dims(threshold_effect)[2], ").");
array[n] int rv;
for(i in 1:n) {
rv[i] = brm_ordinal_rng(
eta[i],
thresholds,
discrimination[i],
structure,
link,
to_vector(threshold_effect[i, :])
);
}
return rv;
}
int brm_ordinal_rng(
real eta,
vector thresholds,
real discrimination,
int structure,
int link,
vector threshold_effect
) {
int k = dims(thresholds)[1];
if (dims(threshold_effect)[1] != k) reject("brm_ordinal_rng: dim mismatch — `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ") does not match `k` (= ", k, "), inferred from `thresholds` dim 1. `k` sizes: `thresholds` dim 1 (= ", dims(thresholds)[1], "), `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ").");
int K = (k + 1);
int rv = K;
if((structure == 1)) {
real u = uniform_rng(0.0, 1.0);
for(j in 1:k) {
if((rv == K)) {
real z_cumulative = (discrimination * (thresholds[j] - eta));
if((u <= brm_ordinal_cdf(z_cumulative | link))) {
rv += (j - rv);
}
}
}
} else {
for(j in 1:k) {
if((rv == K)) {
real z_stopping = (discrimination * ((thresholds[j] - eta) - threshold_effect[j]));
if((bernoulli_rng(brm_ordinal_cdf(z_stopping | link)) == 1)) {
rv += (j - rv);
}
}
}
}
return rv;
}
real brm_ordinal_cdf(
real z,
int link
) {
if((link == 1)) {
return inv_logit(z);
} else {
if((link == 2)) {
return Phi(z);
} else {
return (-expm1((-exp(z))));
}
}
}
}
data {
int x_n;
vector[x_n] x;
int y_n;
array[y_n] int y;
}
transformed data {
matrix[x_n, 1] X_eta = hcat(x);
int pop_eta_n_covariates = 1;
matrix[num_elements(y), 2] y_threshold_effect = rep_matrix(0.0, num_elements(y), 2);
}
parameters {
vector[pop_eta_n_covariates] pop_eta_beta_pop;
ordered[2] y_thresholds;
}
transformed parameters {
vector[x_n] pop_eta = (X_eta * pop_eta_beta_pop);
vector[x_n] eta = pop_eta;
}
model {
pop_eta_beta_pop ~ std_normal();
y_thresholds ~ std_normal();
y ~ brm_ordinal(eta, y_thresholds, 1.0, 1, 2, y_threshold_effect);
}
generated quantities {
vector[num_elements(y)] y_likelihood = brm_ordinal_lpmfs(y, eta, y_thresholds, 1.0, 1, 2, y_threshold_effect);
array[num_elements(y)] int y_gen = brm_ordinal_int_rng(y_n, eta, y_thresholds, 1.0, 1, 2, y_threshold_effect);
}julia
#= line 0 =# Turing.@model(function brm_model(y, X_eta, callable_1)
beta_pop ~ Distributions.product_distribution(Base.fill(Distributions.Normal(), 1))
eta_eta = X_eta * beta_pop
eta = eta_eta
y_thresholds ~ callable_1(2)
begin
for i = Base.eachindex(y)
y[i] ~ BayesianRegressionModels.Ordinal($(QuoteNode(Cumulative())), $(QuoteNode(ProbitLink())), eta[i], y_thresholds; discrimination = 1.0)
end
end
(; eta = eta, y_thresholds = y_thresholds, response = y)
end)The accepted structures are Cumulative() and StoppingRatio(). Each composes with LogitLink(), ProbitLink(), or CloglogLink(). This is intentionally a Julia-native typed surface: BRM does not copy R formula helper names or encode every structure/link pair in a new family type.
For Cumulative(), BRM estimates strictly ordered thresholds
For StoppingRatio(), the estimated stage intercepts need not be ordered and
with the final category equal to the probability of continuing through every stage. F is logistic, standard normal, or complementary-log-log according to the link tag. Both threshold vectors currently receive element-wise standard normal priors.
The thresholds already supply the model location, so the composed surface requires an intercept-free common predictor (eta ~ 0 + ...). They count as that predictor's intercept: a categorical term of eta keeps its K−1 treatment contrasts rather than the cell means an intercept-free predictor otherwise gets ("Cell means" on the overview page). A positive discrimination parameter is explicit, and its predictor is not threshold-located — switch its cell-mean coding off with cmc=false (brms' name for it) so that one group's discrimination stays pinned at exp(0) = 1, which is what identifies the scale against the free thresholds:
brm-comparison
Ordinal discrimination modeljulia
ordinal_disc_data = (;
x=[-1.0, -0.5, 0.0, 0.5, 1.0, 1.5],
group=[1, 2, 1, 2, 1, 2], y=[1, 1, 2, 2, 3, 3],
)
ordinal_disc = @brm ordinal_disc_data begin
eta ~ 0 + x
log(disc) ~ 0 + factor(group; cmc=false)
y ~ Ordinal(Cumulative(), ProbitLink(), eta;
discrimination=disc)
endjulia
BRMI:
x: data (eltype=Float64, n=6)
eta ~ 0 + x
group: data (eltype=Int64, n=6)
log(disc) ~ 0 + factor(group; cmc=false)
y ~ Ordinal(Cumulative(), ProbitLink(), eta; discrimination=disc)julia
SBBRMI with data keys = [:group, :group_idx, :group_n_levels, :x, :y]
emitted @slic body:
begin
X_eta = hcat(x)
pop_eta ~ popefs(; X = X_eta)
eta = pop_eta
cat_log_disc_group ~ _sb_cat(; x = group_idx, n_levels = group_n_levels)
log_disc = cat_log_disc_group
disc = exp(log_disc)
y_thresholds::ordered[2] ~ std_normal()
y_threshold_effect = rep_matrix(0.0, num_elements(y), 2)
y ~ brm_ordinal(eta, y_thresholds, disc, 1, 2, y_threshold_effect)
endstan
functions {
matrix hcat(vector x) {
int n = dims(x)[1];
return to_matrix(x, n, 1);
}
real brm_ordinal_lpmf(
array[] int y,
vector eta,
vector thresholds,
vector discrimination,
int structure,
int link,
matrix threshold_effect
) {
int n = dims(y)[1];
int k = dims(thresholds)[1];
if (dims(eta)[1] != n) reject("brm_ordinal_lpmf: dim mismatch — `eta` dim 1 (= ", dims(eta)[1], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `eta` dim 1 (= ", dims(eta)[1], "), `discrimination` dim 1 (= ", dims(discrimination)[1], "), `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ").");
if (dims(discrimination)[1] != n) reject("brm_ordinal_lpmf: dim mismatch — `discrimination` dim 1 (= ", dims(discrimination)[1], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `eta` dim 1 (= ", dims(eta)[1], "), `discrimination` dim 1 (= ", dims(discrimination)[1], "), `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ").");
if (dims(threshold_effect)[1] != n) reject("brm_ordinal_lpmf: dim mismatch — `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `eta` dim 1 (= ", dims(eta)[1], "), `discrimination` dim 1 (= ", dims(discrimination)[1], "), `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ").");
if (dims(threshold_effect)[2] != k) reject("brm_ordinal_lpmf: dim mismatch — `threshold_effect` dim 2 (= ", dims(threshold_effect)[2], ") does not match `k` (= ", k, "), inferred from `thresholds` dim 1. `k` sizes: `thresholds` dim 1 (= ", dims(thresholds)[1], "), `threshold_effect` dim 2 (= ", dims(threshold_effect)[2], ").");
real rv = 0.0;
for(i in 1:n) {
rv += brm_ordinal_lpmf(y[i] |
eta[i],
thresholds,
discrimination[i],
structure,
link,
to_vector(threshold_effect[i, :])
);
}
return rv;
}
real brm_ordinal_lpmf(
int y,
real eta,
vector thresholds,
real discrimination,
int structure,
int link,
vector threshold_effect
) {
int k = dims(thresholds)[1];
if (dims(threshold_effect)[1] != k) reject("brm_ordinal_lpmf: dim mismatch — `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ") does not match `k` (= ", k, "), inferred from `thresholds` dim 1. `k` sizes: `thresholds` dim 1 (= ", dims(thresholds)[1], "), `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ").");
int K = (k + 1);
if((discrimination <= 0.0)) {
return negative_infinity();
} else {
if((y < 1)) {
return negative_infinity();
} else {
if((y > K)) {
return negative_infinity();
} else {
if((structure == 1)) {
if((link == 1)) {
return ordered_logistic_lpmf(y | (discrimination * eta), (discrimination .* thresholds));
} else {
if((y == 1)) {
real z_first = (discrimination * (thresholds[1] - eta));
return brm_ordinal_logcdf(z_first, link);
} else {
if((y == K)) {
real z_last = (discrimination * (thresholds[k] - eta));
return brm_ordinal_logccdf(z_last, link);
} else {
real z_hi = (discrimination * (thresholds[y] - eta));
real z_lo = (discrimination * (thresholds[(y - 1)] - eta));
return log_diff_exp(brm_ordinal_logcdf(z_hi, link), brm_ordinal_logcdf(z_lo, link));
}
}
}
} else {
real rv = 0.0;
for(j in 1:k) {
real z_stage = (discrimination * ((thresholds[j] - eta) - threshold_effect[j]));
if((j < y)) {
rv += brm_ordinal_logccdf(z_stage, link);
} else {
if((j == y)) {
rv += brm_ordinal_logcdf(z_stage, link);
}
}
}
return rv;
}
}
}
}
}
real brm_ordinal_logcdf(
real z,
int link
) {
if((link == 1)) {
return log_inv_logit(z);
} else {
if((link == 2)) {
return normal_lcdf(z | 0.0, 1.0);
} else {
return log1m_exp((-exp(z)));
}
}
}
real brm_ordinal_logccdf(
real z,
int link
) {
if((link == 1)) {
return log_inv_logit((-z));
} else {
if((link == 2)) {
return normal_lccdf(z | 0.0, 1.0);
} else {
return (-exp(z));
}
}
}
vector brm_ordinal_lpmfs(
array[] int y,
vector eta,
vector thresholds,
vector discrimination,
int structure,
int link,
matrix threshold_effect
) {
int n = dims(y)[1];
int k = dims(thresholds)[1];
if (dims(eta)[1] != n) reject("brm_ordinal_lpmfs: dim mismatch — `eta` dim 1 (= ", dims(eta)[1], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `eta` dim 1 (= ", dims(eta)[1], "), `discrimination` dim 1 (= ", dims(discrimination)[1], "), `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ").");
if (dims(discrimination)[1] != n) reject("brm_ordinal_lpmfs: dim mismatch — `discrimination` dim 1 (= ", dims(discrimination)[1], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `eta` dim 1 (= ", dims(eta)[1], "), `discrimination` dim 1 (= ", dims(discrimination)[1], "), `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ").");
if (dims(threshold_effect)[1] != n) reject("brm_ordinal_lpmfs: dim mismatch — `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `eta` dim 1 (= ", dims(eta)[1], "), `discrimination` dim 1 (= ", dims(discrimination)[1], "), `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ").");
if (dims(threshold_effect)[2] != k) reject("brm_ordinal_lpmfs: dim mismatch — `threshold_effect` dim 2 (= ", dims(threshold_effect)[2], ") does not match `k` (= ", k, "), inferred from `thresholds` dim 1. `k` sizes: `thresholds` dim 1 (= ", dims(thresholds)[1], "), `threshold_effect` dim 2 (= ", dims(threshold_effect)[2], ").");
vector[n] rv;
for(i in 1:n) {
rv[i] = brm_ordinal_lpmf(y[i] |
eta[i],
thresholds,
discrimination[i],
structure,
link,
to_vector(threshold_effect[i, :])
);
}
return rv;
}
array[] int brm_ordinal_int_rng(
int anontok__1,
vector eta,
vector thresholds,
vector discrimination,
int structure,
int link,
matrix threshold_effect
) {
int n = anontok__1;
int k = dims(thresholds)[1];
if (dims(eta)[1] != n) reject("brm_ordinal_rng: dim mismatch — `eta` dim 1 (= ", dims(eta)[1], ") does not match `n` (= ", n, "), inferred from `anontok__1` dim 1. `n` sizes: `anontok__1` dim 1 (= ", anontok__1, "), `eta` dim 1 (= ", dims(eta)[1], "), `discrimination` dim 1 (= ", dims(discrimination)[1], "), `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ").");
if (dims(discrimination)[1] != n) reject("brm_ordinal_rng: dim mismatch — `discrimination` dim 1 (= ", dims(discrimination)[1], ") does not match `n` (= ", n, "), inferred from `anontok__1` dim 1. `n` sizes: `anontok__1` dim 1 (= ", anontok__1, "), `eta` dim 1 (= ", dims(eta)[1], "), `discrimination` dim 1 (= ", dims(discrimination)[1], "), `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ").");
if (dims(threshold_effect)[1] != n) reject("brm_ordinal_rng: dim mismatch — `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ") does not match `n` (= ", n, "), inferred from `anontok__1` dim 1. `n` sizes: `anontok__1` dim 1 (= ", anontok__1, "), `eta` dim 1 (= ", dims(eta)[1], "), `discrimination` dim 1 (= ", dims(discrimination)[1], "), `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ").");
if (dims(threshold_effect)[2] != k) reject("brm_ordinal_rng: dim mismatch — `threshold_effect` dim 2 (= ", dims(threshold_effect)[2], ") does not match `k` (= ", k, "), inferred from `thresholds` dim 1. `k` sizes: `thresholds` dim 1 (= ", dims(thresholds)[1], "), `threshold_effect` dim 2 (= ", dims(threshold_effect)[2], ").");
array[n] int rv;
for(i in 1:n) {
rv[i] = brm_ordinal_rng(
eta[i],
thresholds,
discrimination[i],
structure,
link,
to_vector(threshold_effect[i, :])
);
}
return rv;
}
int brm_ordinal_rng(
real eta,
vector thresholds,
real discrimination,
int structure,
int link,
vector threshold_effect
) {
int k = dims(thresholds)[1];
if (dims(threshold_effect)[1] != k) reject("brm_ordinal_rng: dim mismatch — `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ") does not match `k` (= ", k, "), inferred from `thresholds` dim 1. `k` sizes: `thresholds` dim 1 (= ", dims(thresholds)[1], "), `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ").");
int K = (k + 1);
int rv = K;
if((structure == 1)) {
real u = uniform_rng(0.0, 1.0);
for(j in 1:k) {
if((rv == K)) {
real z_cumulative = (discrimination * (thresholds[j] - eta));
if((u <= brm_ordinal_cdf(z_cumulative | link))) {
rv += (j - rv);
}
}
}
} else {
for(j in 1:k) {
if((rv == K)) {
real z_stopping = (discrimination * ((thresholds[j] - eta) - threshold_effect[j]));
if((bernoulli_rng(brm_ordinal_cdf(z_stopping | link)) == 1)) {
rv += (j - rv);
}
}
}
}
return rv;
}
real brm_ordinal_cdf(
real z,
int link
) {
if((link == 1)) {
return inv_logit(z);
} else {
if((link == 2)) {
return Phi(z);
} else {
return (-expm1((-exp(z))));
}
}
}
}
data {
int x_n;
vector[x_n] x;
int group_n_levels;
int group_idx_n;
array[group_idx_n] int group_idx;
int y_n;
array[y_n] int y;
}
transformed data {
matrix[x_n, 1] X_eta = hcat(x);
int pop_eta_n_covariates = 1;
matrix[num_elements(y), 2] y_threshold_effect = rep_matrix(0.0, num_elements(y), 2);
}
parameters {
vector[pop_eta_n_covariates] pop_eta_beta_pop;
vector[(group_n_levels - 1)] cat_log_disc_group_beta;
ordered[2] y_thresholds;
}
transformed parameters {
vector[x_n] pop_eta = (X_eta * pop_eta_beta_pop);
vector[x_n] eta = pop_eta;
vector[group_idx_n] cat_log_disc_group = append_row(0.0, cat_log_disc_group_beta)[group_idx];
vector[group_idx_n] log_disc = cat_log_disc_group;
vector[group_idx_n] disc = exp(log_disc);
}
model {
pop_eta_beta_pop ~ std_normal();
cat_log_disc_group_beta ~ std_normal();
y_thresholds ~ std_normal();
y ~ brm_ordinal(eta, y_thresholds, disc, 1, 2, y_threshold_effect);
}
generated quantities {
vector[num_elements(y)] y_likelihood = brm_ordinal_lpmfs(y, eta, y_thresholds, disc, 1, 2, y_threshold_effect);
array[num_elements(y)] int y_gen = brm_ordinal_int_rng(y_n, eta, y_thresholds, disc, 1, 2, y_threshold_effect);
}julia
#= line 0 =# Turing.@model(function brm_model(y, X_eta, X_disc, callable_1)
beta_pop ~ Distributions.product_distribution(Base.fill(Distributions.Normal(), 1))
eta_eta = X_eta * beta_pop
eta = eta_eta
beta_pop_disc ~ Distributions.product_distribution(Base.fill(Distributions.Normal(), 1))
eta_disc = X_disc * beta_pop_disc
disc = Base.exp.(eta_disc)
y_thresholds ~ callable_1(2)
begin
for i = Base.eachindex(y)
y[i] ~ BayesianRegressionModels.Ordinal($(QuoteNode(Cumulative())), $(QuoteNode(ProbitLink())), eta[i], y_thresholds; discrimination = disc[i])
end
end
(; disc = disc, eta = eta, y_thresholds = y_thresholds, response = y)
end)discrimination defaults to one. Literal or data-supplied values are checked for finiteness and strict positivity; a modeled value should use a positive-support prior or a link such as log(disc) ~ ....
Stopping-ratio models may add non-proportional effects with a tuple of raw numeric predictors:
brm-comparison
Sequential ordinal modeljulia
sequential_data = (;
period=[1, 2, 3, 1, 2, 3], carry=[0, 0, 1, 0, 1, 1],
treat=[1, 2, 1, 2, 1, 2], y=[1, 1, 2, 2, 3, 3],
)
sequential = @brm sequential_data begin
eta ~ 0 + period + carry
y ~ Ordinal(StoppingRatio(), CloglogLink(), eta;
per_threshold=(treat,))
endjulia
BRMI:
period: data (eltype=Int64, n=6)
carry: data (eltype=Int64, n=6)
eta ~ 0 + period + carry
treat: data (eltype=Int64, n=6)
y ~ Ordinal(StoppingRatio(), CloglogLink(), eta; per_threshold=(treat,))julia
SBBRMI with data keys = [:carry, :carry_idx, :carry_n_levels, :period, :period_idx, :period_n_levels, :treat, :y]
emitted @slic body:
begin
cat_eta_period ~ _sb_cat(; x = period_idx, n_levels = period_n_levels)
cat_eta_carry ~ _sb_cat(; x = carry_idx, n_levels = carry_n_levels)
eta = cat_eta_period + cat_eta_carry
y_thresholds::vector[2] ~ std_normal()
y_threshold_X = hcat(treat)
y_threshold_beta::vector[2, 1] ~ multi_std_normal()
y_threshold_effect = y_threshold_X * adjoint(ranef_b_matrix(y_threshold_beta))
y ~ brm_ordinal(eta, y_thresholds, 1.0, 2, 3, y_threshold_effect)
endstan
functions {
matrix hcat(vector x) {
int n = dims(x)[1];
return to_matrix(x, n, 1);
}
real multi_std_normal_lpdf(
array[] vector x
) {
int m = dims(x)[1];
real rv = 0.0;
for(i in 1:m) {
rv += std_normal_lpdf(x[i, :]);
}
return rv;
}
matrix ranef_b_matrix(
array[] vector b
) {
int m = dims(b)[1];
int n = dims(b)[2];
matrix[m, n] rv;
for(i in 1:m) {
rv[i, :] = (b[i]');
}
return rv;
}
real brm_ordinal_lpmf(
array[] int y,
vector eta,
vector thresholds,
real discrimination,
int structure,
int link,
matrix threshold_effect
) {
int n = dims(y)[1];
int k = dims(thresholds)[1];
if (dims(eta)[1] != n) reject("brm_ordinal_lpmf: dim mismatch — `eta` dim 1 (= ", dims(eta)[1], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `eta` dim 1 (= ", dims(eta)[1], "), `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ").");
if (dims(threshold_effect)[1] != n) reject("brm_ordinal_lpmf: dim mismatch — `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `eta` dim 1 (= ", dims(eta)[1], "), `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ").");
if (dims(threshold_effect)[2] != k) reject("brm_ordinal_lpmf: dim mismatch — `threshold_effect` dim 2 (= ", dims(threshold_effect)[2], ") does not match `k` (= ", k, "), inferred from `thresholds` dim 1. `k` sizes: `thresholds` dim 1 (= ", dims(thresholds)[1], "), `threshold_effect` dim 2 (= ", dims(threshold_effect)[2], ").");
return brm_ordinal_lpmf(y |
eta,
thresholds,
rep_vector(discrimination, n),
structure,
link,
threshold_effect
);
}
real brm_ordinal_lpmf(
array[] int y,
vector eta,
vector thresholds,
vector discrimination,
int structure,
int link,
matrix threshold_effect
) {
int n = dims(y)[1];
int k = dims(thresholds)[1];
if (dims(eta)[1] != n) reject("brm_ordinal_lpmf: dim mismatch — `eta` dim 1 (= ", dims(eta)[1], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `eta` dim 1 (= ", dims(eta)[1], "), `discrimination` dim 1 (= ", dims(discrimination)[1], "), `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ").");
if (dims(discrimination)[1] != n) reject("brm_ordinal_lpmf: dim mismatch — `discrimination` dim 1 (= ", dims(discrimination)[1], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `eta` dim 1 (= ", dims(eta)[1], "), `discrimination` dim 1 (= ", dims(discrimination)[1], "), `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ").");
if (dims(threshold_effect)[1] != n) reject("brm_ordinal_lpmf: dim mismatch — `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `eta` dim 1 (= ", dims(eta)[1], "), `discrimination` dim 1 (= ", dims(discrimination)[1], "), `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ").");
if (dims(threshold_effect)[2] != k) reject("brm_ordinal_lpmf: dim mismatch — `threshold_effect` dim 2 (= ", dims(threshold_effect)[2], ") does not match `k` (= ", k, "), inferred from `thresholds` dim 1. `k` sizes: `thresholds` dim 1 (= ", dims(thresholds)[1], "), `threshold_effect` dim 2 (= ", dims(threshold_effect)[2], ").");
real rv = 0.0;
for(i in 1:n) {
rv += brm_ordinal_lpmf(y[i] |
eta[i],
thresholds,
discrimination[i],
structure,
link,
to_vector(threshold_effect[i, :])
);
}
return rv;
}
real brm_ordinal_lpmf(
int y,
real eta,
vector thresholds,
real discrimination,
int structure,
int link,
vector threshold_effect
) {
int k = dims(thresholds)[1];
if (dims(threshold_effect)[1] != k) reject("brm_ordinal_lpmf: dim mismatch — `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ") does not match `k` (= ", k, "), inferred from `thresholds` dim 1. `k` sizes: `thresholds` dim 1 (= ", dims(thresholds)[1], "), `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ").");
int K = (k + 1);
if((discrimination <= 0.0)) {
return negative_infinity();
} else {
if((y < 1)) {
return negative_infinity();
} else {
if((y > K)) {
return negative_infinity();
} else {
if((structure == 1)) {
if((link == 1)) {
return ordered_logistic_lpmf(y | (discrimination * eta), (discrimination .* thresholds));
} else {
if((y == 1)) {
real z_first = (discrimination * (thresholds[1] - eta));
return brm_ordinal_logcdf(z_first, link);
} else {
if((y == K)) {
real z_last = (discrimination * (thresholds[k] - eta));
return brm_ordinal_logccdf(z_last, link);
} else {
real z_hi = (discrimination * (thresholds[y] - eta));
real z_lo = (discrimination * (thresholds[(y - 1)] - eta));
return log_diff_exp(brm_ordinal_logcdf(z_hi, link), brm_ordinal_logcdf(z_lo, link));
}
}
}
} else {
real rv = 0.0;
for(j in 1:k) {
real z_stage = (discrimination * ((thresholds[j] - eta) - threshold_effect[j]));
if((j < y)) {
rv += brm_ordinal_logccdf(z_stage, link);
} else {
if((j == y)) {
rv += brm_ordinal_logcdf(z_stage, link);
}
}
}
return rv;
}
}
}
}
}
real brm_ordinal_logcdf(
real z,
int link
) {
if((link == 1)) {
return log_inv_logit(z);
} else {
if((link == 2)) {
return normal_lcdf(z | 0.0, 1.0);
} else {
return log1m_exp((-exp(z)));
}
}
}
real brm_ordinal_logccdf(
real z,
int link
) {
if((link == 1)) {
return log_inv_logit((-z));
} else {
if((link == 2)) {
return normal_lccdf(z | 0.0, 1.0);
} else {
return (-exp(z));
}
}
}
vector brm_ordinal_lpmfs(
array[] int y,
vector eta,
vector thresholds,
real discrimination,
int structure,
int link,
matrix threshold_effect
) {
int n = dims(y)[1];
int k = dims(thresholds)[1];
if (dims(eta)[1] != n) reject("brm_ordinal_lpmfs: dim mismatch — `eta` dim 1 (= ", dims(eta)[1], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `eta` dim 1 (= ", dims(eta)[1], "), `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ").");
if (dims(threshold_effect)[1] != n) reject("brm_ordinal_lpmfs: dim mismatch — `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `eta` dim 1 (= ", dims(eta)[1], "), `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ").");
if (dims(threshold_effect)[2] != k) reject("brm_ordinal_lpmfs: dim mismatch — `threshold_effect` dim 2 (= ", dims(threshold_effect)[2], ") does not match `k` (= ", k, "), inferred from `thresholds` dim 1. `k` sizes: `thresholds` dim 1 (= ", dims(thresholds)[1], "), `threshold_effect` dim 2 (= ", dims(threshold_effect)[2], ").");
return brm_ordinal_lpmfs(
y,
eta,
thresholds,
rep_vector(discrimination, n),
structure,
link,
threshold_effect
);
}
vector brm_ordinal_lpmfs(
array[] int y,
vector eta,
vector thresholds,
vector discrimination,
int structure,
int link,
matrix threshold_effect
) {
int n = dims(y)[1];
int k = dims(thresholds)[1];
if (dims(eta)[1] != n) reject("brm_ordinal_lpmfs: dim mismatch — `eta` dim 1 (= ", dims(eta)[1], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `eta` dim 1 (= ", dims(eta)[1], "), `discrimination` dim 1 (= ", dims(discrimination)[1], "), `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ").");
if (dims(discrimination)[1] != n) reject("brm_ordinal_lpmfs: dim mismatch — `discrimination` dim 1 (= ", dims(discrimination)[1], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `eta` dim 1 (= ", dims(eta)[1], "), `discrimination` dim 1 (= ", dims(discrimination)[1], "), `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ").");
if (dims(threshold_effect)[1] != n) reject("brm_ordinal_lpmfs: dim mismatch — `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `eta` dim 1 (= ", dims(eta)[1], "), `discrimination` dim 1 (= ", dims(discrimination)[1], "), `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ").");
if (dims(threshold_effect)[2] != k) reject("brm_ordinal_lpmfs: dim mismatch — `threshold_effect` dim 2 (= ", dims(threshold_effect)[2], ") does not match `k` (= ", k, "), inferred from `thresholds` dim 1. `k` sizes: `thresholds` dim 1 (= ", dims(thresholds)[1], "), `threshold_effect` dim 2 (= ", dims(threshold_effect)[2], ").");
vector[n] rv;
for(i in 1:n) {
rv[i] = brm_ordinal_lpmf(y[i] |
eta[i],
thresholds,
discrimination[i],
structure,
link,
to_vector(threshold_effect[i, :])
);
}
return rv;
}
array[] int brm_ordinal_int_rng(
int anontok__1,
vector eta,
vector thresholds,
real discrimination,
int structure,
int link,
matrix threshold_effect
) {
int n = anontok__1;
int k = dims(thresholds)[1];
if (dims(eta)[1] != n) reject("brm_ordinal_rng: dim mismatch — `eta` dim 1 (= ", dims(eta)[1], ") does not match `n` (= ", n, "), inferred from `anontok__1` dim 1. `n` sizes: `anontok__1` dim 1 (= ", anontok__1, "), `eta` dim 1 (= ", dims(eta)[1], "), `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ").");
if (dims(threshold_effect)[1] != n) reject("brm_ordinal_rng: dim mismatch — `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ") does not match `n` (= ", n, "), inferred from `anontok__1` dim 1. `n` sizes: `anontok__1` dim 1 (= ", anontok__1, "), `eta` dim 1 (= ", dims(eta)[1], "), `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ").");
if (dims(threshold_effect)[2] != k) reject("brm_ordinal_rng: dim mismatch — `threshold_effect` dim 2 (= ", dims(threshold_effect)[2], ") does not match `k` (= ", k, "), inferred from `thresholds` dim 1. `k` sizes: `thresholds` dim 1 (= ", dims(thresholds)[1], "), `threshold_effect` dim 2 (= ", dims(threshold_effect)[2], ").");
return brm_ordinal_int_rng(
n,
eta,
thresholds,
rep_vector(discrimination, n),
structure,
link,
threshold_effect
);
}
array[] int brm_ordinal_int_rng(
int anontok__1,
vector eta,
vector thresholds,
vector discrimination,
int structure,
int link,
matrix threshold_effect
) {
int n = anontok__1;
int k = dims(thresholds)[1];
if (dims(eta)[1] != n) reject("brm_ordinal_rng: dim mismatch — `eta` dim 1 (= ", dims(eta)[1], ") does not match `n` (= ", n, "), inferred from `anontok__1` dim 1. `n` sizes: `anontok__1` dim 1 (= ", anontok__1, "), `eta` dim 1 (= ", dims(eta)[1], "), `discrimination` dim 1 (= ", dims(discrimination)[1], "), `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ").");
if (dims(discrimination)[1] != n) reject("brm_ordinal_rng: dim mismatch — `discrimination` dim 1 (= ", dims(discrimination)[1], ") does not match `n` (= ", n, "), inferred from `anontok__1` dim 1. `n` sizes: `anontok__1` dim 1 (= ", anontok__1, "), `eta` dim 1 (= ", dims(eta)[1], "), `discrimination` dim 1 (= ", dims(discrimination)[1], "), `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ").");
if (dims(threshold_effect)[1] != n) reject("brm_ordinal_rng: dim mismatch — `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ") does not match `n` (= ", n, "), inferred from `anontok__1` dim 1. `n` sizes: `anontok__1` dim 1 (= ", anontok__1, "), `eta` dim 1 (= ", dims(eta)[1], "), `discrimination` dim 1 (= ", dims(discrimination)[1], "), `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ").");
if (dims(threshold_effect)[2] != k) reject("brm_ordinal_rng: dim mismatch — `threshold_effect` dim 2 (= ", dims(threshold_effect)[2], ") does not match `k` (= ", k, "), inferred from `thresholds` dim 1. `k` sizes: `thresholds` dim 1 (= ", dims(thresholds)[1], "), `threshold_effect` dim 2 (= ", dims(threshold_effect)[2], ").");
array[n] int rv;
for(i in 1:n) {
rv[i] = brm_ordinal_rng(
eta[i],
thresholds,
discrimination[i],
structure,
link,
to_vector(threshold_effect[i, :])
);
}
return rv;
}
int brm_ordinal_rng(
real eta,
vector thresholds,
real discrimination,
int structure,
int link,
vector threshold_effect
) {
int k = dims(thresholds)[1];
if (dims(threshold_effect)[1] != k) reject("brm_ordinal_rng: dim mismatch — `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ") does not match `k` (= ", k, "), inferred from `thresholds` dim 1. `k` sizes: `thresholds` dim 1 (= ", dims(thresholds)[1], "), `threshold_effect` dim 1 (= ", dims(threshold_effect)[1], ").");
int K = (k + 1);
int rv = K;
if((structure == 1)) {
real u = uniform_rng(0.0, 1.0);
for(j in 1:k) {
if((rv == K)) {
real z_cumulative = (discrimination * (thresholds[j] - eta));
if((u <= brm_ordinal_cdf(z_cumulative | link))) {
rv += (j - rv);
}
}
}
} else {
for(j in 1:k) {
if((rv == K)) {
real z_stopping = (discrimination * ((thresholds[j] - eta) - threshold_effect[j]));
if((bernoulli_rng(brm_ordinal_cdf(z_stopping | link)) == 1)) {
rv += (j - rv);
}
}
}
}
return rv;
}
real brm_ordinal_cdf(
real z,
int link
) {
if((link == 1)) {
return inv_logit(z);
} else {
if((link == 2)) {
return Phi(z);
} else {
return (-expm1((-exp(z))));
}
}
}
}
data {
int period_n_levels;
int period_idx_n;
array[period_idx_n] int period_idx;
int carry_n_levels;
int carry_idx_n;
array[carry_idx_n] int carry_idx;
int treat_n;
vector[treat_n] treat;
int y_n;
array[y_n] int y;
}
transformed data {
matrix[treat_n, 1] y_threshold_X = hcat(treat);
}
parameters {
vector[(period_n_levels - 1)] cat_eta_period_beta;
vector[(carry_n_levels - 1)] cat_eta_carry_beta;
vector[2] y_thresholds;
array[2] vector[1] y_threshold_beta;
}
transformed parameters {
vector[period_idx_n] cat_eta_period = append_row(0.0, cat_eta_period_beta)[period_idx];
vector[carry_idx_n] cat_eta_carry = append_row(0.0, cat_eta_carry_beta)[carry_idx];
vector[period_idx_n] eta = (cat_eta_period + cat_eta_carry);
matrix[treat_n, 2] y_threshold_effect = (y_threshold_X * (ranef_b_matrix(y_threshold_beta)'));
}
model {
cat_eta_period_beta ~ std_normal();
cat_eta_carry_beta ~ std_normal();
y_thresholds ~ std_normal();
y_threshold_beta ~ multi_std_normal();
y ~ brm_ordinal(eta, y_thresholds, 1.0, 2, 3, y_threshold_effect);
}
generated quantities {
vector[treat_n] y_likelihood = brm_ordinal_lpmfs(y, eta, y_thresholds, 1.0, 2, 3, y_threshold_effect);
array[treat_n] int y_gen = brm_ordinal_int_rng(y_n, eta, y_thresholds, 1.0, 2, 3, y_threshold_effect);
}julia
#= line 0 =# Turing.@model(function brm_model(y, X_eta, callable_1, callable_2, treat)
beta_pop ~ Distributions.product_distribution(Base.fill(Distributions.Normal(), 3))
eta_eta = X_eta * beta_pop
eta = eta_eta
y_threshold_beta ~ callable_1(2)
y_thresholds ~ callable_2(2)
begin
for i = Base.eachindex(y)
y[i] ~ BayesianRegressionModels.Ordinal($(QuoteNode(StoppingRatio())), $(QuoteNode(CloglogLink())), BayesianRegressionModels._brm_threshold_eta(eta[i], (treat[i],), y_threshold_beta, 2), y_thresholds; discrimination = 1.0)
end
end
(; eta = eta, y_threshold_beta = y_threshold_beta, y_thresholds = y_thresholds, response = y)
end)BRM estimates one coefficient per predictor and non-terminal stage, so here per_threshold is deliberately restricted to stopping-ratio models for now: unrestricted cumulative category-specific effects can make cumulative probabilities non-monotone. The predictors must currently be raw numeric data columns.
Outcome categories follow the declared order of a CategoricalVector; plain vectors use sorted unique values. That fitted order is frozen for replay and prediction. The legacy OrderedLogistic(eta) spelling remains supported and continues to lower directly to Stan's native ordered-logistic distribution. The composed cumulative-logit kernel also delegates its scalar density to that native primitive; the other links use Stan's native stable CDF/log-CDF functions. Stopping ratio has no native Stan distribution, so BRM supplies the matching stable lpmf, pointwise log-likelihood, and RNG.
Outside @brm, Ordinal(structure, link, eta, thresholds; discrimination=1) is an executable DiscreteUnivariateDistribution with params, probs, logpdf, and rand. A stopping-ratio eta may be scalar or a vector with one stage-specific value per threshold.
For neutral comparison, brms exposes the same statistical axes through families such as cumulative-probit and stopping-ratio complementary-log-log, and calls threshold-varying terms category-specific effects. BRM preserves that statistical contract while using the typed composition and per_threshold=(...) tuple above rather than importing brms's R formula helpers.
Median regression with Laplace
The StanBlocks backend accepts Distributions.Laplace as an ordinary likelihood. For example, a robust regression for the conditional median can be written as:
brm-comparison
Median regression with Laplacejulia
laplace_data = (;
x = [-1.0, -0.5, 0.0, 0.5, 1.0],
y = [-1.1, -0.2, 0.1, 0.6, 1.4],
)
median_model = @brm laplace_data begin
median_y ~ 1 + x
log(laplace_scale) ~ 1
y ~ Laplace(median_y, laplace_scale)
endjulia
BRMI:
x: data (eltype=Float64, n=5)
median_y ~ 1 + x
log(laplace_scale) ~ 1
y ~ Laplace(median_y, laplace_scale)julia
SBBRMI with data keys = [:x, :y]
emitted @slic body:
begin
X_median_y = hcat(rep_vector(1.0, num_elements(x)), x)
pop_median_y ~ popefs(; X = X_median_y)
median_y = pop_median_y
X_log_laplace_scale = hcat(rep_vector(1.0, num_elements(y)))
pop_log_laplace_scale ~ popefs(; X = X_log_laplace_scale)
log_laplace_scale = pop_log_laplace_scale
laplace_scale = exp(log_laplace_scale)
y ~ double_exponential(median_y, laplace_scale)
endstan
functions {
matrix hcat(
vector x,
vector y
) {
int n = dims(x)[1];
if (dims(y)[1] != n) reject("hcat: dim mismatch — `y` dim 1 (= ", dims(y)[1], ") does not match `n` (= ", n, "), inferred from `x` dim 1. `n` sizes: `x` dim 1 (= ", dims(x)[1], "), `y` dim 1 (= ", dims(y)[1], ").");
return append_col(x, y);
}
matrix hcat(vector x) {
int n = dims(x)[1];
return to_matrix(x, n, 1);
}
vector double_exponential_lpdfs(
vector obs,
vector mu,
vector sigma
) {
return jbroadcasted_double_exponential_lpdfs(obs, mu, sigma);
}
vector jbroadcasted_double_exponential_lpdfs(
vector x1,
vector x2,
vector x3
) {
int n = dims(x1)[1];
vector[n] rv;
for(i in 1:n) {
rv[i] = double_exponential_lpdfs(
broadcasted_getindex(x1, i),
broadcasted_getindex(x2, i),
broadcasted_getindex(x3, i)
);
}
return rv;
}
real double_exponential_lpdfs(
real args1,
real args2,
real args3
) {
return double_exponential_lpdf(args1 | args2, args3);
}
real broadcasted_getindex(vector x, int i) {
return x[i];
}
vector double_exponential_vector_rng(
int anontok__1,
vector a,
vector b
) {
int n = anontok__1;
if((n == 0)) {
vector[n] rv;
return rv;
} else {
return to_vector(double_exponential_rng(a, b));
}
}
}
data {
int x_n;
vector[x_n] x;
int y_n;
vector[y_n] y;
}
transformed data {
matrix[x_n, 2] X_median_y = hcat(rep_vector(1.0, num_elements(x)), x);
int pop_median_y_n_covariates = 2;
matrix[num_elements(y), 1] X_log_laplace_scale = hcat(rep_vector(1.0, num_elements(y)));
int pop_log_laplace_scale_n_covariates = 1;
}
parameters {
vector[pop_median_y_n_covariates] pop_median_y_beta_pop;
vector[pop_log_laplace_scale_n_covariates] pop_log_laplace_scale_beta_pop;
}
transformed parameters {
vector[x_n] pop_median_y = (X_median_y * pop_median_y_beta_pop);
vector[x_n] median_y = pop_median_y;
vector[num_elements(y)] pop_log_laplace_scale = (X_log_laplace_scale * pop_log_laplace_scale_beta_pop);
vector[num_elements(y)] log_laplace_scale = pop_log_laplace_scale;
vector[num_elements(y)] laplace_scale = exp(log_laplace_scale);
}
model {
pop_median_y_beta_pop ~ std_normal();
pop_log_laplace_scale_beta_pop ~ std_normal();
y ~ double_exponential(median_y, laplace_scale);
}
generated quantities {
vector[y_n] y_likelihood = double_exponential_lpdfs(y, median_y, laplace_scale);
vector[y_n] y_gen = double_exponential_vector_rng(y_n, median_y, laplace_scale);
}julia
#= line 0 =# Turing.@model(function brm_model(y, X_median_y, X_laplace_scale)
beta_pop ~ Distributions.product_distribution(Base.fill(Distributions.Normal(), 2))
eta_median_y = X_median_y * beta_pop
median_y = eta_median_y
beta_pop_laplace_scale ~ Distributions.product_distribution(Base.fill(Distributions.Normal(), 1))
eta_laplace_scale = X_laplace_scale * beta_pop_laplace_scale
laplace_scale = Base.exp.(eta_laplace_scale)
begin
for i = Base.eachindex(y)
y[i] ~ Distributions.Laplace(median_y[i], laplace_scale[i])
end
end
(; median_y = median_y, laplace_scale = laplace_scale, response = y)
end)This lowers to Stan's native y ~ double_exponential(median_y, laplace_scale). The second argument is a Laplace scale, not a standard deviation or rate:
This symmetric likelihood is exactly the
In the check-loss convention used by
brms, with scaleand density , use at . In the Bambi/PyMC convention
AsymmetricLaplace(mu, b, kappa), where, use and at .
The Laplace spelling covers median regression only. It does not express an asymmetric likelihood for bs(age, knots=...) term: that basis mapping is a separate unresolved formula-semantic question. Thus the response-family component of the catalogue's quantile_p50 model is available, while the complete historical model remains unsupported.
Quantile regression with SkewDoubleExponential
For a non-median quantile, BRM exposes the executable distribution SkewDoubleExponential(mu, sigma, tau). Its arguments and scale exactly match Stan's native skew_double_exponential family:
Thus cdf(SkewDoubleExponential(mu, sigma, tau), mu) == tau, and SkewDoubleExponential(mu, sigma, 0.5) is exactly Laplace(mu, sigma). There is no hidden brms-scale conversion on this primary Julia surface.
brm-comparison
Non-median quantile regressionjulia
quantile_data = (;
x=[-1.0, -0.5, 0.0, 0.5, 1.0],
y=[-1.1, -0.2, 0.1, 0.6, 1.4],
)
quantile_model = @brm quantile_data begin
q25_y ~ 1 + x
log(native_scale) ~ 1
y ~ SkewDoubleExponential(q25_y, native_scale, 0.25)
endjulia
BRMI:
x: data (eltype=Float64, n=5)
q25_y ~ 1 + x
log(native_scale) ~ 1
y ~ SkewDoubleExponential(q25_y, native_scale, 0.25)julia
SBBRMI with data keys = [:x, :y]
emitted @slic body:
begin
X_q25_y = hcat(rep_vector(1.0, num_elements(x)), x)
pop_q25_y ~ popefs(; X = X_q25_y)
q25_y = pop_q25_y
X_log_native_scale = hcat(rep_vector(1.0, num_elements(y)))
pop_log_native_scale ~ popefs(; X = X_log_native_scale)
log_native_scale = pop_log_native_scale
native_scale = exp(log_native_scale)
y ~ skew_double_exponential(q25_y, native_scale, 0.25)
endstan
functions {
matrix hcat(
vector x,
vector y
) {
int n = dims(x)[1];
if (dims(y)[1] != n) reject("hcat: dim mismatch — `y` dim 1 (= ", dims(y)[1], ") does not match `n` (= ", n, "), inferred from `x` dim 1. `n` sizes: `x` dim 1 (= ", dims(x)[1], "), `y` dim 1 (= ", dims(y)[1], ").");
return append_col(x, y);
}
matrix hcat(vector x) {
int n = dims(x)[1];
return to_matrix(x, n, 1);
}
vector skew_double_exponential_lpdfs(
vector obs,
vector mu,
vector sigma,
real tau
) {
return jbroadcasted_skew_double_exponential_lpdfs(obs, mu, sigma, tau);
}
vector jbroadcasted_skew_double_exponential_lpdfs(
vector x1,
vector x2,
vector x3,
real x4
) {
int n = dims(x1)[1];
vector[n] rv;
for(i in 1:n) {
rv[i] = skew_double_exponential_lpdfs(
broadcasted_getindex(x1, i),
broadcasted_getindex(x2, i),
broadcasted_getindex(x3, i),
x4
);
}
return rv;
}
real skew_double_exponential_lpdfs(
real args1,
real args2,
real args3,
real args4
) {
return skew_double_exponential_lpdf(args1 | args2, args3, args4);
}
real broadcasted_getindex(vector x, int i) {
return x[i];
}
vector skew_double_exponential_vector_rng(
int anontok__1,
vector loc,
vector scale,
real b
) {
int n = anontok__1;
if (dims(loc)[1] != n) reject("skew_double_exponential_rng: dim mismatch — `loc` dim 1 (= ", dims(loc)[1], ") does not match `n` (= ", n, "), inferred from `anontok__1` dim 1. `n` sizes: `anontok__1` dim 1 (= ", anontok__1, "), `loc` dim 1 (= ", dims(loc)[1], "), `scale` dim 1 (= ", dims(scale)[1], ").");
if (dims(scale)[1] != n) reject("skew_double_exponential_rng: dim mismatch — `scale` dim 1 (= ", dims(scale)[1], ") does not match `n` (= ", n, "), inferred from `anontok__1` dim 1. `n` sizes: `anontok__1` dim 1 (= ", anontok__1, "), `loc` dim 1 (= ", dims(loc)[1], "), `scale` dim 1 (= ", dims(scale)[1], ").");
if((n == 0)) {
vector[n] rv;
return rv;
} else {
return to_vector(skew_double_exponential_rng(loc, scale, b));
}
}
}
data {
int x_n;
vector[x_n] x;
int y_n;
vector[y_n] y;
}
transformed data {
matrix[x_n, 2] X_q25_y = hcat(rep_vector(1.0, num_elements(x)), x);
int pop_q25_y_n_covariates = 2;
matrix[num_elements(y), 1] X_log_native_scale = hcat(rep_vector(1.0, num_elements(y)));
int pop_log_native_scale_n_covariates = 1;
}
parameters {
vector[pop_q25_y_n_covariates] pop_q25_y_beta_pop;
vector[pop_log_native_scale_n_covariates] pop_log_native_scale_beta_pop;
}
transformed parameters {
vector[x_n] pop_q25_y = (X_q25_y * pop_q25_y_beta_pop);
vector[x_n] q25_y = pop_q25_y;
vector[num_elements(y)] pop_log_native_scale = (X_log_native_scale * pop_log_native_scale_beta_pop);
vector[num_elements(y)] log_native_scale = pop_log_native_scale;
vector[num_elements(y)] native_scale = exp(log_native_scale);
}
model {
pop_q25_y_beta_pop ~ std_normal();
pop_log_native_scale_beta_pop ~ std_normal();
y ~ skew_double_exponential(q25_y, native_scale, 0.25);
}
generated quantities {
vector[y_n] y_likelihood = skew_double_exponential_lpdfs(y, q25_y, native_scale, 0.25);
vector[num_elements(y)] y_gen = skew_double_exponential_vector_rng(y_n, q25_y, native_scale, 0.25);
}julia
#= line 0 =# Turing.@model(function brm_model(y, X_q25_y, X_native_scale)
beta_pop ~ Distributions.product_distribution(Base.fill(Distributions.Normal(), 2))
eta_q25_y = X_q25_y * beta_pop
q25_y = eta_q25_y
beta_pop_native_scale ~ Distributions.product_distribution(Base.fill(Distributions.Normal(), 1))
eta_native_scale = X_native_scale * beta_pop_native_scale
native_scale = Base.exp.(eta_native_scale)
begin
for i = Base.eachindex(y)
y[i] ~ BayesianRegressionModels.SkewDoubleExponential(q25_y[i], native_scale[i], 0.25)
end
end
(; q25_y = q25_y, native_scale = native_scale, response = y)
end)The brms/check-loss scale
julia
using Distributions: SkewedExponentialPower
y ~ SkewedExponentialPower(mu, sigma_sepd, 1, tau)BRM lowers that spelling with 1; other SkewedExponentialPower shapes are rejected because Stan's asymmetric double-exponential family is not a native implementation of the general SEPD. Density, pointwise log likelihood, and predictive RNG all use the same translation.
Circular regression with VonMises
BRM exposes two deliberately different von-Mises likelihoods. Use Distributions.jl's VonMises when its exact Julia contract is intended:
brm-comparison
Moving-support von Misesjulia
von_mises_data = (;
x=[-1.0, -0.5, 0.0, 0.5, 1.0],
direction=[-0.8, -0.2, 0.1, 0.5, 0.9],
)
exact_model = @brm von_mises_data begin
mu ~ 1 + x
log(kappa) ~ 1
direction ~ VonMises(mu, kappa)
endjulia
BRMI:
x: data (eltype=Float64, n=5)
mu ~ 1 + x
log(kappa) ~ 1
direction ~ VonMises(mu, kappa)julia
SBBRMI with data keys = [:direction, :x]
emitted @slic body:
begin
X_mu = hcat(rep_vector(1.0, num_elements(x)), x)
pop_mu ~ popefs(; X = X_mu)
mu = pop_mu
X_log_kappa = hcat(rep_vector(1.0, num_elements(direction)))
pop_log_kappa ~ popefs(; X = X_log_kappa)
log_kappa = pop_log_kappa
kappa = exp(log_kappa)
direction ~ brm_von_mises(mu, kappa, 0.0, 0.0, 0)
endstan
functions {
matrix hcat(
vector x,
vector y
) {
int n = dims(x)[1];
if (dims(y)[1] != n) reject("hcat: dim mismatch — `y` dim 1 (= ", dims(y)[1], ") does not match `n` (= ", n, "), inferred from `x` dim 1. `n` sizes: `x` dim 1 (= ", dims(x)[1], "), `y` dim 1 (= ", dims(y)[1], ").");
return append_col(x, y);
}
matrix hcat(vector x) {
int n = dims(x)[1];
return to_matrix(x, n, 1);
}
real brm_von_mises_lpdf(
vector y,
vector mu,
vector kappa,
real lo,
real hi,
int principal
) {
int n = dims(y)[1];
if (dims(mu)[1] != n) reject("brm_von_mises_lpdf: dim mismatch — `mu` dim 1 (= ", dims(mu)[1], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `mu` dim 1 (= ", dims(mu)[1], "), `kappa` dim 1 (= ", dims(kappa)[1], ").");
if (dims(kappa)[1] != n) reject("brm_von_mises_lpdf: dim mismatch — `kappa` dim 1 (= ", dims(kappa)[1], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `mu` dim 1 (= ", dims(mu)[1], "), `kappa` dim 1 (= ", dims(kappa)[1], ").");
real rv = 0.0;
for(i in 1:n) {
rv += brm_von_mises_lpdf(y[i] | mu[i], kappa[i], lo, hi, principal);
}
return rv;
}
real brm_von_mises_lpdf(
real y,
real mu,
real kappa,
real lo,
real hi,
int principal
) {
if((kappa <= 0.0)) {
return negative_infinity();
} else {
if((principal == 1)) {
if((y < lo)) {
return negative_infinity();
} else {
if((y >= hi)) {
return negative_infinity();
} else {
real wrapped_mu = (lo + fmod(((fmod((mu - lo), (hi - lo)) + hi) - lo), (hi - lo)));
return von_mises_lpdf(y | wrapped_mu, kappa);
}
}
} else {
if((y < (mu - 3.141592653589793))) {
return negative_infinity();
} else {
if((y > (mu + 3.141592653589793))) {
return negative_infinity();
} else {
return von_mises_lpdf(y | mu, kappa);
}
}
}
}
}
vector brm_von_mises_lpdfs(
vector y,
vector mu,
vector kappa,
real lo,
real hi,
int principal
) {
int n = dims(y)[1];
if (dims(mu)[1] != n) reject("brm_von_mises_lpdfs: dim mismatch — `mu` dim 1 (= ", dims(mu)[1], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `mu` dim 1 (= ", dims(mu)[1], "), `kappa` dim 1 (= ", dims(kappa)[1], ").");
if (dims(kappa)[1] != n) reject("brm_von_mises_lpdfs: dim mismatch — `kappa` dim 1 (= ", dims(kappa)[1], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `mu` dim 1 (= ", dims(mu)[1], "), `kappa` dim 1 (= ", dims(kappa)[1], ").");
vector[n] rv;
for(i in 1:n) {
rv[i] = brm_von_mises_lpdf(y[i] | mu[i], kappa[i], lo, hi, principal);
}
return rv;
}
vector brm_von_mises_vector_rng(
int anontok__1,
vector mu,
vector kappa,
real lo,
real hi,
int principal
) {
int n = anontok__1;
if (dims(mu)[1] != n) reject("brm_von_mises_rng: dim mismatch — `mu` dim 1 (= ", dims(mu)[1], ") does not match `n` (= ", n, "), inferred from `anontok__1` dim 1. `n` sizes: `anontok__1` dim 1 (= ", anontok__1, "), `mu` dim 1 (= ", dims(mu)[1], "), `kappa` dim 1 (= ", dims(kappa)[1], ").");
if (dims(kappa)[1] != n) reject("brm_von_mises_rng: dim mismatch — `kappa` dim 1 (= ", dims(kappa)[1], ") does not match `n` (= ", n, "), inferred from `anontok__1` dim 1. `n` sizes: `anontok__1` dim 1 (= ", anontok__1, "), `mu` dim 1 (= ", dims(mu)[1], "), `kappa` dim 1 (= ", dims(kappa)[1], ").");
vector[n] rv;
for(i in 1:n) {
rv[i] = brm_von_mises_rng(mu[i], kappa[i], lo, hi, principal);
}
return rv;
}
real brm_von_mises_rng(
real mu,
real kappa,
real lo,
real hi,
int principal
) {
if((kappa <= 0.0)) {
reject("brm_von_mises_rng: kappa must be strictly positive");
return 0.0;
} else {
real draw = von_mises_rng(mu, kappa);
if((principal == 1)) {
return (lo + fmod(((fmod((draw - lo), (hi - lo)) + hi) - lo), (hi - lo)));
} else {
real support_lo = (mu - 3.141592653589793);
return (
support_lo +
fmod((fmod((draw - support_lo), 6.283185307179586) + 6.283185307179586), 6.283185307179586)
);
}
}
}
}
data {
int x_n;
vector[x_n] x;
int direction_n;
vector[direction_n] direction;
}
transformed data {
matrix[x_n, 2] X_mu = hcat(rep_vector(1.0, num_elements(x)), x);
int pop_mu_n_covariates = 2;
matrix[num_elements(direction), 1] X_log_kappa = hcat(rep_vector(1.0, num_elements(direction)));
int pop_log_kappa_n_covariates = 1;
}
parameters {
vector[pop_mu_n_covariates] pop_mu_beta_pop;
vector[pop_log_kappa_n_covariates] pop_log_kappa_beta_pop;
}
transformed parameters {
vector[x_n] pop_mu = (X_mu * pop_mu_beta_pop);
vector[x_n] mu = pop_mu;
vector[num_elements(direction)] pop_log_kappa = (X_log_kappa * pop_log_kappa_beta_pop);
vector[num_elements(direction)] log_kappa = pop_log_kappa;
vector[num_elements(direction)] kappa = exp(log_kappa);
}
model {
pop_mu_beta_pop ~ std_normal();
pop_log_kappa_beta_pop ~ std_normal();
direction ~ brm_von_mises(mu, kappa, 0.0, 0.0, 0);
}
generated quantities {
vector[num_elements(direction)] direction_likelihood = brm_von_mises_lpdfs(direction, mu, kappa, 0.0, 0.0, 0);
vector[num_elements(direction)] direction_gen = brm_von_mises_vector_rng(direction_n, mu, kappa, 0.0, 0.0, 0);
}julia
#= line 0 =# Turing.@model(function brm_model(y, X_mu, X_kappa)
beta_pop ~ Distributions.product_distribution(Base.fill(Distributions.Normal(), 2))
eta_mu = X_mu * beta_pop
mu = eta_mu
beta_pop_kappa ~ Distributions.product_distribution(Base.fill(Distributions.Normal(), 1))
eta_kappa = X_kappa * beta_pop_kappa
kappa = Base.exp.(eta_kappa)
begin
for i = Base.eachindex(y)
y[i] ~ Distributions.VonMises(mu[i], kappa[i])
end
end
(; mu = mu, kappa = kappa, response = y)
end)This preserves the constructor order (mu, kappa), the shorthand VonMises(kappa) == VonMises(0, kappa), strict kappa > 0, and the moving closed support [mu - pi, mu + pi]. The backend adds those support/domain guards around Stan's native von_mises_lpdf, and recenters native predictive draws onto the same moving interval. Because the support moves with mu, this surface is usually not the right choice for observations encoded once on a fixed principal interval.
For conventional circular regression on a fixed interval, use the distinct BRM distribution CircularVonMises:
brm-comparison
Fixed-interval circular regressionjulia
circular_data = (;
x=[-1.0, -0.5, 0.0, 0.5, 1.0],
direction=[-0.8, -0.2, 0.1, 0.5, 0.9],
)
circular_model = @brm circular_data begin
mu ~ 1 + x
log(kappa) ~ 1
direction ~ CircularVonMises(mu, kappa; interval=(-pi, pi))
endjulia
BRMI:
x: data (eltype=Float64, n=5)
mu ~ 1 + x
log(kappa) ~ 1
pi: MissingColumn()
direction ~ CircularVonMises(mu, kappa; interval=(-(pi), pi))julia
SBBRMI with data keys = [:direction, :x]
emitted @slic body:
begin
X_mu = hcat(rep_vector(1.0, num_elements(x)), x)
pop_mu ~ popefs(; X = X_mu)
mu = pop_mu
X_log_kappa = hcat(rep_vector(1.0, num_elements(direction)))
pop_log_kappa ~ popefs(; X = X_log_kappa)
log_kappa = pop_log_kappa
kappa = exp(log_kappa)
direction ~ brm_von_mises(mu, kappa, -3.141592653589793, 3.141592653589793, 1)
endstan
functions {
matrix hcat(
vector x,
vector y
) {
int n = dims(x)[1];
if (dims(y)[1] != n) reject("hcat: dim mismatch — `y` dim 1 (= ", dims(y)[1], ") does not match `n` (= ", n, "), inferred from `x` dim 1. `n` sizes: `x` dim 1 (= ", dims(x)[1], "), `y` dim 1 (= ", dims(y)[1], ").");
return append_col(x, y);
}
matrix hcat(vector x) {
int n = dims(x)[1];
return to_matrix(x, n, 1);
}
real brm_von_mises_lpdf(
vector y,
vector mu,
vector kappa,
real lo,
real hi,
int principal
) {
int n = dims(y)[1];
if (dims(mu)[1] != n) reject("brm_von_mises_lpdf: dim mismatch — `mu` dim 1 (= ", dims(mu)[1], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `mu` dim 1 (= ", dims(mu)[1], "), `kappa` dim 1 (= ", dims(kappa)[1], ").");
if (dims(kappa)[1] != n) reject("brm_von_mises_lpdf: dim mismatch — `kappa` dim 1 (= ", dims(kappa)[1], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `mu` dim 1 (= ", dims(mu)[1], "), `kappa` dim 1 (= ", dims(kappa)[1], ").");
real rv = 0.0;
for(i in 1:n) {
rv += brm_von_mises_lpdf(y[i] | mu[i], kappa[i], lo, hi, principal);
}
return rv;
}
real brm_von_mises_lpdf(
real y,
real mu,
real kappa,
real lo,
real hi,
int principal
) {
if((kappa <= 0.0)) {
return negative_infinity();
} else {
if((principal == 1)) {
if((y < lo)) {
return negative_infinity();
} else {
if((y >= hi)) {
return negative_infinity();
} else {
real wrapped_mu = (lo + fmod(((fmod((mu - lo), (hi - lo)) + hi) - lo), (hi - lo)));
return von_mises_lpdf(y | wrapped_mu, kappa);
}
}
} else {
if((y < (mu - 3.141592653589793))) {
return negative_infinity();
} else {
if((y > (mu + 3.141592653589793))) {
return negative_infinity();
} else {
return von_mises_lpdf(y | mu, kappa);
}
}
}
}
}
vector brm_von_mises_lpdfs(
vector y,
vector mu,
vector kappa,
real lo,
real hi,
int principal
) {
int n = dims(y)[1];
if (dims(mu)[1] != n) reject("brm_von_mises_lpdfs: dim mismatch — `mu` dim 1 (= ", dims(mu)[1], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `mu` dim 1 (= ", dims(mu)[1], "), `kappa` dim 1 (= ", dims(kappa)[1], ").");
if (dims(kappa)[1] != n) reject("brm_von_mises_lpdfs: dim mismatch — `kappa` dim 1 (= ", dims(kappa)[1], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `mu` dim 1 (= ", dims(mu)[1], "), `kappa` dim 1 (= ", dims(kappa)[1], ").");
vector[n] rv;
for(i in 1:n) {
rv[i] = brm_von_mises_lpdf(y[i] | mu[i], kappa[i], lo, hi, principal);
}
return rv;
}
vector brm_von_mises_vector_rng(
int anontok__1,
vector mu,
vector kappa,
real lo,
real hi,
int principal
) {
int n = anontok__1;
if (dims(mu)[1] != n) reject("brm_von_mises_rng: dim mismatch — `mu` dim 1 (= ", dims(mu)[1], ") does not match `n` (= ", n, "), inferred from `anontok__1` dim 1. `n` sizes: `anontok__1` dim 1 (= ", anontok__1, "), `mu` dim 1 (= ", dims(mu)[1], "), `kappa` dim 1 (= ", dims(kappa)[1], ").");
if (dims(kappa)[1] != n) reject("brm_von_mises_rng: dim mismatch — `kappa` dim 1 (= ", dims(kappa)[1], ") does not match `n` (= ", n, "), inferred from `anontok__1` dim 1. `n` sizes: `anontok__1` dim 1 (= ", anontok__1, "), `mu` dim 1 (= ", dims(mu)[1], "), `kappa` dim 1 (= ", dims(kappa)[1], ").");
vector[n] rv;
for(i in 1:n) {
rv[i] = brm_von_mises_rng(mu[i], kappa[i], lo, hi, principal);
}
return rv;
}
real brm_von_mises_rng(
real mu,
real kappa,
real lo,
real hi,
int principal
) {
if((kappa <= 0.0)) {
reject("brm_von_mises_rng: kappa must be strictly positive");
return 0.0;
} else {
real draw = von_mises_rng(mu, kappa);
if((principal == 1)) {
return (lo + fmod(((fmod((draw - lo), (hi - lo)) + hi) - lo), (hi - lo)));
} else {
real support_lo = (mu - 3.141592653589793);
return (
support_lo +
fmod((fmod((draw - support_lo), 6.283185307179586) + 6.283185307179586), 6.283185307179586)
);
}
}
}
}
data {
int x_n;
vector[x_n] x;
int direction_n;
vector[direction_n] direction;
}
transformed data {
matrix[x_n, 2] X_mu = hcat(rep_vector(1.0, num_elements(x)), x);
int pop_mu_n_covariates = 2;
matrix[num_elements(direction), 1] X_log_kappa = hcat(rep_vector(1.0, num_elements(direction)));
int pop_log_kappa_n_covariates = 1;
}
parameters {
vector[pop_mu_n_covariates] pop_mu_beta_pop;
vector[pop_log_kappa_n_covariates] pop_log_kappa_beta_pop;
}
transformed parameters {
vector[x_n] pop_mu = (X_mu * pop_mu_beta_pop);
vector[x_n] mu = pop_mu;
vector[num_elements(direction)] pop_log_kappa = (X_log_kappa * pop_log_kappa_beta_pop);
vector[num_elements(direction)] log_kappa = pop_log_kappa;
vector[num_elements(direction)] kappa = exp(log_kappa);
}
model {
pop_mu_beta_pop ~ std_normal();
pop_log_kappa_beta_pop ~ std_normal();
direction ~ brm_von_mises(mu, kappa, -3.141592653589793, 3.141592653589793, 1);
}
generated quantities {
vector[num_elements(direction)] direction_likelihood = brm_von_mises_lpdfs(direction, mu, kappa, -3.141592653589793, 3.141592653589793, 1);
vector[num_elements(direction)] direction_gen = brm_von_mises_vector_rng(direction_n, mu, kappa, -3.141592653589793, 3.141592653589793, 1);
}julia
#= line 0 =# Turing.@model(function brm_model(y, X_mu, X_kappa)
beta_pop ~ Distributions.product_distribution(Base.fill(Distributions.Normal(), 2))
eta_mu = X_mu * beta_pop
mu = eta_mu
beta_pop_kappa ~ Distributions.product_distribution(Base.fill(Distributions.Normal(), 1))
eta_kappa = X_kappa * beta_pop_kappa
kappa = Base.exp.(eta_kappa)
begin
for i = Base.eachindex(y)
y[i] ~ BayesianRegressionModels.CircularVonMises(mu[i], kappa[i]; interval = (-pi, pi))
end
end
(; mu = mu, kappa = kappa, response = y)
end)interval is a compile-time pair of finite numbers with length 2pi; it defaults to (-pi, pi). Observations must lie in the half-open interval [lo, hi). BRM wraps mu and generated draws into that interval, while the density itself remains Stan's native von_mises_lpdf. Both arguments are ordinary distributional parameters: BRM supplies no implicit link or prior. Outside a formula, CircularVonMises(mu, kappa; interval=...) is an executable ContinuousUnivariateDistribution: params, logpdf, and rand preserve the same fixed-interval contract used by the Stan lowering.
For comparison, brms uses a fixed (-pi, pi) response convention and makes both mu and kappa distributional, with default tan_half and log links. Those are a useful neutral baseline, but BRM requires links and priors to be written explicitly rather than silently changing Distributions.jl semantics. In particular, Distributions.Gamma takes a scale, whereas Stan/brms gamma syntax takes a rate: the brms prior gamma(2, 0.01) is spelled Gamma(2, 100.0) in a BRM formula.
Typed observation weights
Observation weights live in the @brm model beside the observation distribution:
brm-comparison
Analytic observation weightsjulia
weighted_model = (@brm begin
y ~ weighted(Normal(mu, sigma), aweights(replicate_k))
mu ~ 1 + x
log(sigma) ~ 1
end)((;
x=[-1.0, 0.0, 1.0], y=[-0.2, 0.3, 1.1],
replicate_k=[1.0, 2.0, 4.0],
))julia
BRMI:
mu ~ 1 + x
log(sigma) ~ 1
replicate_k: data (eltype=Float64, n=3)
y ~ weighted(Normal(mu, sigma), aweights(replicate_k))
x: data (eltype=Float64, n=3)julia
SBBRMI with data keys = [:brm_weight_y, :replicate_k, :x, :y]
emitted @slic body:
begin
X_log_sigma = hcat(rep_vector(1.0, num_elements(y)))
pop_log_sigma ~ popefs(; X = X_log_sigma)
log_sigma = pop_log_sigma
sigma = exp(log_sigma)
X_mu = hcat(rep_vector(1.0, num_elements(x)), x)
pop_mu ~ _popefs_coefs(; X = X_mu)
y ~ normal_id_glm(X_mu, 0.0, pop_mu, sigma ./ sqrt(brm_weight_y))
mu = X_mu * pop_mu
endstan
functions {
matrix hcat(vector x) {
int n = dims(x)[1];
return to_matrix(x, n, 1);
}
matrix hcat(
vector x,
vector y
) {
int n = dims(x)[1];
if (dims(y)[1] != n) reject("hcat: dim mismatch — `y` dim 1 (= ", dims(y)[1], ") does not match `n` (= ", n, "), inferred from `x` dim 1. `n` sizes: `x` dim 1 (= ", dims(x)[1], "), `y` dim 1 (= ", dims(y)[1], ").");
return append_col(x, y);
}
vector normal_id_glm_lpdfs(
vector y,
matrix X,
real alpha,
vector beta,
vector 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,
vector 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], ").");
if((m == 0)) {
vector[m] rv;
return rv;
} else {
return normal_id_glm_rng(X, alpha, beta, sigma);
}
}
vector normal_id_glm_rng(
matrix X,
real alpha,
vector beta,
vector sigma
) {
int m = dims(X)[1];
if((m == 0)) {
vector[m] rv;
return rv;
} else {
return to_vector(normal_rng((rep_vector(alpha, m) + (X * beta)), sigma));
}
}
}
data {
int y_n;
vector[y_n] y;
int x_n;
vector[x_n] x;
int brm_weight_y_n;
vector[brm_weight_y_n] brm_weight_y;
}
transformed data {
matrix[num_elements(y), 1] X_log_sigma = hcat(rep_vector(1.0, num_elements(y)));
int pop_log_sigma_n_covariates = 1;
matrix[x_n, 2] X_mu = hcat(rep_vector(1.0, num_elements(x)), x);
int pop_mu_n_covariates = 2;
}
parameters {
vector[pop_log_sigma_n_covariates] pop_log_sigma_beta_pop;
vector[pop_mu_n_covariates] pop_mu_beta_pop;
}
transformed parameters {
vector[num_elements(y)] pop_log_sigma = (X_log_sigma * pop_log_sigma_beta_pop);
vector[num_elements(y)] log_sigma = pop_log_sigma;
vector[num_elements(y)] sigma = exp(log_sigma);
vector[pop_mu_n_covariates] pop_mu = pop_mu_beta_pop;
}
model {
pop_log_sigma_beta_pop ~ std_normal();
pop_mu_beta_pop ~ std_normal();
y ~ normal_id_glm(X_mu, 0.0, pop_mu, (sigma ./ sqrt(brm_weight_y)));
}
generated quantities {
vector[x_n] y_likelihood = normal_id_glm_lpdfs(y, X_mu, 0.0, pop_mu, (sigma ./ sqrt(brm_weight_y)));
vector[x_n] y_gen = normal_id_glm_vector_rng(y_n, X_mu, 0.0, pop_mu, (sigma ./ sqrt(brm_weight_y)));
vector[x_n] mu = (X_mu * pop_mu);
}julia
#= line 0 =# Turing.@model(function brm_model(y, X_sigma, X_mu, weights_y)
beta_pop_sigma ~ Distributions.product_distribution(Base.fill(Distributions.Normal(), 1))
eta_sigma = X_sigma * beta_pop_sigma
sigma = Base.exp.(eta_sigma)
beta_pop ~ Distributions.product_distribution(Base.fill(Distributions.Normal(), 2))
eta_mu = X_mu * beta_pop
mu = eta_mu
begin
for i = Base.eachindex(y)
y[i] ~ Distributions.Normal(mu[i], sigma[i] / Base.sqrt(weights_y[i]))
end
end
(; mu = mu, sigma = sigma, response = y)
end)The StatsBase constructor determines the statistical meaning:
aweights(k)uses analytic/precision semantics. For a Normal response BRM emitsNormal(mu, sigma / sqrt(k)); model density, pointwise likelihood, and predictive draws all use that adjusted scale.fweights(n)uses frequency/repeat semantics. BRM multiplies each model and pointwise log-likelihood contribution byn; predictive draws remain from the original distribution.weights(w)opts into a power likelihood with the same density/pointwise scaling and unchanged predictive distribution.
The current analytic-weight implementation supports Normal observations. Frequency and power weights support likelihood families that lower through BRM's native Distributions.jl-to-Stan family map. Probability weights and unsupported family/type combinations error instead of silently changing meaning. Weight columns are rebuilt from each dataframe by the reusable @brm builder and by reprocess.