R Under development (unstable) (2023-11-26 r85638 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2023 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

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 an HTML browser interface to help.
Type 'q()' to quit R.

> library("mvtnorm")
> 
> chk <- function(...) isTRUE(all.equal(...))
> 
> ## Showing the TVPACK() gives *NON*-random results:
> (cor1 <- toeplitz(c(1, 1/4, -1/8)))
       [,1] [,2]   [,3]
[1,]  1.000 0.25 -0.125
[2,]  0.250 1.00  0.250
[3,] -0.125 0.25  1.000
> (up1  <- c(1/4, 7/4, 5/8))
[1] 0.250 1.750 0.625
> d <- length(up1) # = 3
> pmvt.. <- function(df, algorithm)
+     vapply(df, function(df) pmvt(upper=up1, corr=cor1, df=df, algorithm=algorithm),
+            numeric(1))
> 
> dfs <- 1:9
> pmvt_TV.7 <- replicate(7, pmvt..(dfs, TVPACK()))
> 
> stopifnot(identical(unique(c(pmvt_TV.7 - pmvt_TV.7[,1])), 0))
> (pmvt.TV. <- pmvt_TV.7[,1])
[1] 0.3554119 0.3817313 0.3923546 0.3980570 0.4016042 0.4040204 0.4057708
[8] 0.4070968 0.4081358
> (pmvt.TV  <- pmvt..(dfs, TVPACK(1e-14)))# has no effect here
[1] 0.3554119 0.3817313 0.3923546 0.3980570 0.4016042 0.4040204 0.4057708
[8] 0.4070968 0.4081358
> chk(max(abs(pmvt.TV - pmvt.TV.)), 0) ## all 0 {unexpectedly ??}
[1] TRUE
> 
> 
> set.seed(47) ## and default algorithm: -> *random* result
> pmvt_7 <- replicate(7, vapply(dfs, function(df) pmvt(df=df, upper=up1, corr=cor1), numeric(1)))
> ## relative errors
> relE <- 1 - pmvt_7 / pmvt.TV
> rng.rE <- range(abs(relE))
> stopifnot(1e-6 < rng.rE[1], rng.rE[2] < 7e-4)
> stopifnot(chk(
+     colMeans(abs(relE)),
+     c(88, 64, 105, 73, 52, 90, 87)*1e-6, tol= 1e-3))
> 
> 
> set.seed(29)
> 
> ########################################################################
> ## 3 dim example
> corr <- cov2cor(crossprod(matrix(runif(9,-1,1),3,3))+diag(3))
> df <- rpois(1,3)+1
> 
> ## central t distribution (-Inf,upper)
> ctrl <- GenzBretz(maxpts = 2500000, abseps = 0.000001, releps = 0)
> upper <- rexp(3,1)
> pmvt(upper=upper, corr=corr, df = df, algorithm = ctrl)
[1] 0.3920567
attr(,"error")
[1] 6.089669e-07
attr(,"msg")
[1] "Normal Completion"
> pmvt(upper=upper, corr=corr, df = df, algorithm = TVPACK())
[1] 0.3920566
attr(,"error")
[1] 1e-06
attr(,"msg")
[1] "Normal Completion"
> 
> ## central t distribution (lower,Inf)
> lower <- -rexp(3,1)
> pmvt(lower=lower, upper=rep(Inf,3), corr=corr, df = df, algorithm = ctrl)
[1] 0.4634843
attr(,"error")
[1] 6.46065e-07
attr(,"msg")
[1] "Normal Completion"
> pmvt(lower=lower, upper=rep(Inf,3), corr=corr, df = df, algorithm = TVPACK())
[1] 0.4634844
attr(,"error")
[1] 1e-06
attr(,"msg")
[1] "Normal Completion"
> 
> ## non-central t (not possible for TVPACK)
> delt <- rexp(3,1/10)
> upper <- delt+runif(3)
> ctrl <- GenzBretz(maxpts = 2500000, abseps = 0.000001, releps = 0)
> pmvt(upper=upper, corr=corr, df = df, algorithm = ctrl, delta = delt)
[1] 0.3235424
attr(,"error")
[1] 9.86407e-07
attr(,"msg")
[1] "Normal Completion"
> tools::assertError(pmvt(upper=upper, corr=corr, df = df, algorithm = TVPACK(), delta = delt))
> 
> ## central mvn (-Inf, upper)
> upper <- rexp(3,1)
> pmvnorm(upper=upper, corr=corr, algorithm = ctrl)
[1] 0.7733949
attr(,"error")
[1] 6.999168e-07
attr(,"msg")
[1] "Normal Completion"
> pmvnorm(upper=upper, corr=corr, algorithm = TVPACK())
[1] 0.7733949
attr(,"error")
[1] 1e-06
attr(,"msg")
[1] "Normal Completion"
> 
> ## central mvn (lower, Inf)
> lower <- rexp(3,5)
> pmvnorm(lower=lower,upper=rep(Inf, 3), corr=corr, algorithm = ctrl)
[1] 0.0827153
attr(,"error")
[1] 4.969585e-07
attr(,"msg")
[1] "Normal Completion"
> pmvnorm(lower=lower,upper=rep(Inf, 3), corr=corr, algorithm = TVPACK())
[1] 0.08271559
attr(,"error")
[1] 1e-06
attr(,"msg")
[1] "Normal Completion"
> 
> ## non-central mvn
> delt <- rexp(3,1/10)
> upper <- delt+rexp(3,1)
> pmvnorm(upper=upper, corr=corr, algorithm = ctrl,     mean = delt)
[1] 0.5824373
attr(,"error")
[1] 9.416391e-07
attr(,"msg")
[1] "Normal Completion"
> pmvnorm(upper=upper, corr=corr, algorithm = TVPACK(), mean = delt) # should not error
[1] 0.582434
attr(,"error")
[1] 1e-06
attr(,"msg")
[1] "Normal Completion"
> 
> ########################################################################
> ## 2 dim example
> corr <- cov2cor(crossprod(matrix(runif(4,-1,1),2,2))+diag(2))
> upper <- rexp(2,1)
> df <- rpois(1, runif(1, 0, 20))
> 
> ## central t (-Inf, upper)
> pmvt(upper=upper, corr=corr, df = df, algorithm = ctrl)
[1] 0.8656102
attr(,"error")
[1] 1e-15
attr(,"msg")
[1] "Normal Completion"
> pmvt(upper=upper, corr=corr, df = df, algorithm = TVPACK())
[1] 0.8656102
attr(,"error")
[1] NA
attr(,"msg")
[1] "Normal Completion"
> 
> ## central t (lower, Inf)
> pmvt(lower=-upper, upper=rep(Inf, 2), corr=corr, df = df, algorithm = ctrl)
[1] 0.8656102
attr(,"error")
[1] 1e-15
attr(,"msg")
[1] "Normal Completion"
> pmvt(lower=-upper, upper=rep(Inf, 2), corr=corr, df = df, algorithm = TVPACK())
[1] 0.8656102
attr(,"error")
[1] NA
attr(,"msg")
[1] "Normal Completion"
> 
> ## non-central t
> delt <- rexp(2,1/5)
> upper <- delt+rexp(2,1)
> pmvnorm(upper=upper, corr=corr, algorithm = ctrl, mean = delt)
[1] 0.6431222
attr(,"error")
[1] 1e-15
attr(,"msg")
[1] "Normal Completion"
> pmvnorm(upper=upper, corr=corr, algorithm = TVPACK(), mean = delt)
[1] 0.6431222
attr(,"error")
[1] NA
attr(,"msg")
[1] "Normal Completion"
> 
> ########################################################################
> ## comparison with Miwa
> ## 2d
> corr <- cov2cor(crossprod(matrix(runif(4,-1,1),2,2))+diag(2))
> upper <- rexp(2, 1)
> 
> pmvnorm(upper=upper, corr=corr, algorithm = Miwa(steps=128))
[1] 0.7331788
attr(,"error")
[1] NA
attr(,"msg")
[1] "Normal Completion"
> pmvnorm(upper=upper, corr=corr, algorithm = TVPACK())
[1] 0.7331788
attr(,"error")
[1] NA
attr(,"msg")
[1] "Normal Completion"
> 
> ## 3d
> corr <- cov2cor(crossprod(matrix(runif(9,-1,1),3,3))+diag(3))
> upper <- rexp(3, 1)
> 
> ctrl <- Miwa(steps=128)
> pmvnorm(upper=upper, corr=corr, algorithm = ctrl)
[1] 0.8981829
attr(,"error")
[1] NA
attr(,"msg")
[1] "Normal Completion"
> pmvnorm(upper=upper, corr=corr, algorithm = TVPACK())
[1] 0.8981829
attr(,"error")
[1] 1e-06
attr(,"msg")
[1] "Normal Completion"
> 
> ##==== Cases where some  (lower[j], upper[j]) == (-Inf, Inf) :
> S <- toeplitz(c(1, 1/2, 1/4))
> 
> set.seed(11)
> P0 <- pmvnorm(lower=c(-Inf, 0, 0), upper=Inf, corr=S)
> P1 <- pmvnorm(lower=c(-Inf, 0, 0), upper=Inf, corr=S, algorithm = TVPACK()) # had failed
> P2 <- pmvnorm(lower=c(-Inf, 0, 0), upper=Inf, corr=S, algorithm = Miwa())
> P2a<- pmvnorm(lower=c(-Inf, 0, 0), upper=Inf, corr=S, algorithm = Miwa(512))
> P2.<- pmvnorm(lower=c(-Inf, 0, 0), upper=Inf, corr=S, algorithm = Miwa(2048))
> 
> stopifnot(chk(1/3, c(P0), tol=1e-14),
+           chk(1/3, c(P1), tol=1e-14),
+           chk(1/3, c(P2), tol=1e-9 ), # 3.765e-10
+           chk(1/3, c(P2a),tol=4e-12), # 8.32 e-13
+           chk(1/3, c(P2.),tol=2e-12) # 5.28 e-13
+ )
> 
> ## t-dist [TVPACK() had failed] :
> set.seed(11)
> Ptdef <- replicate(20, c(pmvt(lower=c(-Inf, 1, 2), upper=Inf, df=2, corr=S)))
> unique(Ptdef)# see length 1; i.e., same result [even though default is Monte-Carlo ??]
[1] 0.0570404
> Pt1 <- pmvt(lower=c(-Inf, 1, 2), upper=Inf, df=2, corr=S, algorithm = TVPACK())
> P. <- 0.0570404044526986
> stopifnot(exprs = {
+     chk(P., c(Pt1), tol = 1e-14)# seen 3.65 e-16
+     abs(P. - Ptdef) < 1e-15 # seen 1.39 e-17
+ })
> 
> proc.time()
   user  system elapsed 
   8.62    0.06    8.65