R Under development (unstable) (2026-02-16 r89426 ucrt) -- "Unsuffered Consequences" Copyright (C) 2026 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 . > > ### --- Khoudraji Copula ----------------------- > > library(copula) > > source(system.file("Rsource", "utils.R", package="copula", mustWork=TRUE)) Loading required package: tools > ## tryCatch.W.E() etc > > (doExtras <- copula:::doExtras()) [1] FALSE > > isExplicit <- copula:::isExplicit > ## if(!dev.interactive(orNone=TRUE)) pdf("khoudraji-ex.pdf") > > ################################################################################ > ## An explicit Khoudraji copula constructed from claytonCopula and indepCoula > ## with one fixed shape parameter > kcf <- khoudrajiCopula(copula2 = claytonCopula(6), + shapes = fixParam(c(0.4, 0.95), c(FALSE, TRUE))) > kcf Khoudraji copula, dim. d = 2, constructed from Independence copula Clayton copula Dimension: 2 Parameters (partly fixed, see ':='): c2.alpha = 6.00 shape1 = 0.40 shape2 := 0.95 > stopifnot( + all.equal(getTheta(kcf, freeOnly = FALSE, named = TRUE) -> th., + c(c2.alpha = 6, shape1 = 0.4, shape2 = 0.95)), + identical(getTheta(kcf, named = TRUE), th.[1:2])) > > ## kcf@exprdist$cdf > stopifnot(isExplicit(kcf)) > > set.seed(123) > u <- rCopula(20, kcf) > > ## cdf/pdf from explicit expression versus from algorithmic implementation > summary(rD <- pCopula(u, kcf) - copula:::pKhoudrajiCopula(u, kcf)) Min. 1st Qu. Median Mean 3rd Qu. Max. 0 0 0 0 0 0 > summary(rD. <- dCopula(u, kcf) - copula:::dKhoudrajiBivCopula(u, kcf)) Min. 1st Qu. Median Mean 3rd Qu. Max. -3.375e-14 -8.327e-16 -2.776e-16 -2.293e-15 0.000e+00 2.220e-16 > stopifnot(max(abs(rD)) < 1e-12, max(abs(rD.)) < 1e-12) > showProc.time() Time (user system elapsed): 0.03 0 0.03 > > ################################################################################ > ## A nested Khoudraji copula from kcf and a gumbelCopula, still explicit > k_kcf_g <- khoudrajiCopula(kcf, gumbelCopula(2), + shapes = fixParam(c(.2, .9), c(TRUE, TRUE))) > k_kcf_g Khoudraji copula, dim. d = 2, constructed from Khoudraji copula constructed from Independence copula Clayton copula Gumbel copula Dimension: 2 Parameters (partly fixed, see ':='): c1.c2.alpha = 6.00 c1.shape1 = 0.40 c1.shape2 := 0.95 c2.alpha = 2.00 shape1 := 0.20 shape2 := 0.90 > stopifnot(isExplicit(k_kcf_g)) > getTheta(k_kcf_g, freeOnly = FALSE, named = TRUE) c1.c2.alpha c1.shape1 c1.shape2 c2.alpha shape1 shape2 6.00 0.40 0.95 2.00 0.20 0.90 > getTheta(k_kcf_g, named = TRUE) c1.c2.alpha c1.shape1 c2.alpha 6.0 0.4 2.0 > > ## k_kcf_g@exprdist$cdf > set.seed(456) > U <- rCopula(100000, k_kcf_g) > > ## check cdf/pdf with C.n and kde2d > require(MASS) Loading required package: MASS > u <- as.matrix(expand.grid((1:9)/10, (1:9)/10)) > ## cdf versus C.n > summary(dp <- pCopula(u, k_kcf_g) - C.n(u, U)) ## |.| <= 9.8e-4) Min. 1st Qu. Median Mean 3rd Qu. Max. -8.705e-04 -1.950e-04 4.744e-05 9.109e-05 3.949e-04 9.825e-04 > stopifnot(max(abs(dp)) < 3e-3) > ## pdf versus kde2d > kde <- kde2d(U[,1], U[,2], n = 9, lims = c(0.1, 0.9, 0.1, 0.9)) > summary(rD <- dCopula(u, k_kcf_g) / c(kde$z) - 1) ## |relative difference| <= 0.123.. Min. 1st Qu. Median Mean 3rd Qu. Max. -0.057814 -0.017136 0.003081 0.003665 0.017643 0.123191 > stopifnot(max(abs(rD)) < 0.13) > showProc.time() Time (user system elapsed): 0.96 0.05 1 > > ################################################################################ > ## one more nesting: kcf and k_kcf_g > monster <- khoudrajiCopula(kcf, k_kcf_g, + shapes = fixParam(c(.9, .1), c(TRUE, FALSE))) > monster; validObject(monster) Khoudraji copula, dim. d = 2, constructed from Khoudraji copula constructed from Independence copula Clayton copula Khoudraji copula constructed from Khoudraji copula constructed from Independence copula Clayton copula Gumbel copula Dimension: 2 Parameters (partly fixed, see ':='): c1.c2.alpha = 6.00 c1.shape1 = 0.40 c1.shape2 := 0.95 c2.c1.c2.alpha = 6.00 c2.c1.shape1 = 0.40 c2.c1.shape2 := 0.95 c2.c2.alpha = 2.00 c2.shape1 := 0.20 c2.shape2 := 0.90 shape1 := 0.90 shape2 = 0.10 [1] TRUE > stopifnot(isExplicit(monster)) > > (th.mA <- getTheta(monster, freeOnly = FALSE, named = TRUE)) c1.c2.alpha c1.shape1 c1.shape2 c2.c1.c2.alpha c2.c1.shape1 6.00 0.40 0.95 6.00 0.40 c2.c1.shape2 c2.c2.alpha c2.shape1 c2.shape2 shape1 0.95 2.00 0.20 0.90 0.90 shape2 0.10 > (th.m <- getTheta(monster, named = TRUE)) c1.c2.alpha c1.shape1 c2.c1.c2.alpha c2.c1.shape1 c2.c2.alpha 6.0 0.4 6.0 0.4 2.0 shape2 0.1 > stopifnot(identical(th.m, th.mA[isFree(monster)])) > > set.seed(488) > U <- rCopula(100000, monster) > ## cdf versus C.n > summary(dp <- pCopula(u, monster) - C.n(u, U)) ## -0.0083 ... 0.00131.. Min. 1st Qu. Median Mean 3rd Qu. Max. -0.0008330 -0.0001364 0.0001537 0.0001981 0.0004075 0.0013149 > ## pdf versus kde2d > kde <- kde2d(U[,1], U[,2], n = 9, lims = c(0.1, 0.9, 0.1, 0.9)) > summary(rdC <- dCopula(u, monster) / c(kde$z) - 1) ## relative difference < 0.09 Min. 1st Qu. Median Mean 3rd Qu. Max. -0.082164 -0.013860 0.001654 0.003298 0.026558 0.057520 > stopifnot(abs(dp) < 0.01, abs(rdC) < 0.09) > showProc.time() Time (user system elapsed): 1.28 0.09 1.38 > > ################################################################################ > ## khoudrajiExplicitCopula with dim 3 > kcd3 <- khoudrajiCopula(copula1 = indepCopula(dim=3), + copula2 = claytonCopula(6, dim=3), + shapes = c(0.4, 0.95, 0.95)) > kcd3 Khoudraji copula, dim. d = 3, constructed from Independence copula Clayton copula Dimension: 3 Parameters: c2.alpha = 6.00 shape1 = 0.40 shape2 = 0.95 shape3 = 0.95 > > set.seed(1712) > U <- rCopula(100000, kcd3) > u <- matrix(runif(15), 5, 3) > ## cdf versus C.n > summary(dP <- pCopula(u, kcd3) - C.n(u, U)) # < 0.00007 Min. 1st Qu. Median Mean 3rd Qu. Max. -6.882e-04 -6.014e-04 -5.745e-04 -4.670e-04 -4.015e-04 -6.929e-05 > ## don't know how to check pdf > (f.v <- dCopula(u, kcd3)) [1] 1.1116773 1.6661734 0.1719080 0.4981574 9.8964415 > stopifnot(expr = { + isExplicit(kcd3) + inherits(kcd3, "khoudrajiExplicitCopula") + all.equal(f.v, tolerance = 1e-7, + c(1.1116773, 1.6661734, 0.1719080, 0.4981574, 9.8964415)) + abs(dP) < 8e-4 + }) > showProc.time() Time (user system elapsed): 0.53 0 0.53 > > kcg <- khoudrajiCopula(claytonCopula(2, dim=3), gumbelCopula(3, dim=3), c(0.2, 0.2, .8)) > kcg Khoudraji copula, dim. d = 3, constructed from Clayton copula Gumbel copula Dimension: 3 Parameters: c1.alpha = 2.0 c2.alpha = 3.0 shape1 = 0.2 shape2 = 0.2 shape3 = 0.8 > u <- rCopula(10, kcg) > dCopula(u, kcg) [1] 2.780102 2.723464 1.041935 2.027283 1.675718 2.307609 1.203909 [8] 10.459996 0.386535 3.480239 > > > kkcgg <- khoudrajiCopula(kcg, gumbelCopula(2, dim = 3), c(.3, .3, .9)) > kkcgg Khoudraji copula, dim. d = 3, constructed from Khoudraji copula constructed from Clayton copula Gumbel copula Gumbel copula Dimension: 3 Parameters: c1.c1.alpha = 2.0 c1.c2.alpha = 3.0 c1.shape1 = 0.2 c1.shape2 = 0.2 c1.shape3 = 0.8 c2.alpha = 2.0 shape1 = 0.3 shape2 = 0.3 shape3 = 0.9 > u <- rCopula(10, kkcgg) > dCopula(u, kkcgg) [1] 0.3665809 1.5981348 1.0638901 1.9684222 1.5384452 1.5214729 1.3522174 [8] 0.9529139 1.2088327 0.4565743 > > > kcgcd3 <- khoudrajiCopula(kcd3, kkcgg, c(.4, .4, .5)) > kcgcd3 Khoudraji copula, dim. d = 3, constructed from Khoudraji copula constructed from Independence copula Clayton copula Khoudraji copula constructed from Khoudraji copula constructed from Clayton copula Gumbel copula Gumbel copula Dimension: 3 Parameters: c1.c2.alpha = 6.00 c1.shape1 = 0.40 c1.shape2 = 0.95 c1.shape3 = 0.95 c2.c1.c1.alpha = 2.00 c2.c1.c2.alpha = 3.00 c2.c1.shape1 = 0.20 c2.c1.shape2 = 0.20 c2.c1.shape3 = 0.80 c2.c2.alpha = 2.00 c2.shape1 = 0.30 c2.shape2 = 0.30 c2.shape3 = 0.90 shape1 = 0.40 shape2 = 0.40 shape3 = 0.50 > u <- rCopula(10, kcgcd3) > dCopula(u, kcgcd3) [1] 0.5406296 4.7736360 0.5877674 0.7449026 2.2931403 0.3452497 1.5496446 [8] 0.3569642 2.9840308 1.2362195 > > showProc.time() Time (user system elapsed): 0.11 0 0.11 > if(doExtras) { + print(compD <- comparederiv(kcgcd3, u), digits = 4) + ## dCdu dCdtheta dlogcdu dlogcdtheta + ## 4.387280e-03 9.791440e-04 3.009563e-08 7.572446e-09 + stopifnot(compD < c(0.01, 3e-3, 1e-7, 8e-8)) + showProc.time() + } > ############################################################################# > ## fitting checking > > do1 <- function(n, cop) { + U <- pobs(rCopula(n, cop)) + fit <- try(fitCopula(cop, U, method = "mpl")) + p <- nParam(cop, freeOnly=TRUE) + if (inherits(fit, "try-error")) rep(NA_real_, 2 * p) + else c(coef(fit), sqrt(diag(vcov(fit)))) + } > > sumsim <- function(sim) { + p <- nrow(sim) / 2 + mm <- apply(sim, 1, mean, na.rm = TRUE) + ss <- apply(sim, 1, sd, na.rm = TRUE) + data.frame(estimate = mm[1:p], ese = ss[1:p], ase = mm[p + 1:p]) + } > > if (doExtras) { + set.seed(123) + nrep <- 50 + n <- 1000 + + sim <- replicate(nrep, do1(n, kcg)) + sumsim(sim) + cbind(true = getTheta(kcg), sumsim(sim)) + ## the average se is greater than empirical for the shape parameters; need more investigation + } > > proc.time() user system elapsed 4.35 0.40 4.76