Simulate from gamBiCop-class object

Simulate from gamBiCop-class object

gamBiCopSimulate(
  object,
  newdata = NULL,
  N = NULL,
  return.calib = FALSE,
  return.par = FALSE,
  return.tau = FALSE
)

Arguments

object

gamBiCop-class object.

newdata

(same as in predict.gam from the mgcv package) A matrix or data frame containing the values of the model covariates at which simulations are required. If this is not provided then simulations corresponding to the original data are returned.

N

sample size.

return.calib

should the calibration function (TRUE) be returned or not (FALSE)?

return.par

should the copula parameter (TRUE) be returned or not (FALSE)?

return.tau

should the Kendall's tau (TRUE) be returned or not (FALSE)?

Value

A list with 1 item data. When N is smaller or larger than the newdata's number of rows (or the number of rows in the original data if newdata is not provided), then N observations are sampled uniformly (with replacement) among the row of newdata (or the rows of the original data if newdata is not provided).

If return.calib = TRUE, return.par = TRUE and/or return.tau = TRUE, then the list also contains respectively items calib, par and/or tau.

Examples

require(copula)
set.seed(1)

## Simulation parameters (sample size, correlation between covariates,
## Gaussian copula family)
n <- 5e2
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.415   1.284   0.292 

## Extract the gamBiCop objects and show various methods
(res <- fit$res)
#> Gaussian copula with tau(z) = (exp(z)-1)/(exp(z)+1) where 
#> z ~ s(x1, k = basis[1], bs = "cr") + s(x2, k = basis[2], bs = "cr") + 
#>     s(x3, k = basis[3], bs = "cr")
EDF(res)
#> [1] 1.000000 1.994559 5.652399 6.866272
sim <- gamBiCopSimulate(fit$res, X)