The R Project SVN R-packages

#include "Wishart.h"

/**
 * Simulate the Cholesky factor of a standardized Wishart variate with
 * dimension p and df degrees of freedom.
 *
 * @param df degrees of freedom
 * @param p dimension of the Wishart distribution
 * @param upper if 0 the result is lower triangular, otherwise upper
                triangular
 * @param ans array of size p * p to hold the result
 *
 * @return ans
 */
double attr_hidden
*std_rWishart_factor(double df, int p, int upper, double ans[])
{
    int i, j, pp1 = p + 1;

    if (df < (double) p || p <= 0)
    error("inconsistent degrees of freedom and dimension");
    for (j = 0; j < p; j++) {   /* jth column */
    ans[j * pp1] = sqrt(rchisq(df - (double) j));
    for (i = 0; i < j; i++) {
        int uind = i + j * p, /* upper triangle index */
        lind = j + i * p; /* lower triangle index */
        ans[(upper ? uind : lind)] = norm_rand();
        ans[(upper ? lind : uind)] = 0;
    }
    }
    return ans;
}

/**
 * Simulate a sample of random matrices from a Wishart distribution
 *
 * @param ns Number of samples to generate
 * @param dfp Degrees of freedom
 * @param scal Positive-definite scale matrix
 *
 * @return
 */
SEXP
lme4_rWishart(SEXP ns, SEXP dfp, SEXP scal)
{
    SEXP ans;
    int *dims = INTEGER(getAttrib(scal, R_DimSymbol)), j,
    n = asInteger(ns), psqr;
    double *scCp, *ansp, *tmp, df = asReal(dfp), one = 1, zero = 0;

    if (!isMatrix(scal) || !isReal(scal) || dims[0] != dims[1])
    error("scal must be a square, real matrix");
    if (n <= 0) n = 1;
    psqr = dims[0] * dims[0];
    tmp = Calloc(psqr, double);
    AZERO(tmp, psqr);
    scCp = Memcpy(Calloc(psqr, double), REAL(scal), psqr);
    F77_CALL(dpotrf)("U", &(dims[0]), scCp, &(dims[0]), &j);
    if (j)
    error("scal matrix is not positive-definite");
    PROTECT(ans = alloc3Darray(REALSXP, dims[0], dims[0], n));
    ansp = REAL(ans);
    GetRNGstate();
    for (j = 0; j < n; j++) {
    double *ansj = ansp + j * psqr;
    std_rWishart_factor(df, dims[0], 1, tmp);
    F77_CALL(dtrmm)("R", "U", "N", "N", dims, dims,
            &one, scCp, dims, tmp, dims);
    F77_CALL(dsyrk)("U", "T", &(dims[1]), &(dims[1]),
            &one, tmp, &(dims[1]),
            &zero, ansj, &(dims[1]));
    internal_symmetrize(ansj, dims[0]);
    }

    PutRNGstate();
    Free(scCp); Free(tmp);
    UNPROTECT(1);
    return ans;
}