[spcopula-commits] r72 - / pkg pkg/R pkg/demo pkg/man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Dec 20 18:32:45 CET 2012
Author: ben_graeler
Date: 2012-12-20 18:32:44 +0100 (Thu, 20 Dec 2012)
New Revision: 72
Added:
pkg/R/spCopula.R
pkg/R/stCopula.R
pkg/demo/spCopula_estimation.R
pkg/man/composeSpCopula.Rd
pkg/man/fitSpCopula.Rd
pkg/man/getNeighbours.Rd
pkg/man/neighbourhood-class.Rd
pkg/man/neighbourhood.Rd
pkg/man/qCopula_u.Rd
pkg/man/stCopula-class.Rd
spcopula_1.0.72.tar.gz
spcopula_1.0.72.zip
Removed:
pkg/R/spcopula.R
pkg/R/stcopula.R
pkg/demo/spcopula_estimation.R
spcopula_1.0.71.tar.gz
spcopula_1.0.71.zip
Modified:
pkg/DESCRIPTION
pkg/NAMESPACE
pkg/R/Classes.R
pkg/R/partialDerivatives.R
pkg/R/returnPeriods.R
pkg/R/spatialPreparation.R
pkg/demo/00Index
pkg/man/invdduCopula-methods.Rd
pkg/man/invddvCopula-methods.Rd
pkg/man/loglikByCopulasLags.Rd
pkg/man/spCopula.Rd
Log:
- few bug fixes
- documentation is now complete!
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2012-12-18 12:10:32 UTC (rev 71)
+++ pkg/DESCRIPTION 2012-12-20 17:32:44 UTC (rev 72)
@@ -1,11 +1,11 @@
Package: spcopula
Type: Package
Title: copula driven spatial analysis
-Version: 1.0.71
-Date: 2012-12-18
+Version: 1.0.72
+Date: 2012-12-20
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 multivariate return periods is implemented.
+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)
@@ -15,8 +15,8 @@
partialDerivatives.R
cqsCopula.R
asCopula.R
- spcopula.R
- stcopula.R
+ spCopula.R
+ stCopula.R
spatialPreparation.R
linkingCDVine.R
BB1copula.R
Modified: pkg/NAMESPACE
===================================================================
--- pkg/NAMESPACE 2012-12-18 12:10:32 UTC (rev 71)
+++ pkg/NAMESPACE 2012-12-20 17:32:44 UTC (rev 72)
@@ -29,8 +29,9 @@
# spatial
export(getNeighbours)
export(calcBins)
+
# fitting
-export(fitCorFun, loglikByCopulasLags, composeSpCop, fitSpCopula)
+export(fitCorFun, loglikByCopulasLags, fitSpCopula, composeSpCopula)
export(spCopula)
# MRP functions
Modified: pkg/R/Classes.R
===================================================================
--- pkg/R/Classes.R 2012-12-18 12:10:32 UTC (rev 71)
+++ pkg/R/Classes.R 2012-12-20 17:32:44 UTC (rev 72)
@@ -57,7 +57,7 @@
# assings valid parameters to the copulas involved
validSpCopula <- function(object) {
- if (length(object at components) != length(object at distances)) return("Length of components + 1 does not equal length of distances. \n Note: The last distance must give the range and it is automatically associated with the indepenence copula.")
+ if (length(object at components) != length(object at distances)) return("Length of components does not equal length of distances. \n Note: The last distance must give the range and it is automatically associated with the indepenence copula.")
check.upper <- NULL
check.lower <- NULL
Modified: pkg/R/partialDerivatives.R
===================================================================
--- pkg/R/partialDerivatives.R 2012-12-18 12:10:32 UTC (rev 71)
+++ pkg/R/partialDerivatives.R 2012-12-20 17:32:44 UTC (rev 72)
@@ -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/>.
-##
-#################################################################################
-
# partial derivatives and their inverse of some copulas from the copula package
# new defined copulas store their partial derivative separately
Modified: pkg/R/returnPeriods.R
===================================================================
--- pkg/R/returnPeriods.R 2012-12-18 12:10:32 UTC (rev 71)
+++ pkg/R/returnPeriods.R 2012-12-20 17:32:44 UTC (rev 72)
@@ -1,6 +1,6 @@
## kendall function (empirical) -> spcopula
genEmpKenFun <- function(copula, sample=NULL) {
- if(is.null(sample)) sample <- rcopula(copula,1e6)
+ 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
@@ -80,21 +80,27 @@
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) { # sample=NULL
dim <- copula at dimension
if(length(p) != length(u)) stop("Length of p and u differ!")
- if(is.null(sample)) sample <- rcopula(copula,1e6)
- empCop <- genEmpCop(sample)
+# if(is.null(sample)) sample <- rCopula(1e5,copula)
+# empCop <- genEmpCop(sample)
params <- NULL
- for(i in 1:length(p)) {
+ for(i in 1:length(p)) { # i <- 1
if (u[i] < p[i]) {
params <- rbind(params,rep(NA,dim-1))
} else {
if (dim == 2) {
- params <- rbind(params,optimize(function(v) (empCop(cbind(u[i],v))-p[i])^2,c(p,1))$minimum)
+ 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) (empCop(c(u[i],vw))-p[i])^2, lower=rep(p[i],dim-1), upper=rep(1,dim-1), method="L-BFGS-B")
+ 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)
}
}
Copied: pkg/R/spCopula.R (from rev 69, pkg/R/spcopula.R)
===================================================================
--- pkg/R/spCopula.R (rev 0)
+++ pkg/R/spCopula.R 2012-12-20 17:32:44 UTC (rev 72)
@@ -0,0 +1,593 @@
+# constructor
+# dimension = "integer" set to 2
+# parameters = "numeric" set of parameters
+# param.names = "character" appropriate names
+# param.lowbnd = "numeric" appropriate lower bounds
+# param.upbnd = "numeric" appropriate upper bounds
+# fullname = "character" messgae printed with "show"
+# components="list" list of copulas (will be automatically supplemented
+# by the independent copula)
+# distances="numeric" the linking distances + the range (will be assigned
+# to the independent copula)
+# unit="character" measurement unit of distance
+# depFun="function" a optional spatial dependence function providing
+# Kendalls tau or Spearman's rho to calib* or exact
+# parameters
+
+
+spCopula <- function(components, distances, spDepFun, unit="m") {
+ if (missing(spDepFun)) {
+ calibMoa <- function(copula, h) return(NULL)
+ } else {
+ if (is.na(match(spDepFun(NULL), c("kendall","spearman","id"))))
+ stop("spDepFun(NULL) must return 'spearman', 'kendall' or 'id'.")
+ cat("The parameters of the components will be recalculated according to the provided spDepFun where possible. \nIn case no 1-1 relation is known, the copula as in components is used. \n")
+ calibMoa <- switch(spDepFun(NULL),
+ kendall=function(copula, h) iTau(copula, spDepFun(h)),
+ spearman=function(copula, h) iRho(copula, spDepFun(h)),
+ id=function(copula, h) return(h))
+
+ for (i in 1:length(components)) {
+ param <- try(calibMoa(components[[i]], distances[i]),T)
+ if (class(param) == "numeric")
+ components[[i]]@parameters[1] <- param # take care of non single parameter copulas
+ }
+ }
+
+# components <- append(components,indepCopula())
+
+ param <- unlist(lapply(components, function(x) x at parameters))
+ param.names <- unlist(lapply(components, function(x) x at param.names))
+ param.low <- unlist(lapply(components, function(x) x at param.lowbnd))
+ param.up <- unlist(lapply(components, function(x) x at param.upbnd))
+
+ new("spCopula", dimension=as.integer(2), parameters=param, param.names=param.names,
+ param.lowbnd=param.low, param.upbnd=param.up,
+ fullname="Spatial Copula: distance dependent convex combination of bivariate copulas",
+ components=components, distances=distances, calibMoa=calibMoa, unit=unit)
+}
+
+## show method
+showCopula <- function(object) {
+ cat(object at fullname, "\n")
+ cat("Dimension: ", object at dimension, "\n")
+ cat("Copulas:\n")
+ for (i in 1:length(object at components)) {
+ cmpCop <- object at components[[i]]
+ cat(" ", cmpCop at fullname, "at", object at distances[i],
+ paste("[",object at unit,"]",sep=""), "\n")
+ if (length(cmpCop at parameters) > 0) {
+ for (i in (1:length(cmpCop at parameters)))
+ cat(" ", cmpCop at param.names[i], " = ", cmpCop at parameters[i], "\n")
+ }
+ }
+ if(!is.null(object at calibMoa(normalCopula(0),0))) cat("A spatial dependence function is used. \n")
+}
+
+setMethod("show", signature("spCopula"), showCopula)
+
+## spatial copula jcdf ##
+
+
+# for spatial copulas with a spatial dependece function
+spDepFunCop <- function(fun, copula, pairs, h) {
+ dists <- copula at distances
+ n.dists <- length(dists)
+ calibPar <- copula at calibMoa
+
+ res <- numeric(0)
+ sel <- which(h < dists[1])
+ if(sum(sel)>0) {
+ 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 (n.dists >= 2) {
+ for ( i in 2:n.dists ) {
+ low <- dists[i-1]
+ high <- dists[i]
+ sel <- which(h >= low & h < high)
+ if (sum(sel)>0) {
+ tmpH <- h[sel]
+ tmpPairs <- pairs[sel,,drop=FALSE]
+ lowerCop <- copula at components[[i-1]]
+ upperCop <- copula at components[[i]]
+ if (class(lowerCop) != class(upperCop)) {
+ 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])
+ 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)
+ } 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)
+ }
+ }
+ }
+ }
+
+ 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]),]))
+ }
+
+ return(res)
+}
+
+# for spatial copulas with a spatial dependece function and a single distance but many pairs
+spDepFunCopSnglDist <- function(fun, copula, pairs, h) {
+ 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 (n.dists >= 2) {
+ for ( i in 2:n.dists ) {
+ low <- dists[i-1]
+ high <- dists[i]
+ if (h >= low & h < high) {
+ 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)
+
+ lowerVals <- fun(pairs, lowerCop)
+ upperVals <- fun(pairs, upperCop)
+
+ res <- (high-h)/(high-low)*lowerVals + (h-low)/(high-low)*upperVals
+ } else {
+ lowerCop at parameters <- calibPar(lowerCop, h)
+ res <- fun(pairs, lowerCop)
+ }
+ }
+ }
+ }
+
+ if(h >= dists[n.dists]) {
+ res <- fun(pairs, copula at components[[n.dists]])
+ }
+
+ return(res)
+}
+
+
+# for static convex combinations of copulas
+spConCop <- function(fun, copula, pairs, h) {
+ 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]])
+ }
+
+ if (n.dists >= 2) {
+ for ( i in 2:n.dists ) {
+ low <- dists[i-1]
+ high <- dists[i]
+ sel <- which(h >= low & h < high)
+ if (sum(sel)>0) {
+ tmpH <- h[sel]
+ tmpPairs <- pairs[sel,,drop=FALSE]
+
+ 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
+ }
+ }
+ }
+
+ sel <- which(h >= dists[n.dists])
+ if(sum(sel)>0) {
+ res[sel] <- fun(pairs[which(h >= dists[n.dists]),],
+ copula at components[[n.dists]])
+ }
+
+ return(res)
+}
+
+
+# u
+# three column matrix providing the transformed pairs and their respective
+# separation distances
+pSpCopula <- function (u, copula) {
+ if (!is.list(u) || !length(u)>=2) stop("Point pairs need to be provided with their separating distance as a list.")
+
+ pairs <- u[[1]]
+ if(!is.matrix(pairs)) pairs <- matrix(pairs,ncol=2)
+ n <- nrow(pairs)
+
+ if(length(u)==3) {
+ block <- u[[3]]
+ if (n%%block != 0) stop("The block size is not a multiple of the data length:",n)
+ } else block <- 1
+
+ h <- u[[2]]
+ if(length(h)>1 && length(h)!=nrow(u[[1]])) {
+ stop("The distance vector must either be of the same length as rows in the data pairs or a single value.")
+ }
+
+
+ if(is.null(copula at calibMoa(normalCopula(0),0))) {
+ res <- spConCop(pCopula, copula, pairs, rep(h,length.out=nrow(pairs)))
+ } else {
+ if(length(h)>1) {
+ if (block == 1){
+ ordering <- order(h)
+
+ # ascending sorted pairs allow for easy evaluation
+ pairs <- pairs[ordering,,drop=FALSE]
+ h <- h[ordering]
+
+ res <- spDepFunCop(pCopula, copula, pairs, h)
+
+ # reordering the values
+ res <- res[order(ordering)]
+ } else {
+ res <- NULL
+ for(i in 1:(n%/%block)) {
+ res <- c(res, spDepFunCopSnglDist(pCopula, copula,
+ pairs[((i-1)*block+1):(i*block),],
+ h[i*block]))
+ }
+ }
+ } else {
+ res <- spDepFunCopSnglDist(pCopula, copula, pairs, h)
+ }
+ }
+
+ return(res)
+}
+
+setMethod("pCopula", signature("list","spCopula"), pSpCopula)
+
+## spatial Copula density ##
+
+# u
+# three column matrix providing the transformed pairs and their respective
+# separation distances
+dSpCopula <- function (u, copula) {
+ if (!is.list(u) || !length(u)>=2) stop("Point pairs need to be provided with their separating distance as a list.")
+
+ pairs <- u[[1]]
+ n <- nrow(pairs)
+
+ if(length(u)==3) {
+ block <- u[[3]]
+ if (n%%block != 0) stop("The block size is not a multiple of the data length:",n)
+ } else block <- 1
+
+ h <- u[[2]]
+ if(length(h)>1 && length(h)!=nrow(u[[1]])) {
+ stop("The distance vector must either be of the same length as rows in the data pairs or a single value.")
+ }
+
+ if(is.null(copula at calibMoa(normalCopula(0),0))){
+ res <- spConCop(dCopula, copula, pairs, rep(h, length.out=nrow(pairs)))
+ }
+ else {
+ if(length(h)>1) {
+ if (block == 1){
+ ordering <- order(h)
+
+ # ascending sorted pairs allow for easy evaluation
+ pairs <- pairs[ordering,,drop=FALSE]
+ h <- h[ordering]
+
+ res <- spDepFunCop(dCopula, copula, pairs, h)
+
+ # reordering the values
+ res <- res[order(ordering)]
+ } else {
+ res <- NULL
+ for(i in 1:(n%/%block)) {
+ res <- c(res, spDepFunCopSnglDist(dCopula, copula,
+ pairs[((i-1)*block+1):(i*block),],
+ h[i*block]))
+ }
+ }
+ } else {
+ res <- spDepFunCopSnglDist(dCopula, copula, pairs, h)
+ }
+ }
+
+ return(res)
+}
+
+setMethod("dCopula", signature("list","spCopula"), dSpCopula)
+
+
+## partial derivatives ##
+
+## dduSpCopula
+###############
+
+dduSpCopula <- function (u, copula) {
+ if (!is.list(u) || !length(u)>=2) stop("Point pairs need to be provided with their separating distance as a list.")
+
+ pairs <- u[[1]]
+ n <- nrow(pairs)
+
+ if(length(u)==3) {
+ block <- u[[3]]
+ if (n%%block != 0) stop("The block size is not a multiple of the data length:",n)
+ } else block <- 1
+
+ h <- u[[2]]
+ if(length(h)>1 && length(h)!=nrow(u[[1]])) {
+ stop("The distance vector must either be of the same length as rows in the data pairs or a single value.")
+ }
+
+ if(is.null(copula at calibMoa(normalCopula(0),0))) res <- spConCop(dduCopula, copula, pairs,
+ rep(h,length.out=nrow(pairs)))
+ else {
+ if(length(h)>1) {
+ if (block == 1){
+ ordering <- order(h)
+
+ # ascending sorted pairs allow for easy evaluation
+ pairs <- pairs[ordering,,drop=FALSE]
+ h <- h[ordering]
+
+ res <- spDepFunCop(dduCopula, copula, pairs, h)
+
+ # reordering the values
+ res <- res[order(ordering)]
+ } else {
+ res <- NULL
+ for(i in 1:(n%/%block)) {
+ res <- c(res, spDepFunCopSnglDist(dduCopula, copula,
+ pairs[((i-1)*block+1):(i*block),],
+ h[i*block]))
+ }
+ }
+ } else {
+ res <- spDepFunCopSnglDist(dduCopula, copula, pairs, h)
+ }
+ }
+
+ return(res)
+}
+
+setMethod("dduCopula", signature("list","spCopula"), dduSpCopula)
+
+## ddvSpCopula
+###############
+
+ddvSpCopula <- function (u, copula) {
+ if (!is.list(u) || !length(u)>=2) stop("Point pairs need to be provided with their separating distance as a list.")
+
+ pairs <- u[[1]]
+ n <- nrow(pairs)
+
+ if(length(u)==3) {
+ block <- u[[3]]
+ if (n%%block != 0) stop("The block size is not a multiple of the data length:",n)
+ } else block <- 1
+
+ h <- u[[2]]
+ if(length(h)>1 && length(h)!=nrow(u[[1]])) {
+ stop("The distance vector must either be of the same length as rows in the data pairs or a single value.")
+ }
+
+
+ if(is.null(copula at calibMoa(normalCopula(0),0))) res <- spConCop(ddvCopula, copula, pairs,
+ rep(h,length.out=nrow(pairs)))
+ else {
+ if(length(h)>1) {
+ if (block == 1){
+ ordering <- order(h)
+
+ # ascending sorted pairs allow for easy evaluation
+ pairs <- pairs[ordering,,drop=FALSE]
+ h <- h[ordering]
+
+ res <- spDepFunCop(ddvCopula, copula, pairs, h)
+
+ # reordering the values
+ res <- res[order(ordering)]
+ } else {
+ res <- NULL
+ for(i in 1:(n%/%block)) {
+ res <- c(res, spDepFunCopSnglDist(ddvCopula, copula,
+ pairs[((i-1)*block+1):(i*block),],
+ h[i*block]))
+ }
+ }
+ } else {
+ res <- spDepFunCopSnglDist(ddvCopula, copula, pairs, h)
+ }
+ }
+
+ return(res)
+}
+
+setMethod("ddvCopula", signature("list","spCopula"), ddvSpCopula)
+
+
+#############
+## ##
+## FITTING ##
+## ##
+#############
+
+# two models:
+# 1) Kendall's tau driven:
+# fit curve through emp. Kendall's tau values, identify validity ranges for
+# copula families deriving parameters from the fit, fade from one family to
+# another at borders
+# 2) convex-linear combination of copulas:
+# fit one per lag, fade from one to another
+
+# towards the first model:
+
+# INPUT: the stBinning
+# steps
+# a) fit a curve
+# b) estimate bivariate copulas per lag (limited to those with some 1-1-relation
+# to Kendall's tau')
+# INTERMEDIATE RESULT
+# c) select best fits based on ... e.g. log-likelihood, visual inspection
+# d) compose bivariate copulas to one spatial copula
+# OUTPUT: a spatial copula parametrised by distance through Kendall's tau
+
+# towards a)
+# bins -> typically output from calcBins
+# degree -> the degree of the polynominal
+# cutoff -> maximal distance that should be considered for fitting
+# bounds -> the bounds of the correlation function (typically c(0,1))
+# method -> the measure of association, either "kendall" or "spearman"
+fitCorFun <- function(bins, degree=3, cutoff=NA, bounds=c(0,1),
+ cor.method=NULL) {
+ if(is.null(cor.method)) {
+ if(is.null(attr(bins,"cor.method")))
+ stop("Neither the bins arguments has an attribute cor.method nor is the parameter cor.method provided.")
+ else
+ cor.method <- attr(bins,"cor.method")
+ } else
+ if(!is.null(attr(bins,"cor.method")) && cor.method != attr(bins,"cor.method"))
+ stop("The cor.method attribute of the bins argument and the argument cor.method do not match.")
+
+ bins <- as.data.frame(bins[1:2])
+ if(!is.na(cutoff)) bins <- bins[which(bins[[1]] <= cutoff),]
+
+ fitCor <- lm(lagCor ~ poly(meanDists, degree), data = bins)
+
+ print(fitCor)
+ cat("Sum of squared residuals:",sum(fitCor$residuals^2),"\n")
+
+ function(x) {
+ if (is.null(x)) return(cor.method)
+ return(pmin(bounds[2], pmax(bounds[1],
+ eval(predict(fitCor, data.frame(meanDists=x))))))
+ }
+}
+
+
+# towards b)
+# bins -> typically output from calcBins
+# calcTau -> a function on distance providing Kendall's tau estimates,
+# families -> a vector of dummy copula objects of each family to be considered
+# DEFAULT: c(normal, t_df=4, clayton, frank, gumbel
+loglikByCopulasLags <- function(bins, calcCor,
+ families=c(normalCopula(0),
+ tCopula(0,dispstr="un"),
+ claytonCopula(0), frankCopula(1),
+ gumbelCopula(1))) {
+ moa <- switch(calcCor(NULL),
+ kendall=function(copula, h) iTau(copula, calcCor(h)),
+ spearman=function(copula, h) iRho(copula, calcCor(h)),
+ id=function(copula, h) calcCor(h))
+ loglik <- NULL
+ for (cop in families) {
+ print(cop)
+ tmploglik <- NULL
+ 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))))
+ }
+ loglik <- cbind(loglik, tmploglik)
+ }
+
+ colnames(loglik) <- sapply(families, function(x) class(x)[1])
+
+ return(loglik)
+}
+
+# towards d)
+composeSpCopula <- function(bestFit, families, bins, calcCor, range=max(bins$meanDists)) {
+ nFits <- length(bestFit)
+ if(nFits > length(bins$meanDists))
+ stop("There may not be less bins than best fits.\n")
+ rangeIndex <- min(nFits, max(which(bins$meanDists <= range)))
+
+ if (missing(calcCor)) {
+ return(spCopula(components = as.list(families[bestFit[1:rangeIndex]]),
+ distances = bins$meanDists[1:rangeIndex],
+ unit = "m"))
+ }
+
+ else {
+ rangeIndex <- min(rangeIndex, which(calcCor(bins$meanDists) <= 0))
+
+ return(spCopula(components = as.list(families[bestFit[1:rangeIndex]]),
+ distances = bins$meanDists[1:rangeIndex],
+ unit = "m", spDepFun = calcCor))
+ }
+}
+
+# older implementation:
+# composeSpCop <- function(bestFit, families, bins, calcCor) {
+# nfits <- length(bestFit)
+# gaps <- which(diff(bestFit)!=0)
+#
+# if(missing(calcCor)) noCor <- nfits
+# else noCor <- min(which(calcCor(bins$meanDists)<=0), nfits)
+#
+# breaks <- sort(c(gaps, gaps+1, noCor))
+# breaks <- breaks[breaks<noCor]
+#
+# cops <- as.list(families[bestFit[breaks]])
+#
+# breaks <- unique(c(breaks, min(nfits,rev(breaks)[1]+1)))
+# distances <- bins$meanDists[breaks]
+#
+# if(length(breaks)==length(cops)) {
+# warning("Lags do not cover the entire range with positive correlation. ",
+# "An artifical one fading towards the indepCopula has been added.")
+# distances <- c(distances,
+# rev(distances)[1]+diff(bins$meanDists[nfits+c(-1,0)]))
+# }
+#
+# if(missing(calcCor)) return(spCopula(components=cops, distances=distances,
+# unit="m"))
+# else return(spCopula(components=cops, distances=distances, unit="m",
+# spDepFun=calcCor))
+# }
+
+# in once
+
+# bins -> typically output from calcBins
+# cutoff -> maximal distance that should be considered for fitting
+# families -> a vector of dummy copula objects of each family to be considered
+# DEFAULT: c(normal, t_df=4, clayton, frank, gumbel
+# ...
+# type -> the type of curve (by now only polynominals are supported)
+# degree -> the degree of the polynominal
+# bounds -> the bounds of the correlation function (typically c(0,1))
+# method -> the measure of association, either "kendall" or "spearman"
+fitSpCopula <- function(bins, cutoff=NA,
+ families=c(normalCopula(0),
+ tCopula(0,dispstr="un"), claytonCopula(0),
+ frankCopula(1), gumbelCopula(1)), ...) {
+ calcCor <- fitCorFun(bins, cutoff=cutoff, ...)
+ loglik <- loglikByCopulasLags(bins, calcCor, families)
+
+ bestFit <- apply(apply(loglik, 1, rank),2,
+ function(x) which(x==length(families)))
+
+ return(composeSpCopula(bestFit, families, bins, calcCor, range=cutoff))
+}
Modified: pkg/R/spatialPreparation.R
===================================================================
--- pkg/R/spatialPreparation.R 2012-12-18 12:10:32 UTC (rev 71)
+++ pkg/R/spatialPreparation.R 2012-12-20 17:32:44 UTC (rev 72)
@@ -1,24 +1,3 @@
-#################################################################################
-##
-## R package spcopula by Benedikt Graeler 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/>
-##
-#################################################################################
-
##################################################################
## ##
## dedicated functions based on sp preparing the use of copulas ##
@@ -33,7 +12,7 @@
# index="matrix" linking the obs. in data to the coordinates
neighbourhood <- function(data, distances, sp, index){
- varNames <- names(data[[1]])
+ varNames <- names(data)[[1]]
sizeN <- ncol(distances)+1
data <- as.data.frame(data)
colnames(data) <- paste(paste("N",rep(0:(sizeN-1),each=length(varNames)),sep=""),rep(varNames,sizeN),sep=".")
@@ -118,7 +97,7 @@
dists <- rbind(dists, c(nbrs[,1]))
}
-return(neighbourhood(lData, dists, SpatialPoints(spData), index))
+return(neighbourhood(as.data.frame(lData), dists, SpatialPoints(spData), index))
}
## testing ##
Deleted: pkg/R/spcopula.R
===================================================================
--- pkg/R/spcopula.R 2012-12-18 12:10:32 UTC (rev 71)
+++ pkg/R/spcopula.R 2012-12-20 17:32:44 UTC (rev 72)
@@ -1,563 +0,0 @@
-# constructor
-# dimension = "integer" set to 2
-# parameters = "numeric" set of parameters
-# param.names = "character" appropriate names
-# param.lowbnd = "numeric" appropriate lower bounds
-# param.upbnd = "numeric" appropriate upper bounds
-# fullname = "character" messgae printed with "show"
-# components="list" list of copulas (will be automatically supplemented
-# by the independent copula)
-# distances="numeric" the linking distances + the range (will be assigned
-# to the independent copula)
-# unit="character" measurement unit of distance
-# depFun="function" a optional spatial dependence function providing
-# Kendalls tau or Spearman's rho to calib* or exact
-# parameters
-
-
-spCopula <- function(components, distances, spDepFun, unit="m") {
- if (missing(spDepFun)) {
- calibMoa <- function(copula, h) return(NULL)
- } else {
- if (is.na(match(spDepFun(NULL),c("kendall","spearman","id","none")))) stop("spDepFun(NULL) must return 'spearman', 'kendall' or 'id'.")
- cat("The parameters of the components will be recalculated according to the provided spDepFun. \n")
- calibMoa <- switch(spDepFun(NULL),
- kendall=function(copula, h) iTau(copula, spDepFun(h)),
- spearman=function(copula, h) iRho(copula, spDepFun(h)),
- id=function(copula, h) return(h))
-
- for (i in 1:length(components)) {
- components[[i]]@parameters[1] <- calibMoa(components[[i]], distances[i]) # take care of non single parameter copulas
- }
- }
-
- components <- append(components,indepCopula())
-
- param <- unlist(lapply(components, function(x) x at parameters))
- param.names <- unlist(lapply(components, function(x) x at param.names))
- param.low <- unlist(lapply(components, function(x) x at param.lowbnd))
- param.up <- unlist(lapply(components, function(x) x at param.upbnd))
-
- new("spCopula", dimension=as.integer(2), parameters=param, param.names=param.names,
- param.lowbnd=param.low, param.upbnd=param.up,
- fullname="Spatial Copula: distance dependent convex combination of bivariate copulas",
- components=components, distances=distances, calibMoa=calibMoa, unit=unit)
-}
-
-## show method
-showCopula <- function(object) {
- cat(object at fullname, "\n")
- cat("Dimension: ", object at dimension, "\n")
- cat("Copulas:\n")
- for (i in 1:length(object at components)) {
- cmpCop <- object at components[[i]]
- cat(" ", cmpCop at fullname, "at", object at distances[i],
- paste("[",object at unit,"]",sep=""), "\n")
- if (length(cmpCop at parameters) > 0) {
- for (i in (1:length(cmpCop at parameters)))
- cat(" ", cmpCop at param.names[i], " = ", cmpCop at parameters[i], "\n")
- }
- }
- if(!is.null(object at calibMoa(normalCopula(0),0))) cat("A spatial dependence function is used. \n")
-}
-
-setMethod("show", signature("spCopula"), showCopula)
-
-## spatial copula jcdf ##
-
-
-# for spatial copulas with a spatial dependece function
-spDepFunCop <- function(fun, copula, pairs, h) {
- dists <- copula at distances
- n.dists <- length(dists)
- calibPar <- copula at calibMoa
-
- res <- numeric(0)
- sel <- which(h < dists[1])
- if(sum(sel)>0) {
- 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 (n.dists >= 2) {
- for ( i in 2:n.dists ) {
- low <- dists[i-1]
- high <- dists[i]
- sel <- which(h >= low & h < high)
- if (sum(sel)>0) {
- tmpH <- h[sel]
- tmpPairs <- pairs[sel,,drop=FALSE]
- lowerCop <- copula at components[[i-1]]
- upperCop <- copula at components[[i]]
- if (class(lowerCop) != class(upperCop)) {
- 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])
- 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)
- } 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)
- }
- }
- }
- }
-
- sel <- which(h >= dists[n.dists])
- if(sum(sel)>0) {
- res <- c(res, fun(copula at components[[n.dists]],
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/spcopula -r 72
More information about the spcopula-commits
mailing list