R/gamBiCopPDF.R
gamBiCopPDF.RdThis function returns the density of a bivariate conditional copula, where either the copula parameter or the Kendall's tau is modeled as a function of the covariates.
gamBiCopPDF(object, newdata = NULL)gamBiCop-class object.
(Same as in predict.gam from the
mgcv package) A matrix or data frame
containing the values of the model covariates at which predictions are
required, along with two columns named `"u1"` and `"u2"`.
If this is not provided then the density corresponding to the
original data are returned. If newdata is provided then it should contain all
the variables needed for prediction: a warning is generated if not.
The conditional density.
gamBiCop and gamBiCopPredict.
require(copula)
set.seed(0)
## Simulation parameters (sample size, correlation between covariates,
## Gaussian copula family)
n <- 2e2
rho <- 0.5
fam <- 1
## A calibration surface depending on three variables
eta0 <- 1
calib.surf <- list(
calib.quad <- function(t, Ti = 0, Tf = 1, b = 8) {
Tm <- (Tf - Ti) / 2
a <- -(b / 3) * (Tf^2 - 3 * Tf * Tm + 3 * Tm^2)
return(a + b * (t - Tm)^2)
},
calib.sin <- function(t, Ti = 0, Tf = 1, b = 1, f = 1) {
a <- b * (1 - 2 * Tf * pi / (f * Tf * pi +
cos(2 * f * pi * (Tf - Ti))
- cos(2 * f * pi * Ti)))
return((a + b) / 2 + (b - a) * sin(2 * f * pi * (t - Ti)) / 2)
},
calib.exp <- function(t, Ti = 0, Tf = 1, b = 2, s = Tf / 8) {
Tm <- (Tf - Ti) / 2
a <- (b * s * sqrt(2 * pi) / Tf) * (pnorm(0, Tm, s) - pnorm(Tf, Tm, s))
return(a + b * exp(-(t - Tm)^2 / (2 * s^2)))
}
)
## 3-dimensional matrix X of covariates
covariates.distr <- mvdc(normalCopula(rho, dim = 3),
c("unif"), list(list(min = 0, max = 1)),
marginsIdentical = TRUE
)
X <- rMvdc(n, covariates.distr)
colnames(X) <- paste("x", 1:3, sep = "")
## U in [0,1]x[0,1] with copula parameter depending on X
U <- condBiCopSim(fam, function(x1, x2, x3) {
eta0 + sum(mapply(function(f, x)
f(x), calib.surf, c(x1, x2, x3)))
}, X[, 1:3], par2 = 6, return.par = TRUE)
## Merge U and X
data <- data.frame(U$data, X)
names(data) <- c(paste("u", 1:2, sep = ""), paste("x", 1:3, sep = ""))
## Model fit with penalized cubic splines (via min GCV)
basis <- c(3, 10, 10)
formula <- ~ s(x1, k = basis[1], bs = "cr") +
s(x2, k = basis[2], bs = "cr") +
s(x3, k = basis[3], bs = "cr")
system.time(fit <- gamBiCopFit(data, formula, fam))
#> user system elapsed
#> 0.162 0.001 0.162
## Evaluate the conditional density
gamBiCopPDF(fit$res)
#> [1] 1.53897409 0.65895255 7.54308919 7.29902794 7.98200466
#> [6] 0.91280983 0.65051146 1.08273677 0.13773165 8.67048283
#> [11] 4.95057615 3.20943703 1.00633039 0.84382926 0.40244706
#> [16] 2.76005256 1.96292286 11.68003966 1.47383017 1.07385649
#> [21] 0.99519798 1.18493994 1.44425297 0.83340597 2.40090338
#> [26] 6.05695932 1.41528126 3.00012259 1.74635340 0.89680036
#> [31] 10.49237658 1.42704608 1.01410397 1.20991114 0.31228006
#> [36] 1.94781172 1.08920172 10.33179260 0.96340114 1.10652903
#> [41] 1.19280224 1.53050049 1.15345886 1.49250150 3.65567450
#> [46] 1.15697132 10.71967518 0.69294323 1.64953687 1.04327879
#> [51] 88.38436240 3.62295213 10.47617993 0.77327325 2.40824831
#> [56] 0.08006449 1.37469193 3.00543370 0.85171208 1.12234724
#> [61] 1.60635642 0.83929746 2.41797654 0.81561166 4.07324027
#> [66] 32.31728518 0.77131922 1.03527579 2.91140932 2.50619363
#> [71] 1.14147003 1.34993051 5.42784075 1.65326684 1.27372418
#> [76] 1.50209872 1.03634247 1.00381308 1.44370479 5.84889100
#> [81] 5.19282176 0.76370453 100.00992066 2.79801297 1.26603622
#> [86] 0.99666029 5.95339495 0.77555019 1.29134500 0.85042384
#> [91] 1.68855016 13.22072485 1.30870031 1.09627773 0.98611410
#> [96] 1.64269684 1.25826086 18.73796411 1.26322251 1.20958580
#> [101] 1.14939531 0.56989490 1.78403902 1.20540159 9.17700886
#> [106] 1.86256652 2.08354774 1.27445574 0.81208319 1.09719242
#> [111] 0.70119548 1.42652282 1.74619315 31.69259035 1.57919837
#> [116] 0.07155906 3.75765970 1.06833619 1.06277686 2.75424876
#> [121] 1.03289717 0.98409803 1.06658308 4.13468165 1.43183856
#> [126] 1.14460591 1.62000488 1.72714023 1.29943373 2.45416520
#> [131] 1.39534538 1.55503242 10.17514164 1.39641877 0.96021880
#> [136] 6.23700638 0.96829628 2.18752443 1.04595150 1.03916368
#> [141] 0.31373842 1.91360632 1.04706871 1.27716309 12.45134113
#> [146] 0.55571760 1.05842882 1.18931050 1.15449172 2.97966281
#> [151] 0.98788715 2.10243551 1.19538471 1.03388254 0.75967887
#> [156] 0.91491599 0.61518975 5.93053002 1.03114386 2.00393186
#> [161] 22.53256040 0.96323476 7.46971594 1.82261251 2.34808682
#> [166] 2.74248586 1.36414169 3.88879034 1.18696751 1.61579038
#> [171] 1.49307845 88.82655883 1.35422005 0.95393911 1.77950413
#> [176] 0.91590531 2.24707243 3.45363709 3.73011797 11.21433096
#> [181] 1.21093188 2.96206743 5.20979379 1.41870963 4.06028494
#> [186] 0.82568046 1.28401255 1.14938890 6.06026205 7.85370444
#> [191] 1.55187793 2.27710006 1.03309736 1.98248935 2.15382610
#> [196] 1.10177792 1.62369088 2.79315172 2.37038226 0.62678583