R/mnormt.R

In mnormt: The Multivariate Normal and t Distributions, and Their Truncated Versions

# R code of package 'mnormt' 
# Author: Adelchi Azzalini (University of Padua, Italy) 

dmnorm <- function(x, mean=rep(0,d), varcov, log=FALSE)
{
  d <- if(is.matrix(varcov)) ncol(varcov) else 1
  if(d==1) return(dnorm(x, mean, sqrt(varcov), log=log))
  x <- if (is.vector(x)) t(matrix(x)) else data.matrix(x)
  if(ncol(x) != d) stop("mismatch of dimensions of 'x' and 'varcov'")
  if(is.matrix(mean)) { if ((nrow(x) != nrow(mean)) || (ncol(mean) != d))
      stop("mismatch of dimensions of 'x' and 'mean'") }
  if(is.vector(mean)) mean <- outer(rep(1, nrow(x)), as.vector(matrix(mean,d)))
  X  <- t(x - mean)
  conc <- pd.solve(varcov, log.det=TRUE)
  Q <- colSums((conc %*% X)* X)
  log.det <- attr(conc, "log.det")
  logPDF <- as.vector(Q + d*logb(2*pi) + log.det)/(-2)
  if(log) logPDF else exp(logPDF)
}

rmnorm <- function(n=1, mean=rep(0,d), varcov, sqrt=NULL)
 {
  sqrt.varcov <- if(is.null(sqrt)) chol(varcov) else sqrt
  d <- if(is.matrix(sqrt.varcov)) ncol(sqrt.varcov) else 1
  mean <- outer(rep(1,n), as.vector(matrix(mean,d)))
  drop(mean + t(matrix(rnorm(n*d), d, n)) %*% sqrt.varcov)
 }


pmnorm <- function(x, mean=rep(0, d), varcov, ...) {
  d <- NCOL(varcov)
  x <- if (is.vector(x)) matrix(x, 1, d) else data.matrix(x)
  n <- nrow(x)
  if(is.vector(mean)) mean <- outer(rep(1, n), as.vector(matrix(mean,d)))
  if(d == 1) p <- as.vector(pnorm(x, mean, sqrt(varcov))) else {
    pv <- numeric(n)
    for (j in 1:n) 
      p <- pv[j] <- {
        if(d == 2)
          biv.nt.prob(Inf, lower=rep(-Inf, 2), upper=x[j,], mean[j,], varcov)   
        else if(d == 3) 
          ptriv.nt(Inf, x=x[j,], mean[j,], varcov)
            else     
             sadmvn(lower=rep(-Inf, d), upper=x[j,], mean[j,], varcov, ...) 
        }     
    if(n > 1) p <- pv 
    }
  return(p)  
  }

sadmvn <- function(lower, upper, mean, varcov,  
                   maxpts=2000*d, abseps=1e-6, releps=0)
{
  if(any(lower > upper)) stop("lower>upper integration limits")
  if(any(lower == upper)) return(0)
  d <- as.integer(if(is.matrix(varcov)) ncol(varcov) else 1)
  varcov <- matrix(varcov, d, d)
  sd  <- sqrt(diag(varcov))
  rho <- cov2cor(varcov)
  lower <- as.double((lower-mean)/sd)
  upper <- as.double((upper-mean)/sd)
  if(d == 1) return(pnorm(upper) - pnorm(lower))
  if(d == 2) return(biv.nt.prob(Inf, lower, upper, rep(0,2), rho))
  infin <- rep(2,d)
  infin <- replace(infin, (upper == Inf) & (lower > -Inf), 1)
  infin <- replace(infin, (upper < Inf) & (lower == -Inf), 0)
  infin <- replace(infin, (upper == Inf) & (lower == -Inf), -1)
  infin <- as.integer(infin)
  if(any(infin == -1)) {
    if(all(infin == -1)) return(1)
    k <- which(infin != -1)
    d <- length(k)
    lower <- lower[k]
    upper <- upper[k]
    if(d == 1) return(pnorm(upper) - pnorm(lower))
    rho <- rho[k, k]
    infin <- infin[k]
    if(d == 2) return(biv.nt.prob(Inf, lower, upper, rep(0,2), rho))
    }
  if(d > 20) # 2022-09-21, overrides a "0 value & error msg" from sadmvn 
     return(NA)
  lower <- replace(lower, lower == -Inf, 0)
  upper <- replace(upper, upper == Inf, 0)
  correl <- as.double(rho[upper.tri(rho, diag=FALSE)])
  maxpts <- as.integer(maxpts)
  abseps <- as.double(abseps)
  releps <- as.double(releps)
  error  <- as.double(0)
  value  <- as.double(0)
  inform <- as.integer(0)
  result <- .Fortran("sadmvn", d, lower, upper, infin, correl, maxpts,
               abseps, releps, error, value, inform, PACKAGE="mnormt")
  prob <- result[[10]]
  attr(prob,"error")  <- result[[9]]
  attr(prob,"status") <- switch(1 + result[[11]], 
                "normal completion", "accuracy non achieved", "oversize")
  return(prob)
}

#----

dmt <- function (x, mean=rep(0,d), S, df = Inf, log = FALSE) 
{
  if (df == Inf)  return(dmnorm(x, mean, S, log = log))
  d <- if(is.matrix(S)) ncol(S) else 1
  if (d==1) {
    y <- dt((x-mean)/sqrt(S), df=df, log=log)
    if(log) y <- (y - 0.5*logb(S)) else y <- y/sqrt(S)
    return(y)
    }
  x <- if (is.vector(x)) t(matrix(x)) else data.matrix(x)
  if (ncol(x) != d) stop("mismatch of dimensions of 'x' and 'varcov'")
  if (is.matrix(mean)) {if ((nrow(x) != nrow(mean)) || (ncol(mean) != d))
      stop("mismatch of dimensions of 'x' and 'mean'") }
  if(is.vector(mean)) mean <- outer(rep(1, nrow(x)), as.vector(matrix(mean,d)))
  X  <- t(x - mean)
  S.inv <- pd.solve(S, log.det=TRUE)
  Q <- colSums((S.inv %*% X) * X)
  logDet <- attr(S.inv, "log.det")
  logPDF <- (lgamma((df + d)/2) - 0.5 * (d * logb(pi * df) + logDet)
             - lgamma(df/2) - 0.5 * (df + d) * logb(1 + Q/df))
  if(log) logPDF else exp(logPDF)
}


rmt <- function(n=1, mean=rep(0,d), S, df=Inf, sqrt=NULL)
{ 
  sqrt.S <- if(is.null(sqrt)) chol(S) else sqrt
  d <- if(is.matrix(sqrt.S)) ncol(sqrt.S) else 1
  # rs <<- .Random.seed 
  v <- if(df==Inf) 1 else rchisq(n, df)/df
  # .Random.seed <<- rs 
  z <- rmnorm(n, rep(0, d), sqrt=sqrt.S)
  mean <- outer(rep(1, n), as.vector(matrix(mean, d)))
  drop(mean + z/sqrt(v))
}
 


pmt <- function(x, mean=rep(0, d), S, df=Inf, ...){
  d <- NCOL(S)
  x <- if(is.vector(x)) matrix(x, 1, d) else data.matrix(x)
  n <- nrow(x)
  if(is.vector(mean)) mean <- outer(rep(1, n), as.vector(matrix(mean,d)))
  if(d == 1) p <- as.vector(pt((x-mean)/sqrt(S), df=df)) else {
    pv <- numeric(n)
    for (j in 1:n) 
       p <- pv[j] <- {
         if(d == 2)
           biv.nt.prob(df, lower=rep(-Inf, 2), upper=x[j,], mean[j,], S)  
          else if(d == 3)
            ptriv.nt(df, x=x[j,], mean=mean[j,], S)      
              else 
                sadmvt(df, lower=rep(-Inf, d), upper=x[j,], mean[j,], S, ...)
         }
     if(n > 1) p <- pv    
     }
  return(p)  
  }
   


sadmvt <- function(df, lower, upper, mean, S, 
                   maxpts=2000*d, abseps=1e-6, releps=0)
{
  if(df == Inf) return(sadmvn(lower, upper, mean, S, maxpts, abseps, releps))
  if(any(lower > upper)) stop("lower>upper integration limits")
  if(any(lower == upper)) return(0)
  if(round(df) != df) warning("non integer df is rounded to integer") 
  df <- as.integer(round(df))
  d  <- as.integer(if(is.matrix(S)) ncol(S) else 1)
  S  <- matrix(S, d, d)
  sd  <- sqrt(diag(S))
  rho <- cov2cor(S)
  lower <- as.double((lower-mean)/sd)
  upper <- as.double((upper-mean)/sd)
  if(d == 1) return(pt(upper, df) - pt(lower, df))
  if(d == 2) return(biv.nt.prob(df, lower, upper, rep(0,2), rho))
  infin <- rep(2,d)
  infin <- replace(infin, (upper == Inf) & (lower > -Inf), 1)
  infin <- replace(infin, (upper < Inf) & (lower == -Inf), 0)
  infin <- replace(infin, (upper == Inf) & (lower == -Inf), -1)
  infin <- as.integer(infin)
  if(any(infin == -1)) {
    if(all(infin == -1)) return(1)
    k <- which(infin != -1)
    d <- length(k)
    lower <- lower[k]
    upper <- upper[k]
    if(d == 1) return(pt(upper, df=df) - pt(lower, df=df))
    rho <- rho[k, k]
    infin <- infin[k]
    if(d == 2) return(biv.nt.prob(df, lower, upper, rep(0,2), rho))
    }
  if(d > 20) # 2022-09-21, overrides a "0 value & error msg" from sadmvn 
     return(NA)  
  lower <- replace(lower, lower == -Inf, 0)
  upper <- replace(upper, upper == Inf, 0)
  correl <- rho[upper.tri(rho, diag=FALSE)]
  maxpts <- as.integer(maxpts)
  abseps <- as.double(abseps)
  releps <- as.double(releps)
  error  <- as.double(0)
  value  <- as.double(0)
  inform <- as.integer(0)
  result <- .Fortran("sadmvt", d, df, lower, upper, infin, correl, maxpts,
                   abseps, releps, error, value, inform, PACKAGE="mnormt")
  prob <- result[[11]]
  attr(prob,"error")  <- result[[10]]
  attr(prob,"status") <- switch(1+result[[12]], 
                "normal completion", "accuracy non achieved", "oversize")
  return(prob)
}


biv.nt.prob <- function(df, lower, upper, mean, S){
  if(any(dim(S) != c(2,2))) stop("dimensions mismatch")
  if(length(mean) != 2) stop("dimensions mismatch") 
  if(round(df) != df) warning("non integer df is rounded to integer") 
  nu <- if(df < Inf) as.integer(round(df)) else 0
  sd <- sqrt(diag(S))
  rho <- cov2cor(S)[1,2]
  lower <- as.double((lower-mean)/sd)
  upper <- as.double((upper-mean)/sd)
  if(any(lower > upper)) stop("lower>upper integration limits")
  if(any(lower == upper)) return(0)
  infin <- c(2,2)
  infin <- replace(infin, (upper == Inf) & (lower > -Inf), 1)
  infin <- replace(infin, (upper < Inf) & (lower == -Inf), 0)
  infin <- replace(infin, (upper == Inf) & (lower == -Inf), -1)
  infin <- as.integer(infin)
  if(any(infin == -1)) {
    if(all(infin == -1)) return(1)
    k <- which(infin != -1)
    return(pt(upper[k], df=df) - pt(lower[k], df=df))
    }
  lower <- replace(lower, lower == -Inf, 0)
  upper <- replace(upper, upper == Inf, 0)
  rho   <- as.double(rho)
  prob  <- as.double(0)
  a <- .Fortran("smvbvt", prob, nu, lower, upper, infin, rho, PACKAGE="mnormt")
  return(a[[1]])
  } 

ptriv.nt <- function(df, x, mean, S){
  if(any(dim(S) != c(3,3))) stop("dimensions mismatch")
  if(length(mean) != 3) stop("dimensions mismatch") 
  if(round(df) != df) warning("non integer df is rounded to integer") 
  nu <- if(df < Inf) as.integer(round(df)) else 0
  if(any(x == -Inf)) return(0)
  ok <- !is.infinite(x)
  p <- if(sum(ok) == 1) 
         pt(df, (x[ok]-mean[ok])/sqrt(S[ok,ok])) 
    else if(sum(ok) == 2) 
         biv.nt.prob(nu, rep(2 -Inf), x[ok], mean[ok], S[ok,ok])
    else {     
      sd <- sqrt(diag(S))
      h <- as.double((x-mean)/sd)
      cor <- cov2cor(S)
      rho  <- as.double(c(cor[2,1], cor[3,1], cor[2,3]))
      prob <- as.double(0)
      epsi <- as.double(1e-14)
      a <- .Fortran("stvtl", prob, nu, h, rho, epsi, PACKAGE="mnormt")
      p <- a[[1]]
      }
    return(p) 
  }  

pd.solve <- function(x, silent=FALSE, log.det=FALSE)
{
  if(is.null(x)) return(NULL)
  if(any(is.na(x)))
    {if(silent) return (NULL) else stop("NA's in x") } 
  if(length(x) == 1) x <- as.matrix(x)
  if(!(inherits(x, "matrix"))) 
    {if(silent) return(NULL) else stop("x is not a matrix")}
  if(max(abs(x - t(x))) > .Machine$double.eps) 
    {if(silent) return (NULL) else stop("x appears to be not symmetric") } 
  x <- (x + t(x))/2
  u <- try(chol(x, pivot = FALSE), silent = silent)
  if(inherits(u, "try-error")) {
     if(silent) return(NULL) else
       stop("x appears to be not positive definite") }
  inv <- chol2inv(u)
  if(log.det) attr(inv, "log.det") <- 2 * sum(log(diag(u)))
  dimnames(inv) <- rev(dimnames(x))
  return(inv)
}

.onLoad <- function(library, pkg)
{ 
   library.dynam("mnormt", pkg, library)
   invisible()
}