[spcopula-commits] r114 - pkg/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Nov 20 10:20:38 CET 2013
Author: ben_graeler
Date: 2013-11-20 10:20:38 +0100 (Wed, 20 Nov 2013)
New Revision: 114
Modified:
pkg/R/asCopula.R
pkg/R/cqsCopula.R
Log:
- taking care of special cases in asCopula and cqsCopula (i.e. equal parameters)
Modified: pkg/R/asCopula.R
===================================================================
--- pkg/R/asCopula.R 2013-11-19 15:49:54 UTC (rev 113)
+++ pkg/R/asCopula.R 2013-11-20 09:20:38 UTC (rev 114)
@@ -79,7 +79,7 @@
u1 <- u[, 1]
u2 <- u[, 2]
- return( u1 + b*(-1 + u1)*u1*(-1 + 2*u2) - (a - b)*(-1 + u1)^2*u1*u2*(-2 + 3*u2))
+ return( u1 + b*(u1-1)*u1*(2*u2-1) - (a - b)*(-1 + u1)^2*u1*u2*(-2 + 3*u2))
}
setMethod("ddvCopula", signature("numeric", "asCopula"),
@@ -94,56 +94,76 @@
## inverse partial derivative
invdduASC2 <- 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]
-# solving the cubic equation: u^3 * c3 + u^2 * c2 + u * c1 + c0 = 0
- usq <- u^2
- c3 <- (a-b)*(-3*usq+4*u-1)
- c2 <- (a-b)*(1-4*u+3*usq)+b*(- 1 + 2*u)
- c1 <- 1+b*(1-2*u)
- c0 <- -y
+ if (a != b) { # the cubic case
+ # solving the cubic equation: u^3 * c3 + u^2 * c2 + u * c1 + c0 = 0
+ usq <- u^2
+ c3 <- (a-b)*(-3*usq+4*u-1)
+ c2 <- (a-b)*(1-4*u+3*usq)+b*(- 1 + 2*u)
+ c1 <- 1+b*(1-2*u)
+ c0 <- -y
- v <- solveCubicEq(c3,c2,c1,c0) # from cqsCopula.R
-
- filter <- function(vec) {
- vec <- vec[!is.na(vec)]
- return(vec[vec >= 0 & vec <= 1])
- }
+ v <- solveCubicEq(c3,c2,c1,c0) # from cqsCopula.R
- return(apply(v,1,filter))
+ filter <- function(vec) {
+ vec <- vec[!is.na(vec)]
+ return(vec[vec >= 0 & vec <= 1])
+ }
+
+ return(apply(v,1,filter))
+ }
+ if(a==0) # and b==0 obvioulsy as well: the independent case
+ return(y)
+
+ # the qudratic cases remain
+ v <- y
+ uR <- u[u != 0.5]
+ v[u != 0.5] <- (sqrt(4*y[u != 0.5]*(2*b*uR-b)+(-2*b*uR+b+1)^2)+2*b*uR-b-1)/(2*b*(2*uR-1))
+
+ return(v)
}
setMethod("invdduCopula", signature("numeric","asCopula","numeric"),invdduASC2)
## inverse partial derivative ddv
invddvASC2 <- 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]
-
-# solving the cubic equation: u^3 * c3 + u^2 * c2 + u * c1 + c0 = 0
+ a <- copula at parameters[1]
+ b <- copula at parameters[2]
+
+ if (a != b) { # the cubic case
+ # solving the cubic equation: u^3 * c3 + u^2 * c2 + u * c1 + c0 = 0
vsq <- v^2
c3 <- (a-b)*(2*v-3*vsq)
c2 <- (a-b)*(-4*v+6*vsq)+b*(-1+2*v)
c1 <- 1+(a-b)*(2*v - 3*vsq)+b*(1-2*v)
c0 <- -y
-u <- solveCubicEq(c3,c2,c1,c0) # from cqsCopula.R
+ u <- solveCubicEq(c3,c2,c1,c0) # from cqsCopula.R
-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))
+ }
+ if(a==0) # and b==0 obvioulsy as well: the independent case
+ return(y)
+
+ # the qudratic cases remain
+ u <- y
+ vR <- v[v != 0.5]
+ u[v != 0.5] <- (sqrt(4*y[v != 0.5]*(2*b*vR-b)+(-2*b*vR+b+1)^2)+2*b*vR-b-1)/(2*b*(2*vR-1))
+
+ return(u)
}
-return(apply(u,1,filter))
-}
-
setMethod("invddvCopula", signature("numeric","asCopula","numeric"),invddvASC2)
## random number generator
Modified: pkg/R/cqsCopula.R
===================================================================
--- pkg/R/cqsCopula.R 2013-11-19 15:49:54 UTC (rev 113)
+++ pkg/R/cqsCopula.R 2013-11-20 09:20:38 UTC (rev 114)
@@ -118,26 +118,37 @@
## 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) {
- stopifnot(length(u)==length(y))
+ stopifnot(length(u) == length(y))
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
-
- v <- solveCubicEq(c3,c2,c1,c0)
+ if (a != b) { # the cubic case
+ # 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
- filter <- function(vec){
- vec <- vec[!is.na(vec)]
- return(vec[vec >= 0 & vec <= 1])
+ v <- solveCubicEq(c3,c2,c1,c0)
+
+ filter <- function(vec){
+ vec <- vec[!is.na(vec)]
+ return(vec[vec >= 0 & vec <= 1])
+ }
+
+ return(apply(v,1,filter))
}
-
- return(apply(v,1,filter))
+ if(a==0) # and b==0 obvioulsy as well: the independent case
+ return(y)
+
+ # the qudratic cases remain
+ v <- y
+ uR <- u[u != 0.5]
+ v[u != 0.5] <- (-sqrt((-2*b*uR+b-1)^2-4*y[u != 0.5]*(2*b*uR-b))+2*b*uR-b+1)/(2*b*(2*uR-1))
+
+ return(v)
}
setMethod("invdduCopula", signature("numeric","cqsCopula","numeric"), invdduCQSec)
@@ -162,28 +173,38 @@
setMethod("ddvCopula", signature("matrix","cqsCopula"), ddvCQSec)
## 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) {
stopifnot(length(v)==length(y))
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
-
- u <- solveCubicEq(c3,c2,c1,c0)
+ if (a != b) { # the cubic case
+ # 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
- filter <- function(vec){
- vec <- vec[!is.na(vec)]
- return(vec[vec >= 0 & vec <= 1])
+ u <- solveCubicEq(c3,c2,c1,c0)
+
+ filter <- function(vec){
+ vec <- vec[!is.na(vec)]
+ return(vec[vec >= 0 & vec <= 1])
+ }
+
+ return(apply(u,1,filter))
}
-
- return(apply(u,1,filter))
+ if(a==0) # and b==0 obvioulsy as well: the independent case
+ return(y)
+
+ # the qudratic cases remain
+ u <- y
+ vR <- v[v != 0.5]
+ u[v != 0.5] <- (-sqrt((-2*b*vR+b-1)^2-4*y[v != 0.5]*(2*b*vR-b))+2*b*vR-b+1)/(2*b*(2*vR-1))
+
+ return(u)
}
setMethod("invddvCopula", signature("numeric","cqsCopula","numeric"), invddvCQSec)
More information about the spcopula-commits
mailing list