\documentclass[a4paper]{report}
\usepackage{a4wide}

%%% DO NOT EDIT THIS FILE
%%% Edit 'lmvnorm_src.w' and run 'nuweb -r lmvnorm_src.w'

%% packages
\usepackage{amsfonts,amstext,amsmath,amssymb,amsthm,nicefrac}

%\VignetteIndexEntry{Multivariate Normal Log-likelihoods in the mvtnorm Package}
%\VignetteDepends{mvtnorm,qrng,numDeriv}
%\VignetteKeywords{multivariate normal distribution}
%\VignettePackage{mvtnorm}

%% hyperlinked dois in jss refstyle
\newcommand{\doi}[1]{\href{https://doi.org/#1}{\normalfont\texttt{\@@doi{#1}}}}


\usepackage[utf8]{inputenc}

\newif\ifshowcode
\showcodetrue

\usepackage{latexsym}
%\usepackage{html}

\usepackage{listings}

\usepackage{color}
\definecolor{linkcolor}{rgb}{0, 0, 0.7}

\usepackage[%
backref,%
raiselinks,%
pdfhighlight=/O,%
pagebackref,%
hyperfigures,%
breaklinks,%
colorlinks,%
pdfpagemode=None,%
pdfstartview=FitBH,%
linkcolor={linkcolor},%
anchorcolor={linkcolor},%
citecolor={linkcolor},%
filecolor={linkcolor},%
menucolor={linkcolor},%
pagecolor={linkcolor},%
urlcolor={linkcolor}%
]{hyperref}

\usepackage[round]{natbib}

%\setlength{\oddsidemargin}{0in}
%\setlength{\evensidemargin}{0in}
%\setlength{\topmargin}{0in}
%\addtolength{\topmargin}{-\headheight}
%\addtolength{\topmargin}{-\headsep}
%\setlength{\textheight}{8.9in}
%\setlength{\textwidth}{6.5in}
%\setlength{\marginparwidth}{0.5in}

\newcommand{\pkg}[1]{\textbf{#1}}
\newcommand{\proglang}[1]{\textsf{#1}}
\newcommand{\code}[1]{\texttt{#1}}
\newcommand{\cmd}[1]{\texttt{#1()}}

\newcommand{\R}{\mathbb{R} }
\newcommand{\Prob}{\mathbb{P} }
\newcommand{\ND}{\mathbb{N} }
\newcommand{\J}{J}
\newcommand{\K}{K}
\newcommand{\V}{\mathbb{V}} %% cal{\mbox{\textnormal{Var}}} }
\newcommand{\E}{\mathbb{E}} %%mathcal{\mbox{\textnormal{E}}} }
\newcommand{\yvec}{\mathbf{y}}
\newcommand{\wvec}{\mathbf{w}}
\newcommand{\avec}{\mathbf{a}}
\newcommand{\bvec}{\mathbf{b}}
\newcommand{\xvec}{\mathbf{x}}
\newcommand{\svec}{\mathbf{s}}
\newcommand{\jvec}{\mathbf{j}}
\newcommand{\muvec}{\boldsymbol{\mu}}
\newcommand{\etavec}{\boldsymbol{\eta}}
\newcommand{\nuvec}{\boldsymbol{\nu}}
\newcommand{\lambdavec}{\boldsymbol{\lambda}}
\newcommand{\betavec}{\boldsymbol{\beta}}
\newcommand{\rY}{\mathbf{Y}}
\newcommand{\rX}{\mathbf{X}}
\newcommand{\rZ}{\mathbf{Z}}
\newcommand{\mC}{\mathbf{C}}
\newcommand{\mL}{\mathbf{L}}
\newcommand{\mP}{\mathbf{P}}
\newcommand{\mR}{\mathbf{R}}
\newcommand{\mT}{\mathbf{T}}
\newcommand{\mB}{\mathbf{B}}
\newcommand{\mI}{\mathbf{I}}
\newcommand{\mS}{\mathbf{S}}
\newcommand{\mA}{\mathbf{A}}
\newcommand{\mX}{\mathbf{X}}
\newcommand{\mD}{\mathbf{D}}
\newcommand{\diag}{\text{diag}}
\newcommand{\mSigma}{\mathbf{\Sigma}}
\newcommand{\argmin}{\operatorname{argmin}\displaylimits}
\newcommand{\argmax}{\operatorname{argmax}\displaylimits}
\newcommand{\vecop}{\text{vec}}


<<mvtnorm-citation, echo = FALSE>>=
year <- substr(packageDescription("mvtnorm")$Date, 1, 4)
version <- packageDescription("mvtnorm")$Version
@@

<<digits, echo = FALSE>>=
options(digits = 4)
@@

\author{Torsten Hothorn}

\date{Version \Sexpr{version}}

\title{Multivariate Normal Log-likelihoods \\ in the \pkg{mvtnorm} Package
\footnote{Please cite this document as: Torsten Hothorn (\Sexpr{year})
Multivariate Normal Log-likelihoods in the \pkg{mvtnorm} Package.
\textsf{R} package vignette version \Sexpr{version}, 
URL \href{https://doi.org/10.32614/CRAN.package.mvtnorm}{DOI:10.32614/CRAN.package.mvtnorm}.}
}

\begin{document}

\pagenumbering{roman}
\maketitle


\tableofcontents


\chapter*{Licence}

{\setlength{\parindent}{0cm}
Copyright (C) 2022-- Torsten Hothorn \\

This file is part of the \pkg{mvtnorm} \proglang{R} add-on package. \\

\pkg{mvtnorm} is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, version 2. \\

\pkg{mvtnorm} is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details. \\

You should have received a copy of the GNU General Public License
along with \pkg{mvtnorm}.  If not, see <http://www.gnu.org/licenses/>.
}

\chapter{Introduction}
\pagenumbering{arabic}

This document describes an implementation of \cite{numerical-:1992} and,
partially, of \cite{Genz_Bretz_2002}, for the  evaluation of
$N$ multivariate $\J$-dimensional normal probabilities
\begin{eqnarray} \label{pmvnorm}
p_i(\mC_i \mid \avec_i, \bvec_i) = \Prob(\avec_i < \rY_i \le \bvec_i \mid \mC_i ) 
  = (2 \pi)^{-\frac{\J}{2}} \text{det}(\mC_i)^{-1} 
    \int_{\avec_i}^{\bvec_i} \exp\left(-\frac{1}{2} \yvec^\top \mC_i^{-\top} \mC_i^{-1} \yvec\right) \, d \yvec
\end{eqnarray}
where $\avec_i = (a^{(i)}_1, \dots, a^{(i)}_\J)^\top \in \R^\J$ and 
$\bvec_i = (b^{(i)}_1, \dots, b^{(i)}_\J)^\top \in \R^\J$ are integration
limits, $\mC_i = (c^{(i)}_{j\jmath}) \in \R^{\J \times
\J}$ is a lower triangular matrix with $c^{(i)}_{j \jmath} = 0$ for $1 \le
j < \jmath < \J$, and thus $\rY_i \sim \ND_\J(\mathbf{0}_\J, \mC_i \mC_i^\top)$ for $i = 1, \dots, N$.

One application of these integrals is the estimation of the Cholesky factor
$\mC$ of a $\J$-dimensional normal distribution based on $N$ interval-censored
observations $\rY_1, \dots, \rY_\J$ (encoded by $\avec$ and $\bvec$) via maximum-likelihood
\begin{eqnarray*}
\hat{\mC} = \argmax_\mC \sum_{i = 1}^N \log(p_i(\mC \mid \avec_i, \bvec_i)).
\end{eqnarray*}
In other applications, the Cholesky factor might also depend on $i$ in some
structured way.

Function \code{pmvnorm} in package \code{mvtnorm} computes $p_i$ based on
the covariance matrix $\mC_i \mC_i^\top$. However, the Cholesky factor $\mC_i$
of the given covariance matrix is computed in \proglang{FORTRAN} first each
time this function is called. Function \code{pmvnorm} is not vectorised
over $i = 1, \dots, N$ and thus separate calls to this function are
necessary in order to compute likelihood contributions.

The implementation described here is a re-implementation (in \proglang{R}
and \proglang{C}) of Alan Genz' original \proglang{FORTRAN} code, focusing 
on efficient computation of the log-likelihood $\sum_{i = 1}^N \log(p_i)$
and the corresponding score function.

The document first describes a class and some useful methods for dealing
with multiple lower triangular matrices $\mC_i, i = 1, \dots, N$ in
Chapter~\ref{ltMatrices}.  The multivariate normal log-likelihood, and the
corresponding score function, is implemented as outlined in
Chapter~\ref{lpmvnorm}.  An example demonstrating maximum-likelihood
estimation of Cholesky factors in the presence of interval-censored
observations is discussed in Chapter~\ref{ML}.  We use the technology
developed here to implement the log-likelihood and score function for
situations where some variables have been observed exactly and others only
in form of interval-censoring in Chapter~\ref{cdl} and for nonparametric
maximum-likelihood estimation in unstructured Gaussian copulae in
Chapter~\ref{copula}. An attempt to provide useRs with a simple and
(hopefully) bullet proof interface is documented in Chapter~\ref{inter}.

The development of this infrastructure was motivated by the necessity to
evaluate probabilities~(\ref{pmvnorm}) arising in the likelihood of
multivariate conditional transformation models \citep{Klein_Hothorn_Barbanti_2020} for
discrete or censored observations. Some forms of the likelihood for such
nonparanormal models are discussed in \cite{Hothorn_2024}.

\chapter{Lower Triangular Matrices} \label{ltMatrices}

@o ltMatrices.R -cp
@{
@<R Header@>
@<.lt@>
@<ltMatrices@>
@<syMatrices@>
@<dim ltMatrices@>
@<dimnames ltMatrices@>
@<names ltMatrices@>
@<is.ltMatrices@>
@<as.ltMatrices@>
@<print ltMatrices@>
@<reorder ltMatrices@>
@<subset ltMatrices@>
@<lower triangular elements@>
@<diagonals ltMatrices@>
@<diagonal matrix@>
@<mult ltMatrices@>
@<mult syMatrices@>
@<solve ltMatrices@>
@<logdet ltMatrices@>
@<tcrossprod ltMatrices@>
@<crossprod ltMatrices@>
@<crossprod tcrossprod methods@>
@<chol syMatrices@>
@<invchol syMatrices@>
@<add diagonal elements@>
@<assign diagonal elements@>
@<kronecker vec trick@>
@<convenience functions@>
@<aperm@>
@<marginal@>
@<conditional@>
@<check obs@>
@<colSumsdnorm ltMatrices@>
@}

@o ltMatrices.c -cc
@{
@<C Header@>
#ifndef  USE_FC_LEN_T
# define USE_FC_LEN_T
#endif
#include <Rconfig.h>
#include <R_ext/Lapack.h> /* for dtptri */
#ifndef FCONE
# define FCONE
#endif
#include <R.h>
#include <Rmath.h>
#include <Rinternals.h>
#include <Rdefines.h>
@<colSumsdnorm@>
@<solve@>
@<solve C@>
@<logdet@>
@<mult@>
@<chol@>
@<invchol workhorse@>
@<invchol@>
@}


We first define and implement infrastructure for dealing with multiple lower triangular matrices
$\mC_i \in \R^{\J \times \J}$ for $i = 1, \dots, N$. We note that each such matrix
$\mC$ can be stored in a vector of length $\J (\J + 1) / 2$. If all
diagonal elements are one (that is, $c^{(i)}_{jj} \equiv 1, j = 1, \dots,
\J$), the length of this vector is $\J (\J - 1) / 2$.

\section{Multiple Lower Triangular Matrices}

We can store $N$ such matrices in an $\J (\J + 1) / 2 \times N$ matrix
(\code{diag = TRUE}) or, for \code{diag = FALSE}, in an $\J (\J
- 1) / 2 \times N$ matrix.

Each vector might define the corresponding lower triangular matrix
either in row or column-major order:

\begin{eqnarray*}
 \mC & = & \begin{pmatrix}
 c_{11} & & & & 0\\
 c_{21} & c_{22} \\
 c_{31} & c_{32} & c_{33} \\
 \vdots & \vdots & & \ddots & \\
 c_{J1} & c_{J2} & \ldots & &  c_{JJ}
 \end{pmatrix}  \text{matrix indexing}\\
& = &  
\begin{pmatrix}
 c_{1} & & & & 0\\
 c_{2} & c_{J + 1} \\
 c_{3} & c_{J + 2} & c_{2J} \\
 \vdots & \vdots & & \ddots & \\
 c_{J} & c_{2J - 1} & \ldots & &  c_{J(J + 1) / 2}
 \end{pmatrix} \text{column-major, \code{byrow = FALSE}} \\
& = & \begin{pmatrix}
 c_{1} & & & & 0\\
 c_{2} & c_{3} \\
 c_{4} & c_{5} & c_{6} \\
 \vdots & \vdots & & \ddots & \\
 c_{J((J + 1) / 2 -1) + 1} & c_{J((J + 1) / 2 -1) + 2} & \ldots & &  c_{J(J + 1) / 2}
 \end{pmatrix} \text{row-major, \code{byrow = TRUE}}
\end{eqnarray*}

Based on some matrix \code{object}, the dimension $\J$ is computed and checked as
@d ltMatrices dim
@{
J <- floor((1 + sqrt(1 + 4 * 2 * nrow(object))) / 2 - diag)
if (nrow(object) != J * (J - 1) / 2 + diag * J)
    stop("Dimension of object does not correspond to lower 
          triangular part of a square matrix")
@}

Before we begin assigning names, we set-up an internal shortcut for
\code{lower.tri}, which is a bit more memory and speed efficient

@d .lt
@{
.lt <- function(J, diag = FALSE) {
    j <- seq_len(J)
    j1 <- rep(c(TRUE, FALSE), J)[- 2 * J]
    j2 <- c(rbind(rev(j), j))[- 2 * J]
    ret <- matrix(rep(j1, j2), nrow = J)
    if (!diag) diag(ret) <- FALSE
    ret
}
.ut <- function(J, diag = FALSE) {
    !.lt(J = J, diag = !diag)
}
@}

Typically the $\J$ dimensions are associated with names, and we therefore
compute identifiers for the vector elements in either column- or row-major
order on request (for later printing)

@d ltMatrices names
@{
nonames <- FALSE
if (!isTRUE(names)) {
    if (is.character(names))
        stopifnot(is.character(names) &&
                  length(unique(names)) == J)
    else
        nonames <- TRUE
} else {
    names <- as.character(seq_len(J))
}

if (!nonames) {
    idx <- matrix(1:J, nrow = J, ncol = J)
    idx <- idx[.lt(J, diag = diag)]
    j <- seq_len(J - !(diag + 0L))
    idx2 <- rep(j, rev(j))
    rn <- paste(names[idx], names[idx2], sep = ".")
    if (byrow)
       rn <- rn[.c2r(J, diag = diag)]
    rownames(object) <- rn
} else {	### add later
    warning("ltMatrices objects should be properly named")
}
@}

If \code{object} is already a classed object representing lower triangular
matrices (we will use the class name \code{ltMatrices}), we might want to
change the storage form from row- to column-major or the other way round.

@d ltMatrices input
@{
if (is.ltMatrices(object)) {
    cls <- class(object)		### keep inheriting classes
    ret <- .reorder(object, byrow = byrow)
    class(ret) <- class(object)
    return(ret)
}
@}

The constructor essentially attaches attributes to a matrix \code{object},
possibly after some reordering / transposing

@d ltMatrices
@{
ltMatrices <- function(object, diag = FALSE, byrow = FALSE, names = TRUE) {

    if (!is.matrix(object)) 
        object <- matrix(object, ncol = 1L)

    @<ltMatrices input@>
    @<ltMatrices dim@>
    @<ltMatrices names@>

    attr(object, "J")       <- J
    attr(object, "diag")    <- diag
    attr(object, "byrow")   <- byrow
    attr(object, "rcnames") <- names

    class(object) <- c("ltMatrices", class(object))
    object
}
@}

For the sake of completeness, we also add a constructor for
multiple symmetric matrices

@d syMatrices
@{
as.syMatrices <- function(x) {
    if (is.syMatrices(x))
        return(x)
    if (is.ltMatrices(x)) {
        class(x) <- gsub("ltMatrices", "syMatrices", class(x))
        return(x)
    }
    as.syMatrices(as.ltMatrices(x * .lt(nrow(x), diag = TRUE)))
}
syMatrices <- function(object, diag = FALSE, byrow = FALSE, names = TRUE) {
    if (inherits(object, "syMatrices"))
        class(object) <- gsub("syMatrices", "ltMatrices", class(object))
    as.syMatrices(ltMatrices(object = object, diag = diag, byrow = byrow, 
                             names = names))
}
@}

The dimensions of such an object are always $N \times \J \times \J$ and are given by

@d dim ltMatrices
@{
dim.ltMatrices <- function(x) {
    J <- attr(x, "J")
    return(c(attr(x, "dim")[2L], J, J)) ### ncol(unclass(x)) may trigger gc
}
dim.syMatrices <- dim.ltMatrices
@}

The corresponding dimnames can be extracted as

@d dimnames ltMatrices
@{
dimnames.ltMatrices <- function(x)
    return(list(attr(x, "dimnames")[[2L]], attr(x, "rcnames"), attr(x, "rcnames")))
dimnames.syMatrices <- dimnames.ltMatrices
@}

The names identifying rows and columns in each $\mC_i$ are

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

Finally, let's add two functions for checking the class and a function for
coersing classes inheriting from \code{ltMatrices} to the latter, the same
for \code{syMatrices}. Furthermode, \code{as.ltMatrices} coerces objects
inheriting from \code{syMatrices} or \code{ltMatrices} to class
\code{ltMatrices} (that is, \code{chol} or \code{invchol} is removed from
the class list, unlike a call to the constructor \code{ltMatrices}). A
\code{default} method is added in Chapter~\ref{inter}.

@d is.ltMatrices
@{
is.ltMatrices <- function(x) inherits(x, "ltMatrices")
is.syMatrices <- function(x) inherits(x, "syMatrices")
as.ltMatrices <- function(x) UseMethod("as.ltMatrices")
as.ltMatrices.syMatrices <- function(x) {
    cls <- class(x)
    class(x) <- cls[which(cls == "syMatrices"):length(cls)]
    class(x)[1L] <- "ltMatrices"
    return(x)
}
as.ltMatrices.ltMatrices <- function(x) {
    cls <- class(x)
    class(x) <- cls[which(cls == "ltMatrices"):length(cls)]
    return(x)
}
@}

Let's set-up an example for illustration. Throughout this document, we will
compare numerical results using
<<chk>>=
chk <- function(...) stopifnot(isTRUE(all.equal(...)))
@@
We start with a a simple example demonstrating how to set-up
\code{ltMatrices} objects

<<example>>=
library("mvtnorm")
set.seed(290875)
N <- 4L
J <- 5L
rn <- paste0("C_", 1:N)
nm <- LETTERS[1:J]
Jn <- J * (J - 1) / 2
## data
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
@@



\section{Printing}

For pretty printing, we coerse objects of class \code{ltMatrices} to
\code{array}. The method has a logical argument called \code{symmetric}, forcing the lower
triangular matrix to by interpreted as a symmetric matrix.

@d extract slots
@{
diag <- attr(x, "diag")
byrow <- attr(x, "byrow")
d <- dim(x)
N <- d[1L]
J <- d[2L]
dn <- dimnames(x)
@}

@d print ltMatrices
@{
as.array.ltMatrices <- function(x, symmetric = FALSE, ...) {

    @<extract slots@>

    x <- unclass(x)

    L <- matrix(1L, nrow = J, ncol = J)
    diag(L) <- 2L
    if (byrow) {
        L[.ut(J, diag = diag)] <- floor(2L + 1:(J * (J - 1) / 2L + diag * J))
        L <- t(L)
    } else {
        L[.lt(J, diag = diag)] <- floor(2L + 1:(J * (J - 1) / 2L + diag * J))
    }
    if (symmetric) {
        L[.ut(J)] <- 0L
        dg <- diag(L)
        L <- L + t(L)
        diag(L) <- dg
    }
    ret <- rbind(0, 1, x)[c(L), , drop = FALSE]
    class(ret) <- "array"
    dim(ret) <- d[3:1]
    if (!isFALSE(dn[[2L]]) && !isFALSE(dn[[3L]]))
        dimnames(ret) <- dn[3:1]
    return(ret)
}

as.array.syMatrices <- function(x, ...)
    return(as.array.ltMatrices(x, symmetric = TRUE))

print.ltMatrices <- function(x, zero.print = ".", ...) {
    if (is.null(attr(x, "dimnames")[[2L]]))
        attr(x, "dimnames")[[2L]] <- as.character(seq_len(dim(x)[1L]))
    print(as.table(as.array(x)), zero.print = zero.print, ...)
}

print.syMatrices <- function(x, ...)
    print(as.array(x))
@}

Symmetric matrices are represented by lower triangular matrix objects, but
we change the class from \code{ltMatrices} to \code{syMatrices} (which
disables all functionality except printing and coersion to arrays).

\section{Reordering}

It is sometimes convenient to have access to lower triangular matrices in
either column- or row-major order and this little helper function switches
between the two forms. This is quite fast but we still try to avoid it, for
example in \code{Mult}, Section~\ref{sec:multiplication}.

@d reorder ltMatrices
@{
.r2c <- function(J, diag = FALSE) {

    rL <- cL <- diag(0L, nrow = J)
    lt <- .lt(J, diag = diag)
    ut <- !lt
    diag(ut) <- !diag(ut)
    x <- seq_len(J * (J + c(-1, 1)[diag + 1L]) / 2)
    cL[ut] <- seq_len(length(x))
    x[t(cL)[lt]]
}

.c2r <- function(J, diag = FALSE) {

    rL <- cL <- diag(0L, nrow = J)
    lt <- .lt(J, diag = diag)
    ut <- !lt
    diag(ut) <- !diag(ut)
    x <- seq_len(J * (J + c(-1, 1)[diag + 1L]) / 2)
    rL[lt] <- seq_len(length(x))
    x[t(rL)[ut]]
}

.reorder <- function(x, byrow = FALSE) {

    ### only we call this function
    ### stopifnot(is.ltMatrices(x))
    if (attr(x, "byrow") == byrow) return(x)

    @<extract slots@>

    x <- unclass(x)

    if (byrow) { ### row -> col order
        return(ltMatrices(x[.r2c(J, diag = diag), , drop = FALSE], 
                          diag = diag, byrow = FALSE, names = dn[[2L]]))
    }
    return(ltMatrices(x[.c2r(J, diag = diag), , drop = FALSE], 
                      diag = diag, byrow = TRUE, names = dn[[2L]]))
}
@}

We can check if this works by switching back and forth between column-major
and row-major order

<<ex-reorder>>=
## constructor + .reorder + as.array
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)
@@

\section{Subsetting}

We might want to select subsets of observations $i \in \{1, \dots, N\}$ or
rows/columns $j \in \{1, \dots, \J\}$ of the corresponding matrices $\mC_i$. 

@d .subset ltMatrices
@{
.subset_ltMatrices <- function(x, i, j, ..., drop = FALSE) {

    if (drop) warning("argument drop is ignored")
    if (missing(i) && missing(j)) return(x)

    @<extract slots@>

    x <- unclass(x) 

    if (!missing(j)) {

        if (is.character(j)) {
            stopifnot(all(j %in% dn[[2L]]))
            j <- match(j, dn[[2L]])
        }

        j <- (1:J)[j] ### get rid of negative indices

        if (length(j) == 1L && !diag) {
            return(ltMatrices(matrix(1, ncol = ncol(x), nrow = 1), diag = TRUE, 
                              byrow = byrow, names = dn[[2L]][j]))
        }
        L <- diag(0L, nrow = J)
        Jp <- sum(.ut(J, diag = diag))
        if (byrow) {
            L[.ut(J, diag = diag)] <- 1:Jp
            L <- L + t(L)
            diag(L) <- diag(L) / 2
            L <- L[j, j, drop = FALSE]
            L <- L[.ut(nrow(L), diag = diag)]
        } else {
            L[.lt(J, diag = diag)] <- 1:Jp
            L <- L + t(L)
            diag(L) <- diag(L) / 2
            L <- L[j, j, drop = FALSE]
            L <- L[.lt(nrow(L), diag = diag)]
        }
        if (missing(i)) {
            return(ltMatrices(x[c(L), , drop = FALSE], diag = diag, 
                              byrow = byrow, names = dn[[2L]][j]))
        }
        return(ltMatrices(x[c(L), i, drop = FALSE], diag = diag, 
                          byrow = byrow, names = dn[[2L]][j]))
    }
    return(ltMatrices(x[, i, drop = FALSE], diag = diag, 
                      byrow = byrow, names = dn[[2L]]))
}
@}

@d subset ltMatrices
@{
@<.subset ltMatrices@>
### if j is not ordered, result is not a lower triangular matrix
"[.ltMatrices" <- function(x, i, j, ..., drop = FALSE) {
    if (!missing(j)) {
        if (is.character(j)) {
            stopifnot(all(j %in% dimnames(x)[[2L]]))
            j <- match(j, dimnames(x)[[2L]])
        }
        if (all(j > 0)) {
            if (any(diff(j) < 0)) stop("invalid subset argument j")
        }
    }

    return(.subset_ltMatrices(x = x, i = i, j = j, ..., drop = drop))
}

"[.syMatrices" <- function(x, i, j, ..., drop = FALSE) {
    x <- as.syMatrices(x)
    ret <- .subset_ltMatrices(x = x, i = i, j = j, ..., drop = drop)
    class(ret)[1L] <- "syMatrices"
    ret
}
@}

We check if this works by first subsetting the \code{ltMatrices} object.
Second, we coerse the object to an array and do the subset for the latter
object. Both results must agree.

<<ex-subset, eval = FALSE>>=
## subset
a <- as.array(ltMatrices(xn, byrow = FALSE, names = nm)[i, j])
b <- as.array(ltMatrices(xn, byrow = FALSE, names = nm))[j, j, i]
chk(a, b)

a <- as.array(ltMatrices(xn, byrow = TRUE, names = nm)[i, j])
b <- as.array(ltMatrices(xn, byrow = TRUE, names = nm))[j, j, i]
chk(a, b)

a <- as.array(ltMatrices(xd, byrow = FALSE, 
                         diag = TRUE, names = nm)[i, j])
b <- as.array(ltMatrices(xd, byrow = FALSE, 
                         diag = TRUE, names = nm))[j, j, i]
chk(a, b)

a <- as.array(ltMatrices(xd, byrow = TRUE, diag = TRUE, 
                         names = nm)[i, j])
b <- as.array(ltMatrices(xd, byrow = TRUE, diag = TRUE, 
                         names = nm))[j, j, i]
chk(a, b)
@@

We start with both indices being positive integers

<<ex-subset-1>>=
i <- colnames(xn)[1:2]
j <- 2:4
<<ex-subset>>
@@

proceed with characters

<<ex-subset-2>>=
i <- 1:2
j <- nm[2:4]
<<ex-subset>>
@@

a different subset

<<ex-subset-3>>=
j <- c(1, 3, 5)
<<ex-subset>>
@@

and characters again

<<ex-subset-4>>=
j <- nm[c(1, 3, 5)]
<<ex-subset>>
@@

and finally with with negative subsets

<<ex-subset-5>>=
j <- -c(1, 3, 5)
<<ex-subset>>
@@

and with non-increasing argument \code{j} (this won't work for lower
triangular matrices, only for symmetric matrices)

<<ex-subset-6>>=
## subset
j <- nm[sample(1:J)]
ltM <- ltMatrices(xn, byrow = FALSE, names = nm)
print(try(ltM[i, j], silent = TRUE))
ltM <- as.syMatrices(ltM)
a <- as.array(ltM[i, j])
b <- as.array(ltM)[j, j, i]
chk(a, b)
@@

Extracting the lower triangular elements from an \code{ltMatrices} object
(or from an object of class \code{syMatrices}) returns a matrix with $N$
columns, undoing the effect of \code{ltMatrices}. Note that ordering of the
rows of this matrix depend on the \code{byrow} attribute of \code{x}, unless
the \code{byrow} to this function is used to overwrite it explicitly

@d lower triangular elements
@{
Lower_tri <- function(x, diag = FALSE, byrow = attr(x, "byrow")) {

    if (is.syMatrices(x))
        x <- as.ltMatrices(x)
    adiag <- diag
    x <- ltMatrices(x, byrow = byrow)

    @<extract slots@>

    if (diag == adiag)
        return(unclass(x)[,,drop = FALSE]) ### remove attributes

    if (!diag && adiag) {
        diagonals(x) <- 1
        return(unclass(x)[,,drop = FALSE]) ### remove attributes
    }

    x <- unclass(x)
    if (J == 1) {
        idx <- 1L
    } else {
      if (byrow)
          idx <- cumsum(c(1, 2:J))
      else
          idx <- cumsum(c(1, J:2))
    }
    return(x[-idx,,drop = FALSE])
}
@}

<<ex-Lower_tri>>=
## 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)
## multiple symmetric matrices
Lower_tri(invchol2cor(M))
@@

\section{Diagonal Elements}

The diagonal elements of each matrix $\mC_i$ can be extracted and are
always returned as an $\J \times N$ matrix.

@d diagonals ltMatrices
@{
diagonals <- function(x, ...)
    UseMethod("diagonals")

diagonals.ltMatrices <- function(x, ...) {

    @<extract slots@>

    x <- unclass(x)

    if (!diag) {
        ret <- matrix(1, nrow = J, ncol = ncol(x))
        colnames(ret) <- dn[[1L]]
        rownames(ret) <- dn[[2L]]
        return(ret)
    } else {
        if (J == 1L) return(x)
        if (byrow)
            idx <- cumsum(c(1, 2:J))
        else
            idx <- cumsum(c(1, J:2))
        ret <- x[idx, , drop = FALSE]
        rownames(ret) <- dn[[2L]]
        return(ret)
    }
}

diagonals.syMatrices <- diagonals.ltMatrices

diagonals.matrix <- function(x, ...) diag(x)
@}

<<ex-diag>>=
all(diagonals(ltMatrices(xn, byrow = TRUE)) == 1L)
@@

Sometimes we need to add diagonal elements to an \code{ltMatrices} object
which was set-up with constant $c_{jj} = 1$ diagonal elements.

@d add diagonal elements
@{
.adddiag <- function(x) {

    ### only we call this function
    # stopifnot(is.ltMatrices(x))

    if (attr(x, "diag")) return(x)

    byrow_orig <- attr(x, "byrow")

    x <- ltMatrices(x, byrow = FALSE)

    N <- dim(x)[1L]
    J <- dim(x)[2L]
    nm <- dimnames(x)[[2L]]

    L <- D <- diag(J)
    lt <- .lt(J, diag = TRUE)
    L[lt] <- seq_len(J * (J + 1) / 2)

    ### by default new diagonal elements are 1
    ret <- matrix(D[lt], nrow = J * (J + 1) / 2, ncol = N)
    colnames(ret) <- dimnames(x)[[1L]]
    diag(lt) <- FALSE
    ret[L[lt],] <- unclass(x)

    ret <- ltMatrices(ret, diag = TRUE, byrow = FALSE, names = nm)
    ret <- ltMatrices(ret, byrow = byrow_orig)

    ret
}
@}


@d assign diagonal elements
@{
"diagonals<-" <- function(x, value)
    UseMethod("diagonals<-")

"diagonals<-.ltMatrices" <- function(x, value) {

    @<extract slots@>

    if (byrow)
        idx <- cumsum(c(1, 2:J))
    else
        idx <- cumsum(c(1, J:2))

    ### diagonals(x) <- NULL returns ltMatrices(..., diag = FALSE)
    if (is.null(value)) {
        if (!attr(x, "diag")) return(x)
        if (J == 1L) {
            x[] <- 1
            return(x)
        }
        return(ltMatrices(unclass(x)[-idx,,drop = FALSE], diag = FALSE, 
                          byrow = byrow, names = dn[[2L]]))
    }

    x <- .adddiag(x)

    if (!is.matrix(value))
        value <- matrix(value, nrow = J, ncol = d[1L])

    stopifnot(is.matrix(value) && nrow(value) == J 
                               && ncol(value) == d[1L])

    if (J == 1L) {
        x[] <- value
        return(x)
    }

    x[idx, ] <- value

    return(x)
}

"diagonals<-.syMatrices" <- function(x, value) {

    x <- as.ltMatrices(x)
    diagonals(x) <- value
    class(x)[1L] <- "syMatrices"

    return(x)
}
@}

<<ex-addiag>>=
lxd2 <- lxn
diagonals(lxd2) <- 1
chk(as.array(lxd2), as.array(lxn))
@@

A unit diagonal matrix is not treated as a special case but as an
\code{ltMatrices} object with all lower triangular elements being zero

@d diagonal matrix
@{
diagonals.integer <- function(x, ...)
    ltMatrices(rep(0, x * (x - 1) / 2), diag = FALSE, ...)
@}

<<ex-diagJ>>=
(I5 <- diagonals(5L))
diagonals(I5) <- 1:5
I5
@@


\section{Multiplication} \label{sec:multiplication}

Products $\mC_i \yvec_i$ or $\mC^\top_i \yvec_i$ with $\yvec_i \in \R^\J$ for $i = 1, \dots,
N$ can be computed with $\code{y}$ being an $J \times N$ matrix of
columns-wise stacked vectors $(\yvec_1 \mid \yvec_2 \mid \dots \mid
\yvec_N)$. If \code{y} is a single vector, it is recycled $N$ times.

If the number of columns of the matrix \code{y} is neither one nor $N$, 
we compute $\mC_i \yvec_j$ for all $i = 1, \dots, N$ and $j$. This is
dangerous but needed in Section~\ref{sec:margcond} for defining \code{cond\_mvnorm} later on.

We have three use cases for this matrix-vector multiplication: (1)
If the number of rows of the matrix \code{y} is $\J$, each column is
interpreted as $\yvec_i \in \R^\J$. (2) If the number of rows is $\J^2$, the
returned matrix is $\J^2 \times N$ with each columns representing
$\text{vec}(\mC_i \rY_i)$ for the $\J \times \J$ matrix $\rY_i = (\yvec_1 \mid
\yvec_2 \mid \dots \mid \yvec_\J)$. (3) If the number of rows of the matrix \code{y} is $\J
(\J + 1) / 2$, we interpret the $i$th column of \code{y} as the lower triangular 
elements of $\mC_i$ (in column-major storage) for computing $\mC_i
\mC_i^\top$ (\code{transpose = FALSE}) or $\mC_i^\top \mC_i$
(\code{transpose = TRUE}) later on in Section~\ref{sec:crossproducts}.
Again, this is dangerous but helps to speed-up crossproducts considerably.

For $\mC_i \yvec_i$, we call \proglang{C} code computing the product
efficiently without copying data by leveraging the lower triangular structure of
\code{x}$=\mC_i$

@d mult ltMatrices
@{
### C %*% y
Mult <- function(x, y, ...)
    UseMethod("Mult")
Mult.default <- function(x, y = NULL, transpose = FALSE, ...) {
    if (!transpose) return(x %*% y)
    return(crossprod(x, y = NULL, ...))
}
Mult.ltMatrices <- function(x, y = NULL, transpose = FALSE, ...) {

    @<extract slots@>
    ### dtpmv requires diagonal elements being present
    if (!diag) diagonals(x) <- diagonals(x)

    ### for crossprod(x), use case (3)    
    if (is.null(y)) {
        stopifnot(!byrow)
        stopifnot(transpose)
        y <- Lower_tri(x, diag = TRUE) 
    } else { 
        stopifnot(NROW(y) == J || NROW(y) == J^2)
    }

    stopifnot(is.numeric(y))
    if (!is.matrix(y)) y <- matrix(y, nrow = J, ncol = N)
    N <- ifelse(d[1L] == 1, ncol(y), N)
    if (ncol(y) != N)
        return(sapply(1:ncol(y), function(i) Mult(x, y[,i], 
                                                  transpose = transpose)))

    if (!is.double(x)) storage.mode(x) <- "double"
    if (!is.double(y)) storage.mode(y) <- "double"

    ret <- .Call(mvtnorm_R_ltMatrices_Mult, x, y, as.integer(N), 
                 as.integer(J), as.logical(diag), as.logical(transpose))

    if (NROW(y) == d[2L] && !isFALSE(dn[[2L]]))    
        rownames(ret) <- dn[[2L]]
    if (length(dn[[1L]]) == N && !isFALSE(dn[[1L]]))
        colnames(ret) <- dn[[1L]]
    return(ret)
}
@}


We start with some options for the \proglang{LAPACK} workhorses, regarding
the presence of non-unit diagonal elements (which, however, have to be part
of $\mC$) and whether $\mC$ comes in row- or column-major order.

@d lapack options
@{
char di, uplo;
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';
}

SEXP abyrow = PROTECT(getAttrib(C, Rf_install("byrow")));
if (TYPEOF(abyrow) != LGLSXP) error("non-logical byrow attribute");
if (LOGICAL(abyrow)[0]) {
    uplo = 'U';
} else {
    uplo = 'L';
}

@}


The underlying \proglang{C} code allows $\mC_i$ (here called \code{C}) to
be in row-major or column-major order. In column-major $\mC \yvec$, we have
\code{uplo = "L"} (lower) and \code{tr = "U"} (take $\mC$ at face value). If
$\mC$ comes in row-major order, we pretent it is an upper triangular matrix
(\code{uplo = "U"}) and and compute $\mC^\top \yvec$ (\code{tr = "T"}).

@d row column
@{
char tr = 'N';
Rboolean Rtranspose = asLogical(transpose);
/* t(C) instead of C for column-major*/
if (Rtranspose && !LOGICAL(abyrow)[0])
    tr = 'T';
/* C instead for row-major*/
if (!Rtranspose && LOGICAL(abyrow)[0])
    tr = 'T';
@}

We extract dimensions next

@d RC input
@{
/* pointer to C matrices */
double *dC = REAL(C);
/* number of matrices */
int iN = INTEGER(N)[0];
/* dimension of matrices */
int iJ = INTEGER(J)[0];
/* C contains diagonal elements */
Rboolean Rdiag = asLogical(diag);
/* p = J * (J - 1) / 2 + diag * J */
int len = iJ * (iJ - 1) / 2 + Rdiag * iJ;
@}

We also allow $\mC_i$ to be constant ($N$ is then determined from
\code{ncol(y)}). The following fragment ensures that we only loop over
$\mC_i$ if \code{dim(x)[1L] > 1}

@d C length
@{
int p;
if (LENGTH(C) == len)
    /* C is constant for i = 1, ..., N */
    p = 0;
else 
    /* C contains C_1, ...., C_N */
    p = len;
@}

We can now implement the main part of the multiplication. We first
differentiate between the three use cases. For (1), we loop over the columns
of \code{C} and / or the columns of \code{y}. For (2), we additionally loop
over the $\J$ vectors of length $\J$ defining the columns of \code{y}. For
(3), we only loop over the vectors of lengths $\J, \J - 1, \dots, 1$ and
also corresponding sub-matrices of $\mC$ such that the lower triangular
elements of the symmetric crossproduct $\mC \mC^\top$ emerge.

@d mult workhorse
@{
/* use case (1) */
if (LENGTH(y) == iJ * iN) {
    nry = 1;
    lty = 0;
/* use case (2) */
} else if (LENGTH(y) == iJ * iJ * iN) {
    nry = iJ;
    lty = 0;
/* use case (3) */
} else if (LENGTH(y) == iJ * (iJ + 1) / 2 * iN) {
    if (!(Rtranspose && !LOGICAL(abyrow)[0]))
        error("matrix multiplication not implemented");
    nry = iJ;
    lty = 1;
} else {
    error("incorrect dimension of y argument");
}

thisJ = iJ;
for (i = 0; i < iN; i++) {
    incr = 0;
    for (j = 0; j < nry; j++) {       
        if (lty) thisJ = iJ - j;
        F77_CALL(dtpmv)(&uplo, &tr, &di, &thisJ,
                        dC + incr, dans, &iONE FCONE FCONE FCONE);
        dans += (lty ? iJ - j : iJ);
        if (lty) 
            incr += thisJ;
    }
    dC += p;
}
@}

We are now ready to implement the \proglang{C} workhorse, using
\proglang{LAPACK} function \code{dtrmv}.

@d mult
@{
SEXP R_ltMatrices_Mult (SEXP C, SEXP y, SEXP N, SEXP J, SEXP diag, SEXP transpose) {

    SEXP ans;
    double *dans;
    int i, j, nry, lty, thisJ, incr;
    int iONE = 1;

    @<RC input@>
    /* diagonal elements are always present */
    if (!Rdiag) len += iJ;
    @<C length@>
    @<lapack options@>
    @<row column@>

    PROTECT(ans = duplicate(y));
    dans = REAL(ans);

    @<mult workhorse@>

    UNPROTECT(2);
    return(ans);
}
@}

Some checks for $\mC_i \yvec_i$

<<ex-mult>>=
### matrix-vector
lxn <- ltMatrices(xn, byrow = TRUE)
lxd <- ltMatrices(xd, byrow = TRUE, diag = TRUE)
y <- matrix(runif(N * J), nrow = J)
a <- Mult(lxn, y)
A <- as.array(lxn)
b <- do.call("rbind", lapply(1:ncol(y), 
    function(i) t(A[,,i] %*% y[,i,drop = FALSE])))
chk(a, t(b), check.attributes = FALSE)

### matrix-matrix
Y <- matrix(runif(N * J^2), nrow = J^2)
a <- Mult(lxn, Y)
b <- do.call("rbind", lapply(1:ncol(y), 
    function(i) t(A[,,i] %*% matrix(Y[,i,drop = FALSE], nrow = J))))
chk(a, t(b), check.attributes = FALSE)

### matrix-vector
a <- Mult(lxd, y)
A <- as.array(lxd)
b <- do.call("rbind", lapply(1:ncol(y), 
    function(i) t(A[,,i] %*% y[,i,drop = FALSE])))
chk(a, t(b), check.attributes = FALSE)

### matrix-matrix
a <- Mult(lxd, Y)
b <- do.call("rbind", lapply(1:ncol(y), 
    function(i) t(A[,,i] %*% matrix(Y[,i,drop = FALSE], nrow = J))))
chk(a, t(b), check.attributes = FALSE)

### recycle C
chk(Mult(lxn[rep(1, N),], y), Mult(lxn[1,], y), check.attributes = FALSE)

### recycle y
chk(Mult(lxn, y[,1]), Mult(lxn, y[,rep(1, N)]))

### tcrossprod as multiplication
i <- sample(1:N)[1]
M <- t(as.array(lxn)[,,i])
a <- sapply(1:J, function(j) Mult(lxn[i,], M[,j,drop = FALSE]))
rownames(a) <- colnames(a) <- dimnames(lxn)[[2L]]
b <- as.array(Tcrossprod(lxn[i,]))[,,1]
chk(a, b, check.attributes = FALSE)
@@

and wrap-up with some tests for computing $\mC^\top_i \yvec_i$

<<ex-tmult>>=
a <- Mult(lxn, y, transpose = TRUE)
A <- as.array(lxn)
b <- do.call("rbind", lapply(1:ncol(y), 
    function(i) t(t(A[,,i]) %*% y[,i,drop = FALSE])))
chk(a, t(b), check.attributes = FALSE)

a <- Mult(lxd, y, transpose = TRUE)
A <- as.array(lxd)
b <- do.call("rbind", lapply(1:ncol(y), 
    function(i) t(t(A[,,i]) %*% y[,i,drop = FALSE])))
chk(a, t(b), check.attributes = FALSE)

### recycle C
chk(Mult(lxn[rep(1, N),], y, transpose = TRUE), 
    Mult(lxn[1,], y, transpose = TRUE), check.attributes = FALSE)

### recycle y
chk(Mult(lxn, y[,1], transpose = TRUE), 
    Mult(lxn, y[,rep(1, N)], transpose = TRUE))
@@

Now we can add a \code{Mult} method for multiple symmetric matrices, noting
that for a symmetric matrix $\mA = \mC + \mC^\top - \text{diag}(\mC)$ with lower triangular
part $\mC$ (including the diagonal) we can compute $\mA \yvec = \mC \yvec +
\mC^\top \yvec - \text{diag}(\mC) \yvec$ using \code{Mult} applied to the lower trianular part:

@d mult syMatrices
@{
Mult.syMatrices <- function(x, y, ...) {

    @<extract slots@>

    x <- as.ltMatrices(x)
    stopifnot(is.numeric(y))
    if (!is.matrix(y)) y <- matrix(y, nrow = J, ncol = N)
    N <- ifelse(d[1L] == 1, ncol(y), N)
    ### <NOTE> not yet implemented for J x J matrices in y, see
    ### Mult.ltMatrices </NOTE>
    stopifnot(nrow(y) == d[2L])
    stopifnot(ncol(y) == N)

    ret <- Mult(x, y) + Mult(x, y, transpose = TRUE) - y * c(diagonals(x))
    return(ret)
}
@}

<<ex-symult>>=
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, check.attributes = FALSE)
})

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

\section{Solving Linear Systems}

Computing $\mC_i^{-1}$ or solving $\mC_i \xvec_i = \yvec_i$ for $\xvec_i$ for
all $i = 1, \dots, N$ is another important task. We sometimes also need $\mC^\top_i \xvec_i =
\yvec_i$ triggered by \code{transpose = TRUE}.

\code{C} is $\mC_i, i = 1, \dots, N$ in column-major order
(matrix of dimension $\J (\J - 1) / 2 + \J \text{diag} \times N$), and
\code{y} is the $\J \times N$ matrix $(\yvec_1 \mid \yvec_2 \mid \dots \mid
\yvec_N)$. This function returns the $\J \times N$ matrix $(\xvec_1 \mid \xvec_2 \mid \dots \mid
\xvec_N)$ of solutions.

If \code{y} is not given, $\mC_i^{-1}$ is returned in 
the same order as the orginal matrix $\mC_i$. If
all $\mC_i$ have unit diagonals, so will $\mC_i^{-1}$.

We set-up a dedicated \proglang{C} function for computing $\mC_i \xvec_i = \yvec_i$

@d solve
@{
SEXP R_ltMatrices_solve (SEXP C, SEXP y, SEXP N, SEXP J, SEXP diag, SEXP transpose)
{

    SEXP ans;
    double *dans;
    int i, ONE = 1;

    @<RC input@>
    /* diagonal elements are always present */
    if (!Rdiag) len += iJ;
    @<C length@>
    @<lapack options@>


    char tr = 'N';
    Rboolean Rtranspose = asLogical(transpose);
    /* t(C) instead of C for column-major*/
    if (Rtranspose && !LOGICAL(abyrow)[0])
        tr = 'T';
    /* C instead for row-major*/
    if (!Rtranspose && LOGICAL(abyrow)[0])
        tr = 'T';

    PROTECT(ans = duplicate(y));
    dans = REAL(ans);
    
    /* loop over matrices, ie columns of C  / y */    
    for (i = 0; i < iN; i++) {

        /* solve linear system */
        F77_CALL(dtpsv)(&uplo, &tr, &di, &iJ, dC, dans, &ONE FCONE FCONE FCONE);
        dans += iJ;
        dC += p;
    }

    UNPROTECT(1);
    return(ans);
}
@}

and then for computing $\mC_i^{-1}$ explicitly

@d solve C
@{
SEXP R_ltMatrices_solve_C (SEXP C, SEXP N, SEXP J, SEXP diag, SEXP transpose)
{

    SEXP ans;
    int i, info;

    @<RC input@>
    /* diagonal elements are always present */
    if (!Rdiag) len += iJ;
    @<lapack options@>

    PROTECT(ans = duplicate(C));
    /* dC is defined in RC input but otherwise not used */
    dC = REAL(ans);
    
    /* loop over matrices, ie columns of C  / y */    
    for (i = 0; i < iN; i++) {

        /* compute inverse */
        F77_CALL(dtptri)(&uplo, &di, &iJ, dC, &info FCONE FCONE);
        if (info != 0)
            error("Cannot solve ltMatrices");

        dC += len;
    }

    UNPROTECT(1);
    /* note: ans always includes diagonal elements */
    return(ans);
}
@}

with \proglang{R} interface

@d solve ltMatrices
@{
solve.ltMatrices <- function(a, b, transpose = FALSE, ...) {

    diag <- attr(a, "diag")
    ### dtptri and dtpsv require diagonal elements being present
    if (!diag) diagonals(a) <- diagonals(a)
    d <- dim(a)
    J <- d[2L]
    dn <- dimnames(a)
    if (!is.double(a)) storage.mode(a) <- "double"

    if (!missing(b)) {
        if (!is.matrix(b)) b <- matrix(b, nrow = J, ncol = d[1L])
        stopifnot(nrow(b) == J)
        N <- ifelse(d[1L] == 1, ncol(b), d[1L])
        stopifnot(ncol(b) == N)
        if (!is.double(b)) storage.mode(b) <- "double"
        ret <- .Call(mvtnorm_R_ltMatrices_solve, a, b, 
                     as.integer(N), as.integer(J), as.logical(diag),
                     as.logical(transpose))
        if (d[1L] == N) {
            colnames(ret) <- dn[[1L]]
        } else {
            colnames(ret) <- colnames(b)
        }
        rownames(ret) <- dn[[2L]]
        return(ret)
    }

    if (transpose) stop("cannot compute inverse of t(a)")
    ret <- .Call(mvtnorm_R_ltMatrices_solve_C, a, 
                 as.integer(d[1L]), as.integer(J), as.logical(diag),
                 as.logical(FALSE))
    ### chol = solve(invchol); invchol = solve(chol)
    if (is.chol(ret)) {
        class(ret) <- gsub("chol", "invchol", class(ret))
    } else if (is.invchol(ret)) {
        class(ret) <- gsub("invchol", "chol", class(ret))
    }

    if (diag)
        return(ret)

    ### ret always includes diagonal elements, remove here
    ret <- ltMatrices(Lower_tri(ret, diag = FALSE), 
                      byrow = attr(a, "byrow"), names = dn[[2L]])
    return(ret)
}
@}

and some checks

<<ex-solve>>=
## 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)

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)

chk(solve(lxn, y), Mult(solve(lxn), y))
chk(solve(lxd, y), Mult(solve(lxd), y))

### recycle C
chk(solve(lxn[1,], y), as.array(solve(lxn[1,]))[,,1] %*% y)
chk(solve(lxn[rep(1, N),], y), solve(lxn[1,], y), check.attributes = FALSE)

### recycle y
chk(solve(lxn, y[,1]), solve(lxn, y[,rep(1, N)]))
@@

also for $\mC^\top_i \xvec_i = \yvec_i$

<<ex-tsolve>>=
chk(solve(lxn[1,], y, transpose = TRUE), 
    t(as.array(solve(lxn[1,]))[,,1]) %*% y)
@@

We run the same checks again switching between row- and column-major order.

<<ex-solve-byrow>>=
tn <- lxn
td <- lxd
lxn <- ltMatrices(lxn, byrow = !attr(lxn, "byrow"))
lxd <- ltMatrices(lxd, byrow = !attr(lxd, "byrow"))
<<ex-solve>>
<<ex-tsolve>>
lxn <- tn
lxd <- td
@@

\section{Log-determinants}

For computing the log-determinant $\log(\text{det}(\mC_i)) = \sum_{j = 1}^J \log(\diag(\mC_i)_j)$
we sum over the log-diagonal entries of a lower triangular matrix in
\proglang{C}, both when the data are stored in row- and column-major order:

@d logdet
@{
SEXP R_ltMatrices_logdet (SEXP C, SEXP N, SEXP J, SEXP diag, SEXP byrow) {

    SEXP ans;
    double *dans;
    int i, j, k;

    @<RC input@>
    Rboolean Rbyrow = asLogical(byrow);
    @<C length@>

    PROTECT(ans = allocVector(REALSXP, iN));
    dans = REAL(ans);

    for (i = 0; i < iN; i++) {
        dans[i] = 0.0;
        if (Rdiag) {
            k = 1;
            for (j = 0; j < iJ; j++) {
                dans[i] += log(dC[k - 1]);
                k += (Rbyrow ? j + 2 : iJ - j);
            }
            dC += p;
        }
    }
    
    UNPROTECT(1);
    return(ans);
}
@}

The \proglang{R} interface now simply calls this low-level function

@d logdet ltMatrices
@{
logdet <- function(x) {

    if (!is.ltMatrices(x))
        stop("x is not an ltMatrices object")

    byrow <- attr(x, "byrow")
    diag <- attr(x, "diag")
    d <- dim(x)
    J <- d[2L]
    dn <- dimnames(x)
    if (!is.double(x)) storage.mode(x) <- "double"

    ret <- .Call(mvtnorm_R_ltMatrices_logdet, x, 
                 as.integer(d[1L]), as.integer(J), as.logical(diag), 
                 as.logical(byrow))
    names(ret) <- dn[[1L]]
    return(ret) 
}
@}

We test the functionality by extracting the diagonal elements from different
matrices and summing over their logarithms

<<ex-logdet>>=
chk(logdet(lxn), colSums(log(diagonals(lxn))))
chk(logdet(lxd[1,]), colSums(log(diagonals(lxd[1,]))))
chk(logdet(lxd), colSums(log(diagonals(lxd))))
lxd2 <- ltMatrices(lxd, byrow = !attr(lxd, "byrow"))
chk(logdet(lxd2), colSums(log(diagonals(lxd2))))
@@

\section{Crossproducts} \label{sec:crossproducts}

We want to compute $\mC_i \mC_i^\top$ or $\text{diag}(\mC_i \mC_i^\top)$
(\code{diag\_only = TRUE}) for $i = 1, \dots, N$. These are symmetric
matrices, so we store them as a lower triangular matrix using a different
class name \code{syMatrices}. We reuse \code{Mult} 
computing $\mC_i \mC_i^\top$ or $\mC_i^\top \mC_i$ (\code{Rtranspose} being
\code{TRUE}).

We differentiate between computation of the diagonal elements of the
crossproduct and computation of the full $\J \times \J$ crossproduct matrix

@d tcrossprod ltMatrices
@{
### C %*% t(C) => returns object of class syMatrices
### diag(C %*% t(C)) => returns matrix of diagonal elements
.Tcrossprod <- function(x, diag_only = FALSE, transpose = FALSE) {

    ### crossprod(x) = t(x) %*% x
    if (transpose) {
        if (diag_only)
            return(marginSums(as.array(x^2), c(2L, 3L)))
        byrow_orig <- attr(x, "byrow")
        x <- ltMatrices(x, byrow = FALSE)
        ret <- syMatrices(Mult(x, transpose = TRUE), 
                          diag = TRUE, byrow = FALSE, 
                          names = dimnames(x)[[2L]])
        return(syMatrices(ret, byrow = byrow_orig))
    }

    ### tcrossprod(x) = x %*% t(x)
    if (diag_only)
        return(marginSums(as.array(x^2), c(1L, 3L)))
    ### t(x) for all i = 1, ..., N
    ax <- matrix(aperm(as.array(x),  c(2L, 1L, 3L)), ncol = dim(x)[1L])
    byrow_orig <- attr(x, "byrow")
    x <- ltMatrices(x, byrow = FALSE) ### not needed but Mult() is faster
    ret <- syMatrices(Mult(x, ax)[.lt(dim(x)[2L], diag = TRUE),,drop = FALSE], 
                      diag = TRUE, byrow = FALSE, 
                      names = dimnames(x)[[2L]])
    return(syMatrices(ret, byrow = byrow_orig))
}

### <FIXME> remove Tcrossprod / Crossprod from interface (mlt/tram may use
### it) </FIXME>
Tcrossprod <- function(x, diag_only = FALSE)
    .Tcrossprod(x = x, diag_only = diag_only, transpose = FALSE)
@}

We could have created yet another generic \code{tcrossprod}, but
\code{base::tcrossprod} is more general and, because speed is an issue, we
don't want to waste time on methods dispatch.

<<ex-tcrossprod>>=
## 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)

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

We also add \code{Crossprod}, which is a call to \code{Tcrossprod} with the
\code{transpose} switch turned on

@d crossprod ltMatrices
@{
Crossprod <- function(x, diag_only = FALSE)
    .Tcrossprod(x, diag_only = diag_only, transpose = TRUE)
@}

and run some checks

<<ex-crossprod>>=
## 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)

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

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)

# diagonal elements only
d <- Crossprod(lxd, diag_only = TRUE)
chk(d, apply(a, 3, diag))
chk(d, diagonals(Crossprod(lxd)))
@@

We now also add methods implementing a standard user interface

@d crossprod tcrossprod methods
@{
crossprod.ltMatrices <- function(x, y = NULL, ...) {

    if (is.null(y))
        return(Crossprod(x = x, ...))
    return(Mult(x, y, transpose = TRUE))
}

crossprod.syMatrices <- crossprod.ltMatrices

tcrossprod.ltMatrices <- function(x, y = NULL, ...) {

    if (is.null(y))
        return(Tcrossprod(x = x, ...))
    return(Mult(x, y, transpose = FALSE))
}

tcrossprod.syMatrices <- tcrossprod.ltMatrices

'%*%.ltMatrices' <- function(x, y)
    Mult(x, y)

'%*%.syMatrices' <- function(x, y)
    Mult(x, y)
@}

which allows coding closer to the standard

<<ex-tcrossprod-2>>=
### 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)
@@

Multiplication also works in this simple way
<<ex-Mult-2>>=
a <- lxn %*% y
A <- as.array(lxn)
b <- do.call("rbind", lapply(1:ncol(y), 
    function(i) t(A[,,i] %*% y[,i,drop = FALSE])))
chk(a, t(b), check.attributes = FALSE)

a <- lxd %*% y
A <- as.array(lxd)
b <- do.call("rbind", lapply(1:ncol(y), 
    function(i) t(A[,,i] %*% y[,i,drop = FALSE])))
chk(a, t(b), check.attributes = FALSE)

### recycle C
chk(lxn[rep(1, N),] %*% y, lxn[1,] %*% y, check.attributes = FALSE)

### recycle y
chk(lxn %*% y[,1], lxn %*% y[,rep(1, N)])

### tcrossprod as multiplication
i <- sample(1:N)[1]
M <- t(as.array(lxn)[,,i])
a <- sapply(1:J, function(j) lxn[i,] %*% M[,j,drop = FALSE])
rownames(a) <- colnames(a) <- dimnames(lxn)[[2L]]
b <- as.array(tcrossprod(lxn[i,]))[,,1]
chk(a, b, check.attributes = FALSE)

a <- crossprod(lxn, y)
A <- as.array(lxn)
b <- do.call("rbind", lapply(1:ncol(y), 
    function(i) t(t(A[,,i]) %*% y[,i,drop = FALSE])))
chk(a, t(b), check.attributes = FALSE)

a <- crossprod(lxd, y)
A <- as.array(lxd)
b <- do.call("rbind", lapply(1:ncol(y), 
    function(i) t(t(A[,,i]) %*% y[,i,drop = FALSE])))
chk(a, t(b), check.attributes = FALSE)

### recycle C
chk(crossprod(lxn[rep(1, N),], y), 
    crossprod(lxn[1,], y), check.attributes = FALSE)

### recycle y
chk(crossprod(lxn, y[,1]), 
    crossprod(lxn, y[,rep(1, N)]))
@@

\section{Cholesky Factorisation}

One might want to compute the Cholesky factorisations $\mSigma_i = \mC_i
\mC_i^\top$ for multiple symmetric matrices $\mSigma_i$, stored as a matrix
in class \code{syMatrices}.

@d chol syMatrices
@{
chol.syMatrices <- function(x, ...) {

    d <- dim(x)
    x <- .adddiag(x)
    if (!is.double(x)) storage.mode(x) <- "double"

    ret <- .Call(mvtnorm_R_syMatrices_chol, x, 
                 as.integer(d[1L]), as.integer(d[2L]))

    ### chol is lower triangular
    class(ret) <- gsub("syMatrices", "ltMatrices", class(ret))

    return(ret)
}
@}

Luckily, we already have the data in the correct packed colum-major storage,
so we swiftly loop over $i = 1, \dots, N$ in \proglang{C} and hand over to
\code{LAPACK}, which requires matrices in column-major order, however,
allows to switch between lower and upper triangular matrices. Thus the
row-major order for our lower triangular matrices can be computed by
pretending to be upper triangular.

@d chol
@{
SEXP R_syMatrices_chol (SEXP Sigma, SEXP N, SEXP J) {

    SEXP ans;
    double *dans;
    int iJ = INTEGER(J)[0];
    int pJ = iJ * (iJ + 1) / 2;
    int iN = INTEGER(N)[0];
    int i, info = 0;
    char uplo = 'L';

    SEXP abyrow = getAttrib(Sigma, Rf_install("byrow"));
    if (TYPEOF(abyrow) != LGLSXP) error("non-logical byrow attribute");
    if (LOGICAL(abyrow)[0])
        uplo = 'U';

    /* duplicate preserves classes, so ltMatrices is returned */
    ans = PROTECT(duplicate(Sigma));
    dans = REAL(ans);

    for (i = 0; i < iN; i++) {

        F77_CALL(dpptrf)(&uplo, &iJ, dans, &info FCONE);

        if (info != 0) {
            if (info > 0)
                error("the leading minor of order %d is not positive definite",
                      info);
            error("argument %d of Lapack routine %s had invalid value",
                  -info, "dpptrf");
        }

        dans += pJ;
    }
    UNPROTECT(1);
    return(ans);
}
@}

This new \code{chol} method can be used to revert \code{Tcrossprod} for
\code{ltMatrices} with and without unit diagonals:

<<chol>>=
Sigma <- tcrossprod(lxd)
chk(chol(Sigma), lxd)
chk(as.array(chol(Sigma)), 
    as.array(chol(syMatrices(Sigma, byrow = !attr(Sigma, "byrow")))))
Sigma <- tcrossprod(lxn)
## Sigma and chol(Sigma) always have diagonal, lxn doesn't
chk(as.array(chol(Sigma)), as.array(lxn))
### check singular matrices
x <- matrix(runif(21), ncol = 7)
S <- crossprod(x)
@@
with meaningful error messages
<<echo = FALSE>>=
cat("> ## IGNORE_RDIFF_BEGIN\n")
@@
<<cholerr>>
print(try(chol(S), silent = TRUE))
print(try(chol(as.syMatrices(S)), silent = TRUE))
@@
<<echo = FALSE>>=
cat("> ## IGNORE_RDIFF_END\n")
@@

Alternatively, we might want to compute the decomposition $\mSigma_i =
\mL_i^{-1} \mL_i^{-\top}$ for multiple symmetric matrices $\mSigma_i$, stored as a matrix
in class \code{syMatrices}. We begin, in analogy to \code{chol}, setting up
a new generic \code{invchol} with methods for \code{syMatrices} and also a
default method handling matrices.

We first recall the interpretation of the diagonal and
off-diagonal elements of $\mL$ (without the $i$ index here) from the
log-likelihood~(\ref{ll_mL_d})
\begin{eqnarray*}
\ell(\nuvec, \mL) & = & -\frac{J}{2} \log(2\pi) + \sum_{j = 1}^\J \log
  \diag(\mL)_j - \frac{1}{2} (\mL \yvec - \nuvec)^\top (\mL \yvec - \nuvec) \\
& \propto & \sum_{j = 1}^\J \log\diag(\mL)_j - 
  \frac{1}{2} (\lambdavec_j^\top \yvec - \nu_j)^2 \\
& = & \sum_{j = 1}^J \log(\lambda_{jj}) - 
  \frac{1}{2}\left(\sum_{\jmath = 1}^{j - 1} \lambda_{j\jmath} y_\jmath + \lambda_{jj} y_j - \nu_j\right)^2.
\end{eqnarray*}
where $\lambda_{jj} = \diag(\mL)_j$ and $\lambdavec_j$ is the $j$th row of
$\mL$ of length $j$.
That is, the joint log-likelihood is the sum the $\J$ log-likelihoods
corresponding to the linear regressions of $Y_j$ on $Y_1, \dots, Y_{j - 1}$,
with residual standard error $\diag(\mL)_j^{-1} = \lambda_{jj}^{-1}$ and
least-square regression coefficients $\beta_\jmath = \lambda_{j\jmath} \lambda_{jj}$
(NB: Veccia approximations use the same setup, but constrain some elements of
$\beta$ to zero).

Thus, in order to obtain $\mL$, 
we only have to compute the whole series of linear regressions from
$\mSigma$. We proceed by induction: For $j = 1$, we simply get $\lambda_{11} =
\diag(\mSigma)_1^{-\nicefrac{1}{2}}$. For $j$,
with $\mS = \sum_{i = }^N (\yvec_i - \bar{\yvec})(\yvec_i - \bar{\yvec})^\top$ we
can write $\betavec_{j} = \mS_{1:(j-1),1:(j-1)}^{-1} \mS_j$ for $j = 2,
\dots, J$ based on the partioned $j \times j $ submatrix
\begin{eqnarray*}
\mS_{\jvec,\jvec} = \begin{pmatrix}
\mS_{1:(j-1),1:(j-1)}^{-1} & \mS_j \\
\mS_j^\top & \sigma^2_j
\end{pmatrix}
\end{eqnarray*}
of $\mSigma$. This matrix corresponds to the crossproduct $\mS_{1:(j-1),1:(j-1)}^{-1}
= \mX^\top \mX$ of the
$j-1$ regressors and $\mX^\top Y = \mS_j$ when using the classical $Y =
\mX^\top \betavec + \sigma \varepsilon$ notation with $\varepsilon \sim \ND_1(0,
1)$.

We already have that $\mS_{1:(j-1),1:(j-1)} = \mL_{1:(j-1),1:(j-1)}^{-1}
\mL_{1:(j-1),1:(j-1)}^{-\top}$, so 
\begin{eqnarray*}
\betavec_j = \mL_{1:(j-1),1:(j-1)}^\top \mL_{1:(j-1),1:(j-1)} \mS_j. 
\end{eqnarray*}
We also know that the variance of 
$Y_j$ was estimated as $\sigma_j^2 = \diag(\mS)_j$ and thus
\begin{eqnarray*}
\sigma_j^2 & = & \betavec_j^\top \mS_{1:(j-1),1:(j-1),1:(j-1),1:(j-1)}
\betavec_j + \lambda_{jj}^{-2} \\
& = & \mS_j^\top \mL_{1:(j-1),1:(j-1)}^\top \mL_{1:(j-1),1:(j-1)} \mS_j + \lambda_{jj}^{-2} \\
& \iff & \lambda_{jj} = (\sigma^2_j - \|\mL_{1:(j-1),1:(j-1)} \mS_j\|^2_2)^{-\nicefrac{1}{2}}
\end{eqnarray*}
and we can update $\mL$ as
\begin{eqnarray*}
\mL_{\jvec, \jvec} = \begin{pmatrix}
\mL_{1:(j-1),1:(j-1)} & \betavec_{j} \lambda_{jj} \\
\betavec_j \lambda_{jj} & \lambda_{jj}
\end{pmatrix}.
\end{eqnarray*}
We proceed until $j = J$ and return $\mL$.

@d invchol syMatrices
@{
invchol <- function(x, ...)
    UseMethod("invchol")

invchol.default <- function(x, ...)
    invchol(as.syMatrices(as.matrix(x)))

invchol.syMatrices <- function(x, ...) {

    byrow_orig <- attr(x, "byrow")
    d <- dim(x)
    x <- .adddiag(x)

    ### x is of class syMatrices, coerse to ltMatrices first and re-arrange
    ### second
    class(x) <- gsub("syMatrices", "ltMatrices", class(x))
    x <- ltMatrices(x, byrow = TRUE)

    if (!is.double(x)) storage.mode(x) <- "double"

    ret <- .Call(mvtnorm_R_syMatrices_invchol, x, 
                 as.integer(d[1L]), as.integer(d[2L]))

    ret <- ltMatrices(ret, byrow = byrow_orig)

    return(ret)
}
@}

The heavy lifting requires some \proglang{C} code

@d invchol
@{
SEXP R_syMatrices_invchol (SEXP Sigma, SEXP N, SEXP J) {

    SEXP ans;
    double *dans;
    int iJ = INTEGER(J)[0];
    int pJ = iJ * (iJ + 1) / 2;
    int iN = INTEGER(N)[0];
    int i, info = 0;

    ans = PROTECT(duplicate(Sigma));
    dans = REAL(ans);

    for (i = 0; i < iN; i++) {

        C_invchol(iJ, dans, &info);

        if (info != 0) {
            if (info > 0)
                error("the leading minor of order %d is not positive definite",
                      info);
        }

        dans += pJ;
    }
    UNPROTECT(1);
    return(ans);
}
@}

@d invchol workhorse
@{
void C_invchol (int J, double* ans, int* info) {
	
    int i, j, k;
    int start = 1, str;
    int end = 0;
    double sd = 0.0;
    double *sigma, *x;

    info[0] = 0;

    if (ans[0] < DBL_EPSILON) {
       info[0] = 1;
    } else {
        ans[0] = 1 / sqrt(ans[0]);

        for (j = 1; j < J; j++) {
            sigma = ans + start;
            sd = 0.0;
            /* compute L mS + its 2norm */
            /* <FIXME> use dtpmv <FIXME> */
            for (i = j - 1; i >= 0; i--) {
                sigma[i] *= ans[end];
                for (k = 1; k <= i; k++)
                    sigma[i] += ans[end - k] * sigma[i - k];
                sd += pow(sigma[i], 2);
                end -= i + 1;
            }
            /* compute lambda_jj^-1 */
            sd = sigma[j] - sd;
            if (sd < DBL_EPSILON) {
                info[0] = j;
                break;
            }
            sd = sqrt(sd);
            str = 0;
            /* compute L^\top (L mS) */
            /* <FIXME> use dtpmv <FIXME> */
            for (i = 0; i < j; i++) {
                x = ans + str + i;
                sigma[i] *= x[0]; 
                for (k = i + 1; k < j; k++) {
                    x += k;
                    sigma[i] += x[0] * sigma[k];
                }
                str += i + 1;
            }
            /* jth row of L */
            for (i = 0; i < j; i++)
                sigma[i] = - sigma[i] / sd;
            sigma[j] = 1 / sd;
            start += j + 1;
            end = start - 1;
        }
    }
}
@}

<<invchol>>=
Sigma <- tcrossprod(lxd)
chk(invchol(Sigma), solve(lxd))
Sigma <- tcrossprod(lxn)
## Sigma and chol(Sigma) always have diagonal, lxn doesn't
chk(as.array(invchol(Sigma)), as.array(solve(lxn)))
### check singular matrices
x <- matrix(runif(21), ncol = 7)
S <- crossprod(x)
@@
with meaningful error messages
<<echo = FALSE>>=
# condition number is different for different systems
cat("> ## IGNORE_RDIFF_BEGIN\n")
@@
<<invcholerr>>=
print(try(solve(S), silent = TRUE))
print(try(invchol(S), silent = TRUE))
@@
<<echo = FALSE>>=
cat("> ## IGNORE_RDIFF_END\n")
@@

\section{Kronecker Products} \label{sec:vectrick}

We sometimes need to compute $\text{vec}(\mS)^\top (\mA^\top \otimes \mC)$,
where $\mS$ is a lower triangular or other $\J \times \J$ matrix and
$\mA$ and $\mC$ are lower triangular $\J \times \J$ matrices. With the ``vec
trick'', we have $\text{vec}(\mS)^\top (\mA^\top \otimes \mC) = 
\text{vec}(\mC^\top \mS \mA^\top)^\top$.  We compute $\mC^\top \mS \mA^\top$ by default or $\mC \mS \mA^\top$ or
$\mC^\top \mS \mA$ or $\mC^\top \mS \mA^\top$ by using the \code{trans}
argument in \code{vectrick}.  Argument \code{C} is an \code{ltMatrices}
object

@d check C argument
@{
C <- as.ltMatrices(C)
C_byrow_orig <- attr(C, "byrow")
dC <- dim(C)
nm <- attr(C, "rcnames")
N <- dC[1L]
J <- dC[2L]
@}

\code{S} can be an \code{ltMatrices} object, an \code{syMatrices} object, or a 
$\J^2 \times N$ matrix featuring columns of vectorised $\J \times \J$, but
not necessarily symmetric, matrices. In case \code{S} is missing, we
interpret it as \code{diag(J)} later on

@d check S argument
@{
SltM <- is.ltMatrices(S)
if (SltM || is.syMatrices(S)) {
    S_byrow_orig <- attr(S, "byrow")
    stopifnot(S_byrow_orig == C_byrow_orig)
    dS <- dim(S)
    stopifnot(dC[2L] == dS[2L])
    if (dC[1] != 1L) {
        stopifnot(dC[1L] == dS[1L])
    } else {
        N <- dS[1L]
    }
    ### we'll use matrix-matrix multiplication in Mult
    S <- matrix(as.array(S), ncol = dS[1L])
} else {
    stopifnot(is.matrix(S))
    stopifnot(nrow(S) == J^2)
    if (dC[1] != 1L) {
        stopifnot(dC[1L] == ncol(S))
    } else {
        N <- ncol(S)
    }
}
@}

\code{A} is an \code{ltMatrices} object

@d check A argument
@{
if (missing(A)) {
    A <- C
} else {
    A <- as.ltMatrices(A)
    A_byrow_orig <- attr(A, "byrow")
    stopifnot(C_byrow_orig == A_byrow_orig)
    dA <- dim(A)
    stopifnot(dC[2L] == dA[2L])
}
@}

We put everything together in function \code{vectrick} by leveraging the
matrix-matrix multiplication in \code{Mult}

@d kronecker vec trick
@{
vectrick <- function(C, S, A, transpose = c(TRUE, TRUE)) {

    stopifnot(all(is.logical(transpose)))
    stopifnot(length(transpose) == 2L)

    if (missing(S) && missing(A)) {
        if (transpose[1L] && !transpose[2L])
            return(crossprod(C)) ### t(C) %*% C
        if (!transpose[1L] && transpose[2L])
            return(tcrossprod(C)) ### C %*% t(C)
    }

    @<check C argument@>

    if (missing(S)) {  ### S = diag(J), interpreted as syMatrices obj
        SltM <- TRUE
        CS <- as.array(C)
        if (transpose[1L]) CS <- aperm(CS, perm = c(2L, 1L, 3L))
        CS <- matrix(CS, nrow = J^2)
    } else {
        @<check S argument@>
        ### CS = C %*% S or CS = t(C) %*% S
        CS <- Mult(C, S, transpose = transpose[1L])
    }

    @<check A argument@>

    ### CS %*% t(A) = t(A %*% t(CS)) or
    ### CS %*% A  = t(t(A) %*% t(CS)) 
    CS <- CS[TR <- c(t(L <- matrix(seq_len(J^2), nrow = J))),,drop = FALSE]
    ret <- Mult(A, CS, transpose = !transpose[2L])[TR,,drop = FALSE]

    if (!SltM) return(matrix(c(ret), ncol = N))

    ret <- ltMatrices(ret[L[.lt(J, diag = TRUE)],,drop = FALSE], 
                      diag = TRUE, byrow = FALSE, names = nm)
    ret <- ltMatrices(ret, byrow = attr(C, "byrow"))
    return(ret)
}
@}

Here is a small example

<<kronecker>>=
J <- 10

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

C <- solve(L)

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

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

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

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

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

SD2 <- -vectrick(a, b, a)
SD2a <- -vectrick(a, b)
chk(SD2, SD2a)

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

### same; but SD2 is vec(SD0)
S <- t(matrix(as.array(b), byrow = FALSE, nrow = 1))
SD2 <- -vectrick(a, S, a)
SD2a <- -vectrick(a, S)
chk(SD2, SD2a)

chk(c(SD0), c(SD2))

### N > 1
N <- 4L
prm <- runif(J * (J - 1) / 2)
C <- ltMatrices(prm)
S <- matrix(runif(J^2 * N), ncol = N)
A <- vectrick(C, S, C)
Cx <- as.array(C)[,,1]
B <- apply(S, 2, function(x) t(Cx) %*% matrix(x, ncol = J) %*% t(Cx))
chk(A, B)

A <- vectrick(C, S, C, transpose = c(FALSE, FALSE))
Cx <- as.array(C)[,,1]
B <- apply(S, 2, function(x) Cx %*% matrix(x, ncol = J) %*% Cx)
chk(A, B)
@@


\section{Convenience Functions} \label{sec:conv}


We add a few convenience functions for computing covariance matrices
$\mSigma_i = \mC_i \mC_i^\top$, precision matrices $\mP_i = \mL_i^\top \mL_i$, 
correlation matrices $\mR_i = \tilde{\mC}_i \tilde{\mC_i}^\top$ 
(where $\tilde{\mC}_i = \text{diag}(\mC_i \mC_i^\top)^{-\frac{1}{2}}
\mC_i)$, or matrices of partial correlations $\mA_i = -\tilde{\mL}_i^\top \tilde{\mL}_i$ with 
$\tilde{\mL}_i = \mL_i \text{diag}(\mL_i^\top \mL_i)^{-\frac{1}{2}}$
from $\mL_i$ (\code{invchol}) or $\mC_i =
\mL_i^{-1}$ (\code{chol}) for $i = 1, \dots, N$. 

Before we start, let us put a label on lower triangular matrices, such that
we can differentiate between $\mC$ and $\mL$.

@d chol classes
@{
is.chol <- function(x) inherits(x, "chol")
as.chol <- function(x) {
    stopifnot(is.ltMatrices(x))
    if (is.chol(x)) return(x)
    if (is.invchol(x))
        return(invchol2chol(x))
    class(x) <- c("chol", class(x))
    return(x)
}
is.invchol <- function(x) inherits(x, "invchol")
as.invchol <- function(x) {
    stopifnot(is.ltMatrices(x))
    if (is.invchol(x)) return(x)
    if (is.chol(x))
        return(chol2invchol(x))
    class(x) <- c("invchol", class(x))
    return(x)
}
@}

First, we set-up functions for computing $\tilde{\mC}_i$
@d D times C
@{
Dchol <- function(x, D = 1 / sqrt(Tcrossprod(x, diag_only = TRUE))) {

    if (is.invchol(x)) stop("Dchol cannot work with invchol objects")

    x <- .adddiag(x)

    N <- dim(x)[1L]
    J <- dim(x)[2L]
    nm <- dimnames(x)[[2L]]

    ### for some parameter configurations logdet(ret) would
    ### be -Inf; make sure this does't happen
    if (any(D < .Machine$double.eps))
        D[D < .Machine$double.eps] <- 2 * .Machine$double.eps
    if (any(D > 1 / .Machine$double.eps))
        D[D > 1 / .Machine$double.eps] <-  (1 / .Machine$double.eps) / 2

    if (attr(x, "byrow")) {
       idx <- rep(1:J, 1:J)
    } else {
       idx <- do.call("c", sapply(seq_len(J), 
           function(j) seq(from = j, to = J, by = 1L)))
    }
    x[] <- unclass(x) * D[idx,,drop = FALSE]

    ret <- as.chol(x)
    return(ret)
}
@}

and $\tilde{\mC}_i^{-1} = \mL_i \text{diag}(\mL_i^{-1} \mL_i^{-\top})^{\frac{1}{2}}$

@d L times D
@{
### invcholD = solve(Dchol)
invcholD <- function(x, D = sqrt(Tcrossprod(solve(x), diag_only = TRUE))) {

    if (is.chol(x)) stop("invcholD cannot work with chol objects")

    x <- .adddiag(x)

    N <- dim(x)[1L]
    J <- dim(x)[2L]
    nm <- dimnames(x)[[2L]]

    ### for some parameter configurations logdet(ret) would
    ### be -Inf; make sure this does't happen
    if (any(D < .Machine$double.eps))
        D[D < .Machine$double.eps] <- 2 * .Machine$double.eps
    if (any(D > 1 / .Machine$double.eps))
        D[D > 1 / .Machine$double.eps] <-  (1 / .Machine$double.eps) / 2

    if (attr(x, "byrow")) {
       idx <- do.call("c",sapply(seq_len(J), function(j) seq_len(j)))
    } else {
       idx <- rep(1:J, J:1)
    }
    x[] <- unclass(x) * D[idx,,drop = FALSE]

    ret <- as.invchol(x)
    return(ret)
}
@}

and now the convenience functions are one-liners:

@d convenience functions
@{
@<chol classes@>
@<D times C@>
@<L times D@>

### C -> Sigma
chol2cov <- function(x)
    Tcrossprod(x)

### Sigma -> C
cov2chol <- function(x)
    as.chol(chol(x))

### L -> C
invchol2chol <- function(x)
    as.chol(solve(x))

### C -> L
chol2invchol <- function(x)
    as.invchol(solve(x))

### L -> Sigma
invchol2cov <- function(x)
    chol2cov(invchol2chol(x))

### Sigma -> L
cov2invchol <- function(x)
    as.invchol(invchol(x))

### L -> Precision
invchol2pre <- function(x)
    Crossprod(x)

### C -> Precision
chol2pre <- function(x)
    Crossprod(chol2invchol(x))

### C -> R
chol2cor <- function(x) {
    ret <- Tcrossprod(Dchol(x))
    diagonals(ret) <- NULL
    return(ret)
}

### L -> R
invchol2cor <- function(x) {
    ret <- chol2cor(invchol2chol(x))
    diagonals(ret) <- NULL
    return(ret)
}

### L -> A
invchol2pc <- function(x) {
    ret <- -Crossprod(invcholD(x, D = 1 / sqrt(Crossprod(x, diag_only = TRUE))))
    diagonals(ret) <- 0
    ret
}

### C -> A
chol2pc <- function(x)
    invchol2pc(solve(x))
@}

Here are some tests

<<conv-ex-1>>=
prec2pc <- function(x) {
    ret <- -cov2cor(x)
    diag(ret) <- 0
    ret
}
L <- lxn
Sigma <- apply(as.array(L), 3, 
               function(x) tcrossprod(solve(x)), simplify = FALSE)
Prec <- lapply(Sigma, solve)
Corr <- lapply(Sigma, cov2cor)
CP <- lapply(Corr, solve)
PC <- lapply(Prec, function(x) prec2pc(x))
chk(unlist(Sigma), c(as.array(invchol2cov(L))), 
    check.attributes = FALSE)
chk(unlist(Prec), c(as.array(invchol2pre(L))), 
    check.attributes = FALSE)
chk(unlist(Corr), c(as.array(invchol2cor(L))), 
    check.attributes = FALSE)
chk(unlist(CP), c(as.array(crossprod(invcholD(L)))), 
    check.attributes = FALSE)
chk(unlist(PC), c(as.array(invchol2pc(L))), 
    check.attributes = FALSE)
@@

<<conv-ex-2>>=
C <- lxn
Sigma <- apply(as.array(C), 3, 
               function(x) tcrossprod(x), simplify = FALSE)
Prec <- lapply(Sigma, solve)
Corr <- lapply(Sigma, cov2cor)
CP <- lapply(Corr, solve)
PC <- lapply(Prec, function(x) prec2pc(x))
chk(unlist(Sigma), c(as.array(chol2cov(C))), 
    check.attributes = FALSE)
chk(unlist(Prec), c(as.array(chol2pre(C))), 
    check.attributes = FALSE)
chk(unlist(Corr), c(as.array(chol2cor(C))), 
    check.attributes = FALSE)
chk(unlist(CP), c(as.array(crossprod(solve(Dchol(C))))), 
    check.attributes = FALSE)
chk(unlist(PC), c(as.array(chol2pc(C))), 
    check.attributes = FALSE)
@@

<<conv-ex-3>>=
L <- lxd
Sigma <- apply(as.array(L), 3, 
               function(x) tcrossprod(solve(x)), simplify = FALSE)
Prec <- lapply(Sigma, solve)
Corr <- lapply(Sigma, cov2cor)
CP <- lapply(Corr, solve)
PC <- lapply(Prec, function(x) prec2pc(x))
chk(unlist(Sigma), c(as.array(invchol2cov(L))), 
    check.attributes = FALSE)
chk(unlist(Prec), c(as.array(invchol2pre(L))), 
    check.attributes = FALSE)
chk(unlist(Corr), c(as.array(invchol2cor(L))), 
    check.attributes = FALSE)
chk(unlist(CP), c(as.array(crossprod(invcholD(L)))), 
    check.attributes = FALSE)
chk(unlist(PC), c(as.array(invchol2pc(L))), 
    check.attributes = FALSE)
@@

<<conv-ex-4>>=
C <- lxd
Sigma <- apply(as.array(C), 3, 
               function(x) tcrossprod(x), simplify = FALSE)
Prec <- lapply(Sigma, solve)
Corr <- lapply(Sigma, cov2cor)
CP <- lapply(Corr, solve)
PC <- lapply(Prec, function(x) prec2pc(x))
chk(unlist(Sigma), c(as.array(chol2cov(C))), 
    check.attributes = FALSE)
chk(unlist(Prec), c(as.array(chol2pre(C))), 
    check.attributes = FALSE)
chk(unlist(Corr), c(as.array(chol2cor(C))), 
    check.attributes = FALSE)
chk(unlist(CP), c(as.array(crossprod(solve(Dchol(C))))), 
    check.attributes = FALSE)
chk(unlist(PC), c(as.array(chol2pc(C))), 
    check.attributes = FALSE)
@@

We also add an \code{aperm} method for class \code{ltMatrices}, 
implementing the parameters ($\tilde{\mC}_i$ or $\tilde{\mL}_i$) 
for permuted versions of the
random vectors $\rY_i$. Let $\pi$ denote a permutation of $1, \dots, J$ and
$\Pi$ the corresponding permutation matrix. Then, we have 
$\Pi \rY_i \sim \ND_\J(\mathbf{0}_\J, \Pi \mC_i \mC_i^\top \Pi^\top)$.
Unfortunately, $\Pi \mC_i$ is no longer lower triangular, so we have to find
the Cholesky decompositon $\tilde{\mC}_i \tilde{\mC}_i^\top$ of $\Pi \mC_i \mC_i^\top
\Pi^\top$. Of course, $\tilde{\mL}_i = \tilde{\mC}_i^{-1}$.

The function \code{aperm}, with argument \code{perm} $=\pi$, 
now computes the Cholesky factor $\tilde{\mC}_i$
of the permuted covariance matrix, or the inverse thereof (in case
\code{x} is of class \code{invchol}). We start with some tests

@d aperm checks
@{
J <- dim(a)[2L]
if (missing(perm)) return(a)
if (is.character(perm)) 
    perm <- match(perm, dimnames(a)[[2L]])
stopifnot(all(perm %in% 1:J))

args <- list(...)
if (length(args) > 0L)
    warning("Additional arguments", names(args), "ignored")
@}

and then implement the two methods

@d aperm
@{
aperm.chol <- function(a, perm, ...) {

    @<aperm checks@>

    return(as.chol(chol(chol2cov(a)[,perm])))
}
aperm.invchol <- function(a, perm, ...) {

    @<aperm checks@>

    return(as.invchol(invchol(invchol2cov(a)[,perm])))
}
@}

<<aperm-tests, eval= TRUE>>=
L <- as.invchol(lxn)
J <- dim(L)[2L]
Lp <- aperm(a = L, perm = p <- sample(1:J))
chk(invchol2cov(L)[,p], invchol2cov(Lp))

C <- as.chol(lxn)
J <- dim(C)[2L]
Cp <- aperm(a = C, perm = p <- sample(1:J))
chk(chol2cov(C)[,p], chol2cov(Cp))
@@

We finally add a method for class \code{ltMatrices}, for which we actually cannot
provide a reasonable result, and for symmetric matrices, where we simply
fall-back on subsetting

@d aperm
@{
aperm.ltMatrices <- function(a, perm, ...)
    stop("Cannot permute objects of class ltMatrices, 
          consider calling as.chol() or as.invchol() first")

aperm.syMatrices <- function(a, perm, ...)
    return(a[,perm])
@}

\chapter{Marginal and Conditional Normal Distributions} \label{sec:margcond}

Marginal and conditional distributions from distributions $\rY_i \sim
\ND_\J\left(\muvec_i, \mC_i \mC_i^\top\right)$
(\code{chol} argument for $\mC_i$ for $i = 1, \dots, N$) or $\rY_i \sim
\ND_\J\left(\muvec_i, \mL_i^{-1} \mL_i^{-\top}\right)$
(\code{invchol} argument for $\mL_i$ for $i = 1, \dots, N$) shall be
computed, where $\muvec_i$ is either the mean vector or is given by 
$\muvec_i = \mL_i^{-1} \nuvec_i = \mC_i \nuvec_i$.

The different parameterisations are useful in different contexts. Dropping
the subscript $i$ for notational convenience, we can write the marginal mean
vector and the covariance matrix as
\begin{eqnarray*}
\E(\rY) & = & \muvec \\
\V(\rY) & = & \mSigma = \mC \mC^\top = \mL^{-1} \mL^{-\top}.
\end{eqnarray*}
Of course, $Y_j \sim \ND_1\left(\mu_j, \sigma_j^2\right)$, where $\sigma_j^2$ is the
$j$th diagonal element of $\mSigma$. In contrast to this ``marginal
parameterisation'' where parameters have a direct interpretation in terms of
univariate marginal distributions, a ``conditional parameterisation'' in
terms of $\nuvec$ and $\mL$ is helpful when one is interested in the
conditional distributions
\begin{eqnarray*}
Y_j \mid Y_1 = y_1, \dots, Y_{j-1} = y_{j-1} \sim 
    \ND_1\left(\mL_{jj}^{-1} \nu_j - \sum_{\jmath = 1}^{j - 1} \mL_{jj}^{-1} \mL_{j\jmath} y_\jmath, 
               \mL_{jj}^{-2}\right)
\end{eqnarray*}
for each $j = 2, \dots, J$, for a proof see (\ref{fm:cdistr}). Here, the intercept is given by $\mL_{jj}^{-1} \nu_j$
and the regression coefficients by $\mL_{jj}^{-1} \mL_{j\jmath}$. Based on
these parameters, the joint density is given as
\begin{eqnarray} \label{fm:jd}
f(\yvec \mid \nuvec, \mL) & = & (2 \pi)^{-\frac{J}{2}} \det(\mL) \exp\left(-\frac{1}{2} \| \mL \yvec - \nuvec \|_2^2\right)
\\ \nonumber
& = & (2 \pi)^{-\frac{J}{2}} \prod_{j = 1}^J \mL_{jj} \exp\left(-\frac{1}{2} \| \mL \yvec - \nuvec \|_2^2\right)
\\
& = & f(y_1 \mid \nu_1, \mL_{11}) \times f(y_2 | y_1, \nu_2, \mL_{2\cdot})
\times \cdots \times f(y_J \mid y_1, \dots, y_{J - 1}, \nu_J, \mL_{J\cdot})
\nonumber
\end{eqnarray}
which nicely decomposes into a product of conditional densities.

Before we begin, we implement three convenience functions for extracting
marginal means from any parameterisation using infrastructure later to be 
developed in Chapter~\ref{inter}. We start with marginal means

@d mvnorm marginal means
@{
mean.mvnorm <- function(x, ...) {
    if (!is.null(x$mean)) return(x$mean)
    if (!is.null(x$invcholmean)) return(.i2m(x))
    return(0)
}
@}

and covariances

@d mvnorm vcov
@{
vcov.mvnorm <- function(object, ...) {
    if (is.invchol(object$scale)) return(invchol2cov(object$scale))
    return(chol2cov(object$scale))
}
@}

and also add a function for computing intercepts, regression coefficients and
residual standard errors 

@d mvnorm regression coefs
@{
coef.mvnorm <- function(object, which = dim(object$scale)[2L], ...) {

    if (dim(object$scale)[[2L]] == 1) {
        if (!is.invchol(object$scale)) object$scale <- chol2invchol(object$scale)
        if (!is.null(object$invcholmean)) object$mean <- mean(object)
        if (is.null(object$mean)) return(NULL)
        ret <- object$mean
        attr(ret, "sigma") <- 1 / c(unclass(object$scale))
        return(ret)
    }

    if (is.character(which)) 
        which <- which(dimnames(object$scale)[[2L]] == which)

    stopifnot((length(which) == 1L) && 
              (which > 1) && 
              (which <= dim(object$scale)[2L]))

    if (!is.invchol(object$scale)) object$scale <- chol2invchol(object$scale)

    diag <- attr(object$scale, "diag")
    j <- which - 1L
    start <- j * (j + 1) / 2
    idx <- start + seq_len(which)
    ret <- Lower_tri(ltMatrices(object$scale, byrow = TRUE, diag = diag), 
                     diag = diag)[idx,,drop = FALSE]
    rownames(ret) <- dimnames(object$scale)[[2L]][seq_len(which)]
    if (!is.null(nm <- dimnames(object$scale)[[1L]]))
        colnames(ret) <- nm
    dg <- ret[nrow(ret),]
    cf <- - ret[-nrow(ret),,drop = FALSE]
    if (!is.null(object$mean)) object$invcholmean <- .m2i(object)
    if (!is.null(object$invcholmean))
        cf <- rbind("(Intercept)" = object$invcholmean[which,,drop = TRUE], cf)
    cf <- cf / matrix(dg, nrow = nrow(cf), ncol = ncol(cf), byrow = TRUE)
    cf <- cf[,,drop = TRUE]
    attr(cf, "sigma") <- 1 / dg
    return(cf)
}
@}

We begin computations for marginal and conditional distributions with some simple argument checks

@d mc input checks
@{
stopifnot(xor(missing(chol), missing(invchol)))
x <- if (missing(chol)) invchol else chol

stopifnot(is.ltMatrices(x))

N <- dim(x)[1L]
J <- dim(x)[2L]

if (missing(which)) return(x)

if (is.character(which)) which <- match(which, dimnames(x)[[2L]])
stopifnot(all(which %in% 1:J))
@}

The first $j$ marginal distributions can be obtained from subsetting $\mC$
or $\mL$ directly. Arbitrary marginal distributions are based on the
corresponding subset of the covariance matrix for which we compute a
corresponding Cholesky factor (such that we can use \code{lpmvnorm} later
on).

@d marginal
@{
marg_mvnorm <- function(chol, invchol, which = 1L) {

    @<mc input checks@>

    if (which[1] == 1L && (length(which) == 1L || 
                           all(diff(which) == 1L))) {
        ### which is 1:j
        tmp <- x[,which]
    } else {
        if (missing(chol)) { 
            cv <- invchol2cov(x)
            tmp <- cov2invchol(cv[,which])
        } else {
            cv <- chol2cov(x)
            tmp <- cov2chol(cv[,which])
        }
    }

    if (missing(chol))
        ret <- list(invchol = tmp)
    else
        ret <- list(chol = tmp)

    ret
}
@}

We compute conditional distributions from the precision matrices
$\mSigma^{-1}_i = \mP_i = \mL_i^\top \mL_i$ (we omit the $i$ index from now
on). For an arbitrary subset $\jvec
\subset \{1, \dots, \J\}$, the conditional distribution of $\rY_{-\jvec}$
given $\rY_{\jvec} = \yvec_{\jvec}$ is
\begin{eqnarray*}
\rY_{-\jvec} \mid \rY_{\jvec} = \yvec_{\jvec} \sim 
  \ND_{|\jvec|}\left(\muvec_{-\jvec} - \mP^{-1}_{-\jvec,-\jvec} \mP_{-\jvec, \jvec} 
        (\yvec_{\jvec} - \muvec_{\jvec}), \mP^{-1}_{-\jvec,-\jvec}\right)
\end{eqnarray*}
and we return a Cholesky factor $\tilde{\mC}$ such that
$\mP^{-1}_{-\jvec,-\jvec} = \tilde{\mC} \tilde{\mC}^\top$ (if \code{chol} was
given) or $\tilde{\mL} = \tilde{\mC}^{-1}$ (if \code{invchol} was given).

We can implement this as
@d cond general
@{
stopifnot(!center)

if (!missing(chol)) ### chol is C = Cholesky of covariance
    P <- Crossprod(solve(chol)) ### P = t(L) %*% L with L = C^-1
else                ### invcol is L = Cholesky of precision
    P <- Crossprod(invchol)

Pw <- P[, -which]
# chol <- solve(A <- base::chol(Pw)) ### Pw = A A^\top
### arg invchol is missing, masking mvtnorm::invchol
### which we can not access bc R CMD check is not happy about it
chol <- getFromNamespace("invchol", "mvtnorm")(Pw) ### Pw = A A^\top
g0 <- matrix(0, nrow = J, ncol = NCOL(given))
g0[which,] <- given
S <- Crossprod(chol) ### P^{-1}_jj = A^-top A^-1
if (ncol(g0) %in% c(1L, N)) {
    mean <- - S %*% (P %*% g0)[-which,,drop = FALSE]
} else { ### compare to Mult() with ncol(y) !%in% (1, N)
    ### <FIXME> check dimension and if t() is needed
    mean <- t(sapply(seq_len(ncol(g0)), function(i)
        - S %*% (P %*% g0[,i])[-which,,drop = FALSE]))
}
@}

If $\jvec = \{1, \dots, j < \J \}$ and $\mL$ is given, computations simplify a lot because
with the partitioned matrix
\begin{eqnarray*}
\mL = & = & \begin{pmatrix}
 \mL_{\jvec,\jvec} & 0 \\
 \mL_{-\jvec,\jvec} & \mL_{-\jvec,-\jvec} 
 \end{pmatrix} 
\end{eqnarray*}
the conditional precision matrix is
\begin{eqnarray*}
\mP_{-\jvec, -\jvec} = (\mL^\top \mL)_{-\jvec, -\jvec} = \mL^\top_{-\jvec, -\jvec} \mL_{-\jvec, -\jvec}
\end{eqnarray*}
and thus we return $\tilde{\mL} = \mL_{-\jvec, -\jvec}$ (if \code{invchol}
was given) or $\tilde{\mC} = \mL^{-1}_{-\jvec, -\jvec}$ (if \code{chol} was
given). The conditional mean is
\begin{eqnarray*}
\muvec_{-\jvec} - \mP^{-1}_{-\jvec,-\jvec} \mP_{-\jvec, \jvec} (\yvec_{\jvec} - \muvec_{\jvec})
& = & 
\muvec_{-\jvec}  -\mL^{-1}_{-\jvec, -\jvec} \mL^{-\top}_{-\jvec, -\jvec} \mL^\top_{-\jvec, -\jvec} \mL_{-\jvec, \jvec} (\yvec_{\jvec} - \muvec_{\jvec}) \\
& = & \muvec_{-\jvec} - \mL^{-1}_{-\jvec, -\jvec} \mL_{-\jvec, \jvec} (\yvec_{\jvec} - \muvec_{\jvec}).
\end{eqnarray*}
We sometimes, for example when scores with respect to $\mL^{-1}_{-\jvec,
-\jvec}$ shall be computed in \code{slpmvnorm}, need the negative rescaled mean 
\begin{eqnarray} \label{fm:center}
\mL_{-\jvec, \jvec} (\yvec_{\jvec} - \muvec_{\jvec})
\end{eqnarray} and the \code{center = TRUE} argument triggers this values to
be returned.

An alternative derivation of this result is based on (\ref{fm:jd}), the
ratio of which and the density of the marginal distribution of $\rY_\jvec \sim
\ND_j\left(\mL_\jvec^{-1} \nuvec_\jvec, \mL_\jvec^{-1}\mL_\jvec^{-\top}\right)$ is
\begin{eqnarray*}
f(\yvec_{-\jvec} \mid \yvec_{\jvec}, \nuvec, \mL) = 
    \frac{f(\yvec \mid \nuvec, \mL)}
         {f(\yvec_{\jvec} \mid \nuvec_\jvec, \mL_{\jvec,\jvec})} 
& = & \frac{(2 \pi)^{-\frac{J}{2}} \prod_{j = 1}^J \mL_{jj} \exp\left(-\frac{1}{2} \| \mL \yvec - \nuvec \|_2^2\right)}
         {(2 \pi)^{-\frac{j}{2}} \prod_{\jmath = 1}^j \mL_{\jmath\jmath} \exp\left(-\frac{1}{2} \| \mL_\jvec \yvec_\jvec - \nuvec_\jvec \|_2^2\right)} \\
& = &       (2 \pi)^{-\frac{J - j}{2}} \prod_{\jmath = j + 1}^J \mL_{\jmath\jmath}
              \exp\left(-\frac{1}{2} ( \| \mL \yvec - \nuvec \|_2^2 - \| \mL_\jvec \yvec_\jvec - \nuvec_\jvec\|_2^2)\right) \\
& = & (2 \pi)^{-\frac{J - j}{2}} \prod_{\jmath = j + 1}^J
\mL_{\jmath\jmath}\exp\left( \| \mL_{-\jvec,-\jvec} \yvec_{-\jvec} - (\nuvec_{-\jvec} - \mL_{-\jvec, \jvec} \yvec_\jvec) \|_2^2\right)
\end{eqnarray*}
which corresponds to the distribution 
\begin{eqnarray} \label{fm:cdistr}
\rY_{-\jvec} \mid \yvec_{\jvec} \sim \ND_{J - j}\left(\mL^{-1}_{-\jvec,-\jvec}(\nuvec_{-\jvec} - \mL_{-\jvec, \jvec}
\yvec_\jvec), \mL^{-1}_{-\jvec,-\jvec} \mL^{-\top}_{-\jvec,-\jvec}\right).
\end{eqnarray}
Noting that $\nuvec = \mL \muvec$ and thus $\nuvec_{-\jvec} = \mL_{-\jvec,-\jvec}
\muvec_\jvec + \mL_{-\jvec, \jvec} \muvec_{-\jvec}$, we obtain for the
conditional mean
\begin{eqnarray*}
& & \mL^{-1}_{-\jvec,-\jvec}(\nuvec_{-\jvec} - \mL_{-\jvec, \jvec} \yvec_\jvec) \\
& = & \mL^{-1}_{-\jvec,-\jvec}(\mL_{-\jvec,-\jvec} \muvec_\jvec + \mL_{-\jvec, \jvec} \muvec_{-\jvec} - 
                \mL_{-\jvec, \jvec} \yvec_\jvec) \\
& = & \muvec_\jvec + \mL^{-1}_{-\jvec,-\jvec} \mL_{-\jvec, \jvec} \muvec_{-\jvec}
    - \mL^{-1}_{-\jvec,-\jvec} \mL_{-\jvec, \jvec} \yvec_\jvec \\
& = & \muvec_\jvec - \mL^{-1}_{-\jvec,-\jvec}\mL_{-\jvec, \jvec} (\yvec_\jvec - \muvec_\jvec)
\end{eqnarray*}



The implementation reads

@d cond simple
@{
if (which[1] == 1L && (length(which) == 1L || 
                       all(diff(which) == 1L))) {
    ### which is 1:j
    L <- if (missing(invchol)) solve(chol) else invchol
    tmp <- matrix(0, ncol = ncol(given), nrow = J - length(which))
    centerm <- Mult(L, rbind(given, tmp))		
    ### if ncol(given) is not N = dim(L)[1L] > 1, then
    ### solve() below won't work and we loop over
    ### columns of centerm
    if (dim(L)[1L] > 1 && ncol(given) != N) {
        centerm <- lapply(1:ncol(centerm), function(j)
            matrix(centerm[,j], nrow = J, ncol = N)[-which,,drop = FALSE]
        )
    } else {
        centerm <- centerm[-which,,drop = FALSE]
    }
    L <- L[,-which]
    ct <- centerm
    if (!is.matrix(ct)) ct <- do.call("rbind", ct)
    if (is.matrix(centerm)) {
        m <- -solve(L, centerm)
    } else {
        m <- do.call("rbind", lapply(centerm, function(cm) -solve(L, cm)))
    }
    if (missing(invchol)) {
        if (center)
            return(list(center = ct, chol = solve(L)))
        return(list(mean = m, chol = solve(L)))
    }
    if (center)
        return(list(center = ct, invchol = L))
    return(list(mean = m, invchol = L))
}
@}

Note that we could have avoided the general case altogether by first
computing a Cholesky decomposition of the permuted covariance matrix (such
that the conditioning variables come first). The code above only
decomposes the marginal (and thus lower-dimensional) covariance. However, we
didn't implement the \code{center = TRUE} case, so we can fall back on the
permuted version if this option is requested. Putting everything together
gives

@d conditional
@{
cond_mvnorm <- function(chol, invchol, which_given = 1L, given, center = FALSE) {

    which <- which_given
    @<mc input checks@>

    if (N == 1) N <- NCOL(given)
    stopifnot(is.matrix(given) && nrow(given) == length(which))

    @<cond simple@>

    ### general with center = TRUE => permute first and go simple
    if (center) {
        perm <- c(which, (1:J)[!(1:J) %in% which])
        if (!missing(chol))
        return(cond_mvnorm(chol = aperm(as.chol(chol), perm = perm),
                           which_given = 1:length(which), given = given,
                           center = center))
        return(cond_mvnorm(invchol = aperm(as.invchol(invchol), perm = perm),
                           which_given = 1:length(which), given = given,
                           center = center))
    }

    @<cond general@>

    if (missing(invchol)) 
        return(list(mean = mean, chol = base::chol(S)))

    return(list(mean = mean, invchol = invchol(S)))
}
@}

Let's check this against the commonly used formula based on the covariance
matrix, first for the marginal distribution

<<marg>>=
Sigma <- tcrossprod(lxd)
j <- 1:3
chk(Sigma[,j], tcrossprod(marg_mvnorm(chol = lxd, which = j)$chol))
j <- 2:4
chk(Sigma[,j], tcrossprod(marg_mvnorm(chol = lxd, which = j)$chol))

Sigma <- tcrossprod(solve(lxd))
j <- 1:3
chk(Sigma[,j], tcrossprod(solve(marg_mvnorm(invchol = lxd, which = j)$invchol)))
j <- 2:4
chk(Sigma[,j], tcrossprod(solve(marg_mvnorm(invchol = lxd, which = j)$invchol)))
@@

and then for conditional distributions. The general case is

<<cond-general>>=
Sigma <- as.array(tcrossprod(lxd))[,,1]
j <- 2:4
y <- matrix(c(-1, 2, 1), nrow = 3)

cm <- Sigma[-j, j,drop = FALSE] %*% solve(Sigma[j,j]) %*%  y
colnames(cm) <- dimnames(lxd)[[1L]][1L]
cS <- Sigma[-j, -j] - Sigma[-j,j,drop = FALSE] %*% 
      solve(Sigma[j,j]) %*% Sigma[j,-j,drop = FALSE]

cmv <- cond_mvnorm(chol = lxd[1,], which_given = j, given = y)

chk(cm, cmv$mean)
chk(cS, as.array(tcrossprod(cmv$chol))[,,1])

Sigma <- as.array(tcrossprod(solve(lxd)))[,,1]
j <- 2:4
y <- matrix(c(-1, 2, 1), nrow = 3)

cm <- Sigma[-j, j,drop = FALSE] %*% solve(Sigma[j,j]) %*%  y
colnames(cm) <- dimnames(lxd)[[1L]][1L]
cS <- Sigma[-j, -j] - Sigma[-j,j,drop = FALSE] %*% 
      solve(Sigma[j,j]) %*% Sigma[j,-j,drop = FALSE]

cmv <- cond_mvnorm(invchol = lxd[1,], which_given = j, given = y)

chk(cm, cmv$mean)
chk(cS, as.array(tcrossprod(solve(cmv$invchol)))[,,1])
@@

and the simple case is

<<cond-simple>>=
Sigma <- as.array(tcrossprod(lxd))[,,1]
j <- 1:3
y <- matrix(c(-1, 2, 1), nrow = 3)

cm <- Sigma[-j, j,drop = FALSE] %*% solve(Sigma[j,j]) %*%  y
colnames(cm) <- dimnames(lxd)[[1L]][1L]
cS <- Sigma[-j, -j] - Sigma[-j,j,drop = FALSE] %*% 
      solve(Sigma[j,j]) %*% Sigma[j,-j,drop = FALSE]

cmv <- cond_mvnorm(chol = lxd[1,], which_given = j, given = y)

chk(c(cm), c(cmv$mean))
chk(cS, as.array(tcrossprod(cmv$chol))[,,1])

Sigma <- as.array(tcrossprod(solve(lxd)))[,,1]
j <- 1:3
y <- matrix(c(-1, 2, 1), nrow = 3)

cm <- Sigma[-j, j,drop = FALSE] %*% solve(Sigma[j,j]) %*%  y
colnames(cm) <- dimnames(lxd)[[1L]][1L]
cS <- Sigma[-j, -j] - Sigma[-j,j,drop = FALSE] %*% 
      solve(Sigma[j,j]) %*% Sigma[j,-j,drop = FALSE]

cmv <- cond_mvnorm(invchol = lxd[1,], which_given = j, given = y)

chk(c(cm), c(cmv$mean))
chk(cS, as.array(tcrossprod(solve(cmv$invchol)))[,,1])
@@

\section{Continuous Log-likelihoods}

With $\rZ \sim \ND_J(0, \mI_J)$ and $\rY =  \mC_i \rZ + \muvec_i \sim
\ND_J(\muvec_i, \mC_i \mC_i^\top)$ we want to evaluate the
log-likelihood contributions for observations $\yvec_1, \dots, \yvec_N$ in a
function called \code{ldmvnorm}

@d ldmvnorm
@{
ldmvnorm <- function(obs, mean, invcholmean, chol, invchol, logLik = TRUE) {

    stopifnot(xor(missing(chol), missing(invchol)))
    if (missing(mean) && missing(invcholmean)) invcholmean <- 0
    stopifnot(xor(missing(mean), missing(invcholmean)))
    if (!is.matrix(obs)) obs <- matrix(obs, ncol = 1L)
    p <- ncol(obs)

    if (!missing(chol)) {
         @<ldmvnorm chol@>
    } else {
         @<ldmvnorm invchol@>
    }

    names(logretval) <- colnames(obs)
    if (logLik) return(sum(logretval))
    return(logretval)
}
@}

We first check if the observations $\yvec_1, \dots, \yvec_N$ are given in an
$\J \times N$ matrix \code{obs} with corresponding means $\muvec_1, \dots,
\muvec_N$ in \code{means} and return the mean-centered observations.
Alternatively, when the scaled mean $\nuvec_i = \mL_i \muvec_i$ (\code{invcholmean}) is
given, we only check the dimensions.

@d check obs
@{
.check_obs_mean <- function(obs, mean, J, N) {

    nr <- nrow(obs)
    nc <- ncol(obs)
    if (nc != N)
        stop("obs and (inv)chol have non-conforming size")
    if (nr != J)
        stop("obs and (inv)chol have non-conforming size")
    if (identical(unique(mean), 0)) return(obs)
    if (length(mean) == J) 
        return(obs - c(mean))
    if (!is.matrix(mean))
        stop("obs and mean have non-conforming size")
    if (nrow(mean) != nr)
        stop("obs and mean have non-conforming size")
    if (ncol(mean) != nc)
        stop("obs and mean have non-conforming size")
    return(obs - mean)
}
.check_obs_invcholmean <- function(obs, invcholmean, J, N) {

    nr <- nrow(obs)
    nc <- ncol(obs)
    if (nc != N)
        stop("obs and (inv)chol have non-conforming size")
    if (nr != J)
        stop("obs and (inv)chol have non-conforming size")

    if (is.null(invcholmean)) 
        return(matrix(0, nrow = J, ncol = N))
    if (identical(unique(invcholmean), 0)) 
        return(matrix(0, nrow = J, ncol = N))

    if (!is.matrix(invcholmean))
        invcholmean <- matrix(invcholmean, nrow = J)
    nr <- nrow(invcholmean)
    nc <- ncol(invcholmean)
    if (!(nc %in% c(1L, N)))
        stop("obs and (inv)chol have non-conforming size")
    if (nr != J)
        stop("invcholmean and (inv)chol have non-conforming size")

    if (ncol(invcholmean) == N) 
        return(invcholmean)
    return(matrix(invcholmean, nrow = J, ncol = N))
}
@}

With $\mSigma_i
= \mC_i \mC_i^\top$ the log-likelihood function for $\rY_i = \yvec_i$ is
\begin{eqnarray*}
\ell_i(\muvec_i, \mC_i) = -\frac{J}{2} \log(2\pi) - \frac{1}{2} \log \mid
\mSigma_i \mid - \frac{1}{2} (\yvec_i - \muvec_i)^\top \mSigma^{-1}_i (\yvec_i - \muvec_i)
\end{eqnarray*}
Because $\log \mid \mSigma_i \mid =  \log \mid \mC_i \mC_i^\top \mid = 2 \log \mid
\mC_i \mid = 2 \sum_{j = 1}^\J \log \diag(\mC_i)_j$ we get the simpler expression
\begin{eqnarray} \label{ll_mC}
\ell_i(\muvec_i, \mC_i) & = & -\frac{J}{2} \log(2\pi) - \sum_{j = 1}^\J \log \diag(\mC_i)_j - \frac{1}{2}
(\yvec_i - \muvec_i)^\top \mC_i^{-\top} \mC_i^{-1} (\yvec_i - \muvec_i).
\end{eqnarray}

With $\nuvec_i = \mC_i^{-1} \muvec_i$ (argument \code{invcholmean}) given, we get
\begin{eqnarray} \label{ll_mC_d}
\ell_i(\nuvec_i, \mC_i) & = & -\frac{J}{2} \log(2\pi) - \sum_{j = 1}^\J \log \diag(\mC_i)_j - \frac{1}{2}
(\mC_i^{-1} \yvec_i - \nuvec_i)^\top (\mC_i^{-1} \yvec_i - \nuvec_i).
\end{eqnarray}

The most attractive expression however is based on the parameterisation
using $\mL_i = \mC_i^{-1}$ and $\nuvec_i$ because we get
\begin{eqnarray} \label{ll_mL_d}
\ell_i(\nuvec_i, \mL_i) & = & -\frac{J}{2} \log(2\pi) + \sum_{j = 1}^\J \log
\diag(\mL_i)_j - \frac{1}{2} (\mL_i \yvec_i - \nuvec_i)^\top (\mL_i \yvec_i - \nuvec_i)
\end{eqnarray}
which is simultaneously log-concave in both parameters
\citep{Barrathh_Boyd_2023}.

We need to compute \code{colSums(dnorm(z, log = TRUE))} quite often. This
turns out to be time-consuming and memory intensive, so we provide a small
internal helper function focusing on the necessary computations.

@d colSumsdnorm
@{
SEXP R_ltMatrices_colSumsdnorm (SEXP z, SEXP N, SEXP J) {
    /* number of columns */
    int iN = INTEGER(N)[0];
    /* number of rows */
    int iJ = INTEGER(J)[0];
    SEXP ans;
    double *dans, Jl2pi, *dz;

    Jl2pi = iJ * log(2 * M_PI);
    PROTECT(ans = allocVector(REALSXP, iN));
    dans = REAL(ans);
    dz = REAL(z);

    for (int i = 0; i < iN; i++) {
        dans[i] = 0.0;
        for (int j = 0; j < iJ; j++)
            dans[i] += pow(dz[j], 2);
        dans[i] = - 0.5 * (Jl2pi + dans[i]);
        dz += iJ;
    }
    
    UNPROTECT(1);
    return(ans);
}
@}

@d colSumsdnorm ltMatrices
@{
.colSumsdnorm <- function(z) {
    stopifnot(is.numeric(z))
    if (!is.matrix(z))
        z <- matrix(z, nrow = 1, ncol = length(z))
    ret <- .Call(mvtnorm_R_ltMatrices_colSumsdnorm, z, ncol(z), nrow(z))
    names(ret) <- colnames(z)
    return(ret)
}
@}

The main part is now

@d ldmvnorm chol
@{
if (missing(chol))
    stop("either chol or invchol must be given")
## chol is given
if (!is.ltMatrices(chol))	### NOTE: replace with is.chol
    stop("chol is not an object of class ltMatrices")
N <- dim(chol)[1L]
N <- ifelse(N == 1, p, N)
J <- dim(chol)[2L]
if (missing(invcholmean)) {
    ## NOTE: obs is now mean-centered
    obs <- .check_obs_mean(obs = obs, mean = mean, J = J, N = N)
    z <- solve(chol, obs)
} else {
    invcholmean <- .check_obs_invcholmean(obs = obs, invcholmean = invcholmean, 
                                          J = J, N = N)
    z <- solve(chol, obs) - invcholmean
}
logretval <- .colSumsdnorm(z)
if (attr(chol, "diag"))
    logretval <- logretval - logdet(chol)
@}

where we can use the efficient implementations of \code{solve} and
\code{logdet}.

If $\mL_i = \mC_i^{-1}$ is given, we obtain
\begin{eqnarray*}
\ell_i(\muvec_i, \mL_i) & = & -\frac{J}{2} \log(2\pi) + \sum_{j = 1}^\J \log \diag(\mL_i)_j - \frac{1}{2}
(\yvec_i - \muvec_i)^\top \mL_i^\top \mL_i (\yvec_i - \muvec_i).
\end{eqnarray*}

@d ldmvnorm invchol
@{
## invchol is given
if (!is.ltMatrices(invchol)) 	### NOTE: replace with is.invchol
    stop("invchol is not an object of class ltMatrices")
N <- dim(invchol)[1L]
N <- ifelse(N == 1, p, N)
J <- dim(invchol)[2L]
if (missing(invcholmean)) {
    ## NOTE: obs is now mean-centered
    obs <- .check_obs_mean(obs = obs, mean = mean, J = J, N = N)
    z <- Mult(invchol, obs)
} else {
    invcholmean <- .check_obs_invcholmean(obs = obs, invcholmean = invcholmean, 
                                          J = J, N = N)
    z <- Mult(invchol, obs) - invcholmean
}
logretval <- .colSumsdnorm(z)
## note that the second summand gets recycled the correct number
## of times in case dim(invchol)[1L] == 1 but ncol(obs) > 1
if (attr(invchol, "diag"))
    logretval <- logretval + logdet(invchol)
@}

The score function with respect to \code{obs} is
\begin{eqnarray*}
\frac{\partial \ell_i(\muvec_i, \mL_i)}{\partial \yvec_i} = - \mL_i^\top \mL_i (\yvec_i - \muvec_i),
\end{eqnarray*}
with respect to the scaled mean $\nuvec_i$ we get
\begin{eqnarray*}
\frac{\partial \ell_i(\nuvec_i, \mL_i)}{\partial \nuvec_i} = \mL_i \yvec_i - \nuvec_i,
\end{eqnarray*}
and with respect to \code{invchol} we have
\begin{eqnarray*}
\frac{\partial \ell_i(\muvec_i, \mL_i)}{\partial \mL_i} = 
- 2 \mL_i (\yvec_i - \muvec_i) (\yvec_i - \muvec_i)^\top + \diag(\mL_i)^{-1}.
\end{eqnarray*}
The score function with respect to \code{chol} post-processes the above
score using the vec-trick~(Section~\ref{sec:vectrick}).
For the log-likelihood~(\ref{ll_mC}), the score with respect to $\mC_i$ is the sum of the score 
functions of the two terms. We start with the simpler first term
\begin{eqnarray*}
\frac{\partial - \sum_{j = 1}^\J \log \diag(\mC_i)_j}{\partial \mC_i} & = & - \diag(\mC_i)^{-1}
\end{eqnarray*}

The second term gives (we omit the mean for the sake of simplicity)
\begin{eqnarray*}
\frac{\partial  - \yvec_i^\top \mC_i^{-\top} \mC_i^{-1} \yvec_i}{\partial \mC_i}
& = & - \left. \frac{\partial \yvec_i^\top \mA^\top \mA \yvec_i}{\partial \mA} \right|_{\mA = \mC^{-1}_i}
        \left. \frac{\partial \mA^{-1}}{\partial \mA} \right|_{\mA = \mC_i} \\
& = & - 2 \vecop(\mC_i^{-1} \yvec_i \yvec_i^\top)^\top (-1) (\mC_i^{-\top} \otimes \mC_i^{-1}) \\
& = & 2 \vecop(\mC_i^{-\top} \mC_i^{-1} \yvec_i \yvec_i^\top \mC_i^{-\top})^\top
\end{eqnarray*}
In \code{sldmvnorm}, we compute the score with respect to $\mL_i$ and use
the above relationship to compute the score with respect to $\mC_i$.

@d sldmvnorm
@{
sldmvnorm <- function(obs, mean, invcholmean, chol, invchol, logLik = TRUE) {

    stopifnot(xor(missing(chol), missing(invchol)))
    if (missing(mean) && missing(invcholmean)) invcholmean <- 0
    stopifnot(xor(missing(mean), missing(invcholmean)))
    if (!is.matrix(obs)) obs <- matrix(obs, ncol = 1L)

    if (!missing(invchol)) {

        N <- dim(invchol)[1L]
        N <- ifelse(N == 1, ncol(obs), N)
        J <- dim(invchol)[2L]
        if (missing(invcholmean)) {
            obs <- .check_obs_mean(obs = obs, mean = mean, J = J, N = N)
            ## NOTE: obs is mean-centered now 
            Mix <- Mult(invchol, obs)
        } else {
            invcholmean <- .check_obs_invcholmean(obs = obs, invcholmean = invcholmean, 
                                                  J = J, N = N)
            Mix <- Mult(invchol, obs) - invcholmean
        }
        sobs <- - Mult(invchol, Mix, transpose = TRUE)

        ### score wrt invchol
        idx <- matrix(1:J, nrow = J, ncol = J)
        idx <- idx[.lt(J, diag  = attr(invchol, "diag"))]
        j <- seq_len(J - !(attr(invchol, "diag") + 0L))
        idx2 <- rep(j, rev(j))
        ret <- - obs[idx2,,drop = FALSE] * Mix[idx,,drop = FALSE]

        if (!is.null(dimnames(invchol)[[1L]]))
            colnames(ret) <- dimnames(invchol)[[1]]
        ret <- ltMatrices(ret,
                          diag = attr(invchol, "diag"), byrow = FALSE,
                          names = dimnames(invchol)[[2L]])
        ret <- ltMatrices(ret, diag = attr(invchol, "diag"), 
                          byrow = attr(invchol, "byrow"))
        if (attr(invchol, "diag")) {
            ### recycle properly
            diagonals(ret) <- diagonals(ret) + c(1 / diagonals(invchol))
        } else {
            diagonals(ret) <- 0
        }
        if (missing(invcholmean)) {
            ret <- list(obs = sobs, invchol = ret)
        } else {
            ret <- list(obs = sobs, invcholmean = Mix, invchol = ret)
        }
        if (logLik) 
            ### note: obs is obs - mean if mean was present
            ret$logLik <- ldmvnorm(obs = obs, invcholmean = invcholmean,
                                   invchol = invchol, logLik = FALSE)
        return(ret)
    }

    invchol <- solve(chol)
    ret <- sldmvnorm(obs = obs, mean = mean, invcholmean = invcholmean, 
                     invchol = invchol)
    ### this means: ret$chol <- - vectrick(invchol, ret$invchol, invchol)
    ret$chol <- as.chol(- vectrick(invchol, ret$invchol))
    ret$invchol <- NULL
    return(ret)
}
@}

\section{Application Example}

Let's say we have $\rY_i \sim \ND_\J(\mathbf{0}_J, \mC_i \mC_i^{\top})$
for $i = 1, \dots, N$ and we know the Cholesky factors $\mL_i = \mC_i^{-1}$ of the $N$
precision matrices $\Sigma^{-1}_i = \mL_i \mL_i^{\top}$. We generate $\rY_i = \mL_i^{-1}
\rZ_i$ from $\rZ_i \sim \ND_\J(\mathbf{0}_\J, \mI_\J)$.
Evaluating the corresponding log-likelihood is now straightforward and fast,
compared to repeated calls to \code{dmvnorm}

<<ex-MV>>=
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)
@@

The \code{ldmvnorm} function now also has \code{chol} and \code{invchol}
arguments such that we can use
<<ex-MV-d>>=
ll3 <- ldmvnorm(obs = Y, invchol = lt)
chk(ll1, ll3)
@@
Note that argument \code{obs} in \code{ldmvnorm} is an $\J \times N$ matrix
whereas the traditional interface in \code{dmvnorm} expects
an $N \times \J$ matrix \code{x}.
The reason is that \code{Mult} or \code{solve} work with $\J \times N$
matrices and we want to avoid matrix transposes.


Sometimes it is preferable to split the joint distribution into a marginal
distribution of some elements and the conditional distribution given these
elements. The joint density is, of course, the product of the marginal and
conditional densities and we can check if this works for our example by

<<ex-MV-mc>>=
## marginal of and conditional on these
(j <- 1:5 * 10)
md <- marg_mvnorm(invchol = lt, which = j)
cd <- cond_mvnorm(invchol = lt, which_given = j, given = Y[j,])

ll3 <- sum(dnorm(md$invchol %*% Y[j,], log = TRUE)) + 
       sum(log(diagonals(md$invchol))) +
       sum(dnorm(cd$invchol %*% (Y[-j,] - cd$mean), log = TRUE)) + 
       sum(log(diagonals(cd$invchol)))
chk(ll1, ll3)
@@


\chapter{Multivariate Normal Log-likelihoods} \label{lpmvnorm}

<<chapterseed, echo = FALSE>>=
set.seed(270312)
@@

We now discuss code for evaluating the log-likelihood
\begin{eqnarray*}
\sum_{i = 1}^N \log(p_i(\mC_i \mid \avec_i, \bvec_i))
\end{eqnarray*}

This is relatively simple to achieve using the existing \code{pmvnorm}
function, so a prototype might look like

@d lpmvnormR
@{
lpmvnormR <- function(lower, upper, mean = 0, invcholmean, center = NULL, chol, logLik = TRUE, ...) {

    ### only needed for input checks
    stopifnot(missing(invcholmean))

    ### get access to internal mvtnorm functions
    .check_obs_mean <- mvtnorm:::.check_obs_mean

    @<input checks@>

    sigma <- Tcrossprod(chol)
    S <- as.array(sigma)
    idx <- 1

    ret <- error <- numeric(N)
    for (i in 1:N) {
        if (dim(sigma)[[1L]] > 1) idx <- i
        tmp <- pmvnorm(lower = lower[,i], upper = upper[,i], sigma = S[,,idx], ...)
        ret[i] <- tmp
        error[i] <- attr(tmp, "error")
    }
    attr(ret, "error") <- error

    if (logLik)
        return(sum(log(pmax(ret, .Machine$double.eps))))

    ret
}
@}

<<fct-lpmvnormR, echo = FALSE>>=
@<lpmvnormR@>
@@

However, the underlying \proglang{FORTRAN} code first computes the Cholesky
factor based on the covariance matrix, which is clearly a waste of time.
Repeated calls to \proglang{FORTRAN} also cost some time. The code \citep[based
on and evaluated in][]{Genz_Bretz_2002} implements a
specific form of quasi-Monte-Carlo integration without allowing the user to
change the scheme (or to fall-back to simple Monte-Carlo). We therefore
implement our own simplified version, with the aim to speed-things up
such that maximum-likelihood estimation becomes a bit faster.

Let's look at an example first. This code estimates $p_1, \dots, p_{10}$ for
a $5$-dimensional normal
<<ex-lpmvnorm_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)))
@@

We want to achieve the same result a bit more general and a bit faster, by
making the code more modular and, most importantly, by providing score
functions for all arguments $\avec_i$, $\bvec_i$, and $\mC_i$.

\section{Algorithm}

@o lpmvnorm.R -cp
@{
@<R Header@>
@<lpmvnorm@>
@<slpmvnorm@>
@<ldmvnorm@>
@<sldmvnorm@>
@<ldpmvnorm@>
@<sldpmvnorm@>
@<deperma@>
@<standardize@>
@<destandardize@>
@<lpRR@>
@<slpRR@>
@}

@o lpmvnorm.c -cc
@{
@<C Header@>
#ifndef  USE_FC_LEN_T
# define USE_FC_LEN_T
#endif
#include <Rconfig.h>
#include <R_ext/BLAS.h> /* for dtrmm */
#ifndef FCONE
# define FCONE
#endif
#include <R.h>
#include <Rmath.h>
#include <Rinternals.h>
#include <Rdefines.h>
@<pnorm fast@>
@<pnorm slow@>
@<R lpmvnorm@>
@<R slpmvnorm@>
@}

We implement the algorithm described by \cite{numerical-:1992}. The key
point here is that the original $\J$-dimensional problem~(\ref{pmvnorm}) is transformed into
an integral over $[0, 1]^{\J - 1}$.

For each $i = 1, \dots, N$, do

\begin{enumerate}
  \item Input $\mC_i$ (\code{chol}), $\avec_i$ (\code{lower}), $\bvec_i$
(\code{upper}), and control parameters $\alpha$, $\epsilon$, and $M_\text{max}$ (\code{M}).

@d input checks
@{
if (!is.matrix(lower)) lower <- matrix(lower, ncol = 1)
if (!is.matrix(upper)) upper <- matrix(upper, ncol = 1)
stopifnot(isTRUE(all.equal(dim(lower), dim(upper))))

stopifnot(is.ltMatrices(chol))		### NOTE: replace with is.chol
byrow_orig <- attr(chol, "byrow")
chol <- ltMatrices(chol, byrow = TRUE)
d <- dim(chol)
### allow single matrix C
N <- ifelse(d[1L] == 1, ncol(lower), d[1L])
J <- d[2L]

stopifnot(nrow(lower) == J && ncol(lower) == N)
stopifnot(nrow(upper) == J && ncol(upper) == N)

if (!missing(mean)) {
    lower <- .check_obs_mean(lower, mean, J = J, N = N)
    upper <- .check_obs_mean(upper, mean, J = J, N = N)
}

if (!missing(invcholmean)) {
    invcholmean <- .check_obs_invcholmean(lower, invcholmean, J = J, N = N)
    center <- - invcholmean
}

if (!is.null(center)) {
    if (!is.matrix(center)) center <- matrix(center, nrow = J, ncol = N)
    stopifnot(nrow(center) == J && ncol(center) == N)
}
@}


  \item Standardise integration limits $a^{(i)}_j / c^{(i)}_{jj}$, $b^{(i)}_j / c^{(i)}_{jj}$, and rows $c^{(i)}_{j\jmath} / c^{(i)}_{jj}$ for $1 \le \jmath < j < \J$.


@d standardise
@{
if (attr(chol, "diag")) {
    ### diagonals returns J x N and lower/upper are J x N, so
    ### elementwise standardisation is simple
    dchol <- diagonals(chol)
    ### zero diagonals not allowed
    stopifnot(all(abs(dchol) > (.Machine$double.eps)))
    ac <- lower / c(dchol)
    bc <- upper / c(dchol)
    C <- Dchol(chol, D = 1 / dchol)
    if (J > 1) { ### else: univariate problem; C is no longer used
        uC <- Lower_tri(C)
    } else {
        uC <- unclass(C)
    }
} else {
    ac <- lower
    bc <- upper
    uC <- Lower_tri(chol)
}
@}


  \item Initialise $\text{intsum} = \text{varsum} = 0$, $M = 0$, $d_1 =
\Phi\left(a^{(i)}_1\right)$, $e_1 = \Phi\left(b^{(i)}_1\right)$ and $f_1 = e_1 - d_1$.

Here, \code{center} is the negative rescaled mean~(\ref{fm:center}) allowing
to evaluate probabilities for conditional multivariate normal distributions
in \code{lpmvnorm} and corresponding scores in \code{slpmvnorm}. This
feature is not present in the orginal formulation of the algorithm.

@d initialisation
@{
x0 = 0.0;
if (LENGTH(center))
    x0 = -dcenter[0];
d0 = pnorm_ptr(da[0], x0);
e0 = pnorm_ptr(db[0], x0);
emd0 = e0 - d0;
f0 = emd0;
intsum = (iJ > 1 ? 0.0 : f0);
@}

  \item Repeat

@d init logLik loop
@{
d = d0;
f = f0;
emd = emd0;
start = 0;
@}

    \begin{enumerate}

      \item Generate uniform $w_1, \dots, w_{\J - 1} \in [0, 1]$.

      \item For $j = 2, \dots, J$ set 
        \begin{eqnarray*}
            y_{j - 1} & = & \Phi^{-1}\left(d_{j - 1} + w_{j - 1} (e_{j - 1} - d_{j - 1})\right)
        \end{eqnarray*}

We either generate $w_{j - 1}$ on the fly or use pre-computed weights
(\code{w}).

@d compute y
@{
Wtmp = (W == R_NilValue ? unif_rand() : dW[j - 1]);
tmp = d + Wtmp * emd;
if (tmp < dtol) {
    y[j - 1] = q0;
} else {
    if (tmp > mdtol)
        y[j - 1] = -q0;
    else
        y[j - 1] = qnorm(tmp, 0.0, 1.0, 1L, 0L);
}
@}

        \begin{eqnarray*}
            x_{j - 1} & = & \sum_{\jmath = 1}^{j - 1} c^{(i)}_{j\jmath} y_j
\end{eqnarray*}

@d compute x
@{
x = 0.0;
if (LENGTH(center)) {
    for (k = 0; k < j; k++)
        x += dC[start + k] * (y[k] - dcenter[k]);
    x -= dcenter[j]; 
} else {
    for (k = 0; k < j; k++)
        x += dC[start + k] * y[k];
}
@}

        \begin{eqnarray*}
            d_j & = & \Phi\left(a^{(i)}_j - x_{j - 1}\right) \\
            e_j & = & \Phi\left(b^{(i)}_j - x_{j - 1}\right)
        \end{eqnarray*}

@d update d, e
@{
d = pnorm_ptr(da[j], x);
e = pnorm_ptr(db[j], x);
emd = e - d;
@}

        \begin{eqnarray*}
            f_j & = & (e_j - d_j) f_{j - 1}.
       \end{eqnarray*}

@d update f
@{
start += j;
f *= emd;
@}

We put everything together in a loop starting with the second dimension

@d inner logLik loop
@{
for (j = 1; j < iJ; j++) {

    @<compute y@>
    @<compute x@>
    @<update d, e@>
    @<update f@>
}
@}


      \item Set $\text{intsum} = \text{intsum} + f_\J$, $\text{varsum} = \text{varsum} + f^2_\J$, $M = M + 1$, 
            and $\text{error} = \sqrt{(\text{varsum}/M - (\text{intsum}/M)^2) / M}$.

@d increment
@{
intsum += f;
@}

We refrain from early stopping and error estimation. 
    
      \item[Until] $\text{error} < \epsilon$ or $M = M_\text{max}$

    \end{enumerate}
  \item Output $\hat{p}_i = \text{intsum} / M$.

We return $\log{\hat{p}_i}$ for each $i$, or we immediately sum-up over $i$.

@d output
@{
dans[0] += (intsum < dtol ? l0 : log(intsum)) - lM;
if (!RlogLik)
    dans += 1L;
@}

and move on to the next observation (note that \code{p} might be $0$ in case
$\mC_i \equiv \mC$).

@d move on
@{
da += iJ;
db += iJ;
dC += p;
if (LENGTH(center)) dcenter += iJ;
@}

\end{enumerate}

It turned out that calls to \code{pnorm} are expensive, so a slightly faster
alternative \citep[suggested by][]{Matic_Radoicic_Stefanica_2018} might
provide an alternative which can be requested from using 
(\code{fast = TRUE} in the calls to \code{lpmvnorm} and \code{slpmvnorm}):

@d pnorm fast
@{
/* see https://doi.org/10.2139/ssrn.2842681  */
const double g2 =  -0.0150234471495426236132;
const double g4 = 0.000666098511701018747289;
const double g6 = 5.07937324518981103694e-06;
const double g8 = -2.92345273673194627762e-06;
const double g10 = 1.34797733516989204361e-07;
const double m2dpi = -2.0 / M_PI; //3.141592653589793115998;

double C_pnorm_fast (double x, double m) {

    double tmp, ret;
    double x2, x4, x6, x8, x10;

    if (R_FINITE(x)) {
        x = x - m;
        x2 = x * x;
        x4 = x2 * x2;
        x6 = x4 * x2;
        x8 = x6 * x2;
        x10 = x8 * x2;
        tmp = 1 + g2 * x2 + g4 * x4 + g6 * x6  + g8 * x8 + g10 * x10;
        tmp = m2dpi * x2 * tmp;
        ret = .5 + ((x > 0) - (x < 0)) * sqrt(1 - exp(tmp)) / 2.0;
    } else {
        ret = (x > 0 ? 1.0 : 0.0);
    }
    return(ret);
}
@}

@d pnorm slow
@{
double C_pnorm_slow (double x, double m) {
    return(pnorm(x, m, 1.0, 1L, 0L));
}
@}

The \code{fast} argument can be used to switch on the faster but less
accurate version of \code{pnorm}

@d pnorm
@{
Rboolean Rfast = asLogical(fast);
double (*pnorm_ptr)(double, double) = C_pnorm_slow;
if (Rfast)
    pnorm_ptr = C_pnorm_fast;
@}

We allow a new set of weights for each observation or one set for all
observations. In the former case, the number of columns is $M \times N$ and
in the latter just $M$.

@d W length
@{
int pW = 0;
if (W != R_NilValue) {
    if (LENGTH(W) == (iJ - 1) * iM) {
        pW = 0;
    } else {
        if (LENGTH(W) != (iJ - 1) * iN * iM)
            error("Length of W incorrect");
        pW = 1;
    }
    dW = REAL(W);
}
@}

@d dimensions
@{
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;
@}

@d setup return object
@{
len = (RlogLik ? 1 : iN);
PROTECT(ans = allocVector(REALSXP, len));
dans = REAL(ans);
for (int i = 0; i < len; i++)
    dans[i] = 0.0;
@}

The case $\J = 1$ does not loop over $M$

@d univariate problem
@{
if (iJ == 1) {
    iM = 0; 
    lM = 0.0;
} else {
    lM = log((double) iM);
}
@}

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

We put the code together in a dedicated \proglang{C} function

@d R slpmvnorm variables
@{
SEXP ans;
double *da, *db, *dC, *dW, *dans, dtol = REAL(tol)[0];
double *dcenter;
double mdtol = 1.0 - dtol;
double d0, e0, emd0, f0, q0;
@}

@d R lpmvnorm
@{
SEXP R_lpmvnorm(SEXP a, SEXP b, SEXP C, SEXP center, SEXP N, SEXP J, 
                SEXP W, SEXP M, SEXP tol, SEXP logLik, SEXP fast) {

    @<R slpmvnorm variables@>
    double l0, lM, x0, intsum;
    int p, len;

    Rboolean RlogLik = asLogical(logLik);

    @<pnorm@>
    @<dimensions@>
    @<W length@>
    @<init center@>

    int start, j, k;
    double tmp, Wtmp, e, d, f, emd, x, y[(iJ > 1 ? iJ - 1 : 1)];

    @<setup return object@>

    q0 = qnorm(dtol, 0.0, 1.0, 1L, 0L);
    l0 = log(dtol);

    @<univariate problem@>

    if (W == R_NilValue)
        GetRNGstate();

    for (int i = 0; i < iN; i++) {

        x0 = 0;
        @<initialisation@>

        if (W != R_NilValue && pW == 0)
            dW = REAL(W);

        for (int m = 0; m < iM; m++) {

            @<init logLik loop@>
            @<inner logLik loop@>
            @<increment@>

            if (W != R_NilValue)
                dW += iJ - 1;
        }

        @<output@>
        @<move on@>
    }

    if (W == R_NilValue)
        PutRNGstate();

    UNPROTECT(1);
    return(ans);
}
@}

The \proglang{R} user interface consists of some checks and a call to
\proglang{C}. Note that we need to specify both \code{w} and \code{M} in
case we want a new set of weights for each observation.

@d init random seed, reset on exit
@{
### from stats:::simulate.lm
if (!exists(".Random.seed", envir = .GlobalEnv, inherits = FALSE)) 
    runif(1)
if (is.null(seed)) 
    RNGstate <- get(".Random.seed", envir = .GlobalEnv)
else {
    R.seed <- get(".Random.seed", envir = .GlobalEnv)
    set.seed(seed)
    RNGstate <- structure(seed, kind = as.list(RNGkind()))
    on.exit(assign(".Random.seed", R.seed, envir = .GlobalEnv))
}
@}

@d check and / or set integration weights
@{
if (!is.null(w) && J > 1) {
    stopifnot(is.matrix(w))
    stopifnot(nrow(w) == J - 1)
    if (is.null(M))
        M <- ncol(w)
    stopifnot(ncol(w) %in% c(M, M * N))
    if (!is.double(w)) storage.mode(w) <- "double"
} else {
    if (J > 1) {
        if (is.null(M)) stop("either w or M must be specified")
    } else {
        M <- 1L
    }
}
@}

Sometimes we want to evaluate the log-likelihood based on $\mL = \mC^{-1}$,
the inverse Cholesky factor of the covariance matrix. In this
case, we explicitly invert $\mL$ to give $\mC$ (both matrices are lower
triangular, so this is fast).

@d Cholesky of precision
@{
stopifnot(xor(missing(chol), missing(invchol)))
if (missing(chol)) chol <- solve(invchol)
@}

FIXME: \code{center} is obsolate as \code{center = x} gives the same result
as \code{invcholmean = - x}. Also in \code{slpmvnorm}.

@d lpmvnorm
@{
lpmvnorm <- function(lower, upper, mean, invcholmean, center = NULL, chol, invchol, 
                     logLik = TRUE, M = NULL, w = NULL, seed = NULL, 
                     tol = .Machine$double.eps, fast = FALSE) {

    @<init random seed, reset on exit@>
    @<Cholesky of precision@>
    @<input checks@>
    @<standardise@>
    @<check and / or set integration weights@>

    ret <- .Call(mvtnorm_R_lpmvnorm, ac, bc, uC, as.double(center), 
                 as.integer(N), as.integer(J), w, as.integer(M), as.double(tol), 
                 as.logical(logLik), as.logical(fast));
    return(ret)
}
@}


Coming back to our simple example, we get (with $25000$ simple Monte-Carlo
iterations)
<<ex-again>>=
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))
@@

Next we generate some data and compare our implementation to \code{pmvnorm}
using quasi-Monte-Carlo integration. The \code{pmvnorm}
function uses randomised Korobov rules.
The experiment here applies generalised Halton sequences. Plain Monte-Carlo
(\code{w = NULL}) will also work but produces more variable results. Results
will depend a lot on appropriate choices and it is the user's
responsibility to make sure things work as intended. If you are unsure, you
should use \code{pmvnorm} which provides a well-tested configuration.

<<ex-lpmvnorm>>= )
M <- 10000L
if (require("qrng", quietly = TRUE)) {
    ### quasi-Monte-Carlo
    W <- t(ghalton(M, d = J - 1))
} else {
    ### Monte-Carlo
    ### byrow = TRUE because adding / removing dimensions
    ### keeps the MC points for the remaining dimensions constant
    W <- matrix(runif(M * (J - 1)), nrow = J - 1, byrow = TRUE)
}

### Genz & Bretz, 2002, without early stopping (really?)
pGB <- lpmvnormR(a, b, chol = lx, logLik = FALSE, 
                algorithm = GenzBretz(maxpts = M, abseps = 0, releps = 0))
### Genz 1992 with quasi-Monte-Carlo, fast pnorm
pGqf <- exp(lpmvnorm(a, b, chol = lx, w = W, M = M, logLik = FALSE, 
                     fast = TRUE))
### Genz 1992, original Monte-Carlo, fast pnorm
pGf <- exp(lpmvnorm(a, b, chol = lx, w = NULL, M = M, logLik = FALSE, 
                    fast = TRUE))
### Genz 1992 with quasi-Monte-Carlo, R::pnorm
pGqs <- exp(lpmvnorm(a, b, chol = lx, w = W, M = M, logLik = FALSE, 
                     fast = FALSE))
### Genz 1992, original Monte-Carlo, R::pnorm
pGs <- exp(lpmvnorm(a, b, chol = lx, w = NULL, M = M, logLik = FALSE, 
                    fast = FALSE))

cbind(pGB, pGqf, pGf, pGqs, pGs)
@@

The three versions agree nicely. We now check if the code also works for
univariate problems

<<ex-uni>>=
### test univariate problem
### call pmvnorm
pGB <- lpmvnormR(a[1,,drop = FALSE], b[1,,drop = FALSE], chol = lx[,1], 
                logLik = FALSE, 
                algorithm = GenzBretz(maxpts = M, abseps = 0, releps = 0))
### call lpmvnorm
pGq <- exp(lpmvnorm(a[1,,drop = FALSE], b[1,,drop = FALSE], chol = lx[,1], 
                   logLik = FALSE))
### ground truth
ptr <- pnorm(b[1,] / c(unclass(lx[,1]))) - pnorm(a[1,] / c(unclass(lx[,1])))

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

Because the default \code{fast = FALSE} was used here, all results are
identical.

\section{Score Function}

In addition to the log-likelihood, we would also like to have access to the
scores with respect to $\mC_i$. Because every element of $\mC_i$ only enters
once, the chain rule rules, so to speak.

We need the derivatives of $d$, $e$, $y$, and $f$ with respect to the $c$
parameters
@d chol scores
@{
double dp_c[Jp], ep_c[Jp], fp_c[Jp], yp_c[(iJ > 1 ? iJ - 1 : 1) * Jp];
@}

and the derivates with respect to the mean

@d mean scores
@{
double dp_m[Jp], ep_m[Jp], fp_m[Jp], yp_m[(iJ > 1 ? iJ - 1 : 1) * Jp];
@}

and the derivates with respect to lower (\code{a})

@d lower scores
@{
double dp_l[Jp], ep_l[Jp], fp_l[Jp], yp_l[(iJ > 1 ? iJ - 1 : 1) * Jp];
@}

and the derivates with respect to upper (\code{b})

@d upper scores
@{
double dp_u[Jp], ep_u[Jp], fp_u[Jp], yp_u[(iJ > 1 ? iJ - 1 : 1) * Jp];
@}

and we start allocating the necessary memory. The output object contains the
likelihood contributions (first row), the scores with respect to the mean
(next $\J$ rows), with respect to the lower integration limits (next $\J$
rows), with respect to the upper integration limits (next $\J$ rows) and
finally with respect to the off-diagonal elements of the Cholesky factor
(last $\J (\J - 1) / 2$ rows).

@d score output object
@{
int Jp = iJ * (iJ + 1) / 2;
@<chol scores@>
@<mean scores@>
@<lower scores@>
@<upper scores@>
double dtmp, ndtmp, etmp, netmp, Wtmp, ytmp, xx;

PROTECT(ans = allocMatrix(REALSXP, Jp + 1 + 3 * iJ, iN));
dans = REAL(ans);
for (j = 0; j < LENGTH(ans); j++) dans[j] = 0.0;
@}

For each $i = 1, \dots, N$, do

\begin{enumerate}
  \item Input $\mC_i$ (\code{chol}), $\avec_i$ (\code{lower}), $\bvec_i$
(\code{upper}), and control parameters $\alpha$, $\epsilon$, and $M_\text{max}$ (\code{M}).

  \item Standardise integration limits $a^{(i)}_j / c^{(i)}_{jj}$, $b^{(i)}_j / c^{(i)}_{jj}$, and rows $c^{(i)}_{j\jmath} / c^{(i)}_{jj}$ for $1 \le \jmath < j < \J$.

Note: We later need derivatives wrt $c^{(i)}_{jj}$, so we compute derivates
wrt $a^{(i)}_j$ and $b^{(i)}_j$ and post-differentiate later.

  \item Initialise $\text{intsum} = \text{varsum} = 0$, $M = 0$, $d_1 =
\Phi\left(a^{(i)}_1\right)$, $e_1 = \Phi\left(b^{(i)}_1\right)$ and $f_1 = e_1 - d_1$.

We start initialised the score wrt to $c^{(i)}_{11}$ (the parameter is non-existent
here due to standardisation)

@d score c11
@{
if (LENGTH(center)) {
    dp_c[0] = (R_FINITE(da[0]) ? dnorm(da[0], x0, 1.0, 0L) * (da[0] - x0 - dcenter[0]) : 0);
    ep_c[0] = (R_FINITE(db[0]) ? dnorm(db[0], x0, 1.0, 0L) * (db[0] - x0 - dcenter[0]) : 0);
} else {
    dp_c[0] = (R_FINITE(da[0]) ? dnorm(da[0], x0, 1.0, 0L) * (da[0] - x0) : 0);
    ep_c[0] = (R_FINITE(db[0]) ? dnorm(db[0], x0, 1.0, 0L) * (db[0] - x0) : 0);
}
fp_c[0] = ep_c[0] - dp_c[0];
@}

@d score a, b
@{
dp_m[0] = (R_FINITE(da[0]) ? dnorm(da[0], x0, 1.0, 0L) : 0);
ep_m[0] = (R_FINITE(db[0]) ? dnorm(db[0], x0, 1.0, 0L) : 0);
dp_l[0] = dp_m[0];
ep_u[0] = ep_m[0];
dp_u[0] = 0;
ep_l[0] = 0;
fp_m[0] = ep_m[0] - dp_m[0];
fp_l[0] = -dp_m[0];
fp_u[0] = ep_m[0];
@}

  \item Repeat

@d init score loop
@{
@<init logLik loop@>
@<score c11@>
@<score a, b@>
@}

    \begin{enumerate}

      \item Generate uniform $w_1, \dots, w_{\J - 1} \in [0, 1]$.

      \item For $j = 2, \dots, J$ set 
        \begin{eqnarray*}
            y_{j - 1} & = & \Phi^{-1}\left(d_{j - 1} + w_{j - 1} (e_{j - 1} - d_{j - 1})\right)
        \end{eqnarray*}

We again either generate $w_{j - 1}$ on the fly or use pre-computed weights
(\code{w}). We first compute the scores with respect to the already existing
parameters.

@d update yp for chol
@{
ytmp = exp(- dnorm(y[j - 1], 0.0, 1.0, 1L)); // = 1 / dnorm(y[j - 1], 0.0, 1.0, 0L)

for (k = 0; k < Jp; k++) yp_c[k * (iJ - 1) + (j - 1)] = 0.0;

for (idx = 0; idx < (j + 1) * j / 2; idx++) {
    yp_c[idx * (iJ - 1) + (j - 1)] = ytmp;
    yp_c[idx * (iJ - 1) + (j - 1)] *= (dp_c[idx] + Wtmp * (ep_c[idx] - dp_c[idx]));
}
@}

@d update yp for means, lower and upper
@{
for (k = 0; k < iJ; k++)
    yp_m[k * (iJ - 1) + (j - 1)] = 0.0;

for (idx = 0; idx < j; idx++) {
    yp_m[idx * (iJ - 1) + (j - 1)] = ytmp;
    yp_m[idx * (iJ - 1) + (j - 1)] *= (dp_m[idx] + Wtmp * (ep_m[idx] - dp_m[idx]));
}
for (k = 0; k < iJ; k++)
    yp_l[k * (iJ - 1) + (j - 1)] = 0.0;

for (idx = 0; idx < j; idx++) {
    yp_l[idx * (iJ - 1) + (j - 1)] = ytmp;
    yp_l[idx * (iJ - 1) + (j - 1)] *= (dp_l[idx] + Wtmp * (dp_u[idx] - dp_l[idx]));
}
for (k = 0; k < iJ; k++)
    yp_u[k * (iJ - 1) + (j - 1)] = 0.0;

for (idx = 0; idx < j; idx++) {
    yp_u[idx * (iJ - 1) + (j - 1)] = ytmp;
    yp_u[idx * (iJ - 1) + (j - 1)] *= (ep_l[idx] + Wtmp * (ep_u[idx] - ep_l[idx]));
}
@}

        \begin{eqnarray*}
            x_{j - 1} & = & \sum_{\jmath = 1}^{j - 1} c^{(i)}_{j\jmath} y_j
\end{eqnarray*}

        \begin{eqnarray*}
            d_j & = & \Phi\left(a^{(i)}_j - x_{j - 1}\right) \\
            e_j & = & \Phi\left(b^{(i)}_j - x_{j - 1}\right)
        \end{eqnarray*}

        \begin{eqnarray*}
            f_j & = & (e_j - d_j) f_{j - 1}.
       \end{eqnarray*}

The scores with respect to $c^{(i)}_{j\jmath}, \jmath = 1, \dots, j - 1$ are

@d score wrt new chol off-diagonals
@{
dtmp = dnorm(da[j], x, 1.0, 0L);
ndtmp = dtmp * (-1.0);
etmp = dnorm(db[j], x, 1.0, 0L);
netmp = etmp * (-1.0);

for (k = 0; k < j; k++) {
    idx = start + j + k;
    if (LENGTH(center)) {    
        dp_c[idx] = ndtmp * (y[k] - dcenter[k]);
        ep_c[idx] = netmp * (y[k] - dcenter[k]);
    } else {
        dp_c[idx] = ndtmp * y[k];
        ep_c[idx] = netmp * y[k];
    }
    fp_c[idx] = (ep_c[idx] - dp_c[idx]) * f;
}
@}

and the score with respect to (the here non-existing) $c^{(i)}_{jj}$ is

@d score wrt new chol diagonal
@{
idx = (j + 1) * (j + 2) / 2 - 1;
if (LENGTH(center)) {
    dp_c[idx] = (R_FINITE(da[j]) ? dtmp * (da[j] - x - dcenter[j]) : 0);
    ep_c[idx] = (R_FINITE(db[j]) ? etmp * (db[j] - x - dcenter[j]) : 0);
} else {
    dp_c[idx] = (R_FINITE(da[j]) ? dtmp * (da[j] - x) : 0);
    ep_c[idx] = (R_FINITE(db[j]) ? etmp * (db[j] - x) : 0);
}
fp_c[idx] = (ep_c[idx] - dp_c[idx]) * f;
@}

@d new score means, lower and upper
@{
dp_m[j] = (R_FINITE(da[j]) ? dtmp : 0);
ep_m[j] = (R_FINITE(db[j]) ? etmp : 0);
dp_l[j] = dp_m[j];
ep_u[j] = ep_m[j];
dp_u[j] = 0;
ep_l[j] = 0;
fp_l[j] = - dp_m[j] * f;
fp_u[j] = ep_m[j] * f;
fp_m[j] = fp_u[j] + fp_l[j];
@}


We next update scores for parameters introduced for smaller $j$

@d update score for chol
@{
for (idx = 0; idx < j * (j + 1) / 2; idx++) {
    idxJ1 = idx * (iJ - 1);
    xx = 0.0;
    for (k = 0; k < j; k++)
        xx += dC[start + k] * yp_c[idxJ1 + k];

    dp_c[idx] = ndtmp * xx;
    ep_c[idx] = netmp * xx;
    fp_c[idx] = (ep_c[idx] - dp_c[idx]) * f + emd * fp_c[idx];
}
@}

@d update score means, lower and upper
@{
for (idx = 0; idx < j; idx++) {
    idxJ1 = idx * (iJ - 1);
    xx = 0.0;
    for (k = 0; k < j; k++)
        xx += dC[start + k] * yp_m[idxJ1 + k];

    dp_m[idx] = ndtmp * xx;
    ep_m[idx] = netmp * xx;
    fp_m[idx] = (ep_m[idx] - dp_m[idx]) * f + emd * fp_m[idx];

    xx = 0.0;
    for (k = 0; k < j; k++)
        xx += dC[start + k] * yp_l[idxJ1 + k];

    dp_l[idx] = ndtmp * xx;
    dp_u[idx] = netmp * xx;
    fp_l[idx] = (dp_u[idx] - dp_l[idx]) * f + emd * fp_l[idx];

    xx = 0.0;
    for (k = 0; k < j; k++)
        xx += dC[start + k] * yp_u[idxJ1 + k];

    ep_l[idx] = ndtmp * xx;
    ep_u[idx] = netmp * xx;
    fp_u[idx] = (ep_u[idx] - ep_l[idx]) * f + emd * fp_u[idx];
}
@}

We put everything together in a loop starting with the second dimension

@d inner score loop
@{
for (j = 1; j < iJ; j++) {

    @<compute y@>
    @<compute x@>
    @<update d, e@>
    @<update yp for chol@>
    @<update yp for means, lower and upper@>
    @<score wrt new chol off-diagonals@>
    @<score wrt new chol diagonal@>
    @<new score means, lower and upper@>
    @<update score for chol@>
    @<update score means, lower and upper@>
    @<update f@>

}
@}

      \item Set $\text{intsum} = \text{intsum} + f_\J$, $\text{varsum} = \text{varsum} + f^2_\J$, $M = M + 1$, 
            and $\text{error} = \sqrt{(\text{varsum}/M - (\text{intsum}/M)^2) / M}$.

We refrain from early stopping and error estimation. 
    
      \item[Until] $\text{error} < \epsilon$ or $M = M_\text{max}$

    \end{enumerate}
  \item Output $\hat{p}_i = \text{intsum} / M$.

We return $\log{\hat{p}_i}$ for each $i$, or we immediately sum-up over $i$.


@d score output
@{
dans[0] += f;
for (j = 0; j < Jp; j++)
    dans[j + 1] += fp_c[j];
for (j = 0; j < iJ; j++) {
    idx = Jp + j + 1;
    dans[idx] += fp_m[j];
    dans[idx + iJ] += fp_l[j];
    dans[idx + 2 * iJ] += fp_u[j];
}
@}

\end{enumerate}

@d init dans
@{
if (iM == 0) {
    dans[0] = intsum;
    dans[1] = fp_c[0];
    dans[2] = fp_m[0];
    dans[3] = fp_l[0];
    dans[4] = fp_u[0];
}
@}

We put everything together in \proglang{C}

@d R slpmvnorm
@{
SEXP R_slpmvnorm(SEXP a, SEXP b, SEXP C, SEXP center, SEXP N, SEXP J, SEXP W, 
               SEXP M, SEXP tol, SEXP fast) {

    @<R slpmvnorm variables@>
    double intsum;
    int p, idx, idxJ1;
    @<dimensions@>
    @<pnorm@>
    @<W length@>
    @<init center@>
    int start, j, k;
    double tmp, e, d, f, emd, x, x0, y[(iJ > 1 ? iJ - 1 : 1)];

    @<score output object@>

    q0 = qnorm(dtol, 0.0, 1.0, 1L, 0L);

    /* univariate problem */
    if (iJ == 1) iM = 0; 

    if (W == R_NilValue)
        GetRNGstate();

    for (int i = 0; i < iN; i++) {

        @<initialisation@>
        @<score c11@>
        @<score a, b@>
        @<init dans@>

        if (W != R_NilValue && pW == 0)
            dW = REAL(W);

        for (int m = 0; m < iM; m++) {
            @<init score loop@>
            @<inner score loop@>
            @<score output@>
            if (W != R_NilValue)
                dW += iJ - 1;
        }

        @<move on@>
        dans += Jp + 1 + 3 * iJ;
    }

    if (W == R_NilValue)
        PutRNGstate();

    UNPROTECT(1);
    return(ans);
}
@}

The \proglang{R} code is now essentially identical to \code{lpmvnorm},
however, we need to undo the effect of standardisation once the scores have
been computed

@d post differentiate mean score
@{
Jp <- J * (J + 1) / 2;
smean <- - ret[Jp + 1:J, , drop = FALSE]
if (attr(chol, "diag"))
    smean <- smean / c(dchol)
@}

@d post differentiate lower score
@{
slower <- ret[Jp + J + 1:J, , drop = FALSE]
if (attr(chol, "diag"))
    slower <- slower / c(dchol)
@}

@d post differentiate upper score
@{
supper <- ret[Jp + 2 * J + 1:J, , drop = FALSE]
if (attr(chol, "diag"))
    supper <- supper / c(dchol)
@}

@d post differentiate chol score
@{
if (J == 1) {
    idx <- 1L
} else {
    idx <- cumsum(c(1, 2:J))
}
if (attr(chol, "diag")) {
    ret <- ret / c(dchol[rep(1:J, 1:J),]) ### because 1 / dchol already there
    ret[idx,] <- -ret[idx,]
}
@}


We sometimes parameterise models in terms of $\mL = \mC^{-1}$, the Cholesky
factor of the precision matrix. The log-likelihood operates on $\mC$, so we
need to post-differentiate the score function. We have
\begin{eqnarray*}
\mA = \frac{\partial \mL^{-1}}{\partial \mL} = - \mL^{-\top} \otimes \mL^{-1}
\end{eqnarray*}
and computing $\svec \mA$ for a score vector $\svec$ with respect to $\mL$ can be
implemented by the ``vec-trick''~(Section~\ref{sec:vectrick})
\begin{eqnarray*}
\svec \mA = \mL^{-\top} \mS \mL^{-\top}
\end{eqnarray*}
where $\svec = \text{vec}(\mS)$.

@d post differentiate invchol score
@{
if (!missing(invchol)) {
    ret <- ltMatrices(ret, diag = TRUE, byrow = TRUE,
                      names = dimnames(chol)[[2L]])
    ### this means vectrick(chol, ret, chol)
    ret <- - unclass(vectrick(chol, ret))
}
@}

If the diagonal elements are constants, we set them to zero. The function
always returns an object of class \code{ltMatrices} with explicit diagonal
elements (use \code{Lower\_tri(, diag = FALSE)} to extract the lower
triangular elements such that the scores match the input)

@d post process score
@{
if (!attr(chol, "diag"))
    ### remove scores for constant diagonal elements
    ret[idx,] <- 0
ret <- ltMatrices(ret, diag = TRUE, byrow = TRUE, 
                  names = dimnames(chol)[[2L]])
@}

If the model was parameterised in terms of $\nuvec_i = \mC_i^{-1} \muvec_i$,
we compute the score for $\nuvec_i$ based on the score for $\muvec_i$

@d invcholmean score
@{
if (!missing(invcholmean))
    sinvcholmean <- Mult(chol, smean, transpose = TRUE)
@}

We can now finally put everything together in a single score function.

@d slpmvnorm
@{
slpmvnorm <- function(lower, upper, mean, invcholmean, center = NULL, 
                      chol, invchol, logLik = TRUE, M = NULL, 
                      w = NULL, seed = NULL, tol = .Machine$double.eps, 
                      fast = FALSE) {

    @<init random seed, reset on exit@>
    @<Cholesky of precision@>
    @<input checks@>
    @<standardise@>
    @<check and / or set integration weights@>

    ret <- .Call(mvtnorm_R_slpmvnorm, ac, bc, uC, as.double(center), as.integer(N), 
                 as.integer(J), w, as.integer(M), as.double(tol), as.logical(fast));

    ll <- log(pmax(ret[1L,], tol)) - log(M)
    intsum <- ret[1L,]
    m <- matrix(intsum, nrow = nrow(ret) - 1, ncol = ncol(ret), byrow = TRUE)
    ret <- ret[-1L,,drop = FALSE] / m
    ### NOTE: division by zero may have happened for observations with
    ### probability < tol; the log-lik is constant in these cases
    ### and thus the derivative = 0
    ret[m < tol] <- 0

    @<post differentiate mean score@>
    @<post differentiate lower score@>
    @<post differentiate upper score@>

    ret <- ret[1:Jp, , drop = FALSE]

    @<post differentiate chol score@>
    @<post differentiate invchol score@>
    @<post process score@>
    @<invcholmean score@>

    ret <- ltMatrices(ret, byrow = byrow_orig)

    rownames(smean) <- rownames(slower) <- 
        rownames(supper) <- dimnames(chol)[[2L]]

    if (logLik) {
        ret <- list(logLik = ll, 
                    mean = smean, 
                    lower = slower,
                    upper = supper,
                    chol = ret)
        if (!missing(invchol)) names(ret)[names(ret) == "chol"] <- "invchol"
        if (!missing(invcholmean)) ret$invcholmean <- sinvcholmean
        return(ret)
    }
    
    return(ret)
}
@}

Let's look at an example, where we use \code{numDeriv::grad} to check the
results (this functionality from package \pkg{numDeriv} was absolutely 
instrumental for this project)

<<ex-score>>=
J <- 5L
N <- 4L

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

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

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

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

chk(lli, sC$logLik)

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

We can do the same when $\mL$ (and not $\mC$) is given
<<ex-Lscore>>=
mL <- solve(mC)

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

chk(lli, lliL)

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

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

chk(lliL, sL$logLik)

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

The score function also works for univariate problems
<<ex-uni-score>>=
ptr <- pnorm(b[1,] / c(unclass(mC[,1]))) - pnorm(a[1,] / c(unclass(mC[,1])))
log(ptr)
lpmvnorm(a[1,,drop = FALSE], b[1,,drop = FALSE], chol = mC[,1], logLik = FALSE)
lapply(slpmvnorm(a[1,,drop = FALSE], b[1,,drop = FALSE], chol = mC[,1], 
                 logLik = TRUE), unclass)
sd1 <- c(unclass(mC[,1]))
(dnorm(b[1,] / sd1) * b[1,] - dnorm(a[1,] / sd1) * a[1,]) * (-1) / sd1^2 / ptr
@@

\chapter{Maximum-likelihood Example} \label{ML}

<<chapterseed, echo = FALSE>>=
set.seed(110515)
@@

We now discuss how this infrastructure can be used to estimate the Cholesky
factor of a multivariate normal in the presence of interval-censored
observations.

We first generate a covariance matrix $\Sigma = \mC \mC^\top$ and extract the Cholesky factor
$\mC$
<<ex-ML-dgp>>=
J <- 4
R <- diag(J)
R[1,2] <- R[2,1] <- .25
R[1,3] <- R[3,1] <- .5
R[2,4] <- R[4,2] <- .75
Sigma <- diag(sqrt(1:J / 2)) %*% R %*% diag(sqrt(1:J / 2))
C <- t(chol(Sigma))
@@

We now represent this matrix as \code{ltMatrices} object
<<ex-ML-C>>=
prm <- C[lower.tri(C, diag = TRUE)]
lt <- ltMatrices(matrix(prm, ncol = 1L), 
                 diag = TRUE,    ### has diagonal elements
                 byrow = FALSE)  ### prm is column-major
BYROW <- FALSE   ### later checks
lt <- ltMatrices(lt, 
                 byrow = BYROW)   ### convert to row-major
chk(C, as.array(lt)[,,1], check.attributes = FALSE)
chk(Sigma, as.array(tcrossprod(lt))[,,1], check.attributes = FALSE)
@@

We generate some data from $\ND_\J(\mathbf{0}_\J, \Sigma)$ by first sampling
from $\rZ \sim \ND_\J(\mathbf{0}_\J, \mI_\J)$ and then computing $\rY = \mC \rZ +
\muvec \sim \ND_\J(\muvec, \mC \mC^\top)$

<<ex-ML-data>>=
N <- 100L
Z <- matrix(rnorm(N * J), nrow = J)
Y <- lt %*% Z + (mn <- 1:J)
@@

Before we add some interval-censoring to the data, let's estimate the
Cholesky factor $\mC$ (here called \code{lt}) from the raw continuous data.
The true mean $\muvec$ and the true covariance matrix $\Sigma$ can be estimated 
from the uncensored data via maximum likelihood as

<<ex-ML-mu-vcov>>=
rowMeans(Y)
(Shat <- var(t(Y)) * (N - 1) / N)
@@

We first check if we can obtain the same results by numerial optimisation
using \code{dmvnorm} and the scores \code{sldmvnorm}. The log-likelihood and
the score function (for the centered means) in terms of $\mC$ are

<<ex-ML-clogLik>>=
Yc <- Y - rowMeans(Y)

ll <- function(parm) {
    C <- ltMatrices(parm, diag = TRUE, byrow = BYROW)
    -ldmvnorm(obs = Yc, chol = C)
}

sc <- function(parm) {
    C <- ltMatrices(parm, diag = TRUE, byrow = BYROW)
    -rowSums(unclass(sldmvnorm(obs = Yc, chol = C)$chol))
}
@@

The diagonal elements of $\mC$ are positive, so we need box constraints
<<ex-ML-const>>=
llim <- rep(-Inf, J * (J + 1) / 2)
llim[which(rownames(unclass(lt)) %in% paste(1:J, 1:J, sep = "."))] <- 1e-4
@@

The ML-estimate of $\mC \mC^\top$ is now used to obtain an estimate of $\mC$
and we check the score function for some random starting values
<<ex-ML-c>>=
if (BYROW) {
  cML <- chol(Shat)[upper.tri(Shat, diag = TRUE)]
} else {
  cML <- t(chol(Shat))[lower.tri(Shat, diag = TRUE)]
}
ll(cML)
start <- runif(length(cML))
if (require("numDeriv", quietly = TRUE))
    chk(grad(ll, start), sc(start), check.attributes = FALSE)
@@

Finally, we hand over to \code{optim} and compare the results of the
analytically and numerically obtained ML estimates

<<ex-ML-coptim>>=
op <- optim(start, fn = ll, gr = sc, method = "L-BFGS-B", 
            lower = llim, control = list(trace = FALSE))
## ML numerically
ltMatrices(op$par, diag = TRUE, byrow = BYROW)
ll(op$par)
## ML analytically
t(chol(Shat))
ll(cML)
## true C matrix
lt
@@

Under interval-censoring, the mean and $\mC$ are no longer orthogonal and
there is no analytic solution to the ML estimation problem. So, 
we add some interval-censoring represented by \code{lwr} and \code{upr} and
try to estimate the model parameters via \code{lpmvnorm} and corresponding
scores \code{slpmvnorm}.

<<ex-ML-cens>>=
prb <- 1:9 / 10
sds <- sqrt(diag(Sigma))
ct <- sapply(1:J, function(j) qnorm(prb, mean = mn[j], sd = sds[j])) 
lwr <- upr <- Y
for (j in 1:J) {
    f <- cut(Y[j,], breaks = c(-Inf, ct[,j], Inf))
    lwr[j,] <- c(-Inf, ct[,j])[f]
    upr[j,] <- c(ct[,j], Inf)[f]
}
@@


Let's do some sanity and performance checks first. For different values of
$M$, we evaluate the log-likelihood using \code{pmvnorm} (called in
\code{lpmvnormR}) and the simplified implementation (fast and slow). The
comparison is a bit
unfair, because we do not add the time needed to setup Halton sequences, but
we would do this only once and use the stored values for repeated
evaluations of a log-likelihood (because the optimiser expects a
deterministic function to be optimised)

<<ex-ML-chk, eval = FALSE>>=
M <- floor(exp(0:25/10) * 1000)
lGB <- sapply(M, function(m) {
    st <- system.time(ret <- 
        lpmvnormR(lwr, upr, mean = mn, chol = lt, algorithm = 
                  GenzBretz(maxpts = m, abseps = 0, releps = 0)))
    return(c(st["user.self"], ll = ret))
})
lH <- sapply(M, function(m) {
    W <- NULL
    if (require("qrng", quietly = TRUE))
        W <- t(ghalton(m, d = J - 1))
    st <- system.time(ret <- lpmvnorm(lwr, upr, mean = mn, 
                                      chol = lt, w = W, M = m))
    return(c(st["user.self"], ll = ret))
})
lHf <- sapply(M, function(m) {
    W <- NULL
    if (require("qrng", quietly = TRUE))
        W <- t(ghalton(m, d = J - 1))
    st <- system.time(ret <- lpmvnorm(lwr, upr, mean = mn, chol = lt, 
                                      w = W, M = m, fast = TRUE))
    return(c(st["user.self"], ll = ret))
})
@@
The evaluated log-likelihoods and corresponding timings are given in
Figure~\ref{lleval}. It seems that for $M \ge 3000$, results are reasonably
stable.

\begin{figure}
\begin{center}
<<ex-ML-fig-data, echo = FALSE>>=
### use pre-computed data, otherwise CRAN complains.
M <-
c(1000, 1105, 1221, 1349, 1491, 1648, 1822, 2013, 2225, 2459, 
2718, 3004, 3320, 3669, 4055, 4481, 4953, 5473, 6049, 6685, 7389, 
8166, 9025, 9974, 11023, 12182)
lGB <- matrix(c(0.054, -880.492612, 0.054, -880.492426, 0.054, -880.492996, 0.054,
-880.492629, 0.054, -880.490231, 0.055, -880.492784, 0.054, -880.492632,
0.055, -880.489297, 0.054, -880.492516, 0.054, -880.491339, 0.054,
-880.492091, 0.11, -880.491601, 0.114, -880.493553, 0.111, -880.49125,
0.108, -880.492151, 0.108, -880.492275, 0.109, -880.491879, 0.109,
-880.492008, 0.192, -880.492132, 0.195, -880.491839, 0.194, -880.492139,
0.194, -880.491042, 0.198, -880.492198, 0.328, -880.4916, 0.323,
-880.491941, 0.323, -880.491698), nrow = 2)
rownames(lGB) <- c("user.self", "ll")
lH <- matrix(c(0.023, -880.480296, 0.027, -880.496166, 0.029, -880.488683,
0.032, -880.496171, 0.035, -880.485597, 0.039, -880.491333, 0.043,
-880.494557, 0.048, -880.495429, 0.053, -880.494391, 0.06, -880.485546,
0.064, -880.491455, 0.071, -880.494138, 0.079, -880.491619, 0.087,
-880.493393, 0.095, -880.492541, 0.106, -880.491649, 0.118, -880.492508,
0.129, -880.492558, 0.141, -880.492509, 0.157, -880.490448, 0.173,
-880.491686, 0.193, -880.491178, 0.211, -880.492286, 0.233, -880.491511,
0.258, -880.49153, 0.287, -880.491929), nrow = 2) 
rownames(lH) <- c("user.self", "ll")
lHf <- matrix(c(0.018, -880.487067, 0.019, -880.488639, 0.022, -880.488569,
0.024, -880.49393, 0.026, -880.486029, 0.029, -880.491563, 0.033,
-880.499415, 0.035, -880.494457, 0.038, -880.493954, 0.043, -880.493648,
0.047, -880.492955, 0.052, -880.494667, 0.059, -880.493745, 0.065,
-880.494195, 0.07, -880.49333, 0.078, -880.491451, 0.086, -880.492379,
0.094, -880.490392, 0.106, -880.491061, 0.115, -880.491577, 0.129,
-880.492523, 0.142, -880.491027, 0.158, -880.492086, 0.171, -880.492069,
0.189, -880.492251, 0.208, -880.492347), nrow = 2) 
rownames(lHf) <- c("user.self", "ll")
@@
<<ex-ML-fig, eval = TRUE, echo = FALSE, fig = TRUE, pdf = TRUE, width = 8, height = 5>>=
layout(matrix(1:2, nrow = 1))
plot(M, lGB["ll",], ylim = range(c(lGB["ll",], lH["ll",], lHf["ll",])), ylab = "Log-likelihood")
points(M, lH["ll",], pch = 4)
points(M, lHf["ll",], pch = 5)
plot(M, lGB["user.self",], ylim = c(0, max(lGB["user.self",])), ylab = "Time (in sec)")
points(M, lH["user.self",], pch = 4)
points(M, lHf["user.self",], pch = 5)
legend("bottomright", legend = c("pmvnorm", "lpmvnorm", "lpmvnorm(fast)"), pch = c(1, 4, 5), bty = "n")
@@
\caption{Evaluated log-likelihoods (left) and timings (right).
\label{lleval}}
\end{center}
\end{figure}

We now define the log-likelihood function. It is important to use weights
via the \code{w} argument (or to set the \code{seed}) such that only the
candidate parameters \code{parm} change with repeated calls to \code{ll}. We
use an extremely low number of integration points \code{M}, let's see if
this still works out.

<<ex-ML-ll, eval = TRUE>>=
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, J) {
     m <- parm[1:J]             ### mean parameters
     parm <- parm[-(1:J)]       ### chol parameters
     C <- matrix(c(parm), ncol = 1L)
     C <- ltMatrices(C, diag = TRUE, byrow = BYROW)
     -lpmvnorm(lower = lwr, upper = upr, mean = m, chol = C, 
               w = W, M = M, logLik = TRUE)
}
@@

We can check the correctness of our log-likelihood function
<<ex-ML-check>>=
prm <- c(mn, unclass(lt))
ll(prm, J = J)
### ATLAS gives -880.4908, M1mac gives -880.4911
round(lpmvnormR(lwr, upr, mean = mn, chol = lt, 
                algorithm = GenzBretz(maxpts = M, abseps = 0, releps = 0)), 3)
(llprm <- lpmvnorm(lwr, upr, mean = mn, chol = lt, w = W, M = M))
chk(llprm, sum(lpmvnorm(lwr, upr, mean = mn, chol = lt, w = W, 
                        M = M, logLik = FALSE)))
@@

Before we hand over to the optimiser, we define the score function with
respect to $\muvec$ and $\mC$

<<ex-ML-sc>>=
sc <- function(parm, J) {
    m <- parm[1:J]             ### mean parameters
    parm <- parm[-(1:J)]       ### chol parameters
    C <- matrix(c(parm), ncol = 1L)
    C <- ltMatrices(C, diag = TRUE, byrow = BYROW)
    ret <- slpmvnorm(lower = lwr, upper = upr, mean = m, chol = C, 
                     w = W, M = M, logLik = TRUE)
    return(-c(rowSums(ret$mean), rowSums(unclass(ret$chol))))
}
@@

and check the correctness numerically

<<ex-ML-sc-chk>>=
if (require("numDeriv", quietly = TRUE))
    chk(grad(ll, prm, J = J), sc(prm, J = J), check.attributes = FALSE)
@@


Finally, we can hand-over to \code{optim}. Because we need $\text{diag}(\mC) >
0$, we use box constraints and \code{method = "L-BFGS-B"}. We start with the
estimates obtained from the original continuous data.

<<ex-ML>>=
llim <- rep(-Inf, J + J * (J + 1) / 2)
llim[J + which(rownames(unclass(lt)) %in% paste(1:J, 1:J, sep = "."))] <- 1e-4

if (BYROW) {
  start <- c(rowMeans(Y), chol(Shat)[upper.tri(Shat, diag = TRUE)])
} else {
  start <- c(rowMeans(Y), t(chol(Shat))[lower.tri(Shat, diag = TRUE)])
}

ll(start, J = J)

op <- optim(start, fn = ll, gr = sc, J = J, method = "L-BFGS-B", 
            lower = llim, control = list(trace = FALSE))

op$value ## compare with 
ll(prm, J = J)
@@

We can now compare the true and estimated Cholesky factor $\mC$ of our covariance
matrix $\mSigma = \mC \mC^\top$
<<ex-ML-C>>=
(C <- ltMatrices(matrix(op$par[-(1:J)], ncol = 1), 
                 diag = TRUE, byrow = BYROW))
lt
@@
and the estimated means
<<ex-ML-mu>>=
op$par[1:J]
mn
@@

We can also compare the results on the scale of the covariance matrix

<<ex-ML-Shat>>=
### ATLAS print issues
round(tcrossprod(lt), 4)  ### true Sigma
tcrossprod(C)             ### interval-censored obs
Shat                      ### "exact" obs
@@

This looks reasonably close.

\textbf{Warning:} Do NOT assume the choices made here (especially \code{M}
and \code{W}) to be universally applicable. Make sure to investigate the
accuracy depending on these parameters 
of the log-likelihood and score function in your application.

One could ask what this whole exercise was about statistically. We
estimated a multivariate normal distribution from interval-censored data, so
what? Maybe we were primarily interested in fitting a linear regression 
\begin{eqnarray*}
\E(Y_1 \mid Y_j = y_j, j = 2, \dots, J) = \alpha + \sum_{j = 2}^J \beta_j y_j.
\end{eqnarray*}
Interval-censoring in the response could have been handled by some Tobit model, but
what about interval-censoring in the explanatory variables? Based on the
multivariate distribution just estimated, we can obtain the regression
coefficients $\beta_j$ as
<<regressions>>=
c(cond_mvnorm(chol = C, which_given = 2:J, given = diag(J - 1))$mean)
@@
Alternatively, we can compute these regressions from a permuted Cholesky
factor (this goes into the ``simple'' conditional distribution in Section~\ref{sec:margcond})
<<regressionsC>>=
c(cond_mvnorm(chol = aperm(as.chol(C), perm = c(2:J, 1)),
              which_given = 1:(J - 1), given = diag(J - 1))$mean)
@@
or, as a third option, from the last row of the precision matrix of the 
permuted Cholesky factor
<<regressionsP>>=
x <- as.array(chol2pre(aperm(as.chol(C), perm = c(2:J, 1))))[J,,1]
c(-x[-J] / x[J])
@@
In higher dimensions, the first option is to be preferred, because it
only involves computing the Cholesky decomposition of a $(\J - 1) \times (\J -
1)$ matrix, whereas the latter two options are based on a decomposition of
the full $\J \times \J$ covariance matrix.

We can compare these estimated regression coefficients with those obtained
from a linear model fitted to the exact observations
<<lm-ex>>=
dY <- as.data.frame(t(Y))
colnames(dY) <- paste0("Y", 1:J)
coef(m1 <- lm(Y1 ~ ., data = dY))[-1L]
@@
The estimates are quite close, but what about standard errors?
Interval-censoring means loss of information, so we should see larger
standard errors for the interval-censored data.

Let's obtain the Hessian for all parameters first
<<hessian>>=
H <- optim(op$par, fn = ll, gr = sc, J = J, method = "L-BFGS-B", 
           lower = llim, hessian = TRUE, 
           control = list(trace = FALSE))$hessian
@@
and next we sample from the distribution of the maximum-likelihood
estimators
<<ML-sample>>=
L <- try(t(chol(H)))
### some check on r-oldrel-macos-arm64
if (inherits(L, "try-error"))
    L <- t(chol(H + 1e-4 * diag(nrow(H))))
L <- ltMatrices(L[lower.tri(L, diag = TRUE)], diag = TRUE)
Nsim <- 50000
Z <- matrix(rnorm(Nsim * nrow(H)), ncol = Nsim)
rC <- solve(L, Z)[-(1:J),] + op$par[-(1:J)] ### remove mean parameters
@@
The standard error in this sample should be close to the ones obtained from
the inverse Fisher information
<<ML-check>>=
c(sqrt(rowMeans((rC - rowMeans(rC))^2)))
c(sqrt(diagonals(crossprod(solve(L)))))
@@
We now coerse the matrix \code{rC} to an object of class \code{ltMatrices}
<<rC>>=
rC <- ltMatrices(rC, diag = TRUE)
@@
The object \code{rC} contains all sampled Cholesky factors of the covariance
matrix. From each of these matrices, we compute the regression coefficient,
giving us a sample we can use to compute standard errors from
<<ML-beta>>=
rbeta <- cond_mvnorm(chol = rC, which_given = 2:J, given = diag(J - 1))$mean
sqrt(rowMeans((rbeta - rowMeans(rbeta))^2))
@@
which are, as expected, slightly different from the ones obtained from the more
informative exact observations
<<se-ex>>=
sqrt(diag(vcov(m1)))[-1L]
@@

\chapter{Continuous-discrete Likelihoods} \label{cdl}

We sometimes are faced with outcomes measured at different levels of
precision. Some variables might have been observed very exactly, and
therefore we might want to use the log-Lebesque density for defining the
log-likelihood. Other variables might be available as relatively wide intervals
only, and thus the log-likelihood is a log-probability. We can use the
infrastructure developed so far to compute a joint likelihood. Let's assume 
we are interested in the joint distribution of $(\rY_i, \rX_i)$ and we
observed $\rY_i = \yvec_i$ (that is, exact observations of $\rY$) and 
$\avec_i < \rX_i \le \bvec_i$ (that is, interval-censored observations for
$\rX_i$). We define the log-likelihood based on the joint normal distribution $(\rY_i,
\rX_i) \sim \ND_J((\muvec_i, \etavec_i)^\top, \mC_i \mC_i^\top)$ as
\begin{eqnarray*}
\ell_i(\muvec_i, \etavec_i, \mC_i) = \ell_i(\muvec_i, \mC_{\rY,i}) + 
    \log(\Prob(\avec_i < \rX_i \le \bvec_i \mid \mC_i, \muvec_i, \etavec_i, \rY_i = \yvec_i)).
\end{eqnarray*}
where $\mC_{\rY,i}$ is the upper part of $\mC_i$ corresponding to the
marginal distribution of $\rY_i$. The conditional probability of $\rX$ given
$\rY$ depends on all parameters, as explained in Section~\ref{sec:margcond}.
The trick here is to decompose the joint likelihood into a product of the 
marginal Lebesque density of $\rY_i$ and the conditional probability of
$\rX_i$ given $\rY_i = \yvec_i$.

We first check the data

@d dp input checks
@{
stopifnot(xor(missing(chol), missing(invchol)))
cJ <- nrow(obs)
dJ <- nrow(lower)
N <- ncol(obs)
stopifnot(N == ncol(lower))
stopifnot(N == ncol(upper))
if (missing(mean)) {
    cmean <- 0
    dmean <- 0
} else {
    if (!is.matrix(mean) || NCOL(mean) == 1L) 
        mean <- matrix(mean, nrow = cJ + dJ, ncol = N)
    stopifnot(nrow(mean) == cJ + dJ)
    stopifnot(ncol(mean) == N)
    cmean <- mean[1:cJ,, drop = FALSE]
    dmean <- mean[-(1:cJ),, drop = FALSE]
}
@}

We can use \code{marg\_mvnorm} and \code{cond\_mvnorm} to compute the
marginal and the conditional normal distributions and the joint log-likelihood
is simply the sum of the two corresponding log-likelihoods.

@d ldpmvnorm
@{
ldpmvnorm <- function(obs, lower, upper, mean, invcholmean, chol, invchol, 
                      logLik = TRUE, ...) {

    if (missing(obs) || is.null(obs))
        return(lpmvnorm(lower = lower, upper = upper, 
                        mean = mean, invcholmean = invcholmean,
                        chol = chol, invchol = invchol, logLik = logLik, ...))
    if (missing(lower) && missing(upper) || is.null(lower) && is.null(upper))
        return(ldmvnorm(obs = obs, mean = mean, invcholmean = invcholmean,
                        chol = chol, invchol = invchol, logLik = logLik))

    @<dp input checks@>    

    stopifnot(missing(invcholmean))

    if (!missing(invchol)) {
        J <- dim(invchol)[2L]
        stopifnot(cJ + dJ == J)

        ### check if 1:cJ corresponds to obs
        nm <- dimnames(invchol)[[2L]]
        rno <- rownames(obs)
        rnl <- rownames(lower)
        rnu <- rownames(upper)
        if (!is.null(rnl) && !is.null(rnu))
            stopifnot(isTRUE(all.equal(rnl, rnu)))
        if (!is.null(rnl) && !is.null(rno) && !is.null(nm))
            stopifnot(isTRUE(all.equal(c(rno, rnl), nm)))

        md <- marg_mvnorm(invchol = invchol, which = 1:cJ)
        ret <- ldmvnorm(obs = obs, mean = cmean, invchol = md$invchol, 
                        logLik = logLik)

        cd <- cond_mvnorm(invchol = invchol, which_given = 1:cJ, 
                          given = obs - cmean, center = TRUE)
        ret <- ret + lpmvnorm(lower = lower, upper = upper, mean = dmean, 
                              invchol = cd$invchol, center = cd$center, 
                              logLik = logLik, ...)
        return(ret)
    }

    J <- dim(chol)[2L]
    stopifnot(cJ + dJ == J)

    md <- marg_mvnorm(chol = chol, which = 1:cJ)
    ret <- ldmvnorm(obs = obs, mean = cmean, chol = md$chol, logLik = logLik)

    cd <- cond_mvnorm(chol = chol, which_given = 1:cJ, 
                      given = obs - cmean, center = TRUE)
    ret <- ret + lpmvnorm(lower = lower, upper = upper, mean = dmean, 
                          chol = cd$chol, center = cd$center, 
                          logLik = logLik, ...)
    return(ret)
}
@}

The score function requires a little extra work. We start with the case when
\code{invchol} is given

@d sldpmvnorm invchol
@{
byrow_orig <- attr(invchol, "byrow")
invchol <- ltMatrices(invchol, byrow = TRUE)

J <- dim(invchol)[2L]
stopifnot(cJ + dJ == J)

md <- marg_mvnorm(invchol = invchol, which = 1:cJ)
cs <- sldmvnorm(obs = obs, mean = cmean, invchol = md$invchol, logLik = logLik)

obs_cmean <- obs - cmean
cd <- cond_mvnorm(invchol = invchol, which_given = 1:cJ, 
                  given = obs_cmean, center = TRUE)
ds <- slpmvnorm(lower = lower, upper = upper, mean = dmean, 
                center = cd$center, invchol = cd$invchol, 
                logLik = logLik, ...)

tmp0 <- solve(cd$invchol, ds$mean, transpose = TRUE)
tmp <- - tmp0[rep(1:dJ, each = cJ),,drop = FALSE] * 
         obs_cmean[rep(1:cJ, dJ),,drop = FALSE]

Jp <- nrow(unclass(invchol))
diag <- attr(invchol, "diag")
M <- as.array(ltMatrices(1:Jp, diag = diag, byrow = TRUE))[,,1]
ret <- matrix(0, nrow = Jp, ncol = ncol(obs))
M1 <- M[1:cJ, 1:cJ, drop = FALSE]
idx <- t(M1)[.ut(cJ, diag = diag)]
ret[idx,] <- Lower_tri(cs$invchol, diag = diag)

idx <- c(t(M[-(1:cJ), 1:cJ]))
ret[idx,] <- tmp

M3 <- M[-(1:cJ), -(1:cJ), drop = FALSE]
idx <- t(M3)[.ut(nrow(M3), diag = diag)]
ret[idx,] <- Lower_tri(ds$invchol, diag = diag)

ret <- ltMatrices(ret, diag = diag, byrow = TRUE)
if (!diag) diagonals(ret) <- 0
ret <- ltMatrices(ret, byrow = byrow_orig)

### post differentiate mean 
aL <- as.array(invchol)[-(1:cJ), 1:cJ,,drop = FALSE]
lst <- tmp0[rep(1:dJ, cJ),,drop = FALSE]
if (dim(aL)[3] == 1)
      aL <- aL[,,rep(1, ncol(lst)), drop = FALSE]
dim <- dim(aL)
dobs <- -margin.table(aL * array(lst, dim = dim), 2:3)

ret <- c(list(invchol = ret, obs = cs$obs + dobs), 
         ds[c("lower", "upper")])
if (logLik) ret$logLik <- cs$logLik + ds$logLik
ret$mean <- rbind(-ret$obs, ds$mean)
return(ret)
@}

For \code{chol}, we compute the above code for its inverse and
post-differentiate using the vec-trick

@d sldpmvnorm
@{
sldpmvnorm <- function(obs, lower, upper, mean, invcholmean, chol, invchol, 
                       logLik = TRUE, ...) {

    if (missing(obs) || is.null(obs))
        return(slpmvnorm(lower = lower, upper = upper, 
                         mean = mean, invcholmean = invcholmean, 
                         chol = chol, invchol = invchol, logLik = logLik, ...))
    if (missing(lower) && missing(upper) || is.null(lower) && is.null(upper))
        return(sldmvnorm(obs = obs, mean = mean, invcholmean = invcholmean,
                         chol = chol, invchol = invchol, logLik = logLik))

    @<dp input checks@>    

    ### for mixed data we cannot deal with invcholmean at the moment
    stopifnot(missing(invcholmean))

    if (!missing(invchol)) {
        @<sldpmvnorm invchol@>
    }

    invchol <- solve(chol)
    ret <- sldpmvnorm(obs = obs, lower = lower, upper = upper, 
                      mean = mean, invchol = invchol, logLik = logLik, ...)
    ### this means: ret$chol <- - vectrick(invchol, ret$invchol, invchol)
    ret$chol <- as.chol(- vectrick(invchol, ret$invchol))
    ret$invchol <- NULL
    return(ret)
}
@}

Let's assume we observed the first two dimensions exactly in our small
example, and the remaining two dimensions are only known in intervals. The
log-likelihood and score function for $\muvec$ and $\mC$ are 

<<ex-ML-cd>>=
ic <- 1:2 	### position of continuous variables
ll_cd <- function(parm, J) {
     m <- parm[1:J]             ### mean parameters
     parm <- parm[-(1:J)]       ### chol parameters
     C <- matrix(c(parm), ncol = 1L)
     C <- ltMatrices(C, diag = TRUE, byrow = BYROW)
     -ldpmvnorm(obs = Y[ic,], lower = lwr[-ic,], 
                upper = upr[-ic,], mean = m, chol = C, 
                w = W[-ic,,drop = FALSE], M = M)
}
sc_cd <- function(parm, J) {
    m <- parm[1:J]             ### mean parameters
    parm <- parm[-(1:J)]       ### chol parameters
    C <- matrix(c(parm), ncol = 1L)
    C <- ltMatrices(C, diag = TRUE, byrow = BYROW)
    ret <- sldpmvnorm(obs = Y[ic,], lower = lwr[-ic,],
                      upper = upr[-ic,], mean = m, chol = C, 
                      w = W[-ic,,drop = FALSE], M = M)
    return(-c(rowSums(ret$mean),
              rowSums(Lower_tri(ret$chol, diag = TRUE))))
}
@@
and the score function seems to be correct
<<ex-ML-cd-score>>=
if (require("numDeriv", quietly = TRUE))
    chk(grad(ll_cd, start, J = J), sc_cd(start, J = J), 
        check.attributes = FALSE, tolerance = 1e-6)
@@

We can now jointly estimate all model parameters via
<<ex-ML-cd-optim>>=
op <- optim(start, fn = ll_cd, gr = sc_cd, J = J, 
            method = "L-BFGS-B", lower = llim, 
            control = list(trace = FALSE))
## estimated C
ltMatrices(matrix(op$par[-(1:J)], ncol = 1), 
           diag = TRUE, byrow = BYROW)
## compare with true C
lt
## estimated means
op$par[1:J]
## compare with true means
mn
@@

The one restriction in both \code{ldpmvnorm} and \code{sldpmvnorm} is that the
continuous variables $\rY$ are ranked before the discrete variables $\rX$ in
the observation $(\rY_i, \rX_i)$, and thus also in $(\muvec, \etavec)$ and $\mC$
(the subscript $i$ is dropped from the parameters in the following paragraph
to keep the notational complexity in check).

While the means can be simply permuted, this is not the case for the
Cholesky factor $\mC$ (see function \code{aperm} in
Section~\ref{sec:conv}). Of course, we can simply permute $\hat{\mC}_i$, but
we loose standard errors in this process. Alternatively, we can permute the
order of variables in $\mC$ to our liking in the log-likelihood function (while
keeping the original order of the observations and for the mean parameters)

<<ex-ML-ap>>=
### discrete variables first
perm <- c((1:J)[-ic], ic)
ll_ap <- function(parm, J) {
     m <- parm[1:J]             ### mean parameters; NOT permuted
     parm <- parm[-(1:J)]       ### chol parameters
     C <- matrix(c(parm), ncol = 1L)
     C <- ltMatrices(C, diag = TRUE, byrow = BYROW)
     Ct <- aperm(as.chol(C), perm = perm)
     -ldpmvnorm(obs = Y[ic,], lower = lwr[-ic,], 
                upper = upr[-ic,], mean = m, chol = Ct, 
                w = W[-ic,,drop = FALSE], M = M)
}
@@

Unfortunately, this distorts the score function and we need to
``de-permute'' the scores. We start with $\mSigma = \mC \mC^\top$, the
Cholesky decomposition of a quadratic positive definite $\J \times \J$ covariance
matrix. There are $\J \times (\J + 1) / 2$ parameters in the lower
triagular part (including the diagonal) of $\mC$. Changing the order of the
variables by a permutation $\pi$ with permutation matrix $\Pi$ gives a
covariance $\Pi \mC \mC^\top \Pi^\top$. This is no longer a Cholesky
decomposition, because $\Pi \mC$ is not lower triangular. The new
decomposition is
\begin{eqnarray*}
\Pi \mC \mC^\top \Pi^\top = \tilde{\mC} \tilde{\mC}^\top
\end{eqnarray*}
($\tilde{\mC}$ is what \code{aperm} computes). As $\mC$, the Cholesky factor
$\tilde{\mC}$ is lower triangular with $\J \times (\J + 1) / 2$ parameters.
We could write this operation as a function
\begin{eqnarray*}
& & f_3: \R^{\J \times (\J + 1) / 2} \rightarrow \R^{\J \times (\J + 1) / 2} \\
& & f_3(\mC) = \tilde{\mC},
\end{eqnarray*}
where in fact $f_3 = $\code{aperm}, and we are interested in its gradient. Deriving the gradient of a Cholesky
decomposition might seem hopeless (it certainly did, at least to me, for a
very long time), but there is a trick. Let us define two other functions:
\begin{eqnarray*}
& & f_1: \R^{\J \times (\J + 1) / 2} \rightarrow \R^{\J \times \J} \\
& & f_1(\mC) = \Pi \mC \mC^\top \Pi^\top \\
& & f_2: \R^{\J \times (\J + 1) / 2} \rightarrow \R^{\J \times \J} \\
& & f_2(\tilde{\mC}) = \tilde{\mC} \tilde{\mC}^\top.
\end{eqnarray*}
Exploiting the chain rule for the composition $f_1 = f_2 \circ f_3$, 
we can write the gradient of $f_1$ as the product
of the gradients of $f_2$ and $f_3$:
\begin{eqnarray} \label{fm:chain}
\frac{\partial f_1(\mC)}{\partial \mC} = 
\frac{\partial f_2(\tilde{\mC})}{\partial \tilde{\mC}} \frac{\partial f_3(\mC)}{\partial \mC}.
\end{eqnarray}
The last factor is what we want to compute. It turns out that it is simpler
to compute the first two gradients first and, in a second step, to derive
the last factor. In more detail
\begin{eqnarray*}
\frac{\partial f_1(\mC)}{\partial \mC} & = & \frac{\partial \Pi \mC \mC^\top \Pi^\top}{\partial \mC} \\ 
& = & \frac{\partial \Pi \mC \mC^\top \Pi^\top}{\partial \Pi \mC} \frac{\partial \Pi \mC}{\mC} \\
& = & \left( (\Pi \mC \otimes \mI_\J) + (\mI_\J \otimes \Pi \mC) \frac{\partial \mA^\top}{\partial \mA} \right) (\mI_\J \otimes \Pi).
\end{eqnarray*}
($\mA$ is a quadratic matrix and the gradient of its transpose is a
permutation matrix). This analytic expression only contains known elements
and can be computed. The same applies to
\begin{eqnarray*}
\frac{\partial f_2(\tilde{\mC})}{\partial \tilde{\mC}} & = & \frac{\partial \tilde{\mC} \tilde{\mC}^\top \Pi}{\partial \tilde{\mC}} \\
&= & (\tilde{\mC} \otimes \mI_\J) + (\mI_\J \otimes \tilde{\mC}) \frac{\partial \mA^\top}{\partial \mA}
\end{eqnarray*}
Both expressions treat $\mC$ or $\tilde{\mC}$ as full matrices, we are only
interested in the score contributions by the $\J \times (\J + 1) / 2$ lower
triangular elements. Using sloppy notation, we collect the relevant columns
in matrices $\mB_1 = \frac{\partial f_1(\mC)}{\partial \mC} \in \R^{\J^2 \times \J \times (\J + 1) / 2}$
and $\mB_2 = \frac{\partial f_2(\tilde{\mC})}{\partial \tilde{\mC}} \in \R^{\J^2 \times \J \times (\J + 1) /
2}$. For the last, unknown, factor, we write $\mB_3 = \frac{\partial f_3(\tilde{\mC})}{\partial \tilde{\mC}} \in
\R^{\J \times (\J + 1) / 2 \times \J \times (\J + 1) / 2}$ and, with
formula~(\ref{fm:chain}), $\mB_1 = \mB_2 \mB_3$. We can then solve for
$\mB_3$ in the system $\mB_1^\top \mB_1 = \mB_1^\top \mB_2 \mB_3$.

With \code{chol} $ = \mC$, \code{permuted\_chol} $ = \tilde{\mC}$,
\code{perm} $ = \pi$ and score \code{score\_schol} of the log-likelihood $\ell(\tilde{\mC})$
with respect to the parameters in $\tilde{\mC}$, we can now implement this
de-permutation of the scores. Starting with some basic sanity checks, we
require lower triangular matrix objects as inputs, with diagonal elements,
and check if the dimensions match

@d deperma input checks chol
@{
stopifnot(is.ltMatrices(chol))	### NOTE: replace with is.chol
byrow_orig <- attr(chol, "byrow")
chol <- ltMatrices(chol, byrow = FALSE)
stopifnot(is.ltMatrices(permuted_chol))	### NOTE: replace with is.chol
permuted_chol <- ltMatrices(permuted_chol, byrow = FALSE)
stopifnot(max(abs(dim(chol) - dim(permuted_chol))) == 0)
J <- dim(chol)[2L]
stopifnot(attr(chol, "diag"))
INVCHOL <- !missing(invchol)
@}

Regarding \code{perm}, we check if this is an actual permutation 

@d deperma input checks perm	
@{
if (missing(perm)) return(score_schol)
stopifnot(isTRUE(all.equal(sort(perm), 1:J)))
if (max(abs(perm - 1:J)) == 0) return(score_schol)
@}

The scores with respect to $\tilde{\mC}$ have been computed elsewhere, we
just check the dimensions. In case we were given the scores with respect to
$\mL$, we first compute the scores with respect to $\mC$ (as we were lazy
and only derived the results for $\mC$). As in \code{standardize}, the
argument \code{score\_schol} gives the score with respect to $\mC$ and it is
the user's responsibility to provide this quantity (even when \code{invchol}
is given).

@d deperma input checks schol
@{
if (is.ltMatrices(score_schol)) { 
    byrow_orig_s <- attr(score_schol, "byrow")
    score_schol <- ltMatrices(score_schol, byrow = FALSE)
    ### don't do this here!
    ### if (INVCHOL) score_schol <- -vectrick(permuted_invchol, score_schol)
    score_schol <- unclass(score_schol) ### this preserves byrow
}
stopifnot(is.matrix(score_schol))
N <- ncol(score_schol)
stopifnot(J * (J + 1) / 2 == nrow(score_schol))
@}

We'll have to loop over $i = 1, \dots, N$ eventually and therefore coerce
all objects to objects of class \code{array}, there is no need to worry
about row or column storage order. We set-up indices matrices and the
permutation matrix $\Pi$

@d deperma indices
@{
idx <- matrix(1:J^2, nrow = J, ncol = J)	### assuming byrow = TRUE
tidx <- c(t(idx))
ltT <- idx[.lt(J, diag = TRUE)]
P <- matrix(0, nrow = J, ncol = J)
P[cbind(1:J, perm)] <- 1
ID <- diag(J)
IDP <- (ID %x% P)
@}

and are now ready for the main course. We are gentle and also allow
\code{invchol}$ = \mL$ as input, and we clean-up by post-differentiation at
the very end in this case.

@d deperma
@{
deperma <- function(chol = solve(invchol), 
                    permuted_chol = solve(permuted_invchol), 
                    invchol, permuted_invchol, perm, score_schol) {

    @<deperma input checks chol@>
    @<deperma input checks perm@>
    @<deperma input checks schol@>

    @<deperma indices@>

    Nc <- dim(chol)[1L]
    mC <- as.array(chol)[perm,,,drop = FALSE]
    Ct <- as.array(permuted_chol)
    ret <- lapply(1:Nc, function(i) {
        B1 <- (mC[,,i] %x% ID) + (ID %x% mC[,,i])[,tidx]
        #                                        ^^^^^^^ <- d t(A) / d A
        B1 <- B1 %*% IDP
        B1 <- B1[,ltT] ### relevant columns of B1
        B2 <- (Ct[,,i] %x% ID) + (ID %x% Ct[,,i])[,tidx]
        B2 <- B2[,ltT] ### relevant columns of B2
        B3 <- try(solve(crossprod(B2), crossprod(B2, B1)))
        if (inherits(B3, "try-error")) 
            stop("failure computing permutation score")
        if (Nc == 1L)
            return(crossprod(score_schol, B3))
        return(crossprod(score_schol[,i,drop = FALSE], B3))
    })
    ret <- do.call("rbind", ret)
    ret <- ltMatrices(t(ret), diag = TRUE, byrow = FALSE)
    if (INVCHOL)
        ret <- -vectrick(chol, ret)
    ret <- ltMatrices(ret, byrow = byrow_orig_s)
    return(ret)
}
@}

We can now use this function to estimate the Cholesky factor for $(\rX, \rY)$
when the data comes as $(\rY, \rX)$ (which is needed because continuous
variables come first in our implementation of log-likehood and score
function).

<<ex-ML-ap-score>>=
sc_ap <- function(parm, J) {
    m <- parm[1:J]               ### mean parameters; NOT permuted
    parm <- parm[-(1:J)]         ### chol parameters
    C <- matrix(c(parm), ncol = 1L)
    C <- ltMatrices(C, diag = TRUE, byrow = BYROW)
    ### permutation
    Ct <- aperm(as.chol(C), perm = perm)
    ret <- sldpmvnorm(obs = Y[ic,], lower = lwr[-ic,],
                      upper = upr[-ic,], mean = m, chol = Ct, 
                      w = W[-ic,,drop = FALSE], M = M)
    ### undo permutation for chol
    retC <- deperma(chol = C, permuted_chol = Ct, 
                   perm = perm, score_schol = ret$chol)
    return(-c(rowSums(ret$mean),
              rowSums(Lower_tri(retC, diag = TRUE))))
}
@@
and the score function seems to be correct
<<ex-ML-ap-grad>>=
if (require("numDeriv", quietly = TRUE))
    chk(grad(ll_ap, start, J = J), sc_ap(start, J = J), 
        check.attributes = FALSE, tolerance = 1e-6)
@@

We can now jointly estimate all model parameters via
<<ex-ML-ap-optim->>=
op <- optim(start, fn = ll_ap, gr = sc_ap, J = J, 
            method = "L-BFGS-B", lower = llim, 
            control = list(trace = FALSE))
## estimated C for (X, Y)
ltMatrices(matrix(op$par[-(1:J)], ncol = 1), 
           diag = TRUE, byrow = BYROW)
## compare with true _permuted_ C for (X, Y)
round(as.array(aperm(as.chol(lt), perm = perm)), 4)
@@


\chapter{Unstructured Gaussian Copula Estimation} \label{copula}

With $\rZ \sim \ND_\J(0, \mI_\J)$ and $\rY = \tilde{\mC} \rZ \sim \ND_\J(0, \tilde{\mC}
\tilde{\mC}^\top)$ we want to estimate the off-diagonal elements of the
lower triangular unit-diagonal matrix $\mC$. We have $\tilde{\mC}(\mC) := \diag(\mC \mC^\top)^{-\nicefrac{1}{2}} \mC$ 
such that $\mSigma = \tilde{\mC} \tilde{\mC}^\top$
is a correlation matrix ($\diag(\mSigma) = \mI_\J$). Note that directly
estimating $\tilde{\mC}$ requires $\J (\J + 1) / 2$ parameters under
constraints $\diag(\mSigma) = 1$ whereas only $\J (\J - 1) / 2$ parameters are necessary
when estimating the lower triangular part of $\mC$. The standardisation by
$\diag(\mC \mC^\top)^{-\nicefrac{1}{2}}$ ensures that $\diag(\mSigma)
\equiv 1$, that is, unconstained optimisation can be applied.

The function works both for unit and non-unit triangular matrices

@d standardize
@{
standardize <- function(chol, invchol) {
    stopifnot(xor(missing(chol), missing(invchol)))
    if (!missing(invchol)) {
        return(invcholD(invchol))
    }
    return(Dchol(chol))
}
@}

<<ex-stand>>=
### unit triangular
C <- ltMatrices(runif(10))
chk(as.array(chol2cov(standardize(chol = C))),
    as.array(chol2cor(standardize(chol = C))))
L <- solve(C)
chk(as.array(invchol2cov(standardize(invchol = L))),
    as.array(invchol2cor(standardize(invchol = L))))
### non-unit 
diagonals(C) <- 2
chk(as.array(chol2cov(standardize(chol = C))),
    as.array(chol2cor(standardize(chol = C))))
L <- solve(C)
chk(as.array(invchol2cov(standardize(invchol = L))),
    as.array(invchol2cor(standardize(invchol = L))))
@@

The log-likelihood function is $\ell_i(\mC_i)$ (we omit $i$ in the
following) and we assume the score
\begin{eqnarray*}
\frac{\partial \ell(\mC)}{\partial \mC}
\end{eqnarray*}
is already available. We want to compute the score
\begin{eqnarray*}
\frac{\partial \ell(\tilde{\mC})}{\partial \mC}
\end{eqnarray*}
which gives
\begin{eqnarray*}
\frac{\partial \ell(\tilde{\mC})}{\partial \mC} & = & 
\underbrace{\frac{\partial \ell(\tilde{\mC})}{\partial \tilde{\mC}}}_{=: \mT} \times \frac{\partial \tilde{\mC}(\mC)}{\partial \mC}
\end{eqnarray*}

We further have
\begin{eqnarray*}
\frac{\partial \tilde{\mC}(\mC)}{\partial \mC} = (\mC^\top \otimes \mI_\J)
\frac{\partial \diag(\mC \mC^\top)^{-\nicefrac{1}{2}}}{\partial \mC} +
(\mI_\J \otimes \diag(\mC \mC^\top)^{-\nicefrac{1}{2}})
\end{eqnarray*}
and thus
\begin{eqnarray*}
\frac{\partial \ell(\tilde{\mC})}{\partial \mC}
& = & 
\vecop(\mI_\J \mT \mC^\top)^\top \frac{\partial \diag(\mC \mC^\top)^{-\nicefrac{1}{2}}}{\partial \mC} + 
    \vecop(\diag(\mC \mC^\top)^{-\nicefrac{1}{2}} \mT \mI_\J)^\top
\end{eqnarray*}
and with 
\begin{eqnarray*}
\frac{\partial \diag(\mC \mC^\top)^{-\nicefrac{1}{2}}}{\partial \mC} & = & 
  \left. \frac{\partial \diag(\mA)^{-\nicefrac{1}{2}}}{\partial \mA} \right|_{\mA = \mC \mC^\top} \frac{\partial \mC \mC^\top}{\partial \mC} \\
& = & 
  -\frac{1}{2} \diag(\vecop(\diag(\mC \mC^\top)^{-\nicefrac{3}{2}})) \left[ (\mC \otimes \mI_\J) \frac{\partial \mC}{\partial \mC} + (\mI_\J \otimes \mC) \frac{\partial \mC^\top}{\partial \mC}\right]
\end{eqnarray*}
we can write
\begin{eqnarray*}
\vecop(\mI_\J \mT \mC^\top)^\top (-\frac{1}{2}) \diag(\vecop(\diag(\mC \mC^\top)^{-\nicefrac{3}{2}}))
& = & 
  -\frac{1}{2} \times \vecop(\mI_\J \mT \mC^\top)^\top \times \vecop(\diag(\mC \mC^\top)^{-\nicefrac{3}{2}})^\top \\
& =: & \bvec^\top
\end{eqnarray*}
thus
\begin{eqnarray*}
\frac{\partial \ell(\tilde{\mC})}{\partial \mC}
& = & 
\bvec^\top \left[ (\mC \otimes \mI_\J) \frac{\partial \mC}{\partial \mC} + (\mI_\J \otimes \mC) \frac{\partial \mC^\top}{\partial \mC}\right] 
  + \vecop(\diag(\mC \mC^\top)^{-\nicefrac{1}{2}} \mT \mI_\J)^\top \\
& = & \vecop(\mI_\J \mB \mC)^\top + \vecop(\mC^\top \mB \mI_\J)^\top \frac{\partial \mC^\top}{\partial \mC}
  + \vecop(\diag(\mC \mC^\top)^{-\nicefrac{1}{2}} \mT \mI_\J)^\top
\end{eqnarray*}
when $\bvec = \vecop(\mB)$. These scores are implemented in
\code{destandardize} with \code{chol} $ = \mC$ and \code{score\_schol} $= \mT$.
If the model was parameterised in $\mL = \mC^{-1}$, we have \code{invchol} $
= \mL$, however, we would still need to compute $\mT$ (\code{score\_schol}, the score with
respect to $\mC$, and it is the user's responsibility to provide this
quantity).

@d destandardize
@{
destandardize <- function(chol = solve(invchol), invchol, score_schol)
{
    stopifnot(is.ltMatrices(chol))	### NOTE: replace with is.chol
    J <- dim(chol)[2L]
    # stopifnot(!attr(chol, "diag"))
    byrow_orig <- attr(chol, "byrow")
    chol <- ltMatrices(chol, byrow = FALSE)
    
    ### TODO: check byrow in score_schol?

    if (is.ltMatrices(score_schol))
        score_schol <- matrix(as.array(score_schol), 
                              nrow = dim(score_schol)[2L]^2)
    stopifnot(is.matrix(score_schol))
    N <- ncol(score_schol)
    stopifnot(J^2 == nrow(score_schol))

    CCt <- Tcrossprod(chol, diag_only = TRUE)
    DC <- Dchol(chol, D = Dinv <- 1 / sqrt(CCt))
    SDC <- solve(DC)

    IDX <- t(M <- matrix(1:J^2, nrow = J, ncol = J))
    i <- cumsum(c(1, rep(J + 1, J - 1)))
    ID <- diagonals(as.integer(J), byrow = FALSE)
    if (dim(ID)[1L] != dim(chol)[1L])
        ID <- ID[rep(1, dim(chol)[1L]),]

    B <- vectrick(ID, score_schol, chol)
    B[i,] <- B[i,] * (-.5) * c(CCt)^(-3/2)
    B[-i,] <- 0

    Dtmp <- Dchol(ID, D = Dinv)

    ret <- vectrick(ID, B, chol, transpose = c(TRUE, FALSE)) +
           vectrick(chol, B, ID)[IDX,] +
           vectrick(Dtmp, score_schol, ID)

    if (!missing(invchol)) {
        ### this means: ret <- - vectrick(chol, ret, chol)
        ret <- - vectrick(chol, ret)
    }
    ret <- ret[M[.lt(J, diag = FALSE)],,drop = FALSE]
    if (!is.null(dimnames(chol)[[1L]]))
        colnames(ret) <- dimnames(chol)[[1L]]
    ret <- ltMatrices(ret,
                      diag = FALSE, byrow = FALSE, 
                      names = dimnames(chol)[[2L]])
    ret <- ltMatrices(ret, byrow = byrow_orig)
    diagonals(ret) <- 0
    return(ret)
}
@}

We can now set-up the log-likelihood and score functions for a Gaussian
copula model.  We start with the classical approach of generating the
marginal observations $\rY$ from the ECDF with denominator $N + 1$ and
subsequent use of the Lebesque density as likelihood.  Because no stats text
on multivariate problems is complete without a reference to Edgar Anderson's
iris data, let's set up a model for these four classical variables

<<gc-classical>>=
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)))
@@

This approach is of course a bit strange, because we estimate the marginal
distributions by nonparametric maximum likelihood whereas the joint
distribution is estimated by plain maximum likelihood. For the latter, we
can define the likelihood by boxes given by intervals obtained from the
marginale ECDFs and estimate the Copula parameters by maximisation of this
nonparametric likelihood.

<<gc-NPML>>=
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)
    -lpmvnorm(lower = lwr, upper = upr, chol = Cs, M = M, w = W)
}
sc <- function(parm) {
    C <- ltMatrices(parm)
    Cs <- standardize(chol = C)
    -rowSums(Lower_tri(destandardize(chol = C, 
        score_schol = slpmvnorm(lower = lwr, upper = upr, chol = Cs, 
                               M = M, w = W)$chol)))
}
if (require("numDeriv", quietly = TRUE))
    chk(grad(ll, start), sc(start), check.attributes = FALSE)
op2 <- optim(start, fn = ll, gr = sc, method = "BFGS", 
             control = list(trace = FALSE), hessian = TRUE)
S_NPML <- chol2cov(standardize(chol = ltMatrices(op2$par)))
@@

For $N = \Sexpr{nrow(iris)}$, the difference is (as expected) marginal:
<<gc>>=
S_ML
S_NPML
@@
with relatively close standard errors
<<gc-se>>=
sd_ML <- ltMatrices(sqrt(diag(solve(op$hessian))))
diagonals(sd_ML) <- 0
sd_NPML <- try(ltMatrices(sqrt(diag(solve(op2$hessian)))))
if (!inherits(sd_NPML, "try-error")) {
    diagonals(sd_NPML) <- 0
    print(sd_ML)
    print(sd_NPML)
}
@@

However, we are still in the wrong, because we pretended we know the
marginal distribution functions. A holistic approach would estimate the
marginal empirical distributions and the copula parameters simultaneously.

\chapter{Joint Mean and Covariance Estimation} \label{ch:joint}

From a computational points of view, the parameterisation of the
multivariate normal distribution in terms of the inverse Cholesky factor
$\mL$ and the scaled mean $\nuvec = \mL \muvec$ instead of the mean
$\muvec$ is attractive, 
because the log-density is then jointly concave in $\nuvec$ and $\mL$ and
thus a convex optimisation problem would emerge \citep{Barrathh_Boyd_2023}.
This also carries over to the interval-censored likelihood because
probabilities of log-concave densities are again log-concave
\citep{Prekopa_1973}.

The package implements the log-likelihood contributions as $\ell_i(\muvec_i,
\mL_i)$ and derives the scores with respect to $\muvec_i$ and $\mL_i$ (FIXME:
\code{ldmvnorm/sldmvnorm} now directly work with $\nuvec$, but the other
procedures still need adaptation). With
$\muvec_i = \mL_i^{-1} \nuvec$ it is simple to rewrite the log-likelihood
contribution as a concave function of both parameters as
\begin{eqnarray*}
\ell^\text{concave}_i(\nuvec_i, \mL_i) = \ell_i(\mL_i^{-1} \nuvec, \mL_i).
\end{eqnarray*}
The implementation of the corresponding convex negative log-likelihood is
also easy:

<<concave-L>>=
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; FIXME: update example using invcholmean args
d <- L %*% m

nll <- function(parm) {
    d <- parm[seq_len(J)]
    L <- ltMatrices(parm[-seq_len(J)], diag = TRUE)
    -ldmvnorm(obs = Y, mean = solve(L, d), invchol = L)
}

start <- c(d, prm)

nll(start)
### identical
-ldmvnorm(obs = Y, mean = m, invchol = L)
@@

However, computing the scores is a bit more work. Let's first consider the
scores with respect to $\nuvec_i$. We have
\begin{eqnarray*}
\frac{\partial \ell^\text{concave}_i(\nuvec_i, \mL_i)}{\partial \nuvec_i} & = & 
  \frac{\partial \ell_i(\mL_i^{-1} \nuvec_i, \mL_i)}{\partial \nuvec_i} \\
& = & \left.\frac{\partial \ell_i(\muvec, \mL_i)}{\partial
\muvec}\right|_{\muvec = \mL_i^{-1} \nuvec_i} \frac{\partial \mL_i^{-1}
\nuvec_i}{\partial \nuvec_i} \\
& = & \left.\frac{\partial \ell_i(\muvec, \mL_i)}{\partial
\muvec}\right|_{\muvec = \mL_i^{-1} \nuvec_i} \mL^{-1}_i.
\end{eqnarray*}
The first term is the score with respect to $\muvec_i$ and thus the score with
respect to $\nuvec_i$ is straightforward to compute.
The score with respect to $\mL_i$ is more complex. With the multivariate
chain rule we have
\begin{eqnarray*}
\frac{\partial \ell^\text{concave}_i(\nuvec_i, \mL_i)}{\partial \mL_i} & = & 
   \frac{\partial \ell_i(\mL_i^{-1} \nuvec_i, \mL_i)}{\partial \mL_i} \\
& = & \left.\frac{\partial \ell_i(\muvec, \mL_i)}{\partial
\muvec}\right|_{\muvec = \mL_i^{-1} \nuvec_i} \frac{\partial \mL_i^{-1}
\nuvec_i}{\partial \mL_i} + \left.\frac{\partial \ell_i(\muvec, \mL_i)}{\partial
\mL_i}\right|_{\muvec = \mL_i^{-1} \nuvec_i} \\
& = & \left.\frac{\partial \ell_i(\muvec, \mL_i)}{\partial
\muvec}\right|_{\muvec = \mL_i^{-1} \nuvec_i} \frac{\partial \mL_i^{-1}
\nuvec_i}{\partial \mL_i^{-1}} \frac{\partial \mL_i^{-1}}{\partial \mL_i} + 
\left.\frac{\partial \ell_i(\muvec, \mL_i)}{\partial
\mL_i}\right|_{\muvec = \mL_i^{-1} \nuvec_i} \\
& = & \underbrace{\left.\frac{\partial \ell_i(\muvec, \mL_i)}{\partial
\muvec}\right|_{\muvec = \mL_i^{-1} \nuvec_i}}_{=: \muvec_i^\prime} (\nuvec_i^\top \otimes
\mI_J)(-1)(\mL_i^{-\top} \otimes \mL_i^{-1}) + 
\left.\frac{\partial \ell_i(\muvec, \mL_i)}{\partial
\mL_i}\right|_{\muvec = \mL_i^{-1} \nuvec_i} \\
& = & \vecop(\muvec_i^\prime \nuvec_i^\top)(-1)(\mL_i^{-\top} \otimes \mL_i^{-1})
+ \left.\frac{\partial \ell_i(\muvec, \mL_i)}{\partial
\mL_i}\right|_{\muvec = \mL_i^{-1} \nuvec_i} 
\end{eqnarray*}
where we apply the vec-trick twice in the last line.

The negative gradient is now
<<scores-L>>=
nsc <- function(parm) {
    d <- parm[seq_len(J)]
    L <- ltMatrices(parm[-seq_len(J)], diag = TRUE)
    ret <- sldmvnorm(obs = Y, mean = solve(L, d), invchol = L)
    C <- solve(L)

    J <- dim(L)[2L]
    M <- matrix(seq_len(J^2), nrow = J, byrow = FALSE)
    idx <- M[lower.tri(M, diag = TRUE)]
   
    X <- -ret$obs
    Y <- matrix(d, nrow = nrow(X), ncol = ncol(X))
    A <- X[rep(1:nrow(X), times = nrow(X)),,drop = FALSE] * 
         Y[rep(1:nrow(Y), each = nrow(X)),,drop = FALSE]

    scL <- - vectrick(C, A)
    scL <- scL[idx,,drop = FALSE]
    scL <- rowSums(unclass(scL) + unclass(ret$invchol))
    - c(rowSums(solve(L, -ret$obs, transpose = TRUE)), 
        scL)
}
chk(unname(nsc(start)), grad(nll, start))
@@

Similarily, when the model is parameterized by $\nuvec_i$ and $\mC_i =
\mL_i^{-1}$, we can compute the log-likelihood as

<<concave-C>>=
C <- ltMatrices(prm <- runif(J * (J + 1) / 2), diag = TRUE)
Z <- matrix(rnorm(N * J), nrow = J)
Y <- C %*% Z + m
### scaled mean
d <- solve(C, m)

nll <- function(parm) {
    d <- parm[seq_len(J)]
    C <- ltMatrices(parm[-seq_len(J)], diag = TRUE)
    -ldmvnorm(obs = Y, mean = C %*% d, chol = C)
}

start <- c(d, prm)

nll(start)
### identical
-ldmvnorm(obs = Y, mean = m, chol = C)
@@

By the same arguments as above, the score with respect to $\nuvec_i$ is
\begin{eqnarray*}
\frac{\partial \ell^\text{concave}_i(\nuvec_i, \mC_i)}{\partial \nuvec_i} & = & 
  \frac{\partial \ell_i(\mC_i \nuvec_i, \mC_i)}{\partial \nuvec_i} \\
& = & \left.\frac{\partial \ell_i(\muvec, \mC_i)}{\partial
\muvec}\right|_{\muvec = \mC_i \nuvec_i} \mC_i.
\end{eqnarray*}
and the score with respect to $\mC_i$ is
\begin{eqnarray*}
\frac{\partial \ell^\text{concave}_i(\nuvec_i, \mC_i)}{\partial \mC_i} & = & 
   \frac{\partial \ell_i(\mC_i \nuvec_i, \mC_i)}{\partial \mC_i} \\
& = & \muvec_i^\prime \nuvec_i^\top + \left.\frac{\partial \ell_i(\muvec,
\mC_i)}{\partial \mC_i}\right|_{\muvec = \mC_i \nuvec_i} 
\end{eqnarray*}

<<scores-C>>=
nsc <- function(parm) {
    d <- parm[seq_len(J)]
    C <- ltMatrices(parm[-seq_len(J)], diag = TRUE)
    ret <- sldmvnorm(obs = Y, mean = C %*% d, chol = C)

    J <- dim(C)[2L]
    M <- matrix(seq_len(J^2), nrow = J, byrow = FALSE)
    idx <- M[lower.tri(M, diag = TRUE)]

    X <- -ret$obs
    Y <- matrix(d, nrow = nrow(X), ncol = ncol(X))
    A <- X[rep(1:nrow(X), times = nrow(X)),,drop = FALSE] * 
         Y[rep(1:nrow(Y), each = nrow(X)),,drop = FALSE]

    scC <- A[idx,,drop = FALSE]
    scC <- rowSums(unclass(scC) + unclass(ret$chol))
    - c(rowSums(crossprod(C, -ret$obs)), 
        scC)
}
chk(unname(nsc(start)), grad(nll, start))
@@



\chapter{User Interface} \label{inter}

@o interface.R -cp
@{
@<mvnorm@>
@<mean2invcholmean@>
@<invcholmean2mean@>
@<mvnorm methods@>
@<mvnorm simulate@>
@<mvnorm margDist@>
@<mvnorm condDist@>
@<mvnorm logLik@>
@<mvnorm lLgrad@>
@<mvnorm marginal means@>
@<mvnorm vcov@>
@<mvnorm regression coefs@>
@}

The tools provided in the previous chapters are rather low-level, so we will
invest some time into setting-up a more high-level interface for
representing normal models, either as $\ND_\J(\muvec, \mC \mC^\top)$, or
$\ND_\J(\muvec, \mL^{-1} \mL^{-\top})$, or $\ND_\J(\mL^{-1} \nuvec, \mL^{-1}
\mL^{-\top})$, or $\ND_\J(\mC \nuvec, \mC \mC^\top)$, for simulating from such models, and
for evaluating the log-likelihood and corresponding score functions. The
latter functionality shall also work when only incomplete (variables are
missing) or censored (observations are only known as intervals) data is
available.

We start with the conversion of a lower triangular matrix \code{x} to an
\code{ltMatrices} object

@d as.ltMatrices
@{
as.ltMatrices.default <- function(x) {
    stopifnot(is.numeric(x))
    if (!is.matrix(x)) x <- matrix(x)    
    DIAG <- max(abs(diag(x) - 1)) > .Machine$double.eps
    DIAG <- DIAG & (nrow(x) > 1)
    lt <- x[.lt(nrow(x), diag = DIAG)]
    up <- x[.ut(nrow(x), diag = FALSE)]
    stopifnot(max(abs(up)) < .Machine$double.eps)
    nm <- rownames(x)
    if (!is.null(nm))
        return(ltMatrices(lt, diag = DIAG, names = nm))
    return(ltMatrices(lt, diag = DIAG))
}
@}

and proceed defining a constructor for object respresenting, potentially
multiple, multivariate normal distributions. If the Cholesky factor $\mC$
(or multiple Cholesky factors $\mC_1, \dots, \mC_N$) are given as
\code{chol} argument, we label them as being such objects using \code{as.chol}. If 
only a matrix is given, we convert it (if possible) to a single Cholesky
factor $\mC$. The same is done when $\mL$ is given as \code{invchol}
argument. Of course, only one of these arguments must be specified.

@d mvnorm chol invchol
@{
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)
@}

The mean, or multiple means, is stored as a $\J \times 1$ or $\J \times N$
matrix, and we check if dimensions and, possibly, names are in line with
what was specified as \code{chol} or \code{invchol}

@d mvnorm mean
@{
if (!missing(mean)) {
    stopifnot(is.numeric(mean))
    stopifnot(NROW(mean) == dim(scale)[2L])
    if (!is.matrix(mean)) {
        mean <- matrix(mean, nrow = NROW(mean))
        rownames(mean) <- names(mean)
    }
    nm <- dimnames(scale)[[2L]]
    if (is.null(rownames(mean)))
        rownames(mean) <- nm
    if (!isTRUE(all.equal(rownames(mean), nm)))
        stop("rownames of mean do not match")        
    nm <- dimnames(scale)[[1L]]
    if (!is.null(nm) && dim(scale)[[2L]] == ncol(mean)) {
        if (is.null(colnames(mean)))
            colnames(mean) <- nm
        if (!isTRUE(all.equal(colnames(mean), nm)))
            stop("colnames of mean do not match")        
    }
    ret$mean <- mean
}
@}

@d mvnorm invcholmean
@{
if (!missing(invcholmean)) {
    stopifnot(is.numeric(invcholmean))
    stopifnot(NROW(invcholmean) == dim(scale)[2L])
    if (!is.matrix(invcholmean)) {
        invcholmean <- matrix(invcholmean, nrow = NROW(invcholmean))
        rownames(invcholmean) <- names(invcholmean)
    }
    nm <- dimnames(scale)[[2L]]
    if (is.null(rownames(invcholmean)))
        rownames(invcholmean) <- nm
    if (!isTRUE(all.equal(rownames(invcholmean), nm)))
        stop("rownames of invcholmean do not match")        
    nm <- dimnames(scale)[[1L]]
    if (!is.null(nm) && dim(scale)[[2L]] == ncol(invcholmean)) {
        if (is.null(colnames(invcholmean)))
            colnames(invcholmean) <- nm
        if (!isTRUE(all.equal(colnames(invcholmean), nm)))
            stop("colnames of invcholmean do not match")        
    }
    ret$invcholmean <- invcholmean
}
@}


Finally, we put everything together and return an object of class
\code{mvnorm}, featuring \code{mean} and \code{scale}. The class of the
latter slot carries the information how this object is to be interpreted (as
Cholesky factor or inverse thereof). We also allow the specification of the
scaled mean $\nuvec = \mL \muvec$ as argument \code{invcholmean}:

@d mvnorm
@{
### allow more than one distribution
mvnorm <- function(mean, invcholmean, chol, invchol) {

  @<mvnorm chol invchol@>
  if (!missing(invcholmean)) {
      stopifnot(missing(mean))
      ### note: this is not really necessary once all ls functions
      ### are able to deal with invcholmean
      @<mvnorm invcholmean@>
  } else {
      @<mvnorm mean@>
  }
  class(ret) <- "mvnorm"
  return(ret)
}
@}

After updating the object, it might be necessary to update $\nuvec$ or
$\muvec$ as well:

@d mean2invcholmean
@{
.m2i <- function(object) {
    if (is.chol(object$scale))
        return(solve(object$scale, object$mean))
    return(object$scale %*% object$mean)
}
@}

@d invcholmean2mean
@{
.i2m <- function(object) {
    if (is.chol(object$scale))
        return(object$scale %*% object$invcholmean)
    return(solve(object$scale, object$invcholmean))
}
@}

We add a \code{names} and \code{aperm} method. The latter returns a
multivariate normal distribution with permuted order of the variables

@d mvnorm methods
@{
names.mvnorm <- function(x)
    dimnames(x$scale)[[2L]]

aperm.mvnorm <- function(a, perm, ...) {

    if (!is.null(a$invcholmean))
        a$mean <- .i2m(a)
    ret <- list(scale = aperm(a$scale, perm = perm, ...))
    if (!is.null(a$mean))
        ret$mean <- a$mean[perm,,drop = FALSE]
    if (!is.null(a$invcholmean)) {
        ret$invcholmean <- .m2i(ret)
        ret$mean <- NULL
    }
    class(ret) <- "mvnorm"
    ret
}
@}

We are now ready to draw samples from such an object. If multiple normal
distributions are contained in \code{object}, we return one sample each,
otherwise, \code{nsim} samples are returned. Because most tools in this
package expect data as $\J \times N$ matrices, we return the data in this
format. If a classical \code{data.frame} is preferred, \code{as.data.frame =
TRUE} we provide one

@d mvnorm simulate
@{
simulate.mvnorm <- function(object, nsim = dim(object$scale)[1L], seed = NULL, 
                            standardize = FALSE, as.data.frame = FALSE, ...) {

    @<init random seed, reset on exit@>

    J <- dim(object$scale)[2L]
    N <- dim(object$scale)[1L]
    if (N > 1)
        stopifnot(nsim == N)
    if (standardize) {
        if (is.chol(object$scale)) {
            object$scale <- standardize(chol = object$scale)
        } else {
            object$scale <- standardize(invchol = object$scale)
        }
    }
    Z <- matrix(rnorm(nsim * J), nrow = J)
    if (!is.null(object$invcholmean))
        Z <- Z + c(object$invcholmean)
    if (is.chol(object$scale)) {
        Y <- Mult(object$scale, Z)
    } else {
        Y <- solve(object$scale, Z)
    }
    ret <- Y
    if (!is.null(object$mean))
        ret <- ret + c(object$mean)
    rownames(ret) <- dimnames(object$scale)[[2L]]
    if (!as.data.frame)
        return(ret)
    return(as.data.frame(t(ret)))
}
@}

It is maybe time for a first example, and we return to the iris dataset,
ignoring the iris' species for the time being. We set-up a model 
in terms of the sample maximum-likelihood estimates
<<iris-model>>=
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))
@@

Marginal and conditional distributions might be of interest, the
\code{margDist} and \code{condDist} methods are simple wrappers to
\code{marg\_mvnorm} and \code{cond\_mvnorm}

@d mvnorm margDist
@{
margDist <- function(object, which, ...)
    UseMethod("margDist")

margDist.mvnorm <- function(object, which, ...) {

    if (is.chol(object$scale)) {
        ret <- list(scale = as.chol(marg_mvnorm(chol = object$scale, 
                                                which = which)$chol))
    } else {
        ret <- list(scale = as.invchol(marg_mvnorm(invchol = object$scale, 
                                                   which = which)$invchol))
    }

    if (!is.null(object$invcholmean))
        object$mean <- .i2m(object)
    if (!is.null(object$mean))
        ret$mean <- object$mean[which,,drop = FALSE]
    if (!is.null(object$invcholmean)) {
        ret$invcholmean <- .m2i(ret)
        ret$mean <- NULL
    }

    class(ret) <- "mvnorm"
    return(ret)
}
@}

@d mvnorm condDist
@{
condDist <- function(object, which_given, given, ...)
    UseMethod("condDist")

condDist.mvnorm <- function(object, which_given = 1L, given, ...) {

    if (!is.null(object$invcholmean))
        object$mean <- .i2m(object)
    if (!is.null(object$mean)) {
        if (ncol(object$mean) == 1L && NCOL(given) > 1L) {
            given <- given - object$mean[which_given,,drop = TRUE]
        } else {
            given <- c(given) - object$mean[which_given,,drop = FALSE]
        }
    }

    if (is.chol(object$scale)) {
        ret <- cond_mvnorm(chol = object$scale, which_given = which_given, 
                           given = given, ...)
        ret$scale <- as.chol(ret$chol)
        ret$chol <- NULL
    } else {
        ret <- cond_mvnorm(invchol = object$scale, which_given = which_given, 
                           given = given, ...)
        ret$scale <- as.invchol(ret$invchol)
        ret$invchol <- NULL
    }
    if (!is.null(object$mean)) {
        if (is.character(which_given)) 
            which_given <- match(which_given, dimnames(object$scale)[[2L]])
        if (ncol(object$mean) > 1L && ncol(ret$mean) > 1) {
            ret$mean <- object$mean[-which_given,,drop = FALSE] + ret$mean
        } else if (ncol(object$mean) == 1L && ncol(ret$mean) > 1L) {
            ret$mean <- object$mean[-which_given,,drop = TRUE] + ret$mean
        } else {
            ret$mean <- object$mean[-which_given,,drop = FALSE] + c(ret$mean)
        }
    }
    if (!is.null(object$invcholmean)) {
        ret$invcholmean <- .m2i(ret)
        ret$mean <- NULL
    }
    class(ret) <- "mvnorm"
    return(ret)
}
@}

We could now compute the marginal distribution of two Petal variables
or the bivariate regressions of the two Petal variables given the observed
Sepal variables. Note that the last object contains $N = \Sexpr{nrow(iris)}$
different distributions

<<iris-mc>>=
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)
@@

We now work towards implementating the corresponding log-likelihood
function. This is a trivial task as long as all variables for all
observations have been observed exactly (that is, we can interpret 
the data as being continuous). Here, we also want to allow imprecise, that
is, interval-censored, measurements. The one constraint in \code{ldpmvnorm}
is that the continuous variables come first, followed by the censored ones.
This of course might not be in line with the variable ordering we have in
mind for our model. Our log-likelihood function shall be able to evaluate
the log-likelihood for arbitrary permutations of the variables and,
optionally, also based on marginal distributions in case observations are
missing.

The following \code{logLik} method for objects of class \code{mvnorm} is
essentially a wrapper for \code{ldpmvnorm}, handling permutations,
marginalisation, and standardisation. We begin with some sanity checks

@d argchecks
@{
args <- c(object, list(...))
nargs <- missing(obs) + missing(lower) + missing(upper)
stopifnot(nargs < 3L)

nmobs <- NULL
if (!missing(obs)) {
    if (!is.null(obs)) {
        stopifnot(is.matrix(obs))
        stopifnot(!is.null(rownames(obs)))
        nmobs <- rownames(obs)
    }
}
nmlower <- nmupper <- nmlu <- NULL
if (!missing(lower)) {
    if (!is.null(lower)) {
        stopifnot(is.matrix(lower))
        stopifnot(!is.null(rownames(lower)))
        nmlu <- nmlower <- rownames(lower)
    }
}
if (!missing(upper)) {
    if (!is.null(lower)) {
        stopifnot(is.matrix(upper))
        stopifnot(!is.null(rownames(upper)))
        nmupper <- rownames(upper)
        if (!missing(lower)) {
            stopifnot(isTRUE(all.equal(nmlower, nmupper)))
        } else {
            nmlu <- nmupper
        }
    }
}

nm <- c(nmobs, nmlu)
no <- names(object)
### we allow dimensions w/o data
stopifnot(all(nm %in% no))
perm <- NULL
if (!isTRUE(all.equal(nm, no)))
    perm <- c(nm, no[!no %in% nm])

if (!missing(obs)) args$obs <- obs
if (!missing(lower)) args$lower <- lower
if (!missing(upper)) args$upper <- upper
@}

and handle standardization for the two possible parameterisations, first
when $\mC$ was given and second when the inverse Cholesky factors $\mL$ was
given

@d logLik standardize
@{
if (is.chol(object$scale)) {
    names(args)[names(args) == "scale"] <- "chol"
    if (standardize)
       args$chol <- object$scale <- standardize(chol = args$chol)
} else {
    names(args)[names(args) == "scale"] <- "invchol"
    if (standardize)
        args$invchol <- object$scale <- standardize(invchol = args$invchol)
}
@}

Putting everything together in a corresponding \code{logLik} method finally
calls \code{ldpmvnorm} when there are (i) only continuous observations for
all dimensions, (ii) only intervals for all dimensions, (iii) continuous
observations for at least the first (maybe more) dimension followed by
intervals for the remaining dimensions. If this is not the case, \code{perm}
is not \code{NULL} and the appropriate marginal (for continuous) and conditional
distributions (intervals given continuous) distributions are computed first.
This also handles the case when data is missing for some observations.

In contrast to \code{lLgrad.mvnorm}, which first permutes the dimensions via
\code{aperm} and thus requires a Cholesky decomposition of the permuted $\J
\times \J$ covariance matrix, Cholesky decompositions of lower dimensional
matrices are required when computing the log-likelihood. However, it is best
to avoid this reshuffeling by using an appropriate order of the dimensions
(continuous first, intervals second, no missing data).

@d mvnorm logLik
@{
logLik.mvnorm <- function(object, obs, lower, upper, standardize = FALSE, 
                          ...) {
    @<argchecks@>

    @<logLik standardize@>

    if (!is.null(perm)) {
        ll <- if (!is.null(args$logLik)) args$logLik else TRUE
        ret <- 0
        ### integrate out dimensions w/o data
        if (length(nm) < length(no) && !is.null(nmlu))
            object <- margDist(object, which = nm)
        ### continuous
        if (!is.null(nmobs)) {
            ret <- ret + logLik(margDist(object, which = nmobs), 
                                obs = obs, logLik = ll) 
            ### set-up conditional distribution if there are any intervals
            if (!is.null(nmlu))
                object <- condDist(object, which_given = nmobs, given = obs)
        }
        ### interval given continuous
        if (!is.null(nmlu))
            ret <- ret + logLik(object, lower = lower, upper = upper, ...)
        return(ret)
    }

    ### base likelihood on mean for mixed data
    if (!is.null(args$invcholmean) &&
        (!missing(obs)) + (!missing(lower) || !missing(upper)) == 2L  ### mixed data
        ) {
        args$mean <- .i2m(object)
        args$invcholmean <- NULL
    }
    return(do.call("ldpmvnorm", args))
}
@}

allows us to evaluate the log-likelihood of the conditional models for iris

<<iris-ll>>=
logLik(object = iris_cmvn, obs = t(iris[,vars[-j]]))
@@

This implementation of the log-likelihood silently handles the case when
variables have been specified in a different order than hard-wired into the
model

<<iris-ll-perm>>=
logLik(object = iris_cmvn, obs = t(iris[,rev(vars[-j])]))
@@

The hardest task is the implementation of a score function which features
the same options as the log-likelihood function and provides the gradients
with respect not only to the parameters ($\muvec$ or $\nuvec$ and $\mC$ or $\mL$), but also
with respect to the data objects \code{obs}, \code{lower}, and \code{upper}.

In essence, we have to repair the damage imposed by a series of
transformations in \code{logLik.mvnorm}, that is, by standardisation,
permutation, and marginalisation. We start with the case when $\mC$ was
given. First, we repeat all the steps performed in \code{logLik}, but call 
the score function \code{sldpmvnorm} instead of the log-likelihood function
\code{ldpmvnorm}. A downside of the approach taken here is that
post-processing of the scores takes place on the level of observations, as
we do not compute the score function but the score contributions. If one is
only interested in the former, major speedups could be achieved by summing
over observations before post-processing.

@d lLgrad chol
@{
names(args)[names(args) == "scale"] <- "chol"
sc <- args$chol
if (standardize)
    args$chol <- sc <- standardize(chol = args$chol)
if (!is.null(perm)) {
    if (!attr(args$chol, "diag")) {
        diagonals(args$chol) <- 1
        sc <- args$chol
    }
    args$chol <- pc <- aperm(as.chol(args$chol), perm = perm)
    if (length(nm) < length(no))
        args$chol <- marg_mvnorm(chol = args$chol, which = nm)$chol
    args$mean <- args$mean[nm,,drop = FALSE]
}
ret <- do.call("sldpmvnorm", args)
@<lLgrad mean@>
@<lLgrad marginalisation@>
@<lLgrad deperma@>
@<lLgrad destandarized@>
@<lLgrad diagonals@>
@<lLgrad return@>
@}

The next task is to post-differentiate all scores such that the gradients
with respect to the original arguments of \code{logLik} are obtained. We
start with the gradient with respect to $\muvec$, in case it was not given

@d lLgrad mean
@{
### sldmvnorm returns mean score as -obs
if (is.null(ret$mean)) ret$mean <- - ret$obs
@}

In case we marginalised over some variables, we have to set the omitted
parameters to zero

@d lLgrad marginalisation
@{
om <- length(no) - length(nm)
if (om > 0) {
    am <- matrix(0, nrow = om, ncol = ncol(ret$mean))
    rownames(am) <- no[!no %in% nm]
    ret$mean <- rbind(ret$mean, am)
    Jo <- dim(object$scale)[[2L]]
    pJ <- dim(args$invchol)[[2L]]
    am <- matrix(0, nrow = Jo * (Jo + 1) / 2 - pJ * (pJ + 1) / 2, 
                 ncol = dim(ret$invchol)[1L])
    byrow_orig <- attr(ret$chol, "byrow")
    ret$chol <- ltMatrices(ret$chol, byrow = TRUE)
    ### rbind only works for byrow = TRUE
    ret$chol <- ltMatrices(rbind(unclass(ret$chol), am), 
                           byrow = TRUE, 
                           diag = TRUE,
                           names = perm)
    ret$chol <- ltMatrices(ret$chol, byrow = byrow_orig)
}
@}

If the order of the variables was permuted, we compute the scores for the
original ordering of the variables, as explained in Chapter~\ref{cdl}

@d lLgrad deperma
@{
if (!is.null(perm))
    ret$chol <- deperma(chol = sc, permuted_chol = pc, 
                        perm = match(perm, no), 
                        score_schol = ret$chol)
@}

The effect of standardization can be removed as discussed in
Chapter~\ref{copula}

@d lLgrad destandarized
@{
if (standardize)
    ret$chol <- destandardize(chol = object$scale, 
                              score_schol = ret$chol)
@}

and it remains to remove fix diagonal elements

@d lLgrad diagonals
@{
if (!attr(sc, "diag"))
    ret$chol <- ltMatrices(Lower_tri(ret$chol, diag = FALSE),
                           diag = FALSE, 
                           byrow = attr(ret$chol, "byrow"), 
                           names = dimnames(ret$chol)[[2L]])
@}

and to return the results, with mean scores in the correct ordering.
Eventually, we have to perform the post-processing steps explained in
Chapter~\ref{ch:joint} when the scaled mean $\nuvec = \mL \muvec$ was
specified and the score is not already available from the score workhorse
function \code{sldpmvnorm}

@d lLgrad return
@{
ret$scale <- ret$chol
ret$chol <- NULL
if (!is.null(ret$mean)) ret$mean <- ret$mean[no,, drop = FALSE]
### sldpmvnorm for mixed data returns mean score only
if (!is.null(object$invcholmean) && is.null(ret$invcholmean)) {
    J <- dim(sc)[2L]
    M <- matrix(seq_len(J^2), nrow = J, byrow = FALSE)
    idx <- M[.lt(J, diag = TRUE)]

    X <- ret$mean
    Y <- matrix(object$invcholmean, nrow = nrow(X), ncol = ncol(X))
    A <- X[rep(1:nrow(X), times = nrow(X)),,drop = FALSE] * 
         Y[rep(1:nrow(Y), each = nrow(X)),,drop = FALSE]

    scC <- ltMatrices(A[idx,,drop = FALSE], diag = TRUE, byrow = FALSE)
    if (!attr(sc, "diag"))
        scC <- ltMatrices(Lower_tri(scC, diag = FALSE), diag = FALSE, byrow = FALSE)
    scC <- ltMatrices(scC, byrow = attr(sc, "byrow"))
    ret$scale <- ltMatrices(unclass(scC) + unclass(ret$scale), 
                            diag = attr(sc, "diag"),
                            byrow = attr(sc, "byrow"),
                            names = dimnames(sc)[[2L]])
    ret$invcholmean <- Mult(sc, X, transpose = TRUE)
    ret$mean <- NULL
}
return(ret)
@}

The steps are essentially the same when $\mL$ was given, but we have to
post-differentiate $\mC = \mL^{-1}$ with respect to $\mL$. We start with
marginalisation

@d lLgrad invchol marginalisation
@{
om <- length(no) - length(nm)
if (om > 0) {
    am <- matrix(0, nrow = om, ncol = ncol(ret$mean))
    rownames(am) <- no[!no %in% nm]
    ret$mean <- rbind(ret$mean, am)
    Jo <- dim(object$scale)[[2L]]
    pJ <- dim(args$invchol)[[2L]]
    am <- matrix(0, nrow = Jo * (Jo + 1) / 2 - pJ * (pJ + 1) / 2, 
                 ncol = dim(ret$invchol)[1L])
    byrow_orig <- attr(ret$invchol, "byrow")
    ret$invchol <- ltMatrices(ret$invchol, byrow = TRUE)
    ### rbind only works for byrow = TRUE
    ret$invchol <- ltMatrices(rbind(unclass(ret$invchol), am), 
                              byrow = TRUE,
                              diag = TRUE,
                              names = perm)
    ret$invchol <- ltMatrices(ret$invchol, byrow = byrow_orig)
}
@}

turn to post-processing if necessary

@d lLgrad invchol post
@{
### sldpmvnorm for mixed data returns mean score only
if (!is.null(object$invcholmean) && is.null(ret$invcholmean)) {
    J <- dim(si)[2L]
    M <- matrix(seq_len(J^2), nrow = J, byrow = FALSE)
    idx <- M[.lt(J, diag = TRUE)]

    X <- ret$mean
    Y <- matrix(object$invcholmean, nrow = nrow(X), ncol = ncol(X))
    A <- X[rep(1:nrow(X), times = nrow(X)),,drop = FALSE] * 
         Y[rep(1:nrow(Y), each = nrow(X)),,drop = FALSE]

    scL <- - vectrick(solve(si), A)
    scL <- ltMatrices(scL[idx,,drop = FALSE], diag = TRUE, byrow = FALSE)
    if (!attr(si, "diag"))
        scL <- ltMatrices(Lower_tri(scL, diag = FALSE), diag = FALSE, byrow = FALSE)
    scL <- ltMatrices(scL, byrow = attr(si, "byrow"))
    ret$scale <- ltMatrices(unclass(scL) + unclass(ret$scale), 
                              diag = attr(si, "diag"),
                              byrow = attr(si, "byrow"),
                              names = dimnames(si)[[2L]])
    ret$invcholmean <- solve(si, X, transpose = TRUE)
    ret$mean <- FALSE
}
@}

and finally put everything together

@d lLgrad invchol
@{
names(args)[names(args) == "scale"] <- "invchol"
si <- args$invchol
if (standardize)
    args$invchol <- si <- standardize(invchol = args$invchol)
if (!is.null(perm)) {
    if (!attr(args$invchol, "diag")) {
        diagonals(args$invchol) <- 1
        si <- args$invchol
    }
    args$invchol <- pi <- aperm(as.invchol(args$invchol), perm = perm)
    if (length(nm) < length(no))
        args$invchol <- marg_mvnorm(invchol = args$invchol,
                                    which = nm)$invchol
    args$mean <- args$mean[nm,,drop = FALSE]
}
ret <- do.call("sldpmvnorm", args)
### sldmvnorm returns mean score as -obs
### return only if object had mean specified
if (is.null(ret$mean) && !is.null(args$mean)) 
    ret$mean <- - ret$obs

@<lLgrad invchol marginalisation@>

if (!is.null(perm))
    ret$invchol <- deperma(invchol = si, permuted_invchol = pi, 
                           perm = match(perm, no), 
                           score_schol = -vectrick(pi, ret$invchol))
if (standardize)
    ret$invchol <- destandardize(invchol = object$scale, 
                                 score_schol = -vectrick(si, ret$invchol))
if (!attr(si, "diag"))
    ret$invchol <- ltMatrices(Lower_tri(ret$invchol, diag = FALSE),
                              diag = FALSE, 
                              byrow = attr(ret$invchol, "byrow"), 
                              names = dimnames(ret$invchol)[[2L]])
ret$scale <- ret$invchol
ret$invchol <- NULL
if (!is.null(ret$mean)) ret$mean <- ret$mean[no,, drop = FALSE]
@<lLgrad invchol post@>
return(ret)
@}

We can now provide the log-likelihood gradients

@d mvnorm lLgrad
@{
lLgrad <- function(object, ...)
    UseMethod("lLgrad")

lLgrad.mvnorm <- function(object, obs, lower, upper, standardize = FALSE, 
                          ...) {
    @<argchecks@>

    ### base likelihood on mean when necessary
    ### note that post processing the mean score is time consuming
    if (!is.null(args$invcholmean) &&
        (!is.null(perm) || 	### permutations	
         (!missing(obs)) + (!missing(lower) || !missing(upper)) == 2L  ### mixed data
        )) {
        args$mean <- .i2m(object)
        args$invcholmean <- NULL
    }

    if (is.chol(object$scale)) {
        @<lLgrad chol@>
    }
    @<lLgrad invchol@>
}
@}

Let's use this infrastructure to set-up maximum-likelihood estimation
procedures. We start implementing the log-likelihood and score functions for
the iris dataset

<<iris-lLgrad>>=
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)))
}
@@

and can now estimate the mean and Cholesky factor of the covariance matrix.
Before we start, we check if the gradient, evaluated at the sample
maximum-likelihood estimates, is actually zero.

<<iris-ML>>=
### don't start at the solution
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
op <- optim(start, fn = ll, gr = sc, method = "L-BFGS-B", 
            lower = llim, 
            control = list(trace = FALSE, 
                           factr = 1e-6)) ### noLD machines
Chat <- ltMatrices(op$par[-(1:J)], diag = TRUE, names = vars)
ML <- mvnorm(mean = op$par[1:J], chol = Chat)
@@

Quite unsurprisingly, the results are practically equivalent to the
analytically available maximum-likelihood estimators in this case

<<iris-ML-hat>>=
### covariance (noLD brings up small differences)
round(vcov(ML), 3)
V
### mean
mean(ML)[,,drop = TRUE]
m
@@

We now check if we can obtain corresponding results when jointly optimising
the scaled mean $\nuvec$ and the inverse Cholesky factor $\mL$:

<<iris-lLgrad-nu>>=
ll <- function(parm, logLik = TRUE) {
    L <- ltMatrices(parm[-(1:J)], diag = TRUE, names = vars)
    x <- mvnorm(invcholmean = parm[1:J], invchol = L)
    if (!logLik) return(x)
    -logLik(object = x, obs = obs)
}
sc <- function(parm) {
    x <- ll(parm, logLik = FALSE)
    ret <- lLgrad(object = x, obs = obs)
    -c(rowSums(ret$invcholmean), rowSums(Lower_tri(ret$scale, diag = TRUE)))
}
### note: This is a convex problem now, so (here incorrect) 
### starting values shouldn't matter
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)
@@

where the comparison to the analytic estimates is


<<iris-ML-hat-nu>>=
### log-likelihood
op$value
opL$value
### covariance
round(vcov(MLL), 3)
V
### mean
mean(MLL)[,,drop = TRUE]
m
@@

Now, this was a lot of work to replace \code{mean} and \code{var} with
something more fancy, and we would of course not go down this way in real
life. But how about a more complex situation where one (or more) variables
are only known up to intervals? Let's present the first variable is such a
case

<<iris-interval>>=
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]
@@

The only necessary modification to our code is the specification of
\code{lower} and \code{upper} bounds for these intervals, and the removal of
the first variable from the ``exact continuous'' observations \code{obs}.
The rest of the machinery \emph{doesn't need any update at all}. Note that
the mean and covariance parameters are no longer orthogonal (as in the toy
example above), so we need to optimise over both sets of parameters
simultaneously and therefore opt for the parameterisation in terms of
$\nuvec$ and $\mL$ because this leads to a convex optimisation problem also
for interval-censored observations as mentioned in Chapter~\ref{ch:joint}.

One additional tweak is necessary for dealing with extreme parameter values
assigning a likelihood of numerically zero to some observations. In this
case, the score function is not defined (because we divide the derivative by
the likelihood). The corresponding values should be \code{NA} and removed
from the summation.

<<iris-MLi>>=
ll <- function(parm, logLik = TRUE) {
    L <- ltMatrices(parm[-(1:J)], diag = TRUE, names = vars)
    x <- mvnorm(invcholmean = parm[1:J], invchol = L)
    if (!logLik) return(x)
    -logLik(object = x, obs = obs, lower = lower, upper = upper, 
            tol = 1e-6) ### probs < tol are considered 0
}
sc <- function(parm) {
    x <- ll(parm, logLik = FALSE)
    ret <- lLgrad(object = x, obs = obs, lower = lower, upper = upper, 
                  tol = 1e-6) ### probs < tol are considered 0
    -c(rowSums(ret$invcholmean, na.rm = TRUE), 
       rowSums(Lower_tri(ret$scale, diag = TRUE), na.rm = TRUE))
}
@@

<<iris-MLi-opt>>=
start <- round(opL$par, 2)
if (require("numDeriv", quietly = TRUE))
    chk(grad(ll, start), sc(start), check.attributes = FALSE)
opi <- optim(start, fn = ll, gr = sc, method = "L-BFGS-B", 
             lower = llim, control = list(trace = FALSE, factr = 1e-6))
MLi <- ll(opi$par, logLik = FALSE)
@@

Because the likelihood is a product of a continuous density and a
conditional probability as introduced in Chapter~\ref{cdl}, the two
in-sample log-likelihoods are not comparable. However, the parameters of the
two estimated normal distributions can be compared directly (and are rather
close in our case)

<<iris-MLi-hat>>=
### covariance
round(vcov(MLi), 3)
round(vcov(MLL), 3)
round(vcov(ML), 3)
### mean
mean(MLi)[,,drop = TRUE]
mean(MLL)[,,drop = TRUE]
mean(ML)[,,drop = TRUE]
@@

Again, we can extract the regression coefficients from these multivariate
normal distributions. We start with the $\muvec$ and $\mC$ parameterisation:
<<iris-lm>>=
### least-squares coefficients
coef(irislm <- lm(Petal.Width ~ Sepal.Length + Sepal.Width + Petal.Length, 
                  data = iris))
### residual standard deviation
summary(irislm)$sigma
### compare with 
round(coef(ML, which = "Petal.Width"), 3)
@@
The regression coefficients and the residual variance are almost identical with
those obtained from \code{lm()}. The same exericise, now with the 
$\nuvec$ and $\mL$ parameterisation for exact and interval-censored
observations, gives slightly different results for the latter, 
due to censoring:
<<echo = FALSE>>=
# slightly different results on noLD machines
cat("> ## IGNORE_RDIFF_BEGIN\n")
@@
<<iris-lm-iL>>=
### nu, L for exact observations
round(coef(MLL, which = "Petal.Width"), 3)
### nu, L for censored observations
round(coef(MLi, which = "Petal.Width"), 3)
@@
<<echo = FALSE>>=
## slightly different results on noLD machines
cat("> ## IGNORE_RDIFF_END\n")
@@

The log-likelihood and score function automagically marginalise over
dimensions where all observations are $(-\infty, \infty)$. We can simply
omit these dimensions from the matrices specified as \code{obs}, \code{lower}, and \code{upper}
arguments. Let's say we have four dimensions called $A$ to $D$ and three
observations. All observations have $A = (-1, 1)$ and $B = (-\infty,
\infty)$, so in fact, the likelihood is given by the marginal distribution
of $A, C, D$.
<<marginB>>=
N <- 3
J <- 4
L <- ltMatrices(runif(J * (J + 1) / 2), diag = TRUE, names = LETTERS[1:J])
Z <- matrix(rnorm(J * N), nrow = J)
Y <- solve(L, Z)

lwrA <- matrix(-1, nrow = 1, ncol = N)
uprA <- matrix(1, nrow = 1, ncol = N)
rownames(lwrA) <- rownames(uprA) <- "A"

lwrB <- matrix(-Inf, nrow = 1, ncol = N)
uprB <- matrix(Inf, nrow = 1, ncol = N)
rownames(lwrB) <- rownames(uprB) <- "B"

lwr <- rbind(lwrA, lwrB)
upr <- rbind(uprA, uprB)
obs <- Y[rev(LETTERS[3:J]),]    ### change order of dimensions
@@
With this data, we first compute the log-likelihood and score functions for
the complete data, that is, including the infinite intervals for $B$.
<<marginBllsc>>=
w <- matrix(runif(1000), nrow = 1, byrow = TRUE)
lABCD <- logLik(mvnorm(invchol = L), obs = obs, lower = lwr, upper = upr, w = w)
sABCD <- lLgrad(mvnorm(invchol = L), obs = obs, lower = lwr, upper = upr, w = w)
@@
This is (almost) the same as omitting dimension $B$ from the data, but of
course not from the model
<<marginllsc>>=
lACD <- logLik(mvnorm(invchol = L), obs = obs, lower = lwrA, upper = uprA)
sACD <- lLgrad(mvnorm(invchol = L), obs = obs, lower = lwrA, upper = uprA)
@@
We can compare the results
<<marginchk>>=
chk(lABCD, lACD)
nm <- names(sABCD)
nm <- nm[!nm %in% c("lower", "upper")]
chk(sABCD[nm], sACD[nm])
@@
noting that the scores with respect to the $B$ data in \code{lower} and
\code{upper} are missing from \code{sACD}
<<marginsc>>=
chk(sABCD$lower["A",,drop = FALSE], sACD$lower)
chk(sABCD$upper["A",,drop = FALSE], sACD$upper)
sABCD$lower["B",]	### zero
sABCD$upper["B",]	### zero
@@
Omitting dimensions might be important because \code{lpmvnorm}
introduced in Chapter~\ref{lpmvnorm} does not check if both \code{lower} and
\code{upper} are infinite and omission thus reduces the dimensionality of
the integral we evaluate numerically. 

We close this chapter with a word of warning: If more than one variable is
censored, the \code{M} and \code{w} arguments to \code{lpmvnorm} and
\code{slpmvnorm} have to be specified in \code{logLik} and \code{lLgrad} as
additional arguments (\code{...}) \emph{AND MUST BE IDENTICAL} in both calls.


\chapter{Reduced Rank Covariance Matrices}

We sometimes can write the $\J \times \J$ covariance matrix as $\mSigma = \mB \mB^\top +
\mD$ where $\mB$ is $\J \times \K$ with $\K < \J$ and $\mD$ is a
$\J$-dimensional diagonal matrix. \cite{Marsaglia_1963} and
\cite{Genz_Bretz_2009}, in their Chapter 2.3.1, demonstrated that the
probability $\Prob(\avec < \rY \le \bvec \mid \mSigma)$ 
can be written as
\begin{eqnarray*}
\Prob(\avec < \rY \le \bvec \mid \mSigma ) 
  = \int_{[0,1]^\K} \prod_{j = 1}^J 
      \left(\Phi\left(\frac{b_j - \mB_{jk} \Phi^{-1}(w_k)}{\sqrt{\mD_{jj}}}\right) - 
            \Phi\left(\frac{a_j - \mB_{jk} \Phi^{-1}(w_k)}{\sqrt{\mD_{jj}}}\right)\right) 
    \, d \wvec
\end{eqnarray*}
where the integration takes place with respect to $\wvec = (w_1, \dots,
w_K)^\top$ from the $\K$- (and not $\J$-) dimensional unit hypercube.

We start implementing low-level functionality for computing log-probabilities for such structures with some book
keeping

@d RR input B D
@{
stopifnot(!missing(B))
if (!is.matrix(B)) B <- matrix(B, ncol = 1)
J <- nrow(B)
K <- ncol(B)
Dsqrt <- sqrt(D)
@}

We use $Z = \Phi^{-1}(w)$ and optional weights and compute the products $\mB_{jk}
\Phi^{-1}(w_k)$ for all $j$ and $k$ with standardisation by the diagonal
elements of $\mD$.

@d RR input Z, weights
@{
stopifnot(nrow(Z) == K)
stopifnot(length(weights) == 1 || length(weights) == ncol(Z))
BZ <- (B / Dsqrt) %*% Z
@}

The limits $\avec$

@d RR input lower
@{
if (missing(lower)) {
    pl <- 0
} else {
    stopifnot(length(lower) == J)
    lower <- c(lower)
    lower <- (lower - mean) / Dsqrt
    pl <- pnorm(lBZ <- lower - BZ)
}
@}

and $\bvec$

@d RR input upper
@{
if (missing(upper)) {
    pl <- 0
} else {
    stopifnot(length(upper) == J)
    upper <- c(upper)
    upper <- (upper - mean) / Dsqrt
    pu <- pnorm(uBZ <- upper - BZ)
}
@}

are processed and we finally compute the integrant, making sure to avoid
negative values. We
first compute log-probabilities, compute the sums and exponentiate before
summing up. In contrast to \code{lpmvnorm}, we also allow weights for the
summation, such that sparse grids (for example from add-on package \pkg{SparseGrid})
can be utilised.

@d RR inner
@{
inner <- pu - pl
inner <- pmax(0, inner)
retw <-  weights * exp(.colSums(m = J, n = ncol(Z), 
                                x = log(inner), na.rm = TRUE))
@}

Finally, we wrap everything up in a function called \code{lpRR}

@d lpRR
@{
lpRR <- function(lower, upper, mean = 0, B, D = rep(1, nrow(B)), 
                 Z, weights = 1 / ncol(Z), log.p = TRUE) {

    @<RR input B D@>
    @<RR input Z, weights@>
    @<RR input lower@>
    @<RR input upper@>

    @<RR inner@>
    ret <- sum(retw)
    if (log.p) return(log(max(c(0, ret))))
    return(ret)
}
@}

We test this functionality for dimensions $\J = 6$ and $\K = 3$. That is,
the integration problem reduces from a $\J - 1 = 5$ dimensional one in
\code{lpmvnorm} to a three dimensional one in \code{lpRR}. We first compare the
two log-probabilities, computed with high accuracy
<<RR-ll>>=
J <- 6
K <- 3
B <- matrix(rnorm(J * K), nrow = J)
D <- runif(J)
S <- tcrossprod(B) + diag(D)
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)
@@

The score function with respect to $\avec$, $\bvec$, $\mB$, and $\mD$ is 

@d slpRR
@{
slpRR <- function(lower, upper, mean = 0, B, D = rep(1, nrow(B)), 
                  Z, weights = 1 / ncol(Z), log.p = TRUE) {
    @<RR input B D@>
    @<RR input Z, weights@>
    @<RR input lower@>
    @<RR input upper@>

    @<RR inner@>
 
    dlBZ <- dnorm(lBZ)
    duBZ <- dnorm(uBZ)

    d <- matrix(retw, nrow = nrow(B), ncol = ncol(Z), byrow = TRUE) / inner
    db <- d * (duBZ - dlBZ)
    tdb <- t(db)
    dB <- -1 * do.call("cbind", lapply(1:nrow(Z), 
                       function(r) colSums(tdb * Z[r,], na.rm = TRUE))) / Dsqrt
    Du <- -.5 / D * uBZ
    Dl <- -.5 / D * lBZ
    dD <-  rowSums(d * duBZ * Du, na.rm = TRUE) - 
           rowSums(d * dlBZ * Dl, na.rm = TRUE)
    dl <- -rowSums(d * dlBZ, na.rm = TRUE) / Dsqrt
    du <-  rowSums(d * duBZ, na.rm = TRUE) / Dsqrt
    dm <- -du - dl # Dinb %*% -rowSums(d * (duBZ - dlBZ))
    fct <- 1
    if (log.p) fct <- 1 / sum(retw)
    list(lower = fct * dl, upper = fct * du, mean = fct * dm, 
         B = fct * dB, D = fct * dD)
}
@}

We can now compare the gradients for $\avec$, $\bvec$ and the mean in our
small example. We cannot expect them to be equal but close
<<RR-sc>>=
smv <- slpmvnorm(lower = a, upper = b, chol = Linv, w = w)
sRR <- slpRR(lower = a, upper = b, B = B, D = D, Z = Z)
chk(c(smv$lower), sRR$lower, tolerance = 1e-2)
chk(c(smv$upper), sRR$upper, tolerance = 1e-2)
chk(c(smv$mean), sRR$mean, tolerance = 1e-2 * 2)
@@

The gradient with respect to $\mB$ and $\mD$ are finally checked against their
numerical approximations
<<RR-sc-BD>>=
Z <- matrix(rnorm(K * 1000), nrow = K)
lB <- function(B) lpRR(lower = a, upper = b, B = B, D = D, Z = Z)
gB <- grad(lB, B)
sRR <- slpRR(lower = a, upper = b, B = B, D = D, Z = Z)
chk(gB, c(sRR$B), tolerance = 1e-3)
lD <- function(D) lpRR(lower = a, upper = b, B = B, D = D, Z = Z)
gD <- grad(lD, D)
chk(gD, c(sRR$D), tolerance = 1e-3)
### while we are at it, check lower and again
llwr <- function(a) lpRR(lower = a, upper = b, B = B, D = D, Z = Z)
glwr <- grad(llwr, a)
chk(glwr, c(sRR$lower))
lupr <- function(b) lpRR(lower = a, upper = b, B = B, D = D, Z = Z)
gupr <- grad(lupr, b)
chk(gupr, c(sRR$upper))
@@


\chapter{Package Infrastructure}

@d R Header
@{
###    Copyright (C) 2022- Torsten Hothorn
###
###    This file is part of the 'mvtnorm' R add-on package.
###
###    'mvtnorm' is free software: you can redistribute it and/or modify
###    it under the terms of the GNU General Public License as published by
###    the Free Software Foundation, version 2.
###
###    'mvtnorm' is distributed in the hope that it will be useful,
###    but WITHOUT ANY WARRANTY; without even the implied warranty of
###    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
###    GNU General Public License for more details.
###
###    You should have received a copy of the GNU General Public License
###    along with 'mvtnorm'.  If not, see <http://www.gnu.org/licenses/>.
###
###
###    DO NOT EDIT THIS FILE
###
###    Edit 'lmvnorm_src.w' and run 'nuweb -r lmvnorm_src.w'
@}

@d C Header
@{
/*
    Copyright (C) 2022- Torsten Hothorn

    This file is part of the 'mvtnorm' R add-on package.

    'mvtnorm' is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, version 2.

    'mvtnorm' is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with 'mvtnorm'.  If not, see <http://www.gnu.org/licenses/>.


    DO NOT EDIT THIS FILE

    Edit 'lmvnorm_src.w' and run 'nuweb -r lmvnorm_src.w'
*/
@}

\chapter*{Appendix}

This document uses the following matrix derivatives for vectors $\yvec \in
\R^J$ and matrices $\mA \in \R^{J \times J}$, $\mB \in \R^{K \times J}$, and 
$\rX \in \R^{J \times K}$:
\begin{eqnarray*}
\frac{\partial \mA \yvec}{\partial \yvec} & = & \mA \\
\frac{\partial \mA \yvec}{\partial \mA} & = & \yvec \otimes \mI_J \\
\frac{\partial \yvec^\top \mA^\top \mA \yvec}{\partial \mA} & = & 2 \mA \yvec \yvec^\top \\
\frac{\partial \mA^{-1}}{\partial \mA} & = & -(\mA^{-\top} \otimes \mA^{-1}) \\
\frac{\partial \mA \mA^\top}{\partial \mA} & = & (\mA \otimes \mI_J) \frac{\partial \mA}{\partial \mA} + (\mI_J \otimes \mA) \frac{\partial \mA^\top}{\partial \mA}
\\
& = & (\mA \otimes \mI_J) + (\mI_J \otimes \mA) \frac{\partial \mA^\top}{\partial \mA} \\
\frac{\partial \diag(\mA)}{\partial \mA} & = & \diag(\vecop(\mI_J)) \\
\frac{\partial \mA}{\partial \mA} & = & \diag(\mI_{J^2}) \\
\frac{\partial \yvec^\top \mA \yvec}{\partial \yvec} & = & \yvec^\top (\mA + \mA^\top) \\
\frac{\partial \mB \mA}{\partial \mA} & = & (\mI_J \otimes \mB)
\end{eqnarray*}
and the ``vec-trick'' $\vecop(\rX)^\top (\mB \otimes \mA^\top) = \vecop(\mA
\rX \mB)^\top$.


\chapter*{Index}

\section*{Files}

@f

\section*{Fragments}

@m

%\section*{Identifiers}
%
%@u

\bibliographystyle{jss}

<<bib, echo = FALSE>>=
thisdir <- getwd()
bibfile <- system.file("REFERENCES.bib", package = "mvtnorm")
### bibfile may contain spaces LaTeX is unable to deal with on MacOS it seems
if (file.copy(bibfile, to = thisdir, overwrite = TRUE)) {
    bibfile <- "REFERENCES.bib"
} else {
    ### hope for the best
    bibfile <- file.path("..", "inst", "REFERENCES.bib")
}
@@
\bibliography{\Sexpr{bibfile}}

\end{document}
