Multivariate Normal Log-likelihood and Score Functions

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, invcholmean, center = NULL, chol, invchol, logLik = TRUE, 
         M = NULL, w = NULL, seed = NULL, tol = .Machine$double.eps, fast = FALSE)
slpmvnorm(lower, upper, mean, invcholmean, center = NULL, chol, invchol, logLik = TRUE, 
          M = NULL, w = NULL, seed = NULL, tol = .Machine$double.eps, fast = FALSE)
ldmvnorm(obs, mean, invcholmean, chol, invchol, logLik = TRUE) 
sldmvnorm(obs, mean, invcholmean, chol, invchol, logLik = TRUE) 
ldpmvnorm(obs, lower, upper, mean, invcholmean, chol, invchol, logLik = TRUE, ...) 
sldpmvnorm(obs, lower, upper, mean, invcholmean, 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).
invcholmean: matrix of means left-multiplied with inverse Cholesky factor (invchol %*% mean, 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 mvtnorm::numerical-:1992 with a quasi-Monte-Carlo approach mvtnorm::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 |mvtnorm::numerical-:1992|the default and quasi-Monte-Carlo mvtnorm::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

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))
#> , , 1
#> 
#>            1          2          3          4
#> 1 0.10241835          .          .          .
#> 2 0.09583977 0.04904382          .          .
#> 3 0.67145689 0.36110032 0.10186426          .
#> 4 0.53815421 0.99842188 0.96296355 0.08205368
#> 
  ### corresponding covariance matrix
  (S <- as.array(Tcrossprod(C))[,,1])
#>             1           2          3          4
#> 1 0.010489519 0.009815752 0.06876951 0.05511687
#> 2 0.009815752 0.011590559 0.08206202 0.10054300
#> 3 0.068769509 0.082062016 0.59162412 0.81996938
#> 4 0.055116869 0.100543003 0.81996938 2.22048781

  ### 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))
#> Loading required package: qrng

  ### log-likelihood
  lpmvnorm(lower = lwr, upper = upr, chol = C, w = w, M = M)
#> [1] -13.11619

  ### compare with pmvnorm
  exp(lpmvnorm(lower = lwr, upper = upr, chol = C, logLik = FALSE, w = w, M = M))
#> [1] 0.10183476 0.02499148 0.10141961 0.06908707 0.11285044
  sapply(1:N, function(i) pmvnorm(lower = lwr[,i], upper = upr[,i], sigma = S))
#> [1] 0.10164231 0.02494297 0.10141648 0.06936283 0.11265149

  ### log-lik contributions and score matrix
  slpmvnorm(lower = lwr, upper = upr, chol = C, w = w, M = M, logLik = TRUE)
#> $logLik
#> [1] -2.284404 -3.689220 -2.288489 -2.672388 -2.181692
#> 
#> $mean
#>            [,1]          [,2]          [,3]          [,4]          [,5]
#> 1 -3.553069e-01 -1.064168e+01  1.301111e-02  3.415473e+00 -3.010876e-08
#> 2  5.789007e-19  8.644205e-06  3.400694e-07  5.369940e-20  8.310564e-03
#> 3 -1.103413e+00  2.226243e+00 -2.504700e-01 -4.709515e-01  1.239629e+00
#> 4  5.854073e-01 -6.551919e-01  1.411593e-02  1.417701e-01 -5.795261e-01
#> 
#> $lower
#>            [,1]          [,2]          [,3]          [,4]          [,5]
#> 1 -4.058539e-05 -1.998014e-03 -1.301111e-02 -3.415473e+00  9.077472e-16
#> 2 -7.519498e-19 -8.644205e-06 -3.400694e-07  1.048279e-50 -8.310564e-03
#> 3 -7.548140e-01 -2.233533e+00 -2.046217e+00 -2.078990e-01 -1.772633e+00
#> 4 -1.192508e+00 -2.819646e+00 -6.842197e-01 -2.525375e+00 -6.164127e-01
#> 
#> $upper
#>           [,1]         [,2]          [,3]          [,4]         [,5]
#> 1 3.553474e-01 1.064368e+01 -1.302781e-18  7.033283e-15 3.010876e-08
#> 2 1.730491e-19 1.599524e-74  5.737794e-48 -5.369940e-20 4.545062e-14
#> 3 1.858227e+00 7.289595e-03  2.296687e+00  6.788505e-01 5.330038e-01
#> 4 6.071002e-01 3.474837e+00  6.701038e-01  2.383605e+00 1.195939e+00
#> 
#> $chol
#> , , 1
#> 
#>               1             2             3             4
#> 1 -2.254212e-01             .             .             .
#> 2 -9.874419e-19 -8.938869e-18             .             .
#> 3 -8.218999e-01 -3.743512e-01  2.536208e-01             .
#> 4 -2.272663e-02 -2.768887e-01 -4.170372e-01 -4.504786e-02
#> 
#> , , 2
#> 
#>               1             2             3             4
#> 1 -6.227313e+00             .             .             .
#> 2 -6.397053e-06 -3.511640e-05             .             .
#> 3 -2.698669e-01 -1.737340e-01 -3.039135e-02             .
#> 4  3.749290e-02 -4.859844e-01 -4.756522e-01 -8.784966e-02
#> 
#> , , 3
#> 
#>               1             2             3             4
#> 1 -3.646912e-02             .             .             .
#> 2 -5.424383e-07 -1.431659e-06             .             .
#> 3 -1.472648e+00  1.387406e-02  6.959499e-01             .
#> 4  2.944348e-01 -4.190606e-01 -6.501086e-01 -4.685569e-02
#> 
#> , , 4
#> 
#>               1             2             3             4
#> 1 -2.671198e+00             .             .             .
#> 2  9.736088e-20  4.626732e-19             .             .
#> 3 -7.571720e-01 -1.985563e-01  5.433594e-01             .
#> 4  8.770663e-02 -3.843693e-01 -6.598660e-01 -2.639458e-02
#> 
#> , , 5
#> 
#>               1             2             3             4
#> 1 -1.679813e-07             .             .             .
#> 2 -1.090520e-02  2.333515e-03             .             .
#> 3 -7.165687e-01 -2.168700e-01  7.003675e-02             .
#> 4 -4.706339e-02 -3.404545e-01 -4.033148e-01 -2.587271e-02
#> 
#>