inst/doc/spcopula.R

In spcopula: Copula Driven Analysis - Multivariate, Spatial, Spatio-Temporal

### R code from vignette source 'spcopula.Rnw'
### Encoding: UTF-8

###################################################
### code chunk number 1: spcopula.Rnw:52-54
###################################################
options(prompt = "R> ", continue = "+  ", width = 70, useFancyQuotes = FALSE)
set.seed(1331)


###################################################
### code chunk number 2: spcopula.Rnw:129-144
###################################################
library("copula")
contourGauss <- function (copula, ...) {
  n <- 51
  xis <- yis <- seq(-3, 3, len = n)
  uis <- vis <- pnorm(xis)
  grids <- as.matrix(expand.grid(uis, vis, KEEP.OUT.ATTRS = FALSE))
  zmat <- matrix(dCopula(grids, copula)*dnorm(qnorm(grids[,1]))*dnorm(qnorm(grids[,2])), n, n)
  contour(xis, yis, zmat,asp=1, ...)
}

pdf(file="figures/copula_densities.pdf", 6, 3.5)
par(mfrow=c(1,2),mai=c(0.5,0.5,0.5,0)+0.1)
contourGauss(normalCopula(iTau(normalCopula(0), 0.5)), main="Gaussian")
contourGauss(gumbelCopula(iTau(gumbelCopula(2), 0.5)), main="Gumbel")
dev.off()


###################################################
### code chunk number 3: spcopula.Rnw:256-262 (eval = FALSE)
###################################################
## library("spcopula")
## coVarCop <- function(stInd) {
##   week <- min(ceiling(stInd[2] / 7), 52)
##   copulaFromFamilyIndex(weekCop[[week]]$family, weekCop[[week]]$par, 
##                         weekCop[[week]]$par2)
## }


###################################################
### code chunk number 4: spcopula.Rnw:274-275 (eval = FALSE)
###################################################
## stBins <- calcBins(EU_RB_2005, "marPM10", nbins = 40, tlags = -(0:4))


###################################################
### code chunk number 5: spcopula.Rnw:278-279 (eval = FALSE)
###################################################
## stDepFun <- fitCorFun(stBins, rep(3, 5), tlags = -(0:4))


###################################################
### code chunk number 6: spcopula.Rnw:290-293 (eval = FALSE)
###################################################
## families <- c(normalCopula(0), tCopula(0),
##               claytonCopula(0), frankCopula(1), gumbelCopula(1), 
##               joeBiCopula())


###################################################
### code chunk number 7: spcopula.Rnw:296-298 (eval = FALSE)
###################################################
## loglikTau <- loglikByCopulasStLags(stBins, EU_RB_2005, families,
##                                    stDepFun)


###################################################
### code chunk number 8: spcopula.Rnw:305-308 (eval = FALSE)
###################################################
## bestFitTau <- lapply(loglikTau, 
##                      function(x) apply(apply(x$loglik, 1, rank),
##                                        2, which.max))


###################################################
### code chunk number 9: spcopula.Rnw:328-333 (eval = FALSE)
###################################################
## distSelect <- function(x) {
##   stBins$meanDists[sort(unique(c(which(diff(x) != 0),
##                                  which(diff(x) != 0) + 1, 1, 40)))]
## }
## listDists <- lapply(bestFitTau, distSelect)


###################################################
### code chunk number 10: spcopula.Rnw:336-341 (eval = FALSE)
###################################################
## famSelect <- function(x) {
##   families[x[sort(unique(c(which(diff(x) != 0),
##                            which(diff(x) != 0) + 1, 1, 40)))]]
## }
## listCops <- lapply(bestFitTau, famSelect)


###################################################
### code chunk number 11: spcopula.Rnw:346-348 (eval = FALSE)
###################################################
## stBiCop <- stCopula(components = listCops, distances = listDists, 
##                     tlags = -c(0:4), stDepFun = stDepFun)


###################################################
### code chunk number 12: spcopula.Rnw:354-357 (eval = FALSE)
###################################################
## stNeigh <- getStNeighbours(EU_RB_2005, var = "marPM10", spSize = 9, 
##                            tlags = -(0:4), timeSteps = 90, min.dist = 10)
## stRedNeigh <- reduceNeighbours(stNeigh, stDepFun, 9)


###################################################
### code chunk number 13: spcopula.Rnw:363-364 (eval = FALSE)
###################################################
## condData <- dropStTree(stRedNeigh, EU_RB_2005, stBiCop)


###################################################
### code chunk number 14: spcopula.Rnw:369-370 (eval = FALSE)
###################################################
## condCoVa <- condCovariate(stNeigh, coVarCop)


###################################################
### code chunk number 15: spcopula.Rnw:375-377 (eval = FALSE)
###################################################
## secTreeData <- cbind(condCoVa, as.matrix(condData@data))
## vineFit <- fitCopula(vineCopula(10L), secTreeData)


###################################################
### code chunk number 16: spcopula.Rnw:382-383 (eval = FALSE)
###################################################
## stCVVC <- stCoVarVineCopula(coVarCop, stBiCop, vineFit@copula)


###################################################
### code chunk number 17: spcopula.Rnw:388-395 (eval = FALSE)
###################################################
## predNeigh <- getStNeighbours(EU_RB_2005, targetGeom, 
##                              spSize = 9, tlags = -(0:4),
##                              var = "marPM10", coVar = "marEMEP",
##                              prediction = TRUE, min.dist = 10)
## predNeigh <- reduceNeighbours(predNeigh, stDepFun, 9)
## stVinePred <- stCopPredict(predNeigh, stCVVC, list(q = qFun), 
##                            method = "expectation")