[spcopula-commits] r137 - in pkg: . R demo
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Sep 19 13:15:54 CEST 2014
Author: ben_graeler
Date: 2014-09-19 13:15:54 +0200 (Fri, 19 Sep 2014)
New Revision: 137
Modified:
pkg/DESCRIPTION
pkg/R/spCopula.R
pkg/R/spVineCopula.R
pkg/R/spatialPreparation.R
pkg/demo/pureSpVineCopula.R
Log:
- updated demo on pure spatial vine copulas
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2014-07-16 10:00:26 UTC (rev 136)
+++ pkg/DESCRIPTION 2014-09-19 11:15:54 UTC (rev 137)
@@ -1,8 +1,8 @@
Package: spcopula
Type: Package
Title: copula driven spatial analysis
-Version: 0.2-0
-Date: 2014-07-16
+Version: 0.2-1
+Date: 2014-09-19
Authors at R: c(person("Benedikt", "Graeler", role = c("aut", "cre"),
email = "ben.graeler at uni-muenster.de"),
person("Marius", "Appel",role = "ctb"))
Modified: pkg/R/spCopula.R
===================================================================
--- pkg/R/spCopula.R 2014-07-16 10:00:26 UTC (rev 136)
+++ pkg/R/spCopula.R 2014-09-19 11:15:54 UTC (rev 137)
@@ -638,11 +638,16 @@
calcCor, lagSub=1:length(bins$meanDists)) {
var <- attr(bins, "variable")
- lagData <- lapply(bins$lags[lagSub],
- function(x) {
- as.matrix((cbind(data[x[,1], var]@data,
- data[x[,2], var]@data)))
- })
+ if(missing(data)) {
+ lagData <- bins$lagData
+ }
+ else {
+ lagData <- lapply(bins$lags[lagSub],
+ function(x) {
+ as.matrix((cbind(data[x[, 1], var]@data,
+ data[x[, 2], var]@data)))
+ })
+ }
lagData <- lapply(lagData,
function(pairs) {
Modified: pkg/R/spVineCopula.R
===================================================================
--- pkg/R/spVineCopula.R 2014-07-16 10:00:26 UTC (rev 136)
+++ pkg/R/spVineCopula.R 2014-09-19 11:15:54 UTC (rev 137)
@@ -214,7 +214,7 @@
close(pb)
if ("data" %in% slotNames(predLocs)) {
- res <- predNeigh at predLocs
+ res <- predLocs
res at data[["expect"]] <- predMean
return(res)
} else {
Modified: pkg/R/spatialPreparation.R
===================================================================
--- pkg/R/spatialPreparation.R 2014-07-16 10:00:26 UTC (rev 136)
+++ pkg/R/spatialPreparation.R 2014-09-19 11:15:54 UTC (rev 137)
@@ -204,6 +204,8 @@
data <- as.matrix(data at data)
+ lagData <- list()
+
for (i in 1:nbins) {
bools <- (dists <= boundaries[i+1] & dists > boundaries[i])
@@ -212,7 +214,7 @@
pairs <- rbind(pairs, data[bools[,col],c(1,1+col)])
}
-# lagData <- append(lagData, list(pairs))
+ lagData[[i]] <- pairs
moa[i] <- corFun(pairs)
meanDists[i] <- mean(dists[bools])
np[i] <- sum(bools)
@@ -224,8 +226,7 @@
abline(h=c(-min(moa),0,min(moa)),col="grey")
}
-# res <- list(np=np, meanDists = meanDists, lagCor=moa, lagData=lagData)
- res <- list(np=np, meanDists = meanDists, lagCor=moa)
+ res <- list(np=np, meanDists = meanDists, lagCor=moa, lagData=lagData)
attr(res,"cor.method") <- switch(cor.method, fasttau="kendall", cor.method)
attr(res,"variable") <- var
Modified: pkg/demo/pureSpVineCopula.R
===================================================================
--- pkg/demo/pureSpVineCopula.R 2014-07-16 10:00:26 UTC (rev 136)
+++ pkg/demo/pureSpVineCopula.R 2014-09-19 11:15:54 UTC (rev 137)
@@ -2,7 +2,7 @@
library("spcopula")
# library("evd")
library("sp")
-
+par(mfrow=c(1,1))
## meuse - spatial poionts data.frame ##
data("meuse")
coordinates(meuse) = ~x+y
@@ -41,7 +41,7 @@
surClaytonCopula(1), surGumbelCopula(1), surJoeBiCopula(1.5))
## find best fitting copula per lag class
-loglikTau <- loglikByCopulasLags(bins, families, calcKTau)
+loglikTau <- loglikByCopulasLags(bins, meuse, families, calcKTau)
bestFitTau <- apply(apply(loglikTau$loglik, 1, rank, na.last=T), 2,
function(x) which(x==9))
colnames(loglikTau$loglik)[bestFitTau]
@@ -59,16 +59,15 @@
## second spatial tree
#######################
-meuseNeigh2 <- dropSpTree(meuseNeigh1, spCop)
-bins2 <- calcBins(meuseNeigh2, boundaries=c(0,2:15)*50, plot=F)
+meuseNeigh2 <- dropSpTree(meuseNeigh1, meuse, spCop)
+bins2 <- calcBins(meuseNeigh2, "rtZinc", boundaries=c(0,2:15)*50, plot=F)
points(bins2$meanDists, bins2$lagCor, pch=2)
calcKTau2 <- fitCorFun(bins2, degree=3,cutoff=500)
curve(calcKTau2,0, 800, col="green",add=TRUE)
-loglikTau2 <- loglikByCopulasLags(bins2, families, calcKTau2)
+loglikTau2 <- loglikByCopulasLags(bins2, families = families, calcCor = calcKTau2)
bestFitTau2 <- apply(apply(loglikTau2$loglik, 1, rank, na.last=T), 2,
function(x) which(x==9))
-
colnames(loglikTau2$loglik)[bestFitTau2]
## set up the second bivariate spatial Copula
@@ -79,13 +78,13 @@
## third spatial tree
######################
-meuseNeigh3 <- dropSpTree(meuseNeigh2, spCop2)
-bins3 <- calcBins(meuseNeigh3, plot=F)
+meuseNeigh3 <- dropSpTree(meuseNeigh2, meuse, spCop2)
+bins3 <- calcBins(meuseNeigh3, "rtZinc", plot=F)
points(bins3$meanDists, bins3$lagCor, pch=3)
-calcKTau3 <- fitCorFun(bins3, degree=1, cutoff=400)
+calcKTau3 <- fitCorFun(bins3, degree=1, cutoff=500)
curve(calcKTau3, 0, 500, col="red", add=TRUE)
-loglikTau3 <- loglikByCopulasLags(bins3, families, calcKTau3)
+loglikTau3 <- loglikByCopulasLags(bins3, families = families, calcCor = calcKTau3)
bestFitTau3 <- apply(apply(loglikTau3$loglik, 1, rank, na.last=T), 2,
function(x) which(x==9))
colnames(loglikTau3$loglik)[bestFitTau3]
@@ -98,8 +97,8 @@
## fourth spatial tree
#######################
-meuseNeigh4 <- dropSpTree(meuseNeigh3, spCop3)
-bins4 <- calcBins(meuseNeigh4, plot=F)
+meuseNeigh4 <- dropSpTree(meuseNeigh3, meuse, spCop3)
+bins4 <- calcBins(meuseNeigh4, "rtZinc", plot=F)
points(bins4$meanDists, bins4$lagCor, pch=4)
calcKTau4 <- fitCorFun(bins4, degree=1,cutoff=400)
curve(calcKTau4,0, 500, col="blue",add=TRUE)
@@ -108,7 +107,7 @@
"3. spatial cop.", "4. spatial cop."),
pch=1:4,col=c("purple","green","red","blue"),lty=1)
-loglikTau4 <- loglikByCopulasLags(bins4, families, calcKTau4)
+loglikTau4 <- loglikByCopulasLags(bins4, families = families,calcCor = calcKTau4)
bestFitTau4 <- apply(apply(loglikTau4$loglik, 1, rank, na.last=T), 2,
function(x) which(x==9))
colnames(loglikTau4$loglik)[bestFitTau4]
@@ -135,24 +134,21 @@
# meuseNeigh <- getNeighbours(dataLocs=meuse, predLocs=meuse, prediction=T,
# min.dist=10, var="evZinc", size=vineDim)
-# meuseNeigh <- getNeighbours(dataLocs=meuse, predLocs=meuse, prediction=T,
-# min.dist=10, var="rtZinc", size=vineDim)
-
## leave-one-out x-validation
##############################
time <- proc.time() # ~160 s
-predMedian <- spCopPredict(meuseNeigh, meuseSpVine, list(q=qMar), "quantile")
-predMean <- spCopPredict(meuseNeigh, meuseSpVine, list(q=qMar), "expectation")
+predMedian <- spCopPredict(meuseNeigh, meuse, meuse, meuseSpVine, margin=list(q=qMar), "quantile")
+predMean <- spCopPredict(meuseNeigh, meuse, meuse, meuseSpVine, margin=list(q=qMar), "expectation")
proc.time() - time
c(mean(abs(predMean$expect-meuse$zinc)),
mean(predMean$expect-meuse$zinc),
sqrt(mean((predMean$expect-meuse$zinc)^2)))
-c(mean(abs(predMedian$quantile.0.5-meuse$zinc)),
- mean(predMedian$quantile.0.5-meuse$zinc),
- sqrt(mean((predMedian$quantile.0.5-meuse$zinc)^2)))
+c(mean(abs(predMedian$quantile.0.5-meuse$zinc),na.rm = T),
+ mean(predMedian$quantile.0.5-meuse$zinc, na.rm = T),
+ sqrt(mean((predMedian$quantile.0.5-meuse$zinc)^2, na.rm = T)))
par(mfrow=c(1,2))
plot(predMean$expect, meuse$zinc,asp=1)
@@ -161,6 +157,9 @@
plot(predMedian$quantile.0.5, meuse$zinc,asp=1)
abline(0,1)
+boxplot(predMean at data[c("zinc","expect")])
+boxplot(predMedian at data[c("zinc","quantile.0.5")])
+
## kriging results:
# same neighbourhood size 5L:
# MAE: 158.61
More information about the spcopula-commits
mailing list