R/condbandwidth.R

In np: Nonparametric Kernel Smoothing Methods for Mixed Data Types

condbandwidth <-
  function(xbw,
           ybw,
           bwmethod = c("cv.ls","normal-reference", "manual"),
           bwscaling = FALSE,
           bwtype = c("fixed","generalized_nn","adaptive_nn"),
           cxkertype = c("gaussian","truncated gaussian","epanechnikov","uniform"), 
           cxkerorder = c(2,4,6,8),
           cxkerbound = c("none","range","fixed"),
           cxkerlb = NULL,
           cxkerub = NULL,
           uxkertype = c("aitchisonaitken","liracine"),
           oxkertype = c("liracine","wangvanryzin","racineliyan"),
           cykertype = c("gaussian","truncated gaussian","epanechnikov","uniform"), 
           cykerorder = c(2,4,6,8),
           cykerbound = c("none","range","fixed"),
           cykerlb = NULL,
           cykerub = NULL,
           uykertype = c("aitchisonaitken","liracine"),
           oykertype = c("liracine","wangvanryzin","racineliyan"),
           fval = NA,
           ifval = NA,
           num.feval = NA,
           num.feval.fast = NA,
           nn.cache = NULL,
           fval.history = NA,
           eval.history = NA,
           invalid.history = NA,
           nobs = NA,
           xdati, ydati,
           xnames = character(length(xbw)),
           ynames = character(length(ybw)),
           sfactor = NA, bandwidth = NA,
           rows.omit = NA,
           nconfac = NA,
           ncatfac = NA,
           sdev = NA,
           bandwidth.compute = TRUE,
           timing = NA,
           total.time = NA,
           regtype = "lc",
           pregtype = "Local-Constant",
           basis = "glp",
           degree = integer(0),
           bernstein.basis = FALSE,
           regtype.engine = "lc",
           basis.engine = "glp",
           degree.engine = integer(0),
           bernstein.basis.engine = FALSE,
           ...){

  if (missing(xbw) || missing(ybw))
    stop("improper invocation of condbandwidth constructor: 'bw' or i[cuo]* missing")
  
  xndim = length(xbw)
  yndim = length(ybw)
  
  bwmethod = match.arg(bwmethod)
  bwtype = match.arg(bwtype)

  cxkertype = match.arg(cxkertype)
  cykertype = match.arg(cykertype)
  cxkerbound = match.arg(cxkerbound)
  cykerbound = match.arg(cykerbound)

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

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

  if (bwmethod == "normal-reference" && (cxkertype != "gaussian" || bwtype != "fixed")){    
    .np_warning("normal-reference bandwidth selection assumes gaussian kernel with fixed bandwidth")
    bwtype = "fixed"
    cxkertype = "gaussian"
  }

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

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

  if (bwmethod == "normal-reference" && (cykertype != "gaussian" || bwtype != "fixed")){    
    .np_warning("normal-reference bandwidth selection assumes gaussian kernel with fixed bandwidth")
    bwtype = "fixed"
    cykertype = "gaussian"
  }

  if (cxkerorder != cykerorder && bwscaling)
    stop("scale factors with different order kernels for dependent and explanatory variables is unsupported")
  
  uxkertype = match.arg(uxkertype)
  uykertype = match.arg(uykertype)
  
  oxkertype = match.arg(oxkertype)
  oykertype = match.arg(oykertype)
  cxbounds <- npKernelBoundsResolve(
    dati = xdati,
    varnames = xnames,
    kerbound = cxkerbound,
    kerlb = cxkerlb,
    kerub = cxkerub,
    argprefix = "cxker")
  cybounds <- npKernelBoundsResolve(
    dati = ydati,
    varnames = ynames,
    kerbound = cykerbound,
    kerlb = cykerlb,
    kerub = cykerub,
    argprefix = "cyker")
  bounded_nonfixed_supported <- bwtype %in% c("generalized_nn", "adaptive_nn")
  if (bwtype != "fixed" &&
      (cxbounds$bound != "none" || cybounds$bound != "none") &&
      !bounded_nonfixed_supported)
    stop("finite continuous kernel bounds require bwtype = \"fixed\"")

  pxorder = switch( cxkerorder/2, "Second-Order", "Fourth-Order", "Sixth-Order", "Eighth-Order" )
  pyorder = switch( cykerorder/2, "Second-Order", "Fourth-Order", "Sixth-Order", "Eighth-Order" )

  dati <- list(x = xdati, y = ydati)
  
  if (!identical(sfactor,NA)){
    ## using the new model for generically accessing bandwidth objects

    okertype <- list(x = oxkertype, y = oykertype)
    ukertype <- list(x = uxkertype, y = uykertype)

    scaleOrMax <- function(i, j) {
      if (dati[[j]]$icon[i])
        return((sfactor[[j]])[i])

      if (dati[[j]]$iord[i])
        return(oMaxL(dati[[j]]$all.nlev[[i]], kertype = okertype[[j]]))
      
      if (dati[[j]]$iuno[i])
        return(uMaxL(dati[[j]]$all.nlev[[i]], kertype = ukertype[[j]]))
    }

    sumNum <- list(x = NA, y = NA)
    sumNum[] <- lapply(seq_along(dati), function(i) {
      sapply(seq_along(dati[[i]]$icon), scaleOrMax, j = i)
    })
  } else {
    sumNum <- NA
  }

  if (length(rows.omit) == 0)
    rows.omit <- NA

  mybw = list(
    xbw=xbw,
    ybw=ybw,
    method = bwmethod,
    pmethod = bwmToPrint(bwmethod),
    fval = fval,
    ifval = ifval,
    num.feval = num.feval,
    num.feval.fast = num.feval.fast,
    nn.cache = nn.cache,
    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),
    cxkertype = cxkertype,
    cykertype = cykertype,
    cxkerorder = cxkerorder,
    cykerorder = cykerorder,
    cxkerbound = cxbounds$bound,
    cxkerlb = cxbounds$lb,
    cxkerub = cxbounds$ub,
    cykerbound = cybounds$bound,
    cykerlb = cybounds$lb,
    cykerub = cybounds$ub,
    pcxkertype = cktToPrint(cxkertype, order = pxorder, kerbound = cxbounds$bound),
    pcykertype = cktToPrint(cykertype, order = pyorder, kerbound = cybounds$bound),
    uxkertype = uxkertype,
    uykertype = uykertype,
    puxkertype = uktToPrint(uxkertype),
    puykertype = uktToPrint(uykertype),
    oxkertype = oxkertype,
    oykertype = oykertype,
    poxkertype = oktToPrint(oxkertype),
    poykertype = oktToPrint(oykertype, normalized = TRUE),
    nobs = nobs,
    xndim = xndim,
    yndim = yndim,
    ndim = xndim + yndim,
    xncon = sum(xdati$icon),
    xnuno = sum(xdati$iuno),
    xnord = sum(xdati$iord),
    yncon = sum(ydati$icon),
    ynuno = sum(ydati$iuno),
    ynord = sum(ydati$iord),
    ncon = sum(c(xdati$icon, ydati$icon)),
    ixcon = xdati$icon,
    ixuno = xdati$iuno,
    ixord = xdati$iord,
    iycon = ydati$icon,
    iyuno = ydati$iuno,
    iyord = ydati$iord,
    xnames = xnames,
    ynames = ynames,
    xdati = xdati,
    ydati = ydati,
    xmcv = mcvConstruct(xdati),
    ymcv = mcvConstruct(ydati),
    sfactor = sfactor,
    bandwidth = bandwidth,
    nconfac = nconfac,
    ncatfac = ncatfac,
    sdev = sdev,
    sumNum = sumNum,
    dati = dati, 
    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,
    regtype = regtype,
    pregtype = pregtype,
    basis = basis,
    degree = degree,
    bernstein.basis = bernstein.basis,
    regtype.engine = regtype.engine,
    basis.engine = basis.engine,
    degree.engine = degree.engine,
    bernstein.basis.engine = bernstein.basis.engine)

  mybw$klist = list(
    x =
    list(ckertype = cxkertype,
         ckerbound = cxbounds$bound,
         ckerlb = cxbounds$lb,
         ckerub = cxbounds$ub,
         pckertype = mybw$pcxkertype,
         ukertype = uxkertype,
         pukertype = mybw$puxkertype,
         okertype = oxkertype,
         pokertype = mybw$poxkertype),
    y =
    list(ckertype = cykertype,
         ckerbound = cybounds$bound,
         ckerlb = cybounds$lb,
         ckerub = cybounds$ub,
         pckertype = mybw$pcykertype,
         ukertype = uykertype,
         pukertype = mybw$puykertype,
         okertype = oykertype,
         pokertype = mybw$poykertype))

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


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

print.condbandwidth <- function(x, digits=NULL, ...){
  cat("\nConditional distribution data (",x$nobs," observations, ",
      (x$xndim+x$yndim)," variable(s))",
      "\n(", x$yndim, " dependent variable(s), and ", x$xndim, " explanatory variable(s))\n\n",
      sep="")
  print(matrix(x$ybw,ncol=x$yndim,dimnames=list(paste("Dep. Var. ",x$pscaling,":",sep=""),x$ynames)))

  print(matrix(x$xbw,ncol=x$xndim,dimnames=list(paste("Exp. Var. ",x$pscaling,":",sep=""),x$xnames)))

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


summary.condbandwidth <- function(object, ...) {
  cat("\nConditional distribution data (",object$nobs," observations, ",
      (object$xndim+object$yndim)," variable(s))",
      "\n(", object$yndim, " dependent variable(s), and ", object$xndim, " explanatory variable(s))\n",
      sep="")

  cat(genOmitStr(object))
  cat(genBwSelStr(object))

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

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

predict.condbandwidth <- function(object, ...) { npcdist(bws = object, ...) }