[spcopula-commits] r80 - / pkg pkg/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Jan 23 12:37:45 CET 2013
Author: ben_graeler
Date: 2013-01-23 12:37:45 +0100 (Wed, 23 Jan 2013)
New Revision: 80
Added:
pkg/R/leafCopula.R
Modified:
pkg/DESCRIPTION
pkg/NAMESPACE
pkg/R/Classes.R
spcopula_0.1-1.tar.gz
spcopula_0.1-1.zip
Log:
- developed and added a new copula family: leafCopula
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2013-01-09 10:35:02 UTC (rev 79)
+++ pkg/DESCRIPTION 2013-01-23 11:37:45 UTC (rev 80)
@@ -2,7 +2,7 @@
Type: Package
Title: copula driven spatial analysis
Version: 0.1-1
-Date: 2013-01-06
+Date: 2013-01-23
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.
@@ -16,6 +16,7 @@
partialDerivatives.R
cqsCopula.R
asCopula.R
+ leafCopula.R
spCopula.R
stCopula.R
spatialPreparation.R
Modified: pkg/NAMESPACE
===================================================================
--- pkg/NAMESPACE 2013-01-09 10:35:02 UTC (rev 79)
+++ pkg/NAMESPACE 2013-01-23 11:37:45 UTC (rev 80)
@@ -1,7 +1,7 @@
import(copula, spacetime, CDVine)
# constructor
-export(asCopula, cqsCopula)
+export(asCopula, cqsCopula, leafCopula)
export(BB1Copula, surBB1Copula, r90BB1Copula, r270BB1Copula)
export(BB6Copula, surBB6Copula, r90BB6Copula, r270BB6Copula)
export(BB7Copula, surBB7Copula, r90BB7Copula, r270BB7Copula)
@@ -37,7 +37,7 @@
export(criticalPair, criticalTriple)
## classes
-exportClasses(asCopula, cqsCopula, neighbourhood, empiricalCopula)
+exportClasses(asCopula, cqsCopula, neighbourhood, empiricalCopula, leafCopula)
# wrappers to CDVine
exportClasses(BB1Copula, surBB1Copula, r90BB1Copula, r270BB1Copula)
Modified: pkg/R/Classes.R
===================================================================
--- pkg/R/Classes.R 2013-01-09 10:35:02 UTC (rev 79)
+++ pkg/R/Classes.R 2013-01-23 11:37:45 UTC (rev 80)
@@ -53,6 +53,24 @@
contains = list("copula")
)
+####
+## the leaf copula
+
+validLeafCopula <- function(object) {
+ if (object at dimension != 2)
+ return("The leaf copula only supports two dimensions.")
+
+ if (any(is.na(object at parameters)))
+ return("Parameter value is \"NA\".")
+ else return (TRUE)
+}
+
+setClass("leafCopula",
+ representation = representation("copula"),
+ validity = validLeafCopula,
+ contains = list("copula")
+)
+
##
## the spatial copula
##
Added: pkg/R/leafCopula.R
===================================================================
--- pkg/R/leafCopula.R (rev 0)
+++ pkg/R/leafCopula.R 2013-01-23 11:37:45 UTC (rev 80)
@@ -0,0 +1,546 @@
+## constructor
+leafCopula <- function (param=c(1.446923, -1.722742)) {
+ val <- new("leafCopula", dimension = as.integer(2), parameters = param,
+ param.names = c("a3", "a2"), param.lowbnd = c(-Inf, -Inf),
+ param.upbnd = c(Inf, Inf), fullname = "leaf copula")
+ return(val)
+}
+
+# weak lower border, two-place parameter
+weakBorderPoly <- function(x, par) {
+ par[1]*x^3+par[2]*x^2+x
+}
+
+# ##
+# curve(weakBorderPoly,asp=1,ylim=c(0,1),xlim=c(0,1))
+# ##
+
+ddxweakBorderPoly <- function(x, par) {
+ 3*par[1]*x^2+2*par[2]*x+1
+}
+
+invWeakBor <- function(v, par) {
+ optFun <- function(u) {
+ abs(weakBorderPoly(u, par)-v)
+ }
+
+ optimise(optFun,c(0,1),tol=.Machine$double.eps^0.5)$minimum
+}
+
+# strong upper border, no parameter
+strongBorderPoly <- function(x) {
+ -x^3+x^2+x
+}
+
+# ##
+# curve(strongBorderPoly,add=T)
+# ##
+
+ddxstrongBorderPoly <- function(x) {
+ -3*x^2+2*x+1
+}
+
+invStrongBor <- function(v) {
+ optFun <- function(u) {
+ sqrt(mean((strongBorderPoly(u)-v)^2))
+ }
+
+ optimise(optFun,c(0,1),tol=.Machine$double.eps^0.5)$minimum
+}
+
+# precalculate ellipse parameters
+solveQ <- function(u) {
+ sqrt(0.5*(strongBorderPoly(u)-u)^2)
+}
+
+ddxsolveQ <- function(u) {
+ sBor <- strongBorderPoly(u)
+ 1/(2*sqrt(0.5*(sBor-u)^2))*(sBor-u)*(ddxstrongBorderPoly(u)-1)
+}
+
+# ## double check
+# curve(solveQ)
+# curve(ddxsolveQ, add=T)
+# abline(h=0)
+# abline(solveQ(0.5)-ddxsolveQ(0.5)*0.5,ddxsolveQ(0.5))
+# abline(solveQ(0.9)-ddxsolveQ(0.9)*0.9,ddxsolveQ(0.9))
+# abline(v=c(0.5,0.9),col="grey")
+# ##
+
+solveXb <- function(u, par) {
+ sqrt(2*(u-weakBorderPoly(u, par))^2)+solveQ(u)
+}
+
+ddxsolveXb <- function(u, par) {
+ wBor <- weakBorderPoly(u, par)
+ -(sqrt(2)*(u-wBor)*(ddxweakBorderPoly(u, par)-1))/sqrt((u-wBor)^2)+ddxsolveQ(u)
+}
+
+# ## double check
+# curve(solveXb)
+# curve(ddxsolveXb, add=T)
+# abline(h=0)
+# abline(solveXb(0.5)-ddxsolveXb(0.5)*0.5,ddxsolveXb(0.5))
+# abline(solveXb(0.9)-ddxsolveXb(0.9)*0.9,ddxsolveXb(0.9))
+# abline(v=c(0.5,0.9), col="grey")
+# ## okay
+
+###
+# ellipse(xb)=b-q:
+# b = -q/(sqrt((a^2-x^2)/a^2)-1)
+#
+# ellipse'(xb)=1:
+# b = (a^2 sqrt((a^2-x^2)/a^2))/x
+#
+# solve for a:
+# -q/(sqrt((a^2-x^2)/a^2)-1) = (a^2 sqrt((a^2-x^2)/a^2))/x
+###
+
+solveA <- function(u, par) {
+ xb <- solveXb(u, par)
+ q <- solveQ(u)
+ sqrt((-xb^3+2*q*xb^2-q^2*xb)/(-2*xb+2*q+xb))
+}
+
+ddusolveA <- function(u, par) {
+ xb <- solveXb(u, par)
+ q <- solveQ(u)
+ dduXb <- ddxsolveXb(u, par)
+ dduQ <- ddxsolveQ(u)
+ 1/(2*solveA(u, par)) * (-2*(q-xb)*(q*xb*(dduQ-3*dduXb)+q^2*dduXb+xb^2*dduXb))/(-2*q + xb)^2
+}
+
+# ## double check
+# curve(solveA)
+# curve(ddusolveA,add=T)
+# abline(h=0)
+# abline(solveA(0.5)-ddusolveA(0.5)*0.5,ddusolveA(0.5))
+# abline(solveA(0.9)-ddusolveA(0.9)*0.9,ddusolveA(0.9))
+# abline(v=c(0.5, 0.9), col="grey")
+# ## okay
+
+solveB <- function(u, par) {
+ a <- solveA(u, par)
+ xb <- solveXb(u, par)
+ a^2*sqrt(1-(xb/a)^2)/xb
+}
+
+ddusolveB <- function(u, par) {
+ a <- solveA(u, par)
+ xb <- solveXb(u, par)
+ dduA <- ddusolveA(u, par)
+ dduXb <- ddxsolveXb(u, par)
+ (2*a^2*xb*dduA - xb^3*dduA - a^3*dduXb)/(a*xb^2*sqrt(1 - xb^2/a^2))
+}
+
+# ## double check
+# curve(solveB)
+# curve(ddusolveB,add=T)
+# abline(h=0)
+# abline(solveB(0.5)-ddusolveB(0.5)*0.5,ddusolveB(0.5))
+# abline(solveB(0.9)-ddusolveB(0.9)*0.9,ddusolveB(0.9))
+# abline(v=c(0.5,0.9), col="grey")
+# ## okay
+
+ellipse <- function(v,a,b) {
+ b*sqrt(1-(v/a)^2)-b
+}
+
+# ## double check
+# ellipseInV <- function(v) ellipse(v,a,b)
+#
+# curve(ellipseInV,0,0.00001,asp=1)
+# abline(v=solveXb(0.9)+c(0,-solveQ(0.9)))
+# abline(solveXb(0.9)-solveQ(0.9),-1)
+# abline(h=0)
+# ## okay
+
+dduellipse <- function(v, a, b, dduA, dduB) {
+ dduB*sqrt(1-(v/a)^2)+b/(2*sqrt(1-(v/a)^2))*2*v^2/a^3*dduA-dduB
+}
+
+# ## double check
+# dduellipseInU <- function(u) dduellipse(0.5,solveA(u),solveB(u),ddusolveA(u),ddusolveB(u))
+# ellipseInU <- function(u) ellipse(0.5,solveA(u),solveB(u))
+#
+# curve(ellipseInU,0.5,1,ylim=c(-0.2,.5))
+# curve(dduellipseInU,0.5,1,add=T)
+# abline(h=0)
+# abline(ellipseInU(0.6)-dduellipseInU(0.6)*0.6,dduellipseInU(0.6))
+# abline(ellipseInU(0.9)-dduellipseInU(0.9)*0.9,dduellipseInU(0.9))
+# abline(v=c(0.6,0.9), col="grey")
+# ## okay
+
+ddvellipse <- function(v, a, b) {
+ -b*v/sqrt(1-(v/a)^2)/a^2
+}
+
+# ## double check
+# curve(ellipseInV,0,0.6,ylim=c(-.2,0))
+# curve(ddvellipse,0,0.6,asp=2,add=T)
+# abline(h=0)
+# abline(ellipseInV(0.2)-ddvellipse(0.2)*0.2,ddvellipse(0.2))
+# abline(ellipseInV(0.4)-ddvellipse(0.4)*0.4,ddvellipse(0.4))
+# abline(v=c(0.2,0.4), col="grey")
+# ## okay
+
+ddvuEllipse <- function(v, a, b, dduA, dduB) {
+ (2*a^2*dduA*b*v - a^3*dduB*v - dduA*b*v^3 + a*dduB*v^3)/(a^3*(a^2 - v^2)*sqrt(1 - v^2/a^2))
+}
+
+# ## double checking
+# curve(dduellipse,ylim=c(-1,.5))
+# curve(ddvuEllipse,add=T) #,0,0.6,ylim=c(-.2,0))
+# abline(h=0)
+# abline(dduellipse(0.2)-ddvuEllipse(0.2)*0.2,ddvuEllipse(0.2))
+# abline(dduellipse(0.4)-ddvuEllipse(0.4)*0.4,ddvuEllipse(0.4))
+# abline(v=c(0.2,0.4), col="grey")
+# ## okay
+
+ddvvEllipse <- function(v, a, b) {
+ -b/((a^2 - v^2)*sqrt(1 - v^2/a^2))
+}
+
+# ## double checking
+# curve(ddvellipse,ylim=c(-1,.5))
+# curve(ddvvEllipse,add=T) #,0,0.6,ylim=c(-.2,0))
+# abline(h=0)
+# abline(ddvellipse(0.2)-ddvvEllipse(0.2)*0.2,ddvvEllipse(0.2))
+# abline(ddvellipse(0.4)-ddvvEllipse(0.4)*0.4,ddvvEllipse(0.4))
+# abline(v=c(0.2,0.4), col="grey")
+# ## okay
+
+vNorm <- function(vo, a, b) {
+ (-a^2*b + 2*sqrt(0.5)*a^2*vo + 2*a*b*sqrt(0.25*a^2 + sqrt(0.5)*b*vo - 0.5*vo^2))/(a^2 + b^2)
+}
+
+dduvNorm <- function(vo, a, b, dduA, dduB) {
+ cRoot <- sqrt(-0.5*vo^2 + 0.25*a^2 + sqrt(0.5)*vo*b)
+
+ low <- a^2+b^2
+ high <- -a^2*b+sqrt(2)*a^2*vo+2*a*b*cRoot
+
+ dduLow <- 2*a*dduA + 2*b*dduB
+ dduHigh <- (-2*a*dduA*b - a^2*dduB + 2*sqrt(2)*vo*a*dduA
+ + 2*(dduA*b + a*dduB)*cRoot
+ + a*b/cRoot*(0.5*a*dduA+sqrt(0.5)*vo*dduB))
+
+ (low*dduHigh - high*dduLow)/low^2
+}
+
+# ## double check
+# vNormInU <- function(u) vNorm(.5,solveA(u),solveB(u))
+# dduvNormInU <- function(u) dduvNorm(.5,solveA(u),solveB(u),ddusolveA(u),ddusolveB(u))
+#
+# curve(vNormInU,0.6,0.85,ylim=c(0.5,.7))
+# curve(dduvNormInU,add=T)
+# abline(vNormInU(0.65)-dduvNormInU(0.65)*0.65,dduvNormInU(0.65))
+# abline(vNormInU(0.8)-dduvNormInU(0.8)*0.8,dduvNormInU(0.8))
+# abline(v=c(0.65,0.8), col="grey")
+# ## okay
+
+ddvvNorm <- function(vo, a, b) {
+ c <- sqrt(0.5)
+ (2*a^2*(c + (b^2*(0.5*b*c - 0.5*vo))/sqrt(a^2*b^2*(0.25*a^2 + (b*c - 0.5*vo)*vo))))/(a^2 + b^2)
+}
+
+# ## double check
+# curve(vNorm,0,0.5,ylim=c(0,1))
+# curve(ddvvNorm,add=T)
+# abline(vNorm(0.4)-ddvvNorm(0.4)*0.4,ddvvNorm(0.4))
+# abline(vNorm(0.1)-ddvvNorm(0.1)*0.1,ddvvNorm(0.1))
+# abline(v=c(0.1,0.4), col="grey")
+# ## okay
+
+ddvuvNorm <- function(vo, a, b, dduA, dduB) {
+ cRoot <- sqrt(-0.5*vo^2 + 0.25*a^2 + sqrt(0.5)*vo*b)
+
+ low <- a^2+b^2
+
+ dduLow <- 2*a*dduA + 2*b*dduB
+ ddvHigh <- sqrt(2)*a^2+(a*b*(-vo + sqrt(0.5)*b))/cRoot
+
+ ddvuHigh <- (2*sqrt(2)*a*dduA
+ + ((dduA*b + a*dduB)*(-vo+sqrt(0.5)*b)+sqrt(1/2)*a*b*dduB)/cRoot
+ - (a*b*(0.5*a*dduA+sqrt(0.5)*dduB*vo)*(-vo+sqrt(0.5)*b))/(2*cRoot^3))
+
+ # ddv (low*dduHigh - high*dduLow):
+ (low * ddvuHigh - ddvHigh*dduLow)/low^2
+}
+
+# ## double check
+# curve(dduvNorm,0,0.5,ylim=c(-2,1))
+# curve(ddvuvNorm,add=T)
+# abline(h=0)
+# abline(dduvNorm(0.2)-ddvuvNorm(0.2)*0.2,ddvuvNorm(0.2))
+# abline(dduvNorm(0.4)-ddvuvNorm(0.4)*0.4,ddvuvNorm(0.4))
+# abline(v=c(0.2,0.4), col="grey")
+# ## okay
+
+ddvvvNorm <- function(vo, a, b) {
+ (-2*a^4*b^4)/((a^2*b^2*(a^2 + 2*(2*sqrt(0.5)*b - vo)*vo))^(3/2))
+}
+
+# ## double check
+# curve(ddvvNorm,0,0.5)#,ylim=c(0,1))
+# curve(ddvvvNorm,add=T)
+# abline(ddvvNorm(0.4)-ddvvvNorm(0.4)*0.4,ddvvvNorm(0.4))
+# abline(ddvvNorm(0.1)-ddvvvNorm(0.1)*0.1,ddvvvNorm(0.1))
+# abline(v=c(0.1,0.4), col="grey")
+# ## okay
+
+yOrig <- function(vn, a, b) {
+ sqrt(0.5)*(vn+ellipse(vn, a, b))
+}
+
+dduyOrig <- function(vn, a, b, dduA, dduB) {
+ sqrt(0.5)*dduellipse(vn, a, b, dduA, dduB)
+}
+
+# ## double check
+# yOrigInU <- function(u) yOrig(.5,solveA(u),solveB(u))
+# dduyOrigInU <- function(u) dduyOrig(.5,solveA(u),solveB(u),ddusolveA(u),ddusolveB(u))
+#
+# curve(yOrigInU,0.5,1,ylim=c(0.2,.4))
+# curve(dduyOrigInU,add=T)
+# abline(yOrigInU(0.6)-dduyOrigInU(0.6)*0.6,dduyOrigInU(0.6))
+# abline(yOrigInU(0.9)-dduyOrigInU(0.9)*0.9,dduyOrigInU(0.9))
+# abline(v=c(0.6,0.9), col="grey")
+# ## okay
+
+ddvyOrig <- function(vn, a, b) {
+ sqrt(0.5)*(1+ddvellipse(vn, a, b))
+}
+
+# ## double check
+# curve(yOrig,0,0.5,ylim=c(0,1))
+# curve(ddvyOrig,add=T)
+# abline(yOrig(0.475)-ddvyOrig(0.475)*0.475,ddvyOrig(0.475))
+# abline(yOrig(0.1)-ddvyOrig(0.1)*0.1,ddvyOrig(0.1))
+# abline(v=c(0.1,0.475), col="grey")
+# ## okay
+
+ddvuyOrig <- function(vn, a, b, dduA, dduB) {
+ sqrt(0.5)*ddvuEllipse(vn, a, b, dduA, dduB)
+}
+
+# ## double check
+# curve(dduyOrig,0,0.5,ylim=c(-0.2,0.4))
+# curve(ddvuyOrig,add=T)
+# abline(dduyOrig(0.475)-ddvuyOrig(0.475)*0.475,ddvuyOrig(0.475))
+# abline(dduyOrig(0.1)-ddvuyOrig(0.1)*0.1,ddvuyOrig(0.1))
+# abline(v=c(0.1,0.475), col="grey")
+# ## okay
+
+ddvvyOrig <- function(vn, a, b) {
+ sqrt(0.5)*(ddvvEllipse(vn, a, b))
+}
+
+# ## double check
+# curve(ddvyOrig,0,0.5)#,ylim=c(-0.2,0.4))
+# # curve(ddvvyOrig,add=F)
+# abline(ddvyOrig(0.475)-ddvvyOrig(0.475)*0.475,ddvvyOrig(0.475))
+# abline(ddvyOrig(0.1)-ddvvyOrig(0.1)*0.1,ddvvyOrig(0.1))
+# abline(v=c(0.1,0.475), col="grey")
+# ## okay
+
+ellipseOrig <- function(vo, a, b) {
+ yOrig(vNorm(vo, a, b), a, b)
+}
+
+dduellipseOrig <- function(vo, a, b, dduA, dduB) {
+ ddvyOrig(vNorm(vo, a, b), a, b)*dduvNorm(vo, a, b, dduA, dduB)+dduyOrig(vNorm(vo, a, b), a, b, dduA, dduB)
+}
+
+# ## double check
+# ellipseOrigInU <- function(u) ellipseOrig(.5,solveA(u),solveB(u))
+# dduellipseOrigInU <- function(u) dduellipseOrig(.5,solveA(u),solveB(u),ddusolveA(u),ddusolveB(u))
+#
+# curve(ellipseOrigInU,0.6,0.85,ylim=c(0.2,0.4))
+# curve(dduellipseOrigInU,add=T)
+# abline(ellipseOrigInU(0.8)-dduellipseOrigInU(0.8)*0.8,dduellipseOrigInU(0.8))
+# abline(ellipseOrigInU(0.65)-dduellipseOrigInU(0.65)*0.65,dduellipseOrigInU(0.65))
+# abline(v=c(0.8,0.65), col="grey")
+# ## okay
+
+ddvellipseOrig <- function(vo, a, b) {
+ ddvyOrig(vNorm(vo, a, b), a, b)*ddvvNorm(vo, a, b)
+}
+
+# ## double check
+# curve(ellipseOrig,0,0.6,ylim=c(0,1))
+# curve(ddvellipseOrig,add=T,col="green")
+# abline(ellipseOrig(0.4)-ddvellipseOrig(0.4)*0.4,ddvellipseOrig(0.4))
+# abline(ellipseOrig(0.1)-ddvellipseOrig(0.1)*0.1,ddvellipseOrig(0.1))
+# abline(v=c(0.1,0.4), col="grey")
+# abline(h=0)
+# ## okay
+
+ddvuEllipseOrig <- function(vo, a, b, dduA, dduB) {
+ cvNorm <- vNorm(vo, a, b)
+ (ddvuvNorm(vo, a, b, dduA, dduB)*ddvyOrig(cvNorm, a, b)
+ + ddvvNorm(vo, a, b)*(dduvNorm(vo, a, b, dduA, dduB)*ddvvyOrig(cvNorm, a, b)
+ + ddvuyOrig(cvNorm, a, b, dduA, dduB)))
+}
+
+# ## double check
+# curve(dduellipseOrig,0,0.5,ylim=c(-1,1.5))
+# curve(ddvuEllipseOrig,add=T)
+# abline(dduellipseOrig(0.4)-ddvuEllipseOrig(0.4)*0.4,ddvuEllipseOrig(0.4))
+# abline(dduellipseOrig(0.1)-ddvuEllipseOrig(0.1)*0.1,ddvuEllipseOrig(0.1))
+# abline(v=c(0.1,0.4), col="grey")
+# ## okay
+
+ddvvEllipseOrig <- function(vo, a, b) {
+ cvNorm <- vNorm(vo, a, b)
+ (ddvvNorm(vo, a, b)^2*ddvvyOrig(cvNorm, a, b)
+ +ddvvvNorm(vo, a, b)*ddvyOrig(cvNorm, a, b))
+}
+
+# ## double check
+# curve(ddvellipseOrig,0,0.5,ylim=c(-1,1.5))
+# curve(ddvvEllipseOrig,add=T)
+# abline(ddvellipseOrig(0.4)-ddvvEllipseOrig(0.4)*0.4,ddvvEllipseOrig(0.4))
+# abline(ddvellipseOrig(0.1)-ddvvEllipseOrig(0.1)*0.1,ddvvEllipseOrig(0.1))
+# abline(v=c(0.1,0.4), col="grey")
+# ## okay
+
+# cdf
+cdf.leaf <- function(uvabw) {
+ uvabw[5]+ellipseOrig(uvabw[2]-uvabw[5], uvabw[3], uvabw[4])
+}
+
+pLeafCopula <- function(u, copula) {
+ aVec <- solveA(u[,1], copula at parameters)
+ bVec <- solveB(u[,1], copula at parameters)
+ wBor <- weakBorderPoly(u[,1], copula at parameters)
+ sBor <- strongBorderPoly(u[,1])
+
+ ret <- apply(u,1,min)
+ bool <- (u[,2] > wBor & u[,2] < sBor)
+ if(any(bool))
+ ret[bool] <- apply(cbind(u,aVec,bVec, wBor)[bool,,drop=F], 1, cdf.leaf)
+ return(ret)
+}
+
+setMethod("pCopula",signature("matrix","leafCopula"), pLeafCopula)
+setMethod("pCopula",signature("numeric","leafCopula"),
+ function(u, copula) pLeafCopula(matrix(u,ncol=copula at dimension), copula))
+
+# partial derivative d/dv
+ddv.cdf.leaf <- function(uvabw) {
+ ddvellipseOrig(uvabw[2]-uvabw[5], uvabw[3], uvabw[4])
+}
+
+ddvLeafCopula <- function(u, copula) {
+ aVec <- solveA(u[,1], copula at parameters)
+ bVec <- solveB(u[,1], copula at parameters)
+ wBor <- weakBorderPoly(u[,1], copula at parameters)
+ sBor <- strongBorderPoly(u[,1])
+
+ ret <- numeric(nrow(u))
+ ret[u[,2] <= wBor] <- 1
+ bool <- (u[,2] > wBor & u[,2] < sBor)
+ if(any(bool))
+ ret[bool] <- apply(cbind(u,aVec,bVec, wBor)[bool,,drop=F], 1, ddv.cdf.leaf)
+ return(ret)
+}
+
+setMethod("ddvCopula",signature("matrix","leafCopula"), ddvLeafCopula)
+setMethod("ddvCopula",signature("numeric","leafCopula"),
+ function(u, copula) ddvLeafCopula(matrix(u,ncol=copula at dimension), copula))
+
+# partial derivative d/du
+ddu.cdf.leaf <- function(uvabw, par) {
+ (ddxweakBorderPoly(uvabw[1], par)*(1-ddvellipseOrig(uvabw[2]-uvabw[5], uvabw[3], uvabw[4]))
+ + dduellipseOrig(uvabw[2]-uvabw[5], uvabw[3], uvabw[4], ddusolveA(uvabw[1], par), ddusolveB(uvabw[1], par)))
+}
+
+dduLeafCopula <- function(u, copula) {
+ aVec <- solveA(u[,1], copula at parameters)
+ bVec <- solveB(u[,1], copula at parameters)
+ wBor <- weakBorderPoly(u[,1], copula at parameters)
+ sBor <- strongBorderPoly(u[,1])
+
+ ret <- numeric(nrow(u))
+ ret[u[,2] >= sBor] <- 1
+ bool <- (u[,2] > wBor & u[,2] < sBor)
+ if(any(bool))
+ ret[bool] <- apply(cbind(u,aVec,bVec, wBor)[bool,,drop=F], 1, ddu.cdf.leaf,
+ par=copula at parameters)
+ return(ret)
+}
+
+setMethod("dduCopula",signature("matrix","leafCopula"), dduLeafCopula)
+setMethod("dduCopula",signature("numeric","leafCopula"),
+ function(u, copula) dduLeafCopula(matrix(u, ncol=copula at dimension), copula))
+
+# density
+pdf.leaf <- function(uvabw, par) {
+ (ddvuEllipseOrig(uvabw[2]-uvabw[5], uvabw[3], uvabw[4],
+ ddusolveA(uvabw[1], par), ddusolveB(uvabw[1], par))
+ - ddxweakBorderPoly(uvabw[1], par)*ddvvEllipseOrig(uvabw[2]-uvabw[5],
+ uvabw[3], uvabw[4]))
+}
+
+dLeafCopula <- function(u, copula) {
+ aVec <- solveA(u[,1], copula at parameters)
+ bVec <- solveB(u[,1], copula at parameters)
+ wBor <- weakBorderPoly(u[,1], copula at parameters)
+ sBor <- strongBorderPoly(u[,1])
+
+ ret <- numeric(nrow(u))
+ bool <- (u[,2] > wBor & u[,2] < sBor)
+ if(any(bool))
+ ret[bool] <- apply(cbind(u,aVec,bVec,wBor)[bool,,drop=F],1,pdf.leaf,
+ par=copula at parameters)
+ return(ret)
+}
+
+setMethod("dCopula",signature("matrix","leafCopula"), dLeafCopula)
+setMethod("dCopula",signature("numeric","leafCopula"),
+ function(u, copula) dLeafCopula(matrix(u,ncol=copula at dimension), copula))
+
+## random pair generator
+
+# inverse of partial derivative d/dv
+invddvLeafCop <- function(vy, par) {
+
+ optFun <- function(u) {
+ y <- ddv.cdf.leaf(c(u,vy[1], solveA(u, par), solveB(u, par), weakBorderPoly(u, par)))
+ ret <- abs(y-vy[2])
+ if(is.nan(ret) | is.infinite(ret)) {
+ cat("The warning occured for the points: \n")
+ cat("u:", u,"\n")
+ cat("y:", y,"\n")
+ cat("vy", vy,"\n")
+ }
+ return(ret)
+ }
+
+ invWBor <- invWeakBor(vy[1], par)
+ invSBor <- invStrongBor(vy[1])
+
+ if(invSBor == invWBor)
+ return(invSBor)
+ optimise(optFun,c(min(invSBor,invWBor),max(invSBor,invWBor)))$minimum
+}
+
+r.leaf <- function(n, copula) {
+ v <- runif(n, min = 0, max = 1)
+ y <- runif(n, min = 0, max = 1)
+
+ res <- cbind(apply(cbind(v,y), 1, invddvLeafCop, par=copula at parameters), v)
+ colnames(res) <- c("u","v")
+
+ return(res)
+}
+
+setMethod("rCopula",signature("numeric","leafCopula"),r.leaf)
+
+# ## example section
+# leafCop <- leafCopula()
+# plot(rCopula(500,leafCop))
+#
+# persp(leafCop,dCopula,zlim=c(0,10))
+# persp(leafCop,pCopula)
\ No newline at end of file
Modified: spcopula_0.1-1.tar.gz
===================================================================
(Binary files differ)
Modified: spcopula_0.1-1.zip
===================================================================
(Binary files differ)
More information about the spcopula-commits
mailing list