[spcopula-commits] r101 - / pkg/R thoughts-Ben
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Aug 1 12:31:15 CEST 2013
Author: ben_graeler
Date: 2013-08-01 12:31:15 +0200 (Thu, 01 Aug 2013)
New Revision: 101
Added:
thoughts-Ben/
thoughts-Ben/Meuse_spcopula_estimation.R
thoughts-Ben/empiricalCopula_COPULA.R
thoughts-Ben/tawn3pCopula-stub.R
Modified:
pkg/R/BB1copula.R
pkg/R/BB6copula.R
pkg/R/BB7copula.R
pkg/R/BB8copula.R
pkg/R/ClaytonGumbelCopula.R
pkg/R/joeBiCopula.R
pkg/R/linkingVineCopula.R
pkg/R/partialDerivatives.R
spcopula_0.1-1.tar.gz
spcopula_0.1-1.zip
Log:
- some textual updates
- added thought folder
Modified: pkg/R/BB1copula.R
===================================================================
--- pkg/R/BB1copula.R 2013-07-23 11:46:51 UTC (rev 100)
+++ pkg/R/BB1copula.R 2013-08-01 10:31:15 UTC (rev 101)
@@ -30,7 +30,7 @@
stop(paste("Parameter values out of bounds: theta: (0,Inf), delta: [1,Inf)."))
new("BB1Copula", dimension = as.integer(2), parameters = param,
param.names = c("theta", "delta"), param.lowbnd = c(0, 1), param.upbnd = c(Inf, Inf),
- family=7, fullname = "BB1 copula family. Number 7 in CDVine.")
+ family=7, fullname = "BB1 copula family. Number 7 in VineCopula.")
}
## density ##
@@ -101,7 +101,7 @@
stop(paste("Parameter values out of bounds: theta: (0,Inf), delta: [1,Inf)."))
new("surBB1Copula", dimension = as.integer(2), parameters = param,
param.names = c("theta", "delta"), param.lowbnd = c(0, 1), param.upbnd = c(Inf, Inf),
- family=17, fullname = "Survival BB1 copula family. Number 17 in CDVine.")
+ family=17, fullname = "Survival BB1 copula family. Number 17 in VineCopula.")
}
## density ##
@@ -168,7 +168,7 @@
stop(paste("Parameter values out of bounds: theta: (-Inf,0), delta: (-Inf,-1]."))
new("r90BB1Copula", dimension = as.integer(2), parameters = param,
param.names = c("theta", "delta"), param.lowbnd = c(-Inf, -Inf), param.upbnd = c(0, -1),
- family=27, fullname = "90 deg rotated BB1 copula family. Number 27 in CDVine.")
+ family=27, fullname = "90 deg rotated BB1 copula family. Number 27 in VineCopula.")
}
BiCopCDF
## density ##
@@ -222,7 +222,7 @@
stop(paste("Parameter values out of bounds: theta: (-Inf,0), delta: (-Inf,-1]."))
new("r270BB1Copula", dimension = as.integer(2), parameters = param,
param.names = c("theta", "delta"), param.lowbnd = c(-Inf, -Inf), param.upbnd = c(0, -1),
- family=37, fullname = "270 deg rotated BB1 copula family. Number 37 in CDVine.")
+ family=37, fullname = "270 deg rotated BB1 copula family. Number 37 in VineCopula.")
}
## density ##
Modified: pkg/R/BB6copula.R
===================================================================
--- pkg/R/BB6copula.R 2013-07-23 11:46:51 UTC (rev 100)
+++ pkg/R/BB6copula.R 2013-08-01 10:31:15 UTC (rev 101)
@@ -32,7 +32,7 @@
stop("Parameter value(s) out of bound(s): theta: [1,Inf), delta: [1,Inf).")
new("BB6Copula", dimension = as.integer(2), parameters = param,
param.names = c("theta", "delta"), param.lowbnd = c(1, 1), param.upbnd = c(Inf, Inf),
- family=8, fullname = "BB6 copula family. Number 8 in CDVine.")
+ family=8, fullname = "BB6 copula family. Number 8 in VineCopula.")
}
## density ##
@@ -67,7 +67,7 @@
## random number generater ??
setMethod("rCopula", signature("numeric","BB6Copula"), linkVineCop.r)
-## kendall distribution/measure, taken from CDVine:::obs.stat
+## kendall distribution/measure, taken from VineCopula:::obs.stat
kendall.BB6 <- function(copula, t){
theta = copula at parameters[1]
delta = copula at parameters[2]
@@ -103,7 +103,7 @@
stop("Parameter value(s) out of bound(s): theta: [1,Inf), delta: [1,Inf).")
new("surBB6Copula", dimension = as.integer(2), parameters = param,
param.names = c("theta", "delta"), param.lowbnd = c(1, 1), param.upbnd = c(Inf, Inf),
- family=18, fullname = "Survival BB6 copula family. Number 18 in CDVine.")
+ family=18, fullname = "Survival BB6 copula family. Number 18 in VineCopula.")
}
## density ##
@@ -170,7 +170,7 @@
stop("Parameter value out of bound: theta: (-Inf,1], delta: (-Inf,1].")
new("r90BB6Copula", dimension = as.integer(2), parameters = param,
param.names = c("theta", "delta"), param.lowbnd = c(-Inf, -Inf), param.upbnd = c(-1, -1),
- family=28, fullname = "90 deg rotated BB6 copula family. Number 28 in CDVine.")
+ family=28, fullname = "90 deg rotated BB6 copula family. Number 28 in VineCopula.")
}
## density ##
@@ -224,7 +224,7 @@
stop("Parameter value out of bound: theta: (-Inf,1], delta: (-Inf,1].")
new("r270BB6Copula", dimension = as.integer(2), parameters = param,
param.names = c("theta", "delta"), param.lowbnd = c(-Inf, -Inf), param.upbnd = c(-1, -1),
- family=38, fullname = "270 deg rotated BB6 copula family. Number 38 in CDVine.")
+ family=38, fullname = "270 deg rotated BB6 copula family. Number 38 in VineCopula.")
}
## density ##
Modified: pkg/R/BB7copula.R
===================================================================
--- pkg/R/BB7copula.R 2013-07-23 11:46:51 UTC (rev 100)
+++ pkg/R/BB7copula.R 2013-08-01 10:31:15 UTC (rev 101)
@@ -30,7 +30,7 @@
stop(paste("Parameter values out of bounds: theta: [1,Inf), delta: (0,Inf)."))
new("BB7Copula", dimension = as.integer(2), parameters = param,
param.names = c("theta", "delta"), param.lowbnd = c(1, 0), param.upbnd = c(Inf, Inf),
- family=9, fullname = "BB7 copula family. Number 9 in CDVine.")
+ family=9, fullname = "BB7 copula family. Number 9 in VineCopula.")
}
## density ##
@@ -65,7 +65,7 @@
## random number generator
setMethod("rCopula", signature("numeric","BB7Copula"), linkVineCop.r)
-## kendall distribution/measure, taken from CDVine:::obs.stat
+## kendall distribution/measure, taken from VineCopula:::obs.stat
kendall.BB7 <- function(copula, t){
theta = copula at parameters[1]
delta = copula at parameters[2]
@@ -103,7 +103,7 @@
stop(paste("Parameter values out of bounds: theta: [1,Inf), delta: (0,Inf)."))
new("surBB7Copula", dimension = as.integer(2), parameters = param,
param.names = c("theta", "delta"), param.lowbnd = c(1, 0), param.upbnd = c(Inf, Inf),
- family= 19, fullname = "Survival BB7 copula family. Number 19 in CDVine.")
+ family= 19, fullname = "Survival BB7 copula family. Number 19 in VineCopula.")
}
## density ##
@@ -172,7 +172,7 @@
stop(paste("Parameter values out of bounds: theta: (-Inf,-1], delta: (-Inf,0)."))
new("r90BB7Copula", dimension = as.integer(2), parameters = param,
param.names = c("theta", "delta"), param.lowbnd = c(-Inf, -Inf), param.upbnd = c(-1, 0),
- family=29, fullname = "90 deg rotated BB7 copula family. Number 29 in CDVine.")
+ family=29, fullname = "90 deg rotated BB7 copula family. Number 29 in VineCopula.")
}
## density ##
@@ -226,7 +226,7 @@
stop(paste("Parameter values out of bounds: theta: (-Inf,-1], delta: (-Inf,0)."))
new("r270BB7Copula", dimension = as.integer(2), parameters = param,
param.names = c("theta", "delta"), param.lowbnd = c(-Inf, -Inf), param.upbnd = c(-1, -1),
- family=39, fullname = "270 deg rotated BB7 copula family. Number 39 in CDVine.")
+ family=39, fullname = "270 deg rotated BB7 copula family. Number 39 in VineCopula.")
}
## density ##
Modified: pkg/R/BB8copula.R
===================================================================
--- pkg/R/BB8copula.R 2013-07-23 11:46:51 UTC (rev 100)
+++ pkg/R/BB8copula.R 2013-08-01 10:31:15 UTC (rev 101)
@@ -32,7 +32,7 @@
stop("Parameter value out of bound: theta: [1,Inf), delta: (0,1].")
new("BB8Copula", dimension = as.integer(2), parameters = param,
param.names = c("theta", "delta"), param.lowbnd = c(1, 0), param.upbnd = c(Inf, 1),
- family=10, fullname = "BB8 copula family. Number 10 in CDVine.")
+ family=10, fullname = "BB8 copula family. Number 10 in VineCopula.")
}
## density ##
@@ -67,7 +67,7 @@
## random number generator
setMethod("rCopula", signature("numeric","BB8Copula"),linkVineCop.r)
-## kendall distribution/measure, taken from CDVine:::obs.stat
+## kendall distribution/measure, taken from VineCopula:::obs.stat
kendall.BB8 <- function(copula, t){
theta = copula at parameters[1]
delta = copula at parameters[2]
@@ -103,7 +103,7 @@
stop("Parameter value out of bound: theta: [1,Inf), delta: (0,1].")
new("surBB8Copula", dimension = as.integer(2), parameters = param,
param.names = c("theta", "delta"), param.lowbnd = c(1, 0), param.upbnd = c(Inf, 1),
- family=20, fullname = "Survival BB8 copula family. Number 20 in CDVine.")
+ family=20, fullname = "Survival BB8 copula family. Number 20 in VineCopula.")
}
## density ##
@@ -170,7 +170,7 @@
stop("Parameter value out of bound: theta: (-Inf,-1], delta: [-1,0).")
new("r90BB8Copula", dimension = as.integer(2), parameters = param,
param.names = c("theta", "delta"), param.lowbnd = c(-Inf, -1), param.upbnd = c(-1, 0),
- family=30, fullname = "90 deg rotated BB8 copula family. Number 30 in CDVine.")
+ family=30, fullname = "90 deg rotated BB8 copula family. Number 30 in VineCopula.")
}
## density ##
@@ -220,7 +220,7 @@
# constructor
r270BB8Copula <- function (param=c(-1,-1)) {
- val <- new("r270BB8Copula", dimension = as.integer(2), parameters = param, param.names = c("theta", "delta"), param.lowbnd = c(-Inf, -1), param.upbnd = c(-1, 0), family=40, fullname = "270 deg rotated BB8 copula family. Number 40 in CDVine.")
+ val <- new("r270BB8Copula", dimension = as.integer(2), parameters = param, param.names = c("theta", "delta"), param.lowbnd = c(-Inf, -1), param.upbnd = c(-1, 0), family=40, fullname = "270 deg rotated BB8 copula family. Number 40 in VineCopula.")
val
}
Modified: pkg/R/ClaytonGumbelCopula.R
===================================================================
--- pkg/R/ClaytonGumbelCopula.R 2013-07-23 11:46:51 UTC (rev 100)
+++ pkg/R/ClaytonGumbelCopula.R 2013-08-01 10:31:15 UTC (rev 101)
@@ -33,7 +33,7 @@
surClaytonCopula <- function (param=1) {
new("surClaytonCopula", dimension = as.integer(2), parameters = param, param.names = c("theta"),
param.lowbnd = 0, param.upbnd = Inf, family=13,
- fullname = "Survival Clayton copula family. Number 13 in CDVine.")
+ fullname = "Survival Clayton copula family. Number 13 in VineCopula.")
}
## density ##
@@ -107,7 +107,7 @@
r90ClaytonCopula <- function (param=-1) {
new("r90ClaytonCopula", dimension = as.integer(2), parameters = param, param.names = c("theta", "delta"),
param.lowbnd = -Inf, param.upbnd = 0, family=23,
- fullname = "90 deg rotated Clayton copula family. Number 23 in CDVine.")
+ fullname = "90 deg rotated Clayton copula family. Number 23 in VineCopula.")
}
## density ##
@@ -167,7 +167,7 @@
r270ClaytonCopula <- function (param=-1) {
new("r270ClaytonCopula", dimension = as.integer(2), parameters = param, param.names = c("theta", "delta"),
param.lowbnd = -Inf, param.upbnd = 0, family=33,
- fullname = "270 deg rotated Clayton copula family. Number 33 in CDVine.")
+ fullname = "270 deg rotated Clayton copula family. Number 33 in VineCopula.")
}
## density ##
@@ -248,7 +248,7 @@
surGumbelCopula <- function (param=1) {
new("surGumbelCopula", dimension = as.integer(2), parameters = param, param.names = c("theta"),
param.lowbnd = 1, param.upbnd = Inf, family=14,
- fullname = "Survival Gumbel copula family. Number 14 in CDVine.")
+ fullname = "Survival Gumbel copula family. Number 14 in VineCopula.")
}
## density ##
@@ -323,7 +323,7 @@
r90GumbelCopula <- function (param=-1) {
new("r90GumbelCopula", dimension = as.integer(2), parameters = param, param.names = c("theta", "delta"),
param.lowbnd = -Inf, param.upbnd = -1, family=24,
- fullname = "90 deg rotated Gumbel copula family. Number 24 in CDVine.")
+ fullname = "90 deg rotated Gumbel copula family. Number 24 in VineCopula.")
}
## density ##
@@ -383,7 +383,7 @@
r270GumbelCopula <- function (param=-1) {
new("r270GumbelCopula", dimension = as.integer(2), parameters = param, param.names = c("theta", "delta"),
param.lowbnd = -Inf, param.upbnd = -1, family=34,
- fullname = "270 deg rotated Gumbel copula family. Number 34 in CDVine.")
+ fullname = "270 deg rotated Gumbel copula family. Number 34 in VineCopula.")
}
## density ##
Modified: pkg/R/joeBiCopula.R
===================================================================
--- pkg/R/joeBiCopula.R 2013-07-23 11:46:51 UTC (rev 100)
+++ pkg/R/joeBiCopula.R 2013-08-01 10:31:15 UTC (rev 101)
@@ -30,7 +30,7 @@
stop("Parameter is outside of the allowed interval (1,Inf).")
new("joeBiCopula", dimension = as.integer(2), parameters = param, param.names = c("theta"),
param.lowbnd = 1, param.upbnd = Inf, family=6,
- fullname = "Joe copula family. Number 6 in CDVine.")
+ fullname = "Joe copula family. Number 6 in VineCopula.")
}
## density ##
@@ -77,7 +77,7 @@
setMethod("tailIndex",signature("joeBiCopula"),linkVineCop.tailIndex)
-## kendall distribution/measure, taken from CDVine:::obs.stat
+## kendall distribution/measure, taken from VineCopula:::obs.stat
kendall.Joe <- function(copula, t){
par = copula at parameters[1]
@@ -109,7 +109,7 @@
stop("Parameter is outside of the allowed interval (1,Inf).")
new("surJoeBiCopula", dimension = as.integer(2), parameters = param, param.names = c("theta", "delta"),
param.lowbnd = 1, param.upbnd = Inf, family=16,
- fullname = "Survival Joe copula family. Number 16 in CDVine.")
+ fullname = "Survival Joe copula family. Number 16 in VineCopula.")
}
## density ##
@@ -186,7 +186,7 @@
stop("Parameter is outside of the allowed interval (-Inf,-1).")
new("r90JoeBiCopula", dimension = as.integer(2), parameters = param, param.names = c("theta", "delta"),
param.lowbnd = -Inf, param.upbnd = -1, family=26,
- fullname = "90 deg rotated Joe copula family. Number 26 in CDVine.")
+ fullname = "90 deg rotated Joe copula family. Number 26 in VineCopula.")
}
## density ##
@@ -248,7 +248,7 @@
stop("Parameter is outside of the allowed interval (-Inf,-1).")
new("r270JoeBiCopula", dimension = as.integer(2), parameters = param, param.names = c("theta", "delta"),
param.lowbnd = -Inf, param.upbnd = -1, family=36,
- fullname = "270 deg rotated Joe copula family. Number 36 in CDVine.")
+ fullname = "270 deg rotated Joe copula family. Number 36 in VineCopula.")
}
## density ##
Modified: pkg/R/linkingVineCopula.R
===================================================================
--- pkg/R/linkingVineCopula.R 2013-07-23 11:46:51 UTC (rev 100)
+++ pkg/R/linkingVineCopula.R 2013-08-01 10:31:15 UTC (rev 101)
@@ -20,7 +20,7 @@
}
#####################################################
-## generic wrapper functions to the CDVine package ##
+## generic wrapper functions to the VineCopula package ##
#####################################################
# density from BiCopPDF
@@ -33,7 +33,7 @@
coplik = RLL_mod_separate(fam, n, u, param)[[7]]
# coplik = .C("LL_mod_seperate", as.integer(fam), as.integer(n), as.double(u[,1]),
# as.double(u[,2]), as.double(param[1]), as.double(param[2]),
-# as.double(rep(0, n)), PACKAGE = "CDVine")[[7]]
+# as.double(rep(0, n)), PACKAGE = "VineCopula")[[7]]
if(log) return(coplik)
else return(exp(coplik))
}
@@ -49,7 +49,7 @@
res <- RarchCDF(fam, n, u, param)[[6]]
# res <- .C("archCDF", as.double(u[,1]), as.double(u[,2]), as.integer(n), as.double(param),
-# as.integer(fam), as.double(rep(0, n)), PACKAGE = "CDVine")[[6]]
+# as.integer(fam), as.double(rep(0, n)), PACKAGE = "VineCopula")[[6]]
return(res)
}
@@ -65,7 +65,7 @@
res <- u1 + u2 - 1 + RarchCDF(fam-10, n, cbind(1-u1,1-u2), param)[[6]]
# res <- u1 + u2 - 1 + .C("archCDF", as.double(1 - u1), as.double(1 - u2), as.integer(n),
# as.double(param), as.integer(fam - 10), as.double(rep(0, n)),
-# PACKAGE = "CDVine")[[6]]
+# PACKAGE = "VineCopula")[[6]]
return(res)
}
@@ -81,7 +81,7 @@
res <- u2 - RarchCDF(fam - 20, n, cbind(1-u1,u2), -param)[[6]]
# u2 - .C("archCDF", as.double(1 - u1), as.double(u2), as.integer(n),
# as.double(-param), as.integer(fam - 20), as.double(rep(0, n)),
-# PACKAGE = "CDVine")[[6]]
+# PACKAGE = "VineCopula")[[6]]
return(res)
}
@@ -97,7 +97,7 @@
res <- u1 - RarchCDF(fam-30, n, cbind(u1,1-u2), -param)[[6]]
# u1 - .C("archCDF", as.double(u1), as.double(1 - u2), as.integer(n),
# as.double(-param), as.integer(fam - 30), as.double(rep(0, n)),
-# PACKAGE = "CDVine")[[6]]
+# PACKAGE = "VineCopula")[[6]]
return(res)
}
@@ -112,7 +112,7 @@
res <- RHfunc1(fam, n, u, param)[[7]]
# .C("Hfunc1", as.integer(fam), as.integer(n), as.double(u[,2]), as.double(u[,1]),
# as.double(param[1]), as.double(param[2]), as.double(rep(0, n)),
-# PACKAGE = "CDVine")[[7]]
+# PACKAGE = "VineCopula")[[7]]
return(res)
}
@@ -126,12 +126,12 @@
res <- RHfunc2(fam, n, u, param)[[7]]
# .C("Hfunc2", as.integer(fam), as.integer(n), as.double(u[,1]), as.double(u[,2]),
# as.double(param[1]), as.double(param[2]), as.double(rep(0, n)),
-# PACKAGE = "CDVine")[[7]]
+# PACKAGE = "VineCopula")[[7]]
return(res)
}
-## random numbers from CDVineSim
+## random numbers from VineCopulaSim
linkVineCop.r <- function (n, copula){
param <- copula at parameters
fam <- copula at family
@@ -140,7 +140,7 @@
tmp <- Rpcc(fam, n, param)[[7]]
# .C("pcc", as.integer(n), as.integer(2), as.integer(fam), as.integer(1),
# as.double(param[1]), as.double(param[2]), as.double(rep(0, n * 2)),
-# PACKAGE = "CDVine")[[7]]
+# PACKAGE = "VineCopula")[[7]]
return(matrix(tmp, ncol = 2))
}
@@ -150,8 +150,8 @@
-# ## transform a fit from CDVine to a list of copula objects
-# castCDvine <- function(cdvEst) {
+# ## transform a fit from VineCopula to a list of copula objects
+# castVineCopula <- function(cdvEst) {
# copulas <- NULL
# for(i in 1: length(cdvEst$family)) {
# par1 <- cdvEst$par[i]
Modified: pkg/R/partialDerivatives.R
===================================================================
--- pkg/R/partialDerivatives.R 2013-07-23 11:46:51 UTC (rev 100)
+++ pkg/R/partialDerivatives.R 2013-08-01 10:31:15 UTC (rev 101)
@@ -393,7 +393,7 @@
setGeneric("kendallDistribution")
## Clayton
-## kendall distribution/measure, taken from CDVine:::obs.stat
+## kendall distribution/measure, taken from VineCopula:::obs.stat
kendall.Clayton <- function(copula, t){
par = copula at parameters
@@ -410,7 +410,7 @@
function(copula) return(function(t) kendall.Clayton(copula, t)))
## Gumbel
-## kendall distribution/measure, taken from CDVine:::obs.stat
+## kendall distribution/measure, taken from VineCopula:::obs.stat
kendall.Gumbel <- function(copula, t){
par = copula at parameters
@@ -427,7 +427,7 @@
function(copula) return(function(t) kendall.Gumbel(copula, t)))
## Frank
-## kendall distribution/measure, taken from CDVine:::obs.stat
+## kendall distribution/measure, taken from VineCopula:::obs.stat
kendall.Frank <- function(copula, t){
par = copula at parameters
Modified: spcopula_0.1-1.tar.gz
===================================================================
(Binary files differ)
Modified: spcopula_0.1-1.zip
===================================================================
(Binary files differ)
Added: thoughts-Ben/Meuse_spcopula_estimation.R
===================================================================
--- thoughts-Ben/Meuse_spcopula_estimation.R (rev 0)
+++ thoughts-Ben/Meuse_spcopula_estimation.R 2013-08-01 10:31:15 UTC (rev 101)
@@ -0,0 +1,56 @@
+## librarys ##
+library(spcopula)
+library(evd)
+
+## dataset - spatial poionts data.frame ##
+data(meuse)
+coordinates(meuse) = ~x+y
+dataSet <- meuse
+
+spplot(meuse,"zinc")
+
+## margins ##
+hist(dataSet[["zinc"]],freq=F,n=30,ylim=c(0,0.0035))
+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")
+curve(dlnorm(x,meanLog,sdLog),add=T,col="green")
+
+ks.test(dataSet[["zinc"]],pgev,loc,scale,shape) # p: 0.07
+ks.test(dataSet[["zinc"]],plnorm,meanLog,sdLog) # p: 0.03
+
+## lag classes ##
+bins <- calcBins(dataSet,var="zinc",nbins=10,cutoff=800)
+
+# transform data to the unit interval
+bins$lagData <- lapply(bins$lagData, function(x) cbind(rank(x[,1])/(nrow(x)+1),rank(x[,2])/(nrow(x)+1)))
+
+abline(h=0,col="red")
+modelCoeff <- lm(bins$lagCor[1:7]~bins$meanDists[1:7])$coefficients
+abline(modelCoeff,col="blue") # 8 lags
+range <- as.numeric(-modelCoeff[1]/modelCoeff[2])
+
+## calculate parameters for Kendall's tau function ##
+# either linear
+calcKTau <- function(h) pmax(h*modelCoeff[2]+modelCoeff[1],0)
+curve(calcKTau,0, 1000, col="red", add=T)
+
+# or polynomial
+calcKTau <- fitCorFun(bins)
+curve(calcKTau,0, 1000, col="purple",add=T)
+
+## find best fitting copula per lag class
+loglikTau <- loglikByCopulasLags(bins, calcKTau)
+bestFitTau <- apply(apply(loglikTau[1:7,], 1, rank),2,function(x) which(x==5))
+
+## set-up a spatial Copula ##
+spCop <- spCopula(components=list(normalCopula(0), tCopula(0, dispstr = "un"),
+ frankCopula(1), normalCopula(0), claytonCopula(0),
+ claytonCopula(0), claytonCopula(0)),
+ distances=c(bins$meanDists[1:7],range),
+ spDepFun=calcKTau, unit="m")
+
Added: thoughts-Ben/empiricalCopula_COPULA.R
===================================================================
--- thoughts-Ben/empiricalCopula_COPULA.R (rev 0)
+++ thoughts-Ben/empiricalCopula_COPULA.R 2013-08-01 10:31:15 UTC (rev 101)
@@ -0,0 +1,80 @@
+########################################
+## ##
+## an empirical copula representation ##
+## ##
+########################################
+
+validEmpCopula <- function(object) {
+ if(ncol(object at sample) != object at dimension)
+ return("Dimension of the copula and the sample do not match.")
+ else
+ return(TRUE)
+}
+
+setClass("empiricalCopula",
+ representation = representation("copula", sample="matrix"),
+ validity = validEmpCopula,
+ contains = list("copula")
+)
+
+# constructor
+empiricalCopula <- function (sample=NULL, copula) {
+ if(is.null(sample) && missing(copula))
+ stop("At least one parameter of copula or sample must be provided.")
+ if(is.null(sample))
+ return(genEmpCop(copula))
+ if(missing(copula))
+ return(new("empiricalCopula", dimension = as.integer(ncol(sample)),
+ parameters = as.numeric(NA), param.names = "unknown",
+ param.lowbnd = as.numeric(NA), param.upbnd = as.numeric(NA),
+ fullname = "Unkown empirical copula based on a sample.",
+ sample=sample))
+ new("empiricalCopula", dimension = copula at dimension,
+ parameters = copula at parameters, param.names = copula at param.names,
+ param.lowbnd = copula at param.lowbnd, param.upbnd = copula at param.upbnd,
+ fullname = paste("Empirical copula derived from",copula at fullname),
+ sample=sample)
+}
+
+# simplified constructor
+genEmpCop <- function(copula, sample.size=1e5) {
+ cat("Note: the copula will be empirically represented by a sample of size:", sample.size, "\n")
+ empiricalCopula(rCopula(sample.size,copula), copula)
+}
+
+
+## density, not yet needed and hence not implemented ##
+
+## jcdf ##
+# from package copula
+pempCop.C <- function(u, copula) {
+ # r-forge "hack", to be removed after release of copula 0.999-6
+ if(length(formals(Cn))==2) {
+ return(Cn(copula at sample,u))
+ }
+ if (length(formals(Cn))== 5) {
+ return(Cn(u, copula at sample, do.pobs=FALSE, offset=0, method="C"))
+ }
+ stop(length(formals(Cn))) # to trace back changes
+}
+
+setMethod("pCopula", signature("numeric", "empiricalCopula"),
+ function(u, copula, ...) {
+ pempCop.C(matrix(u,ncol=copula at dimension),copula)
+ })
+setMethod("pCopula", signature("matrix", "empiricalCopula"), pempCop.C)
+
+
+tauempCop <- function(copula){
+ # CDVine:::fasttau(copula at sample[,1],copula at sample[,2]) # a fast version from CDVine
+ cor(copula at sample,method="kendall")
+}
+
+setMethod("tau",signature("asCopula"),tauempCop)
+
+
+rhoempCop <- function(copula){
+ cor(copula at sample,method="spearman")
+}
+
+setMethod("rho",signature("asCopula"),rhoempCop)
\ No newline at end of file
Added: thoughts-Ben/tawn3pCopula-stub.R
===================================================================
--- thoughts-Ben/tawn3pCopula-stub.R (rev 0)
+++ thoughts-Ben/tawn3pCopula-stub.R 2013-08-01 10:31:15 UTC (rev 101)
@@ -0,0 +1,476 @@
+#################
+## tawn copula ##
+#################
+
+library(copula)
+
+Atawn3p <- function(t, param = c(0.9302082, 1, 8.355008)) {
+ alpha <- param[1]
+ beta <- param[2]
+ theta <- param[3]
+ (1-beta)*(t) + (1-alpha)*(1-t) + ((alpha*(1-t))^theta+(beta*t)^theta)^(1/theta)
+# 1-beta + (beta-alpha)*t + ((alpha*t)^theta + (beta*(1-t))^theta)^(1/theta)
+}
+
+curve(Atawn3p)
+
+tawn3pCopula <- function (param = c(0.5, 0.5, 2)) {
+ # A(t) = (1-beta)*(1-t) + (1-alpha)*t + ((alpha*(1-t))^theta+(beta*t)^theta)^(1/theta)
+ # C(u1,u2) = exp(log(u1*u2) * A(log(u2)/log(u1*u2)))
+
+ cdf <- expression(exp(log(u1*u2)*((1-beta)*(log(u2)/log(u1*u2)) +
+ (1-alpha)*(1-log(u2)/log(u1*u2)) +
+ ((alpha*(1-log(u2)/log(u1*u2)))^theta+(beta*log(u2)/log(u1*u2))^theta)^(1/theta))))
+ dCdU1 <- D(cdf, "u1")
+ pdf <- D(dCdU1, "u2")
+
+ new("tawn3pCopula", dimension = 2L, exprdist = c(cdf = cdf, pdf = pdf),
+ parameters = param, param.names = c("alpha", "beta", "theta"),
+ param.lowbnd = c(0,0,1), param.upbnd = c(1,1,Inf),
+ fullname = "Tawn copula family with three parameters; Extreme value copula")
+}
+
+setClass("tawn3pCopula", representation(exprdist = "expression"),
+ contains = "evCopula")
+
+
+dtawn3pCopula <- function(u, copula, log=FALSE, ...) {
+ dim <- copula at dimension
+ for (i in 1:dim) {
+ assign(paste("u", i, sep=""), u[,i])
+ }
+ alpha <- copula at parameters[1]
+ beta <- copula at parameters[2]
+ theta <-copula at parameters[3]
+
+ val <- c(eval(copula at exprdist$pdf))
+ ## FIXME: improve log-case
+ if(log) log(val) else val
+}
+
+setMethod("dCopula", signature("matrix", "tawn3pCopula"), dtawn3pCopula)
+setMethod("dCopula", signature("numeric", "tawn3pCopula"),dtawn3pCopula)
+
+ptawn3pCopula <- function(u, copula, ...) {
+ dim <- copula at dimension
+ for (i in 1:dim) {
+ assign(paste("u", i, sep=""), u[,i])
+ }
+ alpha <- copula at parameters[1]
+ beta <- copula at parameters[2]
+ theta <-copula at parameters[3]
+
+ val <- c(eval(copula at exprdist$cdf))
+}
+
+setMethod("pCopula", signature("matrix", "tawn3pCopula"), ptawn3pCopula)
+setMethod("pCopula", signature("numeric", "tawn3pCopula"), ptawn3pCopula)
+
+
+persp(tawn3pCopula(c(0.25, 0.75, 2)), dCopula)
+persp(tawn3pCopula(c(0.5, 1, 20)), pCopula)
+
+plot(1-rtTriples[,c(1,3)])
+
+dCopula(c(0.15,0.95), tawn3pCopula(c(0.5, 0.5, 2)))
+
+
+tawnFit <- fitCopula(tawn3pCopula(), 1-as.matrix(rtTriples[,c(1,3)]), hideWarnings=F,estimate.variance=F,
+ start=c(0.9, 1, 8), method="mpl", lower=c(0,0,1), upper=c(1, 1, 10),
+ optim.method="L-BFGS-B",)
+tawnFit at loglik # 742
+tawnFit at copula
+
+fitCopula(gumbelCopula(5),1-as.matrix(rtTriples[,c(1,3)]))@loglik # 723
+
+dLeaf <- dCopula(as.matrix(rtTriples[,c(1,3)]), spcopula:::leafCopula())
+sum(log(dLeaf[dLeaf>0]))
+
+persp(tawnFit at copula, dCopula)
+contour(tawnFit at copula, dCopula, levels=c(0,0.5,1,2,4,8,100), asp=1)
+
+sum(dCopula(as.matrix(rtTriples[,c(1,3)]), cop13, log=T))
+sum(dCopula(1-as.matrix(rtTriples[,c(1,3)]), tawnFit at copula, log=T))
+
+plot(1-as.matrix(rtTriples[,c(1,3)]),asp=1,cex=0.5)
+curve(x^(tawnFit at copula@parameters[1]),add=T, col="red")
+abline(0,1,col="grey")
+
+copula:::fitCopula.ml
+
+
+tawn3pCopula()
+
+persp(tawn3pCopula(),dCopula)
+
+
+###
+# h(t)
+
+hist(log(1-rtTriples[,3])/log((1-rtTriples[,1])*(1-rtTriples[,3])),n=20,
+ xlim=c(0,1), freq=F, add=T, col="blue")
+
+tSmpl <- log(1-rtTriples[,3])/log((1-rtTriples[,1])*(1-rtTriples[,3]))
+ ((1-rtTriples[,1])-(1-rtTriples[,3]))*(0.5)+0.5
+
+dlogNorm <- function(x) dlnorm(x, mean(log(tSmpl)), sd(log(tSmpl)))
+aGevPar <- fgev(tSmpl)$estimate
+dGev <- function(x) dgev(x, 0.46938, 0.05057, 0.01720)
+dGamma <- function(x) dgamma(x, 60.40556, 120.92807)
+
+hist(tSmpl, freq=F, xlim=c(0,1), n=20, add=F)
+curve(dlogNorm, add=T)
+curve(dGev, add=T, col="red")
+curve(dGamma, add=T, col="purple")
+
+sum(log(dGev(tSmpl))) # 690
+sum(log(dlogNorm(tSmpl))) # 666
+sum(log(dGamma(tSmpl))) # 658
+
+
+Afit <- function(t) {
+ res <- t
+ res[res == 0] <- 1
+ intFun <- function(z) (pgev(z, 0.46938, 0.05057, 0.01720)-z)/(z-z^2)
+ for(i in which(res != 1)) {
+ res[i] <- exp(integrate(intFun,0,t[i])$value)
+ }
+ return(res)
+}
+
+curve(Afit,ylim=c(0.5,1))
+abline(0, 1, col="grey")
+abline(1, -1, col="grey")
+curve(Atawn3p,col="red",add=T)
+
+
+optFun <- function(param) {
+ -sum(log(dgamma(tSmpl,param[1],param[2])))
+}
+
+optim(c(1,0.5),optFun)
+
+
+library(evd)
+
+?dnorm
+
+##
+AFun <- function(u, v, cdf=function(z) pgev(z, 0.46938, 0.05057, 0.01720)) {
+ t <- log(v)/log(u*v)
+
+ # rough limit case corrections
+ t[v == 0 | u == 1 | v == 1 | u == 0] <- 1
+
+ res <- t
+ intFun <- function(z) (cdf(z)-z)/(z-z^2)
+ for(i in which(res < 1 & res > 0)) {
+ res[i] <- exp(integrate(intFun,0,t[i],stop.on.error=F,subdivisions=1000L)$value)
+ }
+ return(res) #pmax(pmax(res,t),1-t))
+}
+
+plot(u,AFun(u),typ="l",ylim=c(0,1))
+points(u,log(0.3)/log(0.3*u),typ="l",col="red") # slope less than this one
+points(u,1-log(0.3)/log(0.3*u),typ="l",col="red") # slope larger than this one
+
+AFunInT <- function(t, cdf=function(z) pgev(z, 0.5,
+ 0.05,
+ 0.017)) {
+ res <- t
+ res[res == 0] <- 1
+ intFun <- function(z) (cdf(z)-z)/(z-z^2)
+
+ for(i in which(res < 1 & res > 0)) {
+ res[i] <- exp(integrate(intFun,0,t[i])$value)
+ }
+ return(res)
+}
+
+curve(AFunInT)
+abline(0,1)
+abline(1,-1)
+
+fooFun <- function(z) pgev(z, 0.5, 0.05, 0.017)
+curve(fooFun,add=T)
+abline(v=0.5)
+
+aGevPar
+newFit at estimate
+
+curve(AFunInT)
+abline(0,1)
+abline(1,-1)
+
+
+dduAFun <- function(u, v, cdf=function(z) pgev(z, 0.46938, 0.05057, 0.01720)) {
+ logV <- log(v)
+ logUV <- log(u*v)
+
+ t <- logV/logUV
+ # limit cases:
+ # v == 0 -> t=1
+ # v == 1 -> t=0
+ # u == 1 -> t=1
+ # u == 0 -> t=0
+ t[v == 0] <- 1
+ t[v == 1] <- 0
+ t[u == 1] <- 1
+ t[u == 0] <- 0
+
+ res <- t
+ res[t == 0] <- -10
+ res[t == 1] <- (-logV[t == 1])/u[t == 1]/logUV[t == 1]^2
+
+ boolDo <- t > 0 & t < 1
+ tDo <- t[boolDo]
+ uDo <- u[boolDo]
+ vDo <- v[boolDo]
+
+ res[boolDo] <- AFun(uDo, vDo, cdf)*(cdf(tDo)-tDo)/(tDo-tDo^2)*(-logV[boolDo])/(u[boolDo]*logUV[boolDo]^2)
+
+ return(res)
+}
+
+plot(u, dduAFun(u),typ="l")
+points(u,-log(0.3)/log(0.3*u)^2/u, typ="l",col="red") # slope less than this one
+points(u, log(0.3)/log(0.3*u)^2/u, typ="l",col="red") # slope larger than this one
+
+curve(AFun,0,1,add=F,ylim=c(-2,1))
+curve(dduAFun, add=T, col="green")
+abline(h=0,col="grey")
+abline(v=uniroot(dduAFun,c(0.2,0.4))$root, col="grey")
+
+ddvAFun <- function(u, v=rep(0.3,length(u)), cdf=function(z) pgev(z, 0.46938, 0.05057, 0.01720)) {
+ logV <- log(v)
+ logUV <- log(u*v)
+
+ t <- logV/logUV
+ # limit cases:
+ # v == 0 -> t=1
+ # v == 1 -> t=0
+ # u == 1 -> t=1
+ # u == 0 -> t=0
+ t[v == 0] <- 1
+ t[v == 1] <- 0
+ t[u == 1] <- 1
+ t[u == 0] <- 0
+
+ res <- t
+ res[t == 0] <- -(logUV[t == 0]-logV[t == 0])/(logUV[t == 0]^2*v[t == 0])
+ res[t == 1] <- (logUV[t == 1]-logV[t == 1])/(logUV[t == 1]^2*v[t == 1])
+
+ boolDo <- t > 0 & t < 1
+ tDo <- t[boolDo]
+ uDo <- u[boolDo]
+ vDo <- v[boolDo]
+
+ res[boolDo] <- AFun(uDo, vDo, cdf)*(cdf(tDo)-tDo)/(tDo-tDo^2)*(logUV[boolDo]-logV[boolDo])/(vDo*logUV[boolDo]^2)
+
+ return(res)
+}
+
+
+plot(v, ddvAFun(v),typ="l")
+points(v,-(log(0.3*v)-log(v))/(log(0.3*v)^2*v), typ="l",col="red") # slope less than this one
+points(v, (log(0.3*v)-log(v))/(log(0.3*v)^2*v), typ="l",col="red") # slope larger than this one
+
+curve(AFun,0,1,add=F,ylim=c(-2,1))
+curve(ddvAFun, add=T, col="green")
+abline(h=0,col="grey")
+abline(v=uniroot(dduAFun,c(0.2,0.4))$root, col="grey")
+
+dduvAFun <- function(u, v=rep(0.3,length(u)),
+ cdf=function(z) pgev(z, 0.46938, 0.05057, 0.01720),
+ pdf=function(z) dgev(z, 0.46938, 0.05057, 0.01720)) {
+ logV <- log(v)
+ logUV <- log(u*v)
+
+ .g <- function(z) {
+ (cdf(z)-z)/(z-z^2)
+ }
+
+ .ddzg <- function(z) {
+ ((2*z-1)*cdf(z)-z*((z-1)*pdf(z)+z))/((z-1)^2*z^2)
+ }
+
+ t <- logV/logUV
+ # limit cases
+ t[v == 0] <- 1
+ t[v == 1] <- 0
+ t[u == 1] <- 1
+ t[u == 0] <- 0
+
+ p1 <- dduAFun(u, v, cdf)*.g(t)*u*(logUV-logV)*logUV^2
+ p2 <- AFun(u, v, cdf)*(logV*(logV-logUV) * .ddzg(t) + (2*logV-logUV)*logUV * .g(t) )
+
+ (p1+p2)/(u*v*logUV^4)
+}
+
+curve(ddvAFun,0,1,ylim=c(-4,2),asp=1, add=T)
+curve(dduvAFun,0.01,0.99,add=F)
+abline(h=0, col="grey")
+
+curve(AFun, add=T, col="red")
+curve(dduAFun)
+# curve(ddvAFun,add=T, col="red")
+# curve(dduvAFun, add=F)
+
+CDFevCopula <- function(u, v=rep(0.3,length(u)),
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/spcopula -r 101
More information about the spcopula-commits
mailing list