# Density estimation functions
# ----------------------------
#' Bivariate density estimation
#'
#' \code{dlg_bivariate} returns the locally Gaussian density estimate of a
#' bivariate distribution on a given grid.
#'
#' This function serves as the backbone in the body of methods concerning local
#' Gaussian correlation. It takes a bivariate data set, \code{x}, and a
#' bivariate set of grid points \code{eval_points}, and returns the bivariate,
#' locally Gaussian density estimate in these points. We also need a vector of
#' bandwidths, \code{bw}, with two elements, and an estimation method
#' \code{est_method}
#'
#' @param x The data matrix (or data frame). Must have exactly 2 columns.
#' @param eval_points The grid where the density should be estimated. Must have
#' exactly 2 columns.
#' @param grid_size If \code{eval_points} is not supplied, then the function
#' will create a suitable grid diagonally through the data, with this many
#' grid points.
#' @param bw The two bandwidths, a numeric vector of length 2.
#' @param est_method The estimation method, must either be "1par" for estimation
#' with just the local correlation, or "5par" for a full locally Gaussian fit
#' with all 5 parameters.
#' @param tol The numerical tolerance to be used in the optimization. Only
#' applicable in the 1-parameter optimization.
#' @param run_checks Logical. Should sanity checks be run on the arguments?
#' Useful to disable this when doing cross-validation for example.
#' @param marginal_estimates Provide the marginal estimates here if estimation
#' method is "\code{5par_marginals_fixed}", and the marginal estimates have
#' already been found. Useful for cross-validation. List with two elements as
#' returned by \code{dlg_marginal_wrapper}.
#' @param bw_marginal Vector of bandwidths used to estimate the marginal
#' distributions.
#'
#' @return A list including the data set \code{$x}, the grid
#' \code{$eval_points}, the bandwidths \code{$bw}, as well as a matrix of the
#' estimated parameter estimates \code{$par_est} and the estimated bivariate
#' density \code{$f_est}.
#'
#' @examples
#' x <- cbind(rnorm(100), rnorm(100))
#' bw <- c(1, 1)
#' eval_points <- cbind(seq(-4, 4, 1), seq(-4, 4, 1))
#'
#' estimate <- dlg_bivariate(x, eval_points = eval_points, bw = bw)
#'
#' @export
dlg_bivariate <- function(x,
eval_points = NA,
grid_size = 15,
bw = c(1, 1),
est_method = "1par",
tol = .Machine$double.eps^0.25/10^4,
run_checks = TRUE,
marginal_estimates = NA,
bw_marginal = NA) {
# Check that everything is the way it should be -----------------------------
if(run_checks) {
x <- check_data(x, dim_check = 2, type = "data")
check_bw_bivariate(bw = bw)
check_est_method(est_method)
if(!identical(eval_points, NA))
eval_points <- check_data(eval_points, dim_check = 2, type = "grid")
if(!identical(bw_marginal, NA))
check_bw_bivariate(bw = bw_marginal)
}
# If eval_points not supplied, create a suitable grid
if(identical(eval_points, NA)) {
eval_points <- apply(x,
2,
function(x) seq(quantile(x,.001),quantile(x,.999),
length.out = grid_size))
}
# Preliminary housekeeping - mostly so that we can re-use code from the previous
# package that is working alright. We split the data matrix and the bandwidth vector
# into their individual components:
x1 <- x[,1]
x2 <- x[,2]
h1 <- bw[1]
h2 <- bw[2]
# OPTION 1: THE ONE-PARAMETER OPTIMIZATION
if((est_method == "1par")) {
# We declare a function that maximizes the local likelihood in one grid point.
maximize_likelihood = function(grid_point) {
x1_0 <- grid_point[1]
x2_0 <- grid_point[2]
# We need weights and some empirical moments in this grid point
W <- dnorm(x1, mean = x1_0, sd = h1)*dnorm(x2, mean = x2_0, sd = h2)
m1 <- mean(W)
m2 <- mean(W*x1^2)
m3 <- mean(W*x2^2)
m4 <- mean(W*x1*x2)
# The likelihood function
lik <- function(rho) {
- log(2*pi*sqrt(1 - rho^2))*m1 - m2/(2*(1 - rho^2)) - m3/(2*(1 - rho^2)) +
rho*m4/(1 - rho^2) - 1/2*exp(-1/2*(x2_0^2*h1^2 + x1_0^2 + x1_0^2*h2^2 -
2*x1_0*rho*x2_0 + x2_0^2)/(-rho^2 + h2^2 + 1 + h1^2 + h1^2*h2^2))/
(pi*(-rho^2 + h2^2 + 1 + h1^2 + h1^2*h2^2)^(1/2))
}
# Return the maximum of the likelihood and the density estimate
opt <- try(optimise(lik,
lower = -1,
upper = 1,
maximum = TRUE,
tol = tol),
silent = TRUE)
# Store the result if the optimization went alright. Return NA if not.
if(class(opt) != "try-error") {
return(c(opt$maximum,
mvtnorm::dmvnorm(c(x1_0, x2_0), mean = c(0,0),
sigma = matrix(c(1, opt$maximum, opt$maximum, 1), 2))))
} else {
return(c(NA, NA))
}
}
## Send the grid points to 'maximize_likelihood'
est <- cbind(do.call(rbind,
lapply(X = split(eval_points, row(eval_points)),
FUN = maximize_likelihood)))
par_est <- matrix(est[,1])
f_est = as.vector(est[,2])
colnames(par_est) <- c('rho')
}
# OPTION 2: THE FIVE PARAMETER OPTIMIZATION
if(est_method == "5par") {
# Use the localgauss-package to find the local parameter estimates
localgauss_result <- localgauss::localgauss(x = x1, y = x2,
b1 = h1, b2 = h2,
xy.mat = eval_points)
# Collect the results
par_est <- localgauss_result$par.est
f_est <- dmvnorm_wrapper(eval_points = eval_points,
mu_1 = par_est[, "mu_1"],
mu_2 = par_est[, "mu_2"],
sig_1 = par_est[, "sig_1"],
sig_2 = par_est[, "sig_2"],
rho = par_est[, "rho"],
run_checks = run_checks)
}
# OPTION 3: FIVE PARAMETERS WITH MARGINALS FIXED
if(est_method == "5par_marginals_fixed") {
# Have the marginals already been estimated? If not, do it
if(identical(marginal_estimates, NA)) {
marginal_estimates <- dlg_marginal_wrapper(data_matrix = x,
eval_matrix = eval_points,
bw_vector = bw_marginal)
}
# We declare a function that maximizes the local likelihood in one grid point,
# with fixed marginal parameter estimates.
maximize_likelihood_fixed_marginals =
function(grid_point_and_marginal_parameters) {
x1_0 <- grid_point_and_marginal_parameters[1]
x2_0 <- grid_point_and_marginal_parameters[2]
mu_1 <- grid_point_and_marginal_parameters[3]
sig_1 <- grid_point_and_marginal_parameters[4]
mu_2 <- grid_point_and_marginal_parameters[5]
sig_2 <- grid_point_and_marginal_parameters[6]
# We need weights and some empirical moments in this grid point
W <- dnorm(x1, mean = x1_0, sd = h1)*dnorm(x2, mean = x2_0, sd = h2)
m1 <- mean(W)
m2 <- mean(W*((x1-mu_1)/sig_1)^2)
m3 <- mean(W*((x2-mu_2)/sig_2)^2)
m4 <- mean(W*((x1-mu_1)/sig_1)*((x2-mu_2)/sig_2))
# The likelihood function
lik <- function(rho) {
- log(2*pi*sqrt(1 - rho^2)*sig_1*sig_2)*m1 - m2/(2*(1 - rho^2)) - m3/(2*(1 - rho^2)) +
rho*m4/(1 - rho^2) -
mvtnorm::dmvnorm(c(x1_0, x2_0),
mean = c(mu_1, mu_2),
sigma = matrix(c(sig_1^2 + h1^2,
sig_1*sig_2*rho,
sig_1*sig_2*rho,
sig_2^2 + h2^2),
byrow = TRUE,
ncol = 2))
}
# Return the maximum of the likelihood and the density estimate
opt <- try(optimise(lik,
lower = -1,
upper = 1,
maximum = TRUE,
tol = tol),
silent = TRUE)
# Store the result if the optimization went alright. Return NA if not.
if(class(opt) != "try-error") {
return(c(opt$maximum,
mvtnorm::dmvnorm(c(x1_0, x2_0), mean = c(mu_1,mu_2),
sigma = matrix(c(sig_1^2,
sig_1*sig_2*opt$maximum,
sig_1*sig_2*opt$maximum,
sig_2^2), 2))))
} else {
return(c(NA, NA))
}
}
## Send the grid points to 'maximize_likelihood_fixed_marginals'
est <- cbind(do.call(rbind,
lapply(X = split(cbind(eval_points,
marginal_estimates[[1]], marginal_estimates[[2]]),
row(eval_points)),
FUN = maximize_likelihood_fixed_marginals)))
par_est <- cbind(marginal_estimates[[1]][,1],
marginal_estimates[[2]][,1],
marginal_estimates[[1]][,2],
marginal_estimates[[2]][,2],
matrix(est[,1]))
f_est = as.vector(est[,2])
colnames(par_est) <- c("mu_1", "mu_2", "sig_1", "sig_2", "rho")
}
return(list(x = x,
eval_points = eval_points,
bw = bw,
par_est = par_est,
f_est = f_est))
}
#' Marginal density estimation
#'
#' Function that estimates a univariate density estimation by local Gaussian
#' approximations, as described in Hufthammer and Tjøstheim (2009).
#'
#' This function is mainly mean to be used as a tool in multivariate analysis as
#' away to obtain the estimate of a univariate (marginal) density function, but
#' it can of course be used in general to estimate univariate densities.
#'
#' @param x The data vector.
#' @param bw The bandwidth (a single number).
#' @param grid_size Number of grid points if grid is not provided.
#' @param eval_points The grid where we want to evaluate the density. Chosen
#' suitably if not provided, with length equal to grid_size.
#'
#' @return A list including the data set \code{$x}, the grid
#' \code{$eval_points}, the bandwidth \code{$bw}, as well as a matrix of the
#' estimated parameter estimates \code{$par_est} and the estimated bivariate
#' density \code{$f_est}.
#'
#' @examples
#' x <- rnorm(100)
#' estimate <- dlg_marginal(x, bw = 1, eval_points = -4:4)
#'
#' @references
#'
#' Hufthammer, Karl Ove, and Dag Tjøstheim. "Local Gaussian Likelihood and Local
#' Gaussian Correlation" PhD Thesis of Karl Ove Hufthammer, University of
#' Bergen, 2009.
#'
#' @export
dlg_marginal <- function(x,
bw = 1,
eval_points = seq(quantile(x, .01), quantile(x, .99), length.out = grid_size),
grid_size = 15) {
# Function that maximizes the local Gaussian likelihood in one grid point
maximize_likelihood_univariate <- function(eval_point) {
W <- dnorm(x, mean = eval_point, sd = bw)
m1 <- mean(W)
m2 <- mean(W*x)
m3 <- mean(W*(x^2))
# The negative likelihood
lik <- function(par) {
mu <- par[1]
sig <- par[2]
- m2*mu/sig^2 + m3/(2*sig^2) + m1*(mu^2/(2*sig^2) + .5*log(2*pi*sig^2)) +
dnorm((eval_point - mu)/sqrt(sig^2 + bw^2))/sqrt(sig^2 + bw^2)
}
# Do the optimization
opt <- try(optim(c(0, 1),
lik))
# Return the result, take absolute value of the sd.
if(class(opt) != "try-error") {
return(c(opt$par[1], abs(opt$par[2]),
dnorm(eval_point,
mean = opt$par[1],
sd = abs(opt$par[2]))))
} else {
return(c(NA, NA, NA))
}
}
# Send the grid points to 'maximize_likelihood_univariate'
est <- cbind(do.call(rbind,
lapply(X = as.list(eval_points),
FUN = maximize_likelihood_univariate)))
# Collect and return the results
if(nrow(est) == 1) {
par_est <- matrix(est[, 1:2], ncol = 2)
} else {
par_est <- est[, 1:2]
}
colnames(par_est) <- c("mu", "sig")
f_est <- as.vector(est[, 3])
return(list(x = x,
eval_points = eval_points,
bw = bw,
par_est = par_est,
f_est = f_est))
}
#' Marginal estimates for multivariate data
#'
#' Estimates the marginal locally Gaussian parameters for a
#' multivariate data set
#'
#' This function takes in a matrix of observations, a matrix of evaluation
#' points and a vector of bandwidths, and does a locally Gaussian fit on each of
#' the marginals using the \code{dlg_bivariate}-function. This function assumes
#' that the data and evaluation points are organized column-wise in matrices,
#' and that the bandwidth is found in the corresponding element in the bandwidth
#' matrix. The primary use for this function is multivariate density estimation
#' using the "5par_marginals_fixed"-method.
#'
#' @param data_matrix The matrix of data points. One column constitutes an
#' observation vector.
#' @param eval_matrix The matrix of evaluation points. One column constitutes a
#' vector of grid points.
#' @param bw_vector The vector of bandwidths, one element per component.
#'
#' @return A list with marginal parameter and density estimates as provided by
#' the \code{dlg_bivariate}-function. One element per column in the data.
#'
#' @examples
#' data_matrix <- cbind(rnorm(100), rnorm(100))
#' eval_matrix <- cbind(seq(-4, 4, 1), seq(-4, 4, 1))
#' bw <- c(1, 1)
#'
#' estimate <- dlg_marginal_wrapper(data_matrix, eval_matrix = eval_matrix, bw = bw)
#'
#' @export
dlg_marginal_wrapper <- function(data_matrix, eval_matrix, bw_vector){
# Return the parameters in a list
ret <- list()
# estimate the parameters for each marginal component
for(i in 1:ncol(data_matrix)) {
ret[[i]] <- dlg_marginal(x = data_matrix[,i],
bw = bw_vector[i],
eval_points = eval_matrix[,i])$par_est
}
return(ret)
}
#' The locally Gaussian density estimator (LGDE)
#'
#' Estimate a multivariate density function using locally Gaussian
#' approximations
#'
#' This function does multivariate density estimation using the locally Gaussian
#' density estimator (LGDE), that was introduced by Otneim & Tjøstheim (2017).
#' The function takes as arguments an \code{lg}-object as produced by the main
#' \code{lg_main}-function (where all the running parameters are specified), and
#' a grid of points where the density estimate should be estimated.
#'
#' @param lg_object An object of type \code{lg}, as produced by the
#' \code{lg_main}-function.
#' @param grid A matrix of grid points, where we want to evaluate the density
#' estimate.
#' @param level Specify a level if asymptotic standard deviations and confidence
#' intervals should be returned.
#' @param normalization_points How many grid points for approximating the integral
#' of the density estimate, to use for normalization?
#' @param bootstrap Calculate bootstrapped confidence intervals instead.
#' @param B Number of bootstrap replications if using bootstrapped confidence
#' intervals.
#'
#' @return A list containing the density estimate as well as all the running
#' parameters that has been used. The elements are:
#'
#' \itemize{
#' \item \code{f_est}: The estimated multivariate density.
#' \item \code{loc_mean}: The estimated local means if \code{est_method}
#' is "5par" or "5par_marginals_fixed", a matrix of zeros if
#' \code{est_method} is "1par".
#' \item \code{loc_sd}: The estimated local st. deviations if
#' \code{est_method} is "5par" or "5par_marginals_fixed", a matrix
#' of ones if \code{est_method} is "1par".
#' \item \code{loc_cor}: Matrix of estimated local correlations, one
#' column for each pair of variables, in the same order as specified
#' in the bandwidth object.
#' \item \code{x}: The data set.
#' \item \code{bw}: The bandwidth object.
#' \item \code{transformed_data}: The data transformed to approximate
#' marginal standard normality.
#' \item \code{normalizing_constants}: The normalizing constants used to
#' transform data and grid back and forth to the marginal standard
#' normality scale, as seen in eq. (8) of Otneim & Tjøstheim (2017).
#' \item \code{grid}: The grid where the estimation was performed, on the
#' original scale.
#' \item \code{transformed_grid}: The grid where the estimation was
#' performed, on the marginal standard normal scale.
#' \item \code{normalization_points} Number of grid points used
#' to approximate the integral of the density estimate, in order to
#' normalize?
#' \item \code{normalization_constant} If approximated, the integral of the
#' non-normalized density estimate. NA if not normalized.
#' \item \code{density_normalized} Logical, indicates whether the final
#' density estimate (contained in f_est) has been approximately
#' normalized to have unit integral.
#' \item \code{loc_cor_sd} Estimated asymptotic standard deviation for the
#' local correlations.
#' \item \code{loc_cor_lower} Lower confidence limit based on the asymptotic
#' standard deviation.
#' \item \code{loc_cor_upper} Upper confidence limit based on the asymptotic
#' standard deviation.
#' }
#'
#' @examples
#' x <- cbind(rnorm(100), rnorm(100), rnorm(100))
#' lg_object <- lg_main(x) # Put all the running parameters in here.
#' grid <- cbind(seq(-4, 4, 1), seq(-4, 4, 1), seq(-4, 4, 1))
#' density_estimate <- dlg(lg_object, grid = grid)
#'
#' @references
#'
#' Otneim, Håkon, and Dag Tjøstheim. "The locally gaussian density estimator for
#' multivariate data." Statistics and Computing 27, no. 6 (2017): 1595-1616.
#'
#' @export
dlg <- function(lg_object, grid, level = 0.95, normalization_points = NULL, bootstrap = F, B = 500) {
# Do some checks first
check_lg(lg_object)
if(!is.null(grid)) {
grid <- check_data(grid,
dim_check = ncol(lg_object$transformed_data),
type = "grid")
}
if(!is.null(normalization_points)) {
if(!is.numeric(normalization_points)) {
stop("normalization_points must either be NA (for no normalization) or an integer.")
}
}
# Sample size and number of variables
n <- nrow(lg_object$x)
d <- ncol(lg_object$x)
# If we are going to normalize the density estimate, we calculate it on a
# grid of size normalization_points, that we generate from a multivariate
# normal having the same mean and covariance matrix as the original data. We
# use this to calculate the integral. This is only temporary, and we will
# remove this in the end.
if(!is.null(normalization_points)) {
extra_grid <- mvtnorm::rmvnorm(n = normalization_points,
mean = colMeans(x),
sigma = cov(x))
grid <- rbind(grid, extra_grid)
}
# We do not have asymptotic standard deviations for the trivariate local
# correlations (yet)
if((lg_object$est_method == "trivariate") & (!is.null(level)) & !bootstrap) {
level <- NULL
message("Asymptotic standard deviations are not available for trivariate fits. Use bootstrap instead.")
}
# Extract the data that we will us in the optimization, transform the grid if needed
if(lg_object$transform_to_marginal_normality) {
x <- lg_object$transformed_data
transformed <- lg_object$trans_new(grid)
x0 <- transformed$trans
normalizing_constants <- transformed$normalizing_constants
} else {
x <- lg_object$x
x0 <- grid
normalizing_constants <- matrix(1, nrow = nrow(grid), ncol = d)
}
# Extract the pairs from the list of bandwidths
pairs <- lg_object$bw$joint[, c("x1", "x2")]
# Initialize the matrices where we eventually are going to store the parameter estimates
loc_mean <- matrix(0, ncol = d, nrow = nrow(grid))
loc_sd <- loc_mean + 1
loc_cor <- matrix(0, ncol = nrow(pairs), nrow = nrow(grid))
# If method "5par_marginals_fixed" we estimate the marginals now.
if(lg_object$est_method == "5par_marginals_fixed") {
marginal_estimates <- dlg_marginal_wrapper(data_matrix = x,
eval_matrix = x0,
bw_vector = lg_object$bw$marginal)
for(i in 1:d) {
loc_mean[,i] <- marginal_estimates[[i]][, "mu"]
loc_sd[,i] <- marginal_estimates[[i]][, "sig"]
}
}
# If method is "5par", then we have only one pair, and we can do all of the estimation now
if(lg_object$est_method == "5par") {
estimate <- dlg_bivariate(x = x,
eval_points = x0,
bw = c(lg_object$bw$joint[1, "bw1"],
lg_object$bw$joint[1, "bw2"]),
est_method = lg_object$est_method)
loc_mean <- cbind(estimate$par_est[, "mu_1"],
estimate$par_est[, "mu_2"])
loc_sd <- cbind(estimate$par_est[, "sig_1"],
estimate$par_est[, "sig_2"])
loc_cor[,1] <- estimate$par_est[, "rho"]
} else if (lg_object$est_method == "trivariate"){
estimate <- dlg_trivariate(x = x,
eval_points = x0,
bw = c(lg_object$bw$joint[1, "bw1"],
lg_object$bw$joint[1, "bw2"],
lg_object$bw$joint[1, "bw3"]),
est_method = lg_object$est_method)
loc_cor <- estimate$par_est
f_est_trivariate <- estimate$f_est
} else {
for(i in 1:nrow(pairs)) {
if(lg_object$est_method == "1par") {
pairwise_marginal_estimates <- NA
} else {
pairwise_marginal_estimates <- list(marginal_estimates[[pairs$x1[i]]],
marginal_estimates[[pairs$x2[i]]])
}
pairwise_estimate <-
dlg_bivariate(x = x[, c(pairs$x1[i], pairs$x2[i])],
eval_points = x0[, c(pairs$x1[i], pairs$x2[i])],
bw = c(lg_object$bw$joint[i, "bw1"],
lg_object$bw$joint[i, "bw2"]),
est_method = lg_object$est_method,
marginal_estimates = pairwise_marginal_estimates)
loc_cor[,i] <- pairwise_estimate$par_est[, "rho"]
}
}
# Evaluate the density estimate
if(lg_object$est_method != "trivariate") {
f_est <- mvnorm_eval(eval_points = x0,
loc_mean = loc_mean,
loc_sd = loc_sd,
loc_cor = loc_cor,
pairs = pairs)*
apply(normalizing_constants, 1, prod)
} else {
f_est <- f_est_trivariate*apply(normalizing_constants, 1, prod)
}
# If we normalize, calculate the constant and remove the extra grid points.
# If not, set the normalizing constant equal to one.
if(!is.null(normalization_points)) {
normalization_grid <- tail(grid, n = normalization_points)
normalization_est <- tail(f_est, n = normalization_points)
normalization_constant <- mean(normalization_est/
mvtnorm::dmvnorm(normalization_grid,
mean = colMeans(x),
sigma = cov(x)))
# Remove the extra grid points that were added to perform normalization of
# the estimated density.
grid <- grid[1:(nrow(grid) - normalization_points), ]
loc_mean <- loc_mean[1:nrow(grid),]
loc_sd <- loc_sd[1:nrow(grid),]
loc_cor <- loc_cor[1:nrow(grid), , drop = FALSE]
normalizing_constants <- normalizing_constants[1:nrow(grid),]
x0 <- x0[1:nrow(grid),]
f_est <- f_est[1:nrow(grid)]/
normalization_constant
density_normalized <- TRUE
} else {
normalization_constant <- NA
density_normalized <- FALSE
}
# Return
ret <- list(f_est = f_est,
loc_mean = loc_mean,
loc_sd = loc_sd,
loc_cor = loc_cor,
x = lg_object$x,
bw = lg_object$bw,
transformed_data = lg_object$transformed_data,
normalizing_constants = normalizing_constants,
grid = grid,
transformed_grid = x0,
pairs = pairs,
normalization_points = normalization_points,
normalization_constant = normalization_constant,
density_normalized = density_normalized)
# If the level-argument is provided, we also calculate the standard
# deviations based on asymptotic formulas with corresponding confidence
# intervals for the local correlations and the density estimates.
if(!is.null(level)) {
if(bootstrap) {
bootstrapped_local_correlations <- list()
for ( i in 1:B) {
lg_object_replicated <-
lg_main(x = lg_object$x[sample(1:nrow(x), replace = T), ],
bw_method = lg_object$bw_method,
est_method = lg_object$est_method,
transform_to_marginal_normality = lg_object$transform_to_marginal_normality,
bw = lg_object$bw)
bootstrapped_local_correlations[[i]] <-
dlg(lg_object_replicated,
grid = grid,
level = NULL,
normalization_points = normalization_points)$loc_cor
}
loc_cor_lower <- loc_cor_upper <- loc_cor_sd <- NA*loc_cor
# Loop over the pairs
for(i in 1:ncol(loc_cor)) {
replicated <- matrix(NA, ncol = B, nrow = nrow(loc_cor))
# Loop over the replicates to extract them from the list
for(j in 1:B) {
replicated[,j] <- bootstrapped_local_correlations[[j]][,i]
}
loc_cor_sd[, i] <- apply(replicated, 1, stats::sd)
loc_cor_lower[,i] <- apply(replicated, 1, quantile, (1 - level)/2)
loc_cor_upper[,i] <- apply(replicated, 1, quantile, 1 - (1 - level)/2)
ret$loc_cor_sd <- loc_cor_sd
ret$loc_cor_lower <- loc_cor_lower
ret$loc_cor_upper <- loc_cor_upper
}
}
if(!bootstrap & (lg_object$est_method != "trivariate")) {
# x is the data, or the transformed data. x0 is the grid
# where we do estimation, transformed or not.
# We need to evaluate the density estimate in the observations, and do
# that by running the dlg-function.
density_object_obs <- dlg(lg_object, grid = lg_object$x, level = NULL)
# We initialize a matrix for the asymptotic standard deviation for
# the local correlations. It will have the same dimension as the matrix of
# estimates:
loc_cor_sd <- NA*loc_cor
# We approximate the J- and M integrals by Monte Carlo an d the LLN, using
# locally Gaussian estimates in the observations. This has to be done for
# each pair of variables.
for(i in 1:nrow(pairs)) {
# We first calculate the J- and M integrals. The result are big
# matrices, with one column per grid point.
integrand_J <- apply(x0[, unlist(c(pairs[i,]))],
1,
function(t) mvtnorm::dmvnorm(x = (x[, unlist(c(pairs[i,]))] - t),
sigma = diag(c(lg_object$bw$joint$bw1[i]^2,
lg_object$bw$joint$bw2[i]^2)))*
u(x[,pairs[i,1]],
x[,pairs[i,2]],
c(density_object_obs$loc_cor[,i]))^2)
integrand_M1 <- apply(x0[, unlist(c(pairs[i,]))],
1,
function(t)
mvtnorm::dmvnorm(x = (x[, unlist(c(pairs[i,]))] - t),
sigma = diag(c(lg_object$bw$joint$bw1[i]^2,
lg_object$bw$joint$bw2[i]^2)))^2*
u(x[,pairs[i,1]],
x[,pairs[i,2]],
c(density_object_obs$loc_cor[,i]))^2)
integrand_M2 <- apply(x0[, unlist(c(pairs[i,]))],
1,
function(t)
mvtnorm::dmvnorm(x = (x[, unlist(c(pairs[i,]))] - t),
sigma = diag(c(lg_object$bw$joint$bw1[i]^2,
lg_object$bw$joint$bw2[i]^2)))*
u(x[,pairs[i,1]],
x[,pairs[i,2]],
c(density_object_obs$loc_cor[,i])))
# We approximate the J and the two parts of the M matrices by a Monte
# Carlo approximation:
J <- colMeans(integrand_J)
M1 <- lg_object$bw$joint$bw1[i] * lg_object$bw$joint$bw2[i] * colMeans(integrand_M1)
M2 <- lg_object$bw$joint$bw1[i] * lg_object$bw$joint$bw2[i] * colMeans(integrand_M2)^2
# The final estimate of the standard deviation will be the the expression below:
loc_cor_sd[, i] <- sqrt( (M1 - M2)/(J^2*
nrow(x)*
lg_object$bw$joint$bw1[i]*
lg_object$bw$joint$bw2[i]) )
}
# Calculate the confidence limits based on the asymptotic normality:
loc_cor_lower <- loc_cor + qnorm((1 - level)/2) * loc_cor_sd
loc_cor_upper <- loc_cor - qnorm((1 - level)/2) * loc_cor_sd
# Add the asymptotic standard deviations and confidence limits to the list that we return:
ret$loc_cor_sd <- loc_cor_sd
ret$loc_cor_lower <- loc_cor_lower
ret$loc_cor_upper <- loc_cor_upper
}
}
class(ret) <- "dlg"
return(ret)
}
#' The locally Gaussian conditional density estimator
#'
#' Estimate a conditional density function using locally Gaussian
#' approximations.
#'
#' This function is the conditional version of the locally Gaussian density
#' estimator (LGDE), described in Otneim & Tjøstheim (2018). The function takes
#' as arguments an \code{lg}-object as produced by the main \code{lg_main}- function,
#' a grid of points where the density estimate should be estimated, and a set of
#' conditions.
#'
#' The variables must be sorted before they are supplied to this function. It
#' will always assume that the free variables come before the conditioning
#' variables.
#'
#' Assume that X is a stochastic vector with two components X1 and X2. This
#' function will thus estimate the conditional density of X1 given a specified
#' value of X2.
#'
#' @param lg_object An object of type \code{lg}, as produced by the
#' \code{lg_main}-function.
#' @param grid A matrix of grid points, where we want to evaluate the density
#' estimate. Number of columns *must* be the same as number of variables in
#' X1.
#' @param condition A vector with conditions for the variables that we condition
#' upon. Length of this vector *must* be the same as the number of variables
#' in X2. The function will throw an error of there is any discrepancy in the
#' dimensions of the \code{grid}, \code{condition} and data set.
#' @param normalization_points How many grid points for approximating the integral
#' of the density estimate, to use for normalization?
#' @param fixed_grid Not used presently.
#'
#' @return A list containing the conditional density estimate as well as all the
#' running parameters that has been used. The elements are:
#'
#' \itemize{
#' \item \code{f_est}: The estimated conditional density.
#' \item \code{c_mean}: The estimated local conditional means as defined in
#' equation (10) of Otneim & Tjøstheim (2017).
#' \item \code{c_cov}: The estimated local conditional covariance matrices
#' as defined in equation (11) of Otneim & Tjøstheim (2017).
#' \item \code{x}: The data set.
#' \item \code{bw}: The bandwidth object.
#' \item \code{transformed_data}: The data transformed to approximate
#' marginal standard normality (if selected).
#' \item \code{normalizing_constants}: The normalizing constants used to
#' transform data and grid back and forth to the marginal standard
#' normality scale, as seen in eq. (8) of Otneim & Tjøstheim (2017)
#' (if selected).
#' \item \code{grid}: The grid where the estimation was performed, on the
#' original scale.
#' \item \code{transformed_grid}: The grid where the estimation was
#' performed, on the marginal standard normal scale.
#' \item \code{normalization_points} Number of grid points used
#' to approximate the integral of the density estimate, in order to
#' normalize?
#' \item \code{normalization_constant} If approximated, the integral of the
#' non-normalized density estimate. NA if not normalized.
#' \item \code{density_normalized} Logical, indicates whether the final
#' density estimate (contained in f_est) has been approximately
#' normalized to have unit integral.
#' }
#'
#' @examples
#' # A 3 variate example
#' x <- cbind(rnorm(100), rnorm(100), rnorm(100))
#'
#' # Generate the lg-object with default settings
#' lg_object <- lg_main(x)
#'
#' # Estimate the conditional density of X1|X2 = 0, X3 = 1 on a small grid
#' cond_dens <- clg(lg_object, grid = matrix(-4:4, ncol = 1), condition = c(0, 1))
#'
#' @references
#'
#' Otneim, Håkon, and Dag Tjøstheim. "Conditional density estimation using
#' the local Gaussian correlation" Statistics and Computing 28, no. 2 (2018):
#' 303-321.
#'
#' @export
clg <- function(lg_object, grid = NULL, condition = NULL,
normalization_points = NULL, fixed_grid = NULL) {
# Extract some basic info
n <- nrow(lg_object$x) # Sample size
d <- ncol(lg_object$x) # Number of variables
x <- lg_object$x
if(is.null(fixed_grid)) {
nc <- length(condition) # Number of conditioning variables
m <- nrow(grid) # Number of grid points
} else {
nc <- d - 2
m <- nrow(fixed_grid)
}
# Do some checks first
check_lg(lg_object)
if(!is.null(grid)) {
grid <- check_data(grid,
dim_check = d - nc, # The grid must have d - c columns,
type = "grid") # that is: the number of free variables
}
# Create the estimation grid, where we need the local correlation.
if(is.null(fixed_grid)) {
estimation_grid <- matrix(NA, ncol = d, nrow = m)
estimation_grid[, 1:(d-nc)] <- grid
for(i in 1:nc) {
estimation_grid[,d-nc+i] <- rep(condition[i], m)
}
} else {
estimation_grid <- fixed_grid
}
# If we are going to normalize the conditional density estimate, we
# calculate it on a grid of size normalization_points, that we generate from
# a multivariate normal having the same mean and covariance matrix as the
# independent variables of the original data. We use this to calculate the
# integral. This is only temporary, and we will remove this in the end.
if(!is.null(normalization_points)) {
extra_grid <- matrix(NA, ncol = d, nrow = normalization_points)
extra_grid[, 1:(d-nc)] <- mvtnorm::rmvnorm(n = normalization_points,
mean = colMeans(x[, 1:(d-nc), drop = FALSE]),
sigma = cov(x[, 1:(d-nc), drop = FALSE]))
for(i in 1:nc) {
extra_grid[,d-nc+i] <- rep(condition[i], normalization_points)
}
estimation_grid <- rbind(estimation_grid, extra_grid)
}
# Estimate the local correlation in these points using the density function
density_object <- dlg(lg_object, grid = estimation_grid, level = NULL)
# In each grid point, we need to calculate the conditional mean vector and
# covariance matrix in order to calculate the conditional density estimate.
# Some useful quantities:
pairs <- lg_object$bw$joint[, c("x1", "x2")] # All pairs of variables
n_pairs <- nrow(pairs) # The number of pairs
z2 <- matrix(density_object$transformed_grid[1, (d - nc + 1):d],
ncol = 1) # The transformed location
# This is a function that does this task in grid point number i.
f_eval <- function(i) {
# This grid point
grid_point <- density_object$transformed_grid[i, 1:(d - nc)]
# First, this is the conditional mean at that point
mu <- density_object$loc_mean[i,]
# First, we build up the unconditional covariance matrix
sigma <- diag(density_object$loc_sd[i,]^2)
for(j in 1:n_pairs) { # For each pair of variables
var1 <- pairs[j,1] #| The variables
var2 <- pairs[j,2] #|
rho <- density_object$loc_cor[i, j] #| The local correlation in *this*
#| grid point, and for *this* variable
#| pair.
# Finally, fill in the covariances
sigma[var1, var2] <- rho*sqrt(sigma[var1, var1]*sigma[var2, var2])
sigma[var2, var1] <- sigma[var1, var2]
}
# Next, we partition the mean vector and covariance matrix to produce the conditional
# parameters
mu1 <- matrix(mu[1:(d - nc)], ncol = 1)
mu2 <- matrix(mu[(d-nc+1):d], ncol = 1)
S11 <- as.matrix(sigma[1:(d - nc), 1:(d - nc)])
S22 <- as.matrix(sigma[(d - nc + 1):d, (d - nc + 1):d])
S12 <- as.matrix(sigma[1:(d - nc), (d - nc + 1):d])
dim(S12) <- c(d - nc, nc)
S21 <- as.matrix(sigma[(d - nc + 1):d, 1:(d - nc)])
dim(S21) <- rev(dim(S12))
# The formulas for conditional mean and conditional covariance matrix:
c_mean <- mu1 + S12%*%solve(S22)%*%(z2 - mu2)
c_cov <- S11 - S12%*%solve(S22)%*%S21
f_est <- mvtnorm::dmvnorm(x = grid_point,
mean = c_mean,
sigma = c_cov)*
prod(density_object$normalizing_constants[i, 1:(d - nc)])
return(list(f_est = f_est,
c_mean = c_mean,
c_cov = c_cov))
}
if(is.null(normalization_points)) {
# Apply the f_eval function to all the grid points
estimate <- lapply(as.list(1:m), f_eval)
} else {
# Apply the f_eval function to all the grid points
estimate <- lapply(as.list(1:(m+normalization_points)), f_eval)
}
# Extract the conditional density estimates into a vector
f_est <- unlist(lapply(estimate, "[[", 1))
# Extract the conditional means and covariance matrices into lists
c_mean <- lapply(estimate, "[[", 2)
c_cov <- lapply(estimate, "[[", 3)
# If we normalize, calculate the constant and remove the extra grid points.
# If not, set the normalizing constant equal to one.
if(!is.null(normalization_points)) {
normalization_grid <- tail(estimation_grid,
n = normalization_points)[, 1:(d-nc), drop = FALSE]
normalization_est <- tail(f_est, n = normalization_points)
normalization_constant <- mean(normalization_est/
mvtnorm::dmvnorm(normalization_grid,
mean = colMeans(x[, 1:(d-nc), drop = FALSE]),
sigma = cov(x[, 1:(d-nc), drop = FALSE])))
# Remove the extra grid points
estimation_grid <-
estimation_grid[1:(nrow(estimation_grid) - normalization_points), ]
c_cov <- c_cov[1:nrow(estimation_grid)]
c_mean <- c_mean[1:nrow(estimation_grid)]
normalizing_constants <-
density_object$normalizing_constants[1:nrow(estimation_grid),]
transformed_grid <- density_object$transformed_grid[1:nrow(estimation_grid),]
f_est <- f_est[1:nrow(grid)]/
normalization_constant
density_normalized <- TRUE
} else {
normalization_constant <- NA
density_normalized <- FALSE
normalizing_constants <-
density_object$normalizing_constants
transformed_grid <- density_object$transformed_grid
}
# Return the result
ret <- list(f_est = f_est,
c_mean = c_mean,
c_cov = c_cov,
x = lg_object$x,
bw = lg_object$bw,
transformed_data = lg_object$transformed_data,
normalizing_constants = normalizing_constants,
grid = grid,
transformed_grid = transformed_grid,
density_object = density_object,
normalization_points = normalization_points,
normalization_constant = normalization_constant,
density_normalized = density_normalized)
class(ret) <- "clg"
return(ret)
}
#' Trivariate density estimation
#'
#' \code{dlg_trivariate} returns the locally Gaussian density estimate of a
#' trivariate distribution on a given grid.
#'
#' In some applications it may be desired to produce a full locally Gaussian fit
#' of a trivariate density function without having to resort to bivariate
#' approximations. This function takes a trivariate data set, \code{x}, and a
#' trivariate set of grid points \code{eval_points}, and returns the trivariate,
#' locally Gaussian density estimate in these points. We also need a vector of
#' bandwidths, \code{bw}, with three elements, and an estimation method
#' \code{est_method}, which in this case is fixed at "trivariate", and
#' included only to be fully compatible with the other methods in this package.
#'
#' This function will only work on the marginally standard normal scale! Please
#' use the wrapper function \code{dlg()} for density estimation. This will
#' ensure that all parameters have proper values.
#'
#' @param x The data matrix (or data frame). Must have exactly 2 columns.
#' @param eval_points The grid where the density should be estimated. Must have
#' exactly 2 columns.
#' @param grid_size If \code{eval_points} is not supplied, then the function
#' will create a suitable grid diagonally through the data, with this many
#' grid points.
#' @param bw The two bandwidths, a numeric vector of length 2.
#' @param est_method The estimation method, must either be "1par" for estimation
#' with just the local correlation, or "5par" for a full locally Gaussian fit
#' with all 5 parameters.
#' @param run_checks Logical. Should sanity checks be run on the arguments?
#' Useful to disable this when doing cross-validation for example.
#'
#' @return A list including the data set \code{$x}, the grid
#' \code{$eval_points}, the bandwidths \code{$bw}, as well as a matrix of the
#' estimated parameter estimates \code{$par_est} and the estimated bivariate
#' density \code{$f_est}.
#'
#' @examples
#' x <- cbind(rnorm(100), rnorm(100), rnorm(100))
#' bw <- c(1, 1, 1)
#' eval_points <- cbind(seq(-4, 4, 1), seq(-4, 4, 1), seq(-4, 4, 1))
#'
#' estimate <- dlg_trivariate(x, eval_points = eval_points, bw = bw)
#'
#' @export
dlg_trivariate <- function(x,
eval_points = NULL,
grid_size = 15,
bw = c(1, 1, 1),
est_method = "trivariate",
run_checks = TRUE) {
# Check that everything is the way it should be -----------------------------
if(run_checks) {
x <- check_data(x, dim_check = 3, type = "data")
check_bw_trivariate(bw = bw)
check_est_method(est_method)
if(!identical(eval_points, NA))
eval_points <- check_data(eval_points, dim_check = 3, type = "grid")
}
# If eval_points not supplied, create a suitable grid
if(identical(eval_points, NULL)) {
eval_points <- apply(x,
2,
function(x) seq(quantile(x,.001),quantile(x,.999),
length.out = grid_size))
}
# We declare a function that maximizes the local likelihood in one grid point.
maximize_likelihood = function(grid_point) {
# Some quantities for the likelihood function that do not need to be
# evaluated at every iteration.
K <- 1/((2*pi)*
sqrt(2*pi)*bw[1]*bw[2]*bw[3])*
exp(-0.5*((x[,1] - grid_point[1])^2/bw[1]^2 +
(x[,2] - grid_point[2])^2/bw[2]^2 +
(x[,3] - grid_point[3])^2/bw[3]^2))
mK <- mean(K)
m1 <- mean(x[,1]^2*K)
m2 <- mean(x[,2]^2*K)
m3 <- mean(x[,3]^2*K)
m4 <- mean(x[,1]*x[,2]*K)
m5 <- mean(x[,1]*x[,3]*K)
m6 <- mean(x[,2]*x[,3]*K)
# The likelihood function
lik <- function(rho) {
w <- x[,1]^2*(rho[3]^2 - 1) +
x[,2]^2*(rho[2]^2 - 1) +
x[,3]^2*(rho[1]^2 - 1) +
2*(x[,1]*x[,2]*(rho[1] - rho[2]*rho[3]) +
x[,1]*x[,3]*(rho[2] - rho[1]*rho[3]) +
x[,2]*x[,3]*(rho[3] - rho[1]*rho[2]))
r <- 1 - rho[1]^2 - rho[2]^2 - rho[3]^2 + 2*rho[1]*rho[2]*rho[3]
d <- (1 + bw[1]^2)*(1 + bw[2]^2)*(1 + bw[3]^2) +
2*prod(rho) -
(1 + bw[1]^2)*rho[3]^2 -
(1 + bw[2]^2)*rho[2]^2 -
(1 + bw[3]^2)*rho[1]^2
exponent <-
-0.5*(
grid_point[1]^2 * ((1 + bw[2]^2)*(1 + bw[3]^2) - rho[3]^2) +
grid_point[2]^2 * ((1 + bw[1]^2)*(1 + bw[3]^2) - rho[2]^2) +
grid_point[3]^2 * ((1 + bw[1]^2)*(1 + bw[2]^2) - rho[1]^2) +
2*grid_point[1]*grid_point[2] * (rho[2]*rho[3] - rho[1]*(1 + bw[3]^2)) +
2*grid_point[2]*grid_point[3] * (rho[1]*rho[2] - rho[3]*(1 + bw[1]^2)) +
2*grid_point[1]*grid_point[3] * (rho[1]*rho[3] - rho[2]*(1 + bw[2]^2))
)/d
-(mK*(-1.5*log(2*pi)) +
mK*(-0.5*suppressWarnings(log(r))) +
0.5/r*(
(rho[3]^2 - 1)*m1 +
(rho[2]^2 - 1)*m2 +
(rho[1]^2 - 1)*m3 +
2*(rho[1] - rho[2]*rho[3])*m4 +
2*(rho[2] - rho[1]*rho[3])*m5 +
2*(rho[3] - rho[1]*rho[2])*m6
)) +
(1/(2*pi*sqrt(2*pi*d)))*exp(exponent)
}
# Return the maximum of the likelihood and the density estimate
opt <- try(optim(par = stats::cor(x)[c(2, 3, 6)],
fn = lik),
silent = TRUE)
# Store the result if the optimization went alright. Return NA if not.
if(class(opt) != "try-error") {
return(c(opt$par,
mvtnorm::dmvnorm(c(grid_point), mean = c(0,0,0),
sigma = matrix(c(1, opt$par[1], opt$par[2],
opt$par[1], 1, opt$par[3],
opt$par[2], opt$par[3], 1), 3))))
} else {
return(rep(NA, 4))
}
}
## Send the grid points to 'maximize_likelihood'
est <- cbind(do.call(rbind,
lapply(X = split(eval_points, row(eval_points)),
FUN = maximize_likelihood)))
par_est <- matrix(est[,1:3], ncol = 3)
f_est = as.vector(est[,4])
colnames(par_est) <- c('rho12', 'rho13', 'rho23')
return(list(x = x,
eval_points = eval_points,
bw = bw,
par_est = par_est,
f_est = f_est))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.