[spcopula-commits] r95 - / pkg pkg/R pkg/man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri May 3 17:23:07 CEST 2013
Author: ben_graeler
Date: 2013-05-03 17:23:07 +0200 (Fri, 03 May 2013)
New Revision: 95
Added:
pkg/man/calcSpTreeDists.Rd
pkg/man/dropSpTree.Rd
Modified:
pkg/DESCRIPTION
pkg/NAMESPACE
pkg/R/BB1copula.R
pkg/R/BB6copula.R
pkg/R/BB7copula.R
pkg/R/BB8copula.R
pkg/R/Classes.R
pkg/R/ClaytonGumbelCopula.R
pkg/R/asCopula.R
pkg/R/cqsCopula.R
pkg/R/joeBiCopula.R
pkg/R/linkingVineCopula.R
pkg/R/partialDerivatives.R
pkg/R/spCopula.R
pkg/R/spVineCopula.R
pkg/R/spatialPreparation.R
pkg/R/utilities.R
pkg/man/BB8Copula-class.Rd
pkg/man/calcBins.Rd
pkg/man/condSpVine.Rd
pkg/man/cqsCopula-class.Rd
pkg/man/cqsCopula.Rd
pkg/man/loglikByCopulasLags.Rd
pkg/man/neighbourhood-class.Rd
pkg/man/neighbourhood.Rd
pkg/man/spCopPredict.Rd
pkg/man/spVineCopula-class.Rd
pkg/man/spVineCopula.Rd
spcopula_0.1-1.tar.gz
spcopula_0.1-1.zip
Log:
- multilevel spatial copulas!
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2013-04-23 14:17:11 UTC (rev 94)
+++ pkg/DESCRIPTION 2013-05-03 15:23:07 UTC (rev 95)
@@ -2,7 +2,7 @@
Type: Package
Title: copula driven spatial analysis
Version: 0.1-1
-Date: 2013-04-23
+Date: 2013-05-03
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/NAMESPACE
===================================================================
--- pkg/NAMESPACE 2013-04-23 14:17:11 UTC (rev 94)
+++ pkg/NAMESPACE 2013-05-03 15:23:07 UTC (rev 95)
@@ -27,6 +27,7 @@
# spatial
export(getNeighbours)
export(calcBins)
+export(calcSpTreeDists, dropSpTree)
# fitting
export(fitCorFun, loglikByCopulasLags, fitSpCopula, composeSpCopula)
Modified: pkg/R/BB1copula.R
===================================================================
--- pkg/R/BB1copula.R 2013-04-23 14:17:11 UTC (rev 94)
+++ pkg/R/BB1copula.R 2013-05-03 15:23:07 UTC (rev 95)
@@ -25,7 +25,7 @@
)
# constructor
-BB1Copula <- function (param) {
+BB1Copula <- function (param=c(1,1)) {
if (any(is.na(param) | param >= c(Inf,Inf) | param[1] <= 0 | param[2] < 1))
stop(paste("Parameter values out of bounds: theta: (0,Inf), delta: [1,Inf)."))
new("BB1Copula", dimension = as.integer(2), parameters = param,
@@ -38,7 +38,8 @@
function(u, copula, log) {
linkVineCop.PDF(matrix(u,ncol=copula at dimension),copula, log)
})
-setMethod("dCopula", signature("matrix","BB1Copula"), function(u, copula, log) linkVineCop.PDF(u, copula, log))
+setMethod("dCopula", signature("matrix","BB1Copula"),
+ function(u, copula, log) linkVineCop.PDF(u, copula, log))
## jcdf ##
setMethod("pCopula", signature("numeric","BB1Copula"),
@@ -95,7 +96,7 @@
)
# constructor
-surBB1Copula <- function (param) {
+surBB1Copula <- function (param=c(1,1)) {
if (any(is.na(param) | param >= c(Inf,Inf) | param[1] <= 0 | param[2] < 1))
stop(paste("Parameter values out of bounds: theta: (0,Inf), delta: [1,Inf)."))
new("surBB1Copula", dimension = as.integer(2), parameters = param,
@@ -162,7 +163,7 @@
)
# constructor
-r90BB1Copula <- function (param) {
+r90BB1Copula <- function (param=c(-1,-1)) {
if (any(is.na(param) | param[1] >= 0 | param[2] > -1 | param <= c(-Inf,-Inf)))
stop(paste("Parameter values out of bounds: theta: (-Inf,0), delta: (-Inf,-1]."))
new("r90BB1Copula", dimension = as.integer(2), parameters = param,
@@ -216,7 +217,7 @@
)
# constructor
-r270BB1Copula <- function (param) {
+r270BB1Copula <- function (param=c(-1,-1)) {
if (any(is.na(param) | param[1] >= 0 | param[2] > -1 | param <= c(-Inf,-Inf)))
stop(paste("Parameter values out of bounds: theta: (-Inf,0), delta: (-Inf,-1]."))
new("r270BB1Copula", dimension = as.integer(2), parameters = param,
Modified: pkg/R/BB6copula.R
===================================================================
--- pkg/R/BB6copula.R 2013-04-23 14:17:11 UTC (rev 94)
+++ pkg/R/BB6copula.R 2013-05-03 15:23:07 UTC (rev 95)
@@ -27,7 +27,7 @@
)
# constructor
-BB6Copula <- function (param) {
+BB6Copula <- function (param=c(1,1)) {
if (any(is.na(param) | param >= c(Inf, Inf) | param < c(1,1)))
stop("Parameter value(s) out of bound(s): theta: [1,Inf), delta: [1,Inf).")
new("BB6Copula", dimension = as.integer(2), parameters = param,
@@ -98,7 +98,7 @@
)
# constructor
-surBB6Copula <- function (param) {
+surBB6Copula <- function (param=c(1,1)) {
if (any(is.na(param) | param >= c(Inf, Inf) | param < c(1,1)))
stop("Parameter value(s) out of bound(s): theta: [1,Inf), delta: [1,Inf).")
new("surBB6Copula", dimension = as.integer(2), parameters = param,
@@ -165,7 +165,7 @@
)
# constructor
-r90BB6Copula <- function (param) {
+r90BB6Copula <- function (param=c(-1,-1)) {
if (any(is.na(param) | param > c(-1,-1) | param <= c(-Inf,-Inf)))
stop("Parameter value out of bound: theta: (-Inf,1], delta: (-Inf,1].")
new("r90BB6Copula", dimension = as.integer(2), parameters = param,
@@ -219,7 +219,7 @@
)
# constructor
-r270BB6Copula <- function (param) {
+r270BB6Copula <- function (param=c(-1,-1)) {
if (any(is.na(param) | param > c(-1,-1) | param <= c(-Inf,-Inf)))
stop("Parameter value out of bound: theta: (-Inf,1], delta: (-Inf,1].")
new("r270BB6Copula", dimension = as.integer(2), parameters = param,
Modified: pkg/R/BB7copula.R
===================================================================
--- pkg/R/BB7copula.R 2013-04-23 14:17:11 UTC (rev 94)
+++ pkg/R/BB7copula.R 2013-05-03 15:23:07 UTC (rev 95)
@@ -25,7 +25,7 @@
)
# constructor
-BB7Copula <- function (param) {
+BB7Copula <- function (param=c(1,1)) {
if (any(is.na(param) | param >= c(Inf, Inf) | param[1] < 1 | param[2] <= 0))
stop(paste("Parameter values out of bounds: theta: [1,Inf), delta: (0,Inf)."))
new("BB7Copula", dimension = as.integer(2), parameters = param,
@@ -98,7 +98,7 @@
)
# constructor
-surBB7Copula <- function (param) {
+surBB7Copula <- function (param=c(1,1)) {
if (any(is.na(param) | param >= c(Inf, Inf) | param[1] < 1 | param[2] <= 0))
stop(paste("Parameter values out of bounds: theta: [1,Inf), delta: (0,Inf)."))
new("surBB7Copula", dimension = as.integer(2), parameters = param,
@@ -167,7 +167,7 @@
)
# constructor
-r90BB7Copula <- function (param) {
+r90BB7Copula <- function (param=c(-1,-1)) {
if (any(is.na(param) | param[1] > -1 | param[2] >= 0 | param <= c(-Inf,-Inf)))
stop(paste("Parameter values out of bounds: theta: (-Inf,-1], delta: (-Inf,0)."))
new("r90BB7Copula", dimension = as.integer(2), parameters = param,
@@ -221,7 +221,7 @@
)
# constructor
-r270BB7Copula <- function (param) {
+r270BB7Copula <- function (param=c(-1,-1)) {
if (any(is.na(param) | param[1] > -1 | param[2] >= 0 | param <= c(-Inf,-Inf)))
stop(paste("Parameter values out of bounds: theta: (-Inf,-1], delta: (-Inf,0)."))
new("r270BB7Copula", dimension = as.integer(2), parameters = param,
Modified: pkg/R/BB8copula.R
===================================================================
--- pkg/R/BB8copula.R 2013-04-23 14:17:11 UTC (rev 94)
+++ pkg/R/BB8copula.R 2013-05-03 15:23:07 UTC (rev 95)
@@ -27,7 +27,7 @@
)
# constructor
-BB8Copula <- function (param) {
+BB8Copula <- function (param=c(1,1)) {
if (any(is.na(param)) | param[1] >= Inf | param[2] > 1 | param[1] < 1 | param[2] <= 0)
stop("Parameter value out of bound: theta: [1,Inf), delta: (0,1].")
new("BB8Copula", dimension = as.integer(2), parameters = param,
@@ -98,7 +98,7 @@
)
# constructor
-surBB8Copula <- function (param) {
+surBB8Copula <- function (param=c(1,1)) {
if (any(is.na(param)) | param[1] >= Inf | param[2] > 1 | param[1] < 1 | param[2] <= 0)
stop("Parameter value out of bound: theta: [1,Inf), delta: (0,1].")
new("surBB8Copula", dimension = as.integer(2), parameters = param,
@@ -165,7 +165,7 @@
)
# constructor
-r90BB8Copula <- function (param) {
+r90BB8Copula <- function (param=c(-1,-1)) {
if (any(is.na(param) | param[1] > -1 | param[2] >= 0 | param[1] <= -Inf | param[2] < -1))
stop("Parameter value out of bound: theta: (-Inf,-1], delta: [-1,0).")
new("r90BB8Copula", dimension = as.integer(2), parameters = param,
@@ -219,7 +219,7 @@
)
# constructor
-r270BB8Copula <- function (param) {
+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
}
@@ -257,4 +257,32 @@
setMethod("rCopula", signature("numeric","r270BB8Copula"), linkVineCop.r)
setMethod("tau",signature("r270BB8Copula"),linkVineCop.tau)
-setMethod("tailIndex",signature("r270BB8Copula"),linkVineCop.tailIndex)
\ No newline at end of file
+setMethod("tailIndex",signature("r270BB8Copula"),linkVineCop.tailIndex)
+
+### et union
+
+setClassUnion("twoParamBiCop",c("BB1Copula","BB6Copula","BB7Copula","BB8Copula",
+ "surBB1Copula","surBB6Copula","surBB7Copula","surBB8Copula",
+ "r90BB1Copula","r90BB6Copula","r90BB7Copula","r90BB8Copula",
+ "r270BB1Copula","r270BB6Copula","r270BB7Copula","r270BB8Copula"))
+
+fitCopula.twoParamBiCop <- function(copula, data, method = "mpl",
+ estimate.variance = FALSE) {
+ stopifnot(method=="mpl")
+ fit <- BiCopEst(data[,1], data[,2], copula at family, "mle",
+ se=estimate.variance)
+
+ if(!estimate.variance) {
+ fit$se <- NA
+ fit$se2 <- NA
+ }
+
+ copFit <- copulaFromFamilyIndex(copula at family, fit$par, fit$par2)
+ new("fitCopula", estimate = c(fit$par, fit$par2), var.est = cbind(fit$se, fit$se2),
+ method = "estimation based on 'maximum pseudo-likelihood' via BiCopEst",
+ loglik = sum(dCopula(data, copFit, log=T)),
+ fitting.stats=list(convergence = as.integer(NA)), nsample = nrow(data),
+ copula=copFit)
+}
+
+setMethod("fitCopula", signature("twoParamBiCop"), fitCopula.twoParamBiCop)
\ No newline at end of file
Modified: pkg/R/Classes.R
===================================================================
--- pkg/R/Classes.R 2013-04-23 14:17:11 UTC (rev 94)
+++ pkg/R/Classes.R 2013-05-03 15:23:07 UTC (rev 95)
@@ -17,9 +17,18 @@
# the lower bound of the parameter a dependening on the parameter b
limA <- function (b) {
- (b-3-sqrt(9+6*b-3*b^2))/2
+ stopifnot(abs(b) <= 1)
+ 0.5*(-sqrt(-3*b^2+6*b+9)+b-3)
}
+# the lower and upper bound of the parameter b dependening on the parameter a
+limB <- function (a) {
+ stopifnot(a <=1 & a >= -3)
+ if(a>-2)
+ return(c(-1,1))
+ pmax(pmin(0.5*(c(-1,1)*(sqrt(3)*sqrt(-a^2-2*a+3))+a+3),1),-1)
+}
+
setClass("asCopula",
representation = representation("copula"),
validity = validAsCopula,
@@ -29,10 +38,29 @@
####
## a symmetric copula with cubic and quadratic sections
-validCqsCopula <- validAsCopula
+validCqsCopula <- function(object) {
+ if (object at dimension != 2)
+ return("Only copulas with cubic quadratic sections of dimension 2 are supported.")
+ param <- object at parameters
+ upper <- object at param.upbnd
+ lower <- object at param.lowbnd
+ if (length(param) != length(upper))
+ return("Parameter and upper bound have non-equal length")
+ if (length(param) != length(lower))
+ return("Parameter and lower bound have non-equal length")
+ if (any(is.na(param) | param > upper | param < lower))
+ return("Parameter value out of bound")
+ if (length(object at fixed >0)){
+ if(!("a" %in% object at fixed | "b" %in% object at fixed))
+ return("The slot fixed may only refer to \"a\" or \"b\".")
+ if ("a" %in% object at fixed & "b" %in% object at fixed)
+ return("Only one of the parameters may be kept fixed.")
+ }
+ else return (TRUE)
+}
setClass("cqsCopula",
- representation = representation("copula"),
+ representation = representation("copula",fixed="character"),
validity = validCqsCopula,
contains = list("copula")
)
@@ -91,7 +119,8 @@
# assings valid parameters to the copulas involved
validSpCopula <- function(object) {
- if (length(object at components) != length(object at distances)) return("Length of components does not equal length of distances. \n Note: The last distance must give the range and it is automatically associated with the indepenence copula.")
+ if (length(object at components) != length(object at distances))
+ return("Length of components does not equal length of distances. \n Note: The last distance must give the range and it is automatically associated with the indepenence copula.")
check.upper <- NULL
check.lower <- NULL
@@ -146,7 +175,7 @@
return("Number of provided copulas does not match given dimension.")
if(!any(unlist(lapply(object at copulas,function(x) is(x,"copula")))))
return("Not all provided copulas in your list are indeed copulas.")
- else return (TRUE)
+ return (TRUE)
}
setOldClass("RVineMatrix")
@@ -163,13 +192,22 @@
## Spatial Vine Copula ##
#########################
-validSpVineCopula <- function(object) {
- return(validObject(object at spCop)&validObject(object at vineCop))
+validMixedSpVineCopula <- function(object) {
+ return(all(sapply(object at spCop,validSpCopula) & validObject(object at topCop)))
}
-setClass("spVineCopula", representation("copula",spCop="spCopula",vineCop="vineCopula"),
- validity = validSpVineCopula, contains=list("copula"))
+setClass("mixedSpVineCopula", representation("copula", spCop="list", topCop="copula"),
+ validity = validMixedSpVineCopula, contains=list("copula"))
+validPureSpVineCopula <- function(object) {
+ return(all(sapply(object at spCop,validSpCopula)))
+}
+
+setClass("pureSpVineCopula", representation("copula", spCop="list"),
+ validity = validPureSpVineCopula, contains=list("copula"))
+
+setClassUnion("spVineCopula",c("mixedSpVineCopula","pureSpVineCopula"))
+
########################################
## spatial classes providing the data ##
########################################
@@ -193,19 +231,30 @@
validNeighbourhood <- function(object) {
sizeN <- ncol(object at distances)+1
nVars <- length(object at var)
- if (nrow(object at data) != nrow(object at distances)) return("Data and distances have unequal number of rows.")
- if (ncol(object at data) %% (sizeN-object at prediction) != 0) return("Data and distances have non matching number of columns.")
- if (nrow(object at data) != nrow(object at index)) return("Data and index have unequal number of rows.")
- if (ncol(object at distances) != ncol(object at index)) return("Data and index have unequal number of columns.")
- if (ncol(object at data) != (sizeN-object at prediction) * nVars) return(paste("Number of columns in data does not equal the product of the neighbourhood's size (",sizeN,") with number of variables (",nVars,").",sep=""))
- else return(TRUE)
+ if (object at prediction & is.null(object at dataLocs))
+ return("The locations of the data have to provided for the estimation procedure.")
+ if (nrow(object at data) != nrow(object at distances))
+ return("Data and distances have unequal number of rows.")
+ if (ncol(object at data) %% (sizeN-object at prediction) != 0)
+ return("Data and distances have non matching number of columns.")
+ if (nrow(object at data) != nrow(object at index))
+ return("Data and index have unequal number of rows.")
+ if (ncol(object at distances) != ncol(object at index))
+ return("Data and index have unequal number of columns.")
+ if (ncol(object at data) != (sizeN-object at prediction) * nVars)
+ return(paste("Number of columns in data does not equal the product of the neighbourhood's size (",sizeN,") with number of variables (",nVars,").",sep=""))
+ else
+ return(TRUE)
}
+setClassUnion("optionalDataLocs",c("NULL","Spatial"))
+
setClass("neighbourhood",
representation = representation(data = "data.frame",
distances="matrix",
index="matrix",
- locations="Spatial",
+ locations="Spatial",
+ dataLocs="optionalDataLocs",
var="character",
prediction="logical"),
validity = validNeighbourhood, contains = list("Spatial"))
Modified: pkg/R/ClaytonGumbelCopula.R
===================================================================
--- pkg/R/ClaytonGumbelCopula.R 2013-04-23 14:17:11 UTC (rev 94)
+++ pkg/R/ClaytonGumbelCopula.R 2013-05-03 15:23:07 UTC (rev 95)
@@ -30,7 +30,7 @@
)
# constructor
-surClaytonCopula <- function (param) {
+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.")
@@ -54,14 +54,14 @@
# ddu
setMethod("dduCopula", signature("numeric","surClaytonCopula"),
function(u, copula, log) {
- linkVineCop.ddu(matrix(u,ncol=copula at dimension),copula, log)
+ linkVineCop.ddu(matrix(u,ncol=copula at dimension),copula)
})
setMethod("dduCopula", signature("matrix","surClaytonCopula"), linkVineCop.ddu)
# ddv
setMethod("ddvCopula", signature("numeric","surClaytonCopula"),
function(u, copula, log) {
- linkVineCop.ddv(matrix(u,ncol=copula at dimension),copula, log)
+ linkVineCop.ddv(matrix(u,ncol=copula at dimension),copula)
})
setMethod("ddvCopula", signature("matrix","surClaytonCopula"), linkVineCop.ddv)
@@ -104,7 +104,7 @@
)
# constructor
-r90ClaytonCopula <- function (param) {
+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.")
@@ -164,7 +164,7 @@
)
# constructor
-r270ClaytonCopula <- function (param) {
+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.")
@@ -245,7 +245,7 @@
)
# constructor
-surGumbelCopula <- function (param) {
+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.")
@@ -320,7 +320,7 @@
)
# constructor
-r90GumbelCopula <- function (param) {
+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.")
@@ -380,7 +380,7 @@
)
# constructor
-r270GumbelCopula <- function (param) {
+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.")
Modified: pkg/R/asCopula.R
===================================================================
--- pkg/R/asCopula.R 2013-04-23 14:17:11 UTC (rev 94)
+++ pkg/R/asCopula.R 2013-05-03 15:23:07 UTC (rev 95)
@@ -6,7 +6,7 @@
# (see Example 3.16 in: Nelsen, Roger B. (2006): An Introduction to Copulas, second edition, Springer)
# constructor
-asCopula <- function (param) {
+asCopula <- function (param=c(0,0)) {
val <- new("asCopula", dimension = as.integer(2), parameters = param,
param.names = c("a", "b"), param.lowbnd = c(limA(param[2]), -1),
param.upbnd = c(1, 1), fullname = "asymmetric copula family with cubic and quadratic sections")
Modified: pkg/R/cqsCopula.R
===================================================================
--- pkg/R/cqsCopula.R 2013-04-23 14:17:11 UTC (rev 94)
+++ pkg/R/cqsCopula.R 2013-05-03 15:23:07 UTC (rev 95)
@@ -1,3 +1,6 @@
+## make fitCopula generic
+setGeneric("fitCopula")
+
######################################################
## ##
## a symmetric copula with cubic quadratic sections ##
@@ -4,16 +7,15 @@
## ##
######################################################
-cqsCopula <- function (param) {
- val <- new("cqsCopula", dimension = as.integer(2), parameters = param,
+cqsCopula <- function (param=c(0,0), fixed=character(0)) {
+ new("cqsCopula", dimension = as.integer(2), parameters = param,
param.names = c("a", "b"), param.lowbnd = c(limA(param[2]),-1),
- param.upbnd = c(1, 1), fullname = "copula family with cubic quadratic sections")
- val
+ param.upbnd = c(1, 1),
+ fullname = "copula family with cubic quadratic sections", fixed=fixed)
}
## density ##
-
-dCQSec <- function (u, copula) {
+dCQSec <- function (u, copula, log=F) {
a <- copula at parameters[1]
b <- copula at parameters[2]
@@ -22,17 +24,19 @@
u1 <- u[, 1]
u2 <- u[, 2]
- return(pmax(1-b*(1-2*u2)*(1-2*u1)+(b-a)*(1-u2)*(1-3*u2)*(1-u1)*(1-3*u1),0))
+ if (log)
+ return(log(pmax(1-b*(1-2*u2)*(1-2*u1)+(b-a)*(1-u2)*(1-3*u2)*(1-u1)*(1-3*u1),0)))
+ else
+ return(pmax(1-b*(1-2*u2)*(1-2*u1)+(b-a)*(1-u2)*(1-3*u2)*(1-u1)*(1-3*u1),0))
}
setMethod("dCopula", signature("numeric", "cqsCopula"),
- function(u, copula, ...) {
- dCQSec(matrix(u,ncol=copula at dimension), copula)
+ function(u, copula, log) {
+ dCQSec(matrix(u,ncol=copula at dimension), copula, log)
})
setMethod("dCopula", signature("matrix", "cqsCopula"), dCQSec)
## jcdf ##
-
pCQSec <- function (u, copula) {
a <- copula at parameters[1]
b <- copula at parameters[2]
@@ -42,6 +46,7 @@
u2 <- u[, 2]
return(u1*u2*(1- b*(1-u1)*(1-u2) + (b-a)*(1-u2)^2*(1-u1)^2))
}
+
setMethod("pCopula", signature("numeric", "cqsCopula"),
function(u, copula, ...) {
pCQSec(matrix(u,ncol=copula at dimension), copula)
@@ -55,7 +60,7 @@
solveCubicEq <- function(a,b,c,d){
eps <- .Machine$double.eps
-# using the reduced equation z^3 + 3 * p * z + q = 0 with:
+ # using the reduced equation z^3 + 3 * p * z + q = 0 with:
p <- 3*a*c-b^2
q <- 2*b^3-9*a*b*c+27*a^2*d
D <- q^2+4*p^3
@@ -113,42 +118,43 @@
## inverse partial derivative ddu
# seems to be accurate (1.4e-05 is the max out of 1000 random CQSec-copulas for 1000 random pairs (u,v) each.)
invdduCQSec <- function (u, copula, y) {
- if (length(u)!=length(y))
- stop("Length of u and y differ!")
+ stopifnot(length(u)==length(y))
+
+ a <- copula at parameters[1]
+ b <- copula at parameters[2]
- a <- copula at parameters[1]
- b <- copula at parameters[2]
+ # solving the cubic equation: u^3 * c3 + u^2 * c2 + u * c1 + c0 = 0
+ usq <- u^2
+ c3 <- (b-a)*(1-4*u+3*usq)
+ c2 <- (b-a)*(-2+8*u-6*u^2)-b*(-1+2*u)
+ c1 <- (b-a)*(1-4*u+3*u^2)-b*(1-2*u)+1
+ c0 <- -y
-# solving the cubic equation: u^3 * c3 + u^2 * c2 + u * c1 + c0 = 0
- usq <- u^2
- c3 <- (b-a)*(1-4*u+3*usq)
- c2 <- (b-a)*(-2+8*u-6*u^2)-b*(-1+2*u)
- c1 <- (b-a)*(1-4*u+3*u^2)-b*(1-2*u)+1
- c0 <- -y
+ v <- solveCubicEq(c3,c2,c1,c0)
-v <- solveCubicEq(c3,c2,c1,c0)
-
-filter <- function(vec){
- vec <- vec[!is.na(vec)]
- return(vec[vec >= 0 & vec <= 1])
+ return(v)
+
+# filter <- function(vec){
+# vec <- vec[!is.na(vec)]
+# return(vec[vec >= 0 & vec <= 1])
+# }
+#
+# return(apply(v,1,filter))
}
-return(apply(v,1,filter))
-}
-
setMethod("invdduCopula", signature("numeric","cqsCopula","numeric"), invdduCQSec)
## partial derivative ddv
ddvCQSec <- function (u, copula) {
- a <- copula at parameters[1]
- b <- copula at parameters[2]
- if (!is.matrix(u)) u <- matrix(u, ncol = 2)
+ a <- copula at parameters[1]
+ b <- copula at parameters[2]
+ if (!is.matrix(u)) u <- matrix(u, ncol = 2)
- u1 <- u[, 1]
- u2 <- u[, 2]
+ u1 <- u[, 1]
+ u2 <- u[, 2]
-return(u1-b*(u1-2*u1*u2-u1^2+2*u1^2*u2)+(b-a)*(u1-2*u1^2+u1^3-4*u1*u2+8*u1^2*u2-4*u1^3*u2+3*u1*u2^2-6*u1^2*u2^2+3*u1^3*u2^2))
+ return(u1-b*(u1-2*u1*u2-u1^2+2*u1^2*u2)+(b-a)*(u1-2*u1^2+u1^3-4*u1*u2+8*u1^2*u2-4*u1^3*u2+3*u1*u2^2-6*u1^2*u2^2+3*u1^3*u2^2))
}
setMethod("ddvCopula", signature("numeric","cqsCopula"),
@@ -160,29 +166,30 @@
## inverse partial derivative ddv
# seems to be accurate (1e-05 is the max out of 5000 random CQSec-copulas for 1000 random pairs (u,v) each. Very most are below 10*.Machine$double.eps)
invddvCQSec <- function (v, copula, y) {
- if (length(v)!=length(y))
- stop("Length of v and y differ!")
+ stopifnot(length(v)==length(y))
- a <- copula at parameters[1]
- b <- copula at parameters[2]
+ a <- copula at parameters[1]
+ b <- copula at parameters[2]
-# solving the cubic equation: u^3 * c3 + u^2 * c2 + u * c1 + c0 = 0
- vsq <- v^2
- c3 <- (b-a)*(1-4*v+3*vsq)
- c2 <- (b-a)*(-2+8*v-6*vsq)-b*(-1+2*v)
- c1 <- (b-a)*(1-4*v+3*vsq)-b*(1-2*v)+1
- c0 <- -y
+ # solving the cubic equation: u^3 * c3 + u^2 * c2 + u * c1 + c0 = 0
+ vsq <- v^2
+ c3 <- (b-a)*(1-4*v+3*vsq)
+ c2 <- (b-a)*(-2+8*v-6*vsq)-b*(-1+2*v)
+ c1 <- (b-a)*(1-4*v+3*vsq)-b*(1-2*v)+1
+ c0 <- -y
-u <- solveCubicEq(c3,c2,c1,c0)
+ u <- solveCubicEq(c3,c2,c1,c0)
+
+ return(u)
-filter <- function(vec){
- vec <- vec[!is.na(vec)]
- return(vec[vec >= 0 & vec <= 1])
+# filter <- function(vec){
+# vec <- vec[!is.na(vec)]
+# return(vec[vec >= 0 & vec <= 1])
+# }
+#
+# return(apply(u,1,filter))
}
-return(apply(u,1,filter))
-}
-
setMethod("invddvCopula", signature("numeric","cqsCopula","numeric"), invddvCQSec)
## random number generator
@@ -227,62 +234,73 @@
# one of kendall or spearman according to the calculation of moa
fitCQSec.itau <- function(copula, data, estimate.variance, tau=NULL) {
-if(is.null(tau))
- tau <- VineCopula:::fasttau(data[,1],data[,2])
-esti <- fitCQSec.moa(tau, data, method="itau")
-copula <- cqsCopula(esti)
-return(new("fitCopula",
- estimate = esti,
- var.est = matrix(NA),
- method = "Inversion of Kendall's tau and MLE",
- loglik = sum(log(dCopula(data, copula))),
- fitting.stats=list(convergence = as.integer(NA)),
- nsample = nrow(data),
- copula=copula
-))
-}
+ if(is.null(tau))
+ tau <- VineCopula:::fasttau(data[,1],data[,2])
+ if(copula at fixed=="a")
+ esti <- c(copula at parameters[1], iTauCQSec.a(copula at parameters[1], tau))
+ if(copula at fixed=="b")
+ esti <- c(iTauCQSec.b(copula at parameters[2], tau),copula at parameters[2])
+ else
+ esti <- fitCQSec.moa(tau, data, method="itau")
+
+ copula <- cqsCopula(esti, fixed=copula at fixed)
-fitCQSec.irho <- function(copula, data, estimate.variance){
-rho <- cor(data,method="spearman")[1,2]
-esti <- fitCQSec.moa(rho, data, method="irho")
-copula <- cqsCopula(esti)
-return(new("fitCopula",
- estimate = esti,
- var.est = matrix(NA),
- method = "Inversion of Spearman's rho and MLE",
- loglik = sum(log(dCopula(data, copula))),
- fitting.stats=list(convergence = as.integer(NA)),
- nsample = nrow(data),
- copula=copula
-))
+ new("fitCopula", estimate = esti, var.est = matrix(NA),
+ method = "Inversion of Kendall's tau and MLE",
+ loglik = sum(dCopula(data, copula, log=T)),
+ fitting.stats=list(convergence = as.integer(NA)), nsample = nrow(data),
+ copula=copula)
}
-fitCQSec.moa <- function(moa, data, method="itau", tol=.Machine$double.eps^.5) {
-smpl <- as.matrix(data)
-iTau <- function(p) {
- iTauCQSec(p,moa)
-}
+fitCQSec.irho <- function(copula, data, estimate.variance, rho=NULL){
+ if(is.null(rho))
+ rho <- cor(data,method="spearman")[1,2]
+ if(copula at fixed=="a")
+ esti <- c(copula at parameters[1], iRhoCQSec.a(copula at parameters[1],rho))
+ if(copula at fixed=="b")
+ esti <- c(iRhoCQSec.b(copula at parameters[2],rho),copula at parameters[2])
+ else
+ esti <- fitCQSec.moa(rho, data, method="irho")
+
+ copula <- cqsCopula(esti, fixed=copula at fixed)
-iRho <- function(p) {
- iRhoCQSec(p,moa)
+ new("fitCopula", estimate = esti, var.est = matrix(NA),
+ method = "Inversion of Spearman's rho and MLE",
+ loglik = sum(dCopula(data, copula, log=T)),
+ fitting.stats=list(convergence = as.integer(NA)), nsample = nrow(data),
+ copula=copula)
}
-iFun <- switch(method, itau=iTau, irho=iRho)
-sec <- function (parameters) {
-res <- NULL
-for(param in parameters) {
- res <- rbind(res, -sum(log( dCQSec(smpl, cqsCopula(c(iFun(param),param))) )))
-}
-return(res)
-}
-b <- optimize(sec,c(-1,1), tol=tol)$minimum
+fitCQSec.moa <- function(moa, data, method="itau", tol=.Machine$double.eps^.5) {
+ smpl <- as.matrix(data)
-param <- c(iFun(b),b)
+ iRho.CQS <- function(p) {
+ iRhoCQSec.b(p,moa)
+ }
+
+ iTau.CQS <- function(p) {
+ iTauCQSec.b(p,moa)
+ }
+
+ iFun <- switch(method, itau=iTau.CQS, irho=iRho.CQS)
-return(param)
+ sec <- function (parameters) {
+ res <- NULL
+ for(param in parameters) {
+ res <- rbind(res, -sum(dCQSec(smpl, cqsCopula(c(iFun(param),param)),log=T)))
+ }
+
+ return(res)
+ }
+
+ b <- optimize(sec,c(-1,1), tol=tol)$minimum
+
+ param <- c(iFun(b),b)
+
+ return(param)
}
# maximum log-likelihood estimation of a and b using optim
@@ -290,26 +308,60 @@
fitCQSec.ml <- function(copula, data, start, lower, upper, optim.control, optim.method) {
if(length(start)!=2) stop("Start values need to have same length as parameters.")
- optFun <- function(param=c(0,0)) {
- if(any(param > 1) | param[2] < -1 | param[1] < limA(param[2])) return(1)
- return(-sum(log( dCQSec(data, cqsCopula(param)) )))
- }
+ if (length(copula at fixed)==0) {
+ optFun <- function(param=c(0,0)) {
+ if(any(param > 1) | param[2] < -1 | param[1] < limA(param[2]))
+ return(100)
+ return(-sum(log( dCQSec(data, cqsCopula(param)) )))
+ }
- optimized <- optim(par=start, fn=optFun, method = optim.method,
- lower=lower, upper=upper, control = optim.control)
+ optimized <- optim(par=start, fn=optFun, method = optim.method,
+ lower=lower, upper=upper, control = optim.control)
- return(new("fitCopula", estimate = optimized$par, var.est = matrix(NA),
- method = "Numerical MLE over the full range.",
- loglik = -optimized$value, fitting.stats= optimized,
- nsample = nrow(data), copula=cqsCopula(optimized$par)))
+ return(new("fitCopula", estimate = optimized$par, var.est = matrix(NA),
+ method = "Numerical MLE over the full range.",
+ loglik = -optimized$value, fitting.stats= optimized,
+ nsample = nrow(data), copula=cqsCopula(optimized$par)))
+ } else {
+ optFunFixed <- function(p) {
+ param <- switch(copula at fixed, a=c(copula at parameters[1],p),
+ b=c(p,copula at parameters[2]))
+ if(any(param > 1) | param[2] < -1 | param[1] < limA(param[2]))
+ return(100)
+ return(-sum(log( dCQSec(data, cqsCopula(param)) )))
+ }
+
+ optimized <- optimise(optFunFixed, c(-3,1))
+
+ optPar <- switch(copula at fixed, a=c(copula at parameters[1],optimized$minimum),
+ b=c(optimized$minimum,copula at parameters[2]))
+
+ return(new("fitCopula", estimate = optimized$minimum, var.est = matrix(NA),
+ method = "Numerical MLE over the full range.",
+ loglik = -optimized$objective, fitting.stats=list(),
+ nsample = nrow(data), copula=cqsCopula(optPar)))
+ }
}
####
-iTauCQSec <- function(b,tau=0) {
-return(min(max(limA(b),(b^2 + 75*b + 450*tau)/(b - 25)),1))
+iTauCQSec.a <- function(a, tau=0) {
+ limits <- limB(a)
+ min(max(limits[1],0.5*(sqrt(a^2-250*a-1800*tau+5626)+a-75)),limits[2])
}
+iTauCQSec.b <- function(b,tau=0) {
+ min(max(limA(b),(b^2 + 75*b + 450*tau)/(b - 25)),1)
+}
+
+setMethod("iTau",signature="cqsCopula",
+ function(copula, tau) {
+ switch(copula at fixed,
+ a=c(copula at parameters[1],iTauCQSec.a(copula at parameters[1],tau)),
+ b=c(iTauCQSec.b(copula at parameters[2],tau),copula at parameters[2]),
+ stop("iTau may only be used for cqsCopula with one parameter fixed."))
+ })
+
####
tauCQSec <- function(copula){
@@ -325,10 +377,23 @@
# find parameter "a" for parameter "b" under a given measure of association "rho"
# it may return a value exceeding the limit of "a" which may result in an invalid copula.
-iRhoCQSec <- function(b, rho=0) {
- return(min(max(limA(b),-3*b - 12*rho),1))
+iRhoCQSec.a <- function(a, rho=0) {
+ limits <- limB(a)
+ min(max(limits[1],-a/3 - 4*rho),limits[2])
}
+iRhoCQSec.b <- function(b, rho=0) {
+ min(max(limA(b),-3*b - 12*rho),1)
+}
+
+setMethod("iRho",signature="cqsCopula",
+ function(copula, rho) {
+ switch(copula at fixed,
+ a=function(copula, rho) c(a,iRhoCQSec.a(copula at parameters[1],rho)),
+ b=function(copula, rho) c(iRhoCQSec.b(copula at parameters[2],rho),b),
+ stop("iRho may only be used for cqsCopula with one parameter fixed."))
+ })
+
####
rhoCQSec <- function(copula){
Modified: pkg/R/joeBiCopula.R
===================================================================
--- pkg/R/joeBiCopula.R 2013-04-23 14:17:11 UTC (rev 94)
+++ pkg/R/joeBiCopula.R 2013-05-03 15:23:07 UTC (rev 95)
@@ -25,7 +25,7 @@
)
# constructor
-joeBiCopula <- function (param) {
+joeBiCopula <- function (param=2) {
if (any(is.na(param) | param >= Inf | param <= 1 ))
stop("Parameter is outside of the allowed interval (1,Inf).")
new("joeBiCopula", dimension = as.integer(2), parameters = param, param.names = c("theta"),
@@ -104,7 +104,7 @@
)
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/spcopula -r 95
More information about the spcopula-commits
mailing list