lpmvnorm: Multivariate Normal Log-likelihood and Score Functions

Description

Computes the log-likelihood (contributions) of multiple exact or interval-censored observations (or a mix thereof) from multivariate normal distributions and evaluates corresponding score functions.

Usage

lpmvnorm(lower, upper, mean = 0, center = NULL, chol, invchol, logLik = TRUE, 
         M = NULL, w = NULL, seed = NULL, tol = .Machine$double.eps, fast = FALSE)
slpmvnorm(lower, upper, mean = 0, center = NULL, chol, invchol, logLik = TRUE, 
          M = NULL, w = NULL, seed = NULL, tol = .Machine$double.eps, fast = FALSE)
ldmvnorm(obs, mean = 0, chol, invchol, logLik = TRUE) 
sldmvnorm(obs, mean = 0, chol, invchol, logLik = TRUE) 
ldpmvnorm(obs, lower, upper, mean = 0, chol, invchol, logLik = TRUE, ...) 
sldpmvnorm(obs, lower, upper, mean = 0, chol, invchol, logLik = TRUE, ...)

Value

The log-likelihood (logLik = TRUE) or the individual contributions to the log-likelihood. slpmvnorm, sldmvnorm, and sldpmvnorm return the score matrices and, optionally (logLik = TRUE), the individual log-likelihood contributions as well as scores for obs, lower, upper, and mean.

Arguments

lower: matrix of lower limits (one column for each observation, $J$ rows).
upper: matrix of upper limits (one column for each observation, $J$ rows).
obs: matrix of exact observations (one column for each observation, $J$ rows).
mean: matrix of means (one column for each observation, length is recycled to length of obs, lower and upper).
center: matrix of negative rescaled means (one column for each observation, length is recycled to length of lower and upper) as returned by cond_mvnorm(..., center = TRUE).
chol: Cholesky factors of covariance matrices as ltMatrices object, length is recylced to length of obs, lower and upper.
invchol: Cholesky factors of precision matrices as ltMatrices object, length is recylced to length of lower and upper. Either chol or invchol must be given.
logLik: logical, if TRUE, the log-likelihood is returned, otherwise the individual contributions to the sum are returned.
M: number of iterations, early stopping based on estimated errors is NOT implemented.
w: an optional matrix of weights with $J - 1$ rows. This allows to replace the default Monte-Carlo procedure (Genz, 1992) with a quasi-Monte-Carlo approach (Genz & Bretz, 2002). Note that the same weights for evaluating the multivariate normal probability are used for all observations when ncol(w) == M is specified. If ncol(w) == ncol(lower) * M, each likelihood contribution is evaluated on the corresponding sub-matrix. If w is NULL, different uniform numbers are drawn for each observation.
seed: an object specifying if and how the random number generator should be initialized, see simulate. Only applied when w is NULL.
tol: tolerance limit, values smaller than tol are interpreted as zero.
fast: logical, if TRUE, a faster but less accurate version of pnorm is used internally.
...: additional arguments to lpmvnorm.

Details

Evaluates the multivariate normal log-likelihood defined by means and chol over boxes defined by lower and upper or for exact observations obs.

Monte-Carlo (Genz, 1992, the default) and quasi-Monte-Carlo (Genz & Bretz, 2002) integration is implemented, the latter with weights obtained, for example, from packages qrng or randtoolbox. It is the responsibility of the user to ensure a meaningful lattice is used. In case of doubt, use plain Monte-Carlo (w = NULL) or pmvnorm.

slpmvnorm computes both the individual log-likelihood contributions and the corresponding score matrix (of dimension $J \times (J + 1) / 2 \times N$ ) if chol contains diagonal elements. Otherwise, the dimension is $J \times (J - 1) / 2 \times N$ . The scores for exact or mixed exact-interval observations are computed by sldmvnorm and sldpmvnorm, respectively.

More details can be found in the lmvnorm_src package vignette.

References

Genz, A. (1992). Numerical computation of multivariate normal probabilities. Journal of Computational and Graphical Statistics, 1, 141--150.

Genz, A. and Bretz, F. (2002), Methods for the computation of multivariate t-probabilities. Journal of Computational and Graphical Statistics, 11, 950--971.

Examples

Run this code


  ### five observations
  N <- 5L
  ### dimension
  J <- 4L

  ### lower and upper bounds, ie interval-censoring
  lwr <- matrix(-runif(N * J), nrow = J)
  upr <- matrix(runif(N * J), nrow = J)

  ### Cholesky factor
  (C <- ltMatrices(runif(J * (J + 1) / 2), diag = TRUE))
  ### corresponding covariance matrix
  (S <- as.array(Tcrossprod(C))[,,1])

  ### plain Monte-Carlo (Genz, 1992)
  w <- NULL
  M <- 25000
  ### quasi-Monte-Carlo (Genz & Bretz, 2002, but with different weights)
  if (require("qrng")) w <- t(ghalton(M * N, J - 1))

  ### log-likelihood
  lpmvnorm(lower = lwr, upper = upr, chol = C, w = w, M = M)

  ### compare with pmvnorm
  exp(lpmvnorm(lower = lwr, upper = upr, chol = C, logLik = FALSE, w = w, M = M))
  sapply(1:N, function(i) pmvnorm(lower = lwr[,i], upper = upr[,i], sigma = S))

  ### log-lik contributions and score matrix
  slpmvnorm(lower = lwr, upper = upr, chol = C, w = w, M = M, logLik = TRUE)

Run the code above in your browser using DataLab

Last chance! 50% off unlimited learning