R Under development (unstable) (2024-08-15 r87022 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2024 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.

> ## Copyright (C) 2016 Marius Hofert, Ivan Kojadinovic, Martin Maechler, and Jun Yan
> ##
> ## This program 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; either version 3 of the License, or (at your option) any later
> ## version.
> ##
> ## This program 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
> ## this program; if not, see <http://www.gnu.org/licenses/>.
> 
> ## Tests for Explicit Copulas in combination with
> ## rotCopula, mixCopula, and khoudrajiCopula
> 
> library(copula)
> source(system.file("Rsource", "utils.R",     package="copula", mustWork=TRUE))
Loading required package: tools
> ##-> assertError(), assert.EQ(), ... showProc.time()
> isExplicit <- copula:::isExplicit
> (doExtras <- copula:::doExtras())
[1] FALSE
> options(warn = 1)# print warnings as they happen
> 
> 
> ## FIXME? add "revealing" column names to the resulting matrices
> exprDerivs <- function(copula, u) {
+     cbind(copula:::dCdu       (copula, u),
+ 	  copula:::dCdtheta   (copula, u),
+ 	  copula:::dlogcdu    (copula, u),
+ 	  copula:::dlogcdtheta(copula, u))
+ }
> 
> numeDerivs <- function(copula, u) {
+     cbind(copula:::dCduNumer       (copula, u, may.warn = FALSE),
+ 	  copula:::dCdthetaNumer   (copula, u, may.warn = FALSE),
+ 	  copula:::dlogcduNumer    (copula, u),
+ 	  copula:::dlogcdthetaNumer(copula, u))
+ }
> 
> set.seed(123)
> showProc.time()
Time (user system elapsed): 0 0 0 
> 
> mC <- mixCopula(list(gumbelCopula(2.5, dim = 2),
+                      claytonCopula(pi, dim = 2),
+                      indepCopula(dim = 2)),
+                 fixParam(c(2,2,4)/8, c(TRUE, TRUE, TRUE)))
> mC
mixCopula from 3 components
(Gmbc, Clyc, Indc)  with weights:
(0.25, 0.25, 0.50) 
Dimension:  2 
Parameters (partly fixed, see ':='):
  m1.alpha   = 2.500000
  m2.alpha   = 3.141593
  w1        := 0.250000
  w2        := 0.250000
  w3        := 0.500000
> 
> u <- rCopula(100, mC)
> u1 <- (0:16)/16 # includes {0,1}  --> boundary corners & edges
> u12 <- as.matrix(expand.grid(u1=u1, u2=u1, KEEP.OUT.ATTRS=FALSE))
> 
> ## verify density and df in comparison with expression based evaluation
> ## dExplicit..algr() is NaN at boundary (u1 or u2 == 0 or both == 1), and is continuous at other boundaries
> ## dC...() is  := 0  there:
> dC12 <- dCopula(u12, mC)
> dE12 <- copula:::dExplicitCopula.algr(u12, mC)
> i.n <- is.na(dE12)
> i.1 <- u12 == 1
> stopifnot(identical(i.n,
+                     u12[,1] == 0 | u12[,2] == 0 | (i.1[,1] & i.1[,2])))
> ii <- !i.n  &  !i.1[,1]  &  !i.1[,2]
> stopifnot(all.equal(dC12[ii], dE12[ii]))
> stopifnot(all.equal(pCopula(u12, mC), copula:::pExplicitCopula.algr(u12, mC)))
> showProc.time()
Time (user system elapsed): 0.06 0 0.06 
> 
> ## derivatives
> derExp <- exprDerivs(mC, u)
> showProc.time()
Time (user system elapsed): 0 0 0 
> 
> if(doExtras) {
+     derNum <- numeDerivs(mC, u)
+     ## FIXME: dlogcdu and dlogcdtheta do not match the numerical derivatives!
+     ## Could the numerical derivates be wrong??? Why???
+     print(cbind(sapply(1:ncol(derExp), function(i) all.equal(derExp[,i], derNum[,i]))),
+           quote=FALSE)
+     ## very bad from [5] on ..
+     showProc.time()
+ }
> 
> ## benchmark
> ## library(microbenchmark)
> ## microbenchmark(dCopula(u, mC), copula:::dExplicitCopula.algr(u, mC))
> ## microbenchmark(pCopula(u, mC), copula:::pExplicitCopula.algr(u, mC))
> 
> dCd <- copula:::dCdu(mC, u12)
> if(dev.interactive(orNone=TRUE)) {
+     image(u1,u1, matrix(dCd[,1], 16+1))
+     image(u1,u1, matrix(dCd[,2], 16+1))
+ }
> head(cbind(copula:::dCdu(mC, u),    copula:::dCdtheta(mC, u),
+            copula:::dlogcdu(mC, u), copula:::dlogcdtheta(mC, u)))
          [,1]       [,2]         [,3]         [,4]       [,5]       [,6]
[1,] 0.4239764 0.33630909 9.224860e-03 0.0040431376  0.2726189 -1.5074644
[2,] 0.6659714 0.51365201 8.369030e-03 0.0046841141  1.2944452 -0.7795907
[3,] 0.2226718 0.71896413 5.801589e-03 0.0023172852 -2.3256517  2.6063056
[4,] 0.4961576 0.37012542 9.797473e-03 0.0046189657  0.6297653 -1.2699945
[5,] 0.9587609 0.09279068 5.896722e-05 0.0000358527  0.3082482 -0.2891975
[6,] 0.1468872 0.92697224 4.158990e-04 0.0002900826 -0.8519112  0.9271734
            [,7]        [,8]
[1,]  0.09351527  0.10152763
[2,]  0.11740830  0.03889536
[3,] -0.04601827 -0.04515029
[4,]  0.09611079  0.07686739
[5,] -0.02990324 -0.01902671
[6,] -0.07855317 -0.05322670
> 
> ## rotate it
> mC.surv <- rotCopula(mC)
> isExplicit(mC.surv)
[1] TRUE
> 
> stopifnot(all.equal(dCopula(u, mC.surv), dCopula(1 - u, mC)))
> ## rotate it back
> stopifnot(all.equal(dCopula(u, rotCopula(mC.surv)), dCopula(u, mC)))
> showProc.time()
Time (user system elapsed): 0.05 0 0.05 
> 
> ## derivatives
> derExpS <- exprDerivs(mC.surv, u)
Warning in copula:::dCdtheta(copula, u) :
  dCdtheta() returned NaN in column(s) 1 for this explicit copula; falling back to numeric derivative for those columns
> ## (The only numerical derivative case we compute when  doExtras is FALSE : )
> derNumS <- numeDerivs(mC.surv, u)
Warning in copula:::dlogcduNumer(copula, u) :
  Function dlogcdu() not implemented for "rotExplicitCopula" copulas; numerical differentiation used
Warning in copula:::dlogcdthetaNumer(copula, u) :
  Function dlogcdtheta() not implemented for "rotExplicitCopula" copulas; numerical differentiation used
> ## tol=0: show the "closeness" achieved:
> sapply(1:ncol(derExpS), function(i) all.equal(derExpS[,i], derNumS[,i], tol=0)) #
[1] "Mean relative difference: 0.0003624512"
[2] "Mean relative difference: 0.000145243" 
[3] "TRUE"                                  
[4] "Mean relative difference: 1.429454e-12"
[5] "Mean relative difference: 5.991181e-10"
[6] "Mean relative difference: 4.844984e-10"
[7] "Mean relative difference: 8.403848e-10"
[8] "Mean relative difference: 2.206548e-09"
> stopifnot(sapply(1:ncol(derExpS), function(i) all.equal(derExpS[,i], derNumS[,i],
+                                                         tol = 1e-3)))
> showProc.time()
Time (user system elapsed): 1.52 0 1.51 
> 
> 
> ## nest the survival copula in a khoudraji Copula
> k.mC.g <- khoudrajiCopula(mC.surv, gumbelCopula(3, dim = 2), c(.2, .9))
> isExplicit(k.mC.g)
[1] TRUE
> k.mC.g
Khoudraji copula, dim. d = 2, constructed from
 Rotated copula constructed from
  mixCopula from 3 components
 Gumbel copula 
Dimension:  2 
Parameters (partly fixed, see ':='):
  c1.m1.alpha   = 2.500000
  c1.m2.alpha   = 3.141593
  c1.w1        := 0.250000
  c1.w2        := 0.250000
  c1.w3        := 0.500000
  c2.alpha      = 3.000000
  shape1        = 0.200000
  shape2        = 0.900000
> U <- rCopula(100000, k.mC.g)
> 
> ## check cdf with C.n
> 
> ## Must work for d=1 [e.g. in 'gofCopula::gofKendallCvM()]:
> u1 <- as.matrix((0:7)/7)
> Cuu <- C.n(u1,u1) ## <- did not work for 1-dim. u[] matrix for a while
> stopifnot(all.equal(Cuu, c(0:3,5:8)/8, tol = 1e-15))
> ## NOTE(MM): From the above result I think the finite-sample definition of C.n()
> ## ----      "sucks": It is not equidistant when u[] is.
> 
> 
> ## cdf versus C.n
> stopifnot(max(abs(pCopula(u, k.mC.g) - C.n(u, U))) < 0.002)
> 
> ## pdf versus kde2d
> require(MASS)
Loading required package: MASS
> kde <- kde2d(U[,1], U[,2], n = 9, lims = c(0.1, 0.9, 0.1, 0.9))
> max(abs(dCopula(u, k.mC.g) / c(kde$z) - 1)) ## relative difference < 0.185
Warning in dCopula(u, k.mC.g)/c(kde$z) :
  longer object length is not a multiple of shorter object length
[1] 6.010969
> showProc.time()
Time (user system elapsed): 0.81 0.06 0.88 
> 
> ## derivatives
> derE.k <- exprDerivs(k.mC.g, u)
Warning in copula:::dCdtheta(copula, u) :
  dCdtheta() returned NaN in column(s) 1 for this explicit copula; falling back to numeric derivative for those columns
Warning in copula:::dlogcdtheta(copula, u) :
  dlogcdtheta() returned NaN in column(s) 1 for this explicit copula; falling back to numeric derivative for those columns
> showProc.time()
Time (user system elapsed): 0.5 0 0.5 
> 
> if(doExtras) {
+     derN.k <- numeDerivs(k.mC.g, u)
+     print(cbind(sapply(1:ncol(derE.k),
+                        function(i) all.equal(derE.k[,i], derN.k[,i], tol = 0))),
+           quote=FALSE)
+     ## portable ? -- why is [7] so bad ?
+     stopifnot(sapply(1:2,
+                      function(i) all.equal(derE.k[,i], derN.k[,i], tol = 2e-3)),
+               local({ i <- 7 ;   all.equal(derE.k[,i], derN.k[,i], tol = 4e-3) }),
+               sapply(c(3:6, 8:ncol(derE.k)),
+                      function(i) all.equal(derE.k[,i], derN.k[,i], tol = 1e-6))
+               )
+     showProc.time()
+ } # only if ..xtra
> 
> ## nest k.mC.g and mC in a mixture copula
> m.k.m <- mixCopula(list(mC, k.mC.g), c(.5, .5))
> m.k.m
mixCopula from 2 components
(mf3c          , KccfRccfmf3cGc)  with weights:
(0.5, 0.5) 
Dimension:  2 
Parameters (partly fixed, see ':='):
  m1.m1.alpha      = 2.500000
  m1.m2.alpha      = 3.141593
  m1.w1           := 0.250000
  m1.w2           := 0.250000
  m1.w3           := 0.500000
  m2.c1.m1.alpha   = 2.500000
  m2.c1.m2.alpha   = 3.141593
  m2.c1.w1        := 0.250000
  m2.c1.w2        := 0.250000
  m2.c1.w3        := 0.500000
  m2.c2.alpha      = 3.000000
  m2.shape1        = 0.200000
  m2.shape2        = 0.900000
  w1               = 0.500000
  w2               = 0.500000
> U <- rCopula(10000, m.k.m)
> stopifnot(max(abs(pCopula(u, m.k.m) - C.n(u, U))) < 0.02)
> 
> ## a monster from m.k.m and mC.surv in khoudraji
> monster <- khoudrajiCopula(m.k.m, mC.surv, c(.2, .8))
> monster
Khoudraji copula, dim. d = 2, constructed from
 mixCopula from 2 components
 Rotated copula constructed from
  mixCopula from 3 components 
Dimension:  2 
Parameters (partly fixed, see ':='):
  c1.m1.m1.alpha      = 2.500000
  c1.m1.m2.alpha      = 3.141593
  c1.m1.w1           := 0.250000
  c1.m1.w2           := 0.250000
  c1.m1.w3           := 0.500000
  c1.m2.c1.m1.alpha   = 2.500000
  c1.m2.c1.m2.alpha   = 3.141593
  c1.m2.c1.w1        := 0.250000
  c1.m2.c1.w2        := 0.250000
  c1.m2.c1.w3        := 0.500000
  c1.m2.c2.alpha      = 3.000000
  c1.m2.shape1        = 0.200000
  c1.m2.shape2        = 0.900000
  c1.w1               = 0.500000
  c1.w2               = 0.500000
  c2.m1.alpha         = 2.500000
  c2.m2.alpha         = 3.141593
  c2.w1              := 0.250000
  c2.w2              := 0.250000
  c2.w3              := 0.500000
  shape1              = 0.200000
  shape2              = 0.800000
> U <- rCopula(10000, monster)
> stopifnot(max(abs(pCopula(u, monster) - C.n(u, U))) < 0.007)
> showProc.time()
Time (user system elapsed): 0.12 0 0.12 
> 
> ## derivatives
> derE.M <- exprDerivs(monster, u)
Warning in copula:::dCdtheta(copula, u) :
  dCdtheta() returned NaN in column(s) 3, 10 for this explicit copula; falling back to numeric derivative for those columns
Warning in copula:::dlogcdtheta(copula, u) :
  dlogcdtheta() returned NaN in column(s) 3, 10 for this explicit copula; falling back to numeric derivative for those columns
> showProc.time()
Time (user system elapsed): 2.13 0 2.13 
> if(doExtras) {
+     derN.M <- numeDerivs(monster, u)
+     print(cbind(sapply(1:ncol(derE.M), function(i) all.equal(derE.M[,i], derN.M[,i], tol=0))),
+           quote = FALSE)
+     showProc.time()
+ } # only if ..xtra
> 
> 
> ## rotate the monster
> rM <- rotCopula(monster, flip=c(TRUE, FALSE))
> isExplicit(rM)
[1] TRUE
> rM
Rotated copula constructed from
Khoudraji copula, dim. d = 2, constructed from
 mixCopula from 2 components
 Rotated copula constructed from
  mixCopula from 3 components 
Dimension:  2 
Parameters (partly fixed, see ':='):
  c1.m1.m1.alpha      = 2.500000
  c1.m1.m2.alpha      = 3.141593
  c1.m1.w1           := 0.250000
  c1.m1.w2           := 0.250000
  c1.m1.w3           := 0.500000
  c1.m2.c1.m1.alpha   = 2.500000
  c1.m2.c1.m2.alpha   = 3.141593
  c1.m2.c1.w1        := 0.250000
  c1.m2.c1.w2        := 0.250000
  c1.m2.c1.w3        := 0.500000
  c1.m2.c2.alpha      = 3.000000
  c1.m2.shape1        = 0.200000
  c1.m2.shape2        = 0.900000
  c1.w1               = 0.500000
  c1.w2               = 0.500000
  c2.m1.alpha         = 2.500000
  c2.m2.alpha         = 3.141593
  c2.w1              := 0.250000
  c2.w2              := 0.250000
  c2.w3              := 0.500000
  shape1              = 0.200000
  shape2              = 0.800000
> 
> U <- rCopula(10000, rM)
> max(abs(pCopula(u, rM) - C.n(u, U))) # < 0.005
[1] 0.005771696
> stopifnot(identical(dCopula(u, rM), dCopula(cbind(1 - u[,1], u[,2]), monster)))
> ## derivatives
> derE.rM <- exprDerivs(rM, u)
Warning in copula:::dCdtheta(copula, u) :
  dCdtheta() returned NaN in column(s) 1, 2, 3, 4, 5, 6, 7, 10, 11, 12, 13 for this explicit copula; falling back to numeric derivative for those columns
Warning in copula:::dlogcdtheta(copula, u) :
  dlogcdtheta() returned NaN in column(s) 3, 10 for this explicit copula; falling back to numeric derivative for those columns
> showProc.time()
Time (user system elapsed): 2.81 0.02 2.82 
> if(doExtras) {
+     derN.rM <- numeDerivs(rM, u)
+     ## show diff's:
+     print(cbind(sapply(1:ncol(derE.rM), function(i) all.equal(derE.rM[,i], derN.rM[,i], tol=0))),
+           quote = FALSE)
+     showProc.time()
+ } # only if ..xtra
> 
> ##########################################################
> ## joeCopula
> 
> jC <- joeCopula(4, dim = 2)
> 
> ## "derNum <- cbind"
> ## pdf/cdf
> stopifnot(all.equal(dCopula(u, jC), copula:::dExplicitCopula.algr(u, jC)))
> stopifnot(all.equal(pCopula(u, jC), copula:::pExplicitCopula.algr(u, jC)))
> showProc.time()
Time (user system elapsed): 0.02 0 0.02 
> 
> ## derivatives  -- are fast here
> derE.j <- exprDerivs(jC, u)
> derN.j <- numeDerivs(jC, u)
Warning in copula:::dlogcduNumer(copula, u) :
  Function dlogcdu() not implemented for "joeCopula" copulas; numerical differentiation used
Warning in copula:::dlogcdthetaNumer(copula, u) :
  Function dlogcdtheta() not implemented for "joeCopula" copulas; numerical differentiation used
> sapply(1:ncol(derE.j), function(i) all.equal(derE.j[,i], derN.j[,i], tol=0))
[1] "Mean relative difference: 0.0003640905"
[2] "Mean relative difference: 9.717875e-05"
[3] "Mean relative difference: 3.033551e-13"
[4] "Mean relative difference: 3.59782e-10" 
[5] "Mean relative difference: 4.718684e-10"
[6] "Mean relative difference: 6.362816e-10"
> showProc.time()
Time (user system elapsed): 0.57 0.01 0.59 
> 
> ## rotate it
> rJ <- rotCopula(jC, flip=c(TRUE, FALSE))
> derE.rJ <- exprDerivs(rJ, u)
Warning in copula:::dCdtheta(copula, u) :
  dCdtheta() returned NaN in column(s) 1 for this explicit copula; falling back to numeric derivative for those columns
> derN.rJ <- numeDerivs(rJ, u)
> sapply(1:ncol(derE.rJ), function(i) all.equal(derE.rJ[,i], derN.rJ[,i], tol=0))
[1] "Mean relative difference: 0.0004290902"
[2] "Mean relative difference: 0.000516896" 
[3] "TRUE"                                  
[4] "Mean relative difference: 4.302875e-09"
[5] "Mean relative difference: 1.161346e-09"
[6] "Mean relative difference: 3.264156e-09"
> showProc.time()
Time (user system elapsed): 1.07 0.02 1.08 
> 
> ## mix it with k.mc.g
> hiro <- mixCopula(list(jC, k.mC.g), c(.5, .5))
> hiro
mixCopula from 2 components
(Jcpl          , KccfRccfmf3cGc)  with weights:
(0.5, 0.5) 
Dimension:  2 
Parameters (partly fixed, see ':='):
  m1.alpha         = 4.000000
  m2.c1.m1.alpha   = 2.500000
  m2.c1.m2.alpha   = 3.141593
  m2.c1.w1        := 0.250000
  m2.c1.w2        := 0.250000
  m2.c1.w3        := 0.500000
  m2.c2.alpha      = 3.000000
  m2.shape1        = 0.200000
  m2.shape2        = 0.900000
  w1               = 0.500000
  w2               = 0.500000
> 
> ##########################################################
> 
> proc.time()
   user  system elapsed 
  10.93    0.23   11.17