R/rbandwidth.R

In np: Nonparametric Kernel Smoothing Methods for Mixed Data Types

rbandwidth <-
  function(bw = stop("rbandwidth:argument 'bw' missing"),
           regtype = c("lc","ll","lp"),
           basis = c("glp","additive","tensor"),
           degree = NULL,
           bernstein.basis = FALSE,
           bwmethod = c("cv.ls","cv.aic"),
           bwscaling = FALSE,
           bwtype = c("fixed","generalized_nn","adaptive_nn"),
           ckertype = c("gaussian","truncated gaussian","epanechnikov","uniform"),
           ckerorder = c(2,4,6,8),
           ckerbound = c("none","range","fixed"),
           ckerlb = NULL,
           ckerub = NULL,
           ukertype = c("aitchisonaitken", "liracine"),
           okertype = c("liracine","wangvanryzin","racineliyan"),
           fval = NA,
           ifval = NA,
           num.feval = NA,
           num.feval.fast = NA,
           fval.history = NA,
           eval.history = NA,
           invalid.history = NA,
           nobs = NA,
           xdati = stop("rbandwidth:argument 'xdati' missing"),
           ydati = stop("rbandwidth:argument 'ydati' missing"),
           xnames = character(length(bw)),
           ynames = character(1),
           sfactor = NA, bandwidth = NA,
           rows.omit = NA,
           nconfac = NA,
           ncatfac = NA,
           sdev = NA,
           bandwidth.compute = TRUE,
           timing = NA,
           total.time = NA,
           ...){

  ndim = length(bw)
  npRejectLegacyLpArgs(names(list(...)), where = "rbandwidth")
  regtype = match.arg(regtype)
  basis <- npValidateLpBasis(regtype = regtype, basis = basis)
  bwmethod = match.arg(bwmethod)
  bwtype = match.arg(bwtype)
  ckertype = match.arg(ckertype)
  ckerbound = match.arg(ckerbound)

  if(missing(ckerorder))
    ckerorder = 2
  else if (ckertype == "uniform")
    .np_warning("ignoring kernel order specified with uniform kernel type")
  else {
    kord = c(2,4,6,8) 
    if (!any(kord == ckerorder))
      stop("ckerorder must be one of ", paste(kord,collapse=" "))
  }

  if (ckertype == "truncated gaussian" && ckerorder != 2)
    .np_warning("using truncated gaussian of order 2, higher orders not yet implemented")

  ukertype = match.arg(ukertype)
  okertype = match.arg(okertype)
  cbounds <- npKernelBoundsResolve(
    dati = xdati,
    varnames = xnames,
    kerbound = ckerbound,
    kerlb = ckerlb,
    kerub = ckerub,
    argprefix = "cker")
  bounded_nonfixed_supported <- bwtype %in% c("generalized_nn", "adaptive_nn")
  if (bwtype != "fixed" &&
      cbounds$bound != "none" &&
      !bounded_nonfixed_supported)
    stop("finite continuous kernel bounds require bwtype = \"fixed\"")

  ncon <- sum(xdati$icon)
  degree <- npValidateGlpDegree(regtype = regtype,
                                    degree = degree,
                                    ncon = ncon)
  bernstein.basis <- npValidateGlpBernstein(regtype = regtype,
                                          bernstein.basis = bernstein.basis)
  reg.spec <- npCanonicalConditionalRegSpec(
    regtype = regtype,
    basis = basis,
    degree = degree,
    bernstein.basis = bernstein.basis,
    ncon = ncon,
    where = "rbandwidth"
  )
  if (identical(regtype, "lp") && ncon > 0L && is.finite(nobs)) {
    lp.dim <- dim_basis(basis = basis,
                        kernel = TRUE,
                        degree = degree,
                        segments = rep.int(1L, ncon))
    if (is.finite(lp.dim) && lp.dim > (nobs - 1.0))
      stop(sprintf("LP basis dimension (%s) exceeds nobs - 1 (%s); reduce degree",
                   format(lp.dim, trim = TRUE, scientific = FALSE),
                   format(nobs - 1.0, trim = TRUE, scientific = FALSE)))
  }

  porder = switch( ckerorder/2, "Second-Order", "Fourth-Order", "Sixth-Order", "Eighth-Order" )
  ## calculate some info to be pretty-printed

  if (!identical(sfactor,NA)){
    sumNum <- sapply(seq_len(ndim), function(i) {
      if (xdati$icon[i])
        return(sfactor[i])

      if (xdati$iord[i])
        return(oMaxL(xdati$all.nlev[[i]], kertype = okertype))
      
      if (xdati$iuno[i])
        return(uMaxL(xdati$all.nlev[[i]], kertype = ukertype))
    })
  } else {
    sumNum <- NA
  }

  if (length(rows.omit) == 0)
    rows.omit <- NA
  
  mybw = list(
    bw=bw,
    regtype = regtype,
    pregtype = switch(regtype,
      lc = "Local-Constant",
      ll = "Local-Linear",
      lp = "Local-Polynomial"),
    basis = basis,
    degree = degree,
    bernstein.basis = bernstein.basis,
    regtype.engine = reg.spec$regtype.engine,
    basis.engine = reg.spec$basis.engine,
    degree.engine = reg.spec$degree.engine,
    bernstein.basis.engine = reg.spec$bernstein.basis.engine,
    method = bwmethod,
    pmethod = bwmToPrint(bwmethod),
    fval = fval,
    ifval = ifval,
    num.feval = num.feval,
    num.feval.fast = num.feval.fast,
    fval.history = fval.history,
    eval.history = eval.history,
    invalid.history = invalid.history,
    scaling = bwscaling,
    pscaling = npBandwidthSummaryLabel(bwtype = bwtype, bwscaling = bwscaling),
    type = bwtype,
    ptype = bwtToPrint(bwtype),
    ckertype = ckertype,    
    ckerorder = ckerorder,
    ckerbound = cbounds$bound,
    ckerlb = cbounds$lb,
    ckerub = cbounds$ub,
    pckertype = cktToPrint(ckertype, order = porder, kerbound = cbounds$bound),
    ukertype = ukertype,
    pukertype = uktToPrint(ukertype),
    okertype = okertype,
    pokertype = oktToPrint(okertype),
    nobs = nobs,
    ndim = ndim,
    ncon = ncon,
    nuno = sum(xdati$iuno),
    nord = sum(xdati$iord),
    icon = xdati$icon,
    iuno = xdati$iuno,
    iord = xdati$iord,
    xnames = xnames,
    ynames = ynames,
    xdati = xdati,
    ydati = ydati,
    xmcv = mcvConstruct(xdati),
    sfactor = list(x = sfactor),
    bandwidth = list(x = bandwidth),
    nconfac = nconfac,
    ncatfac = ncatfac,
    sdev = sdev,
    sumNum = list(x = sumNum),
    dati = list(x = xdati, y = ydati),
    varnames = list(x = xnames, y = ynames),
    vartitle = list(x = "Explanatory", y = "Dependent"),
    vartitleabb = list(x = "Exp.", y = "Dep."),
    rows.omit = rows.omit,
    nobs.omit = if (identical(rows.omit, NA)) 0 else length(rows.omit),
    timing = timing,
    total.time = total.time)

  
  mybw$klist = list(
    x =
    list(ckertype = ckertype,
         pckertype = mybw$pckertype,
         ukertype = ukertype,
         pukertype = mybw$pukertype,
         okertype = okertype,
         pokertype = mybw$pokertype))

  if(!bandwidth.compute)
    mybw$pmethod <- "Manual"

  class(mybw) = "rbandwidth"
  if(!any(is.na(mybw$bandwidth)))
    validateBandwidth(mybw)
  mybw
  
}

as.double.rbandwidth <- function(x, ...){
  x$bw
}

npInsertGlpSummary <- function(txt, basis, degree, bernstein){
  lines <- strsplit(txt, "\n", fixed = TRUE)[[1]]
  reg.idx <- grep("^Regression Type:", lines)

  if (length(reg.idx) != 1L)
    return(txt)
  if (is.null(basis) || !length(basis))
    basis <- "glp"

  basis.family <- npLpBasisFamilyLabel(basis)
  basis.rep <- npLpBasisRepresentationLabel(bernstein)
  basis.dim <- npLpBasisNcol(basis = basis, degree = degree)

  glp.lines <- c(
    paste("LP Basis Family:", basis.family),
    paste("LP Basis Representation:", basis.rep),
    paste("LP Basis Dimension (NCOL(B)):",
          format(basis.dim, scientific = FALSE, trim = TRUE))
  )

  lines <- append(lines, glp.lines, after = reg.idx)
  paste(lines, collapse = "\n")
}

## feature: when using dataframe interface, summary and print methods don't 
## provide info on the dependent variable
print.rbandwidth <- function(x, digits=NULL, ...){
  cat("\nRegression Data (",x$nobs," observations, ",x$ndim," variable(s)):\n\n",sep="")
  print(matrix(x$bw,ncol=x$ndim,dimnames=list(paste(x$pscaling,":",sep=""),x$xnames)))

  bw.sel.str <- genBwSelStr(x)
  if (identical(x$regtype, "lp") && x$ncon > 0)
    bw.sel.str <- npInsertGlpSummary(txt = bw.sel.str,
                                     basis = x$basis,
                                     degree = x$degree,
                                     bernstein = x$bernstein.basis)
  cat(bw.sel.str)
  cat(genBwKerStrs(x))

  cat("\n\n")
  if(!missing(...))
    print(...,digits=digits)
  invisible(x)
}


predict.rbandwidth <- function(object, ...) { npreg(bws = object, ...) }

summary.rbandwidth <- function(object, ...){
  cat("\nRegression Data (",object$nobs," observations, ",object$ndim," variable(s)):\n",sep="")

  cat(genOmitStr(object))
  bw.sel.str <- genBwSelStr(object)
  if (identical(object$regtype, "lp") && object$ncon > 0)
    bw.sel.str <- npInsertGlpSummary(txt = bw.sel.str,
                                     basis = object$basis,
                                     degree = object$degree,
                                     bernstein = object$bernstein.basis)
  cat(bw.sel.str)

  cat('\n')
  cat(genBwScaleStrs(object))
  cat(genBwKerStrs(object))

  cat(genTimingStr(object))
  cat("\n\n")

}