[spcopula-commits] r113 - in pkg: . R demo man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Nov 19 16:49:54 CET 2013
Author: ben_graeler
Date: 2013-11-19 16:49:54 +0100 (Tue, 19 Nov 2013)
New Revision: 113
Added:
pkg/R/spatialGaussianCopula.R
pkg/man/spGaussCopPredict.Rd
pkg/man/spGaussLogLik.Rd
Modified:
pkg/DESCRIPTION
pkg/NAMESPACE
pkg/R/spatialPreparation.R
pkg/demo/spCopula.R
pkg/man/rankTransform.Rd
Log:
- newly implemented spatial Gaussian Copula prediction
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2013-11-19 08:57:27 UTC (rev 112)
+++ pkg/DESCRIPTION 2013-11-19 15:49:54 UTC (rev 113)
@@ -2,7 +2,7 @@
Type: Package
Title: copula driven spatial analysis
Version: 0.1-1
-Date: 2013-10-24
+Date: 2013-11-19
Authors at R: c(person("Benedikt", "Graeler", role = c("aut", "cre"), email =
"ben.graeler at uni-muenster.de"), person("Marius", "Appel",
role = "ctb"))
@@ -37,4 +37,5 @@
stVineCopula.R
utilities.R
returnPeriods.R
+ spatialGaussianCopula.R
zzz.R
Modified: pkg/NAMESPACE
===================================================================
--- pkg/NAMESPACE 2013-11-19 08:57:27 UTC (rev 112)
+++ pkg/NAMESPACE 2013-11-19 15:49:54 UTC (rev 113)
@@ -24,6 +24,7 @@
export(qCopula_u)
export(condSpVine,spCopPredict)
export(condStVine,stCopPredict)
+export(spGaussCopPredict, spGaussLogLik)
# tweaks
# export(setSizeLim)
Added: pkg/R/spatialGaussianCopula.R
===================================================================
--- pkg/R/spatialGaussianCopula.R (rev 0)
+++ pkg/R/spatialGaussianCopula.R 2013-11-19 15:49:54 UTC (rev 113)
@@ -0,0 +1,87 @@
+## spatial Gaussian Copula
+
+# "density" evaluation
+spGaussLogLik <- function(corFun, neigh, log=T) {
+ neighDim <- ncol(neigh at data)
+
+ allDataDists <- spDists(neigh at dataLocs)
+
+ pb <- txtProgressBar(0, nrow(neigh at data), 0, width = getOption("width") - 10, style = 3)
+
+ loglik <- 0
+
+ for(i in 1:nrow(neigh at data)) { # i <- 2
+ setTxtProgressBar(pb, i)
+ tmpDists <- allDataDists[neigh at index[i,], neigh at index[i,]]
+
+ tmpCor <- corFun(tmpDists)
+
+ tmpGaussCop <- normalCopula(tmpCor[lower.tri(tmpCor)], neighDim, dispstr="un")
+
+ loglik <- loglik + dCopula(as.numeric(neigh at data[i,]), tmpGaussCop, log=T)
+ }
+ close(pb)
+
+ if(log)
+ return(loglik)
+ else
+ return(exp(loglik))
+}
+
+# interpolation based on a valid corelogram function
+spGaussCopPredict <- function(corFun, predNeigh, margin, p=0.5, ..., n=1000) {
+ stopifnot(is.list(margin))
+ stopifnot(is(margin$q, "function"))
+
+ neighDim <- ncol(predNeigh at data)
+ allDataDists <- spDists(predNeigh at dataLocs)
+
+ pb <- txtProgressBar(0, nrow(predNeigh at data), 0, width = getOption("width") - 10, style = 3)
+
+ predQuantile <- NULL
+
+ for(i in 1:nrow(predNeigh at data)) { # i <- 2
+ setTxtProgressBar(pb, i)
+ tmpDataDists <- allDataDists[predNeigh at index[i,-1], predNeigh at index[i,-1]]
+
+ tmpDists <- rbind(c(0,predNeigh at distances[i,]),
+ cbind(predNeigh at distances[i,], tmpDataDists))
+ tmpCor <- corFun(tmpDists)
+
+ tmpGaussCop <- normalCopula(tmpCor[lower.tri(tmpCor)], neighDim, dispstr="un")
+ rat <- 50:1 %x% c(1e-06, 1e-05, 1e-04, 0.001)
+ xVals <- unique(sort(c(rat, 1 - rat, 1:(n - 1)/n)))
+ nx <- length(xVals)
+
+ condGausCop <- dCopula(cbind(xVals, matrix(rep(as.numeric(predNeigh at data[i,-1]),
+ length(xVals)), ncol=neighDim-1, byrow=T)),
+ tmpGaussCop)
+
+ condGausCop <- c(max(0, 2 * condGausCop[1] - condGausCop[2]), condGausCop,
+ max(0, 2 * condGausCop[nx] - condGausCop[nx - 1]))
+ int <- sum(diff(c(0, xVals, 1)) * (0.5 * diff(condGausCop) + condGausCop[-(nx + 2)]))
+
+ condVineFun <- approxfun(c(0, xVals, 1), condGausCop/int)
+
+ condGausCop <- condVineFun(xVals)
+ int <- cumsum(c(0, diff(xVals) * (0.5 * diff(condGausCop) + condGausCop[-nx])))
+ lower <- max(which(int <= p))
+ m <- (condGausCop[lower + 1] - condGausCop[lower])/(xVals[lower + 1] - xVals[lower])
+ b <- condGausCop[lower]
+ xRes <- -b/m + sign(m) * sqrt(b^2/m^2 + 2 * (p - int[lower])/m)
+ predQuantile <- c(predQuantile, margin$q(xVals[lower] + xRes))
+ }
+ close(pb)
+
+ if ("data" %in% slotNames(predNeigh at predLocs)) {
+ res <- predNeigh at predLocs
+ res at data[[paste("quantile.", p, sep = "")]] <- predQuantile
+ return(res)
+ }
+ else {
+ predQuantile <- data.frame(predQuantile)
+ colnames(predQuantile) <- paste("quantile.", p, sep = "")
+ return(addAttrToGeom(predNeigh at predLocs, predQuantile,
+ match.ID = FALSE))
+ }
+}
Modified: pkg/R/spatialPreparation.R
===================================================================
--- pkg/R/spatialPreparation.R 2013-11-19 08:57:27 UTC (rev 112)
+++ pkg/R/spatialPreparation.R 2013-11-19 15:49:54 UTC (rev 113)
@@ -209,7 +209,7 @@
if(cor.method == "fasttau")
lagCor <- sapply(lagData, function(x) VineCopula:::fasttau(x[,1], x[,2]))
- if(cor.method %in% c("kendall","spearman","perarson"))
+ if(cor.method %in% c("kendall","spearman","pearson"))
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])
Modified: pkg/demo/spCopula.R
===================================================================
--- pkg/demo/spCopula.R 2013-11-19 08:57:27 UTC (rev 112)
+++ pkg/demo/spCopula.R 2013-11-19 15:49:54 UTC (rev 113)
@@ -85,11 +85,13 @@
meuseSpVine <- fitCopula(spVineCopula(spCop, vineCopula(as.integer(vineDim-1))),
meuseNeigh)
+# log-likelihood:
+meuseSpVine at loglik
+
meuseSpVine <- meuseSpVine at copula
##
# leave-one-out x-validation
-
time <- proc.time() # ~100 s
predMedian <- NULL
predMean <- NULL
Modified: pkg/man/rankTransform.Rd
===================================================================
--- pkg/man/rankTransform.Rd 2013-11-19 08:57:27 UTC (rev 112)
+++ pkg/man/rankTransform.Rd 2013-11-19 15:49:54 UTC (rev 113)
@@ -7,19 +7,14 @@
performs the rank order transformational
}
\usage{
-rankTransform(u, v = NULL, ties.method = "average")
+rankTransform(u, v = NULL, na.last=TRUE, ties.method = "average")
}
\arguments{
- \item{u}{
-a matrix or data.frame with at least two columns holding the data or a vector holding the first column of the data
+ \item{u}{a matrix or data.frame with at least two columns holding the data or a vector holding the first column of the data}
+ \item{v}{a vector holding the second column of the data}
+ \item{na.last}{puts \code{NA}s last by default}
+ \item{ties.method}{How should ties be treated in \code{\link{rank}}?}
}
- \item{v}{
-a vector holding the second column of the data
-}
- \item{ties.method}{
-How should ties be treated in \code{\link{rank}}?
-}
-}
\value{
A matrix or data.frame (as provided) with the transformed data. Rows containing any \code{NA} will be dropped.
}
Added: pkg/man/spGaussCopPredict.Rd
===================================================================
--- pkg/man/spGaussCopPredict.Rd (rev 0)
+++ pkg/man/spGaussCopPredict.Rd 2013-11-19 15:49:54 UTC (rev 113)
@@ -0,0 +1,92 @@
+\name{spGaussCopPredict}
+\alias{spGaussCopPredict}
+\title{
+spatial prediction using a Gaussian Copula
+}
+\description{
+This function allows to do spatial prediction using a Gaussian Copula based on the spatial parametrization.
+}
+\usage{
+spGaussCopPredict(corFun, predNeigh, margin, p = 0.5, ..., n = 1000)
+}
+
+\arguments{
+ \item{corFun}{
+A valid correlogram (i.e. producing a valid correlation matrix; e.g. based on a variogram).
+}
+ \item{predNeigh}{
+A \code{\linkS4class{neighbourhood}} object used for prediction.
+}
+ \item{margin}{
+a list containing the marginal distribution. Currently only the entry \code{q} is required defining the quantile function.
+}
+ \item{p}{
+the fraction the quantile function shall be evaluated for. This can be used to calculate besides the median estimate confidence estimates as well.
+}
+ \item{\dots}{
+currently unused
+}
+ \item{n}{
+the approximate number of points used in the linear approximation of the conditional distribution function.
+}
+}
+\details{
+Based on \code{corFun} provided with a distance matrix a Gaussian copula (\code{\linkS4class{normalCopula}}) is generated and conditioned under the data of the neighbouring locations. The 1-dimensional conditional distribution is approximated and evaluated for the given fraction \code{p}. This conditioned fraction is than passed on to the marginal quantile function and evaluated providing an estimate.
+}
+\value{
+According to the \code{predLocs} slot of the provided \code{\linkS4class{neighbourhood}} a spatial data structure extended with an variable holding the predicted values.
+}
+\author{
+Benedikt Graeler
+}
+
+
+\seealso{
+\code{\link{spCopPredict}}
+}
+\examples{
+# load data from the Meuse demo
+data(spCopDemo)
+
+# calculate the correlation function based on Kendall's tau
+calcKTauPol <- fitCorFun(bins, degree=1)
+
+# translate Kendall's tau correlation function into Gaussian Copula parameters
+# using a linear variogram
+meuseGaussCorFun <- function(h) {
+ res <- pmax(iTau(normalCopula(0), calcKTauPol(0))/658*(658-h),0)
+ res[h ==0] <- 1
+ return(res)
+}
+
+# get some prediction data
+library(sp)
+data(meuse.grid)
+coordinates(meuse.grid) <- ~x+y
+gridded(meuse.grid) <- TRUE
+
+data(meuse)
+coordinates(meuse) <- ~x+y
+
+meuse$rtZinc <- rank(meuse$zinc)/(length(meuse)+1)
+
+# obtain the prediction neighbourhoods
+predMeuseNeigh <- getNeighbours(meuse[1:4,], meuse.grid[c(9:12,16:19,25:28),],
+ "rtZinc", 5L, TRUE, -1)
+
+qMar <- function(x) {
+ qlnorm(x,mean(log(meuse$zinc)),sd(log(meuse$zinc)))
+}
+
+# predict using the Gaussian Copula
+predMedian <- spGaussCopPredict(meuseGaussCorFun, predMeuseNeigh, list(q=qMar))
+
+\dontrun{
+ spplot(predMedian,"quantile.0.5",
+ sp.layout=list("sp.points", meuse, pch = 19, col = "red"),
+ col.regions=bpy.colors())
+}
+}
+
+\keyword{ prediction }
+\keyword{ distribution }
Added: pkg/man/spGaussLogLik.Rd
===================================================================
--- pkg/man/spGaussLogLik.Rd (rev 0)
+++ pkg/man/spGaussLogLik.Rd 2013-11-19 15:49:54 UTC (rev 113)
@@ -0,0 +1,66 @@
+\name{spGaussLogLik}
+\alias{spGaussLogLik}
+
+\title{
+Density evalaution for a spatial Gaussian Copula
+}
+\description{
+Evaluates the density for a spatial Gaussian Copula.
+}
+\usage{
+spGaussLogLik(corFun, neigh, log = T)
+}
+
+\arguments{
+ \item{corFun}{
+A valid correlogram (i.e. producing a valid correlation matrix; e.g. based on a variogram).
+}
+ \item{neigh}{
+A \code{\linkS4class{neighbourhood}} object to be evaluated.
+}
+ \item{log}{
+Should the log-likelihood be returned?
+}
+}
+\details{
+Evaluates the density for all neioghbourhoods in \code{neigh} and returns the (log)-likelihood.
+}
+\value{
+The (log)-likelihood value.
+}
+\author{
+Benedikt Graeler
+}
+
+\examples{
+library(spcopula)
+
+# load data from the Meuse demo
+data(spCopDemo)
+
+# calculate the correlation function based on Kendall's tau
+calcKTauPol <- fitCorFun(bins, degree=1)
+
+# translate Kendall's tau correlation function into Gaussian Copula parameters
+# using a linear variogram
+meuseGaussCorFun <- function(h) {
+ res <- pmax(iTau(normalCopula(0), calcKTauPol(0))/658*(658-h),0)
+ res[h ==0] <- 1
+ return(res)
+}
+
+# get the neighbours
+library(sp)
+data(meuse)
+coordinates(meuse) <- ~x+y
+
+meuse$rtZinc <- rank(meuse$zinc)/(length(meuse)+1)
+
+meuseNeigh <- getNeighbours(meuse, var="rtZinc", size=5L,
+ prediction=FALSE)
+
+# calculate the log-likelihood
+spGaussLogLik(meuseGaussCorFun, meuseNeigh)
+}
+
+\keyword{ distribution }
More information about the spcopula-commits
mailing list