[spcopula-commits] r74 - / pkg pkg/R pkg/demo pkg/man

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Fri Dec 28 17:05:22 CET 2012


Author: ben_graeler
Date: 2012-12-28 17:05:21 +0100 (Fri, 28 Dec 2012)
New Revision: 74

Added:
   pkg/man/criticalPair.Rd
   pkg/man/criticalTriple.Rd
   pkg/man/empiricalCopula-class.Rd
   pkg/man/empiricalCopula.Rd
   spcopula_1.0.74.tar.gz
   spcopula_1.0.74.zip
Removed:
   spcopula_1.0.73.tar.gz
   spcopula_1.0.73.zip
Modified:
   pkg/DESCRIPTION
   pkg/R/asCopula.R
   pkg/R/cqsCopula.R
   pkg/R/empiricalCopula.R
   pkg/R/partialDerivatives.R
   pkg/R/returnPeriods.R
   pkg/R/spCopula.R
   pkg/R/stCopula.R
   pkg/demo/MRP.R
   pkg/demo/spCopula_estimation.R
   pkg/man/genEmpCop.Rd
   pkg/man/genInvKenFun.Rd
Log:
- additional documentation for empiricalCopula, criticalPair and criticalTriple
- empiricalCopula can now be used without a theoretical copula family
- implentation of log support for densities of spCopula

Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION	2012-12-21 20:36:33 UTC (rev 73)
+++ pkg/DESCRIPTION	2012-12-28 16:05:21 UTC (rev 74)
@@ -1,14 +1,14 @@
 Package: spcopula
 Type: Package
 Title: copula driven spatial analysis
-Version: 1.0.73
-Date: 2012-12-21
+Version: 1.0.74
+Date: 2012-12-28
 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.
 License: GPL-2
 LazyLoad: yes
-Depends: copula (>= 0.99-2), spacetime, CDVine, methods, lattice, R (>= 2.13.2)
+Depends: copula (>= 0.999-5), spacetime (>= 1.0-2), CDVine, methods, lattice, R (>= 2.15.0)
 URL: http://r-forge.r-project.org/projects/spcopula/
 Collate:
   Classes.R

Modified: pkg/R/asCopula.R
===================================================================
--- pkg/R/asCopula.R	2012-12-21 20:36:33 UTC (rev 73)
+++ pkg/R/asCopula.R	2012-12-28 16:05:21 UTC (rev 74)
@@ -208,7 +208,7 @@
            copula = copula))
 }
 
-fitASC2.moa <- function(moa, data, method="itau") {
+fitASC2.moa <- function(moa, data, method="itau", tol=.Machine$double.eps^.5) {
   smpl <- as.matrix(data)
 
   iTau <- function(p) {
@@ -229,7 +229,7 @@
     return(res)
   }
 
-  b <- optimize(sec,c(-1,1))$minimum
+  b <- optimize(sec,c(-1,1), tol=tol)$minimum
 
   param <- c(iFun(b),b)
 

Modified: pkg/R/cqsCopula.R
===================================================================
--- pkg/R/cqsCopula.R	2012-12-21 20:36:33 UTC (rev 73)
+++ pkg/R/cqsCopula.R	2012-12-28 16:05:21 UTC (rev 74)
@@ -252,7 +252,7 @@
 ))
 }
 
-fitCQSec.moa <- function(moa, data, method="itau") {
+fitCQSec.moa <- function(moa, data, method="itau", tol=.Machine$double.eps^.5) {
 smpl <- as.matrix(data)
 
 iTau <- function(p) {
@@ -273,7 +273,7 @@
 return(res)
 }
 
-b <- optimize(sec,c(-1,1))$minimum
+b <- optimize(sec,c(-1,1), tol=tol)$minimum
 
 param <- c(iFun(b),b)
 

Modified: pkg/R/empiricalCopula.R
===================================================================
--- pkg/R/empiricalCopula.R	2012-12-21 20:36:33 UTC (rev 73)
+++ pkg/R/empiricalCopula.R	2012-12-28 16:05:21 UTC (rev 74)
@@ -4,9 +4,18 @@
 ##                                    ##
 ########################################
 
-
 # constructor
-empiricalCopula <- function (copula, sample) {
+empiricalCopula <- function (sample=NULL, copula) {
+  if(is.null(sample) && missing(copula))
+    stop("At least one parameter of copula or sample must be provided.")
+  if(is.null(sample))
+    return(genEmpCop(copula))
+  if(missing(copula))
+    return(new("empiricalCopula", dimension = as.integer(ncol(sample)), 
+               parameters = as.numeric(NA), param.names = "unknown", 
+               param.lowbnd = as.numeric(NA), param.upbnd = as.numeric(NA), 
+               fullname = "Unkown empirical copula based on a sample.",
+               sample=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, 
@@ -16,10 +25,11 @@
 
 # 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))
+  cat("Note: the copula will be empirically represented by a sample of size:", sample.size, "\n")
+  empiricalCopula(rCopula(sample.size,copula), copula)
 }
 
+
 ## density, not yet needed and hence not implemented ##
 
 ## jcdf ##

Modified: pkg/R/partialDerivatives.R
===================================================================
--- pkg/R/partialDerivatives.R	2012-12-21 20:36:33 UTC (rev 73)
+++ pkg/R/partialDerivatives.R	2012-12-28 16:05:21 UTC (rev 74)
@@ -6,30 +6,30 @@
 
 ## inverse partial derivatives 
 # numerical standard function
-invdduCopula <- function(u, copula, y) {
+invdduCopula <- function(u, copula, y, tol=.Machine$double.eps^0.5) {
     if (length(u) != length(y)) 
         stop("Length of u and y differ!")
-    warning("Numerical evaluation of invddu takes place.")
+    message("Numerical evaluation of invddu takes place.")
     res <- NULL
     for (i in 1:length(u)) {
         res <- rbind(res, optimize( function(x) 
           (dduCopula(cbind(rep(u[i], length(x)), x),copula) - y[i])^2, 
-            interval = c(0, 1))$minimum)
+            interval = c(0, 1), tol=tol)$minimum)
     }
     return(res)
 }
 
 setGeneric("invdduCopula")
 
-invddvCopula <- function(v, copula, y) {
+invddvCopula <- function(v, copula, y, tol=.Machine$double.eps^0.5) {
     if (length(v) != length(y)) 
         stop("Length of v and y differ!")
-  warning("Numerical evaluation of invddv takes place.")
+  message("Numerical evaluation of invddv takes place.")
     res <- NULL
     for (i in 1:length(v)) {
         res <- rbind(res, optimize(function(x) 
           (ddvCopula(cbind(x, rep(v[i], length(x))),copula) - y[i])^2, 
-            interval = c(0, 1))$minimum)
+            interval = c(0, 1), tol=tol)$minimum)
     }
     return(res)
 }

Modified: pkg/R/returnPeriods.R
===================================================================
--- pkg/R/returnPeriods.R	2012-12-21 20:36:33 UTC (rev 73)
+++ pkg/R/returnPeriods.R	2012-12-28 16:05:21 UTC (rev 74)
@@ -22,11 +22,11 @@
 # <=> K_C(t) = 1 - \mu /KRP
 # <=> t = K_C^{-1}(1 - \mu /KRP)
 
-genInvKenFun <- function(kenFun, ...) {
+genInvKenFun <- function(kenFun, tol=.Machine$double.eps^.5) {
   invKenFun <- function(k){
     res <- NULL
     for(i in 1:length(k)) {
-      res <- c(res, optimize(function(x) (kenFun(x)-k[i])^2,c(0,1))$minimum)
+      res <- c(res, optimize(function(x) (kenFun(x)-k[i])^2, c(0,1), tol=tol)$minimum)
     }
     return(res)
   }
@@ -51,7 +51,7 @@
 
 ## next: calculating critical layer, sampling from the layer, selecting "typical" points
 # calculate critical layer (ONLY 2D by now)
-criticalPair <- function(copula, cl, u, ind) {
+criticalPair <- function(copula, cl, u, ind, tol=sqrt(.Machine$double.eps)) {
   
   optimFun <- function(x, u, ind) {
     pair <- cbind(x,x)
@@ -66,13 +66,13 @@
               if (upper == cl) 
                 return(cl)
               optimize(function(x) optimFun(x, uRow, ind),
-                       interval=c(cl,upper))$minimum
+                       interval=c(cl,upper), tol=tol)$minimum
             })
 }
 
 
 # calculate critical layer (ONLY 3D by now)
-criticalTriple <- function(copula, cl, u, ind) {
+criticalTriple <- function(copula, cl, u, ind, tol=sqrt(.Machine$double.eps)) {
   if(!is.matrix(u)) u <- matrix(u,ncol=2)
     
   optimFun <- function(x, u, ind) {
@@ -90,14 +90,14 @@
           if (upper == cl) 
             return(cl)
           optimize(function(x) optimFun(x, uRow, ind), 
-                   interval=c(cl,upper))$minimum
+                   interval=c(cl,upper),tol=tol)$minimum
         })
 }
 
 
 setGeneric("qCopula_u",function(copula,p,u,...) {standardGeneric("qCopula_u")})
 
-qCopula_u.def <- function(copula,p,u) { # sample=NULL
+qCopula_u.def <- function(copula,p,u, tol=.Machine$double.eps^.5) { # sample=NULL
   dim <- copula at dimension
   if(length(p) != length(u)) stop("Length of p and u differ!")
   
@@ -109,7 +109,7 @@
       if (dim == 2) {
         params <- rbind(params, 
                         optimize(function(v) abs(pCopula(cbind(rep(u[i],length(v)),v),copula)-p[i]),
-                                 c(p,1))$minimum)
+                                 c(p,1), tol=tol)$minimum)
       } else {
         opt <- optim(par=rep(p[i],dim-1), 
                      function(vw) abs(pCopula(c(u[i],vw), copula)-p[i]), 

Modified: pkg/R/spCopula.R
===================================================================
--- pkg/R/spCopula.R	2012-12-21 20:36:33 UTC (rev 73)
+++ pkg/R/spCopula.R	2012-12-28 16:05:21 UTC (rev 74)
@@ -70,7 +70,7 @@
 
 
 # for spatial copulas with a spatial dependece function
-spDepFunCop <- function(fun, copula, pairs, h) {
+spDepFunCop <- function(fun, copula, pairs, h, do.logs=F, ...) {
   dists <- copula at distances
   n.dists <- length(dists)
   calibPar <- copula at calibMoa
@@ -81,9 +81,13 @@
     tmpH <- h[sel]
     tmpCop <- copula at components[[1]]
     tmpPairs <- pairs[sel,,drop=FALSE]
-    for (j in 1:length(tmpH)) {
-      tmpCop at parameters[1] <- calibPar(tmpCop, tmpH[j])
-      res <- c(res, fun(tmpPairs[j,], tmpCop))
+    if(class(tmpCop) == "indepCopula")
+      res <- fun(tmpPairs, tmpCop, ...)
+    else {
+      for (j in 1:length(tmpH)) {
+        tmpCop at parameters[1] <- calibPar(tmpCop, tmpH[j])
+        res <- c(res, fun(tmpPairs[j,], tmpCop, ...))
+      }
     }
   }
   
@@ -101,17 +105,27 @@
           lowerVals <- numeric(0)
           upperVals <- numeric(0)
           for (j in 1:length(tmpH)) {
-            lowerCop at parameters[1] <- calibPar(lowerCop,  tmpH[j])
-            upperCop at parameters[1] <- calibPar(upperCop, tmpH[j])
+            if (class(lowerCop) != "indepCopula")
+              lowerCop at parameters[1] <- calibPar(lowerCop,  tmpH[j])
+            if (class(upperCop) != "indepCopula")
+              upperCop at parameters[1] <- calibPar(upperCop, tmpH[j])
             lowerVals <- c(lowerVals, fun(tmpPairs[j,], lowerCop))
             upperVals <- c(upperVals, fun(tmpPairs[j,], upperCop))
           }
-          res <- c(res,(high-tmpH)/(high-low)*lowerVals+(tmpH-low)/(high-low)*upperVals)
+          if(do.logs)
+            res <- c(res,
+                     log((high-tmpH)/(high-low)*lowerVals+(tmpH-low)/(high-low)*upperVals))
+          else 
+            res <- c(res,(high-tmpH)/(high-low)*lowerVals+(tmpH-low)/(high-low)*upperVals)
         } else {
-          newVals <- numeric(0)
-          for (j in 1:length(tmpH)) {
-            lowerCop at parameters <- calibPar(lowerCop, tmpH[j])
-            newVals <- c(newVals, fun(tmpPairs[j,], lowerCop))
+          if (class(lowerCop) == "indepCopula")
+            newVals <- fun(tmpPairs, lowerCop, ...)
+          else {
+            newVals <- numeric(0)
+            for (j in 1:length(tmpH)) {
+              lowerCop at parameters <- calibPar(lowerCop, tmpH[j])
+              newVals <- c(newVals, fun(tmpPairs[j,], lowerCop, ...))
+            }
           }
           res <- c(res, newVals)
         }
@@ -121,23 +135,24 @@
   
   sel <- which(h >= dists[n.dists])
   if(sum(sel)>0) {
-    res <- c(res, fun(copula at components[[n.dists]], 
-                      pairs[which(h >= dists[n.dists]),]))
+    res <- c(res, fun(pairs[which(h >= dists[n.dists]),],
+                      copula at components[[n.dists]],, ...))
   }
 
   return(res)
 }
 
 # for spatial copulas with a spatial dependece function and a single distance but many pairs
-spDepFunCopSnglDist <- function(fun, copula, pairs, h) {
+spDepFunCopSnglDist <- function(fun, copula, pairs, h, do.logs=F, ...) {
   dists <- copula at distances
   n.dists <- length(dists)
   calibPar <- copula at calibMoa
 
   if(h < dists[1]) {
     tmpCop <- copula at components[[1]]
-    tmpCop at parameters[1] <- calibPar(tmpCop, h)
-    res <- fun(pairs, tmpCop)
+    if(class(tmpCop) != "indepCopula")
+      tmpCop at parameters[1] <- calibPar(tmpCop, h)
+    res <- fun(pairs, tmpCop, ...)
   }
   
   if (n.dists >= 2) {
@@ -148,23 +163,28 @@
         lowerCop <- copula at components[[i-1]]
         upperCop <- copula at components[[i]]
         if (class(lowerCop) != class(upperCop)) {
-          lowerCop at parameters[1] <- calibPar(lowerCop, h)
-          upperCop at parameters[1] <- calibPar(upperCop, h)
+          if(class(lowerCop) != "indepCopula")
+            lowerCop at parameters[1] <- calibPar(lowerCop, h)
+          if(class(upperCop) != "indepCopula")
+            upperCop at parameters[1] <- calibPar(upperCop, h)
 
           lowerVals <- fun(pairs, lowerCop)
           upperVals <- fun(pairs, upperCop)
 
           res <- (high-h)/(high-low)*lowerVals + (h-low)/(high-low)*upperVals
+          if(do.logs)
+            res <- log(res)
         } else {
-          lowerCop at parameters <- calibPar(lowerCop, h)
-          res <- fun(pairs, lowerCop)
+          if(class(tmpCop) != "indepCopula")
+            lowerCop at parameters <- calibPar(lowerCop, h)
+          res <- fun(pairs, lowerCop, ...)
         }
       }
     }
   }
   
   if(h >= dists[n.dists]) {
-    res <- fun(pairs, copula at components[[n.dists]])
+    res <- fun(pairs, copula at components[[n.dists]], ...)
   }
   
   return(res)
@@ -172,14 +192,14 @@
 
 
 # for static convex combinations of copulas
-spConCop <- function(fun, copula, pairs, h) {
+spConCop <- function(fun, copula, pairs, h, do.logs=F, ...) {
   dists <- copula at distances
   n.dists <- length(dists)
   
   res <- numeric(nrow(pairs))
   sel <- which(h < dists[1])
   if(sum(sel)>0) {
-    res[sel] <- fun(pairs[sel,,drop=FALSE],copula at components[[1]])
+    res[sel] <- fun(pairs[sel,,drop=FALSE],copula at components[[1]],...)
   }
   
   if (n.dists >= 2) {
@@ -194,7 +214,10 @@
         lowerVals <- fun(tmpPairs[,], copula at components[[i-1]])
         upperVals <- fun(tmpPairs[,], copula at components[[i]])
 
-        res[sel] <- (high-tmpH)/(high-low)*lowerVals+(tmpH-low)/(high-low)*upperVals
+        if(do.logs)
+          res[sel] <- log((high-tmpH)/(high-low)*lowerVals+(tmpH-low)/(high-low)*upperVals)
+        else
+          res[sel] <- (high-tmpH)/(high-low)*lowerVals+(tmpH-low)/(high-low)*upperVals
       }
     }
   }
@@ -202,7 +225,7 @@
   sel <- which(h >= dists[n.dists])
   if(sum(sel)>0) {
     res[sel] <- fun(pairs[which(h >= dists[n.dists]),],
-                    copula at components[[n.dists]])
+                    copula at components[[n.dists]], ...)
   }
 
   return(res)
@@ -268,7 +291,7 @@
 # u 
 #   three column matrix providing the transformed pairs and their respective 
 #   separation distances
-dSpCopula <- function (u, copula) {
+dSpCopula <- function (u, copula, log=F) {
   if (!is.list(u) || !length(u)>=2) stop("Point pairs need to be provided with their separating distance as a list.")
   
   pairs <- u[[1]]
@@ -285,7 +308,8 @@
   }
   
   if(is.null(copula at calibMoa(normalCopula(0),0))){
-    res <- spConCop(dCopula, copula, pairs, rep(h, length.out=nrow(pairs)))
+    res <- spConCop(dCopula, copula, pairs, rep(h, length.out=nrow(pairs)), 
+                    do.logs=log, log=log)
   }
   else {
     if(length(h)>1) {
@@ -296,7 +320,7 @@
         pairs <- pairs[ordering,,drop=FALSE] 
         h <- h[ordering]
         
-        res <- spDepFunCop(dCopula, copula, pairs, h)
+        res <- spDepFunCop(dCopula, copula, pairs, h, do.logs=log, log=log)
         
         # reordering the values
         res <- res[order(ordering)]
@@ -305,20 +329,19 @@
         for(i in 1:(n%/%block)) {
           res <- c(res, spDepFunCopSnglDist(dCopula, copula, 
                                             pairs[((i-1)*block+1):(i*block),], 
-                                            h[i*block]))
+                                            h[i*block], do.logs=log, log=log))
         }
       }
     } else {
-      res <- spDepFunCopSnglDist(dCopula, copula, pairs, h)
+      res <- spDepFunCopSnglDist(dCopula, copula, pairs, h, do.logs=log, log=log)
     }
   }
   
   return(res)
 }
 
-setMethod("dCopula", signature("list","spCopula"), dSpCopula)
+setMethod(dCopula, signature("list","spCopula"), dSpCopula)
 
-
 ## partial derivatives ##
 
 ## dduSpCopula
@@ -507,7 +530,7 @@
     for(i in 1:length(bins$meanDists)) {
       if(class(cop)!="indepCopula")
         cop at parameters[1] <- moa(cop, bins$meanDists[i])
-      tmploglik <- c(tmploglik, sum(log(dCopula(bins$lagData[[i]],cop))))
+      tmploglik <- c(tmploglik, sum(dCopula(bins$lagData[[i]],cop, log=T)))
     }
     loglik <- cbind(loglik, tmploglik)
   }

Modified: pkg/R/stCopula.R
===================================================================
--- pkg/R/stCopula.R	2012-12-21 20:36:33 UTC (rev 73)
+++ pkg/R/stCopula.R	2012-12-28 16:05:21 UTC (rev 74)
@@ -1,24 +1,3 @@
-#################################################################################
-##
-##   R package spCopula by Benedikt Gräler Copyright (C) 2011
-##
-##   This file is part of the R package spCopula.
-##
-##   The R package spCopula 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.
-##
-##   The R package spCopula 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 the R package spCopula. If not, see <http://www.gnu.org/licenses/>.
-##
-#################################################################################
-
 # constructor
 # dimension = "integer"     set to 2
 # parameters = "numeric"    set of parameters

Modified: pkg/demo/MRP.R
===================================================================
--- pkg/demo/MRP.R	2012-12-21 20:36:33 UTC (rev 73)
+++ pkg/demo/MRP.R	2012-12-28 16:05:21 UTC (rev 74)
@@ -10,31 +10,33 @@
 cor(triples,method="kendall")
 
 # estiamte the BB7 copula by means of maximum likelihood
-copQV <- fitCopula(BB7Copula(param=c(2,14)), peakVol, method="ml",start=c(2,14), estimate.variance=F)@copula
+copQV <- fitCopula(BB7Copula(param=c(2,14)), peakVol, method="ml", 
+                   start=c(2,14), estimate.variance=F)@copula
 copQV
 
-# we use a design return period of 100 years
-# the MAR-case: given the 0.99 quantile of one marginal, what is the 
+# we use a design return period of 10 years
+# the MAR-case: given the 0.9 quantile of one marginal, what is the 
 # corresponding quantile of the second variable with the given joint return 
-# period of 1/(1-0.99)
-v_MAR <- c(0.99,invdduCopula(0.99,copQV,0.99))
+# period of 1/(1-0.9)
+v_MAR <- c(0.9, invdduCopula(0.9, copQV, 0.9))
 v_MAR
+pCopula(v_MAR, copQV)
 
 ## the anlytical kendall distribution
 kendallFunQV <- getKendallDistr(copQV)
 
 # the Kendall distribution value (fraction of pairs having a smaller copula value than "t")
-kendallFunQV(t=0.99)
+kendallFunQV(t=0.9)
 
 curve(kendallFunQV, from=0, to=1, asp=1)
 
-# the critical level of the KEN2-RP for 100 years
-t_KEN2 <- criticalLevel(kendallFunQV,100,mu=1)
+# the critical level of the KEN2-RP for 10 years
+t_KEN2 <- criticalLevel(kendallFunQV, 10, mu=1)
 t_KEN2
 
 # the corresponmding KEN2-RP for the OR-RP
-kendallRP(kendallFunQV,cl=0.99,mu=1,copula=copQV)
+kendallRP(kendallFunQV, cl=0.9, mu=1, copula=copQV)
 
 # illustrating the critical lines (empirically)
-contour(copQV,pCopula,levels=c(0.99,t_KEN2),
-        xlim=c(0.98,1), ylim=c(0.98,1), n=1000, asp=1, col="blue")
\ No newline at end of file
+contour(copQV,pCopula,levels=c(0.9,t_KEN2),
+        xlim=c(0.8,1), ylim=c(0.8,1), n=100, asp=1, col="blue")
\ No newline at end of file

Modified: pkg/demo/spCopula_estimation.R
===================================================================
--- pkg/demo/spCopula_estimation.R	2012-12-21 20:36:33 UTC (rev 73)
+++ pkg/demo/spCopula_estimation.R	2012-12-28 16:05:21 UTC (rev 74)
@@ -9,7 +9,8 @@
 spplot(meuse,"zinc", col.regions=bpy.colors(5))
 
 ## margins ##
-hist(dataSet[["zinc"]],freq=F,n=30,ylim=c(0,0.0035))
+hist(dataSet[["zinc"]],freq=F,n=30,ylim=c(0,0.0035), 
+     main="Histogram of zinc", xlab="zinc concentration")
 gevEsti <- fgev(dataSet[["zinc"]])$estimate
 loc <- gevEsti[1]
 scale <- gevEsti[2]
@@ -26,7 +27,7 @@
 bins <- calcBins(dataSet,var="zinc",nbins=10,cutoff=800)
 
 # transform data to the unit interval
-bins$lagData <- lapply(bins$lagData, function(x) cbind(rank(x[,1])/(nrow(x)+1),rank(x[,2])/(nrow(x)+1)))
+bins$lagData <- lapply(bins$lagData, rankTransform)
 
 ## calculate parameters for Kendall's tau function ##
 # either linear
@@ -39,16 +40,32 @@
 
 ## find best fitting copula per lag class
 loglikTau <- loglikByCopulasLags(bins, calcKTauPol,
-                                 families=c(normalCopula(0), tCopula(0,dispstr = "un"),
+                                 families=c(normalCopula(0), tCopula(0),
                                             claytonCopula(0), frankCopula(1), 
                                             gumbelCopula(1), joeBiCopula(1.5),
                                             indepCopula()))
-bestFitTau <- apply(apply(loglikTau, 1, rank),2,function(x) which(x==7))
+bestFitTau <- apply(apply(loglikTau, 1, rank, na.last=T), 2, function(x) which(x==7))
+bestFitTau
 
 ## set-up a spatial Copula ##
-spCop <- spCopula(components=list(normalCopula(0), tCopula(0, dispstr = "un"),
+spCop <- spCopula(components=list(normalCopula(0), tCopula(0),
                                   frankCopula(1), normalCopula(0), claytonCopula(0),
                                   claytonCopula(0), claytonCopula(0), claytonCopula(0),
                                   claytonCopula(0), indepCopula()),
                   distances=bins$meanDists,
                   spDepFun=calcKTauPol, unit="m")
+
+## compare spatial copula loglik by lag:
+spLoglik <- NULL
+for(i in 1:length(bins$lags)) { # i <- 8
+  spLoglik <- c(spLoglik,
+                sum(dCopula(list(u=bins$lagData[[i]], h=bins$lags[[i]][,3]), 
+                                spCop,log=T)))
+}
+
+plot(spLoglik, ylab="log-likelihood", xlim=c(1,11)) 
+points(loglikTau[cbind(1:10,bestFitTau)], col="green", pch=16)
+points(loglikTau[,1], col="red", pch=5)
+legend(6, 50,c("Spatial Copula", "best copula per lag", "Gaussian Copula","number of pairs"), 
+       pch=c(1,16,5,50), col=c("black", "green", "red"))
+text(x=(1:10+0.5),y=spLoglik,lapply(bins$lagData,length))
\ No newline at end of file

Added: pkg/man/criticalPair.Rd
===================================================================
--- pkg/man/criticalPair.Rd	                        (rev 0)
+++ pkg/man/criticalPair.Rd	2012-12-28 16:05:21 UTC (rev 74)
@@ -0,0 +1,49 @@
+\name{criticalPair}
+\alias{criticalPair}
+\title{
+Calculate Crtitical Pairs
+}
+\description{
+For a given argument, the second argument is calculated for a given critical level (cumulated probability level) and copula.
+}
+\usage{
+criticalPair(copula, cl, u, ind, tol=sqrt(.Machine$double.eps))
+}
+\arguments{
+  \item{copula}{
+The three-dimesnional copula.
+}
+  \item{cl}{
+The critical level (cumulative probability level).
+}
+  \item{u}{
+The two given arguments.
+}
+  \item{ind}{
+The index of the given arguments.
+}
+  \item{tol}{
+  The tolerance value as used by \code{\link{optimise}}
+}
+}
+\value{
+The second argument for which the geiven critical level (cumulative probability level) is achieved.
+}
+\author{
+Benedikt Graeler
+}
+\note{
+\code{\link{optimise}} is used to find the third argument.
+}
+
+\seealso{
+\code{\link{criticalTriple}}
+}
+\examples{
+v <- criticalPair(frankCopula(0.7), 0.9, u=.97, 1)
+pCopula(c(0.97, v),frankCopula(0.7))
+}
+% Add one or more standard keywords, see file 'KEYWORDS' in the
+% R documentation directory.
+\keyword{ ~kwd1 }
+\keyword{ ~kwd2 }% __ONLY ONE__ keyword per line

Added: pkg/man/criticalTriple.Rd
===================================================================
--- pkg/man/criticalTriple.Rd	                        (rev 0)
+++ pkg/man/criticalTriple.Rd	2012-12-28 16:05:21 UTC (rev 74)
@@ -0,0 +1,49 @@
+\name{criticalTriple}
+\alias{criticalTriple}
+\title{
+calculate critical triples
+}
+\description{
+For two given arguments, the third argument is calculated for a given critical level (cumulated probability level) and copula.
+}
+\usage{
+criticalTriple(copula, cl, u, ind, tol=sqrt(.Machine$double.eps))
+}
+\arguments{
+  \item{copula}{
+The three-dimesnional copula.
+}
+  \item{cl}{
+The critical level (cumulative probability level).
+}
+  \item{u}{
+The two given arguments.
+}
+  \item{ind}{
+The index of the given arguments.
+}
+\item{tol}{
+  The tolerance value as used by \code{\link{optimise}}
+}
+}
+\value{
+The third argument for which the geiven critical level (cumulative probability level) is achieved.
+}
+\author{
+Benedikt Graeler
+}
+\note{
+\code{\link{optimise}} is used to find the third argument.
+}
+\seealso{
+\code{\link{criticalPair}}
+}
+\examples{
+w <- criticalTriple(frankCopula(0.7,dim=3), 0.9, c(.97,.97), c(1,2))
+
+# check the triple
+pCopula(c(0.97, 0.97, w), frankCopula(0.7, dim=3))
+
+}
+\keyword{ multivariate }
+\keyword{ distribution }

Added: pkg/man/empiricalCopula-class.Rd
===================================================================
--- pkg/man/empiricalCopula-class.Rd	                        (rev 0)
+++ pkg/man/empiricalCopula-class.Rd	2012-12-28 16:05:21 UTC (rev 74)
@@ -0,0 +1,44 @@
+\name{empiricalCopula-class}
+\Rdversion{1.1}
+\docType{class}
+\alias{empiricalCopula-class}
+
+\title{Class \code{"empiricalCopula"}}
+\description{
+A class representing an empirical version of a given copula.
+}
+\section{Objects from the Class}{
+Objects can be created by calls of the form \code{new("empiricalCopula", ...)}, 
+through the constructor \code{\link{empiricalCopula}} or the simplified constructor \code{\link{genEmpCop}}.
+
+}
+\section{Slots}{
+  \describe{
+    \item{\code{sample}:}{Object of class \code{"matrix"} representing the sample. }
+    \item{\code{dimension}:}{Object of class \code{"integer"} giving the dimension. }
+    \item{\code{parameters}:}{Object of class \code{"numeric"} giving the parameters of the underlying copula family if present, \"NA\" otherwise.}
+    \item{\code{param.names}:}{Object of class \code{"character"} giving the parameter names of the underlying copula family if present, \"unknown\" otherwise. }
+    \item{\code{param.lowbnd}:}{Object of class \code{"numeric"}  giving the lower bound of the underlying copula family if present, \"NA\" otherwise.}
+    \item{\code{param.upbnd}:}{Object of class \code{"numeric"}  giving the upper bound of the underlying copula family if present, \"NA\" otherwise.}
+    \item{\code{fullname}:}{Object of class \code{"character"}  giving a descriptive name including the underlying copula family's name if present. }
+  }
+}
+\section{Extends}{
+Class \code{"\linkS4class{copula}"}, directly.
+Class \code{"\linkS4class{Copula}"}, by class "copula", distance 2.
+}
+\section{Methods}{
+No additional methods defined with class "empiricalCopula" in the signature, but \code{\link{pCopula}} works as expected.
+}
+\author{
+Benedikt Graeler}
+\note{
+Its implementation of \code{\link{pCopula}} is based on C-code from \code{\link{copula-package}}.
+}
+
+
+\seealso{\code{\link{genEmpCop}}}
+\examples{
+showClass("empiricalCopula")
+}
+\keyword{classes}

Added: pkg/man/empiricalCopula.Rd
===================================================================
--- pkg/man/empiricalCopula.Rd	                        (rev 0)
+++ pkg/man/empiricalCopula.Rd	2012-12-28 16:05:21 UTC (rev 74)
@@ -0,0 +1,49 @@
+\name{empiricalCopula}
+\alias{empiricalCopula}
+\title{
+Constructor of an empircal copula class
+}
+\description{
+Constructs an object of the empirical copula class \code{\linkS4class{empiricalCopula}}, see \code{\link{genEmpCop}} for a simplified version.
+}
+\usage{
+empiricalCopula(sample=NULL, copula)
+}
+\arguments{
+  \item{sample}{
+A sample from a provided or unknown copula family.
+}
+  \item{copula}{
+The underlying theoretical copula, in case it is known or a sample should be gnerated.
+}
+}
+\value{
+An object of \code{\linkS4class{empiricalCopula}}.
+}
+\author{
+Benedikt Graeler
+}
+\note{
+Its implementation of \code{\link{pCopula}} is based on C-code from \code{\link{copula-package}}.
+}
+
+\seealso{
+\code{\link{genEmpCop}} for a simplified constructor with sample length controle.
+}
+\examples{
+empCop <- empiricalCopula(rCopula(500,frankCopula(0.7)))
+str(empCop)
+
+empCop <- empiricalCopula(copula=frankCopula(0.7))
+str(empCop)
+
+empCop <- empiricalCopula(rCopula(500,frankCopula(0.7)), frankCopula(0.7))
+str(empCop)
+
+# the empirical value
+pCopula(c(0.3, 0.5), empCop)
+
+# the theoretical value
+pCopula(c(0.3, 0.5), frankCopula(0.7))
+}
+\keyword{ multivariate }

Modified: pkg/man/genEmpCop.Rd
===================================================================
--- pkg/man/genEmpCop.Rd	2012-12-21 20:36:33 UTC (rev 73)
+++ pkg/man/genEmpCop.Rd	2012-12-28 16:05:21 UTC (rev 74)
@@ -24,6 +24,10 @@
 Benedikt Graeler
 }
 
+\note{
+Its implementation of \code{\link{pCopula}} is based on C-code from \code{\link{copula-package}}.
+}
+
 \examples{
 empCop <- genEmpCop(frankCopula(0.7), 500)
 

Modified: pkg/man/genInvKenFun.Rd
===================================================================
--- pkg/man/genInvKenFun.Rd	2012-12-21 20:36:33 UTC (rev 73)
+++ pkg/man/genInvKenFun.Rd	2012-12-28 16:05:21 UTC (rev 74)
@@ -7,14 +7,14 @@
 The inverse of a (empirical) Kendall distribution function is generated based on numerical inversion using optimise.
 }
 \usage{
-genInvKenFun(kenFun, ...)
+genInvKenFun(kenFun, tol)
 }
 \arguments{
   \item{kenFun}{
 The (empirical) Kendall distribution function to be inverted.
 }
-  \item{\dots}{
-Control options passed on to \code{\link{optimise}}.
+  \item{tol}{
+Tolerance passed on to \code{\link{optimise}}.
 }
 }
 \value{
@@ -29,10 +29,10 @@
 }
 \examples{
 frankKenDistrFun <- getKendallDistr(frankCopula(.5))
-frankInvKenDitrFun <- genInvKenFun(frankKenDistrFun)
+frankInvKenDistrFun <- genInvKenFun(frankKenDistrFun)
 
-frankInvKenDitrFun(.8)
-frankKenDistrFun(frankInvKenDitrFun(.8))
+frankInvKenDistrFun(.8)
+frankKenDistrFun(frankInvKenDistrFun(.8))
 }
 \keyword{ multivariate }
 \keyword{ distribution }

Deleted: spcopula_1.0.73.tar.gz
===================================================================
(Binary files differ)

Deleted: spcopula_1.0.73.zip
===================================================================
(Binary files differ)

Added: spcopula_1.0.74.tar.gz
===================================================================
(Binary files differ)


Property changes on: spcopula_1.0.74.tar.gz
___________________________________________________________________
Added: svn:mime-type
   + application/octet-stream

Added: spcopula_1.0.74.zip
===================================================================
(Binary files differ)


Property changes on: spcopula_1.0.74.zip
___________________________________________________________________
Added: svn:mime-type
   + application/octet-stream



More information about the spcopula-commits mailing list