EL: Self-concordant multi-variate empirical likelihood with...

In smoothemplik: Smoothed Empirical Likelihood

Self-concordant multi-variate empirical likelihood with counts

Description

Implements the empirical-likelihood-ratio test for the mean of the coordinates of z (with the hypothesised value mu). The counts need not be integer; in the context of local likelihoods, they can be kernel observation weights.

Usage

EL(
  z,
  mu = NULL,
  ct = NULL,
  lambda.init = NULL,
  SEL = FALSE,
  return.weights = FALSE,
  lower = NULL,
  upper = NULL,
  order = 4L,
  weight.tolerance = NULL,
  thresh = 1e-16,
  itermax = 100L,
  verbose = FALSE,
  alpha = 0.3,
  beta = 0.8,
  backeps = 0
)

Arguments

z	A numeric vector or a matrix with one data vector per column.
mu	Hypothesised mean, default (0 ... 0) in R^ncol(z)
ct	A numeric vector of non-negative counts.
lambda.init	Starting lambda, default (0 ... 0)
SEL	If FALSE, the default weight tolerance is MachEps^(1/3), otherwise it is MachEps^(1/2) of the maximum count.
return.weights	Logical: if TRUE, returns the empirical probabilities. Default is memory-saving (FALSE).
lower	Lower cut-off for [logTaylor()], default 1/nrow(z)
upper	Upper cutoff for [logTaylor()], default Inf
order	Positive integer such that the Taylor approximation of this order to log(x) is self-concordant; usually 4 or higher. Passed to [logTaylor()].
weight.tolerance	Weight tolerance for counts to improve numerical stability
thresh	Convergence threshold for log-likelihood (the default is aggressive)
itermax	Upper bound on number of Newton steps (seems ample)
verbose	Logical: print output diagnostics?
alpha	Backtracking line search parameter: acceptance of a decrease in function value by ALPHA*f of the prediction based on the linear extrapolation.
beta	Backtracking line search reduction factor. 0.1 corresponds to a very crude search, 0.8 corresponds to a less crude search.
backeps	Backtrack threshold: the search can miss by this much. Consider setting it to 1e-10 if backtracking seems to be failing due to round-off.

Details

Negative weights are not allowed. They could be useful in some applications, but they can destroy convexity or even boundedness. They also make the Newton step fail to be of least squares type.

This function relies on the improved computational strategy for the empirical likelihood. The search of the lambda multipliers is carried out via a dampened Newton method with guaranteed convergence owing to the fact that the log-likelihood is replaced by its Taylor approximation of any desired order (default: 4, the minimum value that ensures self-concordance).

Tweak alpha and beta with extreme caution. See \insertCiteboyd2004convexsmoothemplik, pp. 464–466 for details. It is necessary that 0 < alpha < 1/2 and 0 < beta < 1. alpha = 0.3 seems better than 0.01 on some 2-dimensional test data (sometimes fewer iterations).

The argument names, except for lambda.init, are matching the original names in Art B. Owen's implementation. The highly optimised one-dimensional counterpart, EL0, is designed to return a faster and a more accurate solution in the one-dimensional case.

Value

A list with the following values:

logelr

Log of empirical likelihood ratio (equal to 0 if the hypothesised mean is equal to the sample mean)

lam

Vector of Lagrange multipliers

wts

Observation weights/probabilities (vector of length n)

converged

TRUE if algorithm converged. FALSE usually means that mu is not in the convex hull of the data. Then, a very small likelihood is returned (instead of zero).

iter

Number of iterations taken.

ndec

Newton decrement (see Boyd & Vandenberghe).

gradnorm

Norm of the gradient of log empirical likelihood.

Source

This original code was written for \insertCiteowen2013selfsmoothemplik and [published online](https://artowen.su.domains/empirical/) by Art B. Owen (March 2015, February 2017). The present version was rewritten in Rcpp and slightly reworked to contain fewer inner functions and loops.

References

\insertAllCited

See Also

[logTaylor()], [EL0()]

Examples

earth <- c(
  5.5, 5.61, 4.88, 5.07, 5.26, 5.55, 5.36, 5.29, 5.58, 5.65, 5.57, 5.53, 5.62, 5.29,
  5.44, 5.34, 5.79, 5.1, 5.27, 5.39, 5.42, 5.47, 5.63, 5.34, 5.46, 5.3, 5.75, 5.68, 5.85
)
EL(earth, mu = 5.517, verbose = TRUE) # 5.517 is the modern accepted value

# Linear regression through empirical likelihood
coef.lm <- coef(lm(mpg ~ hp + am, data = mtcars))
xmat <- cbind(1, as.matrix(mtcars[, c("hp", "am")]))
yvec <- mtcars$mpg
foc.lm <- function(par, x, y) {  # The sample average of this
  resid <- y - drop(x %*% par)   # must be 0
  resid * x
}
minusEL <- function(par) -EL(foc.lm(par, xmat, yvec), itermax = 10)$logelr
coef.el <- optim(c(mean(yvec), 0, 0), minusEL)$par
abs(coef.el - coef.lm) / coef.lm  # Relative difference

# Likelihood ratio testing without any variance estimation
# Define the profile empirical likelihood for the coefficient on am
minusPEL <- function(par.free, par.am)
  -EL(foc.lm(c(par.free, par.am), xmat, yvec), itermax = 20)$logelr
# Constrained maximisation assuming that the coef on par.am is 3.14
coef.el.constr <- optim(coef.el[1:2], minusPEL, par.am = 3.14)$par
print(-2 * EL(foc.lm(c(coef.el.constr, 3.14), xmat, yvec))$logelr)
# Exceeds the critical value qchisq(0.95, df = 1)

smoothemplik documentation built on Aug. 22, 2025, 1:11 a.m.