Nothing
Bcopula <-
function(mat.xy, m, both.cont = FALSE, tolimit = 0.00001){
#
# Input: mat.xy = 2-column matrix with bivariate observations of a random vector (X,Y)
# ** repeated values are allowed, NA values not allowed
# m = Approximation order (an integer 2,3,...,n = sample size)
# ** it is recommended tu use m = min{50, n^(1/2)}
# both.cont = if TRUE then (X,Y) are considered (both) as continuos random variables,
# and jittering will be applied to repeated values (if any)
# tolimit = tolerance limit in numerical approximation of the inverse of the
# partial derivatives of the estimated Bernstein copula
#
# Output: copula = bivariate function Bernstein Copula (BC) of order m
# du = bivariate function dBCopula(u,v)/du
# dv = bivariate function dBCopula(u,v)/dv
# du.inv = inverse of BC.du with respect to v given u and alpha (numerical aprox with tol)
# dv.inv = inverse of BC.du with respect to u given v and alpha (numerical aprox with tol)
# density = bivariate function Bernstein copula density of order m
# bilinearCopula = bivariate function of bilinear approximation of copula
# bilinearSubcopula = (m+1) x (m+1) matrix with empirical subcopula values
# sample.size of bivariate observations
# order of approximation (m)
# both.cont = logical value, TRUE if both variables considered as continuous
# tolimit = tolerance limit in numerical approximation of du.inv and dv.inv
# subcopemObject = list with the output of subcopem (both.cont = FALSE)
# or subcopemc (both.cont = TRUE)
#
# Dependencies: subcopem and subcopemc functions in subcopem2D package
#
# Checking for appropriate order m
n <- nrow(mat.xy) # sample size
if (!(m %in% (2:n))){
error.msg <- paste("Order m must be an integer value in 2:n, n =", n)
stop(error.msg)
} else{
# Calculate bilinear interpolation of order m:
if (both.cont){
SC <- subcopemc(mat.xy, m)
} else{
SC <- subcopem(mat.xy)
}
um <- (0:m)/m
vm <- um
Cm <- matrix(0, nrow = (m + 1), ncol = (m + 1))
Cm[ , m + 1] <- um # C(i/m, 1) = i/m
Cm[m + 1, ] <- vm # C(1, i/m) = i/m
L <- function(u){
u.inf <- max(SC$part1[SC$part1 <= u])
u.sup <- min(SC$part1[SC$part1 >= u])
valor <- ifelse(u.inf < u.sup, (u - u.inf)/(u.sup - u.inf), 1)
return(valor)
}
M <- function(v){
v.inf <- max(SC$part2[SC$part2 <= v])
v.sup <- min(SC$part2[SC$part2 >= v])
valor <- ifelse(v.inf < v.sup, (v - v.inf)/(v.sup - v.inf), 1)
return(valor)
}
S <- SC$matrix # empirical subcopula
n1 <- nrow(S)
n2 <- ncol(S)
iu.inf <- function(u) sum(SC$part1 <= u)
iu.sup <- function(u) n1 - sum(SC$part1 >= u) + 1
iv.inf <- function(v) sum(SC$part2 <= v)
iv.sup <- function(v) n2 - sum(SC$part2 >= v) + 1
C.bilineal <- function(u, v) (1-L(u))*(1-M(v))*S[iu.inf(u),iv.inf(v)] +
(1-L(u))*M(v)*S[iu.inf(u),iv.sup(v)] +
L(u)*(1-M(v))*S[iu.sup(u),iv.inf(v)] +
L(u)*M(v)*S[iu.sup(u),iv.sup(v)]
for (i in 2:m){
for (j in 2:m){
Cm[i, j] <- C.bilineal(um[i], vm[j])
}
}
# Bernstein copula functions
BCopula <- function(u, v) sum(Cm*(dbinom(0:(dim(Cm)[1] - 1), dim(Cm)[1] - 1, u)%*%t(dbinom(0:(dim(Cm)[1] - 1), dim(Cm)[1] - 1, v))))
BC.du <- function(u, v) (dim(Cm)[1] - 1)*sum(Cm*((dbinom(-1:(dim(Cm)[1] - 2), dim(Cm)[1] - 2, u) - dbinom(0:(dim(Cm)[1] - 1),
dim(Cm)[1] - 2, u)*c(-1, rep(1, dim(Cm)[1] - 1)))%*%t(dbinom(0:(dim(Cm)[1] - 1), dim(Cm)[1] - 1, v))))
BC.du.aux <- function(v, ua.vec) BC.du(ua.vec[1], v) - ua.vec[2]
BC.du.inv <- function(u, a) uniroot(BC.du.aux, interval = c(0, 1), ua.vec = c(u, a), tol = tolimit)$root
tCm <- t(Cm) # transpose u <-> v
BC.dv0 <- function(u, v) (dim(tCm)[1] - 1)*sum(tCm*((dbinom(-1:(dim(tCm)[1] - 2), dim(tCm)[1] - 2, u) - dbinom(0:(dim(tCm)[1] - 1),
dim(tCm)[1] - 2, u)*c(-1, rep(1, dim(tCm)[1] - 1)))%*%t(dbinom(0:(dim(tCm)[1] - 1), dim(tCm)[1] - 1, v))))
BC.dv0.aux <- function(v, ua.vec) BC.dv0(ua.vec[1], v) - ua.vec[2]
BC.dv0.inv <- function(u, a) uniroot(BC.dv0.aux, interval = c(0, 1), ua.vec = c(u, a), tol = tolimit)$root
BC.dv <- function(u, v) BC.dv0(v, u)
BC.dv.inv <- function(v, a) BC.dv0.inv(v, a)
Bdensity <- function(u, v) ((dim(Cm)[1] - 1)^2)*sum(Cm*((dbinom(-1:(dim(Cm)[1] - 2), dim(Cm)[1] - 2, u) - dbinom(0:(dim(Cm)[1] - 1),
dim(Cm)[1] - 2, u)*c(-1, rep(1, dim(Cm)[1] - 1)))%*%t(dbinom(-1:(dim(Cm)[1] - 2), dim(Cm)[1] - 2, v) -
dbinom(0:(dim(Cm)[1] - 1), dim(Cm)[1] - 2, v) * c(-1, rep(1, dim(Cm)[1] - 1)))))
# Output:
lista <- list(copula = BCopula, du = BC.du, du.inv = BC.du.inv, dv = BC.dv, dv.inv = BC.dv.inv,
density = Bdensity, bilinearCopula = C.bilineal, bilinearSubcopula = Cm, sample.size = n,
order = m, both.cont = both.cont, tolerance = tolimit, subcopemObject = SC)
return(lista)
} # cierre de else inicial
}
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.