# R/EL.R In EL: Two-Sample Empirical Likelihood

#### Documented in EL.HuberEL.meansEL.plotEL.smoothEL.statistic

```EL.statistic <- function(method, X, Y, d, t, bw = bw.nrd0)
{
EL.local(ELLR, X, Y, d, t, method=method, bw = bw)
}

EL.plot <- function(method, X, Y, bw = bw.nrd0,
conf.level=NULL, simultaneous=FALSE,
bootstrap.samples=300, more.warnings=FALSE, ...)
{
numpoints <- 50
call <- match.call(expand.dots = FALSE)
optional <- call\$...
switch(match(method, c("fdiff", "qdiff", "qq", "pp", "roc")),
{
dxlab <- expression(x)
dylab <- substitute(expression(F[a]-F[b]), list(a = call\$Y, b = call\$X))
},
{
dxlab <- "P"
dylab <- substitute(expression(F[a]^-1-F[b]^-1), list(a = call\$Y, b = call\$X))
},
{
dxlab <- substitute(expression(F[a]^-1), list(a = call\$Y))
dylab <- substitute(expression(F[a]^-1), list(a = call\$X))
},
{
dxlab <- substitute(expression(F[a]), list(a = call\$Y))
dylab <- substitute(expression(F[a]), list(a = call\$X))
},
{
dxlab <- "False positive rate"
dylab <- "True positive rate"
}
)
tr <- EL.local(function(globals, X, Y) globals\$trange(globals, X, Y), X, Y, method=method, bw=bw)
if (simultaneous)
{
trlen <- tr[2] - tr[1]
tr[1] <- tr[1] + 0.05 * trlen
tr[2] <- tr[2] - 0.05 * trlen
}
t <- seq(tr[1], tr[2], length=numpoints)
zz <- EL.smooth(method=method, X, Y, t, bw=bw, conf.level = conf.level,
simultaneous = simultaneous, bootstrap.samples=bootstrap.samples, more.warnings=more.warnings)
if (!is.null(conf.level))
dylim <- c(min(zz\$conf.int[1,]), max(zz\$conf.int[2,]))
else
dylim <- c(min(zz\$estim), max(zz\$estim))

if (!("xlab" %in% names(optional)))
optional <- c(optional, list(xlab=dxlab))
if (!("ylab" %in% names(optional)))
optional <- c(optional, list(ylab=dylab))
if (!("ylim" %in% names(optional)))
optional <- c(optional, list(ylim=dylim))
eval(as.call(c(as.list(quote(plot(t, zz\$estim, type='l'))), optional)))
if (!is.null(conf.level))
{
eval(as.call(c(as.list(quote(lines(t, zz\$conf.int[1,], lty="dashed"))), optional)))
eval(as.call(c(as.list(quote(lines(t, zz\$conf.int[2,], lty="dashed"))), optional)))
}
}

EL.smooth <- function(method, X, Y, t, bw = bw.nrd0,
conf.level=NULL, simultaneous=FALSE, bootstrap.samples=300,
more.warnings = FALSE)
{
warnings <- FALSE
mm <- function(globals, tt)
{
d <- sapply(tt, function(t) deltasolve(globals, X, Y, t))
bs <- NA
conffun <- function(conf.level)
{
function(k) withCallingHandlers(ELconf(globals, X, Y, tt[k], conf.level, delta=d[k]),
simpleWarning = function(e)
{
if (!more.warnings)
{
warnings <<- TRUE
invokeRestart("muffleWarning")
}
})
}
if (is.null(conf.level))
ci <- NA
else
{
if (!simultaneous)
{
ci <- sapply(1:length(tt), conffun(conf.level))
attr(ci, "conf.level") <- conf.level
}
else
{
bootup <- function()
{
X1 <- sample(X, replace=TRUE)
Y1 <- sample(Y, replace=TRUE)
max(mapply(ELLR, delta=d, t=tt, MoreArgs=list(globals=globals, X=X1, Y=Y1)))
}
bs <- sort(replicate(bootstrap.samples, bootup()))
bs <- bs[ceiling(conf.level * bootstrap.samples)]
ci <- sapply(1:length(tt), conffun(pchisq(bs, 1)))
}
}

list(estimate=d, conf.int=ci, simultaneous.conf.int=simultaneous, bootstrap.crit = bs)
}
res <- EL.local(mm, t, method=method, bw = bw)
if (warnings)
warning("Estimates or confidence bands for some values could not be found. For a detailed report, run this command again with 'more.warnings = TRUE'.", call. = FALSE)
res
}

EL.means <- function(X, Y, mu = 0, conf.level=0.95)
{
call <- match.call()
mm <- function(globals)
{
d <- deltasolve(globals, X, Y, 0)
names(d) <- "Mean difference"
ci <- ELconf(globals, X, Y, 0, conf.level, delta=d)
attr(ci, "conf.level") <- conf.level
ll <- ELLR(globals, X, Y, mu, 0)
names(ll) <- "-2 * LogLikelihood"
res <- list(estimate=d, conf.int=ci, p.value = 1 - pchisq(ll, 1), statistic = ll,
method="Empirical likelihood mean difference test",
null.value=mu, data.name = paste(deparse(call\$X), " and ", deparse(call\$Y)))
class(res) <- "htest"
res
}
EL.local(mm, method="mean")
}

EL.Huber <- function(X, Y, mu=0, conf.level=0.95, scaleX=1, scaleY=1, VX=2.046, VY=2.046, k=1.35)
{
call <- match.call()
mm <- function(globals)
{
globals\$Hsigma1 <- scaleX
globals\$Hsigma2 <- scaleY
globals\$Hconst <- k
globals\$HV1 <- VX
globals\$HV2 <- VY
d <- deltasolve(globals, X, Y, 0)
names(d) <- "Huber estimator difference"
ci <- ELconf(globals, X, Y, 0, conf.level, delta=d)
ll <- ELLR(globals, X, Y, mu, 0)
names(ll) <- "-2 * LogLikelihood"
attr(ci, "conf.level") <- conf.level
res <- list(estimate=d, conf.int=ci, p.value = 1 - pchisq(ll, 1), statistic = ll,
method="Empirical likelihood Huber estimator difference test",
null.value=mu, data.name = paste(deparse(call\$X), " and ", deparse(call\$Y)))
class(res) <- "htest"
res
}
EL.local(mm, method="huber")
}

#### Globals setup

set.globals <- function(globals)
{
globals\$lambda1 <- 0
globals\$lambda2 <- 0
globals\$wval1 <- NULL
globals\$wval2 <- NULL
globals\$lambdarange1 <- 0
globals\$lambdarange2 <- 0
globals\$degenerate <- FALSE
globals\$precision <- 1e-4
globals
}

#### Bandwidth setup

set.bw <- function(globals, bw)
{
globals\$bw.X <- 0
globals\$bw.Y <- 0
globals\$bw <- bw
globals\$initbw <- bw.init
globals
}

bw.init <- function(globals, X, Y)
{
if ((globals\$bw.X != 0) && (globals\$bw.Y != 0))
return(globals)

if (is.function(globals\$bw))
{
globals\$bw.X <- globals\$bw(X)
globals\$bw.Y <- globals\$bw(Y)
}
else
{
globals\$bw.X <- globals\$bw[1]
globals\$bw.Y <- globals\$bw[2]
}
globals
}

#### Kernel setup

set.kernel <- function(globals)
{
globals\$krnl <- kernel.norm
globals\$dkrnl <- kernel.dnorm
globals\$ikrnl <- kernel.inorm
globals
}

kernel.norm <- function(t)
pnorm(t)

kernel.dnorm <- function(t)
dnorm(t)

kernel.inorm <- function(t)
qnorm(t)

#### Method setup

set.method <- function(globals, type)
{
if (!(type %in% c("pp", "qq", "roc", "mean", "fdiff", "qdiff", "huber")))
stop(paste("Unknown EL method: ", type, "\n",
"Recognized EL methods are: pp, qq, roc, mean, fdiff, qdiff."))

globals\$w1 <- get(paste("w1.", type, sep=""))
globals\$w2 <- get(paste("w2.", type, sep=""))
globals\$alpha1 <- get(paste("alpha1.", type, sep=""))
globals\$alpha2 <- get(paste("alpha2.", type, sep=""))
globals\$trange <- get(paste("trange.", type, sep=""))
globals\$deltarange <- get(paste("deltarange.", type, sep=""))
globals\$thetarange <- get(paste("thetarange.", type, sep=""))
globals\$deltarange.opt <- get(paste("deltarange.opt.", type, sep=""))
globals\$thetarange.opt <- get(paste("thetarange.opt.", type, sep=""))
globals\$use.smoothing <- get(paste("use.smoothing.", type, sep=""))
globals
}

### Mean
w1.mean <- function(globals, X, theta, delta, t)
X - theta

w2.mean <- function(globals, X, theta, delta, t)
X - theta + delta

alpha1.mean <- function(globals, X, theta, delta, t)
-1

alpha2.mean <- function(globals, X, theta, delta, t)
-1

trange.mean <- function(globals, X, Y)
stop("Mean differences are not dependent on the parameter t.")

deltarange.mean <- function(globals, X, Y, t)
{
tr <- c(min(X), max(X))
tol <- 0.01 * (tr[2] - tr[1])
c(tr[1] -max(Y) + tol, tr[2]-min(Y)- tol)
}

deltarange.opt.mean <- deltarange.mean

thetarange.mean <- function(globals, X, Y, delta, t)
{
tr <- c(min(X), max(X))
dr <- c(min(Y) + delta, max(Y) + delta)
tol <- 0.01 * (tr[2] - tr[1])
c(max(tr[1], dr[1]) + tol, min(tr[2], dr[2]) - tol)
}

thetarange.opt.mean <- thetarange.mean

use.smoothing.mean <- FALSE

### Huber estimator

w1.huber <- function(globals, X, theta, delta, t)
{
globals\$Hsigma <- globals\$Hsigma1
globals\$HV <- globals\$HV1
Hsmooth(globals, (X - theta) / globals\$Hsigma)
}

w2.huber <- function(globals, X, theta, delta, t)
{
globals\$Hsigma <- globals\$Hsigma2
globals\$HV <- globals\$HV2
Hsmooth(globals, (X - theta + delta) / globals\$Hsigma)
}

alpha1.huber <- function(globals, X, theta, delta, t)
{
globals\$Hsigma <- globals\$Hsigma1
globals\$HV <- globals\$HV1
- DHsmooth(globals, (X - theta) / globals\$Hsigma) / globals\$Hsigma
}

alpha2.huber <- function(globals, X, theta, delta, t)
{
globals\$Hsigma <- globals\$Hsigma2
globals\$HV <- globals\$HV2
- DHsmooth(globals, (X - theta + delta) / globals\$Hsigma) / globals\$Hsigma
}

trange.huber <- function(globals, X, Y)
stop("Huber estimator differences are not dependent on the parameter t.")

deltarange.huber <- function(globals, X, Y, t)
{
tr <- c(min(X), max(X))
c(tr[1] -max(Y) + globals\$precision, tr[2]-min(Y)- globals\$precision)
}

deltarange.huber <- deltarange.mean

deltarange.opt.huber <- deltarange.huber

thetarange.huber <- function(globals, X, Y, delta, t)
{
tr <- c(min(X) + globals\$precision, max(X) - globals\$precision)
dr <- c(min(Y) + delta, max(Y) + delta)
c(max(tr[1], dr[1]-globals\$precision), min(tr[2], dr[2]+globals\$precision))
}

thetarange.huber <- thetarange.mean

thetarange.opt.huber <- thetarange.huber

use.smoothing.huber <- FALSE

Hsmooth <- function(globals, x)
{
x.l <- (x - globals\$Hconst) / globals\$HV
x.r <- (x + globals\$Hconst) / globals\$HV
pnorm(x.l) * (globals\$Hconst - x) + pnorm(x.r) * (globals\$Hconst + x) -
globals\$Hconst + globals\$HV * (dnorm(x.r) - dnorm(x.l))
}

DHsmooth <- function(globals, x)
{
x.l <- (x - globals\$Hconst) / globals\$HV
x.r <- (x + globals\$Hconst) / globals\$HV
pnorm(x.r) - pnorm(x.l)
}

### ROC

w1.roc <- function(globals, X, theta, delta, t)
globals\$krnl((theta - X) / globals\$bw.X) - (1 - delta)

w2.roc <- function(globals, X, theta, delta, t)
globals\$krnl((theta - X) / globals\$bw.Y) - (1 - t)

alpha1.roc <- function(globals, X, theta, delta, t)
globals\$dkrnl((theta - X) / globals\$bw.X) / globals\$bw.X

alpha2.roc <- function(globals, X, theta, delta, t)
globals\$dkrnl((theta - X) / globals\$bw.Y) / globals\$bw.Y

trange.roc <- function(globals, X, Y)
{
c(0.01, 0.99)
}

deltarange.roc <- function(globals, X, Y, t)
{
dY <- globals\$bw.Y * globals\$ikrnl(t)
diffs <- c(max(X) - min(Y) + dY, min(X) - max(Y) + dY)
range <- globals\$krnl(diffs / globals\$bw.X)
c(min(range)+globals\$precision, max(range)-globals\$precision)
}

deltarange.opt.roc <- function(globals, X, Y, t)
{
kk <- 1 - ecdf(X)(c(quantile(Y, 1-t) - 3 * globals\$bw.X, quantile(Y, 1-t) + 3 * globals\$bw.X))
c(max(0.001, kk[2] - 0.001), min(0.999, kk[1] + 0.001))
}

thetarange.roc <- function(globals, X, Y, delta, t)
{
dX <- globals\$bw.X * globals\$ikrnl(delta)
dY <- globals\$bw.Y * globals\$ikrnl(t)
dd <- c(max(min(X) - dX, min(Y) - dY), min(max(X) - dX, max(Y) - dY))
c(dd[1] + 0.01 * (dd[2] - dd[1]), dd[2] - 0.01 * (dd[2] - dd[1]))
}

thetarange.opt.roc <- function(globals, X, Y, delta, t)
{
extreme <- thetarange.roc(globals, X, Y, delta, t)

tol <- 2 * globals\$bw.Y
qY <- quantile(Y, 1-t)
c(max(extreme[1], qY - tol), min(extreme[2], qY + tol))
}

use.smoothing.roc <- TRUE

### P-P plots

w1.pp <- function(globals, X, theta, delta, t)
globals\$krnl((theta - X) / globals\$bw.X) - delta

w2.pp <- function(globals, X, theta, delta, t)
globals\$krnl((theta - X) / globals\$bw.Y) - t

alpha1.pp <- function(globals, X, theta, delta, t)
globals\$dkrnl((theta - X) / globals\$bw.X) / globals\$bw.X

alpha2.pp <- function(globals, X, theta, delta, t)
globals\$dkrnl((theta - X) / globals\$bw.Y) / globals\$bw.Y

trange.pp <- function(globals, X, Y)
{
c(0.01, 0.99)
}

deltarange.pp <- function(globals, X, Y, t)
{
c(0.001, 0.999)
}

deltarange.opt.pp <- function(globals, X, Y, t)
{
kk <- ecdf(X)(c(quantile(Y, t) - 3 * globals\$bw.X, quantile(Y, t) + 3 * globals\$bw.X))
c(max(0.001, kk[1] - 0.001), min(0.999, kk[2] + 0.001))
}

thetarange.pp <- function(globals, X, Y, delta, t)
{
dX <- globals\$bw.X * globals\$ikrnl(delta)
dY <- globals\$bw.Y * globals\$ikrnl(t)
dd <- c(max(min(X) + dX, min(Y) + dY), min(max(X) + dX, max(Y) + dY))
c(dd[1] + 0.01 * (dd[2] - dd[1]), dd[2] - 0.01 * (dd[2] - dd[1]))
}

thetarange.opt.pp <- function(globals, X, Y, delta, t)
{
dX <- globals\$bw.X * globals\$ikrnl(delta)
dY <- globals\$bw.Y * globals\$ikrnl(t)
dd <- c(max(min(X) + dX, min(Y) + dY), min(max(X) + dX, max(Y) + dY))
extreme <- c(dd[1] + 0.01 * (dd[2] - dd[1]), dd[2] - 0.01 * (dd[2] - dd[1]))

tol <- 2 * globals\$bw.Y
qY <- quantile(Y, t)
c(max(extreme[1],  qY - tol), min(extreme[2], qY + tol))
}

use.smoothing.pp <- TRUE

### Q-Q plots
w1.qq <- function(globals, X, theta, delta, t)
globals\$krnl((delta - X) / globals\$bw.X) - theta

w2.qq <- function(globals, X, theta, delta, t)
globals\$krnl((t - X) / globals\$bw.Y) - theta

alpha1.qq <- function(globals, X, theta, delta, t)
-1

alpha2.qq <- function(globals, X, theta, delta, t)
-1

trange.qq <- function(globals, X, Y)
c(min(Y), max(Y))

deltarange.qq <- function(globals, X, Y, t)
{
c(min(X) + globals\$precision, max(X) - globals\$precision)
}

deltarange.opt.qq <- function(globals, X, Y, t)
{
extreme <- deltarange.qq(globals, X, Y, t)

tol <- globals\$bw.Y
tr <- ecdf(Y)(c(t - tol, t + tol))
upper <- quantile(X, tr[2]) +  globals\$bw.X
lower <- quantile(X, tr[1]) -  globals\$bw.X
c(max(lower, extreme[1]), min(upper, extreme[2]))
}

thetarange.qq <- function(globals, X, Y, delta, t)
{
ww <- w2.qq(globals, Y, 0, 0, t)
tr <- c(min(ww), max(ww))
ww <- w1.qq(globals, X, 0, delta, 0)
dr <- c(min(ww), max(ww))
c(max(tr[1], dr[1]) + globals\$precision, min(tr[2], dr[2]) - globals\$precision)
}

thetarange.opt.qq <- thetarange.qq

use.smoothing.qq <- TRUE

### F difference
w1.fdiff <- function(globals, X, theta, delta, t)
globals\$krnl((t - X) / globals\$bw.X) - theta

w2.fdiff <- function(globals, X, theta, delta, t)
globals\$krnl((t - X) / globals\$bw.Y) - theta - delta

alpha1.fdiff <- function(globals, X, theta, delta, t)
-1

alpha2.fdiff <- function(globals, X, theta, delta, t)
-1

trange.fdiff <- function(globals, X, Y)
c(max(c(min(X), min(Y))), min(c(max(X), max(Y))))

deltarange.fdiff.prim <- function(globals, X, Y, t, prec)
{
ww <- w1.fdiff(globals, X, 0, 0, t)
tr <- c(min(ww), max(ww))
ww <- w2.fdiff(globals, Y, 0, 0, t)
c(min(ww) - tr[2] + prec, max(ww) - tr[1] - prec)
}

deltarange.fdiff <- function(globals, X, Y, t)
deltarange.fdiff.prim(globals, X, Y, t, globals\$precision)

deltarange.opt.fdiff <- function(globals, X, Y, t)
deltarange.fdiff.prim(globals, X, Y, t, 1e-7)

thetarange.fdiff <- function(globals, X, Y, delta, t)
{
ww <- w1.fdiff(globals, X, 0, 0, t)
tr <- c(min(ww), max(ww))
ww <- w2.fdiff(globals, Y, 0, delta, t)
tr2 <- c(min(ww), max(ww))
c(max(tr2[1], tr[1]) + 1e-7, min(tr2[2], tr[2]) - 1e-7)
}

thetarange.opt.fdiff <- thetarange.fdiff

use.smoothing.fdiff <- TRUE

### Quantile difference

w1.qdiff <- function(globals, X, theta, delta, t)
globals\$krnl((theta - X) / globals\$bw.X) - t

w2.qdiff <- function(globals, X, theta, delta, t)
globals\$krnl((theta + delta - X) / globals\$bw.Y) - t

alpha1.qdiff <- function(globals, X, theta, delta, t)
globals\$dkrnl((theta - X) / globals\$bw.X) / globals\$bw.X

alpha2.qdiff <- function(globals, X, theta, delta, t)
globals\$dkrnl((theta + delta - X) / globals\$bw.Y) / globals\$bw.Y

trange.qdiff <- function(globals, X, Y)
{
minlen <- min(length(X), length(Y))
c(1/minlen, 1 - 1/minlen)
}

deltarange.qdiff <- function(globals, X, Y, t)
{
c(min(Y) - max(X) + globals\$precision, max(Y) - min(X) - globals\$precision)
}

deltarange.opt.qdiff <- function(globals, X, Y, t)
{
extreme <- deltarange.qdiff(globals, X, Y, t)

tol <- 2 * globals\$bw.Y
qdiff <- quantile(Y, t) - quantile(X, t)
c(max(extreme[1], qdiff - tol), min(extreme[2], qdiff + tol))
}

thetarange.qdiff <- function(globals, X, Y, delta, t)
{
dX <- globals\$bw.X * globals\$ikrnl(1-t)
dY <- globals\$bw.Y * globals\$ikrnl(1-t)
c(max(min(X) - dX, min(Y) - delta - dY) + globals\$precision, min(max(X) - dX, max(Y) - delta - dY) - globals\$precision)
}

thetarange.opt.qdiff <- function(globals, X, Y, delta, t)
{
extreme <- thetarange.qdiff(globals, X, Y, delta, t)

tol <- 2 * globals\$bw.X
qX <- quantile(X, t)
c(max(extreme[1], qX - tol), min(extreme[2], qX + tol))
}

use.smoothing.qdiff <- TRUE

#### Generalized calling

EL.local <- function(fun, ..., bw = bw.nrd0, method)
{
globals <- list()
globals <- set.method(globals, method)
globals <- set.globals(globals)
globals <- set.kernel(globals)
globals <- set.bw(globals, bw)
fun(globals, ...)
}

#### Calculations

ELconf <- function(globals, X, Y, t, p.level=0.95, delta=NULL)
{
globals <- globals\$initbw(globals, X, Y)
if(is.null(delta))
delta <- deltasolve(globals, X, Y, t)

if (is.na(delta))
return(c(NA, NA))

critval <- qchisq(p.level, 1)
ELlim <- function(delt)
{
ELLR(globals, X, Y, delt, t) - critval
}

drange <- globals\$deltarange(globals, X, Y, t)
if (drange[2] < drange[1])
return(NA)

if (ELlim(delta) < 0)
{
if (ELlim(drange[1]) > 0)
{
lo <- uniroot(ELlim, c(drange[1], delta))\$root
}
else
{
warning(paste("Could not find lower confidence bound for t =", t, "."), call. = FALSE)
lo <- drange[1]
}

if (ELlim(drange[2]) > 0)
{
hi <- uniroot(ELlim, c(delta, drange[2]))\$root
}
else
{
warning(paste("Could not find upper confidence bound for t =", t, "."), call. = FALSE)
hi <- drange[2]
}
}
else
{
warning(paste("Could not find estimate at t =", t, "."), call. = FALSE)
hi <- delta
lo <- delta
}
c(lo, hi)
}

deltasolve <- function(globals, X, Y, t)
{
globals <- globals\$initbw(globals, X, Y)
dr <- globals\$deltarange.opt(globals, X, Y, t)
if (dr[2] < dr[1])
return(NA)
if (dr[2] - dr[1] < 1e-7)
return(dr[1])
globals\$tr <- globals\$thetarange
globals\$thetarange <- globals\$thetarange.opt
oo <- suppressWarnings(optimize(function(d) ELLR(globals, X, Y, d, t), dr))
globals\$thetarange <- globals\$tr
if (is.finite(oo\$objective))
oo\$minimum
else
NA
}

ELLR <- function(globals, X, Y, delta, t)
{
globals <- globals\$initbw(globals, X, Y)
ts <- thetasolve(globals, X, Y, delta, t)
theta <- ts\$theta
globals <- ts\$globals
if (!globals\$degenerate)
{
lr <- 2* (sum(log(1 + globals\$lambda1 * globals\$wval1)) +
sum(log(1 + globals\$lambda2 * globals\$wval2)))
if (lr < 1e-7)
0
else
lr
}
else
Inf
}

thetasolve <- function(globals, X, Y, delta, t)
{
globals <- globals\$initbw(globals, X, Y)
len1 <- length(X)
len2 <- length(Y)
{
globals\$degenerate <<- FALSE
globals\$wval1 <<- globals\$w1(globals, X, theta, delta, t)
globals\$wval2 <<- globals\$w2(globals, Y, theta, delta, t)
alphaval1 <- globals\$alpha1(globals, X, theta, delta, t)
if (is.na(alphaval1[2]))
alphaval1 <- rep.int(alphaval1, len1)
alphaval2 <- globals\$alpha2(globals, Y, theta, delta, t)
if (is.na(alphaval2[2]))
alphaval2 <- rep.int(alphaval2, len2)
res <- .C("theta_equation",
as.integer(len1),
as.double(globals\$wval1),
as.double(alphaval1),
as.integer(len2),
as.double(globals\$wval2),
as.double(alphaval2),
l1 = double(1),
l2 = double(1),
res = double(1))
globals\$lambda1 <<- res\$l1
globals\$lambda2 <<- res\$l2
res\$res
}

tr <- globals\$thetarange(globals, X, Y, delta, t)
if (tr[1] > tr[2])
{
globals\$degenerate <- TRUE
return(list(theta = NA, globals=globals))
}
theta <- suppressWarnings(try(uniroot(tgrad, tr)\$root, silent=TRUE))
list(theta = theta, globals = globals)
}
```

## Try the EL package in your browser

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

EL documentation built on Dec. 28, 2022, 2:47 a.m.