PCR sensitivity through time: ten related models

This example ports the model matrix from pcr-sensitivity-vs-time. The binary observations y measure PCR sensitivity over ordered times t. Five structural assumptions for the time curve are each paired with two link-scale conventions, producing ten complete Stan programs from one shared likelihood.

Two likelihood conventions

The function value logit or log is passed as data to the model builder. Julia gives every function its own type, so @deffun can dispatch on ::typeof(logit) and ::typeof(log):

• logit interprets theta as a log-odds value and forwards to bernoulli_logit_lpmf and bernoulli_logit_rng;

• log interprets theta as a log probability and uses bernoulli_lpmf(y, exp(theta)) and bernoulli_rng(exp(theta)).

The log-link parameter must therefore remain non-positive. The heterogeneous and random-walk families achieve that by sampling a positive ordered xi and returning -xi. Under the logit convention, xi is ordered but unbounded. Negation makes the resulting theta sequence descending in both cases.

The five structural families

Family label	Model for the time curve
hetero	One ordered latent value per time point with a common normal(0, 3) density. Apart from the ordering constraint, no temporal smoothing is imposed.
rw1	A first-order Gaussian random walk. sigma controls adjacent changes, so the monotone sequence varies smoothly rather than independently.
rw2	A second-order Gaussian random walk. Its innovation is measured around 2xi[t-1] - xi[t-2], directly penalizing changes in slope.
regression	A two-parameter line alpha + beta * t with non-positive slope. This replaces one latent value per time point with a global trend.
regression_mix	Two independently parameterized regression submodels plus a beta(2, 2) mixing weight. The preserved source applies link_f to both curves, mixes the results, and replaces the shared custom likelihood with bernoulli(theta); the caveat below explains why this historical definition needs review.

The generated-code labels are the Cartesian product of these rows and the logit/log columns. For example, rw2.log means the second-order random walk with the log-probability likelihood, while regression.logit means the linear curve with the log-odds likelihood.

Reuse in the StanBlocks implementation

base_model owns the observed-data likelihood but leaves theta unimplemented. Each structural @slic fragment returns a candidate theta. Base.merge inserts that fragment into the shared base, and mapping over the named tuples bases and link_fs creates the nested posteriors family without duplicating a model body.

The source also exercises several less common StanBlocks features:

• mock_data gives the tracer representative element types and shapes: y is an integer array and t a real vector. It is not embedded as the eventual observed dataset.

• args... defines a variadic pointwise-log-likelihood helper. jbroadcasted applies its scalar method to vectors while forwarding whichever function token and parameter vector follow y.

• @lpxf marks methods usable by sampling syntax. @lhs additionally states the constrained left-hand-side type: ordered[n] for the logit model and positive_ordered[n] for the log model. Those signatures communicate both the parameter transform and the symbolic length n.

• Scalar and sized _rng overloads let the same link-specific definition support a single posterior prediction or an n-element generated array.

• Untyped arguments are deliberately generic. StanBlocks resolves them through its Stan type hierarchy, while the more specific function-token and sized methods keep dispatch unambiguous here.

Historical source caveat

The executable block is preserved for comparison with the original project, including two regression-family lines that should be audited before using these models for inference. It defines upper_alpha(logit) as negative_infinity(), which gives the logit intercept an empty upper range, and regression_mix applies the named link functions directly while combining component curves. The page demonstrates faithful source generation; it does not certify those two historical statistical definitions as a current analysis.

Full Julia source and generated Stan code

The build evaluates the exact displayed Julia source and calls stan_code for every nested member. Choose a family/link label to inspect one complete generated program, or open the side-by-side modal to compare it with the shared Julia definition.

using StanBlocks

import StanBlocks.stan: logit

@deffun begin 
    "Needed for cross validation"
    my_bernoulli_lpmfs(y::int[n], args...) = jbroadcasted(my_bernoulli_lpmfs, y, args...)
    "Needed for cross validation"
    my_bernoulli_lpmfs(y::int, args...) = my_bernoulli_lpmf(y, args...)
    "Needed to compute the joint likelihood"
    @lpxf my_bernoulli_lpmf(y, ::typeof(logit), theta) = bernoulli_logit_lpmf(y, theta)
    "Needed for posterior predictions"
    my_bernoulli_rng(::typeof(logit), theta::real)::int = bernoulli_logit_rng(theta)
    my_bernoulli_rng(int[n], ::typeof(logit), theta::vector[n])::int[n] = bernoulli_logit_rng(theta)
    "Needed to compute the joint likelihood"
    @lpxf my_bernoulli_lpmf(y, ::typeof(log), theta) = bernoulli_lpmf(y, exp(theta))
    "Needed for posterior predictions"
    my_bernoulli_rng(::typeof(log), theta::real)::int = bernoulli_rng(exp(theta))
    my_bernoulli_rng(int[n], ::typeof(log), theta::vector[n])::int[n] = bernoulli_rng(exp(theta))

    "The `Hetero` prior density after constraining - common to both link functions"
    hetero_lpdf(x::vector[n]) = normal_lpdf(x, 0, 3)
    "The `Hetero` prior density after constraining - this function definition gets used to infer the type and shape of the parameter for the logit link"
    @lhs @lpxf hetero_lpdf(x::ordered[n], ::typeof(logit), n) = hetero_lpdf(x)
    "The `Hetero` prior density after constraining - this function definition gets used to infer the type and shape of the parameter for the log link"
    @lhs @lpxf hetero_lpdf(x::positive_ordered[n], ::typeof(log), n) = hetero_lpdf(x)
    "The `RW(1)` prior density after constraining - common to both link functions"
    rw1_lpdf(x::vector[n], sigma) = normal_lpdf(x[1], 0, sigma) + normal_lpdf(x[2:n], x[1:n-1], sigma)
    "The `RW(1)` prior density after constraining - this function definition gets used to infer the type and shape of the parameter for the logit link"
    @lhs @lpxf rw1_lpdf(x::ordered[n], ::typeof(logit), sigma, n) = rw1_lpdf(x, sigma)
    "The `RW(1)` prior density after constraining - this function definition gets used to infer the type and shape of the parameter for the log link"
    @lhs @lpxf rw1_lpdf(x::positive_ordered[n], ::typeof(log), sigma, n) = rw1_lpdf(x, sigma)
    "The `RW(2)` prior density after constraining - common to both link functions"
    rw2_lpdf(x::vector[n], sigma) = rw1_lpdf(x[1:2], sigma) + normal_lpdf(x[3:n], 2x[2:n-1] - x[1:n-2], sigma)
    "The `RW(2)` prior density after constraining - this function definition gets used to infer the type and shape of the parameter for the logit link"
    @lhs @lpxf rw2_lpdf(x::ordered[n], ::typeof(logit), sigma, n) = rw2_lpdf(x, sigma)
    "The `RW(2)` prior density after constraining - this function definition gets used to infer the type and shape of the parameter for the log link"
    @lhs @lpxf rw2_lpdf(x::positive_ordered[n], ::typeof(log), sigma, n) = rw2_lpdf(x, sigma)
    "The upper bound of the alpha parameter (for the logit link function)"
    upper_alpha(::typeof(logit)) = negative_infinity()
    "The upper bound of the alpha parameter (for the log link function)"
    upper_alpha(::typeof(log)) = 0
end
# The main things that StanBlocks.jl has to know is that `y` is an array of `int`s, and `t` is a vector.
mock_data = (;y=[1], t=[1.])

# The base model intentionally leaves theta to the structural fragment merged below.
base_model = @slic mock_data begin 
    n = dims(y)[1]
    "The exact implementation of the likelihood depends on the passed link function `link_f`"
    y ~ my_bernoulli(link_f, theta)
end
# The submodel for the `Hetero` models
centered_hetero = @slic begin 
    "The type of the `xi` parameter depends on the passed link function `link_f`"
    xi ~ hetero(link_f, n)
    "The negation is needed to ensure that `theta` is in descending order"
    return -xi
end
# The submodel for the `RW(1)` models
centered_rw1 = @slic begin 
    sigma ~ std_normal(;lower=0)
    "The type of the `xi` parameter depends on the passed link function `link_f`"
    xi ~ rw1(link_f, sigma, n)
    "The negation is needed to ensure that `theta` is in descending order"
    return -xi
end
# The submodel for the `RW(2)` models
centered_rw2 = @slic begin 
    sigma ~ normal(0, .5; lower=0)
    "The type of the `xi` parameter depends on the passed link function `link_f`"
    xi ~ rw2(link_f, sigma, n)
    "The negation is needed to ensure that `theta` is in descending order"
    return -xi
end
# The submodel for the `regression` models - reused in the `regression_mix` models
regression = @slic begin 
    "The upper bound of alpha `alpha_upper` depends on the link function `link_f`"
    alpha_upper = upper_alpha(link_f)
    # This could be `alpha ~ normal(0, .5; upper=upper_alpha(link_f))`, once https://github.com/nsiccha/StanBlocks.jl/issues/35 is fixed
    alpha ~ normal(0, .5; upper=alpha_upper)
    beta ~ normal(0, .5; upper=0.)
    return alpha + beta * t
end
# The submodel for the `regression_mix` models - reusing the `regression` submodel
regression_mix = @slic begin 
    c1 ~ regression(;link_f, t) 
    c2 ~ regression(;link_f, t)
    lambda ~ beta(2, 2)
    return lambda * link_f(c1) + (1-lambda) * link_f(c2)
end

bases = (;
    hetero=Base.merge(base_model, quote 
        theta ~ centered_hetero(;link_f, n)
    end),
    rw1=Base.merge(base_model, quote 
        theta ~ centered_rw1(;link_f, n)
    end),
    rw2=Base.merge(base_model, quote 
        theta ~ centered_rw2(;link_f, n)
    end),
    regression=Base.merge(base_model, quote 
        theta ~ regression(;link_f, t)
    end),
    regression_mix=Base.merge(base_model, quote 
        theta ~ regression_mix(;link_f, t)
        y ~ bernoulli(theta)
    end)
)
link_fs = (;logit, log)  

# `posteriors` will be a (nested) named tuple - accessing e.g. the `Hetero (log)` model works via `posteriors.hetero.log`
posteriors = map(bases) do base 
    map(link_fs) do link_f 
        base(;link_f)
    end
end

hetero / logit

functions {
// The `Hetero` prior density after constraining - this function definition gets used to infer the type and shape of the parameter for the logit link
real hetero_logit_lpdf(vector x, int n) {
    return hetero_lpdf(x);
}
// The `Hetero` prior density after constraining - common to both link functions
real hetero_lpdf(vector x) {
    return normal_lpdf(x | 0, 3);
}
// Needed to compute the joint likelihood
real my_bernoulli_logit_lpmf(
    array[] int y,
    vector theta
) {
    return bernoulli_logit_lpmf(y | theta);
}
// Needed for cross validation
vector my_bernoulli_logit_lpmfs(
    array[] int y,
    vector args1
) {
    return jbroadcasted_my_bernoulli_lpmfs_logit(y, args1);
}
vector jbroadcasted_my_bernoulli_lpmfs_logit(
    array[] int x1,
    vector x3
) {
    int n = dims(x1)[1];
    vector[n] rv;
    for(i in 1:n) {
        rv[i] = my_bernoulli_logit_lpmfs(broadcasted_getindex(x1, i), broadcasted_getindex(x3, i));
    }
    return rv;
}
// Needed for cross validation
real my_bernoulli_logit_lpmfs(int y, real args1) {
    return my_bernoulli_logit_lpmf(y | args1);
}
// Needed to compute the joint likelihood
real my_bernoulli_logit_lpmf(int y, real theta) {
    return bernoulli_logit_lpmf(y | theta);
}
int broadcasted_getindex(array[] int x, int i) {
    return x[i];
}
real broadcasted_getindex(vector x, int i) {
    return x[i];
}
array[] int my_bernoulli_int_logit_rng(
    int anontok__1,
    vector theta
) {
    int n = anontok__1;
    if (dims(theta)[1] != n) reject("my_bernoulli_rng: dim mismatch — `theta` dim 1 (= ", dims(theta)[1], ") does not match `n` (= ", n, "), inferred from `anontok__1` dim 1. `n` sizes: `anontok__1` dim 1 (= ", anontok__1, "), `theta` dim 1 (= ", dims(theta)[1], ").");
    return bernoulli_logit_rng(theta);
}
}
data {
    int y_n;
    array[y_n] int y;
}
transformed data {
    int n = y_n;
}
parameters {
    ordered[n] theta_xi;
}
transformed parameters {
    // The negation is needed to ensure that `theta` is in descending order
    vector[n] theta = (-theta_xi);
}
model {
    theta_xi ~ hetero_logit(n);
    // The exact implementation of the likelihood depends on the passed link function `link_f`
    y ~ my_bernoulli_logit(theta);
}
generated quantities {
    // The exact implementation of the likelihood depends on the passed link function `link_f`
    vector[y_n] y_likelihood = my_bernoulli_logit_lpmfs(y, theta);
    // The exact implementation of the likelihood depends on the passed link function `link_f`
    array[n] int y_gen = my_bernoulli_int_logit_rng(y_n, theta);
}

hetero / log

functions {
// The `Hetero` prior density after constraining - this function definition gets used to infer the type and shape of the parameter for the log link
real hetero_log_lpdf(vector x, int n) {
    return hetero_lpdf(x);
}
// The `Hetero` prior density after constraining - common to both link functions
real hetero_lpdf(vector x) {
    return normal_lpdf(x | 0, 3);
}
// Needed to compute the joint likelihood
real my_bernoulli_log_lpmf(
    array[] int y,
    vector theta
) {
    return bernoulli_lpmf(y | exp(theta));
}
// Needed for cross validation
vector my_bernoulli_log_lpmfs(
    array[] int y,
    vector args1
) {
    return jbroadcasted_my_bernoulli_lpmfs_log(y, args1);
}
vector jbroadcasted_my_bernoulli_lpmfs_log(
    array[] int x1,
    vector x3
) {
    int n = dims(x1)[1];
    vector[n] rv;
    for(i in 1:n) {
        rv[i] = my_bernoulli_log_lpmfs(broadcasted_getindex(x1, i), broadcasted_getindex(x3, i));
    }
    return rv;
}
// Needed for cross validation
real my_bernoulli_log_lpmfs(int y, real args1) {
    return my_bernoulli_log_lpmf(y | args1);
}
// Needed to compute the joint likelihood
real my_bernoulli_log_lpmf(int y, real theta) {
    return bernoulli_lpmf(y | exp(theta));
}
int broadcasted_getindex(array[] int x, int i) {
    return x[i];
}
real broadcasted_getindex(vector x, int i) {
    return x[i];
}
array[] int my_bernoulli_int_log_rng(
    int anontok__1,
    vector theta
) {
    int n = anontok__1;
    if (dims(theta)[1] != n) reject("my_bernoulli_rng: dim mismatch — `theta` dim 1 (= ", dims(theta)[1], ") does not match `n` (= ", n, "), inferred from `anontok__1` dim 1. `n` sizes: `anontok__1` dim 1 (= ", anontok__1, "), `theta` dim 1 (= ", dims(theta)[1], ").");
    return bernoulli_rng(exp(theta));
}
}
data {
    int y_n;
    array[y_n] int y;
}
transformed data {
    int n = y_n;
}
parameters {
    positive_ordered[n] theta_xi;
}
transformed parameters {
    // The negation is needed to ensure that `theta` is in descending order
    vector[n] theta = (-theta_xi);
}
model {
    theta_xi ~ hetero_log(n);
    // The exact implementation of the likelihood depends on the passed link function `link_f`
    y ~ my_bernoulli_log(theta);
}
generated quantities {
    // The exact implementation of the likelihood depends on the passed link function `link_f`
    vector[y_n] y_likelihood = my_bernoulli_log_lpmfs(y, theta);
    // The exact implementation of the likelihood depends on the passed link function `link_f`
    array[n] int y_gen = my_bernoulli_int_log_rng(y_n, theta);
}

rw1 / logit

functions {
// The `RW(1)` prior density after constraining - this function definition gets used to infer the type and shape of the parameter for the logit link
real rw1_logit_lpdf(vector x, real sigma, int n) {
    return rw1_lpdf(x | sigma);
}
// The `RW(1)` prior density after constraining - common to both link functions
real rw1_lpdf(
    vector x,
    real sigma
) {
    int n = dims(x)[1];
    return (normal_lpdf(x[1] | 0, sigma) + normal_lpdf(x[2:n] | x[1:(n - 1)], sigma));
}
// Needed to compute the joint likelihood
real my_bernoulli_logit_lpmf(
    array[] int y,
    vector theta
) {
    return bernoulli_logit_lpmf(y | theta);
}
// Needed for cross validation
vector my_bernoulli_logit_lpmfs(
    array[] int y,
    vector args1
) {
    return jbroadcasted_my_bernoulli_lpmfs_logit(y, args1);
}
vector jbroadcasted_my_bernoulli_lpmfs_logit(
    array[] int x1,
    vector x3
) {
    int n = dims(x1)[1];
    vector[n] rv;
    for(i in 1:n) {
        rv[i] = my_bernoulli_logit_lpmfs(broadcasted_getindex(x1, i), broadcasted_getindex(x3, i));
    }
    return rv;
}
// Needed for cross validation
real my_bernoulli_logit_lpmfs(int y, real args1) {
    return my_bernoulli_logit_lpmf(y | args1);
}
// Needed to compute the joint likelihood
real my_bernoulli_logit_lpmf(int y, real theta) {
    return bernoulli_logit_lpmf(y | theta);
}
int broadcasted_getindex(array[] int x, int i) {
    return x[i];
}
real broadcasted_getindex(vector x, int i) {
    return x[i];
}
array[] int my_bernoulli_int_logit_rng(
    int anontok__1,
    vector theta
) {
    int n = anontok__1;
    if (dims(theta)[1] != n) reject("my_bernoulli_rng: dim mismatch — `theta` dim 1 (= ", dims(theta)[1], ") does not match `n` (= ", n, "), inferred from `anontok__1` dim 1. `n` sizes: `anontok__1` dim 1 (= ", anontok__1, "), `theta` dim 1 (= ", dims(theta)[1], ").");
    return bernoulli_logit_rng(theta);
}
}
data {
    int y_n;
    array[y_n] int y;
}
transformed data {
    int n = y_n;
}
parameters {
    real<lower=0> theta_sigma;
    // The type of the `xi` parameter depends on the passed link function `link_f`
    ordered[n] theta_xi;
}
transformed parameters {
    // The negation is needed to ensure that `theta` is in descending order
    vector[n] theta = (-theta_xi);
}
model {
    theta_sigma ~ std_normal();
    // The type of the `xi` parameter depends on the passed link function `link_f`
    theta_xi ~ rw1_logit(theta_sigma, n);
    // The exact implementation of the likelihood depends on the passed link function `link_f`
    y ~ my_bernoulli_logit(theta);
}
generated quantities {
    // The exact implementation of the likelihood depends on the passed link function `link_f`
    vector[y_n] y_likelihood = my_bernoulli_logit_lpmfs(y, theta);
    // The exact implementation of the likelihood depends on the passed link function `link_f`
    array[n] int y_gen = my_bernoulli_int_logit_rng(y_n, theta);
}

rw1 / log

functions {
// The `RW(1)` prior density after constraining - this function definition gets used to infer the type and shape of the parameter for the log link
real rw1_log_lpdf(vector x, real sigma, int n) {
    return rw1_lpdf(x | sigma);
}
// The `RW(1)` prior density after constraining - common to both link functions
real rw1_lpdf(
    vector x,
    real sigma
) {
    int n = dims(x)[1];
    return (normal_lpdf(x[1] | 0, sigma) + normal_lpdf(x[2:n] | x[1:(n - 1)], sigma));
}
// Needed to compute the joint likelihood
real my_bernoulli_log_lpmf(
    array[] int y,
    vector theta
) {
    return bernoulli_lpmf(y | exp(theta));
}
// Needed for cross validation
vector my_bernoulli_log_lpmfs(
    array[] int y,
    vector args1
) {
    return jbroadcasted_my_bernoulli_lpmfs_log(y, args1);
}
vector jbroadcasted_my_bernoulli_lpmfs_log(
    array[] int x1,
    vector x3
) {
    int n = dims(x1)[1];
    vector[n] rv;
    for(i in 1:n) {
        rv[i] = my_bernoulli_log_lpmfs(broadcasted_getindex(x1, i), broadcasted_getindex(x3, i));
    }
    return rv;
}
// Needed for cross validation
real my_bernoulli_log_lpmfs(int y, real args1) {
    return my_bernoulli_log_lpmf(y | args1);
}
// Needed to compute the joint likelihood
real my_bernoulli_log_lpmf(int y, real theta) {
    return bernoulli_lpmf(y | exp(theta));
}
int broadcasted_getindex(array[] int x, int i) {
    return x[i];
}
real broadcasted_getindex(vector x, int i) {
    return x[i];
}
array[] int my_bernoulli_int_log_rng(
    int anontok__1,
    vector theta
) {
    int n = anontok__1;
    if (dims(theta)[1] != n) reject("my_bernoulli_rng: dim mismatch — `theta` dim 1 (= ", dims(theta)[1], ") does not match `n` (= ", n, "), inferred from `anontok__1` dim 1. `n` sizes: `anontok__1` dim 1 (= ", anontok__1, "), `theta` dim 1 (= ", dims(theta)[1], ").");
    return bernoulli_rng(exp(theta));
}
}
data {
    int y_n;
    array[y_n] int y;
}
transformed data {
    int n = y_n;
}
parameters {
    real<lower=0> theta_sigma;
    // The type of the `xi` parameter depends on the passed link function `link_f`
    positive_ordered[n] theta_xi;
}
transformed parameters {
    // The negation is needed to ensure that `theta` is in descending order
    vector[n] theta = (-theta_xi);
}
model {
    theta_sigma ~ std_normal();
    // The type of the `xi` parameter depends on the passed link function `link_f`
    theta_xi ~ rw1_log(theta_sigma, n);
    // The exact implementation of the likelihood depends on the passed link function `link_f`
    y ~ my_bernoulli_log(theta);
}
generated quantities {
    // The exact implementation of the likelihood depends on the passed link function `link_f`
    vector[y_n] y_likelihood = my_bernoulli_log_lpmfs(y, theta);
    // The exact implementation of the likelihood depends on the passed link function `link_f`
    array[n] int y_gen = my_bernoulli_int_log_rng(y_n, theta);
}

rw2 / logit

functions {
// The `RW(2)` prior density after constraining - this function definition gets used to infer the type and shape of the parameter for the logit link
real rw2_logit_lpdf(vector x, real sigma, int n) {
    return rw2_lpdf(x | sigma);
}
// The `RW(2)` prior density after constraining - common to both link functions
real rw2_lpdf(
    vector x,
    real sigma
) {
    int n = dims(x)[1];
    return (rw1_lpdf(x[1:2] | sigma) + normal_lpdf(x[3:n] | ((2 * x[2:(n - 1)]) - x[1:(n - 2)]), sigma));
}
// The `RW(1)` prior density after constraining - common to both link functions
real rw1_lpdf(
    vector x,
    real sigma
) {
    int n = dims(x)[1];
    return (normal_lpdf(x[1] | 0, sigma) + normal_lpdf(x[2:n] | x[1:(n - 1)], sigma));
}
// Needed to compute the joint likelihood
real my_bernoulli_logit_lpmf(
    array[] int y,
    vector theta
) {
    return bernoulli_logit_lpmf(y | theta);
}
// Needed for cross validation
vector my_bernoulli_logit_lpmfs(
    array[] int y,
    vector args1
) {
    return jbroadcasted_my_bernoulli_lpmfs_logit(y, args1);
}
vector jbroadcasted_my_bernoulli_lpmfs_logit(
    array[] int x1,
    vector x3
) {
    int n = dims(x1)[1];
    vector[n] rv;
    for(i in 1:n) {
        rv[i] = my_bernoulli_logit_lpmfs(broadcasted_getindex(x1, i), broadcasted_getindex(x3, i));
    }
    return rv;
}
// Needed for cross validation
real my_bernoulli_logit_lpmfs(int y, real args1) {
    return my_bernoulli_logit_lpmf(y | args1);
}
// Needed to compute the joint likelihood
real my_bernoulli_logit_lpmf(int y, real theta) {
    return bernoulli_logit_lpmf(y | theta);
}
int broadcasted_getindex(array[] int x, int i) {
    return x[i];
}
real broadcasted_getindex(vector x, int i) {
    return x[i];
}
array[] int my_bernoulli_int_logit_rng(
    int anontok__1,
    vector theta
) {
    int n = anontok__1;
    if (dims(theta)[1] != n) reject("my_bernoulli_rng: dim mismatch — `theta` dim 1 (= ", dims(theta)[1], ") does not match `n` (= ", n, "), inferred from `anontok__1` dim 1. `n` sizes: `anontok__1` dim 1 (= ", anontok__1, "), `theta` dim 1 (= ", dims(theta)[1], ").");
    return bernoulli_logit_rng(theta);
}
}
data {
    int y_n;
    array[y_n] int y;
}
transformed data {
    int n = y_n;
}
parameters {
    real<lower=0> theta_sigma;
    // The type of the `xi` parameter depends on the passed link function `link_f`
    ordered[n] theta_xi;
}
transformed parameters {
    // The negation is needed to ensure that `theta` is in descending order
    vector[n] theta = (-theta_xi);
}
model {
    theta_sigma ~ normal(0, 0.5);
    // The type of the `xi` parameter depends on the passed link function `link_f`
    theta_xi ~ rw2_logit(theta_sigma, n);
    // The exact implementation of the likelihood depends on the passed link function `link_f`
    y ~ my_bernoulli_logit(theta);
}
generated quantities {
    // The exact implementation of the likelihood depends on the passed link function `link_f`
    vector[y_n] y_likelihood = my_bernoulli_logit_lpmfs(y, theta);
    // The exact implementation of the likelihood depends on the passed link function `link_f`
    array[n] int y_gen = my_bernoulli_int_logit_rng(y_n, theta);
}

rw2 / log

functions {
// The `RW(2)` prior density after constraining - this function definition gets used to infer the type and shape of the parameter for the log link
real rw2_log_lpdf(vector x, real sigma, int n) {
    return rw2_lpdf(x | sigma);
}
// The `RW(2)` prior density after constraining - common to both link functions
real rw2_lpdf(
    vector x,
    real sigma
) {
    int n = dims(x)[1];
    return (rw1_lpdf(x[1:2] | sigma) + normal_lpdf(x[3:n] | ((2 * x[2:(n - 1)]) - x[1:(n - 2)]), sigma));
}
// The `RW(1)` prior density after constraining - common to both link functions
real rw1_lpdf(
    vector x,
    real sigma
) {
    int n = dims(x)[1];
    return (normal_lpdf(x[1] | 0, sigma) + normal_lpdf(x[2:n] | x[1:(n - 1)], sigma));
}
// Needed to compute the joint likelihood
real my_bernoulli_log_lpmf(
    array[] int y,
    vector theta
) {
    return bernoulli_lpmf(y | exp(theta));
}
// Needed for cross validation
vector my_bernoulli_log_lpmfs(
    array[] int y,
    vector args1
) {
    return jbroadcasted_my_bernoulli_lpmfs_log(y, args1);
}
vector jbroadcasted_my_bernoulli_lpmfs_log(
    array[] int x1,
    vector x3
) {
    int n = dims(x1)[1];
    vector[n] rv;
    for(i in 1:n) {
        rv[i] = my_bernoulli_log_lpmfs(broadcasted_getindex(x1, i), broadcasted_getindex(x3, i));
    }
    return rv;
}
// Needed for cross validation
real my_bernoulli_log_lpmfs(int y, real args1) {
    return my_bernoulli_log_lpmf(y | args1);
}
// Needed to compute the joint likelihood
real my_bernoulli_log_lpmf(int y, real theta) {
    return bernoulli_lpmf(y | exp(theta));
}
int broadcasted_getindex(array[] int x, int i) {
    return x[i];
}
real broadcasted_getindex(vector x, int i) {
    return x[i];
}
array[] int my_bernoulli_int_log_rng(
    int anontok__1,
    vector theta
) {
    int n = anontok__1;
    if (dims(theta)[1] != n) reject("my_bernoulli_rng: dim mismatch — `theta` dim 1 (= ", dims(theta)[1], ") does not match `n` (= ", n, "), inferred from `anontok__1` dim 1. `n` sizes: `anontok__1` dim 1 (= ", anontok__1, "), `theta` dim 1 (= ", dims(theta)[1], ").");
    return bernoulli_rng(exp(theta));
}
}
data {
    int y_n;
    array[y_n] int y;
}
transformed data {
    int n = y_n;
}
parameters {
    real<lower=0> theta_sigma;
    // The type of the `xi` parameter depends on the passed link function `link_f`
    positive_ordered[n] theta_xi;
}
transformed parameters {
    // The negation is needed to ensure that `theta` is in descending order
    vector[n] theta = (-theta_xi);
}
model {
    theta_sigma ~ normal(0, 0.5);
    // The type of the `xi` parameter depends on the passed link function `link_f`
    theta_xi ~ rw2_log(theta_sigma, n);
    // The exact implementation of the likelihood depends on the passed link function `link_f`
    y ~ my_bernoulli_log(theta);
}
generated quantities {
    // The exact implementation of the likelihood depends on the passed link function `link_f`
    vector[y_n] y_likelihood = my_bernoulli_log_lpmfs(y, theta);
    // The exact implementation of the likelihood depends on the passed link function `link_f`
    array[n] int y_gen = my_bernoulli_int_log_rng(y_n, theta);
}

regression / logit

functions {
// The upper bound of the alpha parameter (for the logit link function)
real upper_alpha_logit() {
    return negative_infinity();
}
// Needed to compute the joint likelihood
real my_bernoulli_logit_lpmf(
    array[] int y,
    vector theta
) {
    return bernoulli_logit_lpmf(y | theta);
}
// Needed for cross validation
vector my_bernoulli_logit_lpmfs(
    array[] int y,
    vector args1
) {
    return jbroadcasted_my_bernoulli_lpmfs_logit(y, args1);
}
vector jbroadcasted_my_bernoulli_lpmfs_logit(
    array[] int x1,
    vector x3
) {
    int n = dims(x1)[1];
    vector[n] rv;
    for(i in 1:n) {
        rv[i] = my_bernoulli_logit_lpmfs(broadcasted_getindex(x1, i), broadcasted_getindex(x3, i));
    }
    return rv;
}
// Needed for cross validation
real my_bernoulli_logit_lpmfs(int y, real args1) {
    return my_bernoulli_logit_lpmf(y | args1);
}
// Needed to compute the joint likelihood
real my_bernoulli_logit_lpmf(int y, real theta) {
    return bernoulli_logit_lpmf(y | theta);
}
int broadcasted_getindex(array[] int x, int i) {
    return x[i];
}
real broadcasted_getindex(vector x, int i) {
    return x[i];
}
array[] int my_bernoulli_int_logit_rng(
    int anontok__1,
    vector theta
) {
    int n = anontok__1;
    if (dims(theta)[1] != n) reject("my_bernoulli_rng: dim mismatch — `theta` dim 1 (= ", dims(theta)[1], ") does not match `n` (= ", n, "), inferred from `anontok__1` dim 1. `n` sizes: `anontok__1` dim 1 (= ", anontok__1, "), `theta` dim 1 (= ", dims(theta)[1], ").");
    return bernoulli_logit_rng(theta);
}
}
data {
    int y_n;
    int t_n;
    vector[t_n] t;
    array[y_n] int y;
}
transformed data {
    int n = y_n;
    real theta_alpha_upper = upper_alpha_logit();
}
parameters {
    real<upper=theta_alpha_upper> theta_alpha;
    real<upper=0.0> theta_beta;
}
transformed parameters {
    vector[t_n] theta = (theta_alpha + (theta_beta * t));
}
model {
    theta_alpha ~ normal(0, 0.5);
    theta_beta ~ normal(0, 0.5);
    // The exact implementation of the likelihood depends on the passed link function `link_f`
    y ~ my_bernoulli_logit(theta);
}
generated quantities {
    // The exact implementation of the likelihood depends on the passed link function `link_f`
    vector[y_n] y_likelihood = my_bernoulli_logit_lpmfs(y, theta);
    // The exact implementation of the likelihood depends on the passed link function `link_f`
    array[t_n] int y_gen = my_bernoulli_int_logit_rng(y_n, theta);
}

regression / log

functions {
// The upper bound of the alpha parameter (for the log link function)
int upper_alpha_log() {
    return 0;
}
// Needed to compute the joint likelihood
real my_bernoulli_log_lpmf(
    array[] int y,
    vector theta
) {
    return bernoulli_lpmf(y | exp(theta));
}
// Needed for cross validation
vector my_bernoulli_log_lpmfs(
    array[] int y,
    vector args1
) {
    return jbroadcasted_my_bernoulli_lpmfs_log(y, args1);
}
vector jbroadcasted_my_bernoulli_lpmfs_log(
    array[] int x1,
    vector x3
) {
    int n = dims(x1)[1];
    vector[n] rv;
    for(i in 1:n) {
        rv[i] = my_bernoulli_log_lpmfs(broadcasted_getindex(x1, i), broadcasted_getindex(x3, i));
    }
    return rv;
}
// Needed for cross validation
real my_bernoulli_log_lpmfs(int y, real args1) {
    return my_bernoulli_log_lpmf(y | args1);
}
// Needed to compute the joint likelihood
real my_bernoulli_log_lpmf(int y, real theta) {
    return bernoulli_lpmf(y | exp(theta));
}
int broadcasted_getindex(array[] int x, int i) {
    return x[i];
}
real broadcasted_getindex(vector x, int i) {
    return x[i];
}
array[] int my_bernoulli_int_log_rng(
    int anontok__1,
    vector theta
) {
    int n = anontok__1;
    if (dims(theta)[1] != n) reject("my_bernoulli_rng: dim mismatch — `theta` dim 1 (= ", dims(theta)[1], ") does not match `n` (= ", n, "), inferred from `anontok__1` dim 1. `n` sizes: `anontok__1` dim 1 (= ", anontok__1, "), `theta` dim 1 (= ", dims(theta)[1], ").");
    return bernoulli_rng(exp(theta));
}
}
data {
    int y_n;
    int t_n;
    vector[t_n] t;
    array[y_n] int y;
}
transformed data {
    int n = y_n;
    int theta_alpha_upper = upper_alpha_log();
}
parameters {
    real<upper=theta_alpha_upper> theta_alpha;
    real<upper=0.0> theta_beta;
}
transformed parameters {
    vector[t_n] theta = (theta_alpha + (theta_beta * t));
}
model {
    theta_alpha ~ normal(0, 0.5);
    theta_beta ~ normal(0, 0.5);
    // The exact implementation of the likelihood depends on the passed link function `link_f`
    y ~ my_bernoulli_log(theta);
}
generated quantities {
    // The exact implementation of the likelihood depends on the passed link function `link_f`
    vector[y_n] y_likelihood = my_bernoulli_log_lpmfs(y, theta);
    // The exact implementation of the likelihood depends on the passed link function `link_f`
    array[t_n] int y_gen = my_bernoulli_int_log_rng(y_n, theta);
}

regression_mix / logit

functions {
// The upper bound of the alpha parameter (for the logit link function)
real upper_alpha_logit() {
    return negative_infinity();
}
vector bernoulli_lpmfs(
    array[] int obs,
    vector args1
) {
    return jbroadcasted_bernoulli_lpmfs(obs, args1);
}
vector jbroadcasted_bernoulli_lpmfs(
    array[] int x1,
    vector x2
) {
    int n = dims(x1)[1];
    vector[n] rv;
    for(i in 1:n) {
        rv[i] = bernoulli_lpmfs(broadcasted_getindex(x1, i), broadcasted_getindex(x2, i));
    }
    return rv;
}
real bernoulli_lpmfs(int args1, real args2) {
    return bernoulli_lpmf(args1 | args2);
}
int broadcasted_getindex(array[] int x, int i) {
    return x[i];
}
real broadcasted_getindex(vector x, int i) {
    return x[i];
}
array[] int bernoulli_int_rng(
    int anontok__1,
    vector p
) {
    int n = anontok__1;
    return bernoulli_rng(p);
}
// Needed to compute the joint likelihood
real my_bernoulli_logit_lpmf(
    array[] int y,
    vector theta
) {
    return bernoulli_logit_lpmf(y | theta);
}
// Needed for cross validation
vector my_bernoulli_logit_lpmfs(
    array[] int y,
    vector args1
) {
    return jbroadcasted_my_bernoulli_lpmfs_logit(y, args1);
}
vector jbroadcasted_my_bernoulli_lpmfs_logit(
    array[] int x1,
    vector x3
) {
    int n = dims(x1)[1];
    vector[n] rv;
    for(i in 1:n) {
        rv[i] = my_bernoulli_logit_lpmfs(broadcasted_getindex(x1, i), broadcasted_getindex(x3, i));
    }
    return rv;
}
// Needed for cross validation
real my_bernoulli_logit_lpmfs(int y, real args1) {
    return my_bernoulli_logit_lpmf(y | args1);
}
// Needed to compute the joint likelihood
real my_bernoulli_logit_lpmf(int y, real theta) {
    return bernoulli_logit_lpmf(y | theta);
}
array[] int my_bernoulli_int_logit_rng(
    int anontok__1,
    vector theta
) {
    int n = anontok__1;
    if (dims(theta)[1] != n) reject("my_bernoulli_rng: dim mismatch — `theta` dim 1 (= ", dims(theta)[1], ") does not match `n` (= ", n, "), inferred from `anontok__1` dim 1. `n` sizes: `anontok__1` dim 1 (= ", anontok__1, "), `theta` dim 1 (= ", dims(theta)[1], ").");
    return bernoulli_logit_rng(theta);
}
}
data {
    int t_n;
    vector[t_n] t;
    int y_n;
    array[y_n] int y;
}
transformed data {
    real theta_c1_alpha_upper = upper_alpha_logit();
    real theta_c2_alpha_upper = upper_alpha_logit();
    int n = y_n;
}
parameters {
    real<upper=theta_c1_alpha_upper> theta_c1_alpha;
    real<upper=0.0> theta_c1_beta;
    real<upper=theta_c2_alpha_upper> theta_c2_alpha;
    real<upper=0.0> theta_c2_beta;
    real<lower=0, upper=1> theta_lambda;
}
transformed parameters {
    vector[t_n] theta_c1 = (theta_c1_alpha + (theta_c1_beta * t));
    vector[t_n] theta_c2 = (theta_c2_alpha + (theta_c2_beta * t));
    vector[t_n] theta = ((theta_lambda * logit(theta_c1)) + ((1 - theta_lambda) * logit(theta_c2)));
}
model {
    theta_c1_alpha ~ normal(0, 0.5);
    theta_c1_beta ~ normal(0, 0.5);
    theta_c2_alpha ~ normal(0, 0.5);
    theta_c2_beta ~ normal(0, 0.5);
    theta_lambda ~ beta(2, 2);
    y ~ bernoulli(theta);
    // The exact implementation of the likelihood depends on the passed link function `link_f`
    y ~ my_bernoulli_logit(theta);
}
generated quantities {
    vector[y_n] y_likelihood = bernoulli_lpmfs(y, theta);
    array[y_n] int y_gen = bernoulli_int_rng(y_n, theta);
    // The exact implementation of the likelihood depends on the passed link function `link_f`
    vector[y_n] y_likelihood = my_bernoulli_logit_lpmfs(y, theta);
    // The exact implementation of the likelihood depends on the passed link function `link_f`
    array[t_n] int y_gen = my_bernoulli_int_logit_rng(y_n, theta);
}

regression_mix / log

functions {
// The upper bound of the alpha parameter (for the log link function)
int upper_alpha_log() {
    return 0;
}
vector bernoulli_lpmfs(
    array[] int obs,
    vector args1
) {
    return jbroadcasted_bernoulli_lpmfs(obs, args1);
}
vector jbroadcasted_bernoulli_lpmfs(
    array[] int x1,
    vector x2
) {
    int n = dims(x1)[1];
    vector[n] rv;
    for(i in 1:n) {
        rv[i] = bernoulli_lpmfs(broadcasted_getindex(x1, i), broadcasted_getindex(x2, i));
    }
    return rv;
}
real bernoulli_lpmfs(int args1, real args2) {
    return bernoulli_lpmf(args1 | args2);
}
int broadcasted_getindex(array[] int x, int i) {
    return x[i];
}
real broadcasted_getindex(vector x, int i) {
    return x[i];
}
array[] int bernoulli_int_rng(
    int anontok__1,
    vector p
) {
    int n = anontok__1;
    return bernoulli_rng(p);
}
// Needed to compute the joint likelihood
real my_bernoulli_log_lpmf(
    array[] int y,
    vector theta
) {
    return bernoulli_lpmf(y | exp(theta));
}
// Needed for cross validation
vector my_bernoulli_log_lpmfs(
    array[] int y,
    vector args1
) {
    return jbroadcasted_my_bernoulli_lpmfs_log(y, args1);
}
vector jbroadcasted_my_bernoulli_lpmfs_log(
    array[] int x1,
    vector x3
) {
    int n = dims(x1)[1];
    vector[n] rv;
    for(i in 1:n) {
        rv[i] = my_bernoulli_log_lpmfs(broadcasted_getindex(x1, i), broadcasted_getindex(x3, i));
    }
    return rv;
}
// Needed for cross validation
real my_bernoulli_log_lpmfs(int y, real args1) {
    return my_bernoulli_log_lpmf(y | args1);
}
// Needed to compute the joint likelihood
real my_bernoulli_log_lpmf(int y, real theta) {
    return bernoulli_lpmf(y | exp(theta));
}
array[] int my_bernoulli_int_log_rng(
    int anontok__1,
    vector theta
) {
    int n = anontok__1;
    if (dims(theta)[1] != n) reject("my_bernoulli_rng: dim mismatch — `theta` dim 1 (= ", dims(theta)[1], ") does not match `n` (= ", n, "), inferred from `anontok__1` dim 1. `n` sizes: `anontok__1` dim 1 (= ", anontok__1, "), `theta` dim 1 (= ", dims(theta)[1], ").");
    return bernoulli_rng(exp(theta));
}
}
data {
    int t_n;
    vector[t_n] t;
    int y_n;
    array[y_n] int y;
}
transformed data {
    int theta_c1_alpha_upper = upper_alpha_log();
    int theta_c2_alpha_upper = upper_alpha_log();
    int n = y_n;
}
parameters {
    real<upper=theta_c1_alpha_upper> theta_c1_alpha;
    real<upper=0.0> theta_c1_beta;
    real<upper=theta_c2_alpha_upper> theta_c2_alpha;
    real<upper=0.0> theta_c2_beta;
    real<lower=0, upper=1> theta_lambda;
}
transformed parameters {
    vector[t_n] theta_c1 = (theta_c1_alpha + (theta_c1_beta * t));
    vector[t_n] theta_c2 = (theta_c2_alpha + (theta_c2_beta * t));
    vector[t_n] theta = ((theta_lambda * log(theta_c1)) + ((1 - theta_lambda) * log(theta_c2)));
}
model {
    theta_c1_alpha ~ normal(0, 0.5);
    theta_c1_beta ~ normal(0, 0.5);
    theta_c2_alpha ~ normal(0, 0.5);
    theta_c2_beta ~ normal(0, 0.5);
    theta_lambda ~ beta(2, 2);
    y ~ bernoulli(theta);
    // The exact implementation of the likelihood depends on the passed link function `link_f`
    y ~ my_bernoulli_log(theta);
}
generated quantities {
    vector[y_n] y_likelihood = bernoulli_lpmfs(y, theta);
    array[y_n] int y_gen = bernoulli_int_rng(y_n, theta);
    // The exact implementation of the likelihood depends on the passed link function `link_f`
    vector[y_n] y_likelihood = my_bernoulli_log_lpmfs(y, theta);
    // The exact implementation of the likelihood depends on the passed link function `link_f`
    array[t_n] int y_gen = my_bernoulli_int_log_rng(y_n, theta);
}

Where the parameter type and size come from

Sampling syntax needs more than a log-density: StanBlocks must also know the left-hand side's Stan type, constraints, and dimensions. A distribution method can provide that contract symbolically. In this example,

@lhs @lpxf hetero_lpdf(x::ordered[n], ::typeof(logit), n) = hetero_lpdf(x)

says that xi ~ hetero(logit, n) introduces an ordered[n] parameter, while the corresponding typeof(log) method introduces a positive_ordered[n] parameter. The ordinary hetero_lpdf(x::vector[n]) method holds the density calculation shared by both constrained interfaces.

That mechanism is still useful for reusable custom distributions: the caller does not repeat a constraint that is intrinsic to the distribution. It is no longer the only way to communicate intent, however. Current StanBlocks also supports explicit typed left-hand sides such as

theta::vector[n] ~ distribution(foo, bar)

when the declaration belongs at the model site. The preserved PCR source uses @lhs because the same distribution name intentionally selects a different transform for each function token. This is a design choice in the example, not an outstanding type-annotation stub.