[spcopula-commits] r87 - in pkg: . R demo man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Feb 20 12:22:14 CET 2013
Author: ben_graeler
Date: 2013-02-20 12:22:13 +0100 (Wed, 20 Feb 2013)
New Revision: 87
Modified:
pkg/DESCRIPTION
pkg/R/spCopula.R
pkg/R/spVineCopula.R
pkg/R/vineCopulas.R
pkg/demo/spCopula_estimation.R
pkg/man/vineCopula.Rd
Log:
- vineCopula class, automated spatial vine copual fitting
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2013-02-08 08:48:54 UTC (rev 86)
+++ pkg/DESCRIPTION 2013-02-20 11:22:13 UTC (rev 87)
@@ -2,7 +2,7 @@
Type: Package
Title: copula driven spatial analysis
Version: 0.1-1
-Date: 2013-02-08
+Date: 2013-02-20
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/spCopula.R
===================================================================
--- pkg/R/spCopula.R 2013-02-08 08:48:54 UTC (rev 86)
+++ pkg/R/spCopula.R 2013-02-20 11:22:13 UTC (rev 87)
@@ -175,7 +175,7 @@
if(do.logs)
res <- log(res)
} else {
- if(class(tmpCop) != "indepCopula")
+ if(class(lowerCop) != "indepCopula")
lowerCop at parameters <- calibPar(lowerCop, h)
res <- fun(pairs, lowerCop, ...)
}
Modified: pkg/R/spVineCopula.R
===================================================================
--- pkg/R/spVineCopula.R 2013-02-08 08:48:54 UTC (rev 86)
+++ pkg/R/spVineCopula.R 2013-02-20 11:22:13 UTC (rev 87)
@@ -3,7 +3,7 @@
#########################################
# constructor
-spVineCopula <- function(spCop, vineCop) {
+spVineCopula <- function(spCop, vineCop=vineCopula()) {
new("spVineCopula", dimension = as.integer(vineCop at dimension+1), parameters=numeric(),
param.names = character(), param.lowbnd = numeric(),
param.upbnd = numeric(), fullname = "Spatial vine copula family.",
@@ -45,4 +45,24 @@
setMethod("dCopula",signature=signature("numeric","spVineCopula"),
function(u, copula, log, ...) {
dspVine(matrix(u,ncol=copula at dimension), copula at spCop, copula at vineCop, log=log, ...)
- })
\ No newline at end of file
+ })
+
+# fiiting the spatial vine for a given spatial copula
+
+fitSpVine <- function(copula, data) {
+ stopifnot(class(data)=="neighbourhood")
+ stopifnot(copula at dimension == ncol(data at data))
+
+ secLevel <- NULL
+ for (i in 1:(copula at dimension-1)) { # i <- 1
+ secLevel <- cbind(secLevel,
+ dduCopula(u=as.matrix(data at data[,c(1,i+1)]),
+ copula=copula at spCop, h=data at distances[,i]))
+ }
+
+ vineCop <- fitCopula(copula at vineCop, secLevel)
+
+ return(spVineCopula(spCop, vineCop))
+}
+
+setMethod("fitCopula",signature=signature("spVineCopula"),fitSpVine)
\ No newline at end of file
Modified: pkg/R/vineCopulas.R
===================================================================
--- pkg/R/vineCopulas.R 2013-02-08 08:48:54 UTC (rev 86)
+++ pkg/R/vineCopulas.R 2013-02-20 11:22:13 UTC (rev 87)
@@ -5,9 +5,19 @@
####################
# constructor
-vineCopula <- function (RVM) {
- if(class(RVM)=="RVineMatrix") # handling non S4-class as subelement in a S4-class
+vineCopula <- function (RVM) { # RVM <- 4L
+ if (is.integer(RVM)) {# assuming dimension; i <- 1
+ Matrix <- NULL
+ for (i in 1:RVM) {
+ Matrix <- cbind(Matrix,c(rep(0,i-1),(RVM-i+1):1))
+ }
+ RVM <- RVineMatrix(Matrix)
+ }
+
+ # handling non S4-class as sub-element in a S4-class
+ stopifnot(class(RVM)=="RVineMatrix")
class(RVM) <- "list"
+
ltr <- lower.tri(RVM$Matrix)
copDef <- cbind(RVM$family[ltr], RVM$par[ltr], RVM$par2[ltr])
copulas <- apply(copDef,1, function(x) copulaFromFamilyIndex(x[1],x[2],x[3]))
@@ -44,7 +54,7 @@
}
setMethod("dCopula", signature("numeric","vineCopula"),
- function(u, copula, ...) dRVine(matrix(u, ncol=copula at dimension), copula, ...))
+ function(u, copula, log, ...) dRVine(matrix(u, ncol=copula at dimension), copula, log, ...))
setMethod("dCopula", signature("matrix","vineCopula"), dRVine)
# ## d-vine structure
@@ -151,7 +161,7 @@
## jcdf ##
pvineCopula <- function(u, copula) {
- empCop <- genEmpCop(copula,1e5)
+ empCop <- genEmpCop(copula, 1e5)
return(pCopula(u, empCop))
}
@@ -188,4 +198,14 @@
RVineSim(n, RVM)
}
-setMethod("rCopula", signature("numeric","vineCopula"), rRVine)
\ No newline at end of file
+setMethod("rCopula", signature("numeric","vineCopula"), rRVine)
+
+# fitting using RVine
+fitVineCop <- function(copula, data, method) {
+ if("StructureSelect" %in% method)
+ vineCopula(RVineStructureSelect(data, indeptest="indeptest" %in% method))
+ else
+ vineCopula(RVineCopSelect(data, Matrix=copula at RVM$Matrix, indeptest="indeptest" %in% method))
+}
+
+setMethod("fitCopula", signature=signature("vineCopula"), fitVineCop)
\ No newline at end of file
Modified: pkg/demo/spCopula_estimation.R
===================================================================
--- pkg/demo/spCopula_estimation.R 2013-02-08 08:48:54 UTC (rev 86)
+++ pkg/demo/spCopula_estimation.R 2013-02-20 11:22:13 UTC (rev 87)
@@ -1,4 +1,5 @@
## librarys ##
+library(spcopula)
library(evd)
## dataset - spatial poionts data.frame ##
@@ -12,17 +13,22 @@
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]
-shape <- gevEsti[3]
-meanLog <- mean(log(meuse[["zinc"]]))
-sdLog <- sd(log(meuse[["zinc"]]))
-curve(dgev(x,loc,scale,shape),add=T,col="red")
+meanLog <- mean(log(dataSet[["zinc"]]))
+sdLog <- sd(log(dataSet[["zinc"]]))
+curve(dgev(x,gevEsti[1], gevEsti[2], gevEsti[3]),add=T,col="red")
curve(dlnorm(x,meanLog,sdLog),add=T,col="green")
-ks.test(dataSet[["zinc"]],pgev,loc,scale,shape) # p: 0.07
+ks.test(dataSet[["zinc"]],pgev,gevEsti[1], gevEsti[2], gevEsti[3]) # p: 0.07
ks.test(dataSet[["zinc"]],plnorm,meanLog,sdLog) # p: 0.03
+pMar <- function(q) plnorm(q, meanLog, sdLog)
+qMar <- function(p) qlnorm(p, meanLog, sdLog)
+dMar <- function(x) dlnorm(x, meanLog, sdLog)
+
+# pMar <- function(q) pgev(q, gevEsti[1], gevEsti[2], gevEsti[3])
+# qMar <- function(p) qgev(p, gevEsti[1], gevEsti[2], gevEsti[3])
+# dMar <- function(x) dgev(x, gevEsti[1], gevEsti[2], gevEsti[3])
+
## lag classes ##
bins <- calcBins(dataSet,var="zinc",nbins=10,cutoff=800)
@@ -44,13 +50,15 @@
claytonCopula(0), frankCopula(1),
gumbelCopula(1), joeBiCopula(1.5),
indepCopula()))
-bestFitTau <- apply(apply(loglikTau, 1, rank, na.last=T), 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),
- frankCopula(1), normalCopula(0), claytonCopula(0),
- claytonCopula(0), claytonCopula(0), claytonCopula(0),
+ frankCopula(1), normalCopula(0),
+ claytonCopula(0), claytonCopula(0),
+ claytonCopula(0), claytonCopula(0),
claytonCopula(0), indepCopula()),
distances=bins$meanDists,
spDepFun=calcKTauPol, unit="m")
@@ -66,6 +74,77 @@
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"),
+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
+text(x=(1:10+0.5),y=spLoglik,lapply(bins$lagData,length))
+
+##
+# spatial vine
+vineDim <- 5L
+meuseNeigh <- getNeighbours(dataSet,"zinc",vineDim)
+meuseNeigh at data <- rankTransform(meuseNeigh at data)
+
+meuseSpVine <- fitCopula(spVineCopula(spCop, vineCopula(as.integer(vineDim-1))),
+ meuseNeigh)
+
+meuseSpVine at vineCop
+
+##
+# leave-one-out x-validation
+
+condVine <- function(condVar, dists, 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), meuseSpVine, h=dists)
+
+ linAppr <- approxfun(c(0,xVals,1), density[c(1,1:xLength,xLength)] ,yleft=0, yright=0)
+ int <- integrate(linAppr,lower=0, upper=1)$value
+
+ return(function(u) linAppr(u)/int)
+}
+
+time <- proc.time() # ~30 s
+predMedian <- NULL
+predMean <- NULL
+for(loc in 1:nrow(meuseNeigh at data)) { # loc <- 429 predNeigh$data[loc,1]
+ cat("Location:",loc,"\n")
+ condSecVine <- condVine(condVar=as.numeric(meuseNeigh at data[loc,-1]),
+ dists=meuseNeigh at distances[loc,,drop=F])
+
+ predMedian <- c(predMedian, qMar(optimise(function(x) abs(integrate(condSecVine,0,x)$value-0.5),c(0,1))$minimum))
+
+ condExp <- function(x) {
+ condSecVine(pMar(x))*dMar(x)*x
+ }
+
+ predMean <- c(predMean, integrate(condExp,0,3000,subdivisions=1e6)$value)
+}
+proc.time()-time
+
+mean(abs(predMean-dataSet$zinc))
+mean(predMean-dataSet$zinc)
+sqrt(mean((predMean-dataSet$zinc)^2))
+
+mean(abs(predMedian-dataSet$zinc))
+mean(predMedian-dataSet$zinc)
+sqrt(mean((predMedian-dataSet$zinc)^2))
+
+plot(predMean,dataSet$zinc)
+abline(0,1)
+
+plot(predMedian,dataSet$zinc)
+abline(0,1)
+
+## kriging results:
+# same neighbourhood size:
+# MAE: 158.61
+# BIAS: -4.24
+# RMSE: 239.85
+#
+# global kriging:
+# MAE: 148.85
+# BIAS: -3.05
+# RMSE: 226.15
\ No newline at end of file
Modified: pkg/man/vineCopula.Rd
===================================================================
--- pkg/man/vineCopula.Rd 2013-02-08 08:48:54 UTC (rev 86)
+++ pkg/man/vineCopula.Rd 2013-02-20 11:22:13 UTC (rev 87)
@@ -11,7 +11,7 @@
}
\arguments{
\item{RVM}{
-An object of class \code{RVineMatrix} generated from \code{\link{RVineMatrix}} in the package \code{\link{VineCopula-package}}.
+An object of class \code{RVineMatrix} generated from \code{\link{RVineMatrix}} in the package \code{\link{VineCopula-package}} or an integer (e.g. \code{4L}) defining the dimension (an independent C-vine of this dimension will be constructed).
}
}
\value{
More information about the spcopula-commits
mailing list