[spcopula-commits] r156 - in pkg: . R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Sep 1 15:33:20 CEST 2016
Author: ben_graeler
Date: 2016-09-01 15:33:20 +0200 (Thu, 01 Sep 2016)
New Revision: 156
Added:
pkg/R/mixtureCopula.R
Modified:
pkg/DESCRIPTION
pkg/R/asCopula.R
pkg/R/cqsCopula.R
pkg/R/empiricalCopula.R
pkg/R/tawn3pCopula.R
Log:
adds mixtureCopula
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2016-08-30 07:33:12 UTC (rev 155)
+++ pkg/DESCRIPTION 2016-09-01 13:33:20 UTC (rev 156)
@@ -2,7 +2,7 @@
Type: Package
Title: Copula Driven Spatio-Temporal Analysis
Version: 0.2-1
-Date: 2016-08-30
+Date: 2016-09-01
Authors at R: c(person("Benedikt", "Graeler", role = c("aut", "cre"),
email = "ben.graeler at uni-muenster.de"),
person("Marius", "Appel",role = "ctb"))
@@ -21,6 +21,7 @@
cqsCopula.R
asCopula.R
leafCopula.R
+ mixtureCopula.R
spCopula.R
stCopula.R
spatialPreparation.R
Modified: pkg/R/asCopula.R
===================================================================
--- pkg/R/asCopula.R 2016-08-30 07:33:12 UTC (rev 155)
+++ pkg/R/asCopula.R 2016-09-01 13:33:20 UTC (rev 156)
@@ -307,4 +307,7 @@
return((a+3*b)/12)
}
-setMethod("rho",signature("asCopula"), rhoASC2)
\ No newline at end of file
+setMethod("rho", signature("asCopula"), rhoASC2)
+
+setMethod("lambda", signature("asCopula"),
+ function(copula, ...) c(lower = 0, upper = 0))
\ No newline at end of file
Modified: pkg/R/cqsCopula.R
===================================================================
--- pkg/R/cqsCopula.R 2016-08-30 07:33:12 UTC (rev 155)
+++ pkg/R/cqsCopula.R 2016-09-01 13:33:20 UTC (rev 156)
@@ -416,4 +416,7 @@
return( -(a+3*b)/12 )
}
-setMethod("rho",signature("cqsCopula"),rhoCQSec)
\ No newline at end of file
+setMethod("rho",signature("cqsCopula"),rhoCQSec)
+
+setMethod("lambda", signature("cqsCopula"),
+ function(copula, ...) c(lower = 0, upper = 0))
\ No newline at end of file
Modified: pkg/R/empiricalCopula.R
===================================================================
--- pkg/R/empiricalCopula.R 2016-08-30 07:33:12 UTC (rev 155)
+++ pkg/R/empiricalCopula.R 2016-09-01 13:33:20 UTC (rev 156)
@@ -26,7 +26,7 @@
# 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)
+ empiricalCopula(rCopula(sample.size, copula), copula)
}
@@ -35,10 +35,6 @@
## jcdf ##
# from package copula
pempCop.C <- function(u, copula) {
- # r-forge hack, to be removed after release of copula 0.999-6
- if(exists("C.n")) {
- return(C.n(u, copula at sample, offset=0, method="C"))
- } else
return(Cn(copula at sample,u))
}
@@ -53,17 +49,19 @@
TauMatrix(copula at sample)[1,2]
}
-setMethod("tau",signature("asCopula"),tauempCop)
+setMethod("tau",signature("empiricalCopula"), tauempCop)
rhoempCop <- function(copula){
cor(copula at sample,method="spearman")
}
-setMethod("rho",signature("asCopula"),rhoempCop)
+setMethod("rho",signature("empiricalCopula"), rhoempCop)
+setMethod("lambda", signature("empiricalCopula"),
+ function(copula, ...) stop("No evaluation possible, try to plot 'empBivJointDepFun' for a visual assessment."))
-# Vine Copula
+# Vine Copula - empirical evaluation
## jcdf ##
pvineCopula <- function(u, copula) {
empCop <- genEmpCop(copula, 1e5)
@@ -72,9 +70,9 @@
}
setMethod("pCopula", signature("numeric","vineCopula"),
- function(u,copula) {
- pvineCopula(matrix(u, ncol=copula at dimension),copula)
+ function(u, copula) {
+ pvineCopula(matrix(u, ncol=copula at dimension), copula)
})
setMethod("pCopula", signature("data.frame","vineCopula"),
- function(u,copula) pvineCopula(as.matrix(u),copula))
-setMethod("pCopula", signature("matrix","vineCopula"), pvineCopula)
+ function(u, copula) pvineCopula(as.matrix(u), copula))
+setMethod("pCopula", signature("matrix","vineCopula"), pvineCopula)
\ No newline at end of file
Added: pkg/R/mixtureCopula.R
===================================================================
--- pkg/R/mixtureCopula.R (rev 0)
+++ pkg/R/mixtureCopula.R 2016-09-01 13:33:20 UTC (rev 156)
@@ -0,0 +1,166 @@
+##############################
+## ##
+## a general mixture copula ##
+## ##
+##############################
+
+# class
+setClass("mixtureCopula", contains = "copula", slots = list(memberCops= "list"))
+
+# constructor
+mixtureCopula <- function (param = c(0.2, 0.2, 0.5), memberCops = c(normalCopula(), claytonCopula())) {
+ stopifnot(length(memberCops) == 2)
+ stopifnot(memberCops[[1]]@dimension == memberCops[[2]]@dimension)
+
+ cop1.nPar <- length(memberCops[[1]]@parameters)
+ cop2.nPar <- length(memberCops[[2]]@parameters)
+
+ if (missing(param))
+ param <- 0.5
+ if (length(param) == 1)
+ param <- c(memberCops[[1]]@parameters, memberCops[[2]]@parameters, 0.5)
+ else {
+ stopifnot(length(param) == cop1.nPar + cop2.nPar + 1)
+
+ memberCops[[1]]@parameters <- param[1:cop1.nPar]
+ memberCops[[2]]@parameters <- param[(1:cop2.nPar)+cop1.nPar]
+ }
+
+ new("mixtureCopula", dimension = memberCops[[1]]@dimension, parameters = param, memberCops = memberCops,
+ param.names = c(memberCops[[1]]@param.names, memberCops[[2]]@param.names, "mixLambda"),
+ param.lowbnd = c(memberCops[[1]]@param.lowbnd, memberCops[[2]]@param.lowbnd, 0),
+ param.upbnd = c(memberCops[[1]]@param.upbnd, memberCops[[2]]@param.upbnd, 1),
+ fullname = paste("mixture of a", memberCops[[1]]@fullname, "and a", memberCops[[2]]@fullname))
+}
+
+## density ##
+setMethod("dCopula", signature(copula = "mixtureCopula"),
+ function(u, copula, log, ...) {
+ mixLambda <- tail(copula at parameters, 1)
+ res <- (1-mixLambda) * dCopula(u, copula at memberCops[[1]], ...) + mixLambda * dCopula(u, copula at memberCops[[2]], ...)
+ if (log)
+ return(log(res))
+ else
+ return(res)
+ })
+
+## jcdf ##
+setMethod("pCopula", signature( copula = "mixtureCopula"),
+ function(u, copula, ...) {
+ mixLambda <- tail(copula at parameters, 1)
+ (1-mixLambda) * pCopula(u, copula at memberCops[[1]]) + mixLambda * pCopula(u, copula at memberCops[[2]])
+ })
+
+## partial derivatives ##
+## ddu
+
+setMethod("dduCopula", signature(copula = "mixtureCopula"),
+ function(u, copula, ...) {
+ mixLambda <- tail(copula at parameters, 1)
+ (1-mixLambda) * dduCopula(u, copula at memberCops[[1]]) + mixLambda * dduCopula(u, copula at memberCops[[2]])
+ })
+
+# ddv
+setMethod("ddvCopula", signature(copula = "mixtureCopula"),
+ function(u, copula, ...) {
+ mixLambda <- tail(copula at parameters, 1)
+ (1-mixLambda) * ddvCopula(u, copula at memberCops[[1]]) + mixLambda * ddvCopula(u, copula at memberCops[[2]])
+ })
+
+## inverse partial derivative
+# invddu
+invdduMixCop <- function (u, copula, y) {
+ stopifnot(length(u) == length(y))
+
+ opti <- function(ind) {
+ optFun <- function(v) {
+ (dduCopula(cbind(u[ind], v), copula) - y[ind])^2
+ }
+ optimise(optFun, c(0,1))$minimum
+ }
+
+ sapply(1:length(y), opti)
+}
+
+setMethod("invdduCopula",
+ signature("numeric", "mixtureCopula", "numeric"),
+ invdduMixCop)
+
+# invddv
+invddvMixCop <- function (v, copula, y) {
+ stopifnot(length(v) == length(y))
+
+ opti <- function(ind) {
+ optFun <- function(u) {
+ (dduCopula(cbind(u, v[ind]), copula) - y[ind])^2
+ }
+ optimise(optFun, c(0,1))$minimum
+ }
+
+ sapply(1:length(y), opti)
+}
+
+setMethod("invddvCopula",
+ signature("numeric", "mixtureCopula", "numeric"),
+ invddvMixCop)
+
+## random number generator
+
+rMixCop <- function(n, copula, ...) {
+ u <- runif(n)
+ y <- runif(n)
+
+ cbind(u, invdduCopula(u, copula, y))
+}
+
+setMethod("rCopula", signature(copula = "mixtureCopula"), rMixCop)
+
+## fitment
+fitMixCop <- function(copula, data, start, method="mpl",
+ lower = NULL, upper = NULL,
+ optim.method = "BFGS", optim.control = list(maxit = 1000),
+ estimate.variance = FALSE, ...){
+ if (missing(start))
+ start <- copula at parameters
+ stopifnot(method %in% c("ml", "mpl"))
+
+ if(is.null(lower))
+ lower <- copula at param.lowbnd
+ if(is.null(upper))
+ upper <- copula at param.lowbnd
+
+ copula:::fitCopula.ml(copula, data, start = start, method = method,
+ lower = lower, upper = upper,
+ optim.method = optim.method,
+ optim.control = optim.control,
+ estimate.variance = estimate.variance , ...)
+}
+
+setMethod(fitCopula,
+ signature = c(copula = "mixtureCopula"),
+ fitMixCop)
+
+mixCop <- mixtureCopula(c(0.2,0.5,0.3))
+fitCopula(mixCop, rCopula(300, mixCop))
+
+fitMixCop(mixCop, rCopula(300, mixCop))
+
+#
+# fitCopulaASC2 <- function (copula, data, method = "ml", start=c(0,0),
+# lower=c(-3,-1), upper=c(1,1),
+# optim.method="L-BFGS-B", optim.control=list(),
+# estimate.variance = FALSE) {
+# fit <- switch(method,
+# ml=fitASC2.ml(copula, data, start, lower, upper, optim.control, optim.method),
+# itau=fitASC2.itau(copula, data, estimate.variance),
+# irho=fitASC2.irho(copula, data, estimate.variance),
+# stop("Implemented methods for copulas in the spCopula package are: ml, itau, and irho."))
+# return(fit)
+# }
+#
+# setMethod("fitCopula", signature("asCopula"), fitCopulaASC2)
+
+# setMethod("tau",signature("asCopula"),tauASC2)
+# setMethod("rho", signature("asCopula"), rhoASC2)
+# setMethod("lambda", signature("asCopula"),
+# function(copula, ...) c(lower = 0, upper = 0))
\ No newline at end of file
Modified: pkg/R/tawn3pCopula.R
===================================================================
--- pkg/R/tawn3pCopula.R 2016-08-30 07:33:12 UTC (rev 155)
+++ pkg/R/tawn3pCopula.R 2016-09-01 13:33:20 UTC (rev 156)
@@ -1,109 +1,144 @@
-#######################################
-## tawn copula with all 3 parameters ##
-#######################################
-
-setClass("tawn3pCopula", representation(exprdist = "expression"),
- contains = "evCopula")
-
-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)
-
-}
-
-ATawn <- function(copula, w) {
- Atawn3p(w,copula at parameters)
-}
-
-setMethod("A",signature("tawn3pCopula"),ATawn)
-
-dAduTawn <- function(copula, w) {
- alpha <- copula at parameters[1]
- beta <- copula at parameters[2]
- theta <- copula at parameters[3]
-
- # 1st derivative
- p1 <- (alpha*(alpha*(-(w-1)))^(theta-1)-beta*(beta*w)^(theta-1))
- p2 <- ((alpha*(-(w-1)))^theta+(beta*w)^theta)^(1/theta-1)
-
- # 2nd derivative
- p3 <- (alpha*(-(w-1)))^(theta-2)
- p4 <- (beta*w)^(theta-2)
- p5 <- ((alpha*(-(w-1)))^theta+(beta*w)^theta)^(1/theta-2)
-
- data.frame(der1=alpha-beta-p1*p2,
- der2=alpha^2*beta^2*(theta-1)*p3*p4*p5)
-}
-
-setMethod("dAdu",signature("tawn3pCopula"),dAduTawn)
-
-tawn3pCopula <- function (param = c(0.5, 0.5, 2)) {
- # A(t) = (1-beta)*t + (1-alpha)*(1-t) + ((alpha*(1-t))^theta+(beta*t)^theta)^(1/theta)
- # C(u1,u2) = exp(log(u1*u2) * A(log(u2)/log(u1*u2)))
- # = u1*u2 + exp(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")
-}
-
-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))
- ## 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)
-
-
-fitTawn3pCop <- function(copula, data, method = c("mpl", "ml", "itau", "irho"),
- start = copula at parameters,
- lower = copula at param.lowbnd,
- upper = copula at param.upbnd,
- optim.method = "L-BFGS-B",
- optim.control = list(maxit = 1000), estimate.variance = FALSE,
- hideWarnings = TRUE) {
-
- fitCopulaAny <- selectMethod(fitCopula, "copula")
- fitCopulaAny(copula, data, method, start, lower, upper,
- optim.method, optim.control, estimate.variance,
- hideWarnings)
-}
-
+#######################################
+## tawn copula with all 3 parameters ##
+#######################################
+
+setClass("tawn3pCopula", representation(exprdist = "expression"),
+ contains = "evCopula")
+
+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)
+
+}
+
+ATawn <- function(copula, w) {
+ Atawn3p(w, copula at parameters)
+}
+
+setMethod("A", signature("tawn3pCopula"), ATawn)
+
+dAduTawn <- function(copula, w) {
+ alpha <- copula at parameters[1]
+ beta <- copula at parameters[2]
+ theta <- copula at parameters[3]
+
+ # 1st derivative
+ p1 <- (alpha*(alpha*(-(w-1)))^(theta-1)-beta*(beta*w)^(theta-1))
+ p2 <- ((alpha*(-(w-1)))^theta+(beta*w)^theta)^(1/theta-1)
+
+ # 2nd derivative
+ p3 <- (alpha*(-(w-1)))^(theta-2)
+ p4 <- (beta*w)^(theta-2)
+ p5 <- ((alpha*(-(w-1)))^theta+(beta*w)^theta)^(1/theta-2)
+
+ data.frame(der1=alpha-beta-p1*p2,
+ der2=alpha^2*beta^2*(theta-1)*p3*p4*p5)
+}
+
+setMethod("dAdu", signature("tawn3pCopula"), dAduTawn)
+
+tawn3pCopula <- function (param = c(0.5, 0.5, 2)) {
+ # A(t) = (1-beta)*t + (1-alpha)*(1-t) + ((alpha*(1-t))^theta+(beta*t)^theta)^(1/theta)
+ # C(u1,u2) = exp(log(u1*u2) * A(log(u2)/log(u1*u2)))
+ # = u1*u2 + exp(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")
+ dCdU2 <- D(cdf, "u2")
+ pdf <- D(dCdU1, "u2")
+
+ new("tawn3pCopula", dimension = 2L, exprdist = c(cdf = cdf, pdf = pdf,
+ dCdU = dCdU1, dCdV = dCdU2),
+ 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")
+}
+
+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))
+ ## improve log-case
+ if(log)
+ return(log(val))
+ else
+ val
+}
+
+setMethod("dCopula", signature(copula = "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(copula = "tawn3pCopula"), ptawn3pCopula)
+
+# partial derivatives
+
+ddutawn3pCopula <- 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]
+
+ return(eval(copula at exprdist$dCdU))
+}
+
+setMethod("dduCopula", signature(copula = "tawn3pCopula"), ddutawn3pCopula)
+
+ddvtawn3pCopula <- 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]
+
+ return(eval(copula at exprdist$dCdV))
+}
+
+setMethod("ddvCopula", signature(copula = "tawn3pCopula"), ddvtawn3pCopula)
+
+# tawn3pCop <- tawn3pCopula()
+# dduCopula(cbind(runif(10), runif(10)), tawn3pCop)
+## fit
+
+fitTawn3pCop <- function(copula, data, method = c("mpl", "ml"),
+ start = copula at parameters,
+ lower = copula at param.lowbnd,
+ upper = copula at param.upbnd,
+ optim.method = "L-BFGS-B",
+ optim.control = list(maxit = 1000), estimate.variance = FALSE,
+ hideWarnings = TRUE) {
+
+ fitCopulaAny <- selectMethod(fitCopula, "copula")
+ fitCopulaAny(copula, data, method, start, lower, upper,
+ optim.method, optim.control, estimate.variance,
+ hideWarnings)
+}
+
setMethod("fitCopula", signature("tawn3pCopula"), fitTawn3pCop)
\ No newline at end of file
More information about the spcopula-commits
mailing list