Nothing
## 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)
}
# param <- c(1.446923, -1.722742)
# 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)
}
ddxsolveQ <- function(u) {
sqrt(0.5)*(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))+solveQ(u)
}
ddxsolveXb <- function(u, par) {
wBor <- weakBorderPoly(u, par)
sqrt(2)*(1-ddxweakBorderPoly(u, par))+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^2 - q^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*sqrt(a^2-xb^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@parameters)
bVec <- solveB(u[,1], copula@parameters)
wBor <- weakBorderPoly(u[,1], copula@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@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@parameters)
bVec <- solveB(u[,1], copula@parameters)
wBor <- weakBorderPoly(u[,1], copula@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@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@parameters)
bVec <- solveB(u[,1], copula@parameters)
wBor <- weakBorderPoly(u[,1], copula@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@parameters)
return(ret)
}
setMethod("dduCopula",signature("matrix","leafCopula"), dduLeafCopula)
setMethod("dduCopula",signature("numeric","leafCopula"),
function(u, copula) dduLeafCopula(matrix(u, ncol=copula@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@parameters)
bVec <- solveB(u[,1], copula@parameters)
wBor <- weakBorderPoly(u[,1], copula@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@parameters)
return(ret)
}
setMethod("dCopula",signature("matrix","leafCopula"), dLeafCopula)
setMethod("dCopula",signature("numeric","leafCopula"),
function(u, copula) dLeafCopula(matrix(u,ncol=copula@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@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)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.