R : Copyright 2005, The R Foundation for Statistical Computing Version 2.1.1 (2005-06-20), ISBN 3-900051-07-0 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. 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 a HTML browser interface to help. Type 'q()' to quit R. > ### *
> ### > attach(NULL, name = "CheckExEnv") > assign(".CheckExEnv", as.environment(2), pos = length(search())) # base > ## add some hooks to label plot pages for base and grid graphics > setHook("plot.new", ".newplot.hook") > setHook("persp", ".newplot.hook") > setHook("grid.newpage", ".gridplot.hook") > > assign("cleanEx", + function(env = .GlobalEnv) { + rm(list = ls(envir = env, all.names = TRUE), envir = env) + RNGkind("default", "default") + set.seed(1) + options(warn = 1) + delayedAssign("T", stop("T used instead of TRUE"), + assign.env = .CheckExEnv) + delayedAssign("F", stop("F used instead of FALSE"), + assign.env = .CheckExEnv) + sch <- search() + newitems <- sch[! sch %in% .oldSearch] + for(item in rev(newitems)) + eval(substitute(detach(item), list(item=item))) + missitems <- .oldSearch[! .oldSearch %in% sch] + if(length(missitems)) + warning("items ", paste(missitems, collapse=", "), + " have been removed from the search path") + }, + env = .CheckExEnv) > assign("..nameEx", "__{must remake R-ex/*.R}__", env = .CheckExEnv) # for now > assign("ptime", proc.time(), env = .CheckExEnv) > grDevices::postscript("mvtnorm-Examples.ps") > assign("par.postscript", graphics::par(no.readonly = TRUE), env = .CheckExEnv) > options(contrasts = c(unordered = "contr.treatment", ordered = "contr.poly")) > options(warn = 1) > library('mvtnorm') > > assign(".oldSearch", search(), env = .CheckExEnv) > assign(".oldNS", loadedNamespaces(), env = .CheckExEnv) > cleanEx(); ..nameEx <- "Mvnorm" > > ### * Mvnorm > > flush(stderr()); flush(stdout()) > > ### Name: Mvnorm > ### Title: The Multivariate Normal Distribution > ### 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 > x <- rmvnorm(n=100, mean=c(1,1)) > plot(x) > > > > cleanEx(); ..nameEx <- "pmvnorm" > > ### * pmvnorm > > flush(stderr()); flush(stdout()) > > ### 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.580005 attr(,"error") [1] 0.0002696830 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.658117e-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))) > > > > > cleanEx(); ..nameEx <- "pmvt" > > ### * pmvt > > flush(stderr()); flush(stdout()) > > ### Name: pmvt > ### Title: Multivariate t Distribution > ### Aliases: pmvt rmvt > ### 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.0002426591 attr(,"msg") [1] "Normal Completion" > > pmvt(lower=-Inf, upper=3, df = 3, sigma = 1) == pt(3, 3) [1] 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) > ### covariance 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 + } > > round(uniroot(myfct, lower=1, upper=5, alpha=0.95)$root, 3) [1] 2.56 > > # 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)) > a [1] 0.4215702 attr(,"error") [1] 0 attr(,"msg") [1] "Normal Completion" > b [1] 0.4211412 attr(,"error") [1] 3.724896e-06 attr(,"msg") [1] "Normal Completion" > > stopifnot(round(a, 2) == round(b, 2)) > > # correlation and covariance 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.5653977 > b [1] 0.5654018 > 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))) > > rmvt(10, sigma=diag(10)) [,1] [,2] [,3] [,4] [,5] [,6] [1,] 0.58045435 6.74720036 2.96689112 3.0938461 1.49953991 0.88888454 [2,] 13.80520566 -4.38898587 -0.28745384 -14.5762276 17.87266953 32.80428259 [3,] -0.45593386 -0.39997886 0.43090188 -1.3412495 0.09600158 2.31202612 [4,] -3.41158745 -2.66099917 -2.70988090 3.3373662 -0.22885390 0.84239967 [5,] 0.05148137 0.04862904 -0.31780111 0.7405152 -1.17442835 0.83023263 [6,] 1.32104845 0.57212821 -0.51193710 0.3850816 0.41036005 -0.01496229 [7,] -6.57501170 10.16639289 -4.17662982 7.6464247 9.63611592 2.23874176 [8,] -2.76137824 3.68638352 -2.01578371 -2.7415699 5.85968389 3.32451220 [9,] -1.71849590 -1.11906542 -0.53352337 2.6352867 0.21839754 4.68139277 [10,] -1.50759723 -2.07738093 -0.02340225 5.0639045 -0.44541817 -0.22297210 [,7] [,8] [,9] [,10] [1,] -3.2633175 -4.50517238 -0.8189665 5.4916770 [2,] 12.7620653 -30.73869223 33.1063071 -23.6953666 [3,] -0.6351312 -0.86446071 -2.0952141 -2.7979406 [4,] 1.3118519 1.21929555 2.2441642 -1.7635202 [5,] 0.2542591 -0.39286512 -0.2612896 -0.4262057 [6,] 0.7273227 0.74151304 0.4060129 -0.4564147 [7,] -4.3170501 1.09680691 -0.2627469 -1.0686472 [8,] -3.3804101 0.54924603 -3.7736340 -1.6692842 [9,] 2.0162871 -0.01996221 0.4462061 -0.5283162 [10,] 0.9301858 -1.22867136 -2.6359182 1.3542486 > > > > > cleanEx(); ..nameEx <- "qmvnorm" > > ### * qmvnorm > > flush(stderr()); flush(stdout()) > > ### Name: qmvnorm > ### Title: Quantiles of the Multvariate Normal Distribution > ### Aliases: qmvnorm > ### Keywords: distribution > > ### ** Examples > > qmvnorm(0.95, sigma = diag(2), tail = "both") $quantile [1] 2.236497 $f.quantile [1] 2.5831e-06 attr(,"error") [1] 1e-15 attr(,"msg") [1] "Normal Completion" $iter [1] 11 $estim.prec [1] 6.103516e-05 > > > > cleanEx(); ..nameEx <- "qmvt" > > ### * qmvt > > flush(stderr()); flush(stdout()) > > ### Name: qmvt > ### Title: Quantiles of the Multivariate t Distribution > ### Aliases: qmvt > ### Keywords: distribution > > ### ** Examples > > qmvt(0.95, df = 16, tail = "both") $quantile [1] 2.119906 $f.quantile [1] 2.467696e-08 attr(,"error") [1] 0 attr(,"msg") [1] "univariate: using pt" $iter [1] 11 $estim.prec [1] 6.103516e-05 > > > > ### *