commonfunc.stan

    // DO NOT EDIT THIS FILE DIRECTLY. It is created by make_commonfunc_stan.py.

    // ========================================================================
    // Common functions
    // ========================================================================
    /*
        Reminders -- Stan's flavour of C++:
        -----------------------------------------------------------------------

        - Disappointingly, you can't modify arguments to Stan user-defined
          functions. (No pass-by-reference.)

        - You can't have templating of user-defined functions, i.e. not this:

            template<T> T somefunc(T x);

        - Two functions with the same name can't even have different
          signatures. So not this:

            real somefunc(real x);
            vector somefunc(vector x);

        - No default values for function parameters. So not this:

            real somefunc(x, y = 0);

        - We can't use a leading "_" prefix on function names (gives a Stan
          syntax error).

        - The addition-assignment (+=) operator generally doesn't work (it
          appears to be reserved for the one variable "target += ...").
          Similarly for all others you might expect.

          - Aha! By Stan 2.19, this has changed. Can use "x += 1"
            (p19 of Stan 2.19 Reference Manual).

        - The ternary (conditional, ?:) operator *is* supported, e.g.:

            x = a ? b : c;

        - Simpler Stan statements (e.g. with the ternary operator) translate
          to fewer C++ statements and faster code (particularly as Stan inserts
          debugging code around the translated C++ statements).

        Reminders -- Stan, other:
        -----------------------------------------------------------------------

        - Array/vector indexing is 1-based.

        - OUTDATED: previously, size() didn't work on a plain "vector" and one
          should have used num_elements(). This is fixed as of Stan ~2.24: see
          https://discourse.mc-stan.org/t/option-to-keep-constant-terms-in-log-probability-via-standard-sampling-syntax/20278/2.
          But remember that size() is "top-level" size (e.g. the first
          dimension of an array), whereas num_elements() counts all elements.

        - Can't define constants in a functions{} block.

    */

    // ------------------------------------------------------------------------
    // Softmax
    // ------------------------------------------------------------------------

    real softmaxNth(vector softmax_inputs, int index)
    {
        /*
            Returns the nth value (at "index") of the softmax of the inputs.
            Assumes an inverse temperature of 1.

            FOR AN EXPLICIT INVERSE TEMPERATURE, see softmaxNthInvTemp().
            FOR A LOGIT (LOG ODDS) VERSION, see logitSoftmaxNth().

            NOTES:

            A softmax function takes several inputs and normalizes them so
            that:
                - the outputs are in the same relative order as the inputs
                - the outputs sum to 1.

            For softmax: see my miscstat.R; the important points for
            optimization are (1) that softmax is invariant to the addition/
            subtraction of a constant, and subtracting the mean makes the
            numbers less likely to fall over computationally; (2) we often only
            need the final part of the computation for a single number
            (preference for one option), so here we don't waste time
            vector-calculating the preference for the left as well [that is:
            we don't have to calculate s_exp_products / sum(s_exp_products)].

            The constant can be the mean, or the max; Stan uses the max, which
            is probably a little more efficient.

            Since Stan 2.0.0, the alternative is to use softmax(); see
            https://github.com/stan-dev/math/blob/develop/stan/math/prim/mat/fun/softmax.hpp
            The exact syntactic equivalence is:

                real result = softmaxNth(inputs, index);  // this
                real result = softmax(inputs)[index];  // Stan

            This "homebrew" version is faster than using Stan's built-in
            softmax() (our speed comparison is in
            rlib/tests/profile_stan_softmax/profile_softmax.stan), presumably
            because Stan calculates the result for all elements of the input,
            and we only bother with the element we care about.
        */

        int length = num_elements(softmax_inputs);
        real constant = max(softmax_inputs);
        vector[length] s_exp_products = exp(softmax_inputs - constant);
        return s_exp_products[index] / sum(s_exp_products);
    }

    real softmaxNthInvTemp(vector softmax_inputs, real inverse_temp, int index)
    {
        /*
            Version of softmaxNth allowing you to specify the inverse temp.

            These are equivalent:

                real result = softmaxNthInvTemp(inputs, invtemp, index);
                real result = softmax(inputs * invtemp)[index];  // Stan

            See softmaxNth() above for speed comparisons.
        */

        return softmaxNth(softmax_inputs * inverse_temp, index);
        // return softmax(softmax_inputs * inverse_temp)[index];
    }

    real logitSoftmaxNth(vector inputs, int index)
    {
        /*
            Returns
                logit(softmax(inputs))[index];
            that is, the log odds for a probability from a softmax function.

            Recall that:

            - odds = p / (1 - p)
            - x = logit(p) = log(odds) = log(p) - log(1 - p) = -log((1/p) - 1)
            - p = logistic(x) = 1 / (1 + exp(-x)) = exp(x) / (exp(x) + 1)
            - softmax(v, i) = exp(v[i]) / sum(exp(v))
            - log_softmax(v, i) = v[i] - log(sum(exp(v))
            - Stan provides log_sum_exp(), log_softmax(), log1m_exp().

            A fully vectorized version in R:

                library(matrixStats)  # for logSumExp
                logitSoftmax <- function(x, debug = FALSE) {
                    log_sum_exp_x <- logSumExp(x)
                    log_p <- x - log_sum_exp_x  # = log(softmax(x))
                    log_1mp = log(1 - exp(log_p))
                    logit <- log_p - log_1mp
                    if (debug) {
                        cat("log_sum_exp_x:\n"); print(log_sum_exp_x)
                        cat("log_p:\n"); print(log_p)
                        p <- exp(log_p)
                        cat("p:\n"); print(p)
                        stopifnot(isTRUE(all.equal(sum(p), 1)))  # check with tolerance
                        cat("log_1mp:\n"); print(log_1mp)
                        cat("logit:\n"); print(logit)
                    }
                    return(logit)
                }
                logitSoftmax(c(1, 2, 3), debug = TRUE)  # demonstration

        */

        // METHOD 1 (fewer calculations involved and empirically faster):
        real log_p = inputs[index] - log_sum_exp(inputs);

        // METHOD 2 (empirically slower):
        // real log_p = log_softmax(inputs)[index];

        // EITHER WAY:
        // Conceptually:
        // (a) log_1mp = log(1 - p)
        //             = log(1 - exp(log_p))
        //             = log1m_exp(log_p)
        // (b) logit   = log(p) - log(1 - p)
        //             = log_p - log_1mp
        // It is very slightly faster (from profiling) to do this in a single
        // step:

        return log_p - log1m_exp(log_p);
    }

    // ------------------------------------------------------------------------
    // Logistic function
    // ------------------------------------------------------------------------

    // - For the logit function, use Stan's built-in logit().
    // - For the standard logistic (with x0 = 0, k = 1, L = 1), use Stan's
    //   inv_logit().

    real logistic(real x, real x0, real k, real L)
    {
        /*
            Returns x transformed through a logistic function, yielding a
            result in the range (0, L).

            Notation as per https://en.wikipedia.org/wiki/Logistic_function:
            - x0: centre
            - k: steepness
            - L: maximum (usually 1)

            The standard logistic function, the inverse of the logit function,
                p = logistic(x) = sigmoid(x) = expit(x) = 1 / (1 + exp(-x))
            where x is a logit (log odds) and p is the resulting probability,
            is a special case where L = 1, k = 1, x0 = 0. However, for that
            you should use Stan's inv_logit().

            Therefore, if you were to transform x so as to be a logit giving
            the same result via the standard logistic function, 1 / (1 +
            exp(-x)), for L = 1, you want this logit:

                x' = k * (x - x0)
        */

        return L / (1 + exp(-k * (x - x0)));
    }

    // ------------------------------------------------------------------------
    // Boundaries (min, max)
    // ------------------------------------------------------------------------

    real bound(real x, real min_value, real max_value)
    {
        // Returns x with minimum/maximum boundaries applied.
        // We would simply do this:
        //     return max(min_value, min(x, max_value));
        // ... but Stan doesn't have max(real, real) or min(real, real) 
        // functions.

        return x < min_value ? min_value : (x > max_value ? max_value : x);
    }

    real boundLower(real x, real min_value)
    {
        // a.k.a. max()

        return x < min_value ? min_value : x;
    }

    real boundUpper(real x, real max_value)
    {
        // a.k.a. min()

        return x > max_value ? max_value : x;
    }

    // ------------------------------------------------------------------------
    // Basic data manipulation
    // ------------------------------------------------------------------------

    vector vector_from_real_array_row(array[,] real x, int row)
    {
        // Given an array
        //      array[nrows, ncols] real x;
        // you can slice the array with
        //      array[ncols] real a = x[row];
        // but not with
        //      vector[ncols] y = x[row];
        // so this function does that.

        int ncols = dims(x)[2];
        vector[ncols] v;
        for (i in 1:ncols) {
            v[i] = x[row, i];
        }
        return v;
    }

    vector vector_from_int_array_row(array[,] int x, int row)
    {
        // As above, but for an int array.

        int ncols = dims(x)[2];
        vector[ncols] v;
        for (i in 1:ncols) {
            v[i] = x[row, i];
        }
        return v;
    }

    vector except_V_V(vector v, int except)
    {
        // Returns a vector that is the original without the element at index
        // "except".

        int n = num_elements(v);
        vector[n - 1] result;
        int r = 1;  // indexes result
        for (i in 1:n) {
            if (i == except) {
                continue;
            }
            result[r] = v[i];
            r += 1;
        }
        return result;
    }

    int except_I_I(int x, int except)
    {
        // The argument is an index to a vector v; the result is the equivalent
        // index to the vector returned by except_V_V(v, except).

        if (x < 1) {
            reject("Argument x is a Stan index so must be >= 1");
        }
        if (except == x) {
            reject("Cannot remove 'except' where except == x");
        }
        if (except < 1 || except > x) {
            return x;
        }
        return x - 1;
    }

    // ------------------------------------------------------------------------
    // Simple functions: matrix calculations
    // ------------------------------------------------------------------------
    // Note that Stan only provides the following versions of dot_product():
    //      dot_product(vector, vector)
    //      dot_product(row vector, row vector)
    //      dot_product(vector, row vector)
    //      dot_product(row vector, vector)
    //      dot_product(array[] real, array[] real)

    vector dot_product_MV_V(matrix x, vector y)
    {
        // Dot product between a matrix (2 dimensions) and a vector (1
        // dimension):
        //
        //      (p, q) matrix ⋅ (q, 1) vector = (p, 1) vector
        //
        // For example:
        //
        //              [a, b]   [g]      [ag + bh]
        //      x ⋅ y = [c, d] ⋅ [h]    = [cg + dh]
        //              [e, f]            [eg + fh]
        //
        //              (3, 2) ⋅ (2, 1) = (3, 1)

        array[2] int x_dimensions = dims(x);
        int p = x_dimensions[1];
        int q = x_dimensions[2];
        vector[p] z;
        real cell;

        if (q != num_elements(y)) {
            reject("Incompatible arguments");
        }
        for (i in 1:p) {  // rows of x
            cell = 0.0;
            for (j in 1:q) {  // columns of x
                cell += x[i, j] * y[j];
            }
            z[i] = cell;
        }
        return z;
    }

    vector dot_product_2A_V(array[,] real x, array[] real y)
    {
        // As dot_product_MV_V, but for array inputs.

        array[2] int x_dimensions = dims(x);
        int p = x_dimensions[1];
        int q = x_dimensions[2];
        vector[p] z;
        real cell;

        if (q != num_elements(y)) {
            reject("Incompatible arguments");
        }
        for (i in 1:p) {  // rows of x
            cell = 0.0;
            for (j in 1:q) {  // columns of x
                cell += x[i, j] * y[j];
            }
            z[i] = cell;
        }
        return z;
    }

    vector dot_product_VM_V(vector x, matrix y)
    {
        // Dot product between a vector (1 dimension) and a matrix (2
        // dimensions):
        //
        //      (1, p) vector ⋅ (p, q) matrix = (1, q) vector
        //
        // For example:
        //
        //                       [a, c, e]
        //      x ⋅ y = [g, h] ⋅ [b, d, f] = [ag + bh, cg + dh, eg + fh]
        //                                 = y' ⋅ x'
        //
        //              (1, 2) ⋅ (2, 3)    = (1, 3) 

        array[2] int y_dimensions = dims(y);
        int p = y_dimensions[1];
        int q = y_dimensions[2];
        vector[q] z;
        real cell;

        if (p != num_elements(x)) {
            reject("Incompatible arguments");
        }
        for (j in 1:q) {  // columns of y
            cell = 0.0;
            for (i in 1:p) {  // rows of y
                cell += x[j] * y[i, j];
            }
            z[j] = cell;
        }
        return z;
    }

    vector dot_product_A2_V(array[] real x, array[,] real y)
    {
        // As dot_product_VM_V(), but for array inputs.

        array[2] int y_dimensions = dims(y);
        int p = y_dimensions[1];
        int q = y_dimensions[2];
        vector[q] z;
        real cell;

        if (p != num_elements(x)) {
            reject("Incompatible arguments");
        }
        for (j in 1:q) {  // columns of y
            cell = 0.0;
            for (i in 1:p) {  // rows of y
                cell += x[j] * y[i, j];
            }
            z[j] = cell;
        }
        return z;
    }

    real dot_product_AA_R(array[] real x, array[] real y)
    {
        // Dot product of two arrays.

        int n = num_elements(x);
        real z = 0.0;
        if (n != num_elements(y)) {
            reject("Incompatible arguments");
        }
        for (i in 1:n) {
            z += x[i] * y[i];
        }
        return z;
    }

    real dot_product_iAV_R(array[] int x, vector y)
    {
        int n = num_elements(x);
        real z = 0.0;
        if (n != num_elements(y)) {
            reject("Incompatible arguments");
        }
        for (i in 1:n) {
            z += x[i] * y[i];
        }
        return z;
    }

    matrix tensordot_A3_M(array[] real x, array[,,] real y)
    {
        // Equivalent to Numpy's tensordot(x, y, axes=1), for:
        //
        //      (1, p) ⋅ (p, q, r) = (q, r)
        //
        // For example:
        //
        //      [a, b] ⋅ [ [c, d, e, f]       = [ac + bc', ad + bd', ...]
        //                 [g, h, i, j]         [ag + bg', ag + bg', ...]
        //                 [k, l, m, n],        [ak + bk', ak + bk', ...]
        //
        //                 [c', d', e', f']
        //                 [g', h', i', j']
        //                 [k', l', m', n'] ]
        //         
        //      (1, 2) ⋅ (2, 3, 4)            = (3, 4)

        array[3] int dimensions = dims(y);
        int p = dimensions[1];
        int q = dimensions[2];
        int r = dimensions[3];
        matrix[q, r] z;
        real cell;

        if (p != num_elements(x)) {
            reject("Incompatible arguments");
        }
        for (j in 1:q) {
            for (k in 1:r) {
                cell = 0.0;
                for (i in 1:p) {
                    cell += x[i] * y[i, j, k];
                }
                z[j, k] = cell;
            }
        }
        return z;
    }

    array[,] real tensordot_A3_2(array[] real x, array[,,] real y)
    {
        // As for tensordot_A3_M(), but returning an array.

        array[3] int dimensions = dims(y);
        int p = dimensions[1];
        int q = dimensions[2];
        int r = dimensions[3];
        array[q, r] real z;
        real cell;

        if (p != num_elements(x)) {
            reject("Incompatible arguments");
        }
        for (j in 1:q) {
            for (k in 1:r) {
                cell = 0.0;
                for (i in 1:p) {
                    cell += x[i] * y[i, j, k];
                }
                z[j, k] = cell;
            }
        }
        return z;
    }

    // ------------------------------------------------------------------------
    // Pairwise differences in matrix format
    // ------------------------------------------------------------------------
    // Two functions with different signatures can't have the same name...

    matrix pairwiseDifferencesSpecifyDiagonal(vector x, vector y, 
                                              real diagonal_value)
    {
        // - Given two vectors of equal length N, returns a matrix[N, N] result
        //   where each element result[i, j] == x[i] - y[j].
        // - Diagonal values, for which i == j, are populated with
        //   diagonal_value. By default this is zero, but if this is to be a
        //   result from e.g. a generated quantities block, Stan will complain
        //   (that the largest value of Rhat is NaN) if diagonal values is unvaryingly
        //   zero. Under those circumstances, you should pass in a small (e.g.
        //   iteration-specific) random number, e.g. like this:
        //          real tiny_random_number = uniform_rng(-1e-16, 1e-16);
        //          group_diffs = pairwiseDifferences(x, y, tiny_random_number);

        int n = num_elements(x);
        matrix[n, n] result;
        real diff_x_minus_y;  // working variable to save a lookup

        if (n != num_elements(y)) {
            reject("Incompatible arguments");
        }
        for (j in 1:n) {  // access matrices in column-major order
            for (i in 1:n) {
                if (i == j) {
                    result[i, j] = diagonal_value;
                } else if (i > j) {
                    // We populate the bottom-left corner [i, j], where i > j,
                    // and simultaneously cross-populate the corresponding cell
                    // in the top-right corner [j, i].
                    diff_x_minus_y = x[i] - y[j];
                    result[i, j] = diff_x_minus_y;
                    result[j, i] = -diff_x_minus_y;
                }
            }
        }
        return result;
    }

    matrix pairwiseDifferences(vector x, vector y)
    {
        // A version of pairwiseDifferences() with diagonal_value = 0.

        return pairwiseDifferencesSpecifyDiagonal(x, y, 0);
    }

    matrix pairwiseDifferencesSelfSpecifyDiagonal(vector x, real diagonal_value)
    {
        // A version of pairwiseDifferences() to compare a vector to itself
        // pairwise.

        return pairwiseDifferencesSpecifyDiagonal(x, x, diagonal_value);
    }

    matrix pairwiseDifferencesSelf(vector x)
    {
        // A version of pairwiseDifferences() to compare a vector to itself
        // pairwise with diagonal_value = 0.

        return pairwiseDifferencesSpecifyDiagonal(x, x, 0);
    }

    // ------------------------------------------------------------------------
    // Pairwise comparisons in vector format
    // ------------------------------------------------------------------------

    int factorial(int x);  // necessary for self-recursion
    int factorial(int x)
    {
        // We could use tgamma(x + 1). But then we run into the unwillingness
        // of Stan to provide functions that round real numbers to integer, and
        // the need for complex workarounds:
        // https://discourse.mc-stan.org/t/real-to-integer-conversion/5622/9 So
        // we could just implement a factorial algorithm; see
        // http://www.luschny.de/math/factorial/FastFactorialFunctions.htm We
        // will just use the simplest version:

        if (x < 0) {
            reject("Factorial undefined for negative numbers. Called for: ", x);
        }
        if (x == 0 || x == 1) {
            return 1;  // 0! = 1, and 1! = 1
        }
        return x * factorial(x - 1);
    }

    int nCombinations(int n, int k)
    {
        // Returns the number of combinations of size k amongst n items.
        //
        // The two-stage approach is entirely because of a wrong warning
        // message from Stan. If you use
        //    return factorial(n) / (factorial(k) * factorial(n - k));
        // then the integer division warning in Stan will print
        //    factorial(n) / factorial(k) * factorial(n - k);
        // ... the removal of the brackets in the warning message may make the
        // reader think the code is wrong.

        int denominator;
        if (n < 1 || k < 1 || n - k < 0) {
            return 0;
        }
        denominator = factorial(k) * factorial(n - k);
        return factorial(n) %/% denominator;
    }

    vector pairwiseDifferencesVec(vector x)
    {
        // Given a vector x of length n (where n > 1), returns a vector of
        // length C(n, 2) of every pairwise comparison.
        //
        // The first pairwise comparisons is x[1] - x[2], then x[1] - x[3],
        // etc., up to x[1] - x[n]. Then it moves to x[2] - x[3], x[2] - x[4],
        // etc. And so on; the last element is x[n - 1] - x[n].
        //
        // The inverse comparisons, e.g. x[2] - x[1], are not performed.

        int n_items = num_elements(x);
        int n_pairs = nCombinations(n_items, 2);
        int pair = 1;
        vector[n_pairs] differences;
        if (n_pairs < 1) {
            return differences;  // empty vector
        }

        for (i in 1:(n_items - 1)) {
            for (j in (i + 1):n_items) {
                differences[pair] = x[i] - x[j];
                pair += 1;
            }
        }
        return differences;
    }

    vector pairwiseDifferencesVecNPairsKnown(vector x, int n_pairs)
    {
        // As for pairwiseDifferencesVec, but with n_pairs precalculated
        // for speed. (The caller will need to have precalculated this to
        // define the size of the return vector...)

        int n_items = num_elements(x);
        int pair = 1;
        vector[n_pairs] differences;
        if (n_pairs < 1) {
            return differences;  // empty vector
        }

        for (i in 1:(n_items - 1)) {
            for (j in (i + 1):n_items) {
                differences[pair] = x[i] - x[j];
                pair += 1;
            }
        }
        return differences;
    }

    // ------------------------------------------------------------------------
    // AUROC (area under the receiver operating characteristic curve)
    // ------------------------------------------------------------------------

    /*

        Calculates AUROC for a binary dependent variable "outcome" from the
        predictor "probability", which is continuous.

        For example, you could use a calculated probability as a predictor,
        or log odds.

        CONCEPT

        See:

        - https://stats.stackexchange.com/questions/145566/how-to-calculate-area-under-the-curve-auc-or-the-c-statistic-by-hand
        - https://www.r-bloggers.com/2016/11/calculating-auc-the-area-under-a-roc-curve/
        - https://blog.revolutionanalytics.com/2016/11/calculating-auc.html

        We will use the following method in principle:

        - For every unique pair of actual values (one is 0, the other is 1):
        - If p_for_outcome_one > p_for_outcome_zero, that's a win (score 1);
          if p_for_outcome_one < p_for_outcome_zero, that's a loss (score 0);
          if p_for_outcome_one = p_for_outcome_zero, that's a tie (score 0.5).
        - Take the mean of those scores; that is the AUROC.

        This follows Hanley & McNeil (1982, PMID 7063747), section III.

        If the outcome doesn't have both ones and zeros, we fail, as in R:
            library(pROC)
            roc(response = c(1, 1, 1, 1), predictor = c(0.1, 0.2, 0.3, 0.4))

        General speedup techniques:
            https://mc-stan.org/docs/2_27/stan-users-guide/vectorization.html

        However, see this algorithm:

        - https://stephanosterburg.gitbook.io/scrapbook/data-science/ds-cheatsheets/machine-learning/fast-computation-of-auc-roc-score

        ALGORITHM

        After:
        - https://stephanosterburg.gitbook.io/scrapbook/data-science/ds-cheatsheets/machine-learning/fast-computation-of-auc-roc-score
        - https://github.com/jfpuget/metrics/blob/master/auc.ipynb

        "Let's first define some entities.

        - pos is the set of examples with target 1. These are the positive
          examples.
        - neg is the set of examples with target 0. These are the negative
          examples.
        - p(i) is the prediction for example i. p(i) is a number between 0
          and 1.
        - A pair of examples (i, j) is labelled the right way if i is a
          positive example, j is a negative example, and the prediction for
          i is higher than the prediction for j.
        - | s | is the number of elements in set s.

        Then AUC-ROC is the count of pairs labelled the right way divided
        by the number of pairs:

            AUC-ROC = | {(i,j), i in pos, j in neg, p(i) > p(j)} | / (| pos | * | neg |)

        A naive code to compute this would be to consider each possible
        pair and count those labelled the right way. A much better way is
        to sort the predictions first, then visit the examples in
        increasing order of predictions. Each time we see a positive
        example we add the number of negative examples we've seen so far."

        ~~~

        RNC: Accuracy verified against R's pROC::roc(); see
        rlib/tests/auroc/test_auroc_algorithm.R.
    */

    real aurocAV(array[] int binary_outcome, vector predictor)
    {
        int n = num_elements(binary_outcome);
        // Sort the binary outcome by ascending predictor:
        array[n] int y = binary_outcome[sort_indices_asc(predictor)];
        int n_false = 0;
        int current_y;
        real total = 0.0;
        for (i in 1:n) {
            current_y = y[i];
            n_false += 1 - current_y;  // add 1 if false; unchanged if true
            total += current_y * n_false;
            // ... if we are seeing a positive example, add the number of
            // negative examples so far.
        }
        return total / (n_false * (n - n_false));
    }

    real aurocAA(array[] int binary_outcome, array[] real predictor)
    {
        // For comments, see aurocAV.
        int n = num_elements(binary_outcome);
        array[n] int y = binary_outcome[sort_indices_asc(predictor)];
        int n_false = 0;
        int current_y;
        real total = 0.0;
        for (i in 1:n) {
            current_y = y[i];
            n_false += 1 - current_y;
            total += current_y * n_false;
        }
        return total / (n_false * (n - n_false));
    }


    // ========================================================================
    // Probability distribution functions not provided by Stan
    // ========================================================================
    // See extra_distribution_functions.stan, which also implements tests.

    /*
        NOTE STAN VERSION/SYNTAX PROBLEMS:

        - "real a, b;" is disallowed by RStan 2.26.1 but permitted by stanc
          2.31.1. The error from the earlier version is:

            ";" or plain assignment expected after variable declaration.

        - gamma_cdf(x, a, b) is warned about in Stan 2.31.1 with a plan to
          remove in Stan 2.32.0. It wants gamma_lcdf(x | a, b). Likewise for
          beta_cdf(). HOWEVER, Stan 2.26.1 prohibits this, with the error:

            Only functions with names ending in _lpdf, _lupdf, _lpmf, _lupmf,
            _lcdf, _lccdf can make use of conditional notation.

          Currently (2023-02-05), we use comma notation, because the latest
          version of RStan is only 2.26.1. This will need changing, though.
    */

    // ------------------------------------------------------------------------
    // qbeta()
    // ------------------------------------------------------------------------

    real qbeta(real p, real alpha, real beta)
    {
        // Quantile, or inverse cumulative distribution function (inverse CDF),
        // for the beta distribution. Equivalent to qbeta() in R, or at least a
        // less capable version of it. Implements the missing Stan function
        // beta_qf().
        //
        // - The first parameter is a cumulative probability.
        // - Distribution parameter shape1 (R) = alpha (Stan).
        // - Distribution parameter shape2 (R) = beta (Stan).
        // - The result is a value from the beta distribution.
        //
        // From
        // https://github.com/SurajGupta/r-source/blob/master/src/nmath/qnbeta.c,
        // modified as per qbeta_notes.txt. We are just implementing the
        // version with lower_tail = false and log_p = false.

        real DBL_EPSILON_X = machine_precision();
        real ONE_M_DBL_EPSILON = 1 - DBL_EPSILON_X;  // precalculate
        real DBL_MIN_ = 1e-323;
        real accu = 1e-15;
        real Eps = 1e-14;  // must be > accu

        real ux;  // Stan 2.26.1 dislikes "real ux, lx, nx, pp;"
        real lx;
        real nx;
        real pp;

        if (p < 0.0 || p > 1.0) {
            reject("qbeta: bad parameter: p < 0 or p > 1");
        }
        if (alpha <= 0.0 || beta <= 0.0) {
            reject("qbeta: bad parameter: alpha <= 0 or beta <= 0");
        }

        // p = R_DT_qIv(p);
        // ... reduces to p for log_p = false and lower_tail = false.

        // Invert pnbeta(.):
        // 1. finding an upper and lower bound
        if (p > ONE_M_DBL_EPSILON) {
            return 1.0;
        }

        // pp = fmin2(ONE_M_DBL_EPSILON, p * (1 + Eps));
        pp = fmin(ONE_M_DBL_EPSILON, p * (1 + Eps));

        // Start ux at 0.5 and work it up (in big steps) while it's too low.
        ux = 0.5;
        while (ux < ONE_M_DBL_EPSILON && beta_cdf(ux, alpha, beta) < pp) {
            ux = 0.5 * (1 + ux);
        }
        // ux is now 0.5 or higher

        pp = p * (1 - Eps);

        // Start lx at 0.5 and work it down (in big steps) while it's too high.
        lx = 0.5;
        while (lx > DBL_MIN_ && beta_cdf(lx, alpha, beta) > pp) {
            lx *= 0.5;
        }
        // lx is now 0.5 or lower

        // 2. interval (lx,ux) halving:
        // Narrow down the gap to find the answer.
        while (1) {
            nx = 0.5 * (lx + ux);  // nx is the mean of lx and ux
            if (beta_cdf(nx, alpha, beta) > p) {
                ux = nx;  // nx too high; move down (shift the upper boundary down)
            } else {
                lx = nx;  // nx too low; move up (shift the lower boundary up)
            }
            if ((ux - lx) / nx <= accu) {
                // Sufficiently accurate!
                break;
            }
        }

        return 0.5 * (ux + lx);
    }

    // ------------------------------------------------------------------------
    // qcauchy()
    // ------------------------------------------------------------------------

    real tanpi(real x)
    {
        // Computes tan(pi * x).
        // R's nmath/nmath.h defines __STDC_WANT_IEC_60559_FUNCS_EXT__ to
        // ensure that tanpi() is defined. There is a declaration in
        // https://github.com/SurajGupta/r-source/blob/master/src/include/Rmath.h0.in
        // Stan has a change request to add tanpi() as of Apr 2021 at
        // https://github.com/stan-dev/math/issues/2376
        //
        // Even Boost's Cauchy distribution uses a simple implementation:
        // https://www.boost.org/doc/libs/1_81_0/boost/random/cauchy_distribution.hpp
        return tan(pi() * x);
    }

    real qcauchy(real p, real location, real scale)
    {
        // location (R) = mu (Stan)
        // scale (R) = sigma (Stan)
        //
        // Implements "lower.tail = TRUE, log.p = FALSE" version.
        // See https://github.com/SurajGupta/r-source/blob/master/src/nmath/qcauchy.c

        int lower_tail = 1;  // no "bool" data type
        real p_working = p;  // can't rewrite argument

        if (p < 0.0 || p > 1.0) {
            reject("qcauchy: bad parameter: p < 0 or p > 1");
        }
        if (scale <= 0) {
            if (scale == 0) {
                // [RNC] Point mass.
                return location;
            }
            reject("qcauchy: bad parameter: sigma (scale) < 0");
        }

        if (p > 0.5) {
            if (p == 1.0) {
                return location
                    + (lower_tail ? scale : -scale) * positive_infinity();
                // RNC: curious! Was the my_INF macro.
            }
            p_working = 1 - p;
            lower_tail = !lower_tail;
        }
        // Use p_working, not p, below here.
        // Stan doesn't let us rewrite our argument.

        if (p_working == 0.5) {
            return location;  // avoid 1/Inf below
        }
        if (p_working == 0.0) {
            return location
                + (lower_tail ? scale : -scale) * negative_infinity();
            // p = 1 is handled above
        }
        return location + (lower_tail ? -scale : scale) / tanpi(p_working);
        // -1/tan(pi * p) = -cot(pi * p) = tan(pi * (p - 1/2))
    }

    // ------------------------------------------------------------------------
    // qgamma(), and its support functions
    // ------------------------------------------------------------------------

    real logcf(real x, real i, real d, real eps)
    {
        // eps: relative tolerance
        // See https://github.com/SurajGupta/r-source/blob/master/src/nmath/pgamma.c
        real c1 = 2 * d;
        real c2 = i + d;
        real c4 = c2 + d;
        real a1 = c2;
        real b1 = i * (c2 - i * x);
        real b2 = d * d * x;
        real a2 = c4 * c2 - b2;
        real scalefactor = pow(2.0, 256.0);
        // = ((4294967296.0^2)^2)^2 = (2^32)^8 = 2^256 = 1.157921e+77
        real c3;

        b2 = c4 * b1 - i * b2;

        while (abs(a2 * b1 - a1 * b2) > abs(eps * b1 * b2)) {
            c3 = c2 * c2 * x;
            c2 += d;
            c4 += d;
            a1 = c4 * a2 - c3 * a1;
            b1 = c4 * b2 - c3 * b1;

            c3 = c1 * c1 * x;
            c1 += d;
            c4 += d;
            a2 = c4 * a1 - c3 * a2;
            b2 = c4 * b1 - c3 * b2;

            if (abs(b2) > scalefactor) {
                a1 /= scalefactor;
                b1 /= scalefactor;
                a2 /= scalefactor;
                b2 /= scalefactor;
            } else if (abs(b2) < 1 / scalefactor) {
                a1 *= scalefactor;
                b1 *= scalefactor;
                a2 *= scalefactor;
                b2 *= scalefactor;
            }
        }

        return a2 / b2;
    }

    real log1pmx(real x)
    {
        // Accurate calculation of log(1+x)-x, particularly for small x.
        // See https://github.com/SurajGupta/r-source/blob/master/src/nmath/pgamma.c
        real minLog1Value = -0.79149064;

        if (x > 1 || x < minLog1Value) {
            return log1p(x) - x;
        } else {
            real r = x / (2 + x);
            real y = r * r;
            if (abs(x) < 1e-2) {
                real two = 2;
                return r * ((((two / 9 * y + two / 7) * y + two / 5) * y +
                        two / 3) * y - x);
            } else {
                real tol_logcf = 1e-14;
                return r * (2 * y * logcf (y, 3, 2, tol_logcf) - x);
            }
        }
    }

    real lgamma1p(real a)
    {
        // Compute  log(gamma(a+1))  accurately also for small a (0 < a < 0.5).
        // See https://github.com/SurajGupta/r-source/blob/master/src/nmath/pgamma.c
        real eulers_const =	 0.5772156649015328606065120900824024;

        // coeffs[i] holds (zeta(i+2)-1)/(i+2) , i = 0:(N-1), N = 40 :
        int N = 40;
        array[N] real coeffs = {
            0.3224670334241132182362075833230126e-0,  // = (zeta(2)-1)/2
            0.6735230105319809513324605383715000e-1,  // = (zeta(3)-1)/3
            0.2058080842778454787900092413529198e-1,
            0.7385551028673985266273097291406834e-2,
            0.2890510330741523285752988298486755e-2,
            0.1192753911703260977113935692828109e-2,
            0.5096695247430424223356548135815582e-3,
            0.2231547584535793797614188036013401e-3,
            0.9945751278180853371459589003190170e-4,
            0.4492623673813314170020750240635786e-4,
            0.2050721277567069155316650397830591e-4,
            0.9439488275268395903987425104415055e-5,
            0.4374866789907487804181793223952411e-5,
            0.2039215753801366236781900709670839e-5,
            0.9551412130407419832857179772951265e-6,
            0.4492469198764566043294290331193655e-6,
            0.2120718480555466586923135901077628e-6,
            0.1004322482396809960872083050053344e-6,
            0.4769810169363980565760193417246730e-7,
            0.2271109460894316491031998116062124e-7,
            0.1083865921489695409107491757968159e-7,
            0.5183475041970046655121248647057669e-8,
            0.2483674543802478317185008663991718e-8,
            0.1192140140586091207442548202774640e-8,
            0.5731367241678862013330194857961011e-9,
            0.2759522885124233145178149692816341e-9,
            0.1330476437424448948149715720858008e-9,
            0.6422964563838100022082448087644648e-10,
            0.3104424774732227276239215783404066e-10,
            0.1502138408075414217093301048780668e-10,
            0.7275974480239079662504549924814047e-11,
            0.3527742476575915083615072228655483e-11,
            0.1711991790559617908601084114443031e-11,
            0.8315385841420284819798357793954418e-12,
            0.4042200525289440065536008957032895e-12,
            0.1966475631096616490411045679010286e-12,
            0.9573630387838555763782200936508615e-13,
            0.4664076026428374224576492565974577e-13,
            0.2273736960065972320633279596737272e-13,
            0.1109139947083452201658320007192334e-13  // = (zeta(40+1)-1)/(40+1)
        };

        real c = 0.2273736845824652515226821577978691e-12;  // zeta(N+2)-1
        real tol_logcf = 1e-14;
        real lgam;
        int i;

        if (abs(a) >= 0.5) {
            // - R C code: lgammafn(x): computes log|gamma(x)|
            //   https://github.com/SurajGupta/r-source/blob/master/src/nmath/lgamma.c
            // - Stan: lgamma(x): natural log of the gamma function applied to x
            //   https://github.com/stan-dev/math/blob/master/stan/math/prim/fun/lgamma.hpp
            return lgamma(a + 1);
        }

        lgam = c * logcf(-a / 2, N + 2, 1, tol_logcf);
        i = N - 1;
        while (i >= 0) {
            lgam = coeffs[i] - a * lgam;
            i -= 1;
        }

        return (a * lgam - eulers_const) * a - log1pmx(a);
    }

    real qchisq_appr(real p, real nu, real g, real tol)
    {
        // An approximation to R's qchisq(p, nu)?
        // See https://github.com/SurajGupta/r-source/blob/master/src/nmath/qgamma.c
        // nu = 'df'
        // g = log Gamma(nu/2) = lgamma(nu/2)
        real C7 = 4.67;
        real C8 = 6.66;
        real C9 = 6.73;
        real C10 = 13.32;
        real M_LN2_ = 0.693147180559945309417232121458;  // ln(2)

        real alpha;
        real a;
        real c;
        real ch;
        real p1;
        real p2;
        real q;
        real t;
        real x;

        // test arguments and initialise

        if (p < 0.0 || p > 1.0) {
            reject("qchisq_appr: bad parameter: p < 0 or p > 1");
        }
        if (nu <= 0.0) {
            reject("qchisq_appr: bad parameter: nu <= 0");
        }

        alpha = 0.5 * nu;  // = [pq]gamma() shape
        c = alpha - 1;

        // p1 = R_DT_log(p);
        // ... R_DT_log(p) reduces to R_D_log(p) for lower_tail = true
        // ... which reduces to log(p) for log_p = false
        p1 = log(p);
        if (nu < -1.24 * p1) {  // for small chi-squared
            // log(alpha) + g = log(alpha) + log(gamma(alpha)) =
            //        = log(alpha*gamma(alpha)) = lgamma(alpha+1) suffers from
            //  catastrophic cancellation when alpha << 1
            real lgam1pa = (alpha < 0.5) ? lgamma1p(alpha) : (log(alpha) + g);
            ch = exp((lgam1pa + p1) / alpha + M_LN2_);

        } else if (nu > 0.32) {  // using Wilson and Hilferty estimate
            // x = qnorm(p, 0, 1, lower_tail, log_p);
            x = inv_Phi(p);
            if (is_inf(x)) {
                // RNC alteration; if p = 1 then x = inf, and then we end up
                // with ch = inf in the next steps, and then via the (ch > 2.2
                // * ...) condition, we get ch = -nan instead of inf. So:
                ch = x;
                // This makes it work.
            } else {

                p1 = 2. / (9 * nu);
                ch = nu * pow(x * sqrt(p1) + 1 - p1, 3);

                // approximation for p tending to 1:
                if (ch > 2.2 * nu + 6) {
                    ch = -2 * (log1p(-p) - c * log(0.5 * ch) + g);
                    // ... R_DT_Clog(p) reduces to R_D_LExp(p) for lower_tail = true
                    // ... R_D_LExp(p) reduces to log1p(-p) for log_p = false
                }
            }

        } else {  // "small nu" : 1.24*(-log(p)) <= nu <= 0.32
            ch = 0.4;
            a = log1p(-p) + g + c * M_LN2_;
            // R_DT_Clog(p) -> log1p(-p), as above
            while (1) {
                q = ch;
                p1 = 1. / (1 + ch * (C7 + ch));
                p2 = ch * (C9 + ch * (C8 + ch));
                t = -0.5 + (C7 + 2 * ch) * p1 - (C9 + ch * (C10 + 3 * ch)) / p2;
                ch -= (1 - exp(a + 0.5 * ch) * p2 * p1) / t;
                if (abs(q - ch) <= tol * abs(ch)) {
                    // converged
                    break;
                }
            }
        }

        return ch;
    }

    real qgamma(real p, real alpha, real beta) {
        // Quantile, or inverse cumulative distribution function (inverse CDF),
        // for the gamma distribution. Equivalent to qgamma() in R, or at least
        // a less capable version of it. Implements the missing Stan function
        // gamma_qf().
        //
        // - The first parameter is a cumulative probability.
        // - Distribution parameter shape (R) = alpha (Stan).
        // - Distribution parameter rate (R) = beta (Stan).
        //   (R also offers scale = 1/rate.)
        // - The result is a value from the gamma distribution.
        //
        // From
        // https://github.com/SurajGupta/r-source/blob/master/src/nmath/qgamma.c,
        // modified as per qbeta_notes.txt. We are just implementing the
        // version with lower_tail = false and log_p = false.

        real scale = 1. / beta;
        real EPS1 = 1e-2;
        real EPS2 = 5e-7;  // final precision of AS 91
        real EPS_N = 1e-15;  // precision of Newton step / iterations
        real LN_EPS = -36.043653389117156;  // = log(.Machine$double.eps) iff IEEE_754
        int MAXIT = 1000;  // was 20
        real pMIN = 1e-100;  // was 0.000002 = 2e-6
        real pMAX = (1 - 1e-14);  // was (1-1e-12) and 0.999998 = 1 - 2e-6
        real i420 = 1. / 420.;
        real i2520 = 1. / 2520.;
        real i5040 = 1. / 5040;
        real p_;
        real a;
        real b;
        real c;
        real g;
        real ch;
        real ch0;
        real p1;
        real p2;
        real q;
        real s1;
        real s2;
        real s3;
        real s4;
        real s5;
        real s6;
        real t;
        real x;
        int max_it_Newton = 1;
        int iterate = 1;  // no boolean type in Stan

        if (p < 0.0 || p > 1.0) {
            reject("qgamma: bad parameter: p < 0 or p > 1");
        }
        if (alpha < 0.0 || scale <= 0.0) {
            reject("qbeta: bad parameter: alpha < 0 or scale (1/rate) <= 0");
        }

        if (alpha == 0.0) {
            // all mass at 0:
            return 0;
        }

        if (alpha < 1e-10) {
            print(
                "qgamma: alpha (", alpha,
                ") is extremely small: results may be unreliable"
            );
            max_it_Newton = 7;  // may still be increased below
        }

        // p_ = R_DT_qIv(p);  // lower_tail prob (in any case)
        p_ = p;

        // g = lgammafn(alpha);  // log Gamma(v/2)
        g = lgamma(alpha);  // log Gamma(v/2)

        // ----- Phase I : Starting Approximation
        ch = qchisq_appr(
            p,
            2 * alpha,  // nu = 'df'
            g,  // = lgamma(nu/2)
            EPS1  // tol
        );
        // if (!R_FINITE(ch)) {
        if (is_inf(ch)) {
            // forget about all iterations!
            max_it_Newton = 0;
            iterate = 0;
        } else if (ch < EPS2) {
            // Corrected according to AS 91; MM, May 25, 1999
            max_it_Newton = 20;
            iterate = 0;  // and do Newton steps
        } else if (p_ > pMAX || p_ < pMIN) {
            max_it_Newton = 20;
            iterate = 0;  // and do Newton steps
        }

        if (iterate) {
            // ----- Phase II: Iteration
            // Call pgamma() [AS 239]	and calculate seven term taylor series
            real M_LN2_ = 0.693147180559945309417232121458;  // ln(2)
            c = alpha - 1;
            s6 = (120 + c * (346 + 127 * c)) * i5040;  // used below, is "const"

            ch0 = ch;  // save initial approx.
            for (i in 1:MAXIT) {
                q = ch;
                p1 = 0.5 * ch;
                p2 = p_ - gamma_cdf(p1, alpha, 1.0);
                // pgamma_raw(p1, alpha, lower_tail=TRUE, log_p=FALSE)
                // is equivalent to pgamma(p1, alpha, scale=1) in R
                // and thus to gamma_cdf(p1, alpha, 1.0).
                // See https://github.com/SurajGupta/r-source/blob/master/src/nmath/pgamma.c

                if (is_inf(p2) || ch <= 0) {
                    ch = ch0;
                    max_it_Newton = 27;
                    break;
                }  // was  return ML_NAN;

                t = p2 * exp(alpha * M_LN2_ + g + p1 - c * log(ch));
                b = t / ch;
                a = 0.5 * t - b * c;
                s1 = (210 + a * (140 + a * (105 + a * (84 + a * (70 + 60 * a)))))
                    * i420;
                s2 = (420 + a * (735 + a * (966 + a * (1141 + 1278 * a)))) * i2520;
                s3 = (210 + a * (462 + a * (707 + 932 * a))) * i2520;
                s4 = (252 + a * (672 + 1182 * a) + c * (294 + a * (889 + 1740 * a)))
                    * i5040;
                s5 = (84 + 2264 * a + c * (1175 + 606 * a)) * i2520;

                ch += t
                    * (1 + 0.5 * t * s1
                       - b * c
                           * (s1 - b * (s2 - b * (s3 - b * (s4 - b * (s5 - b * s6)))))
                    );
                if (abs(q - ch) < EPS2 * ch) {
                    break;
                }
                if (abs(q - ch) > 0.1 * ch) {
                    // diverging? -- also forces ch > 0
                    if (ch < q) {
                        ch = 0.9 * q;
                    } else {
                        ch = 1.1 * q;
                    }
                }
            }
            // no convergence in MAXIT iterations -- but we add Newton now...
        }

        x = 0.5 * scale * ch;
        if (max_it_Newton) {
            // always use log scale
            real NEG_INF = negative_infinity();
            real DBL_MIN_ = 1e-323;  // as above
            real logp = log(p);  // can't reassign to p
            if (x == 0) {
                real one_p = 1. + 1e-7;
                real one_m = 1. - 1e-7;
                x = DBL_MIN_;
                p_ = gamma_lcdf(x | alpha, beta);
                if (p_ > logp * one_p) {
                    return 0;
                }
                // else:  continue, using x = DBL_MIN instead of 0
            } else {
                p_ = gamma_lcdf(x | alpha, beta);
            }
            if (p_ == NEG_INF) {
                return 0.0;  // PR#14710
            }
            for (i in 1:max_it_Newton) {
                p1 = p_ - logp;
                if (abs(p1) < abs(EPS_N * logp)) {
                    break;
                }
                g = gamma_lpdf(x | alpha, beta);
                if (g == NEG_INF) {
                    break;
                }
                t = p1 * exp(p_ - g);  // = "delta x"
                t = x - t;
                p_ = gamma_cdf(t, alpha, beta);
                if (abs(p_ - logp) > abs(p1)
                        || (i > 1 && abs(p_ - logp) == abs(p1))) {
                    // no improvement
                    break;
                }
                x = t;
            }
        }

        return x;
    }

    // ------------------------------------------------------------------------
    // qwiener(), as in R's RWiener package
    // ------------------------------------------------------------------------
    // Notes:
    // - Stan's ceil() does not return an integer. One might think to_int()
    //   would solve the problem, but it needs an input qualified as "data".
    //   If you mark all function arguments "data real" rather than real,
    //   values calculated from them are not automatically considered "data";
    //   nor can you declare a local variable as "data real". This looks like
    //   nannying on the part of Stan and is unhelpful.
    //   We need a to_int_stupid() function.
    // - pnorm(q, mean = 0, sd = 1, lower_tail = 1, log.p = 0) becomes
    //   std_normal_cdf(q).
    // - qnorm(p, 0, 1, 1, 0) becomes std_normal_qf(q) from Stan 2.31.1 or
    //   inv_Phi(p) before that.
    // - sign() requires implementing

    int to_positive_int(real x) {
        // https://discourse.mc-stan.org/t/real-to-integer-conversion/5622/5
        // This assumes that min_val >= 0 is the minimum integer in range,
        // max_val > min_val, and that x has already been rounded. It should
        // find the integer equivalent to x.
        // RNC: By TRUNCATION.
        int min_val = 0;
        int max_val = 2147483646;  // 2 ^ 31 - 2;
        // ... https://mc-stan.org/docs/2_18/reference-manual/numerical-data-types.html
        int range;
        int mid_pt;
        int out;
        int n_for_simple_version = 20;

        if (x < min_val || x > max_val) {
            reject(
                "to_positive_int: ", x,
                " out of range [", min_val, ", ", max_val, "]"
            );
        }

        // The "halving range" method below takes up to about 30-31 interations
        // for our possible range, so is pretty efficient. But for some small
        // number, this is more efficient:
        if (x - min_val < n_for_simple_version) {
            for (q in min_val:min_val + n_for_simple_version) {
                if (q > x) {  // e.g. 6 > 5.1
                    return q - 1;  // e.g. return 5
                }
            }
        }

        range = (max_val - min_val + 1) %/% 2;
        // We add 1 to make sure that truncation doesn't exclude a number.
        // ... and %/% for integer division in Stan, to avoid a warning.
        mid_pt = min_val + range;
        while (range > 0) {
            if (x >= mid_pt && x < mid_pt + 1) {
                out = mid_pt;
                range = 0;
            } else {
                // figure out if range == 1
                range = (range + 1) %/% 2;  // basically, halve the range
                mid_pt = x > mid_pt ? mid_pt + range : mid_pt - range;
            }
        }
        return out;
    }

    int sign(real x) {
        // https://discourse.mc-stan.org/t/signum-function/9156
        // https://github.com/stan-dev/math/issues/1686
        return x < 0 ? -1 : x > 0;
    }

    real prob_upperbound(real v, real a, real w)
    {
        // Used by pwiener. See
        // https://github.com/cran/RWiener/blob/master/src/pwiener.c
        real e = exp(-2.0 * v * a * (1.0 - w));
        if (is_inf(e)) {
            return 1;
        }
        if (v == 0 || w == 1) {
            return 1 - w;
        }
        return ((1 - e) / (exp(2.0 * v * a * w) - e));
    }

    real exp_pnorm(real a, real b)
    {
        // Used by pwiener. See
        // https://github.com/cran/RWiener/blob/master/src/pwiener.c
        // But the original had some daft tests for impossible conditions.
        return exp(a) * std_normal_cdf(b);
    }

    int K_large(real q, real v, real a, real w)
    {
        // Used by pwiener. See
        // https://github.com/cran/RWiener/blob/master/src/pwiener.c
        real err = 1e-10;
        real PI = pi();
        real sqrtL1 = sqrt(1 / q) * a / PI;
        real sqrtL2 = sqrt(fmax(
            1.0,
            -2 / q * a * a / PI / PI * (
                log(err * PI * q / 2 * (v * v + PI * PI / a / a))
                + v * a * w + v * v * q / 2
            )
        ));
        return to_positive_int(ceil(fmax(sqrtL1, sqrtL2)));
    }

    int K_small(real q, real v, real a, real w, real epsilon)
    {
        // Used by pwiener. See
        // https://github.com/cran/RWiener/blob/master/src/pwiener.c
        real S2;
        real S3;
        real S4;
        int K;
        if (v == 0) {
            return to_positive_int(ceil(
                fmax(
                    0.0,
                    w / 2 - sqrt(q) / 2 / a * inv_Phi(
                        fmax(0.0, fmin(1.0, epsilon / (2 - 2 * w)))
                    )
                )
            ));
        }
        if (v > 0) {
            return K_small(q, -v, a, w, exp(-2 * a * w * v) * epsilon);
        }
        S2 = w - 1 + 0.5 / v / a * log(epsilon / 2 * (1 - exp(2 * v * a)));
        S3 = (0.535 * sqrt(2 * q) + v * q + a * w) / 2 / a;
        S4 = w / 2 - sqrt(q) / 2 / a * inv_Phi(
            fmax(
                0.0,
                fmin(
                    1.0,
                    epsilon * a / 0.3 / sqrt(2 * pi() * q)
                        * exp(v * v * q / 2 + v * a * w)
                )
            )
        );
        return to_positive_int(ceil(fmax(fmax(fmax(S2, S3), S4), 0.0)));
    }

    real Fl_lower(real q, real v, real a, real w, int K)
    {
        // Used by pwiener. See
        // https://github.com/cran/RWiener/blob/master/src/pwiener.c
        real F = 0;
        real PI = pi();
        int k = K;
        while (k >= 1) {
            F = F - k
                / (v * v * 1.0 + k * k * PI * PI / (a * 1.0) / a)
                * exp(
                    -v * a * w * 1.0
                    - 0.5 * v * v * q
                    - 0.5 * k * k * PI * PI / (a * 1.0) / a * q
                )
                * sin(PI * k * w);
            // RNC: What's with all the "* 1.0"? Pointless, or something fancy
            // about floating-point handling?
            k -= 1;
        }
        return prob_upperbound(v, a, w) + 2.0 * PI / (a * 1.0) / a * F;
    }

    real Fs0_lower(real q, real a, real w, int K)
    {
        // Used by pwiener. See
        // https://github.com/cran/RWiener/blob/master/src/pwiener.c
        real F = 0;
        int k = K;
        while (k >= 0) {
            F = F
                - std_normal_cdf((-2 * k - 2 + w) * a / sqrt(q))
                + std_normal_cdf((-2 * k - w) * a / sqrt(q));
            k -= 1;
        }
        return 2 * F;
    }

    real Fs_lower(real q, real v, real a, real w, int K)
    {
        // Used by pwiener. See
        // https://github.com/cran/RWiener/blob/master/src/pwiener.c
        real S1 = 0;
        real S2 = 0;
        real sqt = sqrt(q);
        int k = K;
        if (v == 0) {
            return Fs0_lower(q, a, w, K);
        }
        while (k >= 1) {
            S1 = S1
                + exp_pnorm(
                    2 * v * a * k,
                    -sign(v) * (2 * a * k + a * w + v * q) / sqt
                )
                - exp_pnorm(
                    -2 * v * a * k - 2 * v * a * w,
                    sign(v) * (2 * a * k + a * w - v * q) / sqt
                );
            S2 = S2
                + exp_pnorm(
                    -2 * v * a * k,
                    sign(v) * (2 * a * k - a * w - v * q) / sqt
                )
                - exp_pnorm(
                    2 * v * a * k - 2 * v * a * w,
                    -sign(v) * (2 * a * k - a * w + v * q) / sqt
                );
            k -= 1;
        }
        return prob_upperbound(v, a, w)
            + sign(v) * (
                (
                    std_normal_cdf(-sign(v) * (a * w + v * q) / sqt)
                    - exp_pnorm(
                        -2 * v * a * w,
                        sign(v) * (a * w - v * q) / sqt
                    )
                )
                + S1 + S2
            );
    }

    real F_lower(real q, real v, real a, real w)
    {
        // Used by pwiener. See
        // https://github.com/cran/RWiener/blob/master/src/pwiener.c
        real err = 1e-10;
        int K_l = K_large(q, v, a, w);
        int K_s = K_small(q, v, a, w, err);
        if (K_l < 10 * K_s) {
            return Fl_lower(q, v, a, w, K_l);
        } else {
            return Fs_lower(q, v, a, w, K_s);
        }
    }

    real pwiener(real q, real alpha, real tau, real beta, real delta)
    {
        // https://github.com/cran/RWiener/blob/master/src/pwiener.c
        if (is_inf(q)) {
            return positive_infinity();
        }
        if (is_nan(q)) {
            return not_a_number();
        }
        if (abs(q) <= tau) {
            return 0;
        }

        if (q < 0) { // lower boundary 0
            return F_lower(abs(q) - tau, delta, alpha, beta);
        } else {  // upper boundary a
            return F_lower(q - tau, (-delta), alpha, (1 - beta));
        }
    }

    real qwiener(real p, real alpha, real tau, real beta, real delta)
    {
        // Based on https://github.com/cran/RWiener/blob/master/src/qwiener.c
        // Its license: GPL v2
        //
        // R's qwiener() handles upper, lower, or both boundaries.
        // Stan's only does the upper boundary, so we can just implement that.
        // https://mc-stan.org/docs/functions-reference/wiener-first-passage-time-distribution.html
        //
        // R's qwiener() calls C function qwiener_c, which is likely qwiener()
        // in qwiener.c. That is boilerplate code around qwiener_d().
        // We'll restrict it to p in [0, 1]; the original dealt with negative
        // probabilities too. And pmin/pmax not used.
        //
        // For pwiener(), Stan does not provide wiener_cdf, so we have to
        // implement that too.

        real pmid = 0;  // value of p for current q guess
        real qmin = 0;  // lower bound of q
        real qmax = positive_infinity();  // upper bound of q
        real q = 1;  // current guess for q
        int c = 0;  // loop counter

        if (p < 0.0 || p > 1.0) {
            reject("qwiener: bad parameter: p < 0 or p > 1");
        }

        while (1) {
            c += 1;
            pmid = pwiener(q, alpha, tau, beta, delta);
            if (p <= pmid) {  // near lower point
                // RNC: pmid correct or too high; move q down
                qmax = q;
                q = qmin + (qmax - qmin) / 2;
            } else {  // near upper point
                // RNC: pmid too low; move q up
                qmin = q;
                if (!is_inf(qmax)) {
                    q = qmin + (qmax - qmin) / 2;
                } else {
                    q = q * 10;
                }
            }
            if (is_nan(pmid)) {
                return not_a_number();
            }
            if (q >= 1e+10) {
                return positive_infinity();
            }
            if (abs(p - pmid) <= 1e-10 || c >= 1000) {
                // accurate enough, or enough iterations
                break;
            }
        }
        return q;
    }

    // ------------------------------------------------------------------------
    // Upper-half quantile functions
    // ------------------------------------------------------------------------
    // Return the quantile of a random variate from a HALF-distribution, given
    // the corresponding cumulative probability. The transformation is p/2 +
    // 0.5, which compresses p values from the range [0, 1] to [0.5, 1].

    real qupperhalfnormal(real p, real mu, real sigma)
    {
        return inv_Phi(p / 2.0 + 0.5) * sigma + mu;
    }

    real qupperhalfcauchy(real p, real mu, real sigma)
    {
        return qcauchy(p / 2.0 + 0.5, mu, sigma);
    }


    // ------------------------------------------------------------------------
    // LOG PROBABILITY FUNCTIONS FOR BRIDGE SAMPLING
    // ------------------------------------------------------------------------
    /*
    We can have functions that access the log probability accumulator
    if the function name ends in '_lp'; see Stan manual section 23.3.

    RE ARGUMENTS:

    The Stan manual uses notation like
         real normal_lpdf(reals y | reals mu, reals sigma)
    but "reals" isn't something you can actually use in user functions.
    See p495:
        "reals" means:
                real
                array[] real, formerly called real[]
                vector
                row_vector
        "ints" means
                int
                array[] int, formerly called int[]

    Moreover, you can't define two copies of the same function with
    different names (23.6: no overloading of user-defined functions).
    For real arguments, the options are therefore:
         real
         array[] real  // one-dimensional array, formerly real[]
         array[,] real  // two-dimensional array, formerly real[,]
         array[,,] real  // three-dimensional array (... etc.)
         vector  // vector, similar to a one-dimensional array.
         matrix  // matrix, similar to a two-dimensional array.
    See p297 of the 2017 Stan manual, and also p319.
    Which do we use in practice?
    - Firstly, we use single numbers or one-dimensional collections,
      and generally the latter. So that means array[] real or vector.
    - We use both.
    - So let's have "Real", "Arr" and "Vec" versions.
    - Then, to make things worse, we sometimes have constant parameters,
      and sometimes array/vector parameters...
    - For something with two distribution parameters, like the normal
      distribution and many others, that means that we have 3*3*3 combinations
      for each thing. Urgh. Stan should allow user overloading ;).
    - Let's do it and define "R", "A", "2", "3", "V" for the parameters.
      (Also "M" for matrix.)
    - Except we won't be returning R unless it's RRR!
    - Last thing cycles fastest.
    So:
        RRR
        -- nothing else R*

        ARA
        ARV
        AAR
        AAA
        AAV
        AVR
        AVA
        AVV

        2RR
            ...

        3RR
            ...

        VRA
        VRV
        VAR
        VAA
        VAV
        VVR
        VVA
        VVV

    RE SAMPLING TWO-DIMENSIONAL ARRAYS:

    You can't sample an entire matrix or 2D array; you have do to it row-wise.
    - This isn't very clear in the manual, as far as I can see.
    - The definition of e.g. beta_lpdf() is in terms of "reals", which
      probably means a vector or array of real.
    - Section 9.6 ("Multi-logit regression") of the Stan manual v2.16.0
      gives an example where one would use a matrix sampling statement but
      they don't.
    - But it is explicit in the sense that they define what they mean by
      "reals", as above, and that doesn't include 2D arrays.
    - Better to move the boilerplate code here than in user land, though.

    RE TWO-DIMENSIONAL ARRAYS:

        // real thing[N_A, N_B];  // old Stan syntax
        array[N_A, N_B] real thing;  // new Stan syntax, from v2.26

        // One way to iterate through all elements:
        for (a in 1:N_A) {
            for (b in 1:N_B) {
                do_something(thing[a, b]);
            }
        }

        // NOT another way to iterate through all elements:
        for (i in 1:num_elements(thing)) {
            do_something(thing[i]);  // thing[i] is an array[] real, not a real
            // ... and thing[num_elements(thing)] will be an index overflow
        }

    So for some functions we want array[,] real... let's give this the
    one-character notation "2" (for 2D array).

    Now:
        num_elements() gives the total, in this case N_A * N_B;
            ... but when *accessing* a 2D array, my_array[1] gives the first
                row, not the first element; see Stan 2017 manual p323.
        size() gives the size of first dimension, in this case N_A;
        dims() gives all dimensions, in this case an array[] int containing {N_A, N_B}.

    RE ARITHMETIC:

    Note that we cannot do:
            real * array[] real
            vector * vector

    */

    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    // Helper functions for boundary checking
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    // For the use of "target += negative_infinity()" with truncation, see Stan
    // (2017, v2.16.0) manual p82, or e.g.
    // https://mc-stan.org/docs/2_25/reference-manual/sampling-statements-section.html#truncation-with-lower-and-upper-bounds-in-stan
    //
    // These are internal functions that ASSUME size match.

    // Lower

    void enforceLowerBound_R_lp(real y, real lowerbound)
    {
        if (y < lowerbound) {
            target += negative_infinity();
        }
    }
    void enforceLowerBound_A_lp(array[] real y, real lowerbound)
    {
        for (i in 1:num_elements(y)) {
            if (y[i] < lowerbound) {
                target += negative_infinity();
                return;
            }
        }
    }
    void enforceLowerBound_2_lp(array[,] real y, real lowerbound)
    {
        array[2] int dimensions = dims(y);
        for (i in 1:dimensions[1]) {
            for (j in 1:dimensions[2]) {
                if (y[i, j] < lowerbound) {
                    target += negative_infinity();
                    return;
                }
            }
        }
    }
    void enforceLowerBound_3_lp(array[,,] real y, real lowerbound)
    {
        array[3] int dimensions = dims(y);
        for (i in 1:dimensions[1]) {
            for (j in 1:dimensions[2]) {
                for (k in 1:dimensions[3]) {
                    if (y[i, j, k] < lowerbound) {
                        target += negative_infinity();
                        return;
                    }
                }
            }
        }
    }
    void enforceLowerBound_V_lp(vector y, real lowerbound)
    {
        for (i in 1:num_elements(y)) {
            if (y[i] < lowerbound) {
                target += negative_infinity();
                return;
            }
        }
    }

    // Upper

    void enforceUpperBound_R_lp(real y, real upperbound)
    {
        if (y > upperbound) {
            target += negative_infinity();
        }
    }
    void enforceUpperBound_A_lp(array[] real y, real upperbound)
    {
        for (i in 1:num_elements(y)) {
            if (y[i] > upperbound) {
                target += negative_infinity();
                return;
            }
        }
    }
    void enforceUpperBound_2_lp(array[,] real y, real upperbound)
    {
        array[2] int dimensions = dims(y);
        for (i in 1:dimensions[1]) {
            for (j in 1:dimensions[2]) {
                if (y[i, j] > upperbound) {
                    target += negative_infinity();
                    return;
                }
            }
        }
    }
    void enforceUpperBound_3_lp(array[,,] real y, real upperbound)
    {
        array[3] int dimensions = dims(y);
        for (i in 1:dimensions[1]) {
            for (j in 1:dimensions[2]) {
                for (k in 1:dimensions[3]) {
                    if (y[i, j, k] > upperbound) {
                        target += negative_infinity();
                        return;
                    }
                }
            }
        }
    }
    void enforceUpperBound_V_lp(vector y, real upperbound)
    {
        for (i in 1:num_elements(y)) {
            if (y[i] > upperbound) {
                target += negative_infinity();
                return;
            }
        }
    }

    // Range

    void enforceRangeBounds_R_lp(real y, real lowerbound, real upperbound)
    {
        if (y < lowerbound || y > upperbound) {
            target += negative_infinity();
        }
    }
    void enforceRangeBounds_A_lp(array[] real y, real lowerbound, real upperbound)
    {
        real value;
        for (i in 1:num_elements(y)) {
            value = y[i];  // lookup only once
            if (value < lowerbound || value > upperbound) {
                target += negative_infinity();
                return;
            }
        }
    }
    void enforceRangeBounds_2_lp(array[,] real y, real lowerbound, real upperbound)
    {
        array[2] int dimensions = dims(y);
        real value;
        for (i in 1:dimensions[1]) {
            for (j in 1:dimensions[2]) {
                value = y[i, j];  // lookup only once
                if (value < lowerbound || value > upperbound) {
                    target += negative_infinity();
                    return;
                }
            }
        }
    }
    void enforceRangeBounds_3_lp(array[,,] real y, real lowerbound, real upperbound)
    {
        array[3] int dimensions = dims(y);
        real value;
        for (i in 1:dimensions[1]) {
            for (j in 1:dimensions[2]) {
                for (k in 1:dimensions[3]) {
                    value = y[i, j, k];  // lookup only once
                    if (value < lowerbound || value > upperbound) {
                        target += negative_infinity();
                        return;
                    }
                }
            }
        }
    }
    void enforceRangeBounds_V_lp(vector y, real lowerbound, real upperbound)
    {
        real value;
        for (i in 1:num_elements(y)) {
            value = y[i];  // lookup only once
            if (value < lowerbound || value > upperbound) {
                target += negative_infinity();
                return;
            }
        }
    }

    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    // Normal distribution
    // - mu is the mean; sigma is the standard deviation
    // - See Stan 2.16.0 manual p512;
    //   http://mathworld.wolfram.com/NormalDistribution.html
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    
    // Sampling

    void sampleNormal_RRR_lp(real y, real mu, real sigma)
    {
        target += normal_lpdf(y | mu, sigma);
    }

    void sampleNormal_ARR_lp(array[] real y, real mu, real sigma)
    {
        target += normal_lpdf(y | mu, sigma);
    }

    void sampleNormal_2RR_lp(array[,] real y, real mu, real sigma)
    {
        for (i in 1:size(y)) {
            target += normal_lpdf(y[i] | mu, sigma);
            // ... y[i] is a one-dimensional array
        }
    }

    void sampleNormal_3RR_lp(array[,,] real y, real mu, real sigma)
    {
        array[3] int dimensions = dims(y);
        for (i in 1:dimensions[1]) {
            for (j in 1:dimensions[2]) {
                for (k in 1:dimensions[3]) {
                    target += normal_lpdf(y[i, j, k] | mu, sigma);
                }
            }
        }
    }

    void sampleNormal_VRR_lp(vector y, real mu, real sigma)
    {
        target += normal_lpdf(y | mu, sigma);
    }

    // Sampling with lower bound

    void sampleNormalLowerBound_RRR_lp(real y, real mu, real sigma, real lowerbound)
    {
        real correction_per_value = normal_lccdf(lowerbound | mu, sigma);
        target += normal_lpdf(y | mu, sigma) -
                  correction_per_value;
        enforceLowerBound_R_lp(y, lowerbound);
    }

    void sampleNormalLowerBound_ARR_lp(array[] real y, real mu, real sigma, real lowerbound)
    {
        real correction_per_value = normal_lccdf(lowerbound | mu, sigma);
        target += normal_lpdf(y | mu, sigma) -
                  correction_per_value * num_elements(y);
        enforceLowerBound_A_lp(y, lowerbound);
    }

    void sampleNormalLowerBound_2RR_lp(array[,] real y, real mu, real sigma, real lowerbound)
    {
        real correction_per_value = normal_lccdf(lowerbound | mu, sigma);
        array[2] int dimensions = dims(y);
        real correction_per_row = correction_per_value * dimensions[2];
        for (i in 1:dimensions[1]) {
            target += normal_lpdf(y[i] | mu, sigma) -
                      correction_per_row;
            // ... y[i] is a one-dimensional array
        }
        enforceLowerBound_2_lp(y, lowerbound);
    }

    void sampleNormalLowerBound_3RR_lp(array[,,] real y, real mu, real sigma, real lowerbound)
    {
        real correction_per_value = normal_lccdf(lowerbound | mu, sigma);
        array[3] int dimensions = dims(y);
        for (i in 1:dimensions[1]) {
            for (j in 1:dimensions[2]) {
                for (k in 1:dimensions[3]) {
                    target += normal_lpdf(y[i, j, k] | mu, sigma) -
                              correction_per_value;
                }
            }
        }
        enforceLowerBound_3_lp(y, lowerbound);
    }

    void sampleNormalLowerBound_VRR_lp(vector y, real mu, real sigma, real lowerbound)
    {
        real correction_per_value = normal_lccdf(lowerbound | mu, sigma);
        target += normal_lpdf(y | mu, sigma) -
                  correction_per_value * num_elements(y);
        enforceLowerBound_V_lp(y, lowerbound);
    }

    // Sampling with upper bound

    void sampleNormalUpperBound_RRR_lp(real y, real mu, real sigma, real upperbound)
    {
        real correction_per_value = normal_lcdf(upperbound | mu, sigma);
        target += normal_lpdf(y | mu, sigma) -
                  correction_per_value;
        enforceUpperBound_R_lp(y, upperbound);
    }

    void sampleNormalUpperBound_ARR_lp(array[] real y, real mu, real sigma, real upperbound)
    {
        real correction_per_value = normal_lcdf(upperbound | mu, sigma);
        target += normal_lpdf(y | mu, sigma) -
                  correction_per_value * num_elements(y);
        enforceUpperBound_A_lp(y, upperbound);
    }

    void sampleNormalUpperBound_2RR_lp(array[,] real y, real mu, real sigma, real upperbound)
    {
        real correction_per_value = normal_lcdf(upperbound | mu, sigma);
        array[2] int dimensions = dims(y);
        real correction_per_row = correction_per_value * dimensions[2];
        for (i in 1:dimensions[1]) {
            target += normal_lpdf(y[i] | mu, sigma) -
                      correction_per_row;
            // ... y[i] is a one-dimensional array
        }
        enforceUpperBound_2_lp(y, upperbound);
    }

    void sampleNormalUpperBound_3RR_lp(array[,,] real y, real mu, real sigma, real upperbound)
    {
        real correction_per_value = normal_lcdf(upperbound | mu, sigma);
        array[3] int dimensions = dims(y);
        for (i in 1:dimensions[1]) {
            for (j in 1:dimensions[2]) {
                for (k in 1:dimensions[3]) {
                    target += normal_lpdf(y[i, j, k] | mu, sigma) -
                              correction_per_value;
                }
            }
        }
        enforceUpperBound_3_lp(y, upperbound);
    }

    void sampleNormalUpperBound_VRR_lp(vector y, real mu, real sigma, real upperbound)
    {
        real correction_per_value = normal_lcdf(upperbound | mu, sigma);
        target += normal_lpdf(y | mu, sigma) -
                  correction_per_value * num_elements(y);
        enforceUpperBound_V_lp(y, upperbound);
    }

    // Sampling with range (lower and upper) bounds

    void sampleNormalRangeBound_RRR_lp(real y, real mu, real sigma, real lowerbound, real upperbound)
    {
        real correction_per_value = log_diff_exp(
            normal_lcdf(upperbound | mu, sigma),
            normal_lcdf(lowerbound | mu, sigma));
        target += normal_lpdf(y | mu, sigma) -
                  correction_per_value;
        enforceRangeBounds_R_lp(y, lowerbound, upperbound);
    }

    void sampleNormalRangeBound_ARR_lp(array[] real y, real mu, real sigma, real lowerbound, real upperbound)
    {
        real correction_per_value = log_diff_exp(
            normal_lcdf(upperbound | mu, sigma),
            normal_lcdf(lowerbound | mu, sigma));
        target += normal_lpdf(y | mu, sigma) -
                  correction_per_value * num_elements(y);
        enforceRangeBounds_A_lp(y, lowerbound, upperbound);
    }

    void sampleNormalRangeBound_2RR_lp(array[,] real y, real mu, real sigma, real lowerbound, real upperbound)
    {
        real correction_per_value = log_diff_exp(
            normal_lcdf(upperbound | mu, sigma),
            normal_lcdf(lowerbound | mu, sigma));
        array[2] int dimensions = dims(y);
        real correction_per_row = correction_per_value * dimensions[2];
        for (i in 1:dimensions[1]) {
            target += normal_lpdf(y[i] | mu, sigma) -
                      correction_per_row;
            // ... y[i] is a one-dimensional array
        }
        enforceRangeBounds_2_lp(y, lowerbound, upperbound);
    }

    void sampleNormalRangeBound_3RR_lp(array[,,] real y, real mu, real sigma, real lowerbound, real upperbound)
    {
        real correction_per_value = log_diff_exp(
            normal_lcdf(upperbound | mu, sigma),
            normal_lcdf(lowerbound | mu, sigma));
        array[3] int dimensions = dims(y);
        for (i in 1:dimensions[1]) {
            for (j in 1:dimensions[2]) {
                for (k in 1:dimensions[3]) {
                    target += normal_lpdf(y[i, j, k] | mu, sigma) -
                              correction_per_value;
                }
            }
        }
        enforceRangeBounds_3_lp(y, lowerbound, upperbound);
    }

    void sampleNormalRangeBound_VRR_lp(vector y, real mu, real sigma, real lowerbound, real upperbound)
    {
        real correction_per_value = log_diff_exp(
            normal_lcdf(upperbound | mu, sigma),
            normal_lcdf(lowerbound | mu, sigma));
        target += normal_lpdf(y | mu, sigma) -
                  correction_per_value * num_elements(y);
        enforceRangeBounds_V_lp(y, lowerbound, upperbound);
    }

    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    // Standard normal distribution, N(0,1)
    // - Note that we have to use normal_lcdf (etc.) but can use
    //   std_normal_lpdf.
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    
    // Sampling

    void sampleStdNormal_R_lp(real y)
    {
        target += std_normal_lpdf(y);
    }

    void sampleStdNormal_A_lp(array[] real y)
    {
        target += std_normal_lpdf(y);
    }

    void sampleStdNormal_2_lp(array[,] real y)
    {
        for (i in 1:size(y)) {
            target += std_normal_lpdf(y[i]);
            // ... y[i] is a one-dimensional array
        }
    }

    void sampleStdNormal_3_lp(array[,,] real y)
    {
        array[3] int dimensions = dims(y);
        for (i in 1:dimensions[1]) {
            for (j in 1:dimensions[2]) {
                for (k in 1:dimensions[3]) {
                    target += std_normal_lpdf(y[i, j, k]);
                }
            }
        }
    }

    void sampleStdNormal_V_lp(vector y)
    {
        target += std_normal_lpdf(y);
    }

    // Sampling with lower bound

    void sampleStdNormalLowerBound_R_lp(real y, real lowerbound)
    {
        real correction_per_value = normal_lccdf(lowerbound | 0, 1);
        target += std_normal_lpdf(y) -
                  correction_per_value;
        enforceLowerBound_R_lp(y, lowerbound);
    }

    void sampleStdNormalLowerBound_A_lp(array[] real y, real lowerbound)
    {
        real correction_per_value = normal_lccdf(lowerbound | 0, 1);
        target += std_normal_lpdf(y) -
                  correction_per_value * num_elements(y);
        enforceLowerBound_A_lp(y, lowerbound);
    }

    void sampleStdNormalLowerBound_2_lp(array[,] real y, real lowerbound)
    {
        real correction_per_value = normal_lccdf(lowerbound | 0, 1);
        array[2] int dimensions = dims(y);
        real correction_per_row = correction_per_value * dimensions[2];
        for (i in 1:dimensions[1]) {
            target += std_normal_lpdf(y[i]) -
                      correction_per_row;
            // ... y[i] is a one-dimensional array
        }
        enforceLowerBound_2_lp(y, lowerbound);
    }

    void sampleStdNormalLowerBound_3_lp(array[,,] real y, real lowerbound)
    {
        real correction_per_value = normal_lccdf(lowerbound | 0, 1);
        array[3] int dimensions = dims(y);
        for (i in 1:dimensions[1]) {
            for (j in 1:dimensions[2]) {
                for (k in 1:dimensions[3]) {
                    target += std_normal_lpdf(y[i, j, k]) -
                              correction_per_value;
                }
            }
        }
        enforceLowerBound_3_lp(y, lowerbound);
    }

    void sampleStdNormalLowerBound_V_lp(vector y, real lowerbound)
    {
        real correction_per_value = normal_lccdf(lowerbound | 0, 1);
        target += std_normal_lpdf(y) -
                  correction_per_value * num_elements(y);
        enforceLowerBound_V_lp(y, lowerbound);
    }

    // Sampling with upper bound

    void sampleStdNormalUpperBound_R_lp(real y, real upperbound)
    {
        real correction_per_value = normal_lcdf(upperbound | 0, 1);
        target += std_normal_lpdf(y) -
                  correction_per_value;
        enforceUpperBound_R_lp(y, upperbound);
    }

    void sampleStdNormalUpperBound_A_lp(array[] real y, real upperbound)
    {
        real correction_per_value = normal_lcdf(upperbound | 0, 1);
        target += std_normal_lpdf(y) -
                  correction_per_value * num_elements(y);
        enforceUpperBound_A_lp(y, upperbound);
    }

    void sampleStdNormalUpperBound_2_lp(array[,] real y, real upperbound)
    {
        real correction_per_value = normal_lcdf(upperbound | 0, 1);
        array[2] int dimensions = dims(y);
        real correction_per_row = correction_per_value * dimensions[2];
        for (i in 1:dimensions[1]) {
            target += std_normal_lpdf(y[i]) -
                      correction_per_row;
            // ... y[i] is a one-dimensional array
        }
        enforceUpperBound_2_lp(y, upperbound);
    }

    void sampleStdNormalUpperBound_3_lp(array[,,] real y, real upperbound)
    {
        real correction_per_value = normal_lcdf(upperbound | 0, 1);
        array[3] int dimensions = dims(y);
        for (i in 1:dimensions[1]) {
            for (j in 1:dimensions[2]) {
                for (k in 1:dimensions[3]) {
                    target += std_normal_lpdf(y[i, j, k]) -
                              correction_per_value;
                }
            }
        }
        enforceUpperBound_3_lp(y, upperbound);
    }

    void sampleStdNormalUpperBound_V_lp(vector y, real upperbound)
    {
        real correction_per_value = normal_lcdf(upperbound | 0, 1);
        target += std_normal_lpdf(y) -
                  correction_per_value * num_elements(y);
        enforceUpperBound_V_lp(y, upperbound);
    }

    // Sampling with range (lower and upper) bounds

    void sampleStdNormalRangeBound_R_lp(real y, real lowerbound, real upperbound)
    {
        real correction_per_value = log_diff_exp(
            normal_lcdf(upperbound | 0, 1),
            normal_lcdf(lowerbound | 0, 1));
        target += std_normal_lpdf(y) -
                  correction_per_value;
        enforceRangeBounds_R_lp(y, lowerbound, upperbound);
    }

    void sampleStdNormalRangeBound_A_lp(array[] real y, real lowerbound, real upperbound)
    {
        real correction_per_value = log_diff_exp(
            normal_lcdf(upperbound | 0, 1),
            normal_lcdf(lowerbound | 0, 1));
        target += std_normal_lpdf(y) -
                  correction_per_value * num_elements(y);
        enforceRangeBounds_A_lp(y, lowerbound, upperbound);
    }

    void sampleStdNormalRangeBound_2_lp(array[,] real y, real lowerbound, real upperbound)
    {
        real correction_per_value = log_diff_exp(
            normal_lcdf(upperbound | 0, 1),
            normal_lcdf(lowerbound | 0, 1));
        array[2] int dimensions = dims(y);
        real correction_per_row = correction_per_value * dimensions[2];
        for (i in 1:dimensions[1]) {
            target += std_normal_lpdf(y[i]) -
                      correction_per_row;
            // ... y[i] is a one-dimensional array
        }
        enforceRangeBounds_2_lp(y, lowerbound, upperbound);
    }

    void sampleStdNormalRangeBound_3_lp(array[,,] real y, real lowerbound, real upperbound)
    {
        real correction_per_value = log_diff_exp(
            normal_lcdf(upperbound | 0, 1),
            normal_lcdf(lowerbound | 0, 1));
        array[3] int dimensions = dims(y);
        for (i in 1:dimensions[1]) {
            for (j in 1:dimensions[2]) {
                for (k in 1:dimensions[3]) {
                    target += std_normal_lpdf(y[i, j, k]) -
                              correction_per_value;
                }
            }
        }
        enforceRangeBounds_3_lp(y, lowerbound, upperbound);
    }

    void sampleStdNormalRangeBound_V_lp(vector y, real lowerbound, real upperbound)
    {
        real correction_per_value = log_diff_exp(
            normal_lcdf(upperbound | 0, 1),
            normal_lcdf(lowerbound | 0, 1));
        target += std_normal_lpdf(y) -
                  correction_per_value * num_elements(y);
        enforceRangeBounds_V_lp(y, lowerbound, upperbound);
    }


    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    // Specials for half-standard-normal, constrained to be positive
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    void sampleStdNormalPositive_R_lp(real y)
    {
        sampleStdNormalLowerBound_R_lp(y, 0);
    }

    void sampleStdNormalPositive_A_lp(array[] real y, real lowerbound)
    {
        sampleStdNormalLowerBound_A_lp(y, 0);
    }

    void sampleStdNormalPositive_2_lp(array[,] real y, real lowerbound)
    {
        sampleStdNormalLowerBound_2_lp(y, 0);
    }

    void sampleStdNormalPositive_3_lp(array[,,] real y, real lowerbound)
    {
        sampleStdNormalLowerBound_3_lp(y, 0);
    }

    void sampleStdNormalPositive_V_lp(vector y, real lowerbound)
    {
        sampleStdNormalLowerBound_V_lp(y, 0);
    }


    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    // Cauchy distribution
    // - mu is location parameter; sigma is scale parameter
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    
    // Sampling

    void sampleCauchy_RRR_lp(real y, real mu, real sigma)
    {
        target += cauchy_lpdf(y | mu, sigma);
    }

    void sampleCauchy_ARR_lp(array[] real y, real mu, real sigma)
    {
        target += cauchy_lpdf(y | mu, sigma);
    }

    void sampleCauchy_2RR_lp(array[,] real y, real mu, real sigma)
    {
        for (i in 1:size(y)) {
            target += cauchy_lpdf(y[i] | mu, sigma);
            // ... y[i] is a one-dimensional array
        }
    }

    void sampleCauchy_3RR_lp(array[,,] real y, real mu, real sigma)
    {
        array[3] int dimensions = dims(y);
        for (i in 1:dimensions[1]) {
            for (j in 1:dimensions[2]) {
                for (k in 1:dimensions[3]) {
                    target += cauchy_lpdf(y[i, j, k] | mu, sigma);
                }
            }
        }
    }

    void sampleCauchy_VRR_lp(vector y, real mu, real sigma)
    {
        target += cauchy_lpdf(y | mu, sigma);
    }

    // Sampling with lower bound

    void sampleCauchyLowerBound_RRR_lp(real y, real mu, real sigma, real lowerbound)
    {
        real correction_per_value = cauchy_lccdf(lowerbound | mu, sigma);
        target += cauchy_lpdf(y | mu, sigma) -
                  correction_per_value;
        enforceLowerBound_R_lp(y, lowerbound);
    }

    void sampleCauchyLowerBound_ARR_lp(array[] real y, real mu, real sigma, real lowerbound)
    {
        real correction_per_value = cauchy_lccdf(lowerbound | mu, sigma);
        target += cauchy_lpdf(y | mu, sigma) -
                  correction_per_value * num_elements(y);
        enforceLowerBound_A_lp(y, lowerbound);
    }

    void sampleCauchyLowerBound_2RR_lp(array[,] real y, real mu, real sigma, real lowerbound)
    {
        real correction_per_value = cauchy_lccdf(lowerbound | mu, sigma);
        array[2] int dimensions = dims(y);
        real correction_per_row = correction_per_value * dimensions[2];
        for (i in 1:dimensions[1]) {
            target += cauchy_lpdf(y[i] | mu, sigma) -
                      correction_per_row;
            // ... y[i] is a one-dimensional array
        }
        enforceLowerBound_2_lp(y, lowerbound);
    }

    void sampleCauchyLowerBound_3RR_lp(array[,,] real y, real mu, real sigma, real lowerbound)
    {
        real correction_per_value = cauchy_lccdf(lowerbound | mu, sigma);
        array[3] int dimensions = dims(y);
        for (i in 1:dimensions[1]) {
            for (j in 1:dimensions[2]) {
                for (k in 1:dimensions[3]) {
                    target += cauchy_lpdf(y[i, j, k] | mu, sigma) -
                              correction_per_value;
                }
            }
        }
        enforceLowerBound_3_lp(y, lowerbound);
    }

    void sampleCauchyLowerBound_VRR_lp(vector y, real mu, real sigma, real lowerbound)
    {
        real correction_per_value = cauchy_lccdf(lowerbound | mu, sigma);
        target += cauchy_lpdf(y | mu, sigma) -
                  correction_per_value * num_elements(y);
        enforceLowerBound_V_lp(y, lowerbound);
    }

    // Sampling with upper bound

    void sampleCauchyUpperBound_RRR_lp(real y, real mu, real sigma, real upperbound)
    {
        real correction_per_value = cauchy_lcdf(upperbound | mu, sigma);
        target += cauchy_lpdf(y | mu, sigma) -
                  correction_per_value;
        enforceUpperBound_R_lp(y, upperbound);
    }

    void sampleCauchyUpperBound_ARR_lp(array[] real y, real mu, real sigma, real upperbound)
    {
        real correction_per_value = cauchy_lcdf(upperbound | mu, sigma);
        target += cauchy_lpdf(y | mu, sigma) -
                  correction_per_value * num_elements(y);
        enforceUpperBound_A_lp(y, upperbound);
    }

    void sampleCauchyUpperBound_2RR_lp(array[,] real y, real mu, real sigma, real upperbound)
    {
        real correction_per_value = cauchy_lcdf(upperbound | mu, sigma);
        array[2] int dimensions = dims(y);
        real correction_per_row = correction_per_value * dimensions[2];
        for (i in 1:dimensions[1]) {
            target += cauchy_lpdf(y[i] | mu, sigma) -
                      correction_per_row;
            // ... y[i] is a one-dimensional array
        }
        enforceUpperBound_2_lp(y, upperbound);
    }

    void sampleCauchyUpperBound_3RR_lp(array[,,] real y, real mu, real sigma, real upperbound)
    {
        real correction_per_value = cauchy_lcdf(upperbound | mu, sigma);
        array[3] int dimensions = dims(y);
        for (i in 1:dimensions[1]) {
            for (j in 1:dimensions[2]) {
                for (k in 1:dimensions[3]) {
                    target += cauchy_lpdf(y[i, j, k] | mu, sigma) -
                              correction_per_value;
                }
            }
        }
        enforceUpperBound_3_lp(y, upperbound);
    }

    void sampleCauchyUpperBound_VRR_lp(vector y, real mu, real sigma, real upperbound)
    {
        real correction_per_value = cauchy_lcdf(upperbound | mu, sigma);
        target += cauchy_lpdf(y | mu, sigma) -
                  correction_per_value * num_elements(y);
        enforceUpperBound_V_lp(y, upperbound);
    }

    // Sampling with range (lower and upper) bounds

    void sampleCauchyRangeBound_RRR_lp(real y, real mu, real sigma, real lowerbound, real upperbound)
    {
        real correction_per_value = log_diff_exp(
            cauchy_lcdf(upperbound | mu, sigma),
            cauchy_lcdf(lowerbound | mu, sigma));
        target += cauchy_lpdf(y | mu, sigma) -
                  correction_per_value;
        enforceRangeBounds_R_lp(y, lowerbound, upperbound);
    }

    void sampleCauchyRangeBound_ARR_lp(array[] real y, real mu, real sigma, real lowerbound, real upperbound)
    {
        real correction_per_value = log_diff_exp(
            cauchy_lcdf(upperbound | mu, sigma),
            cauchy_lcdf(lowerbound | mu, sigma));
        target += cauchy_lpdf(y | mu, sigma) -
                  correction_per_value * num_elements(y);
        enforceRangeBounds_A_lp(y, lowerbound, upperbound);
    }

    void sampleCauchyRangeBound_2RR_lp(array[,] real y, real mu, real sigma, real lowerbound, real upperbound)
    {
        real correction_per_value = log_diff_exp(
            cauchy_lcdf(upperbound | mu, sigma),
            cauchy_lcdf(lowerbound | mu, sigma));
        array[2] int dimensions = dims(y);
        real correction_per_row = correction_per_value * dimensions[2];
        for (i in 1:dimensions[1]) {
            target += cauchy_lpdf(y[i] | mu, sigma) -
                      correction_per_row;
            // ... y[i] is a one-dimensional array
        }
        enforceRangeBounds_2_lp(y, lowerbound, upperbound);
    }

    void sampleCauchyRangeBound_3RR_lp(array[,,] real y, real mu, real sigma, real lowerbound, real upperbound)
    {
        real correction_per_value = log_diff_exp(
            cauchy_lcdf(upperbound | mu, sigma),
            cauchy_lcdf(lowerbound | mu, sigma));
        array[3] int dimensions = dims(y);
        for (i in 1:dimensions[1]) {
            for (j in 1:dimensions[2]) {
                for (k in 1:dimensions[3]) {
                    target += cauchy_lpdf(y[i, j, k] | mu, sigma) -
                              correction_per_value;
                }
            }
        }
        enforceRangeBounds_3_lp(y, lowerbound, upperbound);
    }

    void sampleCauchyRangeBound_VRR_lp(vector y, real mu, real sigma, real lowerbound, real upperbound)
    {
        real correction_per_value = log_diff_exp(
            cauchy_lcdf(upperbound | mu, sigma),
            cauchy_lcdf(lowerbound | mu, sigma));
        target += cauchy_lpdf(y | mu, sigma) -
                  correction_per_value * num_elements(y);
        enforceRangeBounds_V_lp(y, lowerbound, upperbound);
    }

    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    // Beta distribution
    // - In R, alpha is called shape1, and beta is called shape2.
    // - The beta distribution is confined to the range [0, 1]. See
    //   https://en.wikipedia.org/wiki/Beta_distribution. In R, try:
    //
    //   curve(dbeta(x, shape1 = 1.2, shape2 = 1.2), -0.1, 1.1, ylab = "density")
    //
    // - Stan 2.16.0 manual p532; R ?dbeta;
    //   https://www.rdocumentation.org/packages/visualize/versions/4.3.0/topics/visualize.beta
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    
    // Sampling

    void sampleBeta_RRR_lp(real y, real alpha, real beta)
    {
        target += beta_lpdf(y | alpha, beta);
    }

    void sampleBeta_ARR_lp(array[] real y, real alpha, real beta)
    {
        target += beta_lpdf(y | alpha, beta);
    }

    void sampleBeta_2RR_lp(array[,] real y, real alpha, real beta)
    {
        for (i in 1:size(y)) {
            target += beta_lpdf(y[i] | alpha, beta);
            // ... y[i] is a one-dimensional array
        }
    }

    void sampleBeta_3RR_lp(array[,,] real y, real alpha, real beta)
    {
        array[3] int dimensions = dims(y);
        for (i in 1:dimensions[1]) {
            for (j in 1:dimensions[2]) {
                for (k in 1:dimensions[3]) {
                    target += beta_lpdf(y[i, j, k] | alpha, beta);
                }
            }
        }
    }

    void sampleBeta_VRR_lp(vector y, real alpha, real beta)
    {
        target += beta_lpdf(y | alpha, beta);
    }

    // Sampling with lower bound

    void sampleBetaLowerBound_RRR_lp(real y, real alpha, real beta, real lowerbound)
    {
        real correction_per_value = beta_lccdf(lowerbound | alpha, beta);
        target += beta_lpdf(y | alpha, beta) -
                  correction_per_value;
        enforceLowerBound_R_lp(y, lowerbound);
    }

    void sampleBetaLowerBound_ARR_lp(array[] real y, real alpha, real beta, real lowerbound)
    {
        real correction_per_value = beta_lccdf(lowerbound | alpha, beta);
        target += beta_lpdf(y | alpha, beta) -
                  correction_per_value * num_elements(y);
        enforceLowerBound_A_lp(y, lowerbound);
    }

    void sampleBetaLowerBound_2RR_lp(array[,] real y, real alpha, real beta, real lowerbound)
    {
        real correction_per_value = beta_lccdf(lowerbound | alpha, beta);
        array[2] int dimensions = dims(y);
        real correction_per_row = correction_per_value * dimensions[2];
        for (i in 1:dimensions[1]) {
            target += beta_lpdf(y[i] | alpha, beta) -
                      correction_per_row;
            // ... y[i] is a one-dimensional array
        }
        enforceLowerBound_2_lp(y, lowerbound);
    }

    void sampleBetaLowerBound_3RR_lp(array[,,] real y, real alpha, real beta, real lowerbound)
    {
        real correction_per_value = beta_lccdf(lowerbound | alpha, beta);
        array[3] int dimensions = dims(y);
        for (i in 1:dimensions[1]) {
            for (j in 1:dimensions[2]) {
                for (k in 1:dimensions[3]) {
                    target += beta_lpdf(y[i, j, k] | alpha, beta) -
                              correction_per_value;
                }
            }
        }
        enforceLowerBound_3_lp(y, lowerbound);
    }

    void sampleBetaLowerBound_VRR_lp(vector y, real alpha, real beta, real lowerbound)
    {
        real correction_per_value = beta_lccdf(lowerbound | alpha, beta);
        target += beta_lpdf(y | alpha, beta) -
                  correction_per_value * num_elements(y);
        enforceLowerBound_V_lp(y, lowerbound);
    }

    // Sampling with upper bound

    void sampleBetaUpperBound_RRR_lp(real y, real alpha, real beta, real upperbound)
    {
        real correction_per_value = beta_lcdf(upperbound | alpha, beta);
        target += beta_lpdf(y | alpha, beta) -
                  correction_per_value;
        enforceUpperBound_R_lp(y, upperbound);
    }

    void sampleBetaUpperBound_ARR_lp(array[] real y, real alpha, real beta, real upperbound)
    {
        real correction_per_value = beta_lcdf(upperbound | alpha, beta);
        target += beta_lpdf(y | alpha, beta) -
                  correction_per_value * num_elements(y);
        enforceUpperBound_A_lp(y, upperbound);
    }

    void sampleBetaUpperBound_2RR_lp(array[,] real y, real alpha, real beta, real upperbound)
    {
        real correction_per_value = beta_lcdf(upperbound | alpha, beta);
        array[2] int dimensions = dims(y);
        real correction_per_row = correction_per_value * dimensions[2];
        for (i in 1:dimensions[1]) {
            target += beta_lpdf(y[i] | alpha, beta) -
                      correction_per_row;
            // ... y[i] is a one-dimensional array
        }
        enforceUpperBound_2_lp(y, upperbound);
    }

    void sampleBetaUpperBound_3RR_lp(array[,,] real y, real alpha, real beta, real upperbound)
    {
        real correction_per_value = beta_lcdf(upperbound | alpha, beta);
        array[3] int dimensions = dims(y);
        for (i in 1:dimensions[1]) {
            for (j in 1:dimensions[2]) {
                for (k in 1:dimensions[3]) {
                    target += beta_lpdf(y[i, j, k] | alpha, beta) -
                              correction_per_value;
                }
            }
        }
        enforceUpperBound_3_lp(y, upperbound);
    }

    void sampleBetaUpperBound_VRR_lp(vector y, real alpha, real beta, real upperbound)
    {
        real correction_per_value = beta_lcdf(upperbound | alpha, beta);
        target += beta_lpdf(y | alpha, beta) -
                  correction_per_value * num_elements(y);
        enforceUpperBound_V_lp(y, upperbound);
    }

    // Sampling with range (lower and upper) bounds

    void sampleBetaRangeBound_RRR_lp(real y, real alpha, real beta, real lowerbound, real upperbound)
    {
        real correction_per_value = log_diff_exp(
            beta_lcdf(upperbound | alpha, beta),
            beta_lcdf(lowerbound | alpha, beta));
        target += beta_lpdf(y | alpha, beta) -
                  correction_per_value;
        enforceRangeBounds_R_lp(y, lowerbound, upperbound);
    }

    void sampleBetaRangeBound_ARR_lp(array[] real y, real alpha, real beta, real lowerbound, real upperbound)
    {
        real correction_per_value = log_diff_exp(
            beta_lcdf(upperbound | alpha, beta),
            beta_lcdf(lowerbound | alpha, beta));
        target += beta_lpdf(y | alpha, beta) -
                  correction_per_value * num_elements(y);
        enforceRangeBounds_A_lp(y, lowerbound, upperbound);
    }

    void sampleBetaRangeBound_2RR_lp(array[,] real y, real alpha, real beta, real lowerbound, real upperbound)
    {
        real correction_per_value = log_diff_exp(
            beta_lcdf(upperbound | alpha, beta),
            beta_lcdf(lowerbound | alpha, beta));
        array[2] int dimensions = dims(y);
        real correction_per_row = correction_per_value * dimensions[2];
        for (i in 1:dimensions[1]) {
            target += beta_lpdf(y[i] | alpha, beta) -
                      correction_per_row;
            // ... y[i] is a one-dimensional array
        }
        enforceRangeBounds_2_lp(y, lowerbound, upperbound);
    }

    void sampleBetaRangeBound_3RR_lp(array[,,] real y, real alpha, real beta, real lowerbound, real upperbound)
    {
        real correction_per_value = log_diff_exp(
            beta_lcdf(upperbound | alpha, beta),
            beta_lcdf(lowerbound | alpha, beta));
        array[3] int dimensions = dims(y);
        for (i in 1:dimensions[1]) {
            for (j in 1:dimensions[2]) {
                for (k in 1:dimensions[3]) {
                    target += beta_lpdf(y[i, j, k] | alpha, beta) -
                              correction_per_value;
                }
            }
        }
        enforceRangeBounds_3_lp(y, lowerbound, upperbound);
    }

    void sampleBetaRangeBound_VRR_lp(vector y, real alpha, real beta, real lowerbound, real upperbound)
    {
        real correction_per_value = log_diff_exp(
            beta_lcdf(upperbound | alpha, beta),
            beta_lcdf(lowerbound | alpha, beta));
        target += beta_lpdf(y | alpha, beta) -
                  correction_per_value * num_elements(y);
        enforceRangeBounds_V_lp(y, lowerbound, upperbound);
    }

    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    // Gamma distribution
    // - Stan's alpha is R's shape; Stan's beta is R's rate.
    //   (R also offers scale = 1/rate.)
    // - https://en.wikipedia.org/wiki/Gamma_distribution
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    
    // Sampling

    void sampleGamma_RRR_lp(real y, real alpha, real beta)
    {
        target += gamma_lpdf(y | alpha, beta);
    }

    void sampleGamma_ARR_lp(array[] real y, real alpha, real beta)
    {
        target += gamma_lpdf(y | alpha, beta);
    }

    void sampleGamma_2RR_lp(array[,] real y, real alpha, real beta)
    {
        for (i in 1:size(y)) {
            target += gamma_lpdf(y[i] | alpha, beta);
            // ... y[i] is a one-dimensional array
        }
    }

    void sampleGamma_3RR_lp(array[,,] real y, real alpha, real beta)
    {
        array[3] int dimensions = dims(y);
        for (i in 1:dimensions[1]) {
            for (j in 1:dimensions[2]) {
                for (k in 1:dimensions[3]) {
                    target += gamma_lpdf(y[i, j, k] | alpha, beta);
                }
            }
        }
    }

    void sampleGamma_VRR_lp(vector y, real alpha, real beta)
    {
        target += gamma_lpdf(y | alpha, beta);
    }

    // Sampling with lower bound

    void sampleGammaLowerBound_RRR_lp(real y, real alpha, real beta, real lowerbound)
    {
        real correction_per_value = gamma_lccdf(lowerbound | alpha, beta);
        target += gamma_lpdf(y | alpha, beta) -
                  correction_per_value;
        enforceLowerBound_R_lp(y, lowerbound);
    }

    void sampleGammaLowerBound_ARR_lp(array[] real y, real alpha, real beta, real lowerbound)
    {
        real correction_per_value = gamma_lccdf(lowerbound | alpha, beta);
        target += gamma_lpdf(y | alpha, beta) -
                  correction_per_value * num_elements(y);
        enforceLowerBound_A_lp(y, lowerbound);
    }

    void sampleGammaLowerBound_2RR_lp(array[,] real y, real alpha, real beta, real lowerbound)
    {
        real correction_per_value = gamma_lccdf(lowerbound | alpha, beta);
        array[2] int dimensions = dims(y);
        real correction_per_row = correction_per_value * dimensions[2];
        for (i in 1:dimensions[1]) {
            target += gamma_lpdf(y[i] | alpha, beta) -
                      correction_per_row;
            // ... y[i] is a one-dimensional array
        }
        enforceLowerBound_2_lp(y, lowerbound);
    }

    void sampleGammaLowerBound_3RR_lp(array[,,] real y, real alpha, real beta, real lowerbound)
    {
        real correction_per_value = gamma_lccdf(lowerbound | alpha, beta);
        array[3] int dimensions = dims(y);
        for (i in 1:dimensions[1]) {
            for (j in 1:dimensions[2]) {
                for (k in 1:dimensions[3]) {
                    target += gamma_lpdf(y[i, j, k] | alpha, beta) -
                              correction_per_value;
                }
            }
        }
        enforceLowerBound_3_lp(y, lowerbound);
    }

    void sampleGammaLowerBound_VRR_lp(vector y, real alpha, real beta, real lowerbound)
    {
        real correction_per_value = gamma_lccdf(lowerbound | alpha, beta);
        target += gamma_lpdf(y | alpha, beta) -
                  correction_per_value * num_elements(y);
        enforceLowerBound_V_lp(y, lowerbound);
    }

    // Sampling with upper bound

    void sampleGammaUpperBound_RRR_lp(real y, real alpha, real beta, real upperbound)
    {
        real correction_per_value = gamma_lcdf(upperbound | alpha, beta);
        target += gamma_lpdf(y | alpha, beta) -
                  correction_per_value;
        enforceUpperBound_R_lp(y, upperbound);
    }

    void sampleGammaUpperBound_ARR_lp(array[] real y, real alpha, real beta, real upperbound)
    {
        real correction_per_value = gamma_lcdf(upperbound | alpha, beta);
        target += gamma_lpdf(y | alpha, beta) -
                  correction_per_value * num_elements(y);
        enforceUpperBound_A_lp(y, upperbound);
    }

    void sampleGammaUpperBound_2RR_lp(array[,] real y, real alpha, real beta, real upperbound)
    {
        real correction_per_value = gamma_lcdf(upperbound | alpha, beta);
        array[2] int dimensions = dims(y);
        real correction_per_row = correction_per_value * dimensions[2];
        for (i in 1:dimensions[1]) {
            target += gamma_lpdf(y[i] | alpha, beta) -
                      correction_per_row;
            // ... y[i] is a one-dimensional array
        }
        enforceUpperBound_2_lp(y, upperbound);
    }

    void sampleGammaUpperBound_3RR_lp(array[,,] real y, real alpha, real beta, real upperbound)
    {
        real correction_per_value = gamma_lcdf(upperbound | alpha, beta);
        array[3] int dimensions = dims(y);
        for (i in 1:dimensions[1]) {
            for (j in 1:dimensions[2]) {
                for (k in 1:dimensions[3]) {
                    target += gamma_lpdf(y[i, j, k] | alpha, beta) -
                              correction_per_value;
                }
            }
        }
        enforceUpperBound_3_lp(y, upperbound);
    }

    void sampleGammaUpperBound_VRR_lp(vector y, real alpha, real beta, real upperbound)
    {
        real correction_per_value = gamma_lcdf(upperbound | alpha, beta);
        target += gamma_lpdf(y | alpha, beta) -
                  correction_per_value * num_elements(y);
        enforceUpperBound_V_lp(y, upperbound);
    }

    // Sampling with range (lower and upper) bounds

    void sampleGammaRangeBound_RRR_lp(real y, real alpha, real beta, real lowerbound, real upperbound)
    {
        real correction_per_value = log_diff_exp(
            gamma_lcdf(upperbound | alpha, beta),
            gamma_lcdf(lowerbound | alpha, beta));
        target += gamma_lpdf(y | alpha, beta) -
                  correction_per_value;
        enforceRangeBounds_R_lp(y, lowerbound, upperbound);
    }

    void sampleGammaRangeBound_ARR_lp(array[] real y, real alpha, real beta, real lowerbound, real upperbound)
    {
        real correction_per_value = log_diff_exp(
            gamma_lcdf(upperbound | alpha, beta),
            gamma_lcdf(lowerbound | alpha, beta));
        target += gamma_lpdf(y | alpha, beta) -
                  correction_per_value * num_elements(y);
        enforceRangeBounds_A_lp(y, lowerbound, upperbound);
    }

    void sampleGammaRangeBound_2RR_lp(array[,] real y, real alpha, real beta, real lowerbound, real upperbound)
    {
        real correction_per_value = log_diff_exp(
            gamma_lcdf(upperbound | alpha, beta),
            gamma_lcdf(lowerbound | alpha, beta));
        array[2] int dimensions = dims(y);
        real correction_per_row = correction_per_value * dimensions[2];
        for (i in 1:dimensions[1]) {
            target += gamma_lpdf(y[i] | alpha, beta) -
                      correction_per_row;
            // ... y[i] is a one-dimensional array
        }
        enforceRangeBounds_2_lp(y, lowerbound, upperbound);
    }

    void sampleGammaRangeBound_3RR_lp(array[,,] real y, real alpha, real beta, real lowerbound, real upperbound)
    {
        real correction_per_value = log_diff_exp(
            gamma_lcdf(upperbound | alpha, beta),
            gamma_lcdf(lowerbound | alpha, beta));
        array[3] int dimensions = dims(y);
        for (i in 1:dimensions[1]) {
            for (j in 1:dimensions[2]) {
                for (k in 1:dimensions[3]) {
                    target += gamma_lpdf(y[i, j, k] | alpha, beta) -
                              correction_per_value;
                }
            }
        }
        enforceRangeBounds_3_lp(y, lowerbound, upperbound);
    }

    void sampleGammaRangeBound_VRR_lp(vector y, real alpha, real beta, real lowerbound, real upperbound)
    {
        real correction_per_value = log_diff_exp(
            gamma_lcdf(upperbound | alpha, beta),
            gamma_lcdf(lowerbound | alpha, beta));
        target += gamma_lpdf(y | alpha, beta) -
                  correction_per_value * num_elements(y);
        enforceRangeBounds_V_lp(y, lowerbound, upperbound);
    }

    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    // Uniform distribution
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    // Always constrained with both a lower and an upper bound.
    // Simple; no extra work for the bridge sampler.
    
    // Sampling

    void sampleUniform_RRR_lp(real y, real lowerbound, real upperbound)
    {
        target += uniform_lpdf(y | lowerbound, upperbound);
    }
    
    void sampleUniform_ARR_lp(array[] real y, real lowerbound, real upperbound)
    {
        target += uniform_lpdf(y | lowerbound, upperbound);
    }
    
    void sampleUniform_ARA_lp(array[] real y, real lowerbound, array[] real upperbound)
    {
        target += uniform_lpdf(y | lowerbound, upperbound);
    }
    
    void sampleUniform_ARV_lp(array[] real y, real lowerbound, vector upperbound)
    {
        target += uniform_lpdf(y | lowerbound, upperbound);
    }
    
    void sampleUniform_AAR_lp(array[] real y, array[] real lowerbound, real upperbound)
    {
        target += uniform_lpdf(y | lowerbound, upperbound);
    }
    
    void sampleUniform_AAA_lp(array[] real y, array[] real lowerbound, array[] real upperbound)
    {
        target += uniform_lpdf(y | lowerbound, upperbound);
    }
    
    void sampleUniform_AAV_lp(array[] real y, array[] real lowerbound, vector upperbound)
    {
        target += uniform_lpdf(y | lowerbound, upperbound);
    }
    
    void sampleUniform_AVR_lp(array[] real y, vector lowerbound, real upperbound)
    {
        target += uniform_lpdf(y | lowerbound, upperbound);
    }
    
    void sampleUniform_AVA_lp(array[] real y, vector lowerbound, array[] real upperbound)
    {
        target += uniform_lpdf(y | lowerbound, upperbound);
    }
    
    void sampleUniform_AVV_lp(array[] real y, vector lowerbound, vector upperbound)
    {
        target += uniform_lpdf(y | lowerbound, upperbound);
    }
    
    void sampleUniform_2RR_lp(array[,] real y, real lowerbound, real upperbound)
    {
        for (i in 1:size(y)) {
            target += uniform_lpdf(y[i] | lowerbound, upperbound);
            // ... y[i] is a one-dimensional array
        }
    }
    
    void sampleUniform_3RR_lp(array[,,] real y, real lowerbound, real upperbound)
    {
        array[3] int dimensions = dims(y);
        for (i in 1:dimensions[1]) {
            for (j in 1:dimensions[2]) {
                for (k in 1:dimensions[3]) {
                    target += uniform_lpdf(y[i, j, k] | lowerbound, upperbound);
                }
            }
        }
    }
    
    void sampleUniform_VRR_lp(vector y, real lowerbound, real upperbound)
    {
        target += uniform_lpdf(y | lowerbound, upperbound);
    }
    
    void sampleUniform_VRA_lp(vector y, real lowerbound, array[] real upperbound)
    {
        target += uniform_lpdf(y | lowerbound, upperbound);
    }
    
    void sampleUniform_VRV_lp(vector y, real lowerbound, vector upperbound)
    {
        target += uniform_lpdf(y | lowerbound, upperbound);
    }
    
    void sampleUniform_VAR_lp(vector y, array[] real lowerbound, real upperbound)
    {
        target += uniform_lpdf(y | lowerbound, upperbound);
    }
    
    void sampleUniform_VAA_lp(vector y, array[] real lowerbound, array[] real upperbound)
    {
        target += uniform_lpdf(y | lowerbound, upperbound);
    }
    
    void sampleUniform_VAV_lp(vector y, array[] real lowerbound, vector upperbound)
    {
        target += uniform_lpdf(y | lowerbound, upperbound);
    }
    
    void sampleUniform_VVR_lp(vector y, vector lowerbound, real upperbound)
    {
        target += uniform_lpdf(y | lowerbound, upperbound);
    }
    
    void sampleUniform_VVA_lp(vector y, vector lowerbound, array[] real upperbound)
    {
        target += uniform_lpdf(y | lowerbound, upperbound);
    }
    
    void sampleUniform_VVV_lp(vector y, vector lowerbound, vector upperbound)
    {
        target += uniform_lpdf(y | lowerbound, upperbound);
    }
    
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    // Bernoulli distribution
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    // y is in {0, 1} and theta is a probability in the range [0, 1].

    void sampleBernoulli_IR_lp(int y, real theta)
    {
        target += bernoulli_lpmf(y | theta);
    }
    void sampleBernoulli_AR_lp(array[] int y, real theta)
    {
        target += bernoulli_lpmf(y | theta);
    }
    void sampleBernoulli_AA_lp(array[] int y, array[] real theta)
    {
        target += bernoulli_lpmf(y | theta);
    }
    void sampleBernoulli_AV_lp(array[] int y, vector theta)
    {
        target += bernoulli_lpmf(y | theta);
    }

    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    // Bernoulli logit distribution
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    // y is in {0, 1} and alpha is a logit (log odds) in the range [-inf, +inf].

    void sampleBernoulliLogit_IR_lp(int y, real alpha)
    {
        target += bernoulli_logit_lpmf(y | alpha);
    }
    void sampleBernoulliLogit_AR_lp(array[] int y, real alpha)
    {
        target += bernoulli_logit_lpmf(y | alpha);
    }
    void sampleBernoulliLogit_AA_lp(array[] int y, array[] real alpha)
    {
        target += bernoulli_logit_lpmf(y | alpha);
    }
    void sampleBernoulliLogit_AV_lp(array[] int y, vector alpha)
    {
        target += bernoulli_logit_lpmf(y | alpha);
    }


    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    // Categorical distribution
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    // y is in {1, ..., K} and theta is a K-simplex (of the probabilities of
    // each of the K outcomes). An N-simplex is a vector (etc.) of non-negative
    // numbers that sum to 1.
    // Note that theta must be vector, not "reals".
    // The logit version is such that categorical_logit_lpmf(beta) is the same
    // as categorical_logit(softmax(beta)), i.e. theta = softmax(beta).

    void sampleCategorical_IV_lp(int y, vector theta)
    {
        target += categorical_lpmf(y | theta);
    }
    void sampleCategorical_AV_lp(array[] int y, vector theta)
    {
        target += categorical_lpmf(y | theta);
    }

    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    // Categorical logit distribution
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    // y is in {1, ..., K} and beta is a K-vector of (log odds) in the range
    // [-inf, +inf]; theta = softmax(beta) as above.

    void sampleCategoricalLogit_IV_lp(int y, vector beta)
    {
        target += categorical_logit_lpmf(y | beta);
    }
    void sampleCategoricalLogit_AV_lp(array[] int y, vector beta)
    {
        target += categorical_logit_lpmf(y | beta);
    }


    // ------------------------------------------------------------------------
    // LOG PROBABILITY FUNCTIONS FOR BRIDGE SAMPLING WITH NON-CENTERED
    // REPARAMETERIZATION
    // ------------------------------------------------------------------------

    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    // Normal distribution, reparameterized to the unit normal distribution
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    // Compare Stan (2017) manual p299, but we use a bridgesampling version.

    real reparameterizedNormalBoundary(real boundary, real mu, real sigma)
    {
        // boundary: in real-world N(mu, sigma) space
        // return value: equivalent in the reparameterized N(0, 1) space
        return (boundary - mu) / sigma;
    }
    
    // Plain

    real getReparameterizedNormal_RRR_lp(real y_unit_normal, real mu, real sigma)
    {
        real result;
        sampleNormal_RRR_lp(y_unit_normal, 0, 1);
        result = mu + sigma * y_unit_normal;
        return result;
    }
    
    array[] real getReparameterizedNormal_ARR_lp(array[] real y_unit_normal, real mu, real sigma)
    {
        int length = num_elements(y_unit_normal);
        array[length] real result;
        for (i in 1:length) {
            sampleNormal_RRR_lp(y_unit_normal[i], 0, 1);
            result[i] = mu + sigma * y_unit_normal[i];
        }
        return result;
    }
    
    array[] real getReparameterizedNormal_ARA_lp(array[] real y_unit_normal, real mu, array[] real sigma)
    {
        int length = num_elements(y_unit_normal);
        array[length] real result;
        if (num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            sampleNormal_RRR_lp(y_unit_normal[i], 0, 1);
            result[i] = mu + sigma[i] * y_unit_normal[i];
        }
        return result;
    }
    
    array[] real getReparameterizedNormal_ARV_lp(array[] real y_unit_normal, real mu, vector sigma)
    {
        int length = num_elements(y_unit_normal);
        array[length] real result;
        if (num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            sampleNormal_RRR_lp(y_unit_normal[i], 0, 1);
            result[i] = mu + sigma[i] * y_unit_normal[i];
        }
        return result;
    }
    
    array[] real getReparameterizedNormal_AAR_lp(array[] real y_unit_normal, array[] real mu, real sigma)
    {
        int length = num_elements(y_unit_normal);
        array[length] real result;
        if (num_elements(mu) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            sampleNormal_RRR_lp(y_unit_normal[i], 0, 1);
            result[i] = mu[i] + sigma * y_unit_normal[i];
        }
        return result;
    }
    
    array[] real getReparameterizedNormal_AAA_lp(array[] real y_unit_normal, array[] real mu, array[] real sigma)
    {
        int length = num_elements(y_unit_normal);
        array[length] real result;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            sampleNormal_RRR_lp(y_unit_normal[i], 0, 1);
            result[i] = mu[i] + sigma[i] * y_unit_normal[i];
        }
        return result;
    }
    
    array[] real getReparameterizedNormal_AAV_lp(array[] real y_unit_normal, array[] real mu, vector sigma)
    {
        int length = num_elements(y_unit_normal);
        array[length] real result;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            sampleNormal_RRR_lp(y_unit_normal[i], 0, 1);
            result[i] = mu[i] + sigma[i] * y_unit_normal[i];
        }
        return result;
    }
    
    array[] real getReparameterizedNormal_AVR_lp(array[] real y_unit_normal, vector mu, real sigma)
    {
        int length = num_elements(y_unit_normal);
        array[length] real result;
        if (num_elements(mu) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            sampleNormal_RRR_lp(y_unit_normal[i], 0, 1);
            result[i] = mu[i] + sigma * y_unit_normal[i];
        }
        return result;
    }
    
    array[] real getReparameterizedNormal_AVA_lp(array[] real y_unit_normal, vector mu, array[] real sigma)
    {
        int length = num_elements(y_unit_normal);
        array[length] real result;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            sampleNormal_RRR_lp(y_unit_normal[i], 0, 1);
            result[i] = mu[i] + sigma[i] * y_unit_normal[i];
        }
        return result;
    }
    
    array[] real getReparameterizedNormal_AVV_lp(array[] real y_unit_normal, vector mu, vector sigma)
    {
        int length = num_elements(y_unit_normal);
        array[length] real result;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            sampleNormal_RRR_lp(y_unit_normal[i], 0, 1);
            result[i] = mu[i] + sigma[i] * y_unit_normal[i];
        }
        return result;
    }
    
    array[,] real getReparameterizedNormal_2RR_lp(array[,] real y_unit_normal, real mu, real sigma)
    {
        array[2] int dimensions = dims(y_unit_normal);
        array[dimensions[1], dimensions[2]] real result;
        for (i in 1:dimensions[1]) {
            for (j in 1:dimensions[2]) {
                sampleNormal_RRR_lp(y_unit_normal[i, j], 0, 1);
                result[i, j] = mu + sigma * y_unit_normal[i, j];
            }
        }
        return result;
    }
    
    array[,,] real getReparameterizedNormal_3RR_lp(array[,,] real y_unit_normal, real mu, real sigma)
    {
        array[3] int dimensions = dims(y_unit_normal);
        array[dimensions[1], dimensions[2], dimensions[3]] real result;
        for (i in 1:dimensions[1]) {
            for (j in 1:dimensions[2]) {
                for (k in 1:dimensions[3]) {
                    sampleNormal_RRR_lp(y_unit_normal[i, j, k], 0, 1);
                    result[i, j, k] = mu + sigma * y_unit_normal[i, j, k];
                }
            }
        }
        return result;
    }
    
    vector getReparameterizedNormal_VRR_lp(vector y_unit_normal, real mu, real sigma)
    {
        int length = num_elements(y_unit_normal);
        vector[length] result;
        for (i in 1:length) {
            sampleNormal_RRR_lp(y_unit_normal[i], 0, 1);
            result[i] = mu + sigma * y_unit_normal[i];
        }
        return result;
    }
    
    vector getReparameterizedNormal_VRA_lp(vector y_unit_normal, real mu, array[] real sigma)
    {
        int length = num_elements(y_unit_normal);
        vector[length] result;
        if (num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            sampleNormal_RRR_lp(y_unit_normal[i], 0, 1);
            result[i] = mu + sigma[i] * y_unit_normal[i];
        }
        return result;
    }
    
    vector getReparameterizedNormal_VRV_lp(vector y_unit_normal, real mu, vector sigma)
    {
        int length = num_elements(y_unit_normal);
        vector[length] result;
        if (num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            sampleNormal_RRR_lp(y_unit_normal[i], 0, 1);
            result[i] = mu + sigma[i] * y_unit_normal[i];
        }
        return result;
    }
    
    vector getReparameterizedNormal_VAR_lp(vector y_unit_normal, array[] real mu, real sigma)
    {
        int length = num_elements(y_unit_normal);
        vector[length] result;
        if (num_elements(mu) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            sampleNormal_RRR_lp(y_unit_normal[i], 0, 1);
            result[i] = mu[i] + sigma * y_unit_normal[i];
        }
        return result;
    }
    
    vector getReparameterizedNormal_VAA_lp(vector y_unit_normal, array[] real mu, array[] real sigma)
    {
        int length = num_elements(y_unit_normal);
        vector[length] result;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            sampleNormal_RRR_lp(y_unit_normal[i], 0, 1);
            result[i] = mu[i] + sigma[i] * y_unit_normal[i];
        }
        return result;
    }
    
    vector getReparameterizedNormal_VAV_lp(vector y_unit_normal, array[] real mu, vector sigma)
    {
        int length = num_elements(y_unit_normal);
        vector[length] result;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            sampleNormal_RRR_lp(y_unit_normal[i], 0, 1);
            result[i] = mu[i] + sigma[i] * y_unit_normal[i];
        }
        return result;
    }
    
    vector getReparameterizedNormal_VVR_lp(vector y_unit_normal, vector mu, real sigma)
    {
        int length = num_elements(y_unit_normal);
        vector[length] result;
        if (num_elements(mu) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            sampleNormal_RRR_lp(y_unit_normal[i], 0, 1);
            result[i] = mu[i] + sigma * y_unit_normal[i];
        }
        return result;
    }
    
    vector getReparameterizedNormal_VVA_lp(vector y_unit_normal, vector mu, array[] real sigma)
    {
        int length = num_elements(y_unit_normal);
        vector[length] result;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            sampleNormal_RRR_lp(y_unit_normal[i], 0, 1);
            result[i] = mu[i] + sigma[i] * y_unit_normal[i];
        }
        return result;
    }
    
    vector getReparameterizedNormal_VVV_lp(vector y_unit_normal, vector mu, vector sigma)
    {
        int length = num_elements(y_unit_normal);
        vector[length] result;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            sampleNormal_RRR_lp(y_unit_normal[i], 0, 1);
            result[i] = mu[i] + sigma[i] * y_unit_normal[i];
        }
        return result;
    }
    
    // With lower bound

    real getReparameterizedNormalLowerBound_RRR_lp(real y_unit_normal, real mu, real sigma, real lowerbound)
    {
        real result;
        real lower_transformed;
        lower_transformed = (lowerbound - mu) / sigma;
        sampleNormalLowerBound_RRR_lp(y_unit_normal, 0, 1, lower_transformed);
        result = mu + sigma * y_unit_normal;
        return result;
    }
    
    array[] real getReparameterizedNormalLowerBound_ARR_lp(array[] real y_unit_normal, real mu, real sigma, real lowerbound)
    {
        int length = num_elements(y_unit_normal);
        array[length] real result;
        real lower_transformed;
        for (i in 1:length) {
            lower_transformed = (lowerbound - mu) / sigma;
            sampleNormalLowerBound_RRR_lp(y_unit_normal[i], 0, 1, lower_transformed);
            result[i] = mu + sigma * y_unit_normal[i];
        }
        return result;
    }
    
    array[] real getReparameterizedNormalLowerBound_ARA_lp(array[] real y_unit_normal, real mu, array[] real sigma, real lowerbound)
    {
        int length = num_elements(y_unit_normal);
        array[length] real result;
        real lower_transformed;
        if (num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = (lowerbound - mu) / sigma[i];
            sampleNormalLowerBound_RRR_lp(y_unit_normal[i], 0, 1, lower_transformed);
            result[i] = mu + sigma[i] * y_unit_normal[i];
        }
        return result;
    }
    
    array[] real getReparameterizedNormalLowerBound_ARV_lp(array[] real y_unit_normal, real mu, vector sigma, real lowerbound)
    {
        int length = num_elements(y_unit_normal);
        array[length] real result;
        real lower_transformed;
        if (num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = (lowerbound - mu) / sigma[i];
            sampleNormalLowerBound_RRR_lp(y_unit_normal[i], 0, 1, lower_transformed);
            result[i] = mu + sigma[i] * y_unit_normal[i];
        }
        return result;
    }
    
    array[] real getReparameterizedNormalLowerBound_AAR_lp(array[] real y_unit_normal, array[] real mu, real sigma, real lowerbound)
    {
        int length = num_elements(y_unit_normal);
        array[length] real result;
        real lower_transformed;
        if (num_elements(mu) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = (lowerbound - mu[i]) / sigma;
            sampleNormalLowerBound_RRR_lp(y_unit_normal[i], 0, 1, lower_transformed);
            result[i] = mu[i] + sigma * y_unit_normal[i];
        }
        return result;
    }
    
    array[] real getReparameterizedNormalLowerBound_AAA_lp(array[] real y_unit_normal, array[] real mu, array[] real sigma, real lowerbound)
    {
        int length = num_elements(y_unit_normal);
        array[length] real result;
        real lower_transformed;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = (lowerbound - mu[i]) / sigma[i];
            sampleNormalLowerBound_RRR_lp(y_unit_normal[i], 0, 1, lower_transformed);
            result[i] = mu[i] + sigma[i] * y_unit_normal[i];
        }
        return result;
    }
    
    array[] real getReparameterizedNormalLowerBound_AAV_lp(array[] real y_unit_normal, array[] real mu, vector sigma, real lowerbound)
    {
        int length = num_elements(y_unit_normal);
        array[length] real result;
        real lower_transformed;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = (lowerbound - mu[i]) / sigma[i];
            sampleNormalLowerBound_RRR_lp(y_unit_normal[i], 0, 1, lower_transformed);
            result[i] = mu[i] + sigma[i] * y_unit_normal[i];
        }
        return result;
    }
    
    array[] real getReparameterizedNormalLowerBound_AVR_lp(array[] real y_unit_normal, vector mu, real sigma, real lowerbound)
    {
        int length = num_elements(y_unit_normal);
        array[length] real result;
        real lower_transformed;
        if (num_elements(mu) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = (lowerbound - mu[i]) / sigma;
            sampleNormalLowerBound_RRR_lp(y_unit_normal[i], 0, 1, lower_transformed);
            result[i] = mu[i] + sigma * y_unit_normal[i];
        }
        return result;
    }
    
    array[] real getReparameterizedNormalLowerBound_AVA_lp(array[] real y_unit_normal, vector mu, array[] real sigma, real lowerbound)
    {
        int length = num_elements(y_unit_normal);
        array[length] real result;
        real lower_transformed;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = (lowerbound - mu[i]) / sigma[i];
            sampleNormalLowerBound_RRR_lp(y_unit_normal[i], 0, 1, lower_transformed);
            result[i] = mu[i] + sigma[i] * y_unit_normal[i];
        }
        return result;
    }
    
    array[] real getReparameterizedNormalLowerBound_AVV_lp(array[] real y_unit_normal, vector mu, vector sigma, real lowerbound)
    {
        int length = num_elements(y_unit_normal);
        array[length] real result;
        real lower_transformed;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = (lowerbound - mu[i]) / sigma[i];
            sampleNormalLowerBound_RRR_lp(y_unit_normal[i], 0, 1, lower_transformed);
            result[i] = mu[i] + sigma[i] * y_unit_normal[i];
        }
        return result;
    }
    
    array[,] real getReparameterizedNormalLowerBound_2RR_lp(array[,] real y_unit_normal, real mu, real sigma, real lowerbound)
    {
        array[2] int dimensions = dims(y_unit_normal);
        array[dimensions[1], dimensions[2]] real result;
        real lower_transformed;
        lower_transformed = (lowerbound - mu) / sigma;
        for (i in 1:dimensions[1]) {
            for (j in 1:dimensions[2]) {
                sampleNormalLowerBound_RRR_lp(y_unit_normal[i, j], 0, 1, lower_transformed);
                result[i, j] = mu + sigma * y_unit_normal[i, j];
            }
        }
        return result;
    }
    
    array[,,] real getReparameterizedNormalLowerBound_3RR_lp(array[,,] real y_unit_normal, real mu, real sigma, real lowerbound)
    {
        array[3] int dimensions = dims(y_unit_normal);
        array[dimensions[1], dimensions[2], dimensions[3]] real result;
        real lower_transformed;
        lower_transformed = (lowerbound - mu) / sigma;
        for (i in 1:dimensions[1]) {
            for (j in 1:dimensions[2]) {
                for (k in 1:dimensions[3]) {
                    sampleNormalLowerBound_RRR_lp(y_unit_normal[i, j, k], 0, 1, lower_transformed);
                    result[i, j, k] = mu + sigma * y_unit_normal[i, j, k];
                }
            }
        }
        return result;
    }
    
    vector getReparameterizedNormalLowerBound_VRR_lp(vector y_unit_normal, real mu, real sigma, real lowerbound)
    {
        int length = num_elements(y_unit_normal);
        vector[length] result;
        real lower_transformed;
        for (i in 1:length) {
            lower_transformed = (lowerbound - mu) / sigma;
            sampleNormalLowerBound_RRR_lp(y_unit_normal[i], 0, 1, lower_transformed);
            result[i] = mu + sigma * y_unit_normal[i];
        }
        return result;
    }
    
    vector getReparameterizedNormalLowerBound_VRA_lp(vector y_unit_normal, real mu, array[] real sigma, real lowerbound)
    {
        int length = num_elements(y_unit_normal);
        vector[length] result;
        real lower_transformed;
        if (num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = (lowerbound - mu) / sigma[i];
            sampleNormalLowerBound_RRR_lp(y_unit_normal[i], 0, 1, lower_transformed);
            result[i] = mu + sigma[i] * y_unit_normal[i];
        }
        return result;
    }
    
    vector getReparameterizedNormalLowerBound_VRV_lp(vector y_unit_normal, real mu, vector sigma, real lowerbound)
    {
        int length = num_elements(y_unit_normal);
        vector[length] result;
        real lower_transformed;
        if (num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = (lowerbound - mu) / sigma[i];
            sampleNormalLowerBound_RRR_lp(y_unit_normal[i], 0, 1, lower_transformed);
            result[i] = mu + sigma[i] * y_unit_normal[i];
        }
        return result;
    }
    
    vector getReparameterizedNormalLowerBound_VAR_lp(vector y_unit_normal, array[] real mu, real sigma, real lowerbound)
    {
        int length = num_elements(y_unit_normal);
        vector[length] result;
        real lower_transformed;
        if (num_elements(mu) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = (lowerbound - mu[i]) / sigma;
            sampleNormalLowerBound_RRR_lp(y_unit_normal[i], 0, 1, lower_transformed);
            result[i] = mu[i] + sigma * y_unit_normal[i];
        }
        return result;
    }
    
    vector getReparameterizedNormalLowerBound_VAA_lp(vector y_unit_normal, array[] real mu, array[] real sigma, real lowerbound)
    {
        int length = num_elements(y_unit_normal);
        vector[length] result;
        real lower_transformed;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = (lowerbound - mu[i]) / sigma[i];
            sampleNormalLowerBound_RRR_lp(y_unit_normal[i], 0, 1, lower_transformed);
            result[i] = mu[i] + sigma[i] * y_unit_normal[i];
        }
        return result;
    }
    
    vector getReparameterizedNormalLowerBound_VAV_lp(vector y_unit_normal, array[] real mu, vector sigma, real lowerbound)
    {
        int length = num_elements(y_unit_normal);
        vector[length] result;
        real lower_transformed;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = (lowerbound - mu[i]) / sigma[i];
            sampleNormalLowerBound_RRR_lp(y_unit_normal[i], 0, 1, lower_transformed);
            result[i] = mu[i] + sigma[i] * y_unit_normal[i];
        }
        return result;
    }
    
    vector getReparameterizedNormalLowerBound_VVR_lp(vector y_unit_normal, vector mu, real sigma, real lowerbound)
    {
        int length = num_elements(y_unit_normal);
        vector[length] result;
        real lower_transformed;
        if (num_elements(mu) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = (lowerbound - mu[i]) / sigma;
            sampleNormalLowerBound_RRR_lp(y_unit_normal[i], 0, 1, lower_transformed);
            result[i] = mu[i] + sigma * y_unit_normal[i];
        }
        return result;
    }
    
    vector getReparameterizedNormalLowerBound_VVA_lp(vector y_unit_normal, vector mu, array[] real sigma, real lowerbound)
    {
        int length = num_elements(y_unit_normal);
        vector[length] result;
        real lower_transformed;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = (lowerbound - mu[i]) / sigma[i];
            sampleNormalLowerBound_RRR_lp(y_unit_normal[i], 0, 1, lower_transformed);
            result[i] = mu[i] + sigma[i] * y_unit_normal[i];
        }
        return result;
    }
    
    vector getReparameterizedNormalLowerBound_VVV_lp(vector y_unit_normal, vector mu, vector sigma, real lowerbound)
    {
        int length = num_elements(y_unit_normal);
        vector[length] result;
        real lower_transformed;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = (lowerbound - mu[i]) / sigma[i];
            sampleNormalLowerBound_RRR_lp(y_unit_normal[i], 0, 1, lower_transformed);
            result[i] = mu[i] + sigma[i] * y_unit_normal[i];
        }
        return result;
    }
    
    // With upper bound

    real getReparameterizedNormalUpperBound_RRR_lp(real y_unit_normal, real mu, real sigma, real upperbound)
    {
        real result;
        real upper_transformed;
        upper_transformed = (upperbound - mu) / sigma;
        sampleNormalUpperBound_RRR_lp(y_unit_normal, 0, 1, upper_transformed);
        result = mu + sigma * y_unit_normal;
        return result;
    }
    
    array[] real getReparameterizedNormalUpperBound_ARR_lp(array[] real y_unit_normal, real mu, real sigma, real upperbound)
    {
        int length = num_elements(y_unit_normal);
        array[length] real result;
        real upper_transformed;
        for (i in 1:length) {
            upper_transformed = (upperbound - mu) / sigma;
            sampleNormalUpperBound_RRR_lp(y_unit_normal[i], 0, 1, upper_transformed);
            result[i] = mu + sigma * y_unit_normal[i];
        }
        return result;
    }
    
    array[] real getReparameterizedNormalUpperBound_ARA_lp(array[] real y_unit_normal, real mu, array[] real sigma, real upperbound)
    {
        int length = num_elements(y_unit_normal);
        array[length] real result;
        real upper_transformed;
        if (num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            upper_transformed = (upperbound - mu) / sigma[i];
            sampleNormalUpperBound_RRR_lp(y_unit_normal[i], 0, 1, upper_transformed);
            result[i] = mu + sigma[i] * y_unit_normal[i];
        }
        return result;
    }
    
    array[] real getReparameterizedNormalUpperBound_ARV_lp(array[] real y_unit_normal, real mu, vector sigma, real upperbound)
    {
        int length = num_elements(y_unit_normal);
        array[length] real result;
        real upper_transformed;
        if (num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            upper_transformed = (upperbound - mu) / sigma[i];
            sampleNormalUpperBound_RRR_lp(y_unit_normal[i], 0, 1, upper_transformed);
            result[i] = mu + sigma[i] * y_unit_normal[i];
        }
        return result;
    }
    
    array[] real getReparameterizedNormalUpperBound_AAR_lp(array[] real y_unit_normal, array[] real mu, real sigma, real upperbound)
    {
        int length = num_elements(y_unit_normal);
        array[length] real result;
        real upper_transformed;
        if (num_elements(mu) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            upper_transformed = (upperbound - mu[i]) / sigma;
            sampleNormalUpperBound_RRR_lp(y_unit_normal[i], 0, 1, upper_transformed);
            result[i] = mu[i] + sigma * y_unit_normal[i];
        }
        return result;
    }
    
    array[] real getReparameterizedNormalUpperBound_AAA_lp(array[] real y_unit_normal, array[] real mu, array[] real sigma, real upperbound)
    {
        int length = num_elements(y_unit_normal);
        array[length] real result;
        real upper_transformed;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            upper_transformed = (upperbound - mu[i]) / sigma[i];
            sampleNormalUpperBound_RRR_lp(y_unit_normal[i], 0, 1, upper_transformed);
            result[i] = mu[i] + sigma[i] * y_unit_normal[i];
        }
        return result;
    }
    
    array[] real getReparameterizedNormalUpperBound_AAV_lp(array[] real y_unit_normal, array[] real mu, vector sigma, real upperbound)
    {
        int length = num_elements(y_unit_normal);
        array[length] real result;
        real upper_transformed;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            upper_transformed = (upperbound - mu[i]) / sigma[i];
            sampleNormalUpperBound_RRR_lp(y_unit_normal[i], 0, 1, upper_transformed);
            result[i] = mu[i] + sigma[i] * y_unit_normal[i];
        }
        return result;
    }
    
    array[] real getReparameterizedNormalUpperBound_AVR_lp(array[] real y_unit_normal, vector mu, real sigma, real upperbound)
    {
        int length = num_elements(y_unit_normal);
        array[length] real result;
        real upper_transformed;
        if (num_elements(mu) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            upper_transformed = (upperbound - mu[i]) / sigma;
            sampleNormalUpperBound_RRR_lp(y_unit_normal[i], 0, 1, upper_transformed);
            result[i] = mu[i] + sigma * y_unit_normal[i];
        }
        return result;
    }
    
    array[] real getReparameterizedNormalUpperBound_AVA_lp(array[] real y_unit_normal, vector mu, array[] real sigma, real upperbound)
    {
        int length = num_elements(y_unit_normal);
        array[length] real result;
        real upper_transformed;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            upper_transformed = (upperbound - mu[i]) / sigma[i];
            sampleNormalUpperBound_RRR_lp(y_unit_normal[i], 0, 1, upper_transformed);
            result[i] = mu[i] + sigma[i] * y_unit_normal[i];
        }
        return result;
    }
    
    array[] real getReparameterizedNormalUpperBound_AVV_lp(array[] real y_unit_normal, vector mu, vector sigma, real upperbound)
    {
        int length = num_elements(y_unit_normal);
        array[length] real result;
        real upper_transformed;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            upper_transformed = (upperbound - mu[i]) / sigma[i];
            sampleNormalUpperBound_RRR_lp(y_unit_normal[i], 0, 1, upper_transformed);
            result[i] = mu[i] + sigma[i] * y_unit_normal[i];
        }
        return result;
    }
    
    array[,] real getReparameterizedNormalUpperBound_2RR_lp(array[,] real y_unit_normal, real mu, real sigma, real upperbound)
    {
        array[2] int dimensions = dims(y_unit_normal);
        array[dimensions[1], dimensions[2]] real result;
        real upper_transformed;
        upper_transformed = (upperbound - mu) / sigma;
        for (i in 1:dimensions[1]) {
            for (j in 1:dimensions[2]) {
                sampleNormalUpperBound_RRR_lp(y_unit_normal[i, j], 0, 1, upper_transformed);
                result[i, j] = mu + sigma * y_unit_normal[i, j];
            }
        }
        return result;
    }
    
    array[,,] real getReparameterizedNormalUpperBound_3RR_lp(array[,,] real y_unit_normal, real mu, real sigma, real upperbound)
    {
        array[3] int dimensions = dims(y_unit_normal);
        array[dimensions[1], dimensions[2], dimensions[3]] real result;
        real upper_transformed;
        upper_transformed = (upperbound - mu) / sigma;
        for (i in 1:dimensions[1]) {
            for (j in 1:dimensions[2]) {
                for (k in 1:dimensions[3]) {
                    sampleNormalUpperBound_RRR_lp(y_unit_normal[i, j, k], 0, 1, upper_transformed);
                    result[i, j, k] = mu + sigma * y_unit_normal[i, j, k];
                }
            }
        }
        return result;
    }
    
    vector getReparameterizedNormalUpperBound_VRR_lp(vector y_unit_normal, real mu, real sigma, real upperbound)
    {
        int length = num_elements(y_unit_normal);
        vector[length] result;
        real upper_transformed;
        for (i in 1:length) {
            upper_transformed = (upperbound - mu) / sigma;
            sampleNormalUpperBound_RRR_lp(y_unit_normal[i], 0, 1, upper_transformed);
            result[i] = mu + sigma * y_unit_normal[i];
        }
        return result;
    }
    
    vector getReparameterizedNormalUpperBound_VRA_lp(vector y_unit_normal, real mu, array[] real sigma, real upperbound)
    {
        int length = num_elements(y_unit_normal);
        vector[length] result;
        real upper_transformed;
        if (num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            upper_transformed = (upperbound - mu) / sigma[i];
            sampleNormalUpperBound_RRR_lp(y_unit_normal[i], 0, 1, upper_transformed);
            result[i] = mu + sigma[i] * y_unit_normal[i];
        }
        return result;
    }
    
    vector getReparameterizedNormalUpperBound_VRV_lp(vector y_unit_normal, real mu, vector sigma, real upperbound)
    {
        int length = num_elements(y_unit_normal);
        vector[length] result;
        real upper_transformed;
        if (num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            upper_transformed = (upperbound - mu) / sigma[i];
            sampleNormalUpperBound_RRR_lp(y_unit_normal[i], 0, 1, upper_transformed);
            result[i] = mu + sigma[i] * y_unit_normal[i];
        }
        return result;
    }
    
    vector getReparameterizedNormalUpperBound_VAR_lp(vector y_unit_normal, array[] real mu, real sigma, real upperbound)
    {
        int length = num_elements(y_unit_normal);
        vector[length] result;
        real upper_transformed;
        if (num_elements(mu) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            upper_transformed = (upperbound - mu[i]) / sigma;
            sampleNormalUpperBound_RRR_lp(y_unit_normal[i], 0, 1, upper_transformed);
            result[i] = mu[i] + sigma * y_unit_normal[i];
        }
        return result;
    }
    
    vector getReparameterizedNormalUpperBound_VAA_lp(vector y_unit_normal, array[] real mu, array[] real sigma, real upperbound)
    {
        int length = num_elements(y_unit_normal);
        vector[length] result;
        real upper_transformed;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            upper_transformed = (upperbound - mu[i]) / sigma[i];
            sampleNormalUpperBound_RRR_lp(y_unit_normal[i], 0, 1, upper_transformed);
            result[i] = mu[i] + sigma[i] * y_unit_normal[i];
        }
        return result;
    }
    
    vector getReparameterizedNormalUpperBound_VAV_lp(vector y_unit_normal, array[] real mu, vector sigma, real upperbound)
    {
        int length = num_elements(y_unit_normal);
        vector[length] result;
        real upper_transformed;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            upper_transformed = (upperbound - mu[i]) / sigma[i];
            sampleNormalUpperBound_RRR_lp(y_unit_normal[i], 0, 1, upper_transformed);
            result[i] = mu[i] + sigma[i] * y_unit_normal[i];
        }
        return result;
    }
    
    vector getReparameterizedNormalUpperBound_VVR_lp(vector y_unit_normal, vector mu, real sigma, real upperbound)
    {
        int length = num_elements(y_unit_normal);
        vector[length] result;
        real upper_transformed;
        if (num_elements(mu) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            upper_transformed = (upperbound - mu[i]) / sigma;
            sampleNormalUpperBound_RRR_lp(y_unit_normal[i], 0, 1, upper_transformed);
            result[i] = mu[i] + sigma * y_unit_normal[i];
        }
        return result;
    }
    
    vector getReparameterizedNormalUpperBound_VVA_lp(vector y_unit_normal, vector mu, array[] real sigma, real upperbound)
    {
        int length = num_elements(y_unit_normal);
        vector[length] result;
        real upper_transformed;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            upper_transformed = (upperbound - mu[i]) / sigma[i];
            sampleNormalUpperBound_RRR_lp(y_unit_normal[i], 0, 1, upper_transformed);
            result[i] = mu[i] + sigma[i] * y_unit_normal[i];
        }
        return result;
    }
    
    vector getReparameterizedNormalUpperBound_VVV_lp(vector y_unit_normal, vector mu, vector sigma, real upperbound)
    {
        int length = num_elements(y_unit_normal);
        vector[length] result;
        real upper_transformed;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            upper_transformed = (upperbound - mu[i]) / sigma[i];
            sampleNormalUpperBound_RRR_lp(y_unit_normal[i], 0, 1, upper_transformed);
            result[i] = mu[i] + sigma[i] * y_unit_normal[i];
        }
        return result;
    }
    
    // With range (lower and upper) bounds

    real getReparameterizedNormalRangeBound_RRR_lp(real y_unit_normal, real mu, real sigma, real lowerbound, real upperbound)
    {
        real result;
        real lower_transformed;
        real upper_transformed;
        lower_transformed = (lowerbound - mu) / sigma;
        upper_transformed = (upperbound - mu) / sigma;
        sampleNormalRangeBound_RRR_lp(y_unit_normal, 0, 1, lower_transformed, upper_transformed);
        result = mu + sigma * y_unit_normal;
        return result;
    }
    
    array[] real getReparameterizedNormalRangeBound_ARR_lp(array[] real y_unit_normal, real mu, real sigma, real lowerbound, real upperbound)
    {
        int length = num_elements(y_unit_normal);
        array[length] real result;
        real lower_transformed;
        real upper_transformed;
        for (i in 1:length) {
            lower_transformed = (lowerbound - mu) / sigma;
            upper_transformed = (upperbound - mu) / sigma;
            sampleNormalRangeBound_RRR_lp(y_unit_normal[i], 0, 1, lower_transformed, upper_transformed);
            result[i] = mu + sigma * y_unit_normal[i];
        }
        return result;
    }
    
    array[] real getReparameterizedNormalRangeBound_ARA_lp(array[] real y_unit_normal, real mu, array[] real sigma, real lowerbound, real upperbound)
    {
        int length = num_elements(y_unit_normal);
        array[length] real result;
        real lower_transformed;
        real upper_transformed;
        if (num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = (lowerbound - mu) / sigma[i];
            upper_transformed = (upperbound - mu) / sigma[i];
            sampleNormalRangeBound_RRR_lp(y_unit_normal[i], 0, 1, lower_transformed, upper_transformed);
            result[i] = mu + sigma[i] * y_unit_normal[i];
        }
        return result;
    }
    
    array[] real getReparameterizedNormalRangeBound_ARV_lp(array[] real y_unit_normal, real mu, vector sigma, real lowerbound, real upperbound)
    {
        int length = num_elements(y_unit_normal);
        array[length] real result;
        real lower_transformed;
        real upper_transformed;
        if (num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = (lowerbound - mu) / sigma[i];
            upper_transformed = (upperbound - mu) / sigma[i];
            sampleNormalRangeBound_RRR_lp(y_unit_normal[i], 0, 1, lower_transformed, upper_transformed);
            result[i] = mu + sigma[i] * y_unit_normal[i];
        }
        return result;
    }
    
    array[] real getReparameterizedNormalRangeBound_AAR_lp(array[] real y_unit_normal, array[] real mu, real sigma, real lowerbound, real upperbound)
    {
        int length = num_elements(y_unit_normal);
        array[length] real result;
        real lower_transformed;
        real upper_transformed;
        if (num_elements(mu) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = (lowerbound - mu[i]) / sigma;
            upper_transformed = (upperbound - mu[i]) / sigma;
            sampleNormalRangeBound_RRR_lp(y_unit_normal[i], 0, 1, lower_transformed, upper_transformed);
            result[i] = mu[i] + sigma * y_unit_normal[i];
        }
        return result;
    }
    
    array[] real getReparameterizedNormalRangeBound_AAA_lp(array[] real y_unit_normal, array[] real mu, array[] real sigma, real lowerbound, real upperbound)
    {
        int length = num_elements(y_unit_normal);
        array[length] real result;
        real lower_transformed;
        real upper_transformed;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = (lowerbound - mu[i]) / sigma[i];
            upper_transformed = (upperbound - mu[i]) / sigma[i];
            sampleNormalRangeBound_RRR_lp(y_unit_normal[i], 0, 1, lower_transformed, upper_transformed);
            result[i] = mu[i] + sigma[i] * y_unit_normal[i];
        }
        return result;
    }
    
    array[] real getReparameterizedNormalRangeBound_AAV_lp(array[] real y_unit_normal, array[] real mu, vector sigma, real lowerbound, real upperbound)
    {
        int length = num_elements(y_unit_normal);
        array[length] real result;
        real lower_transformed;
        real upper_transformed;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = (lowerbound - mu[i]) / sigma[i];
            upper_transformed = (upperbound - mu[i]) / sigma[i];
            sampleNormalRangeBound_RRR_lp(y_unit_normal[i], 0, 1, lower_transformed, upper_transformed);
            result[i] = mu[i] + sigma[i] * y_unit_normal[i];
        }
        return result;
    }
    
    array[] real getReparameterizedNormalRangeBound_AVR_lp(array[] real y_unit_normal, vector mu, real sigma, real lowerbound, real upperbound)
    {
        int length = num_elements(y_unit_normal);
        array[length] real result;
        real lower_transformed;
        real upper_transformed;
        if (num_elements(mu) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = (lowerbound - mu[i]) / sigma;
            upper_transformed = (upperbound - mu[i]) / sigma;
            sampleNormalRangeBound_RRR_lp(y_unit_normal[i], 0, 1, lower_transformed, upper_transformed);
            result[i] = mu[i] + sigma * y_unit_normal[i];
        }
        return result;
    }
    
    array[] real getReparameterizedNormalRangeBound_AVA_lp(array[] real y_unit_normal, vector mu, array[] real sigma, real lowerbound, real upperbound)
    {
        int length = num_elements(y_unit_normal);
        array[length] real result;
        real lower_transformed;
        real upper_transformed;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = (lowerbound - mu[i]) / sigma[i];
            upper_transformed = (upperbound - mu[i]) / sigma[i];
            sampleNormalRangeBound_RRR_lp(y_unit_normal[i], 0, 1, lower_transformed, upper_transformed);
            result[i] = mu[i] + sigma[i] * y_unit_normal[i];
        }
        return result;
    }
    
    array[] real getReparameterizedNormalRangeBound_AVV_lp(array[] real y_unit_normal, vector mu, vector sigma, real lowerbound, real upperbound)
    {
        int length = num_elements(y_unit_normal);
        array[length] real result;
        real lower_transformed;
        real upper_transformed;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = (lowerbound - mu[i]) / sigma[i];
            upper_transformed = (upperbound - mu[i]) / sigma[i];
            sampleNormalRangeBound_RRR_lp(y_unit_normal[i], 0, 1, lower_transformed, upper_transformed);
            result[i] = mu[i] + sigma[i] * y_unit_normal[i];
        }
        return result;
    }
    
    array[,] real getReparameterizedNormalRangeBound_2RR_lp(array[,] real y_unit_normal, real mu, real sigma, real lowerbound, real upperbound)
    {
        array[2] int dimensions = dims(y_unit_normal);
        array[dimensions[1], dimensions[2]] real result;
        real lower_transformed;
        real upper_transformed;
        lower_transformed = (lowerbound - mu) / sigma;
        upper_transformed = (upperbound - mu) / sigma;
        for (i in 1:dimensions[1]) {
            for (j in 1:dimensions[2]) {
                sampleNormalRangeBound_RRR_lp(y_unit_normal[i, j], 0, 1, lower_transformed, upper_transformed);
                result[i, j] = mu + sigma * y_unit_normal[i, j];
            }
        }
        return result;
    }
    
    array[,,] real getReparameterizedNormalRangeBound_3RR_lp(array[,,] real y_unit_normal, real mu, real sigma, real lowerbound, real upperbound)
    {
        array[3] int dimensions = dims(y_unit_normal);
        array[dimensions[1], dimensions[2], dimensions[3]] real result;
        real lower_transformed;
        real upper_transformed;
        lower_transformed = (lowerbound - mu) / sigma;
        upper_transformed = (upperbound - mu) / sigma;
        for (i in 1:dimensions[1]) {
            for (j in 1:dimensions[2]) {
                for (k in 1:dimensions[3]) {
                    sampleNormalRangeBound_RRR_lp(y_unit_normal[i, j, k], 0, 1, lower_transformed, upper_transformed);
                    result[i, j, k] = mu + sigma * y_unit_normal[i, j, k];
                }
            }
        }
        return result;
    }
    
    vector getReparameterizedNormalRangeBound_VRR_lp(vector y_unit_normal, real mu, real sigma, real lowerbound, real upperbound)
    {
        int length = num_elements(y_unit_normal);
        vector[length] result;
        real lower_transformed;
        real upper_transformed;
        for (i in 1:length) {
            lower_transformed = (lowerbound - mu) / sigma;
            upper_transformed = (upperbound - mu) / sigma;
            sampleNormalRangeBound_RRR_lp(y_unit_normal[i], 0, 1, lower_transformed, upper_transformed);
            result[i] = mu + sigma * y_unit_normal[i];
        }
        return result;
    }
    
    vector getReparameterizedNormalRangeBound_VRA_lp(vector y_unit_normal, real mu, array[] real sigma, real lowerbound, real upperbound)
    {
        int length = num_elements(y_unit_normal);
        vector[length] result;
        real lower_transformed;
        real upper_transformed;
        if (num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = (lowerbound - mu) / sigma[i];
            upper_transformed = (upperbound - mu) / sigma[i];
            sampleNormalRangeBound_RRR_lp(y_unit_normal[i], 0, 1, lower_transformed, upper_transformed);
            result[i] = mu + sigma[i] * y_unit_normal[i];
        }
        return result;
    }
    
    vector getReparameterizedNormalRangeBound_VRV_lp(vector y_unit_normal, real mu, vector sigma, real lowerbound, real upperbound)
    {
        int length = num_elements(y_unit_normal);
        vector[length] result;
        real lower_transformed;
        real upper_transformed;
        if (num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = (lowerbound - mu) / sigma[i];
            upper_transformed = (upperbound - mu) / sigma[i];
            sampleNormalRangeBound_RRR_lp(y_unit_normal[i], 0, 1, lower_transformed, upper_transformed);
            result[i] = mu + sigma[i] * y_unit_normal[i];
        }
        return result;
    }
    
    vector getReparameterizedNormalRangeBound_VAR_lp(vector y_unit_normal, array[] real mu, real sigma, real lowerbound, real upperbound)
    {
        int length = num_elements(y_unit_normal);
        vector[length] result;
        real lower_transformed;
        real upper_transformed;
        if (num_elements(mu) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = (lowerbound - mu[i]) / sigma;
            upper_transformed = (upperbound - mu[i]) / sigma;
            sampleNormalRangeBound_RRR_lp(y_unit_normal[i], 0, 1, lower_transformed, upper_transformed);
            result[i] = mu[i] + sigma * y_unit_normal[i];
        }
        return result;
    }
    
    vector getReparameterizedNormalRangeBound_VAA_lp(vector y_unit_normal, array[] real mu, array[] real sigma, real lowerbound, real upperbound)
    {
        int length = num_elements(y_unit_normal);
        vector[length] result;
        real lower_transformed;
        real upper_transformed;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = (lowerbound - mu[i]) / sigma[i];
            upper_transformed = (upperbound - mu[i]) / sigma[i];
            sampleNormalRangeBound_RRR_lp(y_unit_normal[i], 0, 1, lower_transformed, upper_transformed);
            result[i] = mu[i] + sigma[i] * y_unit_normal[i];
        }
        return result;
    }
    
    vector getReparameterizedNormalRangeBound_VAV_lp(vector y_unit_normal, array[] real mu, vector sigma, real lowerbound, real upperbound)
    {
        int length = num_elements(y_unit_normal);
        vector[length] result;
        real lower_transformed;
        real upper_transformed;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = (lowerbound - mu[i]) / sigma[i];
            upper_transformed = (upperbound - mu[i]) / sigma[i];
            sampleNormalRangeBound_RRR_lp(y_unit_normal[i], 0, 1, lower_transformed, upper_transformed);
            result[i] = mu[i] + sigma[i] * y_unit_normal[i];
        }
        return result;
    }
    
    vector getReparameterizedNormalRangeBound_VVR_lp(vector y_unit_normal, vector mu, real sigma, real lowerbound, real upperbound)
    {
        int length = num_elements(y_unit_normal);
        vector[length] result;
        real lower_transformed;
        real upper_transformed;
        if (num_elements(mu) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = (lowerbound - mu[i]) / sigma;
            upper_transformed = (upperbound - mu[i]) / sigma;
            sampleNormalRangeBound_RRR_lp(y_unit_normal[i], 0, 1, lower_transformed, upper_transformed);
            result[i] = mu[i] + sigma * y_unit_normal[i];
        }
        return result;
    }
    
    vector getReparameterizedNormalRangeBound_VVA_lp(vector y_unit_normal, vector mu, array[] real sigma, real lowerbound, real upperbound)
    {
        int length = num_elements(y_unit_normal);
        vector[length] result;
        real lower_transformed;
        real upper_transformed;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = (lowerbound - mu[i]) / sigma[i];
            upper_transformed = (upperbound - mu[i]) / sigma[i];
            sampleNormalRangeBound_RRR_lp(y_unit_normal[i], 0, 1, lower_transformed, upper_transformed);
            result[i] = mu[i] + sigma[i] * y_unit_normal[i];
        }
        return result;
    }
    
    vector getReparameterizedNormalRangeBound_VVV_lp(vector y_unit_normal, vector mu, vector sigma, real lowerbound, real upperbound)
    {
        int length = num_elements(y_unit_normal);
        vector[length] result;
        real lower_transformed;
        real upper_transformed;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = (lowerbound - mu[i]) / sigma[i];
            upper_transformed = (upperbound - mu[i]) / sigma[i];
            sampleNormalRangeBound_RRR_lp(y_unit_normal[i], 0, 1, lower_transformed, upper_transformed);
            result[i] = mu[i] + sigma[i] * y_unit_normal[i];
        }
        return result;
    }
    
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    // Cauchy distribution, reparameterized to the uniform distribution
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    /*
    See p339 of the Stan (2017) manual.
    The transformation is

        y_cauchy(mu, sigma) = mu + sigma * y_uniform(-pi/2, pi/2)

    so the boundary transformation is the reverse, namely

        uniform_boundary = arctan((boundary - mu) / sigma)

    Note that
        arctan(-infinity) is -pi/2; arctan(0) is 0; arctan(infinity) is pi/2
        ... so for lower == 0, lower_transformed == 0

    We can do the range constraints like this:
         real<lower=-pi()/2, upper=pi()/2> y_uniform;  // Cauchy
         real<lower=0, upper=pi()/2> y_uniform;  // positive half-Cauchy
    and the sampling statement would be unnecessary, but we're going to
    to the sampling using "target +=" so that bridgesampling works.

    You might think that because of that, the range constraint is unnecessary,
    but it IS STILL NECESSARY or Stan will explore invalid ranges.

    */

    real reparameterizedCauchyBoundary(real boundary, real mu, real sigma)
    {
        // boundary: in real-world Cauchy(mu, sigma) space
        // return value: equivalent in the reparameterized uniform [-pi/2, +pi/2] space
        return atan((boundary - mu) / sigma);
    }
    
    // Plain

    real getReparameterizedCauchy_RRR_lp(real y_uniform, real mu, real sigma)
    {
        real result;
        sampleUniform_RRR_lp(y_uniform, -pi()/2, pi()/2);
        result = mu + sigma * tan(y_uniform);
        return result;
    }
    
    array[] real getReparameterizedCauchy_ARR_lp(array[] real y_uniform, real mu, real sigma)
    {
        int length = num_elements(y_uniform);
        array[length] real result;
        for (i in 1:length) {
            sampleUniform_RRR_lp(y_uniform[i], -pi()/2, pi()/2);
            result[i] = mu + sigma * tan(y_uniform[i]);
        }
        return result;
    }
    
    array[] real getReparameterizedCauchy_ARA_lp(array[] real y_uniform, real mu, array[] real sigma)
    {
        int length = num_elements(y_uniform);
        array[length] real result;
        if (num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            sampleUniform_RRR_lp(y_uniform[i], -pi()/2, pi()/2);
            result[i] = mu + sigma[i] * tan(y_uniform[i]);
        }
        return result;
    }
    
    array[] real getReparameterizedCauchy_ARV_lp(array[] real y_uniform, real mu, vector sigma)
    {
        int length = num_elements(y_uniform);
        array[length] real result;
        if (num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            sampleUniform_RRR_lp(y_uniform[i], -pi()/2, pi()/2);
            result[i] = mu + sigma[i] * tan(y_uniform[i]);
        }
        return result;
    }
    
    array[] real getReparameterizedCauchy_AAR_lp(array[] real y_uniform, array[] real mu, real sigma)
    {
        int length = num_elements(y_uniform);
        array[length] real result;
        if (num_elements(mu) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            sampleUniform_RRR_lp(y_uniform[i], -pi()/2, pi()/2);
            result[i] = mu[i] + sigma * tan(y_uniform[i]);
        }
        return result;
    }
    
    array[] real getReparameterizedCauchy_AAA_lp(array[] real y_uniform, array[] real mu, array[] real sigma)
    {
        int length = num_elements(y_uniform);
        array[length] real result;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            sampleUniform_RRR_lp(y_uniform[i], -pi()/2, pi()/2);
            result[i] = mu[i] + sigma[i] * tan(y_uniform[i]);
        }
        return result;
    }
    
    array[] real getReparameterizedCauchy_AAV_lp(array[] real y_uniform, array[] real mu, vector sigma)
    {
        int length = num_elements(y_uniform);
        array[length] real result;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            sampleUniform_RRR_lp(y_uniform[i], -pi()/2, pi()/2);
            result[i] = mu[i] + sigma[i] * tan(y_uniform[i]);
        }
        return result;
    }
    
    array[] real getReparameterizedCauchy_AVR_lp(array[] real y_uniform, vector mu, real sigma)
    {
        int length = num_elements(y_uniform);
        array[length] real result;
        if (num_elements(mu) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            sampleUniform_RRR_lp(y_uniform[i], -pi()/2, pi()/2);
            result[i] = mu[i] + sigma * tan(y_uniform[i]);
        }
        return result;
    }
    
    array[] real getReparameterizedCauchy_AVA_lp(array[] real y_uniform, vector mu, array[] real sigma)
    {
        int length = num_elements(y_uniform);
        array[length] real result;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            sampleUniform_RRR_lp(y_uniform[i], -pi()/2, pi()/2);
            result[i] = mu[i] + sigma[i] * tan(y_uniform[i]);
        }
        return result;
    }
    
    array[] real getReparameterizedCauchy_AVV_lp(array[] real y_uniform, vector mu, vector sigma)
    {
        int length = num_elements(y_uniform);
        array[length] real result;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            sampleUniform_RRR_lp(y_uniform[i], -pi()/2, pi()/2);
            result[i] = mu[i] + sigma[i] * tan(y_uniform[i]);
        }
        return result;
    }
    
    array[,] real getReparameterizedCauchy_2RR_lp(array[,] real y_uniform, real mu, real sigma)
    {
        array[2] int dimensions = dims(y_uniform);
        array[dimensions[1], dimensions[2]] real result;
        for (i in 1:dimensions[1]) {
            for (j in 1:dimensions[2]) {
                sampleUniform_RRR_lp(y_uniform[i, j], -pi()/2, pi()/2);
                result[i, j] = mu + sigma * tan(y_uniform[i, j]);
            }
        }
        return result;
    }
    
    array[,,] real getReparameterizedCauchy_3RR_lp(array[,,] real y_uniform, real mu, real sigma)
    {
        array[3] int dimensions = dims(y_uniform);
        array[dimensions[1], dimensions[2], dimensions[3]] real result;
        for (i in 1:dimensions[1]) {
            for (j in 1:dimensions[2]) {
                for (k in 1:dimensions[3]) {
                    sampleUniform_RRR_lp(y_uniform[i, j, k], -pi()/2, pi()/2);
                    result[i, j, k] = mu + sigma * tan(y_uniform[i, j, k]);
                }
            }
        }
        return result;
    }
    
    vector getReparameterizedCauchy_VRR_lp(vector y_uniform, real mu, real sigma)
    {
        int length = num_elements(y_uniform);
        vector[length] result;
        for (i in 1:length) {
            sampleUniform_RRR_lp(y_uniform[i], -pi()/2, pi()/2);
            result[i] = mu + sigma * tan(y_uniform[i]);
        }
        return result;
    }
    
    vector getReparameterizedCauchy_VRA_lp(vector y_uniform, real mu, array[] real sigma)
    {
        int length = num_elements(y_uniform);
        vector[length] result;
        if (num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            sampleUniform_RRR_lp(y_uniform[i], -pi()/2, pi()/2);
            result[i] = mu + sigma[i] * tan(y_uniform[i]);
        }
        return result;
    }
    
    vector getReparameterizedCauchy_VRV_lp(vector y_uniform, real mu, vector sigma)
    {
        int length = num_elements(y_uniform);
        vector[length] result;
        if (num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            sampleUniform_RRR_lp(y_uniform[i], -pi()/2, pi()/2);
            result[i] = mu + sigma[i] * tan(y_uniform[i]);
        }
        return result;
    }
    
    vector getReparameterizedCauchy_VAR_lp(vector y_uniform, array[] real mu, real sigma)
    {
        int length = num_elements(y_uniform);
        vector[length] result;
        if (num_elements(mu) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            sampleUniform_RRR_lp(y_uniform[i], -pi()/2, pi()/2);
            result[i] = mu[i] + sigma * tan(y_uniform[i]);
        }
        return result;
    }
    
    vector getReparameterizedCauchy_VAA_lp(vector y_uniform, array[] real mu, array[] real sigma)
    {
        int length = num_elements(y_uniform);
        vector[length] result;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            sampleUniform_RRR_lp(y_uniform[i], -pi()/2, pi()/2);
            result[i] = mu[i] + sigma[i] * tan(y_uniform[i]);
        }
        return result;
    }
    
    vector getReparameterizedCauchy_VAV_lp(vector y_uniform, array[] real mu, vector sigma)
    {
        int length = num_elements(y_uniform);
        vector[length] result;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            sampleUniform_RRR_lp(y_uniform[i], -pi()/2, pi()/2);
            result[i] = mu[i] + sigma[i] * tan(y_uniform[i]);
        }
        return result;
    }
    
    vector getReparameterizedCauchy_VVR_lp(vector y_uniform, vector mu, real sigma)
    {
        int length = num_elements(y_uniform);
        vector[length] result;
        if (num_elements(mu) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            sampleUniform_RRR_lp(y_uniform[i], -pi()/2, pi()/2);
            result[i] = mu[i] + sigma * tan(y_uniform[i]);
        }
        return result;
    }
    
    vector getReparameterizedCauchy_VVA_lp(vector y_uniform, vector mu, array[] real sigma)
    {
        int length = num_elements(y_uniform);
        vector[length] result;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            sampleUniform_RRR_lp(y_uniform[i], -pi()/2, pi()/2);
            result[i] = mu[i] + sigma[i] * tan(y_uniform[i]);
        }
        return result;
    }
    
    vector getReparameterizedCauchy_VVV_lp(vector y_uniform, vector mu, vector sigma)
    {
        int length = num_elements(y_uniform);
        vector[length] result;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            sampleUniform_RRR_lp(y_uniform[i], -pi()/2, pi()/2);
            result[i] = mu[i] + sigma[i] * tan(y_uniform[i]);
        }
        return result;
    }
    
    // With lower bound

    real getReparameterizedCauchyLowerBound_RRR_lp(real y_uniform, real mu, real sigma, real lowerbound)
    {
        real result;
        real lower_transformed;
        lower_transformed = atan((lowerbound - mu) / sigma);
        sampleUniform_RRR_lp(y_uniform, lower_transformed, pi()/2);
        result = mu + sigma * tan(y_uniform);
        return result;
    }
    
    array[] real getReparameterizedCauchyLowerBound_ARR_lp(array[] real y_uniform, real mu, real sigma, real lowerbound)
    {
        int length = num_elements(y_uniform);
        array[length] real result;
        real lower_transformed;
        for (i in 1:length) {
            lower_transformed = atan((lowerbound - mu) / sigma);
            sampleUniform_RRR_lp(y_uniform[i], lower_transformed, pi()/2);
            result[i] = mu + sigma * tan(y_uniform[i]);
        }
        return result;
    }
    
    array[] real getReparameterizedCauchyLowerBound_ARA_lp(array[] real y_uniform, real mu, array[] real sigma, real lowerbound)
    {
        int length = num_elements(y_uniform);
        array[length] real result;
        real lower_transformed;
        if (num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = atan((lowerbound - mu) / sigma[i]);
            sampleUniform_RRR_lp(y_uniform[i], lower_transformed, pi()/2);
            result[i] = mu + sigma[i] * tan(y_uniform[i]);
        }
        return result;
    }
    
    array[] real getReparameterizedCauchyLowerBound_ARV_lp(array[] real y_uniform, real mu, vector sigma, real lowerbound)
    {
        int length = num_elements(y_uniform);
        array[length] real result;
        real lower_transformed;
        if (num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = atan((lowerbound - mu) / sigma[i]);
            sampleUniform_RRR_lp(y_uniform[i], lower_transformed, pi()/2);
            result[i] = mu + sigma[i] * tan(y_uniform[i]);
        }
        return result;
    }
    
    array[] real getReparameterizedCauchyLowerBound_AAR_lp(array[] real y_uniform, array[] real mu, real sigma, real lowerbound)
    {
        int length = num_elements(y_uniform);
        array[length] real result;
        real lower_transformed;
        if (num_elements(mu) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = atan((lowerbound - mu[i]) / sigma);
            sampleUniform_RRR_lp(y_uniform[i], lower_transformed, pi()/2);
            result[i] = mu[i] + sigma * tan(y_uniform[i]);
        }
        return result;
    }
    
    array[] real getReparameterizedCauchyLowerBound_AAA_lp(array[] real y_uniform, array[] real mu, array[] real sigma, real lowerbound)
    {
        int length = num_elements(y_uniform);
        array[length] real result;
        real lower_transformed;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = atan((lowerbound - mu[i]) / sigma[i]);
            sampleUniform_RRR_lp(y_uniform[i], lower_transformed, pi()/2);
            result[i] = mu[i] + sigma[i] * tan(y_uniform[i]);
        }
        return result;
    }
    
    array[] real getReparameterizedCauchyLowerBound_AAV_lp(array[] real y_uniform, array[] real mu, vector sigma, real lowerbound)
    {
        int length = num_elements(y_uniform);
        array[length] real result;
        real lower_transformed;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = atan((lowerbound - mu[i]) / sigma[i]);
            sampleUniform_RRR_lp(y_uniform[i], lower_transformed, pi()/2);
            result[i] = mu[i] + sigma[i] * tan(y_uniform[i]);
        }
        return result;
    }
    
    array[] real getReparameterizedCauchyLowerBound_AVR_lp(array[] real y_uniform, vector mu, real sigma, real lowerbound)
    {
        int length = num_elements(y_uniform);
        array[length] real result;
        real lower_transformed;
        if (num_elements(mu) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = atan((lowerbound - mu[i]) / sigma);
            sampleUniform_RRR_lp(y_uniform[i], lower_transformed, pi()/2);
            result[i] = mu[i] + sigma * tan(y_uniform[i]);
        }
        return result;
    }
    
    array[] real getReparameterizedCauchyLowerBound_AVA_lp(array[] real y_uniform, vector mu, array[] real sigma, real lowerbound)
    {
        int length = num_elements(y_uniform);
        array[length] real result;
        real lower_transformed;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = atan((lowerbound - mu[i]) / sigma[i]);
            sampleUniform_RRR_lp(y_uniform[i], lower_transformed, pi()/2);
            result[i] = mu[i] + sigma[i] * tan(y_uniform[i]);
        }
        return result;
    }
    
    array[] real getReparameterizedCauchyLowerBound_AVV_lp(array[] real y_uniform, vector mu, vector sigma, real lowerbound)
    {
        int length = num_elements(y_uniform);
        array[length] real result;
        real lower_transformed;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = atan((lowerbound - mu[i]) / sigma[i]);
            sampleUniform_RRR_lp(y_uniform[i], lower_transformed, pi()/2);
            result[i] = mu[i] + sigma[i] * tan(y_uniform[i]);
        }
        return result;
    }
    
    array[,] real getReparameterizedCauchyLowerBound_2RR_lp(array[,] real y_uniform, real mu, real sigma, real lowerbound)
    {
        array[2] int dimensions = dims(y_uniform);
        array[dimensions[1], dimensions[2]] real result;
        real lower_transformed;
        lower_transformed = atan((lowerbound - mu) / sigma);
        for (i in 1:dimensions[1]) {
            for (j in 1:dimensions[2]) {
                sampleUniform_RRR_lp(y_uniform[i, j], lower_transformed, pi()/2);
                result[i, j] = mu + sigma * tan(y_uniform[i, j]);
            }
        }
        return result;
    }
    
    array[,,] real getReparameterizedCauchyLowerBound_3RR_lp(array[,,] real y_uniform, real mu, real sigma, real lowerbound)
    {
        array[3] int dimensions = dims(y_uniform);
        array[dimensions[1], dimensions[2], dimensions[3]] real result;
        real lower_transformed;
        lower_transformed = atan((lowerbound - mu) / sigma);
        for (i in 1:dimensions[1]) {
            for (j in 1:dimensions[2]) {
                for (k in 1:dimensions[3]) {
                    sampleUniform_RRR_lp(y_uniform[i, j, k], lower_transformed, pi()/2);
                    result[i, j, k] = mu + sigma * tan(y_uniform[i, j, k]);
                }
            }
        }
        return result;
    }
    
    vector getReparameterizedCauchyLowerBound_VRR_lp(vector y_uniform, real mu, real sigma, real lowerbound)
    {
        int length = num_elements(y_uniform);
        vector[length] result;
        real lower_transformed;
        for (i in 1:length) {
            lower_transformed = atan((lowerbound - mu) / sigma);
            sampleUniform_RRR_lp(y_uniform[i], lower_transformed, pi()/2);
            result[i] = mu + sigma * tan(y_uniform[i]);
        }
        return result;
    }
    
    vector getReparameterizedCauchyLowerBound_VRA_lp(vector y_uniform, real mu, array[] real sigma, real lowerbound)
    {
        int length = num_elements(y_uniform);
        vector[length] result;
        real lower_transformed;
        if (num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = atan((lowerbound - mu) / sigma[i]);
            sampleUniform_RRR_lp(y_uniform[i], lower_transformed, pi()/2);
            result[i] = mu + sigma[i] * tan(y_uniform[i]);
        }
        return result;
    }
    
    vector getReparameterizedCauchyLowerBound_VRV_lp(vector y_uniform, real mu, vector sigma, real lowerbound)
    {
        int length = num_elements(y_uniform);
        vector[length] result;
        real lower_transformed;
        if (num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = atan((lowerbound - mu) / sigma[i]);
            sampleUniform_RRR_lp(y_uniform[i], lower_transformed, pi()/2);
            result[i] = mu + sigma[i] * tan(y_uniform[i]);
        }
        return result;
    }
    
    vector getReparameterizedCauchyLowerBound_VAR_lp(vector y_uniform, array[] real mu, real sigma, real lowerbound)
    {
        int length = num_elements(y_uniform);
        vector[length] result;
        real lower_transformed;
        if (num_elements(mu) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = atan((lowerbound - mu[i]) / sigma);
            sampleUniform_RRR_lp(y_uniform[i], lower_transformed, pi()/2);
            result[i] = mu[i] + sigma * tan(y_uniform[i]);
        }
        return result;
    }
    
    vector getReparameterizedCauchyLowerBound_VAA_lp(vector y_uniform, array[] real mu, array[] real sigma, real lowerbound)
    {
        int length = num_elements(y_uniform);
        vector[length] result;
        real lower_transformed;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = atan((lowerbound - mu[i]) / sigma[i]);
            sampleUniform_RRR_lp(y_uniform[i], lower_transformed, pi()/2);
            result[i] = mu[i] + sigma[i] * tan(y_uniform[i]);
        }
        return result;
    }
    
    vector getReparameterizedCauchyLowerBound_VAV_lp(vector y_uniform, array[] real mu, vector sigma, real lowerbound)
    {
        int length = num_elements(y_uniform);
        vector[length] result;
        real lower_transformed;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = atan((lowerbound - mu[i]) / sigma[i]);
            sampleUniform_RRR_lp(y_uniform[i], lower_transformed, pi()/2);
            result[i] = mu[i] + sigma[i] * tan(y_uniform[i]);
        }
        return result;
    }
    
    vector getReparameterizedCauchyLowerBound_VVR_lp(vector y_uniform, vector mu, real sigma, real lowerbound)
    {
        int length = num_elements(y_uniform);
        vector[length] result;
        real lower_transformed;
        if (num_elements(mu) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = atan((lowerbound - mu[i]) / sigma);
            sampleUniform_RRR_lp(y_uniform[i], lower_transformed, pi()/2);
            result[i] = mu[i] + sigma * tan(y_uniform[i]);
        }
        return result;
    }
    
    vector getReparameterizedCauchyLowerBound_VVA_lp(vector y_uniform, vector mu, array[] real sigma, real lowerbound)
    {
        int length = num_elements(y_uniform);
        vector[length] result;
        real lower_transformed;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = atan((lowerbound - mu[i]) / sigma[i]);
            sampleUniform_RRR_lp(y_uniform[i], lower_transformed, pi()/2);
            result[i] = mu[i] + sigma[i] * tan(y_uniform[i]);
        }
        return result;
    }
    
    vector getReparameterizedCauchyLowerBound_VVV_lp(vector y_uniform, vector mu, vector sigma, real lowerbound)
    {
        int length = num_elements(y_uniform);
        vector[length] result;
        real lower_transformed;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = atan((lowerbound - mu[i]) / sigma[i]);
            sampleUniform_RRR_lp(y_uniform[i], lower_transformed, pi()/2);
            result[i] = mu[i] + sigma[i] * tan(y_uniform[i]);
        }
        return result;
    }
    
    // With upper bound

    real getReparameterizedCauchyUpperBound_RRR_lp(real y_uniform, real mu, real sigma, real upperbound)
    {
        real result;
        real upper_transformed;
        upper_transformed = atan((upperbound - mu) / sigma);
        sampleUniform_RRR_lp(y_uniform, -pi()/2, upper_transformed);
        result = mu + sigma * tan(y_uniform);
        return result;
    }
    
    array[] real getReparameterizedCauchyUpperBound_ARR_lp(array[] real y_uniform, real mu, real sigma, real upperbound)
    {
        int length = num_elements(y_uniform);
        array[length] real result;
        real upper_transformed;
        for (i in 1:length) {
            upper_transformed = atan((upperbound - mu) / sigma);
            sampleUniform_RRR_lp(y_uniform[i], -pi()/2, upper_transformed);
            result[i] = mu + sigma * tan(y_uniform[i]);
        }
        return result;
    }
    
    array[] real getReparameterizedCauchyUpperBound_ARA_lp(array[] real y_uniform, real mu, array[] real sigma, real upperbound)
    {
        int length = num_elements(y_uniform);
        array[length] real result;
        real upper_transformed;
        if (num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            upper_transformed = atan((upperbound - mu) / sigma[i]);
            sampleUniform_RRR_lp(y_uniform[i], -pi()/2, upper_transformed);
            result[i] = mu + sigma[i] * tan(y_uniform[i]);
        }
        return result;
    }
    
    array[] real getReparameterizedCauchyUpperBound_ARV_lp(array[] real y_uniform, real mu, vector sigma, real upperbound)
    {
        int length = num_elements(y_uniform);
        array[length] real result;
        real upper_transformed;
        if (num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            upper_transformed = atan((upperbound - mu) / sigma[i]);
            sampleUniform_RRR_lp(y_uniform[i], -pi()/2, upper_transformed);
            result[i] = mu + sigma[i] * tan(y_uniform[i]);
        }
        return result;
    }
    
    array[] real getReparameterizedCauchyUpperBound_AAR_lp(array[] real y_uniform, array[] real mu, real sigma, real upperbound)
    {
        int length = num_elements(y_uniform);
        array[length] real result;
        real upper_transformed;
        if (num_elements(mu) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            upper_transformed = atan((upperbound - mu[i]) / sigma);
            sampleUniform_RRR_lp(y_uniform[i], -pi()/2, upper_transformed);
            result[i] = mu[i] + sigma * tan(y_uniform[i]);
        }
        return result;
    }
    
    array[] real getReparameterizedCauchyUpperBound_AAA_lp(array[] real y_uniform, array[] real mu, array[] real sigma, real upperbound)
    {
        int length = num_elements(y_uniform);
        array[length] real result;
        real upper_transformed;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            upper_transformed = atan((upperbound - mu[i]) / sigma[i]);
            sampleUniform_RRR_lp(y_uniform[i], -pi()/2, upper_transformed);
            result[i] = mu[i] + sigma[i] * tan(y_uniform[i]);
        }
        return result;
    }
    
    array[] real getReparameterizedCauchyUpperBound_AAV_lp(array[] real y_uniform, array[] real mu, vector sigma, real upperbound)
    {
        int length = num_elements(y_uniform);
        array[length] real result;
        real upper_transformed;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            upper_transformed = atan((upperbound - mu[i]) / sigma[i]);
            sampleUniform_RRR_lp(y_uniform[i], -pi()/2, upper_transformed);
            result[i] = mu[i] + sigma[i] * tan(y_uniform[i]);
        }
        return result;
    }
    
    array[] real getReparameterizedCauchyUpperBound_AVR_lp(array[] real y_uniform, vector mu, real sigma, real upperbound)
    {
        int length = num_elements(y_uniform);
        array[length] real result;
        real upper_transformed;
        if (num_elements(mu) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            upper_transformed = atan((upperbound - mu[i]) / sigma);
            sampleUniform_RRR_lp(y_uniform[i], -pi()/2, upper_transformed);
            result[i] = mu[i] + sigma * tan(y_uniform[i]);
        }
        return result;
    }
    
    array[] real getReparameterizedCauchyUpperBound_AVA_lp(array[] real y_uniform, vector mu, array[] real sigma, real upperbound)
    {
        int length = num_elements(y_uniform);
        array[length] real result;
        real upper_transformed;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            upper_transformed = atan((upperbound - mu[i]) / sigma[i]);
            sampleUniform_RRR_lp(y_uniform[i], -pi()/2, upper_transformed);
            result[i] = mu[i] + sigma[i] * tan(y_uniform[i]);
        }
        return result;
    }
    
    array[] real getReparameterizedCauchyUpperBound_AVV_lp(array[] real y_uniform, vector mu, vector sigma, real upperbound)
    {
        int length = num_elements(y_uniform);
        array[length] real result;
        real upper_transformed;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            upper_transformed = atan((upperbound - mu[i]) / sigma[i]);
            sampleUniform_RRR_lp(y_uniform[i], -pi()/2, upper_transformed);
            result[i] = mu[i] + sigma[i] * tan(y_uniform[i]);
        }
        return result;
    }
    
    array[,] real getReparameterizedCauchyUpperBound_2RR_lp(array[,] real y_uniform, real mu, real sigma, real upperbound)
    {
        array[2] int dimensions = dims(y_uniform);
        array[dimensions[1], dimensions[2]] real result;
        real upper_transformed;
        upper_transformed = atan((upperbound - mu) / sigma);
        for (i in 1:dimensions[1]) {
            for (j in 1:dimensions[2]) {
                sampleUniform_RRR_lp(y_uniform[i, j], -pi()/2, upper_transformed);
                result[i, j] = mu + sigma * tan(y_uniform[i, j]);
            }
        }
        return result;
    }
    
    array[,,] real getReparameterizedCauchyUpperBound_3RR_lp(array[,,] real y_uniform, real mu, real sigma, real upperbound)
    {
        array[3] int dimensions = dims(y_uniform);
        array[dimensions[1], dimensions[2], dimensions[3]] real result;
        real upper_transformed;
        upper_transformed = atan((upperbound - mu) / sigma);
        for (i in 1:dimensions[1]) {
            for (j in 1:dimensions[2]) {
                for (k in 1:dimensions[3]) {
                    sampleUniform_RRR_lp(y_uniform[i, j, k], -pi()/2, upper_transformed);
                    result[i, j, k] = mu + sigma * tan(y_uniform[i, j, k]);
                }
            }
        }
        return result;
    }
    
    vector getReparameterizedCauchyUpperBound_VRR_lp(vector y_uniform, real mu, real sigma, real upperbound)
    {
        int length = num_elements(y_uniform);
        vector[length] result;
        real upper_transformed;
        for (i in 1:length) {
            upper_transformed = atan((upperbound - mu) / sigma);
            sampleUniform_RRR_lp(y_uniform[i], -pi()/2, upper_transformed);
            result[i] = mu + sigma * tan(y_uniform[i]);
        }
        return result;
    }
    
    vector getReparameterizedCauchyUpperBound_VRA_lp(vector y_uniform, real mu, array[] real sigma, real upperbound)
    {
        int length = num_elements(y_uniform);
        vector[length] result;
        real upper_transformed;
        if (num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            upper_transformed = atan((upperbound - mu) / sigma[i]);
            sampleUniform_RRR_lp(y_uniform[i], -pi()/2, upper_transformed);
            result[i] = mu + sigma[i] * tan(y_uniform[i]);
        }
        return result;
    }
    
    vector getReparameterizedCauchyUpperBound_VRV_lp(vector y_uniform, real mu, vector sigma, real upperbound)
    {
        int length = num_elements(y_uniform);
        vector[length] result;
        real upper_transformed;
        if (num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            upper_transformed = atan((upperbound - mu) / sigma[i]);
            sampleUniform_RRR_lp(y_uniform[i], -pi()/2, upper_transformed);
            result[i] = mu + sigma[i] * tan(y_uniform[i]);
        }
        return result;
    }
    
    vector getReparameterizedCauchyUpperBound_VAR_lp(vector y_uniform, array[] real mu, real sigma, real upperbound)
    {
        int length = num_elements(y_uniform);
        vector[length] result;
        real upper_transformed;
        if (num_elements(mu) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            upper_transformed = atan((upperbound - mu[i]) / sigma);
            sampleUniform_RRR_lp(y_uniform[i], -pi()/2, upper_transformed);
            result[i] = mu[i] + sigma * tan(y_uniform[i]);
        }
        return result;
    }
    
    vector getReparameterizedCauchyUpperBound_VAA_lp(vector y_uniform, array[] real mu, array[] real sigma, real upperbound)
    {
        int length = num_elements(y_uniform);
        vector[length] result;
        real upper_transformed;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            upper_transformed = atan((upperbound - mu[i]) / sigma[i]);
            sampleUniform_RRR_lp(y_uniform[i], -pi()/2, upper_transformed);
            result[i] = mu[i] + sigma[i] * tan(y_uniform[i]);
        }
        return result;
    }
    
    vector getReparameterizedCauchyUpperBound_VAV_lp(vector y_uniform, array[] real mu, vector sigma, real upperbound)
    {
        int length = num_elements(y_uniform);
        vector[length] result;
        real upper_transformed;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            upper_transformed = atan((upperbound - mu[i]) / sigma[i]);
            sampleUniform_RRR_lp(y_uniform[i], -pi()/2, upper_transformed);
            result[i] = mu[i] + sigma[i] * tan(y_uniform[i]);
        }
        return result;
    }
    
    vector getReparameterizedCauchyUpperBound_VVR_lp(vector y_uniform, vector mu, real sigma, real upperbound)
    {
        int length = num_elements(y_uniform);
        vector[length] result;
        real upper_transformed;
        if (num_elements(mu) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            upper_transformed = atan((upperbound - mu[i]) / sigma);
            sampleUniform_RRR_lp(y_uniform[i], -pi()/2, upper_transformed);
            result[i] = mu[i] + sigma * tan(y_uniform[i]);
        }
        return result;
    }
    
    vector getReparameterizedCauchyUpperBound_VVA_lp(vector y_uniform, vector mu, array[] real sigma, real upperbound)
    {
        int length = num_elements(y_uniform);
        vector[length] result;
        real upper_transformed;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            upper_transformed = atan((upperbound - mu[i]) / sigma[i]);
            sampleUniform_RRR_lp(y_uniform[i], -pi()/2, upper_transformed);
            result[i] = mu[i] + sigma[i] * tan(y_uniform[i]);
        }
        return result;
    }
    
    vector getReparameterizedCauchyUpperBound_VVV_lp(vector y_uniform, vector mu, vector sigma, real upperbound)
    {
        int length = num_elements(y_uniform);
        vector[length] result;
        real upper_transformed;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            upper_transformed = atan((upperbound - mu[i]) / sigma[i]);
            sampleUniform_RRR_lp(y_uniform[i], -pi()/2, upper_transformed);
            result[i] = mu[i] + sigma[i] * tan(y_uniform[i]);
        }
        return result;
    }
    
    // With range (lower and upper) bounds

    real getReparameterizedCauchyRangeBound_RRR_lp(real y_uniform, real mu, real sigma, real lowerbound, real upperbound)
    {
        real result;
        real lower_transformed;
        real upper_transformed;
        lower_transformed = atan((lowerbound - mu) / sigma);
        upper_transformed = atan((upperbound - mu) / sigma);
        sampleUniform_RRR_lp(y_uniform, lower_transformed, upper_transformed);
        result = mu + sigma * tan(y_uniform);
        return result;
    }
    
    array[] real getReparameterizedCauchyRangeBound_ARR_lp(array[] real y_uniform, real mu, real sigma, real lowerbound, real upperbound)
    {
        int length = num_elements(y_uniform);
        array[length] real result;
        real lower_transformed;
        real upper_transformed;
        for (i in 1:length) {
            lower_transformed = atan((lowerbound - mu) / sigma);
            upper_transformed = atan((upperbound - mu) / sigma);
            sampleUniform_RRR_lp(y_uniform[i], lower_transformed, upper_transformed);
            result[i] = mu + sigma * tan(y_uniform[i]);
        }
        return result;
    }
    
    array[] real getReparameterizedCauchyRangeBound_ARA_lp(array[] real y_uniform, real mu, array[] real sigma, real lowerbound, real upperbound)
    {
        int length = num_elements(y_uniform);
        array[length] real result;
        real lower_transformed;
        real upper_transformed;
        if (num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = atan((lowerbound - mu) / sigma[i]);
            upper_transformed = atan((upperbound - mu) / sigma[i]);
            sampleUniform_RRR_lp(y_uniform[i], lower_transformed, upper_transformed);
            result[i] = mu + sigma[i] * tan(y_uniform[i]);
        }
        return result;
    }
    
    array[] real getReparameterizedCauchyRangeBound_ARV_lp(array[] real y_uniform, real mu, vector sigma, real lowerbound, real upperbound)
    {
        int length = num_elements(y_uniform);
        array[length] real result;
        real lower_transformed;
        real upper_transformed;
        if (num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = atan((lowerbound - mu) / sigma[i]);
            upper_transformed = atan((upperbound - mu) / sigma[i]);
            sampleUniform_RRR_lp(y_uniform[i], lower_transformed, upper_transformed);
            result[i] = mu + sigma[i] * tan(y_uniform[i]);
        }
        return result;
    }
    
    array[] real getReparameterizedCauchyRangeBound_AAR_lp(array[] real y_uniform, array[] real mu, real sigma, real lowerbound, real upperbound)
    {
        int length = num_elements(y_uniform);
        array[length] real result;
        real lower_transformed;
        real upper_transformed;
        if (num_elements(mu) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = atan((lowerbound - mu[i]) / sigma);
            upper_transformed = atan((upperbound - mu[i]) / sigma);
            sampleUniform_RRR_lp(y_uniform[i], lower_transformed, upper_transformed);
            result[i] = mu[i] + sigma * tan(y_uniform[i]);
        }
        return result;
    }
    
    array[] real getReparameterizedCauchyRangeBound_AAA_lp(array[] real y_uniform, array[] real mu, array[] real sigma, real lowerbound, real upperbound)
    {
        int length = num_elements(y_uniform);
        array[length] real result;
        real lower_transformed;
        real upper_transformed;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = atan((lowerbound - mu[i]) / sigma[i]);
            upper_transformed = atan((upperbound - mu[i]) / sigma[i]);
            sampleUniform_RRR_lp(y_uniform[i], lower_transformed, upper_transformed);
            result[i] = mu[i] + sigma[i] * tan(y_uniform[i]);
        }
        return result;
    }
    
    array[] real getReparameterizedCauchyRangeBound_AAV_lp(array[] real y_uniform, array[] real mu, vector sigma, real lowerbound, real upperbound)
    {
        int length = num_elements(y_uniform);
        array[length] real result;
        real lower_transformed;
        real upper_transformed;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = atan((lowerbound - mu[i]) / sigma[i]);
            upper_transformed = atan((upperbound - mu[i]) / sigma[i]);
            sampleUniform_RRR_lp(y_uniform[i], lower_transformed, upper_transformed);
            result[i] = mu[i] + sigma[i] * tan(y_uniform[i]);
        }
        return result;
    }
    
    array[] real getReparameterizedCauchyRangeBound_AVR_lp(array[] real y_uniform, vector mu, real sigma, real lowerbound, real upperbound)
    {
        int length = num_elements(y_uniform);
        array[length] real result;
        real lower_transformed;
        real upper_transformed;
        if (num_elements(mu) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = atan((lowerbound - mu[i]) / sigma);
            upper_transformed = atan((upperbound - mu[i]) / sigma);
            sampleUniform_RRR_lp(y_uniform[i], lower_transformed, upper_transformed);
            result[i] = mu[i] + sigma * tan(y_uniform[i]);
        }
        return result;
    }
    
    array[] real getReparameterizedCauchyRangeBound_AVA_lp(array[] real y_uniform, vector mu, array[] real sigma, real lowerbound, real upperbound)
    {
        int length = num_elements(y_uniform);
        array[length] real result;
        real lower_transformed;
        real upper_transformed;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = atan((lowerbound - mu[i]) / sigma[i]);
            upper_transformed = atan((upperbound - mu[i]) / sigma[i]);
            sampleUniform_RRR_lp(y_uniform[i], lower_transformed, upper_transformed);
            result[i] = mu[i] + sigma[i] * tan(y_uniform[i]);
        }
        return result;
    }
    
    array[] real getReparameterizedCauchyRangeBound_AVV_lp(array[] real y_uniform, vector mu, vector sigma, real lowerbound, real upperbound)
    {
        int length = num_elements(y_uniform);
        array[length] real result;
        real lower_transformed;
        real upper_transformed;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = atan((lowerbound - mu[i]) / sigma[i]);
            upper_transformed = atan((upperbound - mu[i]) / sigma[i]);
            sampleUniform_RRR_lp(y_uniform[i], lower_transformed, upper_transformed);
            result[i] = mu[i] + sigma[i] * tan(y_uniform[i]);
        }
        return result;
    }
    
    array[,] real getReparameterizedCauchyRangeBound_2RR_lp(array[,] real y_uniform, real mu, real sigma, real lowerbound, real upperbound)
    {
        array[2] int dimensions = dims(y_uniform);
        array[dimensions[1], dimensions[2]] real result;
        real lower_transformed;
        real upper_transformed;
        lower_transformed = atan((lowerbound - mu) / sigma);
        upper_transformed = atan((upperbound - mu) / sigma);
        for (i in 1:dimensions[1]) {
            for (j in 1:dimensions[2]) {
                sampleUniform_RRR_lp(y_uniform[i, j], lower_transformed, upper_transformed);
                result[i, j] = mu + sigma * tan(y_uniform[i, j]);
            }
        }
        return result;
    }
    
    array[,,] real getReparameterizedCauchyRangeBound_3RR_lp(array[,,] real y_uniform, real mu, real sigma, real lowerbound, real upperbound)
    {
        array[3] int dimensions = dims(y_uniform);
        array[dimensions[1], dimensions[2], dimensions[3]] real result;
        real lower_transformed;
        real upper_transformed;
        lower_transformed = atan((lowerbound - mu) / sigma);
        upper_transformed = atan((upperbound - mu) / sigma);
        for (i in 1:dimensions[1]) {
            for (j in 1:dimensions[2]) {
                for (k in 1:dimensions[3]) {
                    sampleUniform_RRR_lp(y_uniform[i, j, k], lower_transformed, upper_transformed);
                    result[i, j, k] = mu + sigma * tan(y_uniform[i, j, k]);
                }
            }
        }
        return result;
    }
    
    vector getReparameterizedCauchyRangeBound_VRR_lp(vector y_uniform, real mu, real sigma, real lowerbound, real upperbound)
    {
        int length = num_elements(y_uniform);
        vector[length] result;
        real lower_transformed;
        real upper_transformed;
        for (i in 1:length) {
            lower_transformed = atan((lowerbound - mu) / sigma);
            upper_transformed = atan((upperbound - mu) / sigma);
            sampleUniform_RRR_lp(y_uniform[i], lower_transformed, upper_transformed);
            result[i] = mu + sigma * tan(y_uniform[i]);
        }
        return result;
    }
    
    vector getReparameterizedCauchyRangeBound_VRA_lp(vector y_uniform, real mu, array[] real sigma, real lowerbound, real upperbound)
    {
        int length = num_elements(y_uniform);
        vector[length] result;
        real lower_transformed;
        real upper_transformed;
        if (num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = atan((lowerbound - mu) / sigma[i]);
            upper_transformed = atan((upperbound - mu) / sigma[i]);
            sampleUniform_RRR_lp(y_uniform[i], lower_transformed, upper_transformed);
            result[i] = mu + sigma[i] * tan(y_uniform[i]);
        }
        return result;
    }
    
    vector getReparameterizedCauchyRangeBound_VRV_lp(vector y_uniform, real mu, vector sigma, real lowerbound, real upperbound)
    {
        int length = num_elements(y_uniform);
        vector[length] result;
        real lower_transformed;
        real upper_transformed;
        if (num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = atan((lowerbound - mu) / sigma[i]);
            upper_transformed = atan((upperbound - mu) / sigma[i]);
            sampleUniform_RRR_lp(y_uniform[i], lower_transformed, upper_transformed);
            result[i] = mu + sigma[i] * tan(y_uniform[i]);
        }
        return result;
    }
    
    vector getReparameterizedCauchyRangeBound_VAR_lp(vector y_uniform, array[] real mu, real sigma, real lowerbound, real upperbound)
    {
        int length = num_elements(y_uniform);
        vector[length] result;
        real lower_transformed;
        real upper_transformed;
        if (num_elements(mu) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = atan((lowerbound - mu[i]) / sigma);
            upper_transformed = atan((upperbound - mu[i]) / sigma);
            sampleUniform_RRR_lp(y_uniform[i], lower_transformed, upper_transformed);
            result[i] = mu[i] + sigma * tan(y_uniform[i]);
        }
        return result;
    }
    
    vector getReparameterizedCauchyRangeBound_VAA_lp(vector y_uniform, array[] real mu, array[] real sigma, real lowerbound, real upperbound)
    {
        int length = num_elements(y_uniform);
        vector[length] result;
        real lower_transformed;
        real upper_transformed;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = atan((lowerbound - mu[i]) / sigma[i]);
            upper_transformed = atan((upperbound - mu[i]) / sigma[i]);
            sampleUniform_RRR_lp(y_uniform[i], lower_transformed, upper_transformed);
            result[i] = mu[i] + sigma[i] * tan(y_uniform[i]);
        }
        return result;
    }
    
    vector getReparameterizedCauchyRangeBound_VAV_lp(vector y_uniform, array[] real mu, vector sigma, real lowerbound, real upperbound)
    {
        int length = num_elements(y_uniform);
        vector[length] result;
        real lower_transformed;
        real upper_transformed;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = atan((lowerbound - mu[i]) / sigma[i]);
            upper_transformed = atan((upperbound - mu[i]) / sigma[i]);
            sampleUniform_RRR_lp(y_uniform[i], lower_transformed, upper_transformed);
            result[i] = mu[i] + sigma[i] * tan(y_uniform[i]);
        }
        return result;
    }
    
    vector getReparameterizedCauchyRangeBound_VVR_lp(vector y_uniform, vector mu, real sigma, real lowerbound, real upperbound)
    {
        int length = num_elements(y_uniform);
        vector[length] result;
        real lower_transformed;
        real upper_transformed;
        if (num_elements(mu) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = atan((lowerbound - mu[i]) / sigma);
            upper_transformed = atan((upperbound - mu[i]) / sigma);
            sampleUniform_RRR_lp(y_uniform[i], lower_transformed, upper_transformed);
            result[i] = mu[i] + sigma * tan(y_uniform[i]);
        }
        return result;
    }
    
    vector getReparameterizedCauchyRangeBound_VVA_lp(vector y_uniform, vector mu, array[] real sigma, real lowerbound, real upperbound)
    {
        int length = num_elements(y_uniform);
        vector[length] result;
        real lower_transformed;
        real upper_transformed;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = atan((lowerbound - mu[i]) / sigma[i]);
            upper_transformed = atan((upperbound - mu[i]) / sigma[i]);
            sampleUniform_RRR_lp(y_uniform[i], lower_transformed, upper_transformed);
            result[i] = mu[i] + sigma[i] * tan(y_uniform[i]);
        }
        return result;
    }
    
    vector getReparameterizedCauchyRangeBound_VVV_lp(vector y_uniform, vector mu, vector sigma, real lowerbound, real upperbound)
    {
        int length = num_elements(y_uniform);
        vector[length] result;
        real lower_transformed;
        real upper_transformed;
        if (num_elements(mu) != length || num_elements(sigma) != length) {
            reject("Incompatible arguments");
        }
        for (i in 1:length) {
            lower_transformed = atan((lowerbound - mu[i]) / sigma[i]);
            upper_transformed = atan((upperbound - mu[i]) / sigma[i]);
            sampleUniform_RRR_lp(y_uniform[i], lower_transformed, upper_transformed);
            result[i] = mu[i] + sigma[i] * tan(y_uniform[i]);
        }
        return result;
    }
    
    // ------------------------------------------------------------------------
    // ANOVA-type designs: DEPRECATED APPROACH
    // ------------------------------------------------------------------------
    // ... rather than coding intercept + main effects + interactions (etc.),
    // as here, it's probably best to code individual cells. That makes
    // distributions more sensible (and predictable/easily calculable).

    int interactionIndex2Way(int first_index, int first_max,
                             int second_index, int second_max)
    {
        /*
            Because Stan doesn't support sampling into matrix, we need to
                         ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
            convert matrix-like concepts to vectors. Specifically, it doesn't
            support either
                matrix[A, B] m;
                m ~ normal(0, 0.5);  // error: "no matches for matrix ~ normal(int, real)"
            or
                real a[A, B];
                a ~ normal(0, 0.5);  // error: "no matches for array[,] real ~ normal(int, real)"

            And note that a vectorized sampling statement is strongly preferred
            (for performance reasons) over iterating through a matrix:
                https://groups.google.com/forum/#!topic/stan-users/4gv3fNCqSNk
                    "Do not loop over sampling statements when a vectorized
                    sampling statement is possible"

            So we use a vector of size A*B, and this index lookup function.
            Parameters:
            - first_index is from 1 to first_max
            - second_index is from 1 to second_max
            - We want a consecutive index from 1 to (first_max * second_max)

            In the output, the FIRST will cycle LEAST rapidly, and the
            LAST will cycle MOST rapidly.
        */
        return (
            (first_index - 1) * first_max +     // slow cycling
            second_index                        // fast cycling
        );
    }

    vector setLastForZeroSum(vector params)
    {
        /*
            Makes a vector of parameters sum to zero, by setting the last
            element to the negative sum of the others.
            Used for ANOVA-style effects; e.g. if you have a grand mean, you
            might specify the effects of a three-level factor A as A1, A2, A3;
            then A1 + A2 + A3 must be zero, so A1 and A2 are free parameters
            that are drawn from an appropriate distribution, and then A3 is
            fully constrainted to be -(A1 + A2).

            Because we can't modify the input parameters, we make a new copy.

            Returns a vector of the SAME LENGTH as the original.
            (The last element of the incoming vector is ignored.)
        */
        int length = num_elements(params);
        vector[length] newparams;
        real total = 0.0;
        real value;
        for (i in 1:length - 1) {
            value = params[i];
            newparams[i] = value;
            total = total + value;
        }
        newparams[length] = -total;
        return newparams;
    }

    vector appendElementForZeroSum(vector params)
    {
        /*
            As for setLastForZeroSum(), but uses all the information in the
            incoming vector, and returns a vector that's one element longer.
        */
        int initial_length = num_elements(params);
        int new_length = initial_length + 1;
        vector[new_length] newparams;
        real total = 0.0;
        real value;
        for (i in 1:initial_length) {
            value = params[i];
            newparams[i] = value;
            total = total + value;
        }
        newparams[new_length] = -total;
        return newparams;
    }