-
Notifications
You must be signed in to change notification settings - Fork 2
/
my_dmvnorm.R
45 lines (45 loc) · 1.36 KB
/
my_dmvnorm.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
#'Copy of mvtnorm::dmvnorm for progiling and reference
#'
#'@param x matrix
#'@param mean vector
#'@param sigma covariance matrix
#'@param log logical. Default is \code{FALSE}
#'@param checkSymmetry logical. Default is \code{TRUE}
#'
#'@export
#'
my_dmvnorm <- function(x, mean = rep(0, p), sigma = diag(p), log = FALSE,
checkSymmetry = TRUE){
if (is.vector(x))
x <- matrix(x, ncol = length(x))
p <- ncol(x)
if (!missing(mean)) {
if (!is.null(dim(mean)))
dim(mean) <- NULL
if (length(mean) != p)
stop("x and mean have non-conforming size")
}
if (!missing(sigma)) {
if (p != ncol(sigma))
stop("x and sigma have non-conforming size")
if (checkSymmetry && !isSymmetric(sigma, tol = sqrt(.Machine$double.eps),
check.attributes = FALSE))
stop("sigma must be a symmetric matrix")
}
dec <- tryCatch(base::chol(sigma), error = function(e) e)
if (inherits(dec, "error")) {
x.is.mu <- colSums(t(x) != mean) == 0
logretval <- rep.int(-Inf, nrow(x))
logretval[x.is.mu] <- Inf
}
else {
tmp <- backsolve(dec, t(x) - mean, transpose = TRUE)
rss <- colSums(tmp^2)
logretval <- -sum(log(diag(dec))) - 0.5 * p * log(2 *
pi) - 0.5 * rss
}
names(logretval) <- rownames(x)
if (log)
logretval
else exp(logretval)
}