diff --git a/src/_quarto.yml b/src/_quarto.yml index e54edc4a0..e418df73c 100644 --- a/src/_quarto.yml +++ b/src/_quarto.yml @@ -150,6 +150,7 @@ website: - stan-users-guide/reparameterization.qmd - stan-users-guide/efficiency-tuning.qmd - stan-users-guide/parallelization.qmd + - stan-users-guide/user_defined_constraints.qmd - section: "Posterior Inference & Model Checking" contents: - stan-users-guide/posterior-prediction.qmd diff --git a/src/bibtex/all.bib b/src/bibtex/all.bib index fbbe9e004..c5e0f87c4 100644 --- a/src/bibtex/all.bib +++ b/src/bibtex/all.bib @@ -554,6 +554,17 @@ @MISC{Boost:2011 {http://en.highscore.de/cpp/boost/} } +@article{Bou-Rabee:2025, + title={Incorporating local step-size adaptivity into the no-{U}-turn sampler using {G}ibbs self-tuning}, + author={Bou-Rabee, Nawaf and Carpenter, Bob and Kleppe, Tore Selland and Marsden, Milo}, + journal={The Journal of Chemical Physics}, + volume={163}, + number={8}, + year={2025}, + publisher={AIP Publishing} +} + + @article{BowlingEtAl:2009, author={Bowling, Shannon R. and Khasawneh, Mohammad T. and Kaewkuekool, Sittichai and Cho, Byung Rae}, year={2009}, title={A logistic approximation to @@ -1547,3 +1558,11 @@ @book{nicenboim2025introduction year={2025}, publisher={CRC Press} } + +@manual{spiegelhalter2003bugs, + title = {{BUGS} Examples, Volume 3 (Version 0.5)}, + author = {Spiegelhalter, D. J. and Thomas, A. and Best, N. G. and Gilks, W. R. and Lunn, D. J.}, + year = {2003}, + organization = {MRC Biostatistics Unit, Institute of Public Health}, + address = {Cambridge, UK} +} \ No newline at end of file diff --git a/src/stan-users-guide/_quarto.yml b/src/stan-users-guide/_quarto.yml index 4ee186407..4d3bc188b 100644 --- a/src/stan-users-guide/_quarto.yml +++ b/src/stan-users-guide/_quarto.yml @@ -71,6 +71,7 @@ book: - reparameterization.qmd - efficiency-tuning.qmd - parallelization.qmd + - user_defined_constraints.qmd - part: "Posterior Inference & Model Checking" chapters: - posterior-prediction.qmd diff --git a/src/stan-users-guide/img/disc.png b/src/stan-users-guide/img/disc.png new file mode 100644 index 000000000..ab02d3190 Binary files /dev/null and b/src/stan-users-guide/img/disc.png differ diff --git a/src/stan-users-guide/img/hexagon.png b/src/stan-users-guide/img/hexagon.png new file mode 100644 index 000000000..21e267788 Binary files /dev/null and b/src/stan-users-guide/img/hexagon.png differ diff --git a/src/stan-users-guide/img/hollow-square.png b/src/stan-users-guide/img/hollow-square.png new file mode 100644 index 000000000..d9b030eca Binary files /dev/null and b/src/stan-users-guide/img/hollow-square.png differ diff --git a/src/stan-users-guide/img/parallelogram.png b/src/stan-users-guide/img/parallelogram.png new file mode 100644 index 000000000..5f1fe812b Binary files /dev/null and b/src/stan-users-guide/img/parallelogram.png differ diff --git a/src/stan-users-guide/img/ring.png b/src/stan-users-guide/img/ring.png new file mode 100644 index 000000000..5eb9da150 Binary files /dev/null and b/src/stan-users-guide/img/ring.png differ diff --git a/src/stan-users-guide/img/soft-ring.png b/src/stan-users-guide/img/soft-ring.png new file mode 100644 index 000000000..55e8c60d9 Binary files /dev/null and b/src/stan-users-guide/img/soft-ring.png differ diff --git a/src/stan-users-guide/user_defined_constraints.qmd b/src/stan-users-guide/user_defined_constraints.qmd new file mode 100644 index 000000000..27f8380bb --- /dev/null +++ b/src/stan-users-guide/user_defined_constraints.qmd @@ -0,0 +1,1351 @@ +--- +pagetitle: User-Defined Constrained Parameterizations +--- + +# User-Defined Constrained Parameters + +This chapter will cover how to define parameterizations that satisfy +constraints that are not handled by Stan's built-in types. There +are several general strategies used to do this. The easiest is to build +the type out of Stan's existing constrained types and transforms, as +that will let Stan automatically work out the change-of-variables. +The more general approach is to code the transform function and a +log Jacobian increment to ensure the baseline constrained distribution +is uniform over parameter values that satisfy the constraints. + +The change-of-variables adjustment issue is a vexing one. It is +only strictly necessary when a distribution is going to be placed +on the transformed parameters. If a variable is only being transformed +to satisfy some constraint and is then used on the right-hand side of +a distribution statement (or the direct equivalent), then it does not +need a change-of-variables adjustment. + +## Specialized matrices + +There are a large number of specialized matrix forms, several of which +will be described in this section. + +### Diagonal matrices + +Whenever possible, diagonal matrices should be avoided. Stan builds in +specialized functions to multiply by diagonal matrices that only +requires the diagonal values; see [diagonal matrix +functions](https://mc-stan.org/docs/functions-reference/matrix_operations.html#diagonal-matrix-functions) +in the functions reference manual. + +Stan doesn't provide a diagonal matrix type, but it builds in +the necessary transform as `diag_matrix()`. Given an $N$-vector +`v`, `diag_matrix(v)` is the $N \times N$ matrix with diagonal `v` +and zeros elsewhere. No Jacobian is required as it is just filling +in zero values. + + +### Tridiagonal matrices + +A tridiagonal matrix has zero entries everywhere other than on the +main diagonal and the diagonals below and above it. An $N \times N$ +tridiagonal matrix has $3 \cdot N - 2$ free entries and the rest are +structural zeros. The matrix is filled from the vector in row +major order. + +```stan +matrix tridiagonal_matrix(vector x, int N) { + if (rows(x) != 3 * N - 2) { + reject("rows(x) = ", rows(x), + "; but 3 * N - 2 = ", 3 * N - 2); + } + if (N == 1) return [[x[1]]]; + matrix[N, N] y; + y[1, 1] = x[1]; + y[1, 2] = x[2]; + int next = 3; + for (n in 2:(N - 1)) { + y[n, n - 1] = x[next]; + y[n, n] = x[next + 1]; + y[n, n + 1] = x[next + 2]; + next += 3; + } + y[N, N - 1] = x[next]; + y[N, N] = x[next + 1]; + return y; +} +``` + +### Toeplitz matrices + +A Toeplitz matrix has constant values along the diagonal. Thus +an $M \times N$ Toeplitz matrix has $M + N - 1$ free variables. These +are filled in from the lower left corner $(M, 1)$ to the upper right $(1, N)$. + +```stan +matrix toeplitz_matrix(vector x, int M, int N) { + if (num_elements(x) != M + N - 1) { + reject("x must have M + N - 1 elements", + "found num_elements(x) = ", num_elements(x), + " but M + N - 1 = ", M + N - 1); + } + matrix[M, N] y; + for (i in 1:M) { + for (j in 1:N) { + y[i, j] = x[M + j - i]; + } + } + return y; +} +``` + +To verify that the lower left corner is assigned `x[1]`, note that the +lower left corner has `i` equal to M and `j` equal to 1, so `M + j - +i` evaluates to 1, as intended. To verify the upper right is assigned +`x[M + N - 1]`, note that the upper upper right has `i` equal to 1 and +`j` equal to `N`, so `M + j - i` is `M + N - 1`. To verify the +diagonals are aligned, note that the requirement is that `M[i, j]` is +the same as `M[i - 1, j + 1]`, which it is because `M + j - i` is +equal to `M + (j - 1) - (i - 1)`. + + +### Lower-triangular matrices + +An $N \times N$ matrix $A$ is lower triangular if $A_{i, j} = 0$ for +all $j > i$. There are $N + \binom{N}{2}$ variables required to +construct a lower-triangular matrix. Because Stan cannot convert a +real value to an integer, the dimension of the resulting matrix is +required as an argument. + +```stan +matrix lower_triangular_matrix(vector x, int N) { + if (cols(x) != N + binom(N, 2)) { + reject("cols(x) = ", cols(x), + "; but N + binom(N, 2) = ", N + binom(N, 2)); + } + matrix[N, N] y = rep_matrix(0, N, N); // fills zero entries + int next = 0; + for (i in 1:N) { + for (j in i) { + y[i, j] = x[next]; + next += 1; + } + } + return y; +} +``` + +Upper-triangular matrices can be constructed in the same way. + +### Symmetric matrices + +An $N \times N$ matrix $A$ is symmetric if $A = A^\top$ (i.e., +$A_{i,j} = A_{j, i}$ for all $1 \leq i, j \leq N$). As with a lower +triangular matrix, there are $N + \binom{N}{2}$ free variables in a +symmetric matrix. To build a symmetric matrix, it suffices to +build a lower-triangular matrix, then copy the entries below the +diagonal into their transposed position above the diagonal. + +```stan +matrix symmetric_matrix(vector x, int N) { + matrix[N, N] y = lower_triangular_matrix(x, N); + for (i in 1:N) { + for (j in 1:(i - 1)) { + y[j, i] = y[i, j]; + } + } + return y; +} +``` + +### Skew-symmetric matrices + +An $N \times N$ matrix $A$ is skew symmetric if its negation is equal to its transpose, +$-A = A^\top$, or equivalently, if it's equal to its negated transpose, $A = -A^\top$. +Building a skew symmetric matrix is just like building a symmetric matrix only the +above-diagonal entries are negated. + +```stan +matrix skew_symmetric_matrix(vector x, int N) { + matrix[N, N] y = lower_triangular_matrix(x, N); + for (i in 1:N) { + for (j in 1:(i - 1)) { + y[j, i] = -y[i, j]; + } + } + return y; +} +``` + +The only difference form the symmetric matrix case is the negation. + + +## Ordered simplexes + +There are often requests for defining simplex distributions that are +ordered, so that for $\theta \in \Delta^{N-1}$, not only is +$\theta_n > 0$ and $\sum{n=1}^N \theta_n = 1$, but also $\theta_n +< \theta_{n+1}$ for $n < N$. + +Now that Stan has the underlying transforms built in, an ordered +simplex can be defined by taking an ordered vector and applying the +simplex transform. + +```stan +data { + int N; +} +parameters { + positive_ordered[N - 1] log_theta; +} +transformed parameters { + simplex[N] theta = simplex_jacobian(-log_theta); +} +``` + +The negation on `log_theta` makes sure the inputs to +`simplex_jacobian` are negative. The internal transform in +`simplex_jacobian` converts this to a sum-to-zero vector in increasing +order, then applies `softmax()`, then adds the required log +change-of-variables adjustment to the log Jacobian accumulator. + +Enforcing an increasing sequence pattern on the simplex leads to a +uniform distribution with small means on the early components and +relatively large means on the later ones. For example, with `N` set +to 2, the posterior means are $[0.25 \ 0.75]$, with `N` set to 5, they +are are $[0.05 \ 0.08 \ 0.12 \ 0.23 \ 0.52]$, and set to 10, they are +$[0.02 \ 0.02 \ 0.03 \ 0.04 \ 0.05 \ 0.07 \ 0.10 \ 0.14 \ 0.20 \ +0.33]$. To derive a more uniform distribution, they can be given +a Dirichlet prior that's skewed in the opposite direction. For example, +consider adding the model block for the case where `N` is 5. +```stan +model { + theta ~ dirichlet([16, 8, 4, 2, 1]'); +} +``` +The resulting posterior means are now $[0.17 \ 0.18 \ 0.19 \ 0.21 \ +0.25]$. + +It feels counterintuitive that giving such a skewed Dirichlet prior +leads to an almost uniform ordered distribution, but that's just the +natural result of two densities pulling in opposite direction. The +uniform distribution on ordered simplexes is heavily skewed. This is +handled by the change-of-variables adjustments Stan does under the +hood for the composition of simplex constraints with positive-ordered +constraints. The explicit Dirichlet prior in the `model` block then +pulls it back toward uniformity. This overall point is characteristic +of constrained distribution, which will not behave like uniform +distributions on the unconstrained space when the transforms are +non-linear. + + +## Correlation matrices with structural zeros + +Stan does _not_ provide a built-in type for correlation matrices that +have structural zeros (i.e., correlations known to be zero). Such +constraints arise from graphical models in structural equation +modeling and as priors in varying (aka "random") effects models. This +section shows how to parameterize Cholesky factors of correlation +matrices with specified structural zeros and define matching +restricted LKJ distributions. + +### Stan's correlation matrix types and the LKJ distribution + +Stan provides the built-in constrained types `corr_matrix` and +`cholesky_factor_corr` for correlation matrices and their Cholesky +factors, as in the following Stan program. + +```stan +parameters { + corr_matrix[5] Omega1; + cholesky_factor_corr[5] L2; +} +transformed parameters { + cholesky_factor_corr[5] L1 = cholesky_factor(Omega1); + corr_marrix[5] Omega2 = L * L'; +} +``` + +In this program, `Omega1` will have a uniform distribution over +correlation matrices and `L2` will have a uniform distribution over +lower triangular matrices with positive diagonals and unit-length +rows, i.e., over the Cholesky factors of correlation matrices. This +matches Stan's overall design for constrained types, which by default, +will be uniform over variables satisfying the constraints. Unbounded +constrained variables like covariance matrices will have improper "uniform" +distributions, whereas bounded variables like correlation matrices will +induce proper uniform distributions. + +Given this behavior, it should not be surprising that the Cholesky +factor `L1` of the correlation matrix `Omega1 = L1 * L1'` will _not_ +have a uniform distribution over Cholesky factors of correlation +matrices. It's missing the change-of-variables adjustment resulting +from `L1 * L1'`. For the same reason, `Omega2`, the reconstructed +correlation matrix from the Cholesky factor `L2`, will not have a +uniform distribution over correlation matrices. + +Suppose a model block is now added to the program above giving +`Omega1` an LKJ distribution with concentration 1 and `L2` a matching +LKJ distribution over Cholesky factors. + +```stan +model { + Omega1 ~ lkj_corr(1); + L2 ~ lkj_corr_cholesky(1); +} +``` + +Setting $\eta = 1$ in $\textrm{LKJ}(\eta)$ produces a uniform +distribution over correlation matrices. Thus `Omega1` still has a +uniform distribution over correlation matrices. `L2`, on the other +hand, now has a distribution that matches that of `L1`. In words, `L2` +now has a distribution over Cholesky factors such that when they are +multiplied out have a uniform distribution over correlation matrices. +The `lkj_corr_cholesky` distribution builds in the +change-of-variables adjustment resulting from transforming `L1` to +`Omega1`, then applies the LKJ distribution to `Omega1` (conceptually; +computationally it's more efficient than that). + +With Stan's built-in constraint transformation functions, another way +of specifying a parameter to be the Cholesky factor of a correlation +matrix is as follows. + +```stan +parameters { + vector[choose(5, 2)] L3_raw; +} +transformed parameters { + cholesky_factor_corr[5] L3 + = cholesky_factor_corr_jacobian(L3_raw, 5); +} +``` + +In this program, `L3` will have a uniform distribution over Cholesky +factors of correlation matrices, matching `L1` in the original +program. In particular, it can be given an LKJ-correlation distribution +to make the resulting correlation matrices LKJ distributed. + + +### Specifying structural zeros + +Stan does not support user-defined types. The closest approximation is +user-defined `_jacobian` transform functions that increment the +Jacobian accumulator with a change-of-variables adjustment. The rest +of this section shows how to use this mechanism to develop correlation +matrix Cholesky factors that lead to structural zeros in specified +positions of the resulting correlation matrices. + +The structural zeros will be specified with a simple symmetric +indicator matrix. For example, specifying that a $4 x 4$ matrix does +not have correlations between variables 1 and 2 or between variables 3 +and 4 would be specified as a $5 \times 5$ integer array, + +$$ +\textrm{zero} += +\begin{bmatrix} +0 & 1 & 0 & 0 & 0 \\ +1 & 0 & 0 & 0 & 0 \\ +0 & 0 & 0 & 1 & 0 \\ +0 & 0 & 1 & 0 & 0 \\ +0 & 0 & 0 & 0 & 0 \\ +\end{bmatrix} +$$ + + +Technically, this array should be symmetric with a zero diagonal. In +practice, only the below-diagonal elements will be used, so there +won't be any error checking. This may seem wasteful if there are only +a few structural zeros, but integer arrays are small and the program +will already be generating dense Cholesky factors of the same order +and automatically differentiating through them. In Stan programs, +indicator integer arrays like this can be used as booleans, e.g., `if +(zero[1, 2])` where the condition is true if `zero[1, 2]` evaluates to +a non-zero value and false otherwise. + + +### Constructing Cholesky factors with structural zeros + +There are easier ways to define this transform which are +mathematically correct, but computationally sub-optimal. Rather than +starting with one of those and then moving on, this section introduces +a single, numerically stable version of the transform. The transform +is complicated in realization, but is based on some simple principles. + +The construction of a lower-triangular Cholesky factor of a +correlation matrix must ensure the result is symmetric and positive +definite with a unit diagonal. Positive-definiteness is +straightforward, because for any $D \times D$ lower-triangular matrix +$L$ with a positive diagonal, $\Omega = L \cdot L^\top$ is symmetric +and positive definite. + +Each of the rows of $L$ must be of unit Euclidean length, so that +$$ +L_n \cdot L_n^\top = 1, +$$ +where $L_n$ is read Stan-style as the $n$-th row of $L$ (and hence a +row vector). This is going to be handled by sampling from a unit +hypersphere in $D$ dimensions, because the unit hypersphere is exactly +the set of vectors such that $v^\top \cdot v = 1$. + +Dealing with the structural zeros is trickier. For example, if +$\Omega_{i, j}$ is a structural zero, then $L_{i,j}$ will be defined +as the solution to an equation making $L_i \cdot L_j^\top = 0$, i.e., +keeping the two rows orthogonal. That's it, in a nutshell. + +The Stan function to carry out this transform will be presented +in a sequence of function definitions. These, of course, should +be defined in a `functions` block. The first function is +just a numerically stable log hyperbolic cosine. + +```stan +/* + * Return the numerically stable log of the hyperbolic cosine of + * the argument. + * + * @param x The argument. + * @return The log of the hyperbolic cosine of the argument. + */ +real log_cosh(real x) { + return log_sum_exp(x, -x) - log2(); +} +``` + +The documentation follows the style of the multi-language +documentation generator +[Doxygen](https://en.wikipedia.org/wiki/Doxygen). + +The next helper function extracts the column indexes that match a +given value for a given row of the zero correlation indicator array. + +```stan +/** + * Return the columns `c` in `1:(i - 1)` for which `row_ind[c] == + * val`, in ascending order. Used to split the columns of row `i` + * of a Cholesky factor into structurally zero columns (`val == 1`) + * and free columns (`val == 0`). + * + * @param i Row index. + * @param row_ind Indicator array of size >= i - 1 (row i of the + * zero-indicator matrix). + * @param val Indicator value to match (0 or 1). + * @return Ascending array of matching column indices. + */ +array[] int columns_matching(int i, array[] int row_ind, int val) { + int n = 0; + for (c in 1:(i - 1)) { + n += (row_ind[c] == val); + } + array[n] int idx; + int pos = 0; + for (c in 1:(i - 1)) { + if (row_ind[c] == val) { + pos += 1; + idx[pos] = c; + } + } + return idx; +} +``` + +The next function does all the heavy lifting. It has a long documentation +string explaining what it does. + +```stan + /** + * Transform unconstrained parameters to the Cholesky factor of a + * correlation matrix Omega = L * L' with structural zeros, i.e., + * Omega[i, j] == 0 wherever is_zero[i, j] == 1 (below diagonal; + * diagonal and above are ignored). A change-of-variables + * adjustment is applied to the log Jacobian accumulator so that + * if `raw` has an (improper) uniform distribution, the returned + * Cholesky factor has a uniform distribution over Cholesky + * factors that produce the specified structural zeros. + * + * @param D The number of dimensions. + * @param raw The unconstrained parameters. + * @param is_zero A D x D indicator array where `is_zero[i, j]` + * is an indicator that there is zero correlation between + * variable `i` and `j`. + * @return The Cholesky factor of a correlation matrix satisfying + * the structural zeros. + */ +``` + +The code itself follows. + +```stan +matrix cholesky_corr_zeros_jacobian(int D, vector raw, + array[,] int is_zero) { + matrix[D, D] L = rep_matrix(0, D, D); + int raw_idx = 1; + L[1, 1] = 1; + for (i in 2:D) { + int size_fc = i - 1 - sum(is_zero[i, 1:(i - 1)]); + // free columns: + array[size_fc] int fc = columns_matching(i, is_zero[i], 0); + int size_zc = sum(is_zero[i, 1:(i - 1)]); + // zero columns: + array[size_zc] int zc = columns_matching(i, is_zero[i], 1); + int d = size(fc); // # free dims + int n_z = size(zc); // # structural zeros + row_vector[d] z; // point in d-dimensional unit ball + real r = 1; // leftover radius + // invariant: r^2 == 1 - dot_self(z). + for (k in 1:d) { + real x = raw[raw_idx]; + z[k] = r * tanh(x); + jacobian += -(d - k + 2) * log_cosh(x); + r /= cosh(x); + raw_idx += 1; + } + L[i, i] = r; + if (n_z == 0) { + L[i, 1:d] = z; // unconstrained row uses z directly + } else if (d > 0) { + // B[k, m]: entry at zero column zc[m] of the basis direction for + // free column fc[k] patched so that the direction is orthogonal + // to row zc[m] of L. + matrix[d, n_z] B = rep_matrix(0, d, n_z); + for (k in 1:d) { + for (m in 1:n_z) { + int j = zc[m]; + if (j > fc[k]) { + real b = L[j, fc[k]] + + dot_product(L[j, zc[1:(m - 1)]], + B[k, 1:(m - 1)]); + B[k, m] = -b / L[j, j]; // solve to make correlation 0 + } + } + } + // overlap (Gram) matrix of the basis directions is + // G = I + B * B'; straighten with Cholesky factor C so + // row's length matches length of z + matrix[d, d] C + = cholesky_decompose(add_diag(tcrossprod(B), 1.0)); + jacobian += -sum(log(diagonal(C))); // free entry Jacobians + row_vector[d] t + = mdivide_right_tri_low(z, C); // free entries + L[i, fc] = t; + L[i, zc] = t * B; // forced entries + } + } + return L; +} +``` + +To unpack the solve required for zeros, consider the situation where +$\Omega_{i, j}$ is a structural zero. Then $L_{i, j}$ will be defined +so that $L_i \cdot L_j^\top = 0$. The lower-triangular form means that +the product defining entry $(i, j)$ in the correlation matrix only +involves the first $k = \textrm{min}(i, j)$ entries in $L_i$ and +$L_j$, and $j < i$ because of the lower triangular structure, so the +condition that must be satisfied is +$$ +\Omega_{i, j} +\ = \ L_i \cdot L_j^\top +\ = \ L_{i, 1:j} \cdot L_{j, 1:j}^\top +\ = \ 0. +$$ +Let +$$ +b = L_{j, 1:(j-1)} \cdot L_{j, 1:(j-1)}^\top. +$$ +Solving for $L_{i,j}$ in $\Omega_{i,j} = 0$ yields +$$ +L_{i, j} = \dfrac{-b}{L_{j,j}}. +$$ + +### LKJ distribution for Cholesky factors leading to structural zeros + +As described above, this new transform results in a uniform distribution over +Cholesky factors that produce correlation matrices with the specified structural +zeros. In order to use it, a specialized form of the `lkj_corr_cholesky` distribution +must be defined to deal with the reduction in the number of underlying parameters due +to the structural zeros (each zero removes an underlying dimension). + +```stan +/** + * Return the log density of the LKJ distribution with the + * specified concentration for the specified Cholesky factor + * of a correlation matrix, corrected for the specified + * structural zeros. + * + * @param L Cholesky factor of a correlation matrix. + * @param eta LKJ concentration (eta > 0). + * @param col_zeros Per-column counts of below-diagonal + * structural zeros. + * @return Log density of `L`, corrected for structural zeros. + */ +real lkj_corr_cholesky_zeros_lpdf(matrix L, real eta, + row_vector col_zeros) { + return lkj_corr_cholesky_lpdf(L | eta) + - col_zeros * log(diagonal(L)); + } +} +``` + +### Stan program using the transform + +With all the heavy lifting done in functions, the Stan program that depends +on them is simple. + +```stan +data { + int D; // dimensions + real eta; // LKJ concentration + array[D, D] int zeros; // structural zeros +} +transformed data { + row_vector[D] col_zeros; + for (j in 1:D) { + col_zeros[j] = sum(zeros[(j + 1):D, j]); + } + int N_zero = sum(col_zeros); +} +parameters { + vector[choose(D, 2) - N_zero] raw; +} +transformed parameters { + matrix[D, D] L_Omega + = cholesky_corr_zeros_jacobian(D, raw, zeros); +} +model { + L_Omega ~ lkj_corr_cholesky_zeros(eta, col_zeros); +} +generated quantities { + matrix[D, D] Omega + = multiply_lower_tri_self_transpose(L_Omega); +} +``` + +This just declares the data, including the LKJ concentration and the +structural zeros. The transformed data pulls out the column zero +counts which will be needed by he new Cholesky distribution defined +above. The number of raw parameters is the number of correlation +matrix parameters minus the number of structural zeros. The variable +`L_Omega` will be implicitly given a uniform distribution over +Cholesky factors that lead to the specified structural zeros. With +the LKJ distribution in the `model` block, the result is that `Omega = +L_Omega * L_Omega'` will have an $\textrm{LKJCorr}(\eta)$ +distribution. + +For each pair of non-zero indexes below the diagonal, there is a +Jacobian adjustment for the change of variables involved in the +constraining (inverse) transform due to the lower/upper bounding. This +adjustment is handled implicitly by Stan's built-in lower/upper-bound +transform. + +This defines a variable `L_Omega` that is uniform over Cholesky +factors of correlation matrices with the specified structural zeros. +The final result is an `Omega` with an $\textrm{LKJ}(\nu)$ +distribution restricted to correlation matrices with the specified +structural zeros. This works no matter what $\eta > 0$ is (e.g., +setting `eta` to 1 produces a uniform distribution). + +As part of a larger application, the correlation matrix or its +Cholesky factor are typically paired with a vector of scale parameters +to produce a covariance matrix with the specified pattern of zeros. + +As an example, consider a $D = 5$ example where setting $\Omega_{1, 4} = 0$, $\Omega_{2,3} = 0$ +and $\Omega{4,5} = 0$, as representing by the following `zeros` matrix. +$$ +\textrm{zeros} += +\begin{bmatrix} +0 & 0 & 0 & 1 & 0 +\\ +0 & 0 & 1 & 0 & 0 +\\ +0 & 1 & 0 & 0 & 0 +\\ +1 & 0 & 0 & 0 & 1 +\\ +0 & 0 & 0 & 1 & 0 +\end{bmatrix}. +$$ +Taking $\eta = 1$ and running with default NUTS setting, the result is a minimum +effective sample size of about half the number of iterations. Here's a representative +draw, +$$ +\Omega^{(112)} += +\begin{bmatrix} +1 & -0.224 & 0.024 & 0 & -0.372 \\ +-0.224 & 1 & 0 & 0.214 & 0.544 \\ + 0.024 & 0 & 1 & -0.072 & 0.737 \\ + 0 & 0.214 & -0.072 & 1 & 0 \\ +-0.372 & 0.544 & 0.737 & 0 & 1 +\end{bmatrix}. +$$ +Being floating point, the zero values are not exactly zero, but +on the order of $10^{-17}$, which is roughly within machine precision +for double-precision floating point like Stan uses. + + +## Orthogonal, rotation, and reflection matrices + +Stan has a built-in `unit_vector` type that produces vectors where the +sum of the squared entry is 1. Those points live on a unit +hypersphere. This section generalizes this notation to orthonormal +matrices, rotation matrices, and reflection matrices. + +### Orthogonality and orthonormality + +A pair of $N$-vectors $u, v$ is _orthogonal_ if $u^\top \cdot v = 0$, +which entails a cosine of zero. Orthogonal is the $N$-dimensional +generalization of the opposite and adjacent edges of a right triangle. +An $N$-vector $u$ has Euclidean _length_ $||u||_2 = \sqrt{u^\top \cdot +u}$. Thus it has unit length if $||u||_2 = 1$. A pair of vectors is +_orthonormal_ if they are orthogonal and unit length. + +An _orthonormal matrix_ is an $M \times N$ matrix in which the columns +are orthonormal (waring: in some treatments rows are orthonormal). In +a somewhat confusing terminological convention, if $M = N$, the matrix +is further said to be _orthogonal_ (the confusion arises because the +columns are orthogonal in an orthonormal matrix even if $M \neq N$). + +Stan's built-in unit vector type (`unit_vector`) produces the boundary +case of an orthonormal matrix with only a single column. + +### Rotation and reflection + +If $y$ is a point on a hypersphere (e.g., a unit hypersphere) and $X$ +is orthogonal, then $y \cdot X$ is on the same hypersphere. A point +on the hypersphere can be derived by either rotation or reflection. + +An orthogonal matrix $X$ with $\det X = 1$ is a _rotation matrix_ in +that its action on a vector $y$, $y \cdot X$, is a rotation. You +can get from any point on the sphere to any other point on a sphere +through a rotation. + +An orthogonal $X$ with $\det X = -1$ is said to be a _reflection +matrix_. A reflection matrix maps a point through the origin to the +opposite pole. + +### The special orthogonal group + +The set of rotation matrices of dimension $N \times N$ forms a group +structure which is written as $\textrm{SO}(N)$ for _special orthogonal +group_. It's easy to verify that the rotation matrices meet the +requirements to form a group under multiplication. First, +multiplication is just matrix multiplication, so it is associative. +So if $X_1, X_2$ are rotations, then applying $X_1$ $X_1 \cdot X_2$ is a rotation. + +First, +$\textrm{SO}(N)$ is clearly closed under multiplication, because +applying two rotations is another rotation. The identity rotation is +the identity matrix $\textrm{I}$. Finally, every rotation has an +inverse that is also a rotation, namely its transpose. + +### Generating orthonormal matrices via QR decomposition + +Suppose that $X$ is an $M \times N$ matrix with independent standard +normal entries, i.e., $X_{m,n} \sim \textrm{normal}(0, 1)$. Then the +$Q$ matrix resulting from a QR-decomposition of $X = Q \cdot R$ has a +uniform distribution over orthonormal matrices. + +The Stan program follows the math, as usual. + +```stan +data { + int M, N; +} +parameters { + matrix[M, N] X_raw; +} +transformed parameters { + matrix[M, N] X = qr_thin_Q(X_raw); +} +model { + to_vector(X_raw) ~ normal(0, 1); +} +``` + +This can be used in a program requiring uniformly distributed +orthonormal matrices $X$. If $M = N$, the result is an orthogonal +matrix. + +The QR decomposition is not unique, in general, because of sign +flipping. Stan's QR decomposition takes the canonical form where all +diagonal entries are positive. A little linear algebra shows the +resulting density is rotation invariant and thus uniformly distributed +on orthonormal matrices. + +The `X` that results from this construction is uniformly distributed +over orthonormal matrices. Here's an example draw for the orthogonal +case where $M = N = 5$. It's easy to verify that each row has unit +Euclidean length a zero dot product with each of the other rows. +$$ +X^{(100)} = \begin{bmatrix} +0.207 & -0.697 & -0.388 & 0.295 & -0.483 \\ +-0.631 & -0.494 & -0.224 & -0.138 & 0.538 \\ +0.738 & -0.294 & 0.060 & -0.199 & 0.571 \\ +-0.118 & -0.409 & 0.892 & 0.092 & -0.120 \\ +-0.010 & -0.127 & -0.015 & -0.920 & -0.371 +\end{bmatrix} +$$ + + +### Uniform distributions on SO(_N_) + +To restrict sampling of orthogonal matrices to just the rotations, +the determinant of the constructed variable `X` needs to be pinned +to 1. To do this, it suffices to flip the sign of a single row. +To generate only rotation matrices, the previous `transformed parameters` +block can be replaced with the following. + +```stan +transformed parameters { + matrix[N, N] X = qr_thin_Q(X_raw); + if (determinant(X) < 0) { // force X to be a rotation + X[, N] = -X[, N]; + } +} +``` + +Replacing the condition `determinant(X) < 0` with `determinant(X) > 0` +will force `X` to be a reflection rather than a rotation. + +The induced geometry from this transform with the sign flip may +seem like it would be challenging for sampling due to the discrete +parameter-dependent discontinuity introduced. It shouldn't be that +bad because of properties of the manifold of orthogonal matrices---there +is a zero density gap between the rotation and reflection matrices, so +the branch will likely get stuck with a single evaluation path and +thus leave everything continuous in practice. + +A safe alternative is to initialize with a rotation matrix, such as +the identity matrix. For the same reason as above, it's unlikely to +jump to reflection matrices because of the energy gap. + +Without fiddling with initialization, it would probably +work to just reject in cases where the determinant is not positive, +as follows. + +```stan +transformed parameters { + matrix[N, N] X = qr_thin_Q(X_raw); + if (determinant(X) < 0) { // force X to be a rotation + reject("-1 determinant in Q from QR, but need value 1"); + } +} +``` + +There's a 50% chance a random initialization will get the right +determinant, so this should be work with Stan's default +initialization, which retries up to 100 times, assuming there +aren't other stochastic initialization bottlenecks. + +### Rotations via skew-symmetric matrix exponential + +A matrix $S \in \mathbb{R}^{N \times N}$ is said to be _skew +symmetric_ if it is equal to the negation of its transpose, $-S = +S^\top$. The matrix exponential of a skew-symmetric matrix will +be a rotation. + +An $N \times N$ skew-symmetric matrix only has $\binom{N}{2}$ degrees +of freedom (the diagonal plus points below the diagonal). +Given a vector $v \in \mathbb{R}^{\binom{N}{2}}$, here's a function +to produce a skew-symmetric matrix. + +```stan + matrix skew_symmetric(vector v, int N) { + matrix[N, N] S = rep_matrix(0, N, N); + int idx = 1; + for (i in 1:N) { + for (j in (i + 1):N) { + S[i, j] = v[idx]; + S[j, i] = -v[idx]; + idx += 1; + } + } + return S; + } +``` + +The transform is going to require a Jacobian adjustment, which +will be written like a built-in transform. The transform itself +is rather opaque, as is the following helper function. +The helper function computes $\textrm{sinc}(x / 2) / (x / 2)$ +accurately even if $x$ is near zero. It applies to a whole vector +and then returns 2 times the sum of the log absolute values of that +vector. + +```stan + real two_log_abs_half_sinc(vector x) { + int N = rows(x); + vector[N] h = x / 2; + vector[N] half_sinc = sin(h) ./ h; + // fall back to Taylor expansion for x near 0 + for (n in 1:N) { + if (abs(h[n]) < 1e-4) { + real h_sq = square(h[n]); + half_sinc[n] = 1 - h_sq / 6 + square(h_sq) / 120; + } + } + return 2 * sum(log(abs(half_sinc))); + } +``` + +Depending on the application, this could be simplified. Then there's the transform +itself, which takes a vector $v$ of size $\binom{N}{2}$, creates a skew-symmetric +matrix, applies matrix exponential, then adds the Jacobian adjustment to the Jacobian +accumulator. + +```stan + matrix special_orthogonal_jacobian(vector v, int N) { + matrix[N, N] S = skew_symmetric(v, N); + matrix[N, N] Q_rot = matrix_exp(S); // determinant(Q_rot) = 1 + int m = N / 2; // rounds down (skip last odd eigenvalue) + matrix[N, N] M = -S * S; // symmetric psd + vector[N] eig = eigenvalues_sym(M); + vector[m] theta; // rotation angles + for (i in 1:m) { + real val = eig[n % 2 + 2 * i]; // one from each eq. pair + theta[i] = val > 0 ? sqrt(val) : 0; + } + real log_J = 0; + if (N % 2 == 1) { + jacobian += two_log_abs_half_sinc(theta); + } + for (j in 1:(m - 1)) { + jacobian + += two_log_abs_half_sinc(theta[j] - theta[(j + 1):m]); + jacobian + += two_log_abs_half_sinc(theta[j] + theta[(j + 1):m]); + } + return Q_rot; + } +``` + +Given these function definitions are placed in a `functions` block, the program +to generate random rotation matrices is simple. + +```stan +data { + int N; +} +parameters { + vector[choose(N, 2)] v; +} +transformed parameters { + matrix[N, N] X = special_orthogonal_jacobian(v, N); +} +``` + +The result is a uniform distribution over $X$. + + +## Matrices matching known margins + +Suppose $X$ is an $M \times N$ matrix $X$, the marginals of which are +known, but the entries of which are unknown. A simple outer product +of the marginals scaled by the sum of either of the marginals (they're +the same) provides a baseline matrix with the correct marginals. +Stan's sum-to-zero matrix type can be used to add to that solution +to generate a uniform distribution over matrices that have the specified +marginals. This distribution is improper and thus requires priors +on the matrix entries to determine their scale of variation. + +Let $r \in \mathbb{R}^M$ be a vector of row marginals and $c \in +\mathbb{R}^N$ be a vector of column marginals. There is only a +solution if $\textrm{sum}(r) = \textrm{sum}(c)$. A matrix with +row marginals $r$ and column marginals $c$ is +$$ +E = r \cdot c^\top / \textrm{sum}(c). +$$ +In this matrix, +$$ +E_{m,n} = r_m \cdot c_n / \textrm{sum}(c). +$$ +To see that this is indeed a solution, consider the row $E_m$, +the sum of which is equal to the correct marginal +$r_m$. +\begin{align*} +\textrm{sum}(E_m) +&= E_{m,1} + \cdots + E_{m, N} +\\[4pt] +&= r_m \cdot c_1 / \textrm{sum}(c) + + \cdots + + r_m \cdot c_N / \textrm{sum}(c) +\\ +&= r_m \cdot (c_1 + \cdots + c_N) / \textrm{sum}(c) +\\ +&= r_m. +\end{align*} +To see that the columns $E_{, n}$ have the correct marginals $c_n$, +expand the sum. +\begin{align*} +\textrm{sum}(E_{, n}) +&= E_{1, n} + \cdots + E_{M, n} +\\[4pt] +&= r_1 \cdot c_n / \textrm{sum}(c) ++ \cdots ++ r_M \cdot c_n / \textrm{sum}(c) +\\[4pt] +&= c_n \cdot (r_1 + \cdots + r_M) / \textrm{sum}(c) +\\[4pt] +&= c_n \cdot (r_1 + \cdots + r_M) / \textrm{sum}(r) +\\[4pt] +&= c_n. +\end{align*} +The penultimate step of the derivation exploits the fact that +$\textrm{sum}(c) = \textrm{sum}(r)$ because they are marginals of the +same matrix and thus both equal to the sum of all matrix entries. + +To generate random matrices that satisfy the marginals, the following +Stan programs takes the matrix $E$ defined above and then perturbs it +with a sum-to-zero matrix. Stan's sum-to-zero matrix type is +constrained to have rows that sum to zero and columns that sum to +zero. Thus when adding it to $E$, the result is still a solution. In +the unconstrained domain, there are only $(M - 1) \cdot (N - 1)$ +degrees of freedom. The program takes the more colloquially named +data variables `row_sums` and `col_sums` for the row and column +marginals. The data variable `sigma` determines the prior scale of a +normal distribution on the perturbation values $Z$, which must have zero means due to symmetry considerations. + +An equivalent, but slightly less efficient way to write the prior +would be +```stan +to_vector(X) ~ normal(to_vector(E), sigma); +``` + +equivalently been put on $X$ with slightly higher automatic +differentiation overhead). + +```stan +data { + int M; + int N; + vector[M] row_sums; + row_vector[N] col_sums; + real sigma; +} +transformed data { + matrix[M, N] E = row_sums * (col_sums / sum(col_sums)); +} +parameters { + sum_to_zero_matrix[M, N] Z; +} +transformed parameters { + matrix[M, N] X = E + Z; +} +model { + to_vector(Z) ~ normal(0, sigma); +} +``` + +The definition of `E` is written as `row_sums * (col_sums / +sum(col_sums))`, because it only has to divide `N` elements. Had `E` +been defined as `(row_sums * col_sums) / sum(col_sums)`, the result +would involve `N * M` divisions. If the program will be run with `N` +larger than `M`, then it's more efficient to code as the equivalent +`row_sums / sum(row_sums) * col_sums()`. + +To work through an example, suppose $r = [55 \ 23 \ 22]$ and +$c = [19 \ 18 \ 25 \ 38]$. Then +$$ +r \cdot c / \textrm{sum}(c) += +\begin{bmatrix} +10.45 & 9.90 & 13.75 & 20.90 \\ + 4.37 & 4.14 & 5.75 & 8.74 \\ + 4.18 & 3.96 & 5.50 & 8.36 +\end{bmatrix}. +$$ + +The solutions for $X$ here are allowed to contain negative entries. +There is no good way to constrain them to integer values---rounding +will not guarantee the recovery of exact marginals. Here are a couple +of posterior draws for $x$ given the example marginals and $\sigma = +2.5$. + +$$ +X^{(100)} += +\begin{bmatrix} +12.03 & 11.00 & 12.09 & 19.89 \\ + 2.67 & 4.91 & 5.19 & 10.23 \\ + 4.31 & 2.09 & 7.72 & 7.88 +\end{bmatrix} +\qquad \qquad +X^{(200)} += +\begin{bmatrix} + 8.76 & 8.64 & 14.17 & 23.43 \\ + 3.49 & 6.63 & 5.22 & 7.66 \\ + 6.75 & 2.73 & 5.61 & 6.91 +\end{bmatrix} +$$ + +The posterior standard deviation on the entries is around 1.8, which +is lower than the stipulated 2.5 due to the lowered degrees of +freedom. To be more specific, there are $3 \cdot 4 = 12$ normal +priors when there are only $(3 - 1) \cdot (4 - 1) = 6$ degrees of +freedom due to the sum-to-zero constraint. This is similar to +coding a prior in the following way in Stan, which is unusual +and appears to be a typo, but is legal. + +```stan +parameters { + real y; +} +model { + y ~ normal(0, 1); + y ~ normal(0, 1); +} +``` + +This is equivalent to saying that the prior for $y$ is + +\begin{align*} +p(y) +&\propto \textrm{normal}(y \mid 0, 1) \cdot \textrm{normal}(y \mid 0, 1) +\\[4pt] +&\propto \textrm{normal}(y \mid 0, \sqrt{1 / (1 + 1)}). +\end{align*} + +The second line follows by plugging in the definitions +and reducing. The result with the sum-to-zero constraint +is a doubling of the number of priors in the same way, leading +to a factor of $1 / \sqrt{2}$ lower standard deviation. + + +## Fun shapes + +The title and many of the examples in this section are ported to Stan from +the chapter "Fun Shapes" in the BUGS examples +[@spiegelhalter2003bugs]. + +### The disc + +How could we generate parameters $(x, y)$ that have a uniform distribution +within a disc of radius $r > 0$? This kind of constraint is easy to +parameterize by leaning on Stan's built-in constraints. + +```stan +data { + real r; +} +transformed data { + real r_squared = r^2; +} +parameters { + real x; + real y; +} +``` + +The $x$ value is free to float anywhere in $(-r, r)$. Conditional on +the $x$ value, $y$ must satisfy $x^2 + y^2 < r^2$, or $|y| < \sqrt{r^2 +- x^2}$. Because the constraint involves the parameter `x`, the +combined adjustment does the right thing and assumes that jointly, `x` +and `y` have a uniform distribution within the disc. The scatterplot +of draws shown in @fig-disc-scatterplot verifies that the distribution +is indeed uniform. + +![Draws from the disc distribution of points within a radius $r$ circle for $r=2.5$.](img/disc.png){#fig-disc-scatterplot width=50%} + +The easy way to generate parameters in a disc is to simply scale +the built-in unit vector. This has the nice property of easily +extending to higher dimensions, like points within a sphere and +having built-in change-of-variables adjustments. + +```stan +data { + real r; +} +parameters { + unit_vector[2] xy_unit; +} +transformed parameters { + real x = r * xy_unit[1]; + real y = r * xy_unit[2]; +} +``` + +The hard way to generate parameters within a disc is with an angle and +a radius. We can uniformly generate the angle in $(0, 2\pi)$, then +generate a radius in $(0, r)$. First, we have to map the angle and radius +to $(x, y)$ coordinates by +$$ +x, y = \rho \cdot \cos \theta, \rho \cdot \sin \theta. +$$ +The transform is non-linear, so it requires a change-of-variables +adjustment (absolute value of the determinant of the Jacobian of the +transform). The Jacobian matrix is +$$ +J_{f^{-1}} = \begin{bmatrix} +\cos \theta & - \rho \cdot \sin \theta +\\ +\sin \theta & \rho \cdot \cos \theta +\end{bmatrix}. +$$ +After some algebra, the absolute determinant reduces to $\left|\, \det +J_{f^{-1}} \, \right| = \rho$, so that the Jacobian adjustment on the +log scale is $\log \rho$. + +The Stan model directly encodes the math. + +```stan +data { + real r; +} +parameters { + real rho; + real theta; +} +transformed parameters { + real x = rho * cos(theta); + real y = rho * sin(theta); + jacobian += log(rho); +} +``` + +The plots look indistinguishable for all three ways of parameterizing +the disc. + + +### Hexagons + +A hexagon can be generated in the same way as a circle, by working out +the upper bound function and coding it directly in Stan. + +```stan +functions { + real upper_bound(real x) { + if (x < -0.5) { + return 2 + 2 * x; + } + if (x < 0.5) { + return 1; + } + return 2 - 2 * x; + } +} +parameters { + real x; + real y; +} +``` + +Draws from the hexagon are shown in @fig-hexagon. + +![Uniform draws from a hexagon.](img/hexagon.png){#fig-hexagon width=50%} + +The draws appear to fall off a bit more than they should as $x +\rightarrow -1$ and $x \rightarrow 1$. + + +### Parallelograms + +A parallelogram can be generated by generating uniformly within a +rectangle, such as $x, y^\textrm{raw} \in \textrm{uniform}(0, 1)$ and +then shifting the raw $y$ values by subtracting $x$, $y = +y^\textrm{raw} - x$. + + +```stan +parameters { + real x; + real y_raw; +} +transformed parameters { + real y = y_raw - x; +} +``` + +Draws from the parallelogram are shown in @fig-parallelogram. Draws +from other parallelograms would be similar. For example, to have the +left side lower, add $x$ rather than subtracting it. To have the $x$ +axis be wider, just increase the bounds on $x$. To rotate 90 degrees +so that the flat edges are top and bottom, reverse the definitions of +$x$ and $y$. + +![Uniform draws from a parallelogram.](img/parallelogram.png){#fig-parallelogram width=50%} + +### Hollow square + +The following Stan program defines a uniform distribution on $(0, 1) +\times (0, 1) \setminus (-0.5, 0.5) \times (-0.5, 0.5)$, which looks +like a square with a smaller square taken out of the middle. This one +is trickier to define and the resulting transform is not continuous. +The idea is to take $x^\textrm{raw} \sim \textrm{uniform}(-1.5, 1.5)$ +and $y^\textrm{raw} \sim \textrm{uniform}(-1, 1)$ then transform. If +$x^\textrm{raw} \in (-0.5, 0.5)$, then positive $y^\textrm{raw}$ are +compressed to half their height and shifted up, and the opposite for +positive. To keep the distribution uniform, if $x^\textrm{raw} +\not\in (-0.5, 0.5)$, then it is compressed to half its width inward. +The Stan code follows the description. + +```stan +parameters { + real x_raw; + real y_raw; +} +transformed parameters { + real x = x_raw; + real y = y_raw; + if (x_raw < -0.5) { + x = 0.5 * (x_raw + 0.5) - 0.5; + } else if (x_raw > 0.5) { + x = 0.5 * (x_raw - 0.5) + 0.5; + } else { + y = (x < -0.5 || x > 0.5) + ? y_raw + : (y_raw < 0 ? -0.5 : 0.5) + 0.5 * y_raw; + } +} +``` + +Draws from the hollow square distribution are shown in @fig-hollow-square. + +![Uniform draws from a hollow square.](img/hollow-square.png){#fig-hollow-square width=50%} + +### Ring + +A ring can be parameterized by taking the polar coordinate version of +the disk, but only sampling the radius $\rho \in (0.5, 1)$. + +```stan +parameters { + real rho; + real phi; +} +transformed parameters { + real x = rho * cos(phi); + real y = rho * sin(phi); + jacobian += log(z); +} +``` + +Draws from the ring distribution are shown in @fig-ring. + +![Uniform draws from a ring.](img/ring.png){#fig-ring width=50%} + + +### Soft ring + +Instead of uniform within a ring with a hard boundary, a "soft" +version of a ring distribution can be defined as follows. Rather than +bounding `z`, it is just declared with a `` constraint. By +itself, this would produce an improper distribution because of the +uniformity over an unbounded area. To turn it into a soft ring, `rho` +can be given a gamma distribution with a mean of 1 and a relatively +low variance; a normal would also work. + +```stan +parameters { + real z; + real phi; +} +transformed parameters { + real x = z * cos(phi); + real y = z * sin(phi); + jacobian += log(z); +} +model { + z ~ gamma(20, 20); +} +``` + +Draws from the soft ring distribution are shown in @fig-soft-ring. + +![Non-uniform draws from a soft ring.](img/soft-ring.png){#fig-soft-ring width=50%} +