[spcopula-commits] r93 - / pkg/demo
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Apr 15 13:39:19 CEST 2013
Author: ben_graeler
Date: 2013-04-15 13:39:18 +0200 (Mon, 15 Apr 2013)
New Revision: 93
Added:
pkg/demo/spCopula.R
Removed:
README
pkg/demo/spCopula_estimation.R
Modified:
pkg/demo/
spcopula_0.1-1.tar.gz
spcopula_0.1-1.zip
Log:
- updated demos
Deleted: README
===================================================================
--- README 2013-04-11 15:38:43 UTC (rev 92)
+++ README 2013-04-15 11:39:18 UTC (rev 93)
@@ -1,76 +0,0 @@
- R-Forge SVN README
-
-This file explains the repository structure of your project. A more
-detailed guide to R-Forge is available by
-Theußl and Zeileis (2010) [1] and the R-Forge Administration and
-Development Team (2009) [2].
-
-1. Introduction
------------------------------------------------------------------------
-R is free software distributed under a GNU-style copyleft. R-Forge is
-a central platform for the development of R packages, R-related
-software and further projects. Among many other web-based features it
-provides facilities for collaborative source code management via
-Subversion (SVN) [3].
-
-2. The directory you're in
------------------------------------------------------------------------
-This is the repository of your project. It contains two important
-pre-defined directories namely 'pkg' and 'www'. These directories must
-not be deleted otherwise R-Forge's core functionality will not be
-available (i.e., daily checking and building of your package or the
-project websites).
-'pkg' and 'www' are standardized and therefore are going to be
-described in this README. The rest of your repository can be used as
-you like.
-
-3. 'pkg' directory
------------------------------------------------------------------------
-To make use of the package building and checking feature the package
-source code has to be put into the 'pkg' directory of your repository
-(i.e., 'pkg/DESCRIPTION', 'pkg/R', 'pkg/man', etc.) or, alternatively,
-a subdirectory of 'pkg'. The latter structure allows for having more
-than one package in a single project, e.g., if a project consists of
-the packages foo and bar then the source code will be located in
-'pkg/foo' and 'pkg/bar', respectively.
-
-R-Forge automatically examines the 'pkg' directory of every repository
-and builds the package sources as well as the package binaries on a
-daily basis for Mac OS X and Windows (if applicable). The package builds
-are provided in the 'R Packages' tab for download or can be installed
-directly in R from a CRAN-style repository using
-'install.packages("foo", repos="http://R-Forge.R-project.org")'.
-Furthermore, in the 'R Packages' tab developers can examine logs
-generated on different platforms by the build and check process.
-
-4. 'www' directory
------------------------------------------------------------------------
-Developers may present their project on a subdomain of R-Forge, e.g.,
-'http://foo.R-Forge.R-project.org', or via a link to an external
-website.
-
-This directory contains the project homepage which gets updated hourly
-on R-Forge, so please take into consideration that it will not be
-available right after you commit your changes or additions.
-
-5. Help
------------------------------------------------------------------------
-If you need help don't hesitate to submit a support request at
-https://r-forge.r-project.org/tracker/?func=add&group_id=34&atid=194,
-search the forum
-https://r-forge.r-project.org/forum/forum.php?forum_id=78&group_id=34,
-or contact us at R-Forge at R-project.org.
-
-6. References
------------------------------------------------------------------------
-
-[1] Stefan Theußl and Achim Zeileis. Collaborative software development
-using R-Forge. The R Journal, 1(1):9-14, May 2009. URL
-http://journal.r-project.org/2009-1/RJournal_2009-1_Theussl+Zeileis.pdf
-
-[2] R-Forge Administration and Development Team. RForge User’s Manual,
-2008. URL http://download.R-Forge.R-project.org/R-Forge.pdf
-
-[3] C. M. Pilato, B. Collins-Sussman, and B. W. Fitzpatrick. Version
-Control with Subversion. O’Reilly, 2004. Full book available online at
-http://svnbook.red-bean.com/
Property changes on: pkg/demo
___________________________________________________________________
Added: svn:ignore
+ spCopula_estimation.R
Added: pkg/demo/spCopula.R
===================================================================
--- pkg/demo/spCopula.R (rev 0)
+++ pkg/demo/spCopula.R 2013-04-15 11:39:18 UTC (rev 93)
@@ -0,0 +1,147 @@
+## librarys ##
+library(spcopula)
+library(evd)
+
+## meuse - spatial poionts data.frame ##
+data(meuse)
+coordinates(meuse) = ~x+y
+
+spplot(meuse,"zinc", col.regions=bpy.colors(5))
+
+## margins ##
+hist(meuse[["zinc"]],freq=F,n=30,ylim=c(0,0.0035),
+ main="Histogram of zinc", xlab="zinc concentration")
+
+gevEsti <- fgev(meuse[["zinc"]])$estimate
+meanLog <- mean(log(meuse[["zinc"]]))
+sdLog <- sd(log(meuse[["zinc"]]))
+curve(dgev(x,gevEsti[1], gevEsti[2], gevEsti[3]),add=T,col="red")
+curve(dlnorm(x,meanLog,sdLog),add=T,col="green")
+
+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(meuse,var="zinc",nbins=10,cutoff=800)
+
+# transform data to the unit interval
+bins$lagData <- lapply(bins$lagData, rankTransform)
+
+## calculate parameters for Kendall's tau function ##
+# either linear
+calcKTauLin <- fitCorFun(bins, degree=1, cutoff=600)
+curve(calcKTauLin,0, 1000, col="red",add=TRUE)
+
+# or polynomial (used here)
+calcKTauPol <- fitCorFun(bins, degree=3)
+curve(calcKTauPol,0, 1000, col="purple",add=TRUE)
+
+## find best fitting copula per lag class
+loglikTau <- loglikByCopulasLags(bins, calcKTauPol,
+ families=c(normalCopula(0), tCopula(0),
+ 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
+
+## 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),
+ claytonCopula(0), indepCopula()),
+ distances=bins$meanDists,
+ spDepFun=calcKTauPol, unit="m")
+
+## compare spatial copula loglik by lag:
+spLoglik <- NULL
+for(i in 1:length(bins$lags)) { # i <- 8
+ spLoglik <- c(spLoglik,
+ sum(dCopula(u=bins$lagData[[i]], spCop,log=T,
+ h=bins$lags[[i]][,3])))
+}
+
+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"),
+ pch=c(1,16,5,50), col=c("black", "green", "red"))
+text(x=(1:10+0.5),y=spLoglik,lapply(bins$lagData,length))
+
+##
+# spatial vine
+vineDim <- 5L
+meuseNeigh <- getNeighbours(meuse,var="zinc",size=vineDim)
+meuseNeigh at data <- rankTransform(meuseNeigh at data)
+
+meuseSpVine <- fitCopula(spVineCopula(spCop, vineCopula(as.integer(vineDim-1))),
+ meuseNeigh)
+
+meuseSpVine <- meuseSpVine at copula
+
+##
+# 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-meuse$zinc))
+mean(predMean-meuse$zinc)
+sqrt(mean((predMean-meuse$zinc)^2))
+
+mean(abs(predMedian-meuse$zinc))
+mean(predMedian-meuse$zinc)
+sqrt(mean((predMedian-meuse$zinc)^2))
+
+plot(predMean,meuse$zinc)
+abline(0,1)
+
+plot(predMedian,meuse$zinc)
+abline(0,1)
+
+## kriging results:
+# same neighbourhood size 5L:
+# 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
Deleted: pkg/demo/spCopula_estimation.R
===================================================================
--- pkg/demo/spCopula_estimation.R 2013-04-11 15:38:43 UTC (rev 92)
+++ pkg/demo/spCopula_estimation.R 2013-04-15 11:39:18 UTC (rev 93)
@@ -1,149 +0,0 @@
-## librarys ##
-library(spcopula)
-library(evd)
-
-## meuse - spatial poionts data.frame ##
-data(meuse)
-coordinates(meuse) = ~x+y
-
-spplot(meuse,"zinc", col.regions=bpy.colors(5))
-
-## margins ##
-hist(meuse[["zinc"]],freq=F,n=30,ylim=c(0,0.0035),
- main="Histogram of zinc", xlab="zinc concentration")
-gevEsti <- fgev(meuse[["zinc"]])$estimate
-meanLog <- mean(log(meuse[["zinc"]]))
-sdLog <- sd(log(meuse[["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(meuse[["zinc"]],pgev,gevEsti[1], gevEsti[2], gevEsti[3]) # p: 0.07
-ks.test(meuse[["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(meuse,var="zinc",nbins=10,cutoff=800)
-
-# transform data to the unit interval
-bins$lagData <- lapply(bins$lagData, rankTransform)
-
-## calculate parameters for Kendall's tau function ##
-# either linear
-calcKTauLin <- fitCorFun(bins, degree=1, cutoff=600)
-curve(calcKTauLin,0, 1000, col="red",add=TRUE)
-
-# or polynomial (used here)
-calcKTauPol <- fitCorFun(bins, degree=3)
-curve(calcKTauPol,0, 1000, col="purple",add=TRUE)
-
-## find best fitting copula per lag class
-loglikTau <- loglikByCopulasLags(bins, calcKTauPol,
- families=c(normalCopula(0), tCopula(0),
- 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
-
-## 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),
- claytonCopula(0), indepCopula()),
- distances=bins$meanDists,
- spDepFun=calcKTauPol, unit="m")
-
-## compare spatial copula loglik by lag:
-spLoglik <- NULL
-for(i in 1:length(bins$lags)) { # i <- 8
- spLoglik <- c(spLoglik,
- sum(dCopula(u=bins$lagData[[i]], spCop,log=T,
- h=bins$lags[[i]][,3])))
-}
-
-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"),
- pch=c(1,16,5,50), col=c("black", "green", "red"))
-text(x=(1:10+0.5),y=spLoglik,lapply(bins$lagData,length))
-
-##
-# spatial vine
-vineDim <- 5L
-meuseNeigh <- getNeighbours(meuse,"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-meuse$zinc))
-mean(predMean-meuse$zinc)
-sqrt(mean((predMean-meuse$zinc)^2))
-
-mean(abs(predMedian-meuse$zinc))
-mean(predMedian-meuse$zinc)
-sqrt(mean((predMedian-meuse$zinc)^2))
-
-plot(predMean,meuse$zinc)
-abline(0,1)
-
-plot(predMedian,meuse$zinc)
-abline(0,1)
-
-## kriging results:
-# same neighbourhood size 5L:
-# 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: 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