R/lpmvnorm.R

In mvtnorm: Multivariate Normal and t Distributions

# R Header

###    Copyright (C) 2022- Torsten Hothorn
###
###    This file is part of the 'mvtnorm' R add-on package.
###
###    'mvtnorm' is free software: you can redistribute it and/or modify
###    it under the terms of the GNU General Public License as published by
###    the Free Software Foundation, version 2.
###
###    'mvtnorm' is distributed in the hope that it will be useful,
###    but WITHOUT ANY WARRANTY; without even the implied warranty of
###    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
###    GNU General Public License for more details.
###
###    You should have received a copy of the GNU General Public License
###    along with 'mvtnorm'.  If not, see <http://www.gnu.org/licenses/>.
###
###
###    DO NOT EDIT THIS FILE
###
###    Edit 'lmvnorm_src.w' and run 'nuweb -r lmvnorm_src.w'

# lpmvnorm

lpmvnorm <- function(lower, upper, mean = 0, center = NULL, chol, invchol, 
                     logLik = TRUE, M = NULL, w = NULL, seed = NULL, 
                     tol = .Machine$double.eps, fast = FALSE) {

    # init random seed, reset on exit
    
    ### from stats:::simulate.lm
    if (!exists(".Random.seed", envir = .GlobalEnv, inherits = FALSE)) 
        runif(1)
    if (is.null(seed)) 
        RNGstate <- get(".Random.seed", envir = .GlobalEnv)
    else {
        R.seed <- get(".Random.seed", envir = .GlobalEnv)
        set.seed(seed)
        RNGstate <- structure(seed, kind = as.list(RNGkind()))
        on.exit(assign(".Random.seed", R.seed, envir = .GlobalEnv))
    }
    
    # Cholesky of precision
    
    stopifnot(xor(missing(chol), missing(invchol)))
    if (missing(chol)) chol <- solve(invchol)
    
    # input checks
    
    if (!is.matrix(lower)) lower <- matrix(lower, ncol = 1)
    if (!is.matrix(upper)) upper <- matrix(upper, ncol = 1)
    stopifnot(isTRUE(all.equal(dim(lower), dim(upper))))

    stopifnot(inherits(chol, "ltMatrices"))
    byrow_orig <- attr(chol, "byrow")
    chol <- ltMatrices(chol, byrow = TRUE)
    d <- dim(chol)
    ### allow single matrix C
    N <- ifelse(d[1L] == 1, ncol(lower), d[1L])
    J <- d[2L]

    stopifnot(nrow(lower) == J && ncol(lower) == N)
    stopifnot(nrow(upper) == J && ncol(upper) == N)
    if (is.matrix(mean))
        stopifnot(nrow(mean) == J && ncol(mean) == N)

    lower <- lower - mean
    upper <- upper - mean

    if (!is.null(center)) {
        if (!is.matrix(center)) center <- matrix(center, ncol = 1)
        stopifnot(nrow(center) == J && ncol(center == N))
    }
    
    # standardise
    
    if (attr(chol, "diag")) {
        ### diagonals returns J x N and lower/upper are J x N, so
        ### elementwise standardisation is simple
        dchol <- diagonals(chol)
        ### zero diagonals not allowed
        stopifnot(all(abs(dchol) > (.Machine$double.eps)))
        ac <- lower / c(dchol)
        bc <- upper / c(dchol)
        C <- Dchol(chol, D = 1 / dchol)
        uC <- unclass(C)
        if (J > 1) ### else: univariate problem; C is no longer used
           uC <- Lower_tri(C)
        } else {
            ac <- lower
            bc <- upper
            uC <- Lower_tri(chol)
        }
    
    # check and / or set integration weights
    
    if (!is.null(w) && J > 1) {
        stopifnot(is.matrix(w))
        stopifnot(nrow(w) == J - 1)
        if (is.null(M))
            M <- ncol(w)
        stopifnot(ncol(w) %in% c(M, M * N))
        storage.mode(w) <- "double"
    } else {
        if (J > 1) {
            if (is.null(M)) stop("either w or M must be specified")
        } else {
            M <- 1L
        }
    }
    

    ret <- .Call(mvtnorm_R_lpmvnorm, ac, bc, uC, as.double(center), 
                 as.integer(N), as.integer(J), w, as.integer(M), as.double(tol), 
                 as.logical(logLik), as.logical(fast));
    return(ret)
}

# slpmvnorm

slpmvnorm <- function(lower, upper, mean = 0, center = NULL, chol, invchol, logLik = TRUE, M = NULL, 
                    w = NULL, seed = NULL, tol = .Machine$double.eps, fast = FALSE) {

    # init random seed, reset on exit
    
    ### from stats:::simulate.lm
    if (!exists(".Random.seed", envir = .GlobalEnv, inherits = FALSE)) 
        runif(1)
    if (is.null(seed)) 
        RNGstate <- get(".Random.seed", envir = .GlobalEnv)
    else {
        R.seed <- get(".Random.seed", envir = .GlobalEnv)
        set.seed(seed)
        RNGstate <- structure(seed, kind = as.list(RNGkind()))
        on.exit(assign(".Random.seed", R.seed, envir = .GlobalEnv))
    }
    
    # Cholesky of precision
    
    stopifnot(xor(missing(chol), missing(invchol)))
    if (missing(chol)) chol <- solve(invchol)
    
    # input checks
    
    if (!is.matrix(lower)) lower <- matrix(lower, ncol = 1)
    if (!is.matrix(upper)) upper <- matrix(upper, ncol = 1)
    stopifnot(isTRUE(all.equal(dim(lower), dim(upper))))

    stopifnot(inherits(chol, "ltMatrices"))
    byrow_orig <- attr(chol, "byrow")
    chol <- ltMatrices(chol, byrow = TRUE)
    d <- dim(chol)
    ### allow single matrix C
    N <- ifelse(d[1L] == 1, ncol(lower), d[1L])
    J <- d[2L]

    stopifnot(nrow(lower) == J && ncol(lower) == N)
    stopifnot(nrow(upper) == J && ncol(upper) == N)
    if (is.matrix(mean))
        stopifnot(nrow(mean) == J && ncol(mean) == N)

    lower <- lower - mean
    upper <- upper - mean

    if (!is.null(center)) {
        if (!is.matrix(center)) center <- matrix(center, ncol = 1)
        stopifnot(nrow(center) == J && ncol(center == N))
    }
    
    # standardise
    
    if (attr(chol, "diag")) {
        ### diagonals returns J x N and lower/upper are J x N, so
        ### elementwise standardisation is simple
        dchol <- diagonals(chol)
        ### zero diagonals not allowed
        stopifnot(all(abs(dchol) > (.Machine$double.eps)))
        ac <- lower / c(dchol)
        bc <- upper / c(dchol)
        C <- Dchol(chol, D = 1 / dchol)
        uC <- unclass(C)
        if (J > 1) ### else: univariate problem; C is no longer used
           uC <- Lower_tri(C)
        } else {
            ac <- lower
            bc <- upper
            uC <- Lower_tri(chol)
        }
    
    # check and / or set integration weights
    
    if (!is.null(w) && J > 1) {
        stopifnot(is.matrix(w))
        stopifnot(nrow(w) == J - 1)
        if (is.null(M))
            M <- ncol(w)
        stopifnot(ncol(w) %in% c(M, M * N))
        storage.mode(w) <- "double"
    } else {
        if (J > 1) {
            if (is.null(M)) stop("either w or M must be specified")
        } else {
            M <- 1L
        }
    }
    

    ret <- .Call(mvtnorm_R_slpmvnorm, ac, bc, uC, as.double(center), as.integer(N), 
                 as.integer(J), w, as.integer(M), as.double(tol), as.logical(fast));

    ll <- log(pmax(ret[1L,], tol)) - log(M)
    intsum <- ret[1L,]
    m <- matrix(intsum, nrow = nrow(ret) - 1, ncol = ncol(ret), byrow = TRUE)
    ret <- ret[-1L,,drop = FALSE] / m ### NOTE: division by zero MAY happen,
                                      ### catch outside

    # post differentiate mean score
    
    Jp <- J * (J + 1) / 2;
    smean <- - ret[Jp + 1:J, , drop = FALSE]
    if (attr(chol, "diag"))
        smean <- smean / c(dchol)
    
    # post differentiate lower score
    
    slower <- ret[Jp + J + 1:J, , drop = FALSE]
    if (attr(chol, "diag"))
        slower <- slower / c(dchol)
    
    # post differentiate upper score
    
    supper <- ret[Jp + 2 * J + 1:J, , drop = FALSE]
    if (attr(chol, "diag"))
        supper <- supper / c(dchol)
    

    ret <- ret[1:Jp, , drop = FALSE]

    # post differentiate chol score
    
    if (J == 1) {
        idx <- 1L
    } else {
        idx <- cumsum(c(1, 2:J))
    }
    if (attr(chol, "diag")) {
        ret <- ret / c(dchol[rep(1:J, 1:J),]) ### because 1 / dchol already there
        ret[idx,] <- -ret[idx,]
    }
    
    # post differentiate invchol score
    
    if (!missing(invchol)) {
        ret <- ltMatrices(ret, diag = TRUE, byrow = TRUE)
        ### this means vectrick(chol, ret, chol)
        ret <- - unclass(vectrick(chol, ret))
    }
    
    # post process score
    
    if (!attr(chol, "diag"))
        ### remove scores for constant diagonal elements
        ret[idx,] <- 0
    ret <- ltMatrices(ret, diag = TRUE, byrow = TRUE)
    

    ret <- ltMatrices(ret, byrow = byrow_orig)

    if (logLik) {
        ret <- list(logLik = ll, 
                    mean = smean, 
                    lower = slower,
                    upper = supper,
                    chol = ret)
        if (!missing(invchol)) names(ret)[names(ret) == "chol"] <- "invchol"
        return(ret)
    }
    
    return(ret)
}