chol()

[Jonathan Rougier]

May I suggest the following modification to chol(), allowing calls on
variance matrices with zeros in the diagonal.  This makes the generation
of multivariate gaussian vectors simpler, and consistent with rnorm()
which can take a zero standard deviation argument. 

"Chol" <-
function (x, tol = sqrt(.Machine$double.eps)) 
{
    if (!is.numeric(x)) 
        stop("non-numeric argument to chol")
    x <- as.matrix(x)
    if (length(n <- unique(dim(x))) > 1) 
        stop("non-square matrix in chol")
    nzero <- diag(x) > tol
    if (!all(nzero)) {
        cx <- Recall(x[nzero, nzero])
        x[] <- 0
        x[nzero, nzero] <- cx
        x
    }
    else {
        if (!is.double(x)) 
            storage.mode(x) <- "double"
        z <- .Fortran("chol", x = x, n, n, v = matrix(0, nr = n, 
            nc = n), info = integer(1), DUP = FALSE)
        if (z$info) 
            stop("singular matrix in chol")
        z$v
    }
}

I know that almost everyone has got one of these, but for completeness
here is my multivariate gaussian function which defaults to rnorm()

"rgauss" <-
function (n, mu, sig, ch = all(sig[lower.tri(sig)] == 0)) 
{
    if (missing(mu)) {
        if (missing(sig)) 
            sig <- 1
        sig <- as.matrix(sig)
        k <- unique(dim(sig))
        if (length(k) > 1) 
            stop("Non-square \"sig\" matrix")
        mu <- rep(0, k)
    }
    else {
        k <- length(mu)
        if (missing(sig)) 
            sig <- diag(rep(1, k))
        else {
            sig <- as.matrix(sig)
            if (any(k != dim(sig))) 
                stop("\"mu\" and \"sig\" inconsistent")
        }
    }
    # uses Chol() to handle zero variances
    if (!ch) 
        sig <- Chol(sig)
    rn <- matrix(rnorm(n * k), nrow = n)
    drop(sweep(rn %*% sig, 2, mu, "+"))
}

Regards, Jonathan.

Jonathan Rougier                       Science Laboratories
Department of Mathematical Sciences    South Road
University of Durham                   Durham DH1 3LE

"[B]egin upon the precept ... that the things we see are to be 
 weighed in the scale with what we know"  (Meredith, 1879, The Egoist)

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._