# R/berncop.R In kdecopula: Kernel Smoothing for Bivariate Copula Densities

```#' Calculates a Bernstein polynomial
#'
#' @param x vector of evaluation points
#' @param k index of the Bernstein polynomial
#' @param m order of the Bernstein polynomial
#' @noRd
bern_poly <- function(x, k, m) {
(choose(m, k) * x^k * (1 - x)^(m - k)) * (m + 1)
}

#' Calculates basis coefficients fo the Bernstein estimator
#'
#' The method of Weiss and Scheffer (2016) is used to adjust the coefficients
#' to uniform margins.
#'
#' @param u data
#' @param m order of the Bernstein polynomial
#' @noRd
bern_coefs <- function(u, m) {
# initialize coefficients with empirical frequencies
ucut <- lapply(1:2, function(i) cut(u[, i], 0:(m + 1) / (m + 1)))
cf0 <- table(ucut[[1]], ucut[[2]]) / nrow(u)
as.matrix(cf0)
# # set up quadratic programming problem
# m <- m + 1
# D <- diag(m^2)
# d <- as.vector(cf0)
# A1 <- A2 <- matrix(0, m, m^2)
# for (i in 1:m) {
#     A1[i, (i - 1) * m + 1:m] <- 1
#     A2[i, m * (1:m) - m + i] <- 1
# }
# A3 <- diag(m^2)
# A <- t(rbind(A1, A2, A3))
# b <- c(rep(1 / m, 2 * m), rep(0, m^2))
#
# # solve
# sol <- tryCatch(solve.QP(D, d, A, b, meq = 2 * m)\$solution,
#                 error = function(e) d)
# sol[sol < 0] <- 0

# return as matrix
# matrix(sol, m, m)
}

#' Fit a Bernstein copula to the data
#'
#' @param u data
#' @param m order of the Bernstein copula
#'
#' @return berncop object
#' @noRd
berncop <- function(u, m = 10) {
out <- list(coefs = bern_coefs(u, m),  m = m)
class(out) <- "berncop"
out
}

#' Evaluate the density of a berncop object
#'
#' @param unew data
#' @param object berncop object
#'
#' @return Bernstein copula density evaluated at uev.
#' @noRd
dberncop <- function(uev, object) {
if (NCOL(uev) == 1)
uev <- matrix(uev, ncol = 2)
list2env(object)
summand <- function(i, j) {
P1 <- bern_poly(uev[, 1], i, object\$m)
P2 <- bern_poly(uev[, 2], j, object\$m)
object\$coefs[i + 1, j + 1] * P1 * P2
}

est <- 0
for (i in 0:object\$m) {
phi_i <- sapply(0:object\$m, function(j) summand(i, j))
if (NCOL(phi_i) == 1)
phi_i <- t(phi_i)
est <- est + rowSums(phi_i)
}
est
}
```

## Try the kdecopula package in your browser

Any scripts or data that you put into this service are public.

kdecopula documentation built on May 2, 2019, 1:06 a.m.