[spcopula-commits] r160 - pkg/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Dec 1 13:46:52 CET 2016
Author: ben_graeler
Date: 2016-12-01 13:46:52 +0100 (Thu, 01 Dec 2016)
New Revision: 160
Modified:
pkg/R/partialDerivatives.R
pkg/R/tawn3pCopula.R
Log:
fix error in partial derivatives of frank copula
Modified: pkg/R/partialDerivatives.R
===================================================================
--- pkg/R/partialDerivatives.R 2016-09-14 10:58:12 UTC (rev 159)
+++ pkg/R/partialDerivatives.R 2016-12-01 12:46:52 UTC (rev 160)
@@ -265,13 +265,20 @@
## partial derivative d/du
##########################
+## Wolfram alpha:
+# -e^α/(e^(α u) - e^α) - ((e^α - 1) e^(α + α u))/((e^(α u) - e^α) (e^α - e^(α + α u) + e^(α u + α v) - e^(α + α v)))
+
dduFrank <- function(u, copula){
rho <- copula at parameters
- u1 <- u[,1]
- u2 <- u[,2]
-
- exp(-rho*u1)*(exp(-rho*u2)-1) / ( (exp(-rho)-1) + (exp(-rho*u1)-1)*(exp(-rho*u2)-1) )
+ v <- u[,2]
+ u <- u[,1]
+ E <- exp(1)
+ eRho <- E^rho
+ eRhoU <- E^(rho * u)
+ eRhoV <- E^(rho * v)
+
+ -eRho/(eRhoU - eRho) - ((eRho - 1) * eRho * eRhoU)/((eRhoU - eRho) * (eRho - eRho * eRhoU + eRhoU * eRhoV - eRho * eRhoV))
}
setMethod("dduCopula", signature("numeric","frankCopula"),
@@ -280,30 +287,20 @@
})
setMethod("dduCopula", signature("matrix","frankCopula"), dduFrank)
-
-## inverse of the partial derivative d/du
-#########################################
-
-invdduFrank <- function(u, copula, y){
- rho <- copula at parameters[1]
- if (length(u)!=length(y))
- stop("Length of u and y differ!")
- return( (-1/rho) * log( y*( exp(-rho)-1)/(exp(-rho*u)-y*(exp(-rho*u)-1)) +1) ) # by DL
-}
-
-setMethod("invdduCopula", signature("numeric", "frankCopula", "numeric"), invdduFrank)
-
-
## partial derivative d/dv
##########################
ddvFrank <- function(u, copula){
rho <- copula at parameters
-
- u1 <- u[,1]
- u2 <- u[,2]
-
- exp(-rho*u2)*(exp(-rho*u1)-1) / ( (exp(-rho)-1) + (exp(-rho*u2)-1)*(exp(-rho*u1)-1) )
+
+ v <- u[,2]
+ u <- u[,1]
+ E <- exp(1)
+ eRho <- E^rho
+ eRhoU <- E^(rho * u)
+ eRhoV <- E^(rho * v)
+
+ -eRho/(eRhoV - eRho) - ((eRho - 1) * eRho * eRhoV)/((eRhoV - eRho) * (eRho - eRho * eRhoV + eRhoV * eRhoU - eRho * eRhoU))
}
setMethod("ddvCopula", signature("numeric","frankCopula"),
@@ -312,19 +309,6 @@
})
setMethod("ddvCopula", signature("matrix","frankCopula"), ddvFrank)
-
-## inverse of the partial derivative d/dv
-#########################################
-
-invddvFrank <- function(v, copula, y){
- rho <- copula at parameters[1]
- if (length(v)!=length(y))
- stop("Length of v and y differ!")
- return( (-1/rho) * log( y*( exp(-rho)-1)/(exp(-rho*v)-y*(exp(-rho*v)-1)) +1) )
-}
-
-setMethod("invddvCopula", signature("numeric", "frankCopula", "numeric"), invddvFrank)
-
####################
## student Copula ##
####################
Modified: pkg/R/tawn3pCopula.R
===================================================================
--- pkg/R/tawn3pCopula.R 2016-09-14 10:58:12 UTC (rev 159)
+++ pkg/R/tawn3pCopula.R 2016-12-01 12:46:52 UTC (rev 160)
@@ -1,144 +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")
- 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
+#######################################
+## 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