R version 3.2.3 (2015-12-10) -- "Wooden Christmas-Tree"
Copyright (C) 2015 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> pkgname <- "mvtnorm"
> source(file.path(R.home("share"), "R", "examples-header.R"))
> options(warn = 1)
> base::assign(".ExTimings", "mvtnorm-Ex.timings", pos = 'CheckExEnv')
> base::cat("name\tuser\tsystem\telapsed\n", file=base::get(".ExTimings", pos = 'CheckExEnv'))
> base::assign(".format_ptime",
+ function(x) {
+   if(!is.na(x[4L])) x[1L] <- x[1L] + x[4L]
+   if(!is.na(x[5L])) x[2L] <- x[2L] + x[5L]
+   options(OutDec = '.')
+   format(x[1L:3L], digits = 7L)
+ },
+ pos = 'CheckExEnv')
> 
> ### * </HEADER>
> library('mvtnorm')
> 
> base::assign(".oldSearch", base::search(), pos = 'CheckExEnv')
> cleanEx()
> nameEx("Mvnorm")
> ### * Mvnorm
> 
> flush(stderr()); flush(stdout())
> 
> base::assign(".ptime", proc.time(), pos = "CheckExEnv")
> ### Name: Mvnorm
> ### Title: Multivariate Normal Density and Random Deviates
> ### Aliases: dmvnorm rmvnorm
> ### Keywords: distribution multivariate
> 
> ### ** Examples
> 
> dmvnorm(x=c(0,0))
[1] 0.1591549
> dmvnorm(x=c(0,0), mean=c(1,1))
[1] 0.05854983
> 
> sigma <- matrix(c(4,2,2,3), ncol=2)
> x <- rmvnorm(n=500, mean=c(1,2), sigma=sigma)
> colMeans(x)
[1] 0.9492868 1.9916602
> var(x)
         [,1]     [,2]
[1,] 4.276371 2.105874
[2,] 2.105874 3.186058
> 
> x <- rmvnorm(n=500, mean=c(1,2), sigma=sigma, method="chol")
> colMeans(x)
[1] 0.8724731 1.9804160
> var(x)
         [,1]     [,2]
[1,] 4.804847 2.430442
[2,] 2.430442 3.148344
> 
> plot(x)
> 
> 
> 
> base::assign(".dptime", (proc.time() - get(".ptime", pos = "CheckExEnv")), pos = "CheckExEnv")
> base::cat("Mvnorm", base::get(".format_ptime", pos = 'CheckExEnv')(get(".dptime", pos = "CheckExEnv")), "\n", file=base::get(".ExTimings", pos = 'CheckExEnv'), append=TRUE, sep="\t")
> cleanEx()
> nameEx("Mvt")
> ### * Mvt
> 
> flush(stderr()); flush(stdout())
> 
> base::assign(".ptime", proc.time(), pos = "CheckExEnv")
> ### Name: Mvt
> ### Title: The Multivariate t Distribution
> ### Aliases: dmvt rmvt
> ### Keywords: distribution multivariate
> 
> ### ** Examples
> 
> ## basic evaluation
> dmvt(x = c(0,0), sigma = diag(2))
[1] -1.837877
> 
> ## check behavior for df=0 and df=Inf
> x <- c(1.23, 4.56)
> mu <- 1:2
> Sigma <- diag(2)
> x0 <- dmvt(x, delta = mu, sigma = Sigma, df = 0) # default log = TRUE!
> x8 <- dmvt(x, delta = mu, sigma = Sigma, df = Inf) # default log = TRUE!
> xn <- dmvnorm(x, mean = mu, sigma = Sigma, log = TRUE)
> stopifnot(identical(x0, x8), identical(x0, xn))
> 
> ## X ~ t_3(0, diag(2))
> x <- rmvt(100, sigma = diag(2), df = 3) # t_3(0, diag(2)) sample
> plot(x)
> 
> ## X ~ t_3(mu, Sigma)
> n <- 1000
> mu <- 1:2
> Sigma <- matrix(c(4, 2, 2, 3), ncol=2)
> set.seed(271)
> x <- rep(mu, each=n) + rmvt(n, sigma=Sigma, df=3)
> plot(x)
> 
> ## Note that the call rmvt(n, mean=mu, sigma=Sigma, df=3) does *not*
> ## give a valid sample from t_3(mu, Sigma)! [and thus throws an error]
> try(rmvt(n, mean=mu, sigma=Sigma, df=3))
Error in rmvt(n, mean = mu, sigma = Sigma, df = 3) : 
  Providing 'mean' does *not* sample from a multivariate t distribution!
> 
> ## df=Inf correctly samples from a multivariate normal distribution
> set.seed(271)
> x <- rep(mu, each=n) + rmvt(n, sigma=Sigma, df=Inf)
> set.seed(271)
> x. <- rmvnorm(n, mean=mu, sigma=Sigma)
> stopifnot(identical(x, x.))
> 
> 
> 
> base::assign(".dptime", (proc.time() - get(".ptime", pos = "CheckExEnv")), pos = "CheckExEnv")
> base::cat("Mvt", base::get(".format_ptime", pos = 'CheckExEnv')(get(".dptime", pos = "CheckExEnv")), "\n", file=base::get(".ExTimings", pos = 'CheckExEnv'), append=TRUE, sep="\t")
> cleanEx()
> nameEx("pmvnorm")
> ### * pmvnorm
> 
> flush(stderr()); flush(stdout())
> 
> base::assign(".ptime", proc.time(), pos = "CheckExEnv")
> ### Name: pmvnorm
> ### Title: Multivariate Normal Distribution
> ### Aliases: pmvnorm
> ### Keywords: distribution
> 
> ### ** Examples
> 
> 
> n <- 5
> mean <- rep(0, 5)
> lower <- rep(-1, 5)
> upper <- rep(3, 5)
> corr <- diag(5)
> corr[lower.tri(corr)] <- 0.5
> corr[upper.tri(corr)] <- 0.5
> prob <- pmvnorm(lower, upper, mean, corr)
> print(prob)
[1] 0.5800051
attr(,"error")
[1] 0.0002696831
attr(,"msg")
[1] "Normal Completion"
> 
> stopifnot(pmvnorm(lower=-Inf, upper=3, mean=0, sigma=1) == pnorm(3))
> 
> a <- pmvnorm(lower=-Inf,upper=c(.3,.5),mean=c(2,4),diag(2))
> 
> stopifnot(round(a,16) == round(prod(pnorm(c(.3,.5),c(2,4))),16))
> 
> a <- pmvnorm(lower=-Inf,upper=c(.3,.5,1),mean=c(2,4,1),diag(3))
> 
> stopifnot(round(a,16) == round(prod(pnorm(c(.3,.5,1),c(2,4,1))),16))
> 
> # Example from R News paper (original by Genz, 1992):
> 
> m <- 3
> sigma <- diag(3)
> sigma[2,1] <- 3/5
> sigma[3,1] <- 1/3
> sigma[3,2] <- 11/15
> pmvnorm(lower=rep(-Inf, m), upper=c(1,4,2), mean=rep(0, m), corr=sigma)
[1] 0.8279847
attr(,"error")
[1] 2.658133e-07
attr(,"msg")
[1] "Normal Completion"
> 
> # Correlation and Covariance
> 
> a <- pmvnorm(lower=-Inf, upper=c(2,2), sigma = diag(2)*2)
> b <- pmvnorm(lower=-Inf, upper=c(2,2)/sqrt(2), corr=diag(2))
> stopifnot(all.equal(round(a,5) , round(b, 5)))
> 
> 
> 
> 
> base::assign(".dptime", (proc.time() - get(".ptime", pos = "CheckExEnv")), pos = "CheckExEnv")
> base::cat("pmvnorm", base::get(".format_ptime", pos = 'CheckExEnv')(get(".dptime", pos = "CheckExEnv")), "\n", file=base::get(".ExTimings", pos = 'CheckExEnv'), append=TRUE, sep="\t")
> cleanEx()
> nameEx("pmvt")
> ### * pmvt
> 
> flush(stderr()); flush(stdout())
> 
> base::assign(".ptime", proc.time(), pos = "CheckExEnv")
> ### Name: pmvt
> ### Title: Multivariate t Distribution
> ### Aliases: pmvt
> ### Keywords: distribution
> 
> ### ** Examples
> 
> 
> n <- 5
> lower <- -1
> upper <- 3
> df <- 4
> corr <- diag(5)
> corr[lower.tri(corr)] <- 0.5
> delta <- rep(0, 5)
> prob <- pmvt(lower=lower, upper=upper, delta=delta, df=df, corr=corr)
> print(prob)
[1] 0.5063832
attr(,"error")
[1] 0.0002426557
attr(,"msg")
[1] "Normal Completion"
> 
> pmvt(lower=-Inf, upper=3, df = 3, sigma = 1) == pt(3, 3)
upper 
 TRUE 
> 
> # Example from R News paper (original by Edwards and Berry, 1987)
> 
> n <- c(26, 24, 20, 33, 32)
> V <- diag(1/n)
> df <- 130
> C <- c(1,1,1,0,0,-1,0,0,1,0,0,-1,0,0,1,0,0,0,-1,-1,0,0,-1,0,0)
> C <- matrix(C, ncol=5)
> ### scale matrix
> cv <- C %*% V %*% t(C)
> ### correlation matrix
> dv <- t(1/sqrt(diag(cv)))
> cr <- cv * (t(dv) %*% dv)
> delta <- rep(0,5)
> 
> myfct <- function(q, alpha) {
+   lower <- rep(-q, ncol(cv))
+   upper <- rep(q, ncol(cv))
+   pmvt(lower=lower, upper=upper, delta=delta, df=df,
+        corr=cr, abseps=0.0001) - alpha
+ }
> 
> ### uniroot for this simple problem
> round(uniroot(myfct, lower=1, upper=5, alpha=0.95)$root, 3)
[1] 2.561
> 
> # compare pmvt and pmvnorm for large df:
> 
> a <- pmvnorm(lower=-Inf, upper=1, mean=rep(0, 5), corr=diag(5))
> b <- pmvt(lower=-Inf, upper=1, delta=rep(0, 5), df=rep(300,5),
+           corr=diag(5))
Warning in if (df < 0) stop("cannot compute multivariate t distribution with ",  :
  the condition has length > 1 and only the first element will be used
Warning in if (isInf(df)) df <- 0 :
  the condition has length > 1 and only the first element will be used
> a
[1] 0.4215702
attr(,"error")
[1] 0
attr(,"msg")
[1] "Normal Completion"
> b
[1] 0.4211423
attr(,"error")
[1] 2.31377e-06
attr(,"msg")
[1] "Normal Completion"
> 
> stopifnot(round(a, 2) == round(b, 2))
> 
> # correlation and scale matrix
> 
> a <- pmvt(lower=-Inf, upper=2, delta=rep(0,5), df=3,
+           sigma = diag(5)*2)
> b <- pmvt(lower=-Inf, upper=2/sqrt(2), delta=rep(0,5),
+           df=3, corr=diag(5))
> attributes(a) <- NULL
> attributes(b) <- NULL
> a
[1] 0.5653944
> b
[1] 0.5654
> stopifnot(all.equal(round(a,3) , round(b, 3)))
> 
> a <- pmvt(0, 1,df=10)
> attributes(a) <- NULL
> b <- pt(1, df=10) - pt(0, df=10)
> stopifnot(all.equal(round(a,10) , round(b, 10)))
> 
> 
> 
> 
> base::assign(".dptime", (proc.time() - get(".ptime", pos = "CheckExEnv")), pos = "CheckExEnv")
> base::cat("pmvt", base::get(".format_ptime", pos = 'CheckExEnv')(get(".dptime", pos = "CheckExEnv")), "\n", file=base::get(".ExTimings", pos = 'CheckExEnv'), append=TRUE, sep="\t")
> cleanEx()
> nameEx("qmvnorm")
> ### * qmvnorm
> 
> flush(stderr()); flush(stdout())
> 
> base::assign(".ptime", proc.time(), pos = "CheckExEnv")
> ### Name: qmvnorm
> ### Title: Quantiles of the Multivariate Normal Distribution
> ### Aliases: qmvnorm
> ### Keywords: distribution
> 
> ### ** Examples
> 
> qmvnorm(0.95, sigma = diag(2), tail = "both")
$quantile
[1] 2.236358

$f.quantile
[1] -1.681424e-06

attr(,"message")
[1] "Normal Completion"
> 
> 
> 
> base::assign(".dptime", (proc.time() - get(".ptime", pos = "CheckExEnv")), pos = "CheckExEnv")
> base::cat("qmvnorm", base::get(".format_ptime", pos = 'CheckExEnv')(get(".dptime", pos = "CheckExEnv")), "\n", file=base::get(".ExTimings", pos = 'CheckExEnv'), append=TRUE, sep="\t")
> cleanEx()
> nameEx("qmvt")
> ### * qmvt
> 
> flush(stderr()); flush(stdout())
> 
> base::assign(".ptime", proc.time(), pos = "CheckExEnv")
> ### Name: qmvt
> ### Title: Quantiles of the Multivariate t Distribution
> ### Aliases: qmvt
> ### Keywords: distribution
> 
> ### ** Examples
> 
> ## basic evaluation
> qmvt(0.95, df = 16, tail = "both")
$quantile
[1] 2.119905

$f.quantile
[1] 0.975

> 
> ## check behavior for df=0 and df=Inf
> Sigma <- diag(2)
> set.seed(29)
> q0 <- qmvt(0.95, sigma = Sigma, df = 0,   tail = "both")$quantile
> set.seed(29)
> q8 <- qmvt(0.95, sigma = Sigma, df = Inf, tail = "both")$quantile
> set.seed(29)
> qn <- qmvnorm(0.95, sigma = Sigma, tail = "both")$quantile
> stopifnot(identical(q0, q8),
+           isTRUE(all.equal(q0, qn, tol = (.Machine$double.eps)^(1/3))))
> 
> ## if neither sigma nor corr are provided, corr = 1 is used internally
> df <- 0
> set.seed(29)
> qt95 <- qmvt(0.95, df = df, tail = "both")$quantile
> set.seed(29)
> qt95.c <- qmvt(0.95, df = df, corr  = 1, tail = "both")$quantile
> set.seed(29)
> qt95.s <- qmvt(0.95, df = df, sigma = 1, tail = "both")$quantile
> stopifnot(identical(qt95, qt95.c),
+           identical(qt95, qt95.s))
> 
> df <- 4
> set.seed(29)
> qt95 <- qmvt(0.95, df = df, tail = "both")$quantile
> set.seed(29)
> qt95.c <- qmvt(0.95, df = df, corr  = 1, tail = "both")$quantile
> set.seed(29)
> qt95.s <- qmvt(0.95, df = df, sigma = 1, tail = "both")$quantile
> stopifnot(identical(qt95, qt95.c),
+           identical(qt95, qt95.s))
> 
> 
> 
> base::assign(".dptime", (proc.time() - get(".ptime", pos = "CheckExEnv")), pos = "CheckExEnv")
> base::cat("qmvt", base::get(".format_ptime", pos = 'CheckExEnv')(get(".dptime", pos = "CheckExEnv")), "\n", file=base::get(".ExTimings", pos = 'CheckExEnv'), append=TRUE, sep="\t")
> ### * <FOOTER>
> ###
> options(digits = 7L)
> base::cat("Time elapsed: ", proc.time() - base::get("ptime", pos = 'CheckExEnv'),"\n")
Time elapsed:  1.044 0.012 1.059 0 0 
> grDevices::dev.off()
null device 
          1 
> ###
> ### Local variables: ***
> ### mode: outline-minor ***
> ### outline-regexp: "\\(> \\)?### [*]+" ***
> ### End: ***
> quit('no')