# Example:  Failure times of integrated circuits under accelerated testing.  
# 4156 circuits were tested for 1370 hours at elevated temperature and 
# humidity.  The goal was to estimate the proportion of circuits with a 
# manufacturing defect that leads to relatively early failure.  28 of the 
# 4156 circuits failed, at various times, during the test; the failure time 
# of the remaining circuits was censored at 1370 hours.  

# The data and proposed statistical model are given in section 18.1 of 
# Meeker, Hahn, and Escobar (2017).  The model assumes that a fraction 
# 'dfrac' of circuits are defective, with failure times having a 
# Weibull distribution with parameters 'shape' and 'scale'.  The remaining 
# fraction '1 - dfrac' are assumed to have essentially infinite failure 
# times.  Meeker, Hahn, and Escobar refer to this as the limited failure 
# population (LFP) model.  It is also a version of a "cure model", where 
# '1 - dfrac' is the proportion of patients who are cured and thus never 
# succumb to the disease being studied. 

# Here we examine the 3-dimensional likelihood ratio-based 95% confidence 
# region for the full set of model parameters (dfrac, shape, scale). 
# Therefore the response function is the log-likelihood, and the 
# confidence region is the level set consisting of parameter vectors whose 
# log-likelihood is no more than 'qchisq(0.95, df=3) / 2 = 3.907' below its 
# maximum.

# Meeker, William Q., Hahn, Gerald J., Escobar, Luis A. Statistical Intervals: 
# A Guide for Practitioners and Researchers, 2nd ed.  John Wiley & Sons, 2017.

# This script creates and saves (as 'cicuitFailure_3dEx.rda' in directory 
# '../data') a list with components:
#  fobj : 'fnObj' object with the response function, response and input  
#      space specifications, and feasible region function.
#  thresh : Log-likelihood value defining the 95% asymptotic confidence 
#      region.
#  theta_mle : Three-element vector containing the maximum likelihood estimate 
#      of the parameters.
#  theta_vcov : 3-by-3 matrix containing the approximate variance-covariance 
#      matrix for 'theta_mle'.
#  bsrch : List returned by 'bdrySearch()' with results from a search for 
#      level set boundary points and segments.
#  slsrch : List returned by 'slicedBdrySearch()', with results from a 
#      level set search within six slices of input space, each defined by 
#      a specified value of dimension 'scale'.

# The 'theta_mle' and 'theta_vcov' components were obtained by 'stats::optim'.

# 11 Feb 2026 / Richard Raubertas
# 11 Apr 2026 : Update to be closer to the version used in the confidence 
#               region vignette.  (No empty slices.)

local({
  #----- Specify the data and find the MLE of the model parameters.
  
  # The data from Meeker, Hahn, and Escobar (2017):
  ftime <- c(0.10, 0.15, 0.60, 0.80, 1.20, 2.50, 3.00, 4.00, 
             6.00, 10.00, 12.50, 20, 43, 48, 54, 74, 84, 94, 
             168, 263, 593, 1370)
  count <- c(   2,    1,    1,    2,    1,    1,    1,    2, 
                1,     2,     1,  2,  2,  2,  1,  1,  1,  1, 
                1,   1,   1, 4128)
  fail <- rep(c(TRUE, FALSE), c(21, 1))

  # Log-likelihood function, allowing for censoring, in the form required by 
  # 'stats::optim'.
  llik1 <- function(theta) {
    # Defined by lexical scoping:  ftime, count, fail
    nfail <- sum(count[fail])
    dfrac <- theta[1]
    shape <- theta[2]
    scale <- theta[3]
    llik_fail <- nfail*log(dfrac) +
                 sum(count[fail] * dweibull(ftime[fail], shape=shape, 
                                            scale=scale, log=TRUE))
    llik_cens <- sum(count[!fail] * log(1 - dfrac*pweibull(ftime[!fail], 
                                                           shape=shape, 
                                                           scale=scale, 
                                                           log.p=FALSE, 
                                                           lower.tail=TRUE)))
    rslt <- llik_fail + llik_cens
    rslt
  }

  # NOTE:  parameter scaling ('parscale') helps a lot.
  mle_fit <- stats::optim(par=c(dfrac=sum(count[fail])/sum(count), 
                                shape=1/2, 
                                scale=exp(3)),  # start values close to ref
                          fn=llik1, hessian=TRUE, 
                          method="L-BFGS-B", lower=c(1e-4, 1e-4, 1e-4), 
                          control=list(fnscale=-1, trace=0, REPORT=4, 
                                       parscale=c(0.0067, 0.5, 28)))
  theta_mle <- mle_fit$par  # c(dfrac, shape, scale) = 
                            # c(0.006744423, 0.495980540, 28.366738156)
  theta_vcov <- solve(-1*mle_fit$hessian)
  max_llik <- mle_fit$value  # -293.0329 = llik1(theta_mle)

  #----- Set up the level set problem.

  # Log-likelihood function for the data, in a form suitable for the 
  # 'levelSets' package (handles multiple parameter vectors, as rows of a 
  # matrix).
  llik_lfp <- function(x) {
    # 'x' : matrix, each row contains coordinates of a point in 
    # input space (i.e., a parameter vector for the model).
    # Return a vector with one log-likelihood value per parameter vector.
    # Defined by lexical scoping:  ftime, count, fail
    ntheta <- nrow(x)
    rslt <- rep(NA_real_, ntheta)
    stopifnot(colnames(x) == c("dfrac", "shape", "scale"))
    nfail <- sum(count[fail])
    for (i in seq_len(ntheta)) {
      theta <- x[i, ]
      dfrac <- theta[1]
      shape <- theta[2]
      scale <- theta[3]
      llik_fail <- nfail*log(dfrac) +
                   sum(count[fail] * dweibull(ftime[fail], shape=shape, 
                                              scale=scale, log=TRUE))
      llik_cens <- sum(count[!fail] * log(1 - dfrac*pweibull(ftime[!fail], 
                                                             shape=shape, 
                                                             scale=scale, 
                                                             log.p=FALSE, 
                                                             lower.tail=TRUE)))
      rslt[i] <- llik_fail + llik_cens
    }
    rslt
  }
  
  fobj <- fnObj(names(theta_mle), respfn=llik_lfp, 
                # Require all parameters to be strictly positive:
                feasbnds=rbind(rep(1e-4, 3), rep(NA, 3)))
  thresh <- max_llik - (qchisq(0.95, df=length(theta_mle))/2)
    # (threshold for the log-likelihood that defines asymptotic 95% confidence 
    # region for the parameter vector)
  # Since the confidence region is defined in terms of simple differences 
  # between loglikelihoods, independent of their magnitudes, set 'resptol' 
  # to use only absolute tolerance, not relative tolerance:
  fobj <- update(fobj, resptol=c(0, 1e-3))

  # A natural choice for a ray origin is the MLE.  The initial search region is 
  # set to be +/-3 standard errors around the MLE.  Input space dimensions 
  # are scaled by the asymptotic variance-covariance matrix from the MLE fit:
  rect0 <- inpRect(rbind(theta_mle - 3*sqrt(diag(theta_vcov)), 
                         theta_mle + 3*sqrt(diag(theta_vcov))), 
                   spec=fobj)
  scl0 <- inpScale(theta_vcov, spec=fobj)

  # Function `bdrySearch()` searches for a level set boundary in a series 
  # of steps, where each step adaptively selects a new ray origin in an 
  # attempt to explore new portions of the boundary.  The search region and 
  # input dimension scaling can also be updated at each step. (See 
  # `?srchControl` for the full set of control options.)
  bsrch <- bdrySearch(fobj, lsetthresh=thresh, rect=rect0, scale=scl0, 
                      initOrigin=theta_mle, 
                      control=list(axisDegree=2, updateScale=2))  # list

  # Can also search for the level set boundary within slices of input space. 
  # Consider slices defined by six values of input dimension 'scale'.
  slices <- sliceMat(list(scale=c(10, 30, 50, 70, 90, 110)), spec=fobj)  
  
  # To apply `slicedBdrySearch()` we need to specify a search region that 
  # intersects all the slices, such as 
  rect1 <- boundingRect(bsrch$rect, slices, expand=1.1)
  
  slsrch <- slicedBdrySearch(fobj, lsetthresh=thresh, slices=slices, 
                             rect=rect1, scale=bsrch$scale, 
                             currentBdry=bsrch, 
                             control=list(axisDegree=2, 
                                          maxSteps=6))  # list with one component per slice
  
  circuitFailure_3dEx <- list(fobj=fobj, thresh=thresh, 
                              theta_mle=theta_mle, theta_vcov=theta_vcov, 
                              bsrch=bsrch, slsrch=slsrch)
  objs <- c("circuitFailure_3dEx")
  for (obj in objs) {
    save(list=obj, file=paste0("../../data/", obj, ".rda"), compress=TRUE)
  }

  # Clean up.  (This may not be necessary due to 'local()' wrap.)
  rm(ftime, count, fail, llik1, mle_fit, theta_mle, theta_vcov, max_llik, 
     llik_lfp, fobj, thresh, objs, obj, circuitFailure_3dEx, rect0, scl0, 
     bsrch, slices, rect1, slsrch)
})

if (FALSE) {
  # Code like this can be inserted at the start of example sections of 
  # help pages:
  Ex <- circuitFailure_3dEx  # see '?circuitFailure_3dEx'
  fobj <- Ex$fobj
  thresh <- Ex$thresh
  theta_mle <- Ex$theta_mle
  theta_vcov <- Ex$theta_vcov
  bsrch <- Ex$bsrch
  slsrch <- Ex$slsrch
}
