This is an old version, view current version.

## 20.2 Functions as statements

In some cases, it makes sense to have functions that do not return a value. For example, a routine to print the lower-triangular portion of a matrix can be defined as follows.

functions {
void pretty_print_tri_lower(matrix x) {
if (rows(x) == 0) {
print("empty matrix");
return;
}
print("rows=", rows(x), " cols=", cols(x));
for (m in 1:rows(x)) {
for (n in 1:m) {
print("[", m, ",", n, "]=", x[m, n]);
}
}
}
}

The special symbol void is used as the return type. This is not a type itself in that there are no values of type void; it merely indicates the lack of a value. As such, return statements for void functions are not allowed to have arguments, as in the return statement in the body of the previous example.

Void functions applied to appropriately typed arguments may be used on their own as statements. For example, the pretty-print function defined above may be applied to a covariance matrix being defined in the transformed parameters block.

transformed parameters {
cov_matrix[K] Sigma;
// ... code to set Sigma ...
pretty_print_tri_lower(Sigma);
// ...
}
stan

## Functions accessing the log probability accumulator

Functions whose names end in _lp are allowed to use sampling
statements and target += statements; other
functions are not.  Because of this access, their use is restricted to
the transformed parameters and model blocks.

Here is an example of a function to assign standard normal priors to a
vector of coefficients, along with a center and scale, and return the
translated and scaled coefficients; see the [reparameterization
efficient non-centered parameterizations

stan
functions {
vector center_lp(vector beta_raw, real mu, real sigma) {
beta_raw ~ std_normal();
sigma ~ cauchy(0, 5);
mu ~ cauchy(0, 2.5);
return sigma * beta_raw + mu;
}
// ...
}
parameters {
vector[K] beta_raw;
real mu_beta;
real<lower=0> sigma_beta;
// ...
}
transformed parameters {
vector[K] beta;
// ...
beta = center_lp(beta_raw, mu_beta, sigma_beta);
// ...
}
stan

## Functions acting as random number generators

A user-specified function can be declared to act as a (pseudo) random
number generator (PRNG) by giving it a name that ends in _rng.
Giving a function a name that ends in _rng allows it to access
built-in functions and user-defined functions that end in
_rng, which includes all the built-in PRNG functions.  Only
functions ending in _rng are able access the built-in PRNG
functions.  The use of functions ending in _rng must therefore
be restricted to transformed data and generated quantities blocks like
other PRNG functions; they may also be used in the bodies of other
user-defined functions ending in _rng.

For example, the following function generates an $N \times K$ data
matrix, the first column of which is filled with 1 values for the
intercept and the remaining entries of which have values drawn
from a standard normal PRNG.

stan
matrix predictors_rng(int N, int K) {
matrix[N, K] x;
for (n in 1:N) {
x[n, 1] = 1.0;  // intercept
for (k in 2:K) {
x[n, k] = normal_rng(0, 1);
}
}
return x;
}

The following function defines a simulator for regression outcomes based on a data matrix x, coefficients beta, and noise scale sigma.

vector regression_rng(vector beta, matrix x, real sigma) {
vector[rows(x)] y;
vector[rows(x)] mu;
mu = x * beta;
for (n in 1:rows(x)) {
y[n] = normal_rng(mu[n], sigma);
}
return y;
}

These might be used in a generated quantity block to simulate some fake data from a fitted regression model as follows.

parameters {
vector[K] beta;
real<lower=0> sigma;
// ...
}
generated quantities {
matrix[N_sim, K] x_sim;
vector[N_sim] y_sim;
x_sim = predictors_rng(N_sim, K);
y_sim = regression_rng(beta, x_sim, sigma);
}

A more sophisticated simulation might fit a multivariate normal to the predictors x and use the resulting parameters to generate multivariate normal draws for x_sim.