[spcopula-commits] r94 - / pkg pkg/R pkg/man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Apr 23 16:17:12 CEST 2013
Author: ben_graeler
Date: 2013-04-23 16:17:11 +0200 (Tue, 23 Apr 2013)
New Revision: 94
Modified:
pkg/DESCRIPTION
pkg/R/Classes.R
pkg/R/partialDerivatives.R
pkg/R/returnPeriods.R
pkg/R/spCopula.R
pkg/R/spVineCopula.R
pkg/R/spatialPreparation.R
pkg/man/condSpVine.Rd
pkg/man/fitCorFun.Rd
pkg/man/spCopPredict.Rd
spcopula_0.1-1.tar.gz
spcopula_0.1-1.zip
Log:
- fine tuning of condSpVine
- weighted estimate in fitCorFun
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2013-04-15 11:39:18 UTC (rev 93)
+++ pkg/DESCRIPTION 2013-04-23 14:17:11 UTC (rev 94)
@@ -2,7 +2,7 @@
Type: Package
Title: copula driven spatial analysis
Version: 0.1-1
-Date: 2013-04-11
+Date: 2013-04-23
Author: Benedikt Graeler
Maintainer: Benedikt Graeler <ben.graeler at uni-muenster.de>
Description: This package provides a framework to analyse via copulas spatial and spatio-temporal data provided in the format of the spacetime package. Additionally, support for calculating different multivariate return periods is implemented.
Modified: pkg/R/Classes.R
===================================================================
--- pkg/R/Classes.R 2013-04-15 11:39:18 UTC (rev 93)
+++ pkg/R/Classes.R 2013-04-23 14:17:11 UTC (rev 94)
@@ -95,8 +95,10 @@
check.upper <- NULL
check.lower <- NULL
+ nComp <- length(object at components)
if(!is.null(object at calibMoa(normalCopula(0),0))) {
- for (i in 1:(length(object at components)-1)) {
+ nonIndep <- sapply(object at components[-nComp], function(x) class(x) != "indepCopula")
+ for (i in (1:(nComp-1))[nonIndep]) {
check.upper <- c(check.upper, is.na(object at calibMoa(object at components[[i]], object at distances[i+1])))
check.lower <- c(check.lower, is.na(object at calibMoa(object at components[[i]], c(0,object at distances)[i])))
}
Modified: pkg/R/partialDerivatives.R
===================================================================
--- pkg/R/partialDerivatives.R 2013-04-15 11:39:18 UTC (rev 93)
+++ pkg/R/partialDerivatives.R 2013-04-23 14:17:11 UTC (rev 94)
@@ -218,7 +218,7 @@
u1 <- u[,1]
u2 <- u[,2]
- pcopula(gumbelCopula(rho),u) * ((-log(u1))^rho+(-log(u2))^rho)^(1/rho-1) * (-log(u1))^(rho-1)/u1
+ pCopula(u,gumbelCopula(rho)) * ((-log(u1))^rho+(-log(u2))^rho)^(1/rho-1) * (-log(u1))^(rho-1)/u1
}
setMethod("dduCopula", signature("numeric","gumbelCopula"),
@@ -237,7 +237,7 @@
u1 <- u[,1]
u2 <- u[,2]
- pcopula(gumbelCopula(rho),u) * ((-log(u2))^rho+(-log(u1))^rho)^(1/rho-1) * (-log(u2))^(rho-1)/u2
+ pCopula(u,gumbelCopula(rho)) * ((-log(u2))^rho+(-log(u1))^rho)^(1/rho-1) * (-log(u2))^(rho-1)/u2
}
setMethod("ddvCopula", signature("numeric","gumbelCopula"),
Modified: pkg/R/returnPeriods.R
===================================================================
--- pkg/R/returnPeriods.R 2013-04-15 11:39:18 UTC (rev 93)
+++ pkg/R/returnPeriods.R 2013-04-23 14:17:11 UTC (rev 94)
@@ -123,4 +123,4 @@
return(cbind(u,params))
}
-setMethod("qCopula_u",signature("copula"),qCopula_u.def)
\ No newline at end of file
+setMethod("qCopula_u", signature("copula"), qCopula_u.def)
Modified: pkg/R/spCopula.R
===================================================================
--- pkg/R/spCopula.R 2013-04-15 11:39:18 UTC (rev 93)
+++ pkg/R/spCopula.R 2013-04-23 14:17:11 UTC (rev 94)
@@ -28,7 +28,8 @@
id=function(copula, h) return(h))
for (i in 1:length(components)) {
- param <- try(calibMoa(components[[i]], distances[i]),T)
+ if(class(components[[i]]) != "indepCopula")
+ param <- try(calibMoa(components[[i]], distances[i]),T)
if (class(param) == "numeric")
components[[i]]@parameters[1] <- param # take care of non single parameter copulas
}
@@ -36,10 +37,26 @@
# 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))
+ param <- unlist(lapply(components,
+ function(x) {
+ if(class(x)=="indepCopula")
+ return(NA)
+ x at parameters}))
+ param.names <- unlist(lapply(components,
+ function(x) {
+ if(class(x)=="indepCopula")
+ return(NA)
+ x at param.names}))
+ param.low <- unlist(lapply(components,
+ function(x) {
+ if(class(x)=="indepCopula")
+ return(NA)
+ x at param.lowbnd}))
+ param.up <- unlist(lapply(components,
+ function(x) {
+ if(class(x)=="indepCopula")
+ return(NA)
+ x at param.upbnd}))
new("spCopula", dimension=as.integer(2), parameters=param, param.names=param.names,
param.lowbnd=param.low, param.upbnd=param.up,
@@ -75,6 +92,7 @@
n.dists <- length(dists)
calibPar <- copula at calibMoa
+ # data is sorted to be ascending in distance
res <- numeric(0)
sel <- which(h < dists[1])
if(sum(sel)>0) {
@@ -136,7 +154,7 @@
sel <- which(h >= dists[n.dists])
if(sum(sel)>0) {
res <- c(res, fun(pairs[which(h >= dists[n.dists]),],
- copula at components[[n.dists]],, ...))
+ copula at components[[n.dists]], ...))
}
return(res)
@@ -152,42 +170,39 @@
tmpCop <- copula at components[[1]]
if(class(tmpCop) != "indepCopula")
tmpCop at parameters[1] <- calibPar(tmpCop, h)
- res <- fun(pairs, tmpCop, ...)
+ return(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)) {
- if(class(lowerCop) != "indepCopula")
- lowerCop at parameters[1] <- calibPar(lowerCop, h)
- if(class(upperCop) != "indepCopula")
- upperCop at parameters[1] <- calibPar(upperCop, h)
+ if(h >= dists[n.dists]) {
+ return(fun(pairs, copula at components[[n.dists]], ...))
+ }
+
+ for (i in 2:n.dists) {
+ low <- dists[i-1]
+ high <- dists[i]
+ if(low <= h & h < high) {
+ lowerCop <- copula at components[[i-1]]
+ upperCop <- copula at components[[i]]
+ if (class(lowerCop) != class(upperCop)) {
+ 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)
+ 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 {
- if(class(lowerCop) != "indepCopula")
- lowerCop at parameters <- calibPar(lowerCop, h)
- res <- fun(pairs, lowerCop, ...)
- }
+ res <- (high-h)/(high-low)*lowerVals + (h-low)/(high-low)*upperVals
+ if(do.logs)
+ return(log(res))
+ return(res)
+ } else {
+ if(class(lowerCop) != "indepCopula")
+ lowerCop at parameters <- calibPar(lowerCop, h)
+ return(fun(pairs, lowerCop, ...))
}
}
}
-
- if(h >= dists[n.dists]) {
- res <- fun(pairs, copula at components[[n.dists]], ...)
- }
-
- return(res)
}
@@ -232,106 +247,63 @@
}
-# u
-# two column matrix providing the transformed pairs and their respective
-# separation distances
-# h providing the separating distance(s)
-# block
-# block distances, pairs are assumed to be ordered block wise
+## spatial copula CDF
+######################
+
pSpCopula <- function (u, copula, h, block=1) {
- if (missing(h)) stop("Point pairs need to be provided with their separating distance \"h\".")
-
- n <- nrow(u)
-
- if (n%%block != 0) stop("The block size is not a multiple of the data length:",n)
-
- if(length(h)>1 && length(h)!=n) {
+ if (missing(h))
+ stop("Point pairs need to be provided with their separating distance \"h\".")
+ if(length(h)>1 && length(h)!=nrow(u))
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)))
+ return(spConCop(pCopula, copula, u, rep(h,length.out=nrow(u))))
- if(is.null(copula at calibMoa(normalCopula(0),0))) {
- res <- spConCop(pCopula, copula, u, rep(h,length.out=nrow(pairs)))
- } else {
- if(length(h)>1) {
- if (block == 1){
- ordering <- order(h)
+ if(length(h)>1) {
+ ordering <- order(h)
- # ascending sorted pairs allow for easy evaluation
- u <- u[ordering,,drop=FALSE]
- h <- h[ordering]
+ # ascending sorted pairs allow for easy evaluation
+ u <- u[ordering,,drop=FALSE]
+ h <- h[ordering]
- res <- spDepFunCop(pCopula, copula, u, h)
+ res <- spDepFunCop(pCopula, copula, u, h)
- # reordering the values
- res <- res[order(ordering)]
- } else {
- res <- NULL
- for(i in 1:(n%/%block)) {
- res <- c(res, spDepFunCopSnglDist(pCopula, copula,
- u[((i-1)*block+1):(i*block),],
- h[i*block]))
- }
- }
- } else {
- res <- spDepFunCopSnglDist(pCopula, copula, u, h)
- }
- }
-
- return(res)
+ # reordering the values
+ return(res[order(ordering)])
+ } else
+ return(spDepFunCopSnglDist(pCopula, copula, u, h))
}
setMethod(pCopula, signature("numeric","spCopula"),
function(u, copula, ...) pSpCopula(matrix(u,ncol=2),copula, ...))
setMethod(pCopula, signature("matrix","spCopula"), pSpCopula)
-## spatial Copula density ##
+## spatial Copula density
+##########################
-# u
-# three column matrix providing the transformed pairs and their respective
-# separation distances
-dSpCopula <- function (u, copula, log, h, block=1) {
- if (missing(h)) stop("Point pairs need to be provided with their separating distance \"h\".")
-
- n <- nrow(u)
-
- if (n%%block != 0) stop("The block size is not a multiple of the data length:",n)
-
- if(length(h)>1 && length(h)!=n) {
+dSpCopula <- function (u, copula, log, h) {
+ if (missing(h))
+ stop("Point pairs need to be provided with their separating distance \"h\".")
+ if(length(h)>1 && length(h)!=nrow(u))
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, u, rep(h, length.out=n),
- do.logs=log, log=log)
- }
- else {
- if(length(h)>1) {
- if (block == 1){
- ordering <- order(h)
+ if(is.null(copula at calibMoa(normalCopula(0),0)))
+ return(spConCop(dCopula, copula, u, rep(h, length.out=nrow(u)), do.logs=log,
+ log=log))
+
+ if(length(h)>1) {
+ ordering <- order(h)
+
+ # ascending sorted pairs allow for easy evaluation
+ u <- u[ordering,,drop=FALSE]
+ h <- h[ordering]
- # ascending sorted pairs allow for easy evaluation
- u <- u[ordering,,drop=FALSE]
- h <- h[ordering]
+ res <- spDepFunCop(dCopula, copula, u, h, do.logs=log, log=log)
- res <- spDepFunCop(dCopula, copula, u, h, do.logs=log, log=log)
-
- # reordering the values
- res <- res[order(ordering)]
- } else {
- res <- NULL
- for(i in 1:(n%/%block)) {
- res <- c(res, spDepFunCopSnglDist(dCopula, copula,
- u[((i-1)*block+1):(i*block),],
- h[i*block], do.logs=log, log=log))
- }
- }
- } else {
- res <- spDepFunCopSnglDist(dCopula, copula, u, h, do.logs=log, log=log)
- }
- }
-
- return(res)
+ # reordering the values
+ return(res[order(ordering)])
+ } else
+ return(spDepFunCopSnglDist(dCopula, copula, u, h, do.logs=log, log=log))
}
setMethod(dCopula, signature("numeric","spCopula"),
@@ -339,56 +311,31 @@
setMethod(dCopula, signature("matrix","spCopula"), dSpCopula)
## partial derivatives ##
-
## dduSpCopula
###############
-dduSpCopula <- function (u, copula, h, block=1) {
- if (missing(h)) stop("Point pairs need to be provided with their separating distance h.")
-
- n <- nrow(u)
-
- if(length(h)>1 && length(h)!=n) {
+dduSpCopula <- function (u, copula, h) {
+ if (missing(h))
+ stop("Point pairs need to be provided with their separating distance h.")
+ if(length(h)>1 && length(h)!=nrow(u))
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, u, rep(h, length.out=n))
+ return(spConCop(dduCopula, copula, u, rep(h, length.out=nrow(u))))
- else {
- if(length(h)>1) {
- if (block == 1){
- ordering <- order(h)
+ if(length(h)>1) {
+ ordering <- order(h)
- # ascending sorted pairs allow for easy evaluation
- u <- u[ordering,,drop=FALSE]
- h <- h[ordering]
+ # ascending sorted pairs allow for easy evaluation
+ u <- u[ordering, , drop=FALSE]
+ h <- h[ordering]
- res <- spDepFunCop(dduCopula, copula, u, h)
+ res <- spDepFunCop(dduCopula, copula, u, h)
- # reordering the values
- res <- res[order(ordering)]
- } else {
- res <- NULL
- for(i in 1:(n%/%block)) {
- res <- c(res, spDepFunCopSnglDist(dduCopula, copula,
- u[((i-1)*block+1):(i*block),],
- h[i*block]))
- }
- }
- } else {
- res <- spDepFunCopSnglDist(dduCopula, copula, u, h)
- }
- }
-
- if(any(res < 0) || any(res > 1)) {
- warning("Partial derivative produced values outside of [0,1], corrections will be applied to:\n",
- paste(res[res<0],collapse=" "),paste(res[res>1],collapse=" "))
- res[res<0] <- 0
- res[res>1] <- 1
- }
-
- return(res)
+ # reordering the values
+ return(res[order(ordering)])
+ } else
+ return(spDepFunCopSnglDist(dduCopula, copula, u, h))
}
setMethod("dduCopula", signature("matrix","spCopula"), dduSpCopula)
@@ -398,53 +345,28 @@
## ddvSpCopula
###############
-
ddvSpCopula <- function (u, copula, h, block=1) {
- if (missing(h)) stop("Point pairs need to be provided with their separating distance h.")
-
- n <- nrow(u)
-
- if(length(h)>1 && length(h)!=n) {
+ if (missing(h))
+ stop("Point pairs need to be provided with their separating distance h.")
+ if(length(h)>1 && length(h)!=nrow(u))
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, u, rep(h, length.out=n))
+ return(spConCop(ddvCopula, copula, u, rep(h, length.out=nrow(u))))
- else {
- if(length(h)>1) {
- if (block == 1){
- ordering <- order(h)
+ if(length(h)>1) {
+ ordering <- order(h)
- # ascending sorted pairs allow for easy evaluation
- u <- u[ordering,,drop=FALSE]
- h <- h[ordering]
+ # ascending sorted pairs allow for easy evaluation
+ u <- u[ordering,,drop=FALSE]
+ h <- h[ordering]
- res <- spDepFunCop(ddvCopula, copula, u, h)
+ res <- spDepFunCop(ddvCopula, copula, u, h)
- # reordering the values
- res <- res[order(ordering)]
- } else {
- res <- NULL
- for(i in 1:(n%/%block)) {
- res <- c(res, spDepFunCopSnglDist(ddvCopula, copula,
- u[((i-1)*block+1):(i*block),],
- h[i*block]))
- }
- }
- } else {
- res <- spDepFunCopSnglDist(ddvCopula, copula, u, h)
- }
- }
-
- if(any(res < 0) || any(res > 1)) {
- warning("Partial derivative produced values outside of [0,1], corrections will be applied to:\n",
- paste(res[res<0],collapse=" "),paste(res[res>1],collapse=" "))
- res[res<0] <- 0
- res[res>1] <- 1
- }
-
- return(res)
+ # reordering the values
+ return(res[order(ordering)])
+ } else
+ return(spDepFunCopSnglDist(ddvCopula, copula, u, h))
}
setMethod("ddvCopula", signature("matrix","spCopula"), ddvSpCopula)
@@ -485,21 +407,29 @@
# 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) {
+ cor.method=NULL, weighted=FALSE) {
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
+ } 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),]
+ }
+
+ if (weighted) {
+ bins <- as.data.frame(bins[c("np","meanDists","lagCor")])
+ if(!is.na(cutoff))
+ bins <- bins[bins$meanDists <= cutoff,]
+ fitCor <- lm(lagCor ~ poly(meanDists, degree), data = bins, weights=bins$np)
+ } else {
+ bins <- as.data.frame(bins[c("meanDists","lagCor")])
+ if(!is.na(cutoff))
+ bins <- bins[bins$meanDists <= cutoff,]
+ fitCor <- lm(lagCor ~ poly(meanDists, degree), data = bins)
+ }
- fitCor <- lm(lagCor ~ poly(meanDists, degree), data = bins)
-
print(fitCor)
cat("Sum of squared residuals:",sum(fitCor$residuals^2),"\n")
@@ -564,35 +494,6 @@
}
}
-# 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
Modified: pkg/R/spVineCopula.R
===================================================================
--- pkg/R/spVineCopula.R 2013-04-15 11:39:18 UTC (rev 93)
+++ pkg/R/spVineCopula.R 2013-04-23 14:17:11 UTC (rev 94)
@@ -24,13 +24,14 @@
l0 <- rep(0,nrow(u)) # level 0 (spatial) density
u0 <- NULL # level 0 conditional data
- if(!is.matrix(h)) h <- matrix(h, ncol=length(h))
+ if(!is.matrix(h))
+ h <- matrix(h, ncol=length(h))
for(i in 1:(ncol(u)-1)) { # i <- 1
- l0 <- l0+dCopula(as.matrix(u[,c(1,i+1)]), spCop, h=h[,i], log=T)
- u0 <- cbind(u0, dduCopula(as.matrix(u[,c(1,i+1)]), spCop, h=h[,i]))
+ l0 <- l0 + dCopula(u[,c(1,i+1)], spCop, h=h[,i], log=T)
+ u0 <- cbind(u0, dduCopula(u[,c(1,i+1)], spCop, h=h[,i]))
}
-
+
l1 <- dCopula(u0, vine, log=T)
if(log)
return(l0+l1)
@@ -65,6 +66,8 @@
copula=copula at spCop, h=data at distances[,i]))
}
+ cat(summary(as.data.frame(secLevel)))
+
vineCopFit <- fitCopula(copula at vineCop, secLevel, method)
spVineCop <- spVineCopula(copula at spCop, vineCopFit at copula)
@@ -79,22 +82,29 @@
setMethod("fitCopula",signature=signature("spVineCopula"),fitSpVine)
# conditional spatial vine
-condSpVine <- function(condVar, dists, spVine, n=100) {
- rat <- 0.2/(1:(n/2))-(0.1/((n+1)/2))
- xVals <- unique(sort(c(rat,1-rat,1:(n-1)/(n))))
- xLength <- length(xVals)
- repCondVar <- matrix(condVar, ncol=length(condVar), nrow=xLength, byrow=T)
- density <- dCopula(cbind(xVals, repCondVar), spVine, h=dists)
+condSpVine <- function (condVar, dists, spVine, n = 1000) {
+ # add some points in the tails
+ rat <- 29:1%x%c(1e-6,1e-5,1e-4,1e-3)
+ xVals <- unique(sort(c(rat, 1 - rat, 1:(n - 1)/(n))))
+ nx <- length(xVals)
- linAppr <- approxfun(c(0,xVals,1), density[c(1,1:xLength,xLength)] ,yleft=0, yright=0)
- int <- integrate(linAppr,lower=0, upper=1)$value
+ repCondVar <- matrix(condVar, ncol = length(condVar), nrow = nx, byrow = T)
+ density <- dCopula(cbind(xVals, repCondVar), spVine, h = dists)
+ # the 1-e7 corners linearily to [0,1], but keep non-negative
+ density <- c(max(0,2*density[1]-density[2]),
+ density, max(0,2*density[nx]-density[nx-1]))
+ linAppr <- approxfun(c(0, xVals, 1), density)
+
+ # sum up the denstiy to rescale
+ int <- sum(diff(c(0,xVals,1))*(0.5*diff(density)+density[-(nx+2)]))
return(function(u) linAppr(u)/int)
}
# interpolation
-spCopPredict.expectation <- function(predNeigh, spVine, margin) {
+spCopPredict.expectation <- function(predNeigh, spVine, margin, range) {
+ stopifnot(!is.null(range))
stopifnot(is.function(margin$d))
stopifnot(is.function(margin$p))
@@ -106,7 +116,7 @@
condSecVine(margin$p(x))*margin$d(x)*x
}
- predMean <- c(predMean, integrate(condExp,0,3000,subdivisions=1e6)$value)
+ predMean <- c(predMean, integrate(condExp,range[1],range[2],subdivisions=1e6)$value)
}
if ("data" %in% slotNames(predNeigh at locations)) {
res <- predNeigh at locations
@@ -125,10 +135,16 @@
predQuantile <- NULL
for(i in 1:nrow(predNeigh at data)) { # i <-1
condSecVine <- condSpVine(as.numeric(predNeigh at data[i,]), predNeigh at distances[i,], spVine)
- pPred <- optimise(function(x) abs(integrate(condSecVine,0,x)$value-p),
- c(0,1))$minimum
- predQuantile <- c(predQuantile, margin$q(pPred))
+ pPred <- optimise(function(x) abs(integrate(condSecVine, 0, x,
+ subdivisions=10000L,
+ abs.tol=1e-6)$value-p),
+ c(0,1))
+ if(pPred$objective > 1e-6)
+ warning("Numerical evaluation in predQuantile achieved an obkective of only ",
+ pPred$objective, " where 0 has been sought.")
+ predQuantile <- c(predQuantile, margin$q(pPred$minimum))
}
+
if ("data" %in% slotNames(predNeigh at locations)) {
res <- predNeigh at locations
res at data[[paste("quantile.",p,sep="")]] <- predQuantile
@@ -140,8 +156,8 @@
}
}
-spCopPredict <- function(predNeigh, spVine, margin, method="quantile", p=0.5) {
+spCopPredict <- function(predNeigh, spVine, margin, method="quantile", p=0.5, range=NULL) {
switch(method,
quantile=spCopPredict.quantile(predNeigh, spVine, margin, p),
- expectation=spCopPredict.expectation(predNeigh, spVine, margin))
+ expectation=spCopPredict.expectation(predNeigh, spVine, margin, range))
}
\ No newline at end of file
Modified: pkg/R/spatialPreparation.R
===================================================================
--- pkg/R/spatialPreparation.R 2013-04-15 11:39:18 UTC (rev 93)
+++ pkg/R/spatialPreparation.R 2013-04-23 14:17:11 UTC (rev 94)
@@ -1,15 +1,11 @@
-##################################################################
-## ##
-## dedicated functions based on sp preparing the use of copulas ##
-## ##
-##################################################################
+########################################################
+## ##
+## functions based on sp preparing the use of copulas ##
+## ##
+########################################################
-## neighbourhood
-# constructor
-# data = "matrix" a matrix or array providing the data
-# distances="matrix" a matrix providing the distances
-# sp="SpatialPoints" SpatialPoints object providing the coordinates
-# index="matrix" linking the obs. in data to the coordinates
+## neighbourhood constructor
+############################
neighbourhood <- function(data, distances, sp, index, prediction, var){
sizeN <- ncol(distances)+1
@@ -31,7 +27,6 @@
setMethod(show,signature("neighbourhood"),showNeighbourhood)
## names (from sp)
-
setMethod(names, signature("neighbourhood"), namesNeighbourhood <- function(x) x at var)
## spplot ##
@@ -49,12 +44,7 @@
## calculate neighbourhood from SpatialPointsDataFrame
# returns an neighbourhood object
-# spData spatialPointsDataFrame
-# locations the prediction locations, for fitting, locations=spData
-# var one or multiple variable names, all is the default
-# size the size of the neighbourhood, note that for prediction the size
-# is one less than for the copula estimation (default of 5)
-# min.dist the minimum distance between neighbours, must be positive
+##################################
getNeighbours <- function(spData, locations, var=names(spData)[1], size=5,
prediction=FALSE, min.dist=0.01) {
@@ -94,14 +84,6 @@
allLocs, prediction, var))
}
-# data(meuse)
-# coordinates(meuse) <- ~x+y
-# meuseNeigh <- getNeighbours(meuse,var="zinc",size=5)
-# str(meuseNeigh)
-#
-# meuseNeigh <- addAttrToGeom(meuseNeigh at locations,data.frame(rnd=runif(155)),match.ID=F)
-# str(meuseNeigh)
-
#############
## BINNING ##
#############
@@ -131,40 +113,46 @@
}
# the generic calcBins, calculates bins for spatial and spatio-temporal data
+setGeneric("calcBins", function(data, var, nbins=15, boundaries=NA, cutoff=NA,
+ cor.method="kendall", plot=T, ...) {
+ standardGeneric("calcBins")
+ })
-setGeneric("calcBins", function(data, var, nbins=15, boundaries=NA, cutoff=NA, cor.method="kendall", plot=T, ...) standardGeneric("calcBins") )
+## calculating the spatial bins
+################################
-# calculating the spatial bins
-#
-# data denotes the spatial data object
-# var denotes the only variable name used
-# cor.method is passed on to cor() (default="kendall")
-# if plot=TRUE (default), the correlation measures are plotted agaisnt the mean lag separation distance
-#
-calcSpBins <- function(data, var=names(data), nbins=15, boundaries=NA, cutoff=NA, cor.method="kendall", plot=TRUE) {
+calcSpBins <- function(data, var=names(data), nbins=15, boundaries=NA,
+ cutoff=NA, cor.method="kendall", plot=TRUE) {
+ if(is.na(cutoff)) {
+ cutoff <- spDists(coordinates(t(data at bbox)))[1,2]/3
+ }
if(is.na(boundaries)) {
- diagonal <- spDists(coordinates(t(data at bbox)))[1,2]
- boundaries <- ((1:nbins) * min(cutoff,diagonal/3,na.rm=T) / nbins)
+ boundaries <- ((1:nbins) * cutoff/nbins)
}
lags <- calcSpLagInd(data, boundaries)
- mDists <- sapply(lags,function(x) mean(x[,3]))
+ mDists <- sapply(lags, function(x) mean(x[,3]))
+ np <- sapply(lags, function(x) length(x[,3]))
lagData <- lapply(lags, function(x) as.matrix((cbind(data[x[,1],var]@data, data[x[,2],var]@data))))
if(cor.method == "fasttau")
lagCor <- sapply(lagData, function(x) VineCopula:::fasttau(x[,1], x[,2]))
- else
+ if(cor.method %in% c("kendall","spearman","perarson"))
lagCor <- sapply(lagData, function(x) cor(x,method=cor.method)[1,2])
-
+ if(cor.method == "normVariogram")
+ lagCor <- sapply(lagData, function(x) 1-cor(x,method="pearson")[1,2])
+ if(cor.method == "variogram")
+ lagCor <- sapply(lagData, function(x) 0.5*mean((x[,1]-x[,2])^2,na.rm=T))
+
if(plot) {
plot(mDists, lagCor, xlab="distance",ylab=paste("correlation [",cor.method,"]",sep=""),
ylim=1.05*c(-abs(min(lagCor)),max(lagCor)), xlim=c(0,max(mDists)))
abline(h=c(-min(lagCor),0,min(lagCor)),col="grey")
}
- res <- list(meanDists = mDists, lagCor=lagCor, lagData=lagData, lags=lags)
+ res <- list(np=np, meanDists = mDists, lagCor=lagCor, lagData=lagData, lags=lags)
attr(res,"cor.method") <- cor.method
return(res)
}
Modified: pkg/man/condSpVine.Rd
===================================================================
--- pkg/man/condSpVine.Rd 2013-04-15 11:39:18 UTC (rev 93)
+++ pkg/man/condSpVine.Rd 2013-04-23 14:17:11 UTC (rev 94)
@@ -8,7 +8,7 @@
A spatial vine copula is conditioned under the observations of all but one neighbour generating a conditional uivariate distribution used ofr prediction.
}
\usage{
-condSpVine(condVar, dists, spVine, n = 100)
+condSpVine(condVar, dists, spVine, n = 1000)
}
\arguments{
\item{condVar}{
Modified: pkg/man/fitCorFun.Rd
===================================================================
--- pkg/man/fitCorFun.Rd 2013-04-15 11:39:18 UTC (rev 93)
+++ pkg/man/fitCorFun.Rd 2013-04-23 14:17:11 UTC (rev 94)
@@ -8,7 +8,7 @@
Polynomials of different degrees can be fitted to the correlogram calculated using \code{\link{calcBins}}. This function will be used to adjust the copula parameter in the spatial/spatio-temporal copula.
}
\usage{
-fitCorFun(bins, degree = 3, cutoff = NA, bounds = c(0, 1), cor.method = NULL)
+fitCorFun(bins, degree = 3, cutoff = NA, bounds = c(0, 1), cor.method = NULL, weighted = FALSE)
}
\arguments{
\item{bins}{
@@ -26,6 +26,9 @@
\item{cor.method}{
The output of \code{\link{calcBins}} has an attribute \code{cor.method}, in case this is not present the parameter \code{cor.method} will be used. In case the parameter \code{cor.method} is not \code{NULL} and the attribute \code{cor.method} is present, they will be compared.
}
+ \item{weighted}{
+ shall the residulas be weighted by the number of points in the lag class?
+ }
}
\value{
Returns a function that provides correlation estimate for every separating distance.
Modified: pkg/man/spCopPredict.Rd
===================================================================
--- pkg/man/spCopPredict.Rd 2013-04-15 11:39:18 UTC (rev 93)
+++ pkg/man/spCopPredict.Rd 2013-04-23 14:17:11 UTC (rev 94)
@@ -8,7 +8,7 @@
A spatial vine copula is used to predict values at unobserved locations conditioned on observations of a local neighbourhood.
}
\usage{
-spCopPredict(predNeigh, spVine, margin, method = "quantile", p = 0.5)
+spCopPredict(predNeigh, spVine, margin, method = "quantile", p = 0.5, range = NULL)
}
\arguments{
\item{predNeigh}{
@@ -26,6 +26,9 @@
\item{p}{
only used for the quantile predictor indicating the desired fraction the quantile should correspond to.
}
+ \item{range}{
+ the range of integartion to be used in the numerical integration in \code{"expectation"} by \code{\link{integrate}}.
+ }
}
\details{
Predictions are done based on \code{\link{condSpVine}} through numerical integration/optimisation.
Modified: spcopula_0.1-1.tar.gz
===================================================================
(Binary files differ)
Modified: spcopula_0.1-1.zip
===================================================================
(Binary files differ)
More information about the spcopula-commits
mailing list