### Load mvtnorm package and set up example parameters

Source: https://cran.r-project.org/web/packages/mvtnorm/vignettes/lmvnorm_src.Rnw

Loads the 'mvtnorm' package and sets random seed for reproducibility. It also defines dimensions N and J, and names for rows and columns for subsequent examples.

```R
library("mvtnorm")
set.seed(290875)
N <- 4L
J <- 5L
rn <- paste0("C_", 1:N)
nm <- LETTERS[1:J]
Jn <- J * (J - 1) / 2
```

--------------------------------

### Kronecker Product Example

Source: https://cran.r-project.org/web/packages/mvtnorm/vignettes/lmvnorm_src.Rnw

Demonstrates the use of the Kronecker product and related matrix operations in R. This example involves setting up matrices, calculating their Kronecker product, and performing checks.

```R
J <- 10
```

```R
d <- TRUE
```

```R
L <- diag(J)
```

```R
L[lower.tri(L, diag = d)] <- prm <- runif(J * (J + c(-1, 1)[d + 1]) / 2)
```

```R
C <- solve(L)
```

```R
D <- -kronecker(t(C), C)
```

```R
S <- diag(J)
```

```R
S[lower.tri(S, diag = TRUE)] <- x <- runif(J * (J + 1) / 2)
```

```R
SD0 <- matrix(c(S) %*% D, ncol = J)
```

```R
SD1 <- -crossprod(C, tcrossprod(S, C))
```

```R
a <- ltMatrices(C[lower.tri(C, diag = TRUE)], diag = TRUE, byrow = FALSE)
```

```R
b <- ltMatrices(x, diag = TRUE, byrow = FALSE)
```

```R
SD2 <- -vectrick(a, b, a)
```

```R
SD2a <- -vectrick(a, b)
```

```R
chk(SD2, SD2a)
```

```R
chk(SD0[lower.tri(SD0, diag = d)], 
    SD1[lower.tri(SD1, diag = d)])
```

```R
chk(SD0[lower.tri(SD0, diag = d)],
    c(unclass(SD2)))
```

```R
S <- t(matrix(as.array(b), byrow = FALSE, nrow = 1))
```

```R
SD2 <- -vectrick(a, S, a)
```

```R
SD2a <- -vectrick(a, S)
```

```R
chk(SD2, SD2a)
```

```R
chk(c(SD0), c(SD2))
```

--------------------------------

### Setup Iris Model

Source: https://cran.r-project.org/web/packages/mvtnorm/vignettes/lmvnorm_src.Rnw

Example of initializing an mvnorm model using maximum-likelihood estimates from the iris dataset.

```R
data("iris", package = "datasets")
vars <- names(iris)[-5L]
N <- nrow(iris)
m <- colMeans(iris[,vars])
V <- var(iris[,vars]) * (N - 1) / N
iris_mvn <- mvnorm(mean = m, chol = t(chol(V)))
iris_var <- simulate(iris_mvn, nsim = nrow(iris))
```

--------------------------------

### Univariate Score Function Example

Source: https://cran.r-project.org/web/packages/mvtnorm/vignettes/lmvnorm_src.Rnw

Demonstrates the score function for univariate problems, comparing direct probability calculations with the output of `slpmvnorm`. Requires prior setup of matrix 'mC'.

```R
ptr <- pnorm(b[1,] / c(unclass(mC[,1]))) - pnorm(a[1,] / c(unclass(mC[,1])))
log(ptr)
```

```R
lpmvnorm(a[1,,drop = FALSE], b[1,,drop = FALSE], chol = mC[,1], logLik = FALSE)
```

```R
lapply(slpmvnorm(a[1,,drop = FALSE], b[1,,drop = FALSE], chol = mC[,1], 
                 logLik = TRUE), unclass)
```

```R
sd1 <- c(unclass(mC[,1]))
(dnorm(b[1,] / sd1) * b[1,] - dnorm(a[1,] / sd1) * a[1,]) * (-1) / sd1^2 / ptr
```

--------------------------------

### Example usage of Tcrossprod

Source: https://cran.r-project.org/web/packages/mvtnorm/vignettes/lmvnorm_src.Rnw

Demonstrates how to use the Tcrossprod function and verify results against base R implementations.

```R
## Tcrossprod
a <- as.array(Tcrossprod(lxn))
b <- array(apply(as.array(lxn), 3L, function(x) tcrossprod(x), simplify = TRUE), 
           dim = rev(dim(lxn)))
chk(a, b, check.attributes = FALSE)

# diagonal elements only
d <- Tcrossprod(lxn, diag_only = TRUE)
chk(d, apply(a, 3, diag))
chk(d, diagonals(Tcrossprod(lxn)))

a <- as.array(Tcrossprod(lxd))
b <- array(apply(as.array(lxd), 3L, function(x) tcrossprod(x), simplify = TRUE), 
           dim = rev(dim(lxd)))
chk(a, b, check.attributes = FALSE)
```

--------------------------------

### Estimate Multivariate Normal Probabilities Example

Source: https://cran.r-project.org/web/packages/mvtnorm/vignettes/lmvnorm_src.Rnw

This example demonstrates how to use the `lpmvnormR` function to estimate probabilities for a multivariate normal distribution. It sets up random lower and upper bounds and a Cholesky factor, then calls the function.

```R
J <- 5L
N <- 10L

x <- matrix(runif(N * J * (J + 1) / 2), ncol = N)
lx <- ltMatrices(x, byrow = TRUE, diag = TRUE)

a <- matrix(runif(N * J), nrow = J) - 2
a[sample(J * N)[1:2]] <- -Inf
b <- a + 2 + matrix(runif(N * J), nrow = J)
b[sample(J * N)[1:2]] <- Inf

(phat <- c(lpmvnormR(a, b, chol = lx, logLik = FALSE)))
```

--------------------------------

### Setup Return Vector in C

Source: https://cran.r-project.org/web/packages/mvtnorm/vignettes/lmvnorm_src.Rnw

Allocates and initializes a real vector of a specified length for the return value in R. Sets all elements to 0.0.

```c
len = (RlogLik ? 1 : iN);
PROTECT(ans = allocVector(REALSXP, len));
dans = REAL(ans);
for (int i = 0; i < len; i++)
    dans[i] = 0.0;
```

--------------------------------

### Maximum Likelihood Estimation Setup

Source: https://cran.r-project.org/web/packages/mvtnorm/vignettes/lmvnorm_src.Rnw

Demonstrates setting up maximum likelihood estimation procedures using the log-likelihood and score functions for the iris dataset.

```APIDOC
## Maximum Likelihood Estimation Setup

### Description
This section illustrates how to set up maximum likelihood estimation (MLE) for a multivariate normal distribution using the `mvtnorm` package. It defines functions for the log-likelihood (`ll`) and the score function (`sc`), and demonstrates their use with the iris dataset to estimate the mean and Cholesky factor of the covariance matrix.

### Method
Conceptual setup for MLE.

### Endpoint
N/A (This is a code example for setting up estimation procedures)

### Parameters
None

### Request Example
```R
# Assuming 'vars' is a character vector of variable names
# Assuming 'iris' is the iris dataset

J <- length(vars)
obs <- t(iris[, vars])

# Log-likelihood function
ll <- function(parm) {
    C <- ltMatrices(parm[-(1:J)], diag = TRUE, names = vars)
    x <- mvnorm(mean = parm[1:J], chol = C)
    -logLik(object = x, obs = obs)
}

# Score function
sc <- function(parm) {
    C <- ltMatrices(parm[-(1:J)], diag = TRUE, names = vars)
    x <- mvnorm(mean = parm[1:J], chol = C)
    ret <- lLgrad(object = x, obs = obs)
    -c(rowSums(ret$mean), rowSums(Lower_tri(ret$scale, diag = TRUE)))
}

# To estimate parameters, one would typically use an optimization function:
# optim(par = initial_parm, fn = ll, gr = sc, method = "L-BFGS-B")
```

### Response
#### Success Response (200)
- **ll** (function) - Log-likelihood function.
- **sc** (function) - Score function (gradient of log-likelihood).

#### Response Example
```json
{
  "ll": "function(parm) { ... }",
  "sc": "function(parm) { ... }"
}
```
```

--------------------------------

### Example Usage of margDist and condDist

Source: https://cran.r-project.org/web/packages/mvtnorm/vignettes/lmvnorm_src.Rnw

This example demonstrates computing the marginal distribution for specific variables and then using the result to compute a conditional distribution. It assumes `iris_mvn` is a pre-defined multivariate normal object and `vars` is a vector of variable names.

```R
j <- 3:4
margDist(iris_mvn, which = vars[j])
gm <- t(iris[,vars[-(j)]])
iris_cmvn <- condDist(iris_mvn, which_given = vars[j], given = gm)
```

--------------------------------

### Test Setup for lpRR

Source: https://cran.r-project.org/web/packages/mvtnorm/vignettes/lmvnorm_src.Rnw

Initializes variables for testing the lpRR functionality with specific dimensions.

```R
J <- 6
K <- 3
B <- matrix(rnorm(J * K), nrow = J)
D <- runif(J)
S <- tcrossprod(B) + diag(D)
```

--------------------------------

### Gaussian Copula Classical Model Setup

Source: https://cran.r-project.org/web/packages/mvtnorm/vignettes/lmvnorm_src.Rnw

Sets up log-likelihood and score functions for a Gaussian copula model using the iris dataset.

```R
data("iris", package = "datasets")
J <- 4
Z <- t(qnorm(do.call("cbind", lapply(iris[1:J], rank, ties.method = "max")) / 
       (nrow(iris) + 1)))
(CR <- cor(t(Z)))
ll <- function(parm) {
    C <- ltMatrices(parm)
    Cs <- standardize(chol = C)
    -ldmvnorm(obs = Z, chol = Cs)
}
sc <- function(parm) {
    C <- ltMatrices(parm)
    Cs <- standardize(chol = C)
    -rowSums(Lower_tri(destandardize(chol = C, 
        score_schol = sldmvnorm(obs = Z, chol = Cs)$chol)))
}
start <- t(chol(CR))
start <- start[lower.tri(start)]
if (require("numDeriv", quietly = TRUE))
    chk(grad(ll, start), sc(start), check.attributes = FALSE)
op <- optim(start, fn = ll, gr = sc, method = "BFGS", 
            control = list(trace = FALSE), hessian = TRUE)
op$value
S_ML <- chol2cov(standardize(chol = ltMatrices(op$par)))
```

--------------------------------

### Solve Matrix Equations

Source: https://cran.r-project.org/web/packages/mvtnorm/vignettes/lmvnorm_src.Rnw

Examples for solving linear systems with matrix objects.

```R
chk(solve(lxn, y[,1]), solve(lxn, y[,rep(1, N)]))
```

```R
chk(solve(lxn[1,], y, transpose = TRUE), 
    t(as.array(solve(lxn[1,]))[,,1]) %*% y)
```

--------------------------------

### Test Lower_tri function

Source: https://cran.r-project.org/web/packages/mvtnorm/vignettes/lmvnorm_src.Rnw

Examples demonstrating the usage of the Lower_tri function.

```R
## J <- 4
M <- ltMatrices(matrix(1:10, nrow = 10, ncol = 2), diag = TRUE)
Lower_tri(M, diag = FALSE)
Lower_tri(M, diag = TRUE)
M <- ltMatrices(matrix(1:6, nrow = 6, ncol = 2), diag = FALSE)
Lower_tri(M, diag = FALSE)
Lower_tri(M, diag = TRUE)
```

--------------------------------

### Example: Evaluating Log-Likelihood with mvtnorm in R

Source: https://cran.r-project.org/web/packages/mvtnorm/vignettes/lmvnorm_src.Rnw

Demonstrates evaluating the log-likelihood using the lpmvnorm function with pre-computed Cholesky factors and specified Monte Carlo iterations.

```R
phat
exp(lpmvnorm(a, b, chol = lx, M = 25000, logLik = FALSE, fast = TRUE))
exp(lpmvnorm(a, b, chol = lx, M = 25000, logLik = FALSE, fast = FALSE))
```

--------------------------------

### Example usage of syMatrices and Mult

Source: https://cran.r-project.org/web/packages/mvtnorm/vignettes/lmvnorm_src.Rnw

Demonstrates the usage of syMatrices for creating symmetric matrices and the Mult function for matrix multiplication. Includes checks for different configurations of matrix dimensions and diagonal properties.

```R
J <- 5
N1 <- 10
ex <- expression({
  C <- syMatrices(matrix(runif(N2 * J * (J + c(-1, 1)[DIAG + 1L] ) / 2),
                          ncol = N2),
                  diag = DIAG)
  x <- matrix(runif(N1 * J), nrow = J)
  Ca <- as.array(C)
  p1 <- do.call("cbind", lapply(1:N1, function(i) 
      Ca[,,c(1,i)[(N2 > 1) + 1]] %*% x[,i]))
  p2 <- Mult(C, x)
  chk(p1, p2)
})

N2 <- N1
DIAG <- TRUE
eval(ex)
N2 <- 1
DIAG <- TRUE
eval(ex)
N2 <- 1
DIAG <- FALSE
eval(ex)
N2 <- N1
DIAG <- FALSE
eval(ex)
```

--------------------------------

### Calculate Ground Truth Probabilities

Source: https://cran.r-project.org/web/packages/mvtnorm/vignettes/lmvnorm_src.Rnw

Example of calculating probabilities using pnorm for comparison.

```R
ptr <- pnorm(b[1,] / c(unclass(lx[,1]))) - pnorm(a[1,] / c(unclass(lx[,1])))

cbind(c(ptr), pGB, pGq)
```

--------------------------------

### Univariate Problem Integration with lpmvnorm

Source: https://cran.r-project.org/web/packages/mvtnorm/vignettes/lmvnorm_src.Rnw

Tests univariate problem integration using the `lpmvnorm` function. This example uses default settings for `w` and `fast` parameters.

```R
pGq <- exp(lpmvnorm(a[1,,drop = FALSE], b[1,,drop = FALSE], chol = lx[,1], 
                   logLik = FALSE))
```

--------------------------------

### Create and Modify Lower Triangular Matrix in R

Source: https://cran.r-project.org/web/packages/mvtnorm/vignettes/lmvnorm_src.Rnw

Example demonstrating the creation of a lower triangular matrix and setting its diagonal values.

```R
lxd2 <- lxn
diagonals(lxd2) <- 1
chk(as.array(lxd2), as.array(lxn))
```

--------------------------------

### Calculate Multivariate Normal Probabilities

Source: https://cran.r-project.org/web/packages/mvtnorm/news.html

Examples demonstrating the use of pmvnorm with Miwa and TVPACK algorithms for dimension reduction scenarios.

```R
  pmvnorm(lower = rep(-Inf,3), upper = c(-1, Inf, Inf),
          sigma = diag(3), algorithm = Miwa())
     # or
  pmvnorm(lower = c(-Inf,-Inf),
          upper = c(- 1,  Inf), sigma=diag(2), algorithm = TVPACK())
```

--------------------------------

### Example Usage of lpRR and slpRR

Source: https://cran.r-project.org/web/packages/mvtnorm/refman/mvtnorm.html

Demonstrates the usage of lpRR and slpRR functions for calculating log-likelihood and score functions with specified parameters. Ensure necessary matrices and vectors are defined before use.

```R
J <- 6
K <- 3
B <- matrix(rnorm(J * K), nrow = J)
D <- runif(J)
S <- tcrossprod(B) + diag(D)
a <- -(2 + runif(J))
b <- 2 + runif(J)
M <- 1e4
Z <- matrix(rnorm(K * M), nrow = K)
## log-likelihood
lpRR(lower = a, upper = b, B = B, D = D, Z = Z)
## score wrt all arguments
slpRR(lower = a, upper = b, B = B, D = D, Z = Z)
```

--------------------------------

### Gaussian Copula Nonparametric Likelihood Setup

Source: https://cran.r-project.org/web/packages/mvtnorm/vignettes/lmvnorm_src.Rnw

Defines intervals for nonparametric maximum likelihood estimation of copula parameters.

```R
lwr <- do.call("cbind", lapply(iris[1:J], rank, ties.method = "min")) - 1L
upr <- do.call("cbind", lapply(iris[1:J], rank, ties.method = "max"))
lwr <- t(qnorm(lwr / nrow(iris)))
upr <- t(qnorm(upr / nrow(iris)))

M <- 500 
if (require("qrng", quietly = TRUE)) {
    ### quasi-Monte-Carlo
    W <- t(ghalton(M, d = J - 1))
} else {
    ### Monte-Carlo
    W <- matrix(runif(M * (J - 1)), nrow = J - 1, byrow = TRUE)
}

ll <- function(parm) {
    C <- ltMatrices(parm)
    Cs <- standardize(chol = C)
```

--------------------------------

### Example Usage of Cholesky Factorization

Source: https://cran.r-project.org/web/packages/mvtnorm/vignettes/lmvnorm_src.Rnw

Demonstrates the usage of the chol function to revert tcrossprod for ltMatrices objects, with and without unit diagonals.

```R
Sigma <- tcrossprod(lxd)
chk(chol(Sigma), lxd)
Sigma <- tcrossprod(lxn)

```

--------------------------------

### Calculate Multivariate Normal Probabilities

Source: https://cran.r-project.org/web/packages/mvtnorm/vignettes/lmvnorm_src.Rnw

Example usage of lpmvnorm and lpRR for calculating probabilities and gradients.

```R
Linv <- t(chol(S))
Linv <- ltMatrices(Linv[lower.tri(Linv, diag = TRUE)], diag = TRUE)
a <- -(2 + runif(J))
b <- 2 + runif(J)
M <- 1e6
dim(w <- matrix(runif((J - 1) * M), nrow = J - 1, byrow = TRUE))
lpmvnorm(lower = a, upper = b, chol = Linv, w = w)
dim(Z <- matrix(rnorm(K * M), nrow = K))
lpRR(lower = a, upper = b, B = B, D = D, Z = Z)
```

--------------------------------

### Compute Multivariate t Probabilities

Source: https://cran.r-project.org/web/packages/mvtnorm/refman/mvtnorm.html

Examples demonstrating the calculation of multivariate t probabilities using pmvt and comparing results with univariate t distributions.

```R
n <- 5
lower <- -1
upper <- 3
df <- 4
corr <- diag(5)
corr[lower.tri(corr)] <- 0.5
delta <- rep(0, 5)
prob <- pmvt(lower=lower, upper=upper, delta=delta, df=df, corr=corr)
print(prob)

pmvt(lower=-Inf, upper=3, df = 3, sigma = 1) == pt(3, 3)

# Example from R News paper (original by Edwards and Berry, 1987)

n <- c(26, 24, 20, 33, 32)
V <- diag(1/n)
df <- 130
C <- c(1,1,1,0,0,-1,0,0,1,0,0,-1,0,0,1,0,0,0,-1,-1,0,0,-1,0,0)
C <- matrix(C, ncol=5)
### scale matrix
cv <- C %*% tcrossprod(V, C)
### correlation matrix
cr <- cov2cor(cv)
delta <- rep(0,5)

myfct <- function(q, alpha) {
  lower <- rep(-q, ncol(cv))
  upper <- rep(q, ncol(cv))
  pmvt(lower=lower, upper=upper, delta=delta, df=df,
       corr=cr, abseps=0.0001) - alpha
}

### uniroot for this simple problem
round(uniroot(myfct, lower=1, upper=5, alpha=0.95)$root, 3)

# compare pmvt and pmvnorm for large df:

a <- pmvnorm(lower=-Inf, upper=1, mean=rep(0, 5), corr=diag(5))
b <- pmvt(lower=-Inf, upper=1, delta=rep(0, 5), df=300,
          corr=diag(5))
a
b

stopifnot(round(a, 2) == round(b, 2))

# correlation and scale matrix

a <- pmvt(lower=-Inf, upper=2, delta=rep(0,5), df=3,
          sigma = diag(5)*2)
b <- pmvt(lower=-Inf, upper=2/sqrt(2), delta=rep(0,5),
          df=3, corr=diag(5))
attributes(a) <- NULL
attributes(b) <- NULL
a
b
stopifnot(all.equal(round(a,3) , round(b, 3)))

a <- pmvt(0, 1,df=10)
attributes(a) <- NULL
b <- pt(1, df=10) - pt(0, df=10)
stopifnot(all.equal(round(a,10) , round(b, 10)))
```

--------------------------------

### Setup for numerical comparisons

Source: https://cran.r-project.org/web/packages/mvtnorm/vignettes/lmvnorm_src.Rnw

Defines a utility function 'chk' for stopping execution if 'all.equal' comparison of two objects returns FALSE. This is useful for verifying numerical results.

```R
chk <- function(...) stopifnot(isTRUE(all.equal(...)))
```

--------------------------------

### Initialize Log-Likelihood Loop Variables

Source: https://cran.r-project.org/web/packages/mvtnorm/vignettes/lmvnorm_src.Rnw

Initializes variables for the log-likelihood calculation loop. Sets initial values for d, f, emd, and start.

```R
d = d0;
f = f0;
emd = emd0;
start = 0;
```

--------------------------------

### Maximum Likelihood Estimation Setup for Iris Dataset

Source: https://cran.r-project.org/web/packages/mvtnorm/vignettes/lmvnorm_src.Rnw

Sets up functions for calculating the log-likelihood and score for maximum likelihood estimation using the iris dataset.

```R
J <- length(vars)
obs <- t(iris[, vars])
ll <- function(parm) {
    C <- ltMatrices(parm[-(1:J)], diag = TRUE, names = vars)
    x <- mvnorm(mean = parm[1:J], chol = C)
    -logLik(object = x, obs = obs)
}
sc <- function(parm) {
    C <- ltMatrices(parm[-(1:J)], diag = TRUE, names = vars)
    x <- mvnorm(mean = parm[1:J], chol = C)
    ret <- lLgrad(object = x, obs = obs)
    -c(rowSums(ret$mean), rowSums(Lower_tri(ret$scale, diag = TRUE)))
}
```

--------------------------------

### Initialize parameters for numerical differentiation

Source: https://cran.r-project.org/web/packages/mvtnorm/vignettes/lmvnorm_src.Rnw

Sets up a positive definite matrix and its Cholesky decomposition for testing purposes.

```R
J <- 5L
N <- 4L

S <- crossprod(matrix(runif(J^2), nrow = J))
prm <- t(chol(S))[lower.tri(S, diag = TRUE)]
```

--------------------------------

### Initialize R Environment

Source: https://cran.r-project.org/web/packages/mvtnorm/vignettes/MVT_Rnews.R

Sets the random seed for reproducibility.

```R
set.seed(290875)
```

--------------------------------

### Subset with character vector

Source: https://cran.r-project.org/web/packages/mvtnorm/vignettes/lmvnorm_src.Rnw

Example of subsetting using a character vector.

```R
j <- nm[c(1, 3, 5)]
```

--------------------------------

### Compare qmvt with qmvnorm for df=0 and df=Inf

Source: https://cran.r-project.org/web/packages/mvtnorm/refman/mvtnorm.html

Compares the quantiles computed by `qmvt` with `df = 0` and `df = Inf` against those computed by `qmvnorm`. This is useful for verifying the behavior of `qmvt` in limiting cases.

```r
Sigma <- diag(2)
set.seed(29)
q0 <- qmvt(0.95, sigma = Sigma, df = 0,   tail = "both")$quantile
set.seed(29)
q8 <- qmvt(0.95, sigma = Sigma, df = Inf, tail = "both")$quantile
set.seed(29)
qn <- qmvnorm(0.95, sigma = Sigma, tail = "both")$quantile
stopifnot(identical(q0, q8),
          isTRUE(all.equal(q0, qn, tol = (.Machine$double.eps)^(1/3))))
```

--------------------------------

### Subset with specific indices

Source: https://cran.r-project.org/web/packages/mvtnorm/vignettes/lmvnorm_src.Rnw

Example of subsetting using a specific numeric vector.

```R
j <- c(1, 3, 5)
```

--------------------------------

### Subset with character indices

Source: https://cran.r-project.org/web/packages/mvtnorm/vignettes/lmvnorm_src.Rnw

Example of subsetting using character vector indices.

```R
i <- 1:2
j <- nm[2:4]
```

--------------------------------

### Initialize Parameters for Convex Optimization

Source: https://cran.r-project.org/web/packages/mvtnorm/vignettes/lmvnorm_src.Rnw

Sets up initial parameters for the convex optimization of the multivariate normal log-likelihood. Includes generating a mean vector and an inverse Cholesky factor.

```R
J <- 5
N <- 100
### mean
m <- rnorm(J)
L <- ltMatrices(prm <- runif(J * (J + 1) / 2), diag = TRUE)
Z <- matrix(rnorm(N * J), nrow = J)
Y <- solve(L, Z) + m
### scaled mean
d <- L %*% m

start <- c(d, prm)
```

--------------------------------

### Subset with positive integers

Source: https://cran.r-project.org/web/packages/mvtnorm/vignettes/lmvnorm_src.Rnw

Example of subsetting using positive integer indices.

```R
i <- colnames(xn)[1:2]
j <- 2:4
```

--------------------------------

### Subset with non-increasing indices

Source: https://cran.r-project.org/web/packages/mvtnorm/vignettes/lmvnorm_src.Rnw

Example of subsetting with non-increasing indices, applicable to symmetric matrices.

```R
## subset
j <- nm[sample(1:J)]
ltM <- ltMatrices(xn, byrow = FALSE, names = nm)
try(ltM[i, j])
ltM <- as.syMatrices(ltM)
a <- as.array(ltM[i, j])
b <- as.array(ltM)[j, j, i]
chk(a, b)
```

--------------------------------

### Initialize ltMatrices and syMatrices

Source: https://cran.r-project.org/web/packages/mvtnorm/vignettes/lmvnorm_src.Rnw

Demonstrates creating lower triangular and symmetric matrix objects from random data.

```R
xn <- matrix(runif(N * Jn), ncol = N)
colnames(xn) <- rn
xd <- matrix(runif(N * (Jn + J)), ncol = N)
colnames(xd) <- rn

(lxn <- ltMatrices(xn, byrow = TRUE, names = nm))
dim(lxn)
dimnames(lxn)
lxd <- ltMatrices(xd, byrow = TRUE, diag = TRUE, names = nm)
dim(lxd)
dimnames(lxd)

lxn <- as.syMatrices(lxn)
lxn
```

--------------------------------

### Solving Linear Systems with solve()

Source: https://cran.r-project.org/web/packages/mvtnorm/vignettes/lmvnorm_src.Rnw

Demonstrates the usage of the solve() function for inverting matrices and solving linear systems. Includes checks for correctness against array operations and matrix multiplication.

```R
## solve
A <- as.array(lxn)
a <- solve(lxn)
a <- as.array(a)
b <- array(apply(A, 3L, function(x) solve(x), simplify = TRUE), 
           dim = rev(dim(lxn)))
chk(a, b, check.attributes = FALSE)
```

```R
A <- as.array(lxd)
a <- as.array(solve(lxd))
b <- array(apply(A, 3L, function(x) solve(x), simplify = TRUE), 
           dim = rev(dim(lxd)))
chk(a, b, check.attributes = FALSE)
```

```R
chk(solve(lxn, y), Mult(solve(lxn), y))
```

```R
chk(solve(lxd, y), Mult(solve(lxd), y))
```

```R
### recycle C
chk(solve(lxn[1,], y), as.array(solve(lxn[1,]))[,,1] %*% y)
```

```R
chk(solve(lxn[rep(1, N),], y), solve(lxn[1,], y), check.attributes = FALSE)
```

--------------------------------

### Subset with negative indices

Source: https://cran.r-project.org/web/packages/mvtnorm/vignettes/lmvnorm_src.Rnw

Example of subsetting using negative indices to exclude elements.

```R
j <- -c(1, 3, 5)
```

--------------------------------

### Basic Usage of qmvt

Source: https://cran.r-project.org/web/packages/mvtnorm/refman/mvtnorm.html

Demonstrates the basic usage of the `qmvt` function to compute a quantile for a multivariate t distribution with specified degrees of freedom and tail type.

```r
qmvt(0.95, df = 16, tail = "both")
```

--------------------------------

### Compute Crossproduct First Element in C

Source: https://cran.r-project.org/web/packages/mvtnorm/vignettes/lmvnorm_src.Rnw

C implementation snippet for calculating the first element of a crossproduct matrix.

```C
dans[0] = 1.0;
if (Rdiag)
    dans[0] = pow(dC[0], 2);
if (Rtranspose) { // crossprod
    for (k = 1; k < iJ; k++) 
        dans[0] += pow(dC[IDX(k + 1, 1, iJ, Rdiag)], 2);
}
```

--------------------------------

### Evaluate Log-Likelihood for Iris Data

Source: https://cran.r-project.org/web/packages/mvtnorm/vignettes/lmvnorm_src.Rnw

Examples of evaluating the log-likelihood for conditional multivariate normal models using the iris dataset.

```R
logLik(object = iris_cmvn, obs = t(iris[,vars[-j]]))
```

```R
logLik(object = iris_cmvn, obs = t(iris[,rev(vars[-j])]))
```

--------------------------------

### LAPACK options for solving linear systems

Source: https://cran.r-project.org/web/packages/mvtnorm/vignettes/lmvnorm_src.Rnw

Sets up options for LAPACK routines, specifically determining whether to handle non-unit ('N') or unit ('U') diagonal elements for matrix operations.

```C
char di, lo = 'L';
if (Rdiag) {
    /* non-unit diagonal elements */
    di = 'N';
} else {
    /* unit diagonal elements; NOTE: these diagonals 1s ARE always present but
       ignored in the computations */
    di = 'U';
}
```

--------------------------------

### Initialize Optimization Parameters for MVN

Source: https://cran.r-project.org/web/packages/mvtnorm/vignettes/lmvnorm_src.Rnw

Sets initial values for optimization by rounding the mean and lower triangular elements of the scale matrix. Defines lower bounds for optimization, ensuring positive values for diagonal elements of the Cholesky factor.

```R
start <- round(c(c(iris_mvn$mean), 
                 Lower_tri(iris_mvn$scale, diag = TRUE)), 2)
llim <- rep(-Inf, J + J * (J + 1) / 2)
llim[J + c(diagonals(ltMatrices(seq_len(J * (J + 1) / 2), diag = TRUE)))] <- 1e-4
```

--------------------------------

### Evaluate Multivariate Normal Log-Likelihood

Source: https://cran.r-project.org/web/packages/mvtnorm/vignettes/lmvnorm_src.Rnw

Examples demonstrating efficient log-likelihood evaluation using Cholesky factors compared to traditional methods.

```R
N <- 1000L
J <- 50L
lt <- ltMatrices(matrix(runif(N * J * (J + 1) / 2) + 1, ncol = N), 
                 diag = TRUE, byrow = FALSE)
Z <- matrix(rnorm(N * J), ncol = N)
Y <- solve(lt, Z)
ll1 <- sum(dnorm(lt %*% Y, log = TRUE)) + sum(log(diagonals(lt)))

S <- as.array(tcrossprod(solve(lt)))
ll2 <- sum(sapply(1:N, function(i) 
                           dmvnorm(x = Y[,i], sigma = S[,,i], log = TRUE)))
chk(ll1, ll2)
```

```R
ll3 <- ldmvnorm(obs = Y, invchol = lt)
chk(ll1, ll3)
```

--------------------------------

### Integration Weights

Source: https://cran.r-project.org/web/packages/mvtnorm/vignettes/lmvnorm_src.Rnw

Handles the initialization and checking of integration weights for Monte Carlo or quasi-Monte Carlo integration.

```APIDOC
## Integration Weights

### Description
Checks and sets integration weights for multivariate normal probability calculations. If weights are not provided, it can generate them using quasi-Monte Carlo methods (e.g., Halton sequences) or standard Monte Carlo methods.

### Method
Internal function logic

### Endpoint
N/A (Internal logic within `lpmvnorm`)

### Parameters
#### Path Parameters
None

#### Query Parameters
None

#### Request Body
None

### Request Example
```R
# Example of generating weights (within lpmvnorm context)
if (require("qrng", quietly = TRUE)) {
    W <- t(ghalton(M, d = J - 1)) # Quasi-Monte Carlo
} else {
    W <- matrix(runif(M * (J - 1)), nrow = J - 1, byrow = TRUE) # Monte Carlo
}
```

### Response
#### Success Response (200)
Integration weights matrix `w` is prepared.

#### Response Example
```R
# Example structure of weights matrix W
# [,1] [,2] ...
# [1,]  0.1  0.5 ...
# [2,]  0.3  0.7 ...
# ...
```
```

--------------------------------

### Initialize Center Vector in C

Source: https://cran.r-project.org/web/packages/mvtnorm/vignettes/lmvnorm_src.Rnw

Sets up a double pointer 'dcenter' for the 'center' R object and validates its dimensions against iN * iJ if it's not empty.

```c
dcenter = REAL(center);
if (LENGTH(center)) {
    if (LENGTH(center) != iN * iJ)
        error("incorrect dimensions of center");
}
```

--------------------------------

### Check Gradient Correctness

Source: https://cran.r-project.org/web/packages/mvtnorm/vignettes/lmvnorm_src.Rnw

Numerically checks the correctness of the gradient function using the 'numDeriv' package. Ensure 'numDeriv' is installed and loaded.

```R
if (require("numDeriv", quietly = TRUE))
    chk(grad(ll, prm, J = J), sc(prm, J = J), check.attributes = FALSE)
```

--------------------------------

### Compare MVN Estimation Results

Source: https://cran.r-project.org/web/packages/mvtnorm/vignettes/lmvnorm_src.Rnw

Compares the log-likelihood values and estimated parameters (covariance and mean) from two different optimization approaches: direct Cholesky factor and inverse Cholesky factor parameterization.

```R
op$value
opL$value
invchol2cov(MLL$scale)
V
MLL$mean[,,drop = TRUE]
m
```

--------------------------------

### Get Dimensions and Names of ltMatrices Object

Source: https://cran.r-project.org/web/packages/mvtnorm/refman/mvtnorm.html

Retrieves the dimensions and names from an ltMatrices object. Useful for inspecting the structure of the matrix data.

```R
dim(C)
dimnames(C)
names(C)
```

--------------------------------

### Initialize score loop

Source: https://cran.r-project.org/web/packages/mvtnorm/vignettes/lmvnorm_src.Rnw

Initializes the loop for score calculation by invoking logLik and score macros.

```C
<init logLik loop 69c>
<score c11 80b>
<score a, b 81a>
```

--------------------------------

### Set Digits for Output

Source: https://cran.r-project.org/web/packages/mvtnorm/vignettes/lmvnorm_src.Rnw

Sets the default number of digits for displaying floating-point numbers in R. Useful for consistent output in examples or reports.

```R
options(digits = 4)
```

--------------------------------

### Construct and Reorder ltMatrices

Source: https://cran.r-project.org/web/packages/mvtnorm/vignettes/lmvnorm_src.Rnw

Demonstrates creating ltMatrices objects and reordering them by toggling the byrow argument.

```R
a <- as.array(ltMatrices(xn, byrow = TRUE))
b <- as.array(ltMatrices(ltMatrices(xn, byrow = TRUE), 
                         byrow = FALSE))
chk(a, b)

a <- as.array(ltMatrices(xn, byrow = FALSE))
b <- as.array(ltMatrices(ltMatrices(xn, byrow = FALSE), 
                         byrow = TRUE))
chk(a, b)

a <- as.array(ltMatrices(xd, byrow = TRUE, diag = TRUE))
b <- as.array(ltMatrices(ltMatrices(xd, byrow = TRUE, diag = TRUE), 
                         byrow = FALSE))
chk(a, b)

a <- as.array(ltMatrices(xd, byrow = FALSE, diag = TRUE))
b <- as.array(ltMatrices(ltMatrices(xd, byrow = FALSE, diag = TRUE), 
                         byrow = TRUE))
chk(a, b)
```

--------------------------------

### Univariate Problem Integration with Genz-Bretz

Source: https://cran.r-project.org/web/packages/mvtnorm/vignettes/lmvnorm_src.Rnw

Tests the Genz-Bretz algorithm for a univariate problem using `lpmvnormR`. Note the use of `a[1,,drop = FALSE]` and `b[1,,drop = FALSE]` for univariate input.

```R
pGB <- lpmvnormR(a[1,,drop = FALSE], b[1,,drop = FALSE], chol = lx[,1], 
                logLik = FALSE, 
                algorithm = GenzBretz(maxpts = M, abseps = 0, releps = 0))
```

--------------------------------

### Get Names of ltMatrices in R

Source: https://cran.r-project.org/web/packages/mvtnorm/vignettes/lmvnorm_src.Rnw

Extracts the names identifying rows and columns within each matrix of an ltMatrices object. This function is also used for syMatrices.

```r
names.ltMatrices <- function(x) {
    return(attr(x, "dimnames")[[1L]])
}
names.syMatrices <- names.ltMatrices
```

--------------------------------

### Evaluate Negative Log-Likelihood with Initial Parameters

Source: https://cran.r-project.org/web/packages/mvtnorm/vignettes/lmvnorm_src.Rnw

Evaluates the negative log-likelihood function using the initialized parameters. This is a preliminary check before optimization.

```R
nll(start)
```

--------------------------------

### Test standard interface methods

Source: https://cran.r-project.org/web/packages/mvtnorm/vignettes/lmvnorm_src.Rnw

Verification of standard tcrossprod and crossprod methods.

```R
### tcrossprod
a <- as.array(tcrossprod(lxn))
b <- array(apply(as.array(lxn), 3L, function(x) tcrossprod(x), simplify = TRUE), 
           dim = rev(dim(lxn)))
chk(a, b, check.attributes = FALSE)

a <- as.array(tcrossprod(lxd))
b <- array(apply(as.array(lxd), 3L, function(x) tcrossprod(x), simplify = TRUE), 
           dim = rev(dim(lxd)))
chk(a, b, check.attributes = FALSE)

## crossprod
a <- as.array(crossprod(lxn))
b <- array(apply(as.array(lxn), 3L, function(x) crossprod(x), simplify = TRUE), 
           dim = rev(dim(lxn)))
chk(a, b, check.attributes = FALSE)

a <- as.array(crossprod(lxd))
b <- array(apply(as.array(lxd), 3L, function(x) crossprod(x), simplify = TRUE), 
           dim = rev(dim(lxd)))
chk(a, b, check.attributes = FALSE)
```

--------------------------------

### Get Dimnames of ltMatrices in R

Source: https://cran.r-project.org/web/packages/mvtnorm/vignettes/lmvnorm_src.Rnw

Extracts the dimnames for an ltMatrices object. It returns a list containing the dimnames for the matrices. This function is also used for syMatrices.

```r
dimnames.ltMatrices <- function(x)
    return(list(attr(x, "dimnames")[[2L]], attr(x, "rcnames"), attr(x, "rcnames")))
dimnames.syMatrices <- dimnames.ltMatrices
```

--------------------------------

### Initialize Dimensions and Data Pointers in C

Source: https://cran.r-project.org/web/packages/mvtnorm/vignettes/lmvnorm_src.Rnw

Initializes integer dimensions (iM, iN, iJ) from R objects and sets up double pointers for data arrays (a, b, C, W). It also determines the 'p' value based on the length of C.

```c
int iM = INTEGER(M)[0]; 
int iN = INTEGER(N)[0]; 
int iJ = INTEGER(J)[0]; 

da = REAL(a);
db = REAL(b);
dC = REAL(C);
dW = REAL(C); // make -Wmaybe-uninitialized happy

if (LENGTH(C) == iJ * (iJ - 1) / 2)
    p = 0;
else 
    p = LENGTH(C) / iN;
```

--------------------------------

### Get Lower Triangular Matrix Diagonals in R

Source: https://cran.r-project.org/web/packages/mvtnorm/vignettes/lmvnorm_src.Rnw

Returns a representation of the diagonals of a lower triangular matrix. If input is an integer, it creates a zero matrix of that size.

```R
diagonals.integer <- function(x, ...) 
    ltMatrices(rep(0, x * (x - 1) / 2), diag = FALSE, ...)
```

--------------------------------

### Genz-1992 Quasi-Monte-Carlo Integration with Fast pnorm

Source: https://cran.r-project.org/web/packages/mvtnorm/vignettes/lmvnorm_src.Rnw

Employs the Genz-1992 algorithm with quasi-Monte-Carlo integration and fast probability density function calculation. Requires pre-defined `W` for quasi-Monte-Carlo.

```R
pGqf <- exp(lpmvnorm(a, b, chol = lx, w = W, M = M, logLik = FALSE, 
                     fast = TRUE))
```

--------------------------------

### Prepare Data for Interval-Censored MVN Estimation

Source: https://cran.r-project.org/web/packages/mvtnorm/vignettes/lmvnorm_src.Rnw

Sets up data for interval-censored multivariate normal estimation. It defines interval bounds for the first variable and removes it from the set of exactly observed variables.

```R
v1 <- vars[1]
q1 <- quantile(iris[[v1]], probs = 1:4 / 5)
head(f1 <- cut(iris[[v1]], breaks = c(-Inf, q1, Inf)))
lower <- matrix(c(-Inf, q1)[f1], nrow = 1)
upper <- matrix(c(q1, Inf)[f1], nrow = 1)
rownames(lower) <- rownames(upper) <- v1
obs <- obs[!rownames(obs) %in% v1,,drop = FALSE]
```

--------------------------------

### Multivariate Normal Probability with Custom Correlation Matrix

Source: https://cran.r-project.org/web/packages/mvtnorm/refman/mvtnorm.html

Example from R News paper (original by Genz, 1992) demonstrating the computation of a multivariate normal probability with a custom correlation matrix.

```R
m <- 3
sigma <- diag(3)
sigma[2,1] <- 3/5
sigma[3,1] <- 1/3
sigma[3,2] <- 11/15
pmvnorm(lower=rep(-Inf, m), upper=c(1,4,2), mean=rep(0, m), corr=sigma)

```

--------------------------------

### Algorithm Configuration

Source: https://cran.r-project.org/web/packages/mvtnorm/refman/mvtnorm.html

Endpoints for configuring the hyper-parameters of the integration algorithms used for multivariate distribution evaluation.

```APIDOC
## GenzBretz

### Description
Configures the randomized Quasi-Monte-Carlo procedure for multivariate normal and t-distribution probabilities.

### Parameters
- **maxpts** (integer) - Optional - Maximum number of function values.
- **abseps** (double) - Optional - Absolute error tolerance.
- **releps** (double) - Optional - Relative error tolerance.

## Miwa

### Description
Configures the Miwa, Hayter, and Kuriki algorithm for normal probabilities in smaller dimensions.

### Parameters
- **steps** (integer) - Optional - Number of grid points to be evaluated (max 4097).
- **checkCorr** (logical) - Optional - Whether to check for singularity of the correlation matrix.
- **maxval** (double) - Optional - Replacement for Inf in non-orthant probability computations.

## TVPACK

### Description
Configures the TVPACK interface for two- and three-dimensional problems.

### Parameters
- **abseps** (double) - Optional - Absolute error tolerance.
```

--------------------------------

### Calculate Likelihood Using Inverse Cholesky Factor

Source: https://cran.r-project.org/web/packages/mvtnorm/vignettes/lmvnorm_src.Rnw

Calculates the log-likelihood for multivariate normal probabilities using the inverse of the Cholesky factor. Requires prior setup of matrices 'a', 'b', and 'W'.

```R
mL <- solve(mC)

lliL <- c(lpmvnorm(a, b, invchol = mL, w = W, M = M, logLik = FALSE))
```

```R
chk(lli, lliL)
```

```R
fL <- function(prm) {
    L <- ltMatrices(matrix(prm, ncol = 1), diag = TRUE)
    lpmvnorm(a, b, invchol = L, w = W, M = M)
}
```

```R
sL <- slpmvnorm(a, b, invchol = mL, w = W, M = M)
```

```R
chk(lliL, sL$logLik)
```

```R
if (require("numDeriv", quietly = TRUE)) 
    chk(grad(fL, unclass(mL)), rowSums(unclass(sL$invchol)),
        check.attributes = FALSE)
```

--------------------------------

### Initialize Cholesky or Inverse Cholesky Factors

Source: https://cran.r-project.org/web/packages/mvtnorm/vignettes/lmvnorm_src.Rnw

Handles the conversion of input matrices to Cholesky factors and ensures that only one of chol or invchol is provided.

```R
if (missing(chol) && missing(invchol))
    chol <- as.chol(ltMatrices(1, diag = TRUE))
stopifnot(xor(missing(chol), missing(invchol)))

if (!missing(chol)) {
    if (!is.ltMatrices(chol))
        chol <- as.ltMatrices(chol)
    scale <- as.chol(chol)
}

if (!missing(invchol)) {
    if (!is.ltMatrices(invchol))
        invchol <- as.ltMatrices(invchol)
    scale <- as.invchol(invchol)
}
ret <- list(scale = scale)
```

--------------------------------

### Define Cholesky Factor and Calculate Likelihood

Source: https://cran.r-project.org/web/packages/mvtnorm/vignettes/lmvnorm_src.Rnw

Defines the Cholesky factor of a covariance matrix and calculates the log-likelihood for multivariate normal probabilities with interval censoring. Requires prior setup of matrices 'a', 'b', and 'W'.

```R
mC <- ltMatrices(matrix(prm, ncol = 1), diag = TRUE)

a <- matrix(runif(N * J), nrow = J) - 2
b <- a + 4
a[2,] <- -Inf
b[3,] <- Inf

M <- 10000L
W <- matrix(runif(M * (J - 1)), ncol = M)

lli <- c(lpmvnorm(a, b, chol = mC, w = W, M = M, logLik = FALSE))
```

```R
fC <- function(prm) {
    C <- ltMatrices(matrix(prm, ncol = 1), diag = TRUE)
    lpmvnorm(a, b, chol = C, w = W, M = M)
}
```

```R
sC <- slpmvnorm(a, b, chol = mC, w = W, M = M)
```

```R
chk(lli, sC$logLik)
```

```R
if (require("numDeriv", quietly = TRUE)) 
    chk(grad(fC, unclass(mC)), rowSums(unclass(sC$chol)), 
        check.attributes = FALSE)
```

--------------------------------

### Optimize MVN Parameters with Inverse Cholesky Factor

Source: https://cran.r-project.org/web/packages/mvtnorm/vignettes/lmvnorm_src.Rnw

Optimizes the multivariate normal distribution parameters using the inverse Cholesky factor parameterization. This approach is noted to be a convex problem, making starting values less critical.

```R
opL <- optim(start, fn = ll, gr = sc, method = "L-BFGS-B", 
            lower = llim, control = list(trace = FALSE, factr = 1e-6))
MLL <- ll(opL$par, logLik = FALSE)
```