Skip to content

Commit

Permalink
Allowing p=Inf in BiCopParamDistLp
Browse files Browse the repository at this point in the history
  • Loading branch information
@AlexisDerumigny
AlexisDerumigny committed Apr 12, 2022
1 parent ce387e7 commit 9bcd88b
Show file tree
Hide file tree
Showing 2 changed files with 151 additions and 33 deletions.
18 changes: 16 additions & 2 deletions NEWS.md
View file
Original file line number Diff line number Diff line change
@@ -1,13 +1,27 @@

# MMDCopula 0.2.1

* New functions `BiCopPar2Tau.MO` and `BiCopTau2Par.MO` to convert between the parameter and the Kendall's tau of a Marshall Olkin copula.

## NEW FEATURES

* New functions `BiCopPar2Tau.MO` and `BiCopTau2Par.MO` to convert between the parameter and the Kendall's tau of a Marshall-Olkin copula.

* New argument `truncVal` of the function `BiCopParamDistLp` to allow for computation of the norm on a smaller subset of the unit square [0,1]^2.
`BiCopParamDistLp` now also works for `p=Inf` (i.e., the supremum norm).

* New checks for finiteness of the arguments to `BiCopEstMMD`.


## OTHER IMPROVEMENTS

* Improvement of default values for gamma in `BiCopEstMMD`.
* Updated the references

* Updated the references.

* Fixed the documentation of the parameter `kernel` for the functions `BiCopEstMMD` and `BiCopGradMMD`.



# MMDCopula 0.2.0

* Change email address of maintainer Alexis Derumigny corresponding to new affiliation.
Expand Down
166 changes: 135 additions & 31 deletions R/BiCopParamDistLp.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,13 +22,16 @@
#' @param p determines the \eqn{L_p} distance that is used.
#'
#' @param type type of the functions considered.
#' Can be \code{cdf} for the distance between the two cumulative distribution functions
#' or \code{pdf} for the distance between the two probability density functions.
#' Can be \code{"cdf"} for the distance between the two cumulative distribution functions
#' or \code{"pdf"} for the distance between the two probability density functions.
#'
#' @param maxEval maximum number of evaluation of the function be integrated.
#' If 0, then no maximum limit is given.
#' @param maxEval maximum number of evaluation of the function to be integrated.
#' If 0, then no maximum limit is given. (Only used if \code{p < Inf}).
#'
#' @return a list of four items
#' @param truncVal the distance is computed using the supremum or the integral
#' of the function on \eqn{[truncVal, 1 - truncVal]^2}.
#'
#' @return If \code{p < Inf}, it returns a list of four items
#' \itemize{
#' \item \code{distance} the value of the distance
#' \item \code{integral} the value of the integral, which is
Expand All @@ -38,11 +41,18 @@
#' called by \code{cubature::\link[cubature]{hcubature}()}.
#' This should be 0 if there is no error.
#' }
#' If \code{p = Inf}, it returns a list of two items
#' \itemize{
#' \item \code{distance} the maximum difference between the two copulas
#' (respectively, between the two copula densities).
#' \item \code{u_max} the point at which this difference is attained.
#' }
#'
#' @examples
#' # Distance between the densities of a Gaussian copula with correlation 0.5
#' # and a Gaussian copula with correlation 0.2
#' BiCopParamDistLp(family = 1, par = 0.5, par_p = 0.2, p = 2, type = "cdf", maxEval = 10)
#' BiCopParamDistLp(family = 1, par = 0.5, par_p = 0.2, p = Inf, type = "cdf")
#'
#' # Distance between the cdf of a Student copula
#' # with correlation 0.5 and 4 degrees of freedom
Expand All @@ -58,38 +68,132 @@
#' @export
#'
BiCopParamDistLp <- function(family, par, par_p, par2 = par, par2_p = par_p, family_p = family,
p, type, maxEval = 0)
p, type, maxEval = 0, truncVal = 0)
{
switch(
type,
"cdf" = {
toIntegrate <- function(u) {
return ( matrix(( abs(
VineCopula::BiCopCDF(u1 = u[1,], u2 = u[2,],
family = family, par = par, par2 = par2) -
if (p < Inf) {
switch(
type,
"cdf" = {
toIntegrate <- function(u) {
return ( matrix(( abs(
VineCopula::BiCopCDF(u1 = u[1,], u2 = u[2,],
family = family_p, par = par_p, par2 = par2_p)
))^p, nrow = 1, ncol = length(u[1,])) )
}
},
family = family, par = par, par2 = par2) -
VineCopula::BiCopCDF(u1 = u[1,], u2 = u[2,],
family = family_p, par = par_p, par2 = par2_p)
))^p, nrow = 1, ncol = length(u[1,])) )
}
},

"pdf" = {
toIntegrate <- function(u) {
return ( matrix(( abs(
VineCopula::BiCopPDF(u1 = u[1,], u2 = u[2,],
family = family, par = par, par2 = par2) -
"pdf" = {
toIntegrate <- function(u) {
return ( matrix(( abs(
VineCopula::BiCopPDF(u1 = u[1,], u2 = u[2,],
family = family_p, par = par_p, par2 = par2_p)
))^p, nrow = 1, ncol = length(u[1,])) )
}
})
family = family, par = par, par2 = par2) -
VineCopula::BiCopPDF(u1 = u[1,], u2 = u[2,],
family = family_p, par = par_p, par2 = par2_p)
))^p, nrow = 1, ncol = length(u[1,])) )
}
})

result = cubature::hcubature(f = toIntegrate,
lower = c(truncVal, truncVal),
upper = c(1 - truncVal, 1 - truncVal),
vectorInterface = TRUE, maxEval = maxEval)
result[["distance"]] = result$integral^(1/p)

return (result[c("distance", "integral", "error", "returnCode")])

} else if (p == Inf){
# if (maxEval == 0){maxEval <- 10000}
# nGrid = floor(sqrt(maxEval))
# univGrid = seq(0 , nGrid/(nGrid+1), length = nGrid) + 1/(2 * nGrid)
# u = expand.grid(univGrid , univGrid)
# switch(
# type,
# "cdf" = {
# diff_ =
# VineCopula::BiCopCDF(u1 = u[,1], u2 = u[,2],
# family = family, par = par, par2 = par2) -
# VineCopula::BiCopCDF(u1 = u[,1], u2 = u[,2],
# family = family_p, par = par_p, par2 = par2_p)
# },
#
# "pdf" = {
# diff_ =
# VineCopula::BiCopPDF(u1 = u[,1], u2 = u[,2],
# family = family, par = par, par2 = par2) -
# VineCopula::BiCopPDF(u1 = u[,1], u2 = u[,2],
# family = family_p, par = par_p, par2 = par2_p)
# })
#
# result = max(abs(diff_))
# which_ = which.max(abs(diff_))
# u_max = as.numeric(u[which_, ])


switch(
type,
"cdf" = {
fn <- function(u) {
return ( - (VineCopula::BiCopCDF(u1 = u[1], u2 = u[2],
family = family, par = par, par2 = par2) -
VineCopula::BiCopCDF(u1 = u[1], u2 = u[2],
family = family_p, par = par_p, par2 = par2_p)
)^2 )
}
gr <- function(u) {
return (
- 2 *
(unlist(VineCopula::BiCopHfunc(u1 = u[1], u2 = u[2],
family = family, par = par, par2 = par2)) -
unlist(VineCopula::BiCopHfunc(u1 = u[1], u2 = u[2],
family = family_p, par = par_p, par2 = par2_p))
) * (VineCopula::BiCopCDF(u1 = u[1], u2 = u[2],
family = family, par = par, par2 = par2) -
VineCopula::BiCopCDF(u1 = u[1], u2 = u[2],
family = family_p, par = par_p, par2 = par2_p)
))
}
},

"pdf" = {
fn <- function(u) {
return ( - (
VineCopula::BiCopPDF(u1 = u[1], u2 = u[2],
family = family, par = par, par2 = par2) -
VineCopula::BiCopPDF(u1 = u[1], u2 = u[2],
family = family_p, par = par_p, par2 = par2_p) )^2)
}

gr <- function(u) {
return (
- 2 *
c(VineCopula::BiCopDeriv(u1 = u[1], u2 = u[2],
family = family, par = par, par2 = par2,
deriv = "u1") -
VineCopula::BiCopDeriv(u1 = u[1], u2 = u[2],
family = family_p, par = par_p, par2 = par2_p,
deriv = "u1")
,
VineCopula::BiCopDeriv(u1 = u[1], u2 = u[2],
family = family, par = par, par2 = par2,
deriv = "u2") -
VineCopula::BiCopDeriv(u1 = u[1], u2 = u[2],
family = family_p, par = par_p, par2 = par2_p,
deriv = "u2")
) * (VineCopula::BiCopPDF(u1 = u[1], u2 = u[2],
family = family, par = par, par2 = par2) -
VineCopula::BiCopPDF(u1 = u[1], u2 = u[2],
family = family_p, par = par_p, par2 = par2_p)
))
}
})

result = cubature::hcubature(f = toIntegrate,
lower = c(0,0), upper = c(1,1),
vectorInterface = TRUE, maxEval = maxEval)
result[["distance"]] = result$integral^(1/p)
result = optim(par = c(0.2, 0.2), fn = fn, gr = gr,
lower = truncVal, upper = 1 - truncVal, method = "L-BFGS-B")

return (result[c("distance", "integral", "error", "returnCode")])
return (list(distance = -result$value, u_max = result$par ))
}
}


0 comments on commit 9bcd88b

Please sign in to comment.