[spcopula-commits] r73 - / pkg pkg/R pkg/demo pkg/man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Dec 21 21:36:33 CET 2012
Author: ben_graeler
Date: 2012-12-21 21:36:33 +0100 (Fri, 21 Dec 2012)
New Revision: 73
Added:
pkg/R/empiricalCopula.R
spcopula_1.0.73.tar.gz
spcopula_1.0.73.zip
Removed:
spcopula_1.0.72.tar.gz
spcopula_1.0.72.zip
Modified:
pkg/DESCRIPTION
pkg/NAMESPACE
pkg/R/Classes.R
pkg/R/asCopula.R
pkg/R/returnPeriods.R
pkg/R/vineCopulas.R
pkg/demo/spCopula_estimation.R
pkg/man/genEmpCop.Rd
Log:
- added new class empiricalCopula
- speedup of ~25 for evaluation of an empirical copula
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2012-12-20 17:32:44 UTC (rev 72)
+++ pkg/DESCRIPTION 2012-12-21 20:36:33 UTC (rev 73)
@@ -1,8 +1,8 @@
Package: spcopula
Type: Package
Title: copula driven spatial analysis
-Version: 1.0.72
-Date: 2012-12-20
+Version: 1.0.73
+Date: 2012-12-21
Author: Benedikt Graeler
Maintainer: Benedikt Graeler <ben.graeler at uni-muenster.de>
Description: This package provides a framework to analyse spatial data provided in the format of the spacetime package with copulas. Additionally, support for calculating different multivariate return periods is implemented.
@@ -12,6 +12,7 @@
URL: http://r-forge.r-project.org/projects/spcopula/
Collate:
Classes.R
+ empiricalCopula.R
partialDerivatives.R
cqsCopula.R
asCopula.R
Modified: pkg/NAMESPACE
===================================================================
--- pkg/NAMESPACE 2012-12-20 17:32:44 UTC (rev 72)
+++ pkg/NAMESPACE 2012-12-21 20:36:33 UTC (rev 73)
@@ -14,6 +14,7 @@
export(surGumbelCopula, r90GumbelCopula, r270GumbelCopula)
export(vineCopula)
export(neighbourhood)
+export(empiricalCopula, genEmpCop)
# general functions
export(rankTransform, dependencePlot, unitScatter, univScatter)
@@ -21,7 +22,6 @@
export(dduCopula,ddvCopula)
export(invdduCopula, invddvCopula)
export(qCopula_u)
-export(genEmpCop)
# tweaks
# export(setSizeLim)
@@ -37,9 +37,10 @@
# MRP functions
export(genEmpKenFun, genInvKenFun)
export(kendallRP, criticalLevel, kendallDistribution, getKendallDistr)
+export(criticalPair, criticalTriple)
## classes
-exportClasses(asCopula, cqsCopula, neighbourhood)
+exportClasses(asCopula, cqsCopula, neighbourhood, empiricalCopula)
# wrappers to CDVine
exportClasses(BB1Copula, surBB1Copula, r90BB1Copula, r270BB1Copula)
Modified: pkg/R/Classes.R
===================================================================
--- pkg/R/Classes.R 2012-12-20 17:32:44 UTC (rev 72)
+++ pkg/R/Classes.R 2012-12-21 20:36:33 UTC (rev 73)
@@ -37,6 +37,22 @@
contains = list("copula")
)
+####
+## an empirical copula representation
+
+validEmpCopula <- function(object) {
+ if(ncol(object at sample) != object at dimension)
+ return("Dimension of the copula and the sample do not match.")
+ else
+ return(TRUE)
+}
+
+setClass("empiricalCopula",
+ representation = representation("copula", sample="matrix"),
+ validity = validEmpCopula,
+ contains = list("copula")
+)
+
##
## the spatial copula
##
Modified: pkg/R/asCopula.R
===================================================================
--- pkg/R/asCopula.R 2012-12-20 17:32:44 UTC (rev 72)
+++ pkg/R/asCopula.R 2012-12-21 20:36:33 UTC (rev 73)
@@ -293,4 +293,4 @@
return((a+3*b)/12)
}
-setMethod("rho",signature("asCopula"),tauASC2)
\ No newline at end of file
+setMethod("rho",signature("asCopula"), rhoASC2)
\ No newline at end of file
Added: pkg/R/empiricalCopula.R
===================================================================
--- pkg/R/empiricalCopula.R (rev 0)
+++ pkg/R/empiricalCopula.R 2012-12-21 20:36:33 UTC (rev 73)
@@ -0,0 +1,78 @@
+########################################
+## ##
+## an empirical copula representation ##
+## ##
+########################################
+
+
+# constructor
+empiricalCopula <- function (copula, sample) {
+ new("empiricalCopula", dimension = copula at dimension,
+ parameters = copula at parameters, param.names = copula at param.names,
+ param.lowbnd = copula at param.lowbnd, param.upbnd = copula at param.upbnd,
+ fullname = paste("Empirical copula derived from",copula at fullname),
+ sample=sample)
+}
+
+# simplified constructor
+genEmpCop <- function(copula, sample.size=1e5) {
+ cat("Note: the copula will be empirically evaluated from a sample of size:", sample.size, "\n")
+ empiricalCopula(copula, rCopula(sample.size,copula))
+}
+
+## density, not yet needed and hence not implemented ##
+
+## jcdf ##
+# pempCop <- function(u, copula) {
+# t_data <- t(copula at sample)
+#
+# u <- matrix(u,ncol=copula at dimension)
+#
+# # --/-- make this a C-function?
+# res <- NULL
+# for(i in 1:nrow(u)) {
+# bool <- t_data <= u[i,]
+# for (i in 2:nrow(t_data)) bool[1,] <- bool[1,] * bool[i,]
+# res <- c(res,sum(bool[1,]))
+# }
+# # --//--
+#
+# return(res/ncol(t_data))
+# }
+
+# from package copula
+pempCop.C <- function(u, copula) {
+ .C("RmultCn", as.double(copula at sample), as.integer(nrow(copula at sample)),
+ copula at dimension, as.double(u), as.integer(nrow(u)), as.double(u[,1]),
+ PACKAGE="copula")[[6]]
+}
+
+##
+# us <- matrix(runif(100),ncol=2)
+# copula <- genEmpCop(normalCopula(.3),1e6)
+#
+# hist(pCopula(us,normalCopula(.3)) - pempCop.C(us,copula),main="C")
+# hist(pCopula(us,normalCopula(.3)) - pempCop(us,copula),main="R")
+#
+# system.time(pempCop.C(us,copula))
+# system.time(pempCop(us,copula))
+
+setMethod("pCopula", signature("numeric", "empiricalCopula"),
+ function(u, copula, ...) {
+ pempCop.C(matrix(u,ncol=copula at dimension),copula)
+ })
+setMethod("pCopula", signature("matrix", "empiricalCopula"), pempCop.C)
+
+
+tauempCop <- function(copula){
+ CDVine:::fasttau(copula at sample[,1],copula at sample[,2])
+}
+
+setMethod("tau",signature("asCopula"),tauempCop)
+
+
+rhoempCop <- function(copula){
+ cor(copula at sample,method="spearman")
+}
+
+setMethod("rho",signature("asCopula"),rhoempCop)
\ No newline at end of file
Modified: pkg/R/returnPeriods.R
===================================================================
--- pkg/R/returnPeriods.R 2012-12-20 17:32:44 UTC (rev 72)
+++ pkg/R/returnPeriods.R 2012-12-21 20:36:33 UTC (rev 73)
@@ -1,15 +1,12 @@
## kendall function (empirical) -> spcopula
genEmpKenFun <- function(copula, sample=NULL) {
if(is.null(sample)) sample <- rCopula(1e6,copula)
- empCop <- genEmpCop(sample)
- ken <- empCop(sample) # takes really long, any suggestions? Comparring a 1e6x3/1e6x2 matrix by 1e6 pairs/triplets values
+ # as empirical copula:
+ # copula <- genEmpCop(copula, sample)
+ ken <- pCopula(sample, copula)
empKenFun <- function(tlevel) {
- res <- NULL
- for(t in tlevel) {
- res <- c(res,sum(ken<=t))
- }
- return(res/nrow(sample))
+ sapply(tlevel,function(t) sum(ken<=t))/nrow(sample)
}
return(empKenFun)
}
@@ -53,28 +50,48 @@
}
## next: calculating critical layer, sampling from the layer, selecting "typical" points
+# calculate critical layer (ONLY 2D by now)
+criticalPair <- function(copula, cl, u, ind) {
+
+ optimFun <- function(x, u, ind) {
+ pair <- cbind(x,x)
+ pair[,ind] <- u
+ return(abs(pCopula(pair,copula)-cl))
+ }
+
+ sapply(u, function(uRow) {
+ upper <- cl+(1-uRow)
+ if (upper<cl | uRow < cl)
+ return(NA)
+ if (upper == cl)
+ return(cl)
+ optimize(function(x) optimFun(x, uRow, ind),
+ interval=c(cl,upper))$minimum
+ })
+}
+
# calculate critical layer (ONLY 3D by now)
-criticalTriple <- function(empCop, cl, u, ind, eps=10e-6) {
+criticalTriple <- function(copula, cl, u, ind) {
if(!is.matrix(u)) u <- matrix(u,ncol=2)
optimFun <- function(x, u, ind) {
x <- matrix(rep(x,3),ncol=3)
x[,ind[1]] <- u[1]
x[,ind[2]] <- u[2]
- return((empCop(x)-cl)^2*1e6)
+ return(abs(pCopula(x,copula)-cl))
}
- res <- apply(u, 1,
- function(uRow) {
- upper <- min(1,2+cl-sum(uRow)) # hyperplane in the hypercube
- if (upper < cl | any(uRow < cl)) return(NA)
- if (upper == cl) return(cl)
- opt <- optimize(function(x) optimFun(x, uRow, ind), interval=c(cl,upper))
- if(opt$objective < eps*1e6) return(opt$minimum)
- else return(NA)
- })
- return(res)
+ apply(u, 1,
+ function(uRow) {
+ upper <- min(1,2+cl-sum(uRow)) # hyperplane in the hypercube
+ if (upper < cl | any(uRow < cl))
+ return(NA)
+ if (upper == cl)
+ return(cl)
+ optimize(function(x) optimFun(x, uRow, ind),
+ interval=c(cl,upper))$minimum
+ })
}
@@ -83,10 +100,8 @@
qCopula_u.def <- function(copula,p,u) { # sample=NULL
dim <- copula at dimension
if(length(p) != length(u)) stop("Length of p and u differ!")
-# if(is.null(sample)) sample <- rCopula(1e5,copula)
-# empCop <- genEmpCop(sample)
+
params <- NULL
-
for(i in 1:length(p)) { # i <- 1
if (u[i] < p[i]) {
params <- rbind(params,rep(NA,dim-1))
@@ -95,12 +110,10 @@
params <- rbind(params,
optimize(function(v) abs(pCopula(cbind(rep(u[i],length(v)),v),copula)-p[i]),
c(p,1))$minimum)
- # function (empCop(cbind(u[i],v))-p[i])^2
} else {
opt <- optim(par=rep(p[i],dim-1),
function(vw) abs(pCopula(c(u[i],vw), copula)-p[i]),
lower=rep(p[i],dim-1), upper=rep(1,dim-1), method="L-BFGS-B")
- # function(vw) (empCop(c(u[i],vw))-p[i])^2
params <- rbind(params, opt$par)
}
}
Modified: pkg/R/vineCopulas.R
===================================================================
--- pkg/R/vineCopulas.R 2012-12-20 17:32:44 UTC (rev 72)
+++ pkg/R/vineCopulas.R 2012-12-21 20:36:33 UTC (rev 73)
@@ -159,31 +159,10 @@
setMethod("dCopula", signature("matrix","vineCopula"), dvineCopula)
## jcdf ##
-genEmpCop <- function(data) {
- t_data <- t(data)
-
- empCop <- function(u) {
- u <- matrix(u,ncol=nrow(t_data))
-
-# --/-- make this a C-function?
- res <- NULL
- for(i in 1:nrow(u)) {
- bool <- t_data <= u[i,]
- for (i in 2:nrow(t_data)) bool[1,] <- bool[1,] * bool[i,]
- res <- c(res,sum(bool[1,]))
- }
-# --//--
-
- return(res/ncol(t_data))
- }
- return(empCop)
-}
-
pvineCopula <- function(u, copula) {
- cat("Note: the copula is empirically evaluated from 100.000 samples.")
- empCop <- genEmpCop(rcopula(copula,1e5))
+ empCop <- genEmpCop(copula,1e5)
- return(empCop(u))
+ return(pCopula(u, empCop))
}
setMethod("pCopula", signature("numeric","vineCopula"), pvineCopula)
Modified: pkg/demo/spCopula_estimation.R
===================================================================
--- pkg/demo/spCopula_estimation.R 2012-12-20 17:32:44 UTC (rev 72)
+++ pkg/demo/spCopula_estimation.R 2012-12-21 20:36:33 UTC (rev 73)
@@ -51,4 +51,4 @@
claytonCopula(0), claytonCopula(0), claytonCopula(0),
claytonCopula(0), indepCopula()),
distances=bins$meanDists,
- spDepFun=calcKTauPol, unit="m")
\ No newline at end of file
+ spDepFun=calcKTauPol, unit="m")
Modified: pkg/man/genEmpCop.Rd
===================================================================
--- pkg/man/genEmpCop.Rd 2012-12-20 17:32:44 UTC (rev 72)
+++ pkg/man/genEmpCop.Rd 2012-12-21 20:36:33 UTC (rev 73)
@@ -7,25 +7,28 @@
Generates an empirical copula from a sample and returns the corresponding function.
}
\usage{
-genEmpCop(data)
+genEmpCop(copula, sample.size=1e5)
}
\arguments{
- \item{data}{
-The sample to be used as input for the empirical copula.
+ \item{copula}{
+The theoretical copula to be represnted numerically.
}
+\item{sample.size}{
+The length of the sample to be used. The default is 1e5.
}
+}
\value{
-The empirical copula as a function (the multivariate cdf).
+An empirical copula class \code{\linkS4class{empiricalCopula}}.
}
\author{
Benedikt Graeler
}
\examples{
-empCop <- genEmpCop(rCopula(500, frankCopula(0.7)))
+empCop <- genEmpCop(frankCopula(0.7), 500)
# the empirical value
-empCop(c(0.3, 0.5))
+pCopula(c(0.3, 0.5), empCop)
# the theoretical value
pCopula(c(0.3, 0.5), frankCopula(0.7))
Deleted: spcopula_1.0.72.tar.gz
===================================================================
(Binary files differ)
Deleted: spcopula_1.0.72.zip
===================================================================
(Binary files differ)
Added: spcopula_1.0.73.tar.gz
===================================================================
(Binary files differ)
Property changes on: spcopula_1.0.73.tar.gz
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: spcopula_1.0.73.zip
===================================================================
(Binary files differ)
Property changes on: spcopula_1.0.73.zip
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
More information about the spcopula-commits
mailing list