Skip to content

Add bayesEL part to jss draft #5

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
yeoveryeah wants to merge 6 commits into
base: develop-jss
Choose a base branch
from
Open

Add bayesEL part to jss draft #5

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
ebfeabc
draft to jss-essay
yeoveryeah Nov 28, 2021
c9bb0e6
Merge pull request #4 from mlysy/develop-jss
mlysy Dec 3, 2021
31d4a6d
update with Bayesian EL parth
yeoveryeah Dec 12, 2021
27141dc
jss draft
yeoveryeah Dec 19, 2021
ff245cc
added convergence check to GenEL
mlysy Dec 21, 2021
473dde2
Merge branch 'main' of https://github.com/mlysy/flexEL into develop-jss
yeoveryeah Jun 19, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion DESCRIPTION
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,7 @@ Description: Various tools for implementing and calibrating empirical likelihood
License: GPL (>= 3)
Imports: Rcpp (>= 0.12.13), numDeriv
LinkingTo: Rcpp, RcppEigen
RoxygenNote: 7.1.1
RoxygenNote: 7.1.2
Encoding: UTF-8
Suggests: testthat, MASS, optimCheck, knitr, rmarkdown, RcppEigen
VignetteBuilder: knitr
Expand Down
24 changes: 14 additions & 10 deletions R/RcppExports.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -89,28 +89,32 @@ GenEL_get_supp_adj <- function(pGEL) {
.Call(`_flexEL_GenEL_get_supp_adj`, pGEL)
}

GenEL_lambda_nr <- function(pGEL, G, weights, verbose) {
.Call(`_flexEL_GenEL_lambda_nr`, pGEL, G, weights, verbose)
GenEL_get_diag <- function(pGEL) {
.Call(`_flexEL_GenEL_get_diag`, pGEL)
}

GenEL_lambda_nr <- function(pGEL, G, weights, check_conv) {
.Call(`_flexEL_GenEL_lambda_nr`, pGEL, G, weights, check_conv)
}

GenEL_omega_hat <- function(pGEL, lambda, G, weights) {
.Call(`_flexEL_GenEL_omega_hat`, pGEL, lambda, G, weights)
}

GenEL_logel <- function(pGEL, G, verbose) {
.Call(`_flexEL_GenEL_logel`, pGEL, G, verbose)
GenEL_logel <- function(pGEL, G, check_conv) {
.Call(`_flexEL_GenEL_logel`, pGEL, G, check_conv)
}

GenEL_weighted_logel <- function(pGEL, G, weights, verbose) {
.Call(`_flexEL_GenEL_weighted_logel`, pGEL, G, weights, verbose)
GenEL_weighted_logel <- function(pGEL, G, weights, check_conv) {
.Call(`_flexEL_GenEL_weighted_logel`, pGEL, G, weights, check_conv)
}

GenEL_logel_grad <- function(pGEL, G, verbose) {
.Call(`_flexEL_GenEL_logel_grad`, pGEL, G, verbose)
GenEL_logel_grad <- function(pGEL, G, check_conv) {
.Call(`_flexEL_GenEL_logel_grad`, pGEL, G, check_conv)
}

GenEL_weighted_logel_grad <- function(pGEL, G, weights, verbose) {
.Call(`_flexEL_GenEL_weighted_logel_grad`, pGEL, G, weights, verbose)
GenEL_weighted_logel_grad <- function(pGEL, G, weights, check_conv) {
.Call(`_flexEL_GenEL_weighted_logel_grad`, pGEL, G, weights, check_conv)
}

MeanReg_evalG <- function(y, X, beta) {
Expand Down
121 changes: 67 additions & 54 deletions R/class_GenEL.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,14 +23,25 @@ GenEL <- R6::R6Class(
n_obs <- GenEL_get_n_obs(private$.GEL)
n_eqs <- GenEL_get_n_eqs(private$.GEL)
if (nrow(G) != n_obs | ncol(G) != n_eqs) {
stop("Dimension of G matrix must be `n_obs` x `n_eqs`.")
stop("Dimension of `G` matrix must be `n_obs` x `n_eqs`.")
}
},

#' @description Check the dimension of weights is as expected.
#' @param weights A numeric vector.
check_weights = function(weights) {
if (length(weights) != GenEL_get_n_obs(private$.GEL)) {
stop("Length of `weights` does not equal to the number of obserations.")
}
if (any(weights < 0)) {
stop("`weights` should contain only non-negative values.")
}
}
),

active = list(

#' @field max_iter Access or reset the value of the maximum number of
#' @field max_iter Access or reset the value of the maximum number of
#' iterations. The value must be a positive integer (default to 100).
max_iter = function(value) {
if (missing(value)) private$.max_iter
Expand All @@ -43,8 +54,8 @@ GenEL <- R6::R6Class(
}
},

#' @field rel_tol Access or reset the value of relative tolerance
#' controlling the accuracy at convergence. The value must be a small
#' @field rel_tol Access or reset the value of relative tolerance
#' controlling the accuracy at convergence. The value must be a small
#' positive number (default to 1e-7).
rel_tol = function(value) {
if (missing(value)) private$.rel_tol
Expand All @@ -57,7 +68,7 @@ GenEL <- R6::R6Class(
}
},

#' @field lambda0 Access or reset the initial value of lambda. The value
#' @field lambda0 Access or reset the initial value of lambda. The value
#' must be a vector of length `n_eqs` (default to a vector of 0).
lambda0 = function(value) {
if (missing(value)) private$.lambda0
Expand All @@ -71,8 +82,8 @@ GenEL <- R6::R6Class(
}
},

#' @field supp_adj Access or reset the support correction flag. The value
#' must be a boolean indicating whether to conduct support correction or
#' @field supp_adj Access or reset the support correction flag. The value
#' must be a boolean indicating whether to conduct support correction or
#' not (default to FALSE).
supp_adj = function(value) {
if (missing(value)) private$.supp_adj
Expand All @@ -87,7 +98,7 @@ GenEL <- R6::R6Class(
}
},

#' @field supp_adj_a Access or reset the value of support corection factor.
#' @field supp_adj_a Access or reset the value of support corection factor.
#' The value must be a scalar (defaults to `max(1.0, log(n_obs)/2)`).
supp_adj_a = function(value) {
if (missing(value)) private$.supp_adj_a
Expand All @@ -103,8 +114,8 @@ GenEL <- R6::R6Class(
}
},

#' @field weight_adj Access or reset the value of the weight for the
#' additional fake observation under support correction. The value must
#' @field weight_adj Access or reset the value of the weight for the
#' additional fake observation under support correction. The value must
#' be a scalar (default to NULL).
weight_adj = function(value) {
if (missing(value)) private$.weight_adj
Expand Down Expand Up @@ -150,99 +161,101 @@ GenEL <- R6::R6Class(
#' @param supp_adj A boolean indicating whether to conduct support correction or not.
#' @param supp_adj_a Support adjustment factor. Defaults to `max(1.0, log(n_obs)/2)`.
#' @param weight_adj Weight under weighted log EL (default to 1.0).
#' @param lambda0 Initialization vector of size `n_eqs`.
#' @param lambda0 Initialization vector of size `n_eqs`. Defaults to a vector of zeros.
set_opts = function(max_iter = 100, rel_tol = 1e-7,
supp_adj = FALSE, supp_adj_a = NULL, weight_adj = NULL,
lambda0 = rep(0, GenEL_get_n_eqs(private$.GEL))) {
lambda0 = NULL) {
self$max_iter <- max_iter
self$rel_tol <- rel_tol
self$set_supp_adj(supp_adj = supp_adj, supp_adj_a = supp_adj_a, weight_adj = weight_adj)
if(is.null(lambda0)) lambda0 <- rep(0, GenEL_get_n_eqs(private$.GEL))
self$lambda0 <- lambda0
},

#' @description Calculate the solution of the dual problem of maximum log EL problem.
#' @param G A numeric matrix of dimension `n_obs x n_eqs`.
#' @param weights A numeric vector of length `n_obs` containing non-negative values.
#' @param verbose A boolean indicating whether to print out number of iterations and maximum error at the end of the Newton-Raphson algorithm.
#' @param check_conv If `TRUE`, checks whether the desired `rel_tol` was reached within the given maximum number of Newton-Raphson iterations. If not, returns a vector of `NaN`s.
#' @return A numeric vector of length `n_eqs`.
lambda_nr = function(G, weights, verbose = FALSE) {
lambda_nr = function(G, weights, check_conv = TRUE) {
private$check_G(G)
n_obs <- GenEL_get_n_obs(private$.GEL)
if (missing(weights) || is.null(weights)) {
weights <- rep(1.0, n_obs)
}
if (length(weights) != n_obs) {
stop("Length of `weights` does not equal to the number of obserations.")
}
if (any(weights < 0)) {
stop("`weights` should contain only non-negative values.")
}
GenEL_lambda_nr(private$.GEL, t(G), weights, verbose)
private$check_weights(weights)
GenEL_lambda_nr(private$.GEL, t(G), weights, check_conv)
},

#' @description Calculate the probability vector base on the given G matrix.
#' @param G A numeric matrix of dimension `n_obs x n_eqs`.
#' @param weights A numeric vector of length `n_obs` containing non-negative values.
#' @param verbose A boolean indicating whether to print out number of iterations and maximum error at the end of the Newton-Raphson algorithm.
#' @param check_conv If `TRUE`, checks whether the desired `rel_tol` was reached within the given maximum number of Newton-Raphson iterations. If not, returns a vector of `NaN`s.
#' @return A probability vector of length `n_obs + supp_adj`.
omega_hat = function(G, weights, verbose = FALSE) {
omega_hat = function(G, weights, check_conv = TRUE) {
private$check_G(G)
n_obs <- GenEL_get_n_obs(private$.GEL)
if (missing(weights) || is.null(weights)) {
weights <- rep(1.0, n_obs)
}
if (length(weights) != n_obs) {
stop("Length of `weights` does not equal to the number of obserations.")
}
if (any(weights < 0)) {
stop("`weights` should contain only non-negative values.")
}
lambda <- self$lambda_nr(G = G, weights = weights, verbose = verbose)
## if (length(weights) != n_obs) {
## stop("Length of `weights` does not equal to the number of obserations.")
## }
## if (any(weights < 0)) {
## stop("`weights` should contain only non-negative values.")
## }
private$check_weights(weights)
lambda <- self$lambda_nr(G = G, weights = weights,
check_conv = check_conv)
GenEL_omega_hat(private$.GEL, lambda, t(G), weights)
},

#' @description Calculate the log empirical likelihood base on the given G matrix.
#' @param G A numeric matrix of dimension `n_eqs x n_obs`.
#' @param weights A numeric vector of length `n_obs` containing non-negative values.
#' @param verbose A boolean indicating whether to print out number of iterations and maximum error at the end of the Newton-Raphson algorithm.
#' @param check_conv If `TRUE`, checks whether the desired `rel_tol` was reached within the given maximum number of Newton-Raphson iterations. If not, returns `-Inf`.
#' @return A scalar.
logel = function(G, weights, verbose = FALSE) {
logel = function(G, weights, check_conv = TRUE) {
private$check_G(G)
if (missing(weights) || is.null(weights)) {
GenEL_logel(private$.GEL, t(G), verbose)
ans <- GenEL_logel(private$.GEL, t(G), check_conv)
}
else {
n_obs <- GenEL_get_n_obs(private$.GEL)
if (length(weights) != n_obs) {
stop("Length of `weights` does not equal to the number of obserations.")
}
if (any(weights < 0)) {
stop("`weights` should contain only non-negative values.")
}
GenEL_weighted_logel(private$.GEL, t(G), weights, verbose)
## n_obs <- GenEL_get_n_obs(private$.GEL)
## if (length(weights) != n_obs) {
## stop("Length of `weights` does not equal to the number of obserations.")
## }
## if (any(weights < 0)) {
## stop("`weights` should contain only non-negative values.")
## }
private$check_weights(weights)
GenEL_weighted_logel(private$.GEL, t(G), weights, check_conv)
}
},

#' @description Calculate log EL and the derivative of log EL w.r.t. G evaluated at G.
#' @param G A numeric matrix of dimension `n_obs x n_eqs`.
#' @param weights A numeric vector of length `n_obs` containing non-negative values.
#' @param verbose A boolean indicating whether to print out number of iterations and maximum error at the end of the Newton-Raphson algorithm.
#' @param check_conv If `TRUE`, checks whether the desired `rel_tol` was reached within the given maximum number of Newton-Raphson iterations. If not, sets the log EL to negative infinity and the gradient to a matrix of `NaN`s.
#' @return A list of two elements.
logel_grad = function(G, weights, verbose = FALSE) {
logel_grad = function(G, weights, check_conv = TRUE) {
private$check_G(G)
if (missing(weights) || is.null(weights)) {
GenEL_logel_grad(private$.GEL, t(G), verbose)
ans <- GenEL_logel_grad(private$.GEL, t(G), check_conv)
}
else {
n_obs <- GenEL_get_n_obs(private$.GEL)
if (length(weights) != n_obs) {
stop("Length of `weights` does not equal to the number of obserations.")
}
if (any(weights < 0)) {
stop("`weights` should contain only non-negative values.")
}
GenEL_weighted_logel_grad(private$.GEL, t(G), weights, verbose)
## n_obs <- GenEL_get_n_obs(private$.GEL)
## if (length(weights) != n_obs) {
## stop("Length of `weights` does not equal to the number of obserations.")
## }
## if (any(weights < 0)) {
## stop("`weights` should contain only non-negative values.")
## }
private$check_weights(weights)
ans <- GenEL_weighted_logel_grad(private$.GEL, t(G), weights, check_conv)
}
ans$dldG <- t(ans$dldG)
ans
}
)
)
8 changes: 8 additions & 0 deletions inst/include/flexEL/gen_el.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -126,6 +126,8 @@ namespace flexEL {
const Ref<const MatrixXd>& G,
const Ref<const VectorXd>& weights);
// const Ref<const VectorXd>& norm_weights);
/// Check whether Newton-Raphson algorithm has converged.
bool has_converged_nr();
/// Calculate the profile probability weights.
void omega_hat(Ref<VectorXd> omega,
const Ref<const VectorXd>& lambda,
Expand Down Expand Up @@ -249,6 +251,12 @@ namespace flexEL {
return;
}

/// @return Whether the desired error tolerance `rel_tol` has been reached in less than `max_iter` Newton-Raphson steps.
/// @warning Assumes that lambda_nr has been run at least once.
inline bool GenEL::has_converged_nr() {
return (nr_err_ > rel_tol_) && (nr_iter_ == max_iter_);
}

/// @param[in] G Moment matrix of size `n_eqs x n_obs` or `n_eqs x (n_obs + 1)`. If `supp_adj = false`, the former is required. If `supp_adj = true` and the former is provided, support adjustment is performed. If `supp_adj = true` and `G.cols() == n_obs + 1`, assumes that support has already been corrected.
inline Ref<const MatrixXd> GenEL::supp_G(const Ref<const MatrixXd>& G) {
if(supp_adj_ && G.cols() == n_obs_) {
Expand Down
Loading