lpmvnorm {mvtnorm}

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, ...) 

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.

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.

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.

See Also

dmvnorm, vignette("lmvnorm_src", package = "mvtnorm")

Examples

  ### 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)