# chol()

Jonathan Rougier J.C.Rougier@durham.ac.uk
Thu, 8 Apr 1999 10:26:23 +0100 (BST)

```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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._

```