--- title: "condTruncMVN: Conditional Truncated Multivariate Normal Distribution" author: "Paul M. Hargarten" date: "`r Sys.Date()`" keywords: geometry: margin=1in preamble: > \usepackage{indentfirst} \usepackage{amsmath} \usepackage{graphicx} output: rmarkdown::html_vignette: toc: true vignette: > %\VignetteIndexEntry{condTruncMVN: Conditional Truncated Multivariate Normal Distribution} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} link-citations: true citation_package: biblatex bibliography: mnvignette.bib csl: multidisciplinary-digital-publishing-institute.csl # csl: /Users/Shared/Zotero/styles/multidisciplinary-digital-publishing-institute.csl --- ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, # If TRUE, all output would be in the code chunk. comment = "#>", fig.path = "man/figures/README-", out.width = "100%", results = "markup", prompt = TRUE, strip.white = TRUE, tidy = TRUE, tidy.opts = list(width.cutoff = 90), # options for tidy to remove blank lines [blank = FALSE] and set the approximate line width to be 80. fig.show = "asis", fig.height = 4.5, # inches fig.width = 4.5 # inches ) library("formatR") ``` The goal of condTruncMVN is to find densities, probabilities, and samples from a conditional truncated multivariate normal distribution. Suppose that **Z = (X,Y)** is from a fully-joint multivariate normal distribution of dimension *n* with **mean** $\boldsymbol\mu$ and covariance matrix **sigma** ( _$\Sigma$_ ) truncated between **lower** and **upper**. Then, Z has the density $$ f_Z(\textbf{z}, \boldsymbol\mu, \Sigma, \textbf{lower}, \textbf{upper})= \frac{exp(-\frac{1}{2}*(\textbf{z}-\boldsymbol\mu)^T \Sigma^{-1} (\textbf{z}-\boldsymbol\mu))} {\int_{\textbf{lower}}^{\textbf{upper}} exp(-\frac{1}{2}*(\textbf{z}-\boldsymbol\mu)^T \Sigma^{-1} (\textbf{z}-\boldsymbol\mu)) d\textbf{z}} $$ for all **z** in [**lower**, **upper**] in $\mathbb{R^{n}}$. This package computes the conditional truncated multivariate normal distribution of Y|X. The conditional distribution follows a truncated multivariate normal [@horraceResultsMultivariateTruncated2005]. Specifically, the functions are arranged such that $$ Y = Z[ , dependent.ind] $$ $$ X = Z[ , given.ind] $$ $$ Y|X = X.given \sim MVN(mean, sigma, lower, upper) $$ The [d,p,r]cmvtnorm() functions create a list of parameters used in truncated conditional normal and then passes the parameters to the source function below. | Function Name | Description | Source Function | Univariate Case| Additional Parameters |---------- | ------------| ---------- | ------------ | --------------------- | | condtMVN | List of parameters used in truncated conditional normal.| condMVNorm:: condMVN() | | | | dcmvtnorm | Calculates the density f(Y=y\| X = X.given) up to a constant. The integral of truncated distribution is not computed. | tmvtnorm:: dtmvnorm() | truncnorm:: dtruncnorm() | y, log | | pcmvtnorm | Calculates the probability that Y\|X is between lowerY and upperY given the parameters. | tmvtnorm:: ptmvnorm() | truncnorm:: ptruncnorm() | lowerY, upperY, maxpts, abseps, releps | | rcmvtnorm | Generate random sample. | tmvmixnorm:: rtmvn() | truncnorm:: rtruncnorm() | n, init, burn, thin | ## Installation You can install the released version of condTruncMVN from [CRAN](https://CRAN.R-project.org) with `install.packages("condTruncMVN")`. You can load the package by: ```{r} library("condTruncMVN") ``` And the development version from [GitHub](https://github.com/) with: ## Example Suppose $X2,X3,X5|X1 = 1,X4 = -1 \sim N_3(1, Sigma, -10, 10)$. The following code finds the parameters of the distribution, calculates the density, probability, and finds random variates from this distribution. ```{r example} library(condTruncMVN) d <- 5 rho <- 0.9 Sigma <- matrix(0, nrow = d, ncol = d) Sigma <- rho^abs(row(Sigma) - col(Sigma)) Sigma ``` First, we find the conditional Truncated Normal Parameters. ```{r} condtMVN(mean = rep(1, d), sigma = Sigma, lower = rep(-10, d), upper = rep(10, d), dependent.ind = c(2, 3, 5), given.ind = c(1, 4), X.given = c(1, -1) ) ``` Find the log-density when X2,X3,X5 all equal $0$: ```{r} dcmvtruncnorm( rep(0, 3), mean = rep(1, 5), sigma = Sigma, lower = rep(-10, 5), upper = rep(10, d), dependent.ind = c(2, 3, 5), given.ind = c(1, 4), X.given = c(1, -1), log = TRUE ) ``` Find $P( -0.5 < X2,X3,X5 < 0 | X1 = 1,X4 = -1)$: ```{r} pcmvtruncnorm( rep(-0.5, 3), rep(0, 3), mean = rep(1, d), sigma = Sigma, lower = rep(-10, d), upper = rep(10, d), dependent.ind = c(2, 3, 5), given.ind = c(1, 4), X.given = c(1, -1) ) ``` Generate two random numbers from the distribution. ```{r} set.seed(2342) rcmvtruncnorm(2, mean = rep(1, d), sigma = Sigma, lower = rep(-10, d), upper = rep(10, d), dependent.ind = c(2, 3, 5), given.ind = c(1, 4), X.given = c(1, -1) ) ``` Another Example: To find the probability that $X1|X2, X3, X4, X5 \sim N(**1**, Sigma, **-10**, **10**)$ is between -0.5 and 0: ```{r} pcmvtruncnorm(-0.5, 0, mean = rep(1, d), sigma = Sigma, lower = rep(-10, d), upper = rep(10, d), dependent.ind = 1, given.ind = 2:5, X.given = c(1, -1, 1, -1) ) ``` If I want to generate 2 random variates from $X1|X2, X3, X4, X5 \sim N(**1**, Sigma, **-10**, **10**)$: ```{r} set.seed(2342) rcmvtruncnorm(2, mean = rep(1, d), sigma = Sigma, lower = rep(-10, d), upper = rep(10, d), dependent.ind = 1, given.ind = 2:5, X.given = c(1, -1, 1, -1) ) ``` ## Computational Details This vignette is successfully processed using the following. ```{r echo=FALSE} # sessioninfo::session_info() # makes a mess! Instead cat(" -- Session info ---------------------------------------------------") sessioninfo::platform_info() cat("-- Packages -------------------------------------------------------") tmp.df <- sessioninfo::package_info( c("condMVNorm", "matrixNormal", "tmvmixnorm", "tmvtnorm", "truncnorm"), dependencies = FALSE ) print(tmp.df) ``` ## Final Notes [![Lifecycle::experimental](https://img.shields.io/badge/lifecycle-experimental-orange.svg)](https://www.tidyverse.org/lifecycle/#experimental) [![CRAN status](https://www.r-pkg.org/badges/version/condTruncMVN)](https://CRAN.R-project.org/package=condTruncMVN) ## References