# R/np.plot.R In np: Nonparametric Kernel Smoothing Methods for Mixed Data Types

#### Documented in npplotnpplot.bandwidthnpplot.conbandwidthnpplot.plbandwidthnpplot.rbandwidthnpplot.scbandwidthnpplot.sibandwidthuocquantile

```## the idea is that you have done bandwidth selection
## you just need to supply training data and the bandwidth

gen.label = function(label, altlabel){
paste(ifelse(is.null(label), altlabel, label))
}

gen.tflabel = function(condition, tlabel, flabel){
paste(ifelse(condition,tlabel,flabel))
}

draw.error.bands = function(ex, ely, ehy, lty = 2){
lines(ex,ely,lty=lty)
lines(ex,ehy,lty=lty)
}

draw.error.bars = function(ex, ely, ehy, hbar = TRUE, hbarscale = 0.3, lty = 2){
yy = double(3*length(ex))
jj = 1:length(ex)*3

yy[jj-2] = ely
yy[jj-1] = ehy
yy[jj] = NA

xx = double(3*length(ex))
xx[jj-2] = ex
xx[jj-1] = ex
xx[jj] = NA

lines(xx,yy,lty=lty)

if (hbar){
## hbars look silly if they are too wide in relation to their height
## this only matters in the limit of few points, since that is when
## hbardist may get relatively large

golden = (1+sqrt(5))/2
hbardist = abs(max(ex) - min(ex))/length(ex)*hbarscale

yg = abs(yy[jj-2]-yy[jj-1])/golden
htest = (hbardist >= yg)

xx[jj-2] = ex - ifelse(htest, yg/2, hbardist/2)
xx[jj-1] = ex + ifelse(htest, yg/2, hbardist/2)

ty = yy[jj-1]
yy[jj-1] = yy[jj-2]

lines(xx,yy)

yy[jj-2] = ty
yy[jj-1] = ty

lines(xx,yy)
}
}

draw.errors =
function(ex, ely, ehy,
plot.errors.style,
plot.errors.bar,
plot.errors.bar.num,
lty){
if (plot.errors.style == "bar"){
ei = seq(1,length(ex),length.out = min(length(ex),plot.errors.bar.num))
draw.error.bars(ex = ex[ei],
ely = ely[ei],
ehy = ehy[ei],
hbar = (plot.errors.bar == "I"),
lty = lty)
} else if (plot.errors.style == "band") {
draw.error.bands(ex = ex,
ely = ely,
ehy = ehy,
lty = lty)
}
}

plotFactor <- function(f, y, ...){
plot(x = f, y = y, lty = "blank", ...)

l.f = rep(f, each=3)
l.f[3*(1:length(f))] = NA

l.y = unlist(lapply(y, function (p) { c(0,p,NA) }))

lines(x = l.f, y = l.y, lty = 2)
points(x = f, y = y)
}

compute.bootstrap.errors = function(...,bws){
UseMethod("compute.bootstrap.errors",bws)
}

compute.bootstrap.errors.rbandwidth =
function(xdat, ydat,
exdat,
slice.index,
plot.errors.boot.method,
plot.errors.boot.blocklen,
plot.errors.boot.num,
plot.errors.center,
plot.errors.type,
plot.errors.quantiles,
bws){
boot.err = matrix(data = NA, nrow = dim(exdat)[1], ncol = 3)

is.inid = plot.errors.boot.method=="inid"

strf = ifelse(is.inid, "function(data,indices){", "function(tsb){")
strtx = ifelse(is.inid, "txdat = xdat[indices,],",
"txdat = tsb[,1:(ncol(tsb)-1),drop=FALSE],")
strty = ifelse(is.inid, "tydat = ydat[indices],",
"tydat = tsb[,ncol(tsb)],")

boofun = eval(parse(text=paste(strf, "npreg(", strtx, strty,
"exdat = exdat, bws = bws,",

if (is.inid){
boot.out = boot(data = data.frame(xdat,ydat), statistic = boofun,
R = plot.errors.boot.num)
} else {
boot.out = tsboot(tseries = data.frame(xdat,ydat), statistic = boofun,
R = plot.errors.boot.num,
l = plot.errors.boot.blocklen,
sim = plot.errors.boot.method)
}

all.bp <- list()

if (slice.index > 0 && (bws\$xdati\$iord | bws\$xdati\$iuno)[slice.index]){
boot.frame <- as.data.frame(boot.out\$t)
u.lev <- bws\$xdati\$all.ulev[[slice.index]]

## if we are bootstrapping a factor, there should be one
## set of replications for each level
stopifnot(length(u.lev)==ncol(boot.frame))

all.bp\$stats <- matrix(data = NA, nrow = 5, ncol = length(u.lev))
all.bp\$conf <- matrix(data = NA, nrow = 2, ncol = length(u.lev))

for (i in 1:length(u.lev)){
t.bp <- boxplot.stats(boot.frame[,i])
all.bp\$stats[,i] <- t.bp\$stats
all.bp\$conf[,i] <- t.bp\$conf
all.bp\$out <- c(all.bp\$out,t.bp\$out)
all.bp\$group <- c(all.bp\$group, rep.int(i,length(t.bp\$out)))
}
all.bp\$n <- rep.int(plot.errors.boot.num, length(u.lev))
all.bp\$names <- bws\$xdati\$all.lev[[slice.index]]
rm(boot.frame)
}

if (plot.errors.type == "standard") {
boot.err[,1:2] = 2.0*sqrt(diag(cov(boot.out\$t)))
}
else if (plot.errors.type == "quantiles") {
boot.err[,1:2] = t(sapply(as.data.frame(boot.out\$t),
function (y) {
quantile(y,probs = plot.errors.quantiles)
}))
boot.err[,1] = boot.out\$t0 - boot.err[,1]
boot.err[,2] = boot.err[,2] - boot.out\$t0
}
if (plot.errors.center == "bias-corrected")
boot.err[,3] <- 2*boot.out\$t0-colMeans(boot.out\$t)
list(boot.err = boot.err, bxp = all.bp)
}

compute.bootstrap.errors.scbandwidth =
function(xdat, ydat, zdat,
exdat, ezdat,
slice.index,
plot.errors.boot.method,
plot.errors.boot.blocklen,
plot.errors.boot.num,
plot.errors.center,
plot.errors.type,
plot.errors.quantiles,
bws){
miss.z <- missing(zdat)
boot.err = matrix(data = NA, nrow = dim(exdat)[1], ncol = 3)

is.inid = plot.errors.boot.method=="inid"

xi <- 1:ncol(xdat)
yi <- ncol(xdat)+1
if (!miss.z)
zi <- yi+1:ncol(zdat)

strf = ifelse(is.inid, "function(data,indices){", "function(tsb){")
strtx = ifelse(is.inid, "txdat = xdat[indices,, drop = FALSE],",
"txdat = tsb[,xi,drop=FALSE],")
strty = ifelse(is.inid, "tydat = ydat[indices],",
"tydat = tsb[,yi],")
strtz <-
ifelse(miss.z, '',
ifelse(is.inid, 'tzdat = zdat[indices,, drop = FALSE],',
'tzdat = tsb[,zi, drop = FALSE],'))

boofun = eval(parse(text=paste(strf, "npscoef(", strtx, strty, strtz,
"exdat = exdat,", ifelse(miss.z,"", "ezdat = ezdat,"),
"bws = bws, iterate = FALSE)\$",
"mean", "}", sep="")))

boot.out <-
eval(parse(text = paste(ifelse(is.inid, 'boot(data = ',
'tsboot(tseries ='), 'data.frame(xdat,ydat',
ifelse(miss.z,'', ',zdat'),
'), statistic = boofun, R = plot.errors.boot.num',
ifelse(is.inid,'','l = plot.errors.boot.blocklen, sim = plot.errors.boot.method'),')')))

all.bp <- list()

if ((slice.index > 0) && (((slice.index <= ncol(xdat)) && (bws\$xdati\$iord | bws\$xdati\$iuno)[slice.index]) ||
((slice.index > ncol(xdat)) && (bws\$zdati\$iord | bws\$zdati\$iuno)[slice.index-ncol(xdat)]))) {
boot.frame <- as.data.frame(boot.out\$t)

if(slice.index <= ncol(xdat))
u.lev <- bws\$xdati\$all.ulev[[slice.index]]
else
u.lev <- bws\$zdati\$all.ulev[[slice.index-ncol(xdat)]]

## if we are bootstrapping a factor, there should be one
## set of replications for each level
stopifnot(length(u.lev)==ncol(boot.frame))

all.bp\$stats <- matrix(data = NA, nrow = 5, ncol = length(u.lev))
all.bp\$conf <- matrix(data = NA, nrow = 2, ncol = length(u.lev))

for (i in 1:length(u.lev)){
t.bp <- boxplot.stats(boot.frame[,i])
all.bp\$stats[,i] <- t.bp\$stats
all.bp\$conf[,i] <- t.bp\$conf
all.bp\$out <- c(all.bp\$out,t.bp\$out)
all.bp\$group <- c(all.bp\$group, rep.int(i,length(t.bp\$out)))
}
all.bp\$n <- rep.int(plot.errors.boot.num, length(u.lev))

if(slice.index <= ncol(xdat))
all.bp\$names <- bws\$xdati\$all.lev[[slice.index]]
else
all.bp\$names <- bws\$zdati\$all.lev[[slice.index-ncol(xdat)]]
rm(boot.frame)
}

if (plot.errors.type == "standard") {
boot.err[,1:2] = 2.0*sqrt(diag(cov(boot.out\$t)))
}
else if (plot.errors.type == "quantiles") {
boot.err[,1:2] = t(sapply(as.data.frame(boot.out\$t),
function (y) {
quantile(y,probs = plot.errors.quantiles)
}))
boot.err[,1] = boot.out\$t0 - boot.err[,1]
boot.err[,2] = boot.err[,2] - boot.out\$t0
}
if (plot.errors.center == "bias-corrected")
boot.err[,3] <- 2*boot.out\$t0-colMeans(boot.out\$t)
list(boot.err = boot.err, bxp = all.bp)
}

compute.bootstrap.errors.plbandwidth =
function(xdat, ydat, zdat,
exdat, ezdat,
slice.index,
plot.errors.boot.method,
plot.errors.boot.blocklen,
plot.errors.boot.num,
plot.errors.center,
plot.errors.type,
plot.errors.quantiles,
bws){
boot.err = matrix(data = NA, nrow = dim(exdat)[1], ncol = 3)

is.inid = plot.errors.boot.method=="inid"

strf = ifelse(is.inid, "function(data,indices){", "function(tsb){")
strtx = ifelse(is.inid, "txdat = xdat[indices,],",
"txdat = tsb[,1:ncol(xdat),drop=FALSE],")
strty = ifelse(is.inid, "tydat = ydat[indices],",
"tydat = tsb[,ncol(xdat)+1],")
strtz = ifelse(is.inid, "tzdat = zdat[indices,],",
"tzdat = tsb[,(ncol(xdat)+2):ncol(tsb), drop=FALSE],")

boofun = eval(parse(text=paste(strf, "npplreg(", strtx, strty, strtz,
"exdat = exdat, ezdat = ezdat, bws = bws)\$mean}", sep="")))

if (is.inid){
boot.out = boot(data = data.frame(xdat,ydat,zdat), statistic = boofun,
R = plot.errors.boot.num)
} else {
boot.out = tsboot(tseries = data.frame(xdat,ydat,zdat), statistic = boofun,
R = plot.errors.boot.num,
l = plot.errors.boot.blocklen,
sim = plot.errors.boot.method)
}

all.bp <- list()

if (slice.index <= bws\$xndim){
tdati <- bws\$xdati
ti <- slice.index
} else {
tdati <- bws\$zdati
ti <- slice.index - bws\$xndim
}

if (slice.index > 0 && (tdati\$iord | tdati\$iuno)[ti]){
boot.frame <- as.data.frame(boot.out\$t)
u.lev <- tdati\$all.ulev[[ti]]

## if we are bootstrapping a factor, there should be one
## set of replications for each level
stopifnot(length(u.lev)==ncol(boot.frame))

all.bp\$stats <- matrix(data = NA, nrow = 5, ncol = length(u.lev))
all.bp\$conf <- matrix(data = NA, nrow = 2, ncol = length(u.lev))

for (i in 1:length(u.lev)){
t.bp <- boxplot.stats(boot.frame[,i])
all.bp\$stats[,i] <- t.bp\$stats
all.bp\$conf[,i] <- t.bp\$conf
all.bp\$out <- c(all.bp\$out,t.bp\$out)
all.bp\$group <- c(all.bp\$group, rep.int(i,length(t.bp\$out)))
}
all.bp\$n <- rep.int(plot.errors.boot.num, length(u.lev))
all.bp\$names <- tdati\$all.lev[[ti]]
rm(boot.frame)
}

if (plot.errors.type == "standard") {
boot.err[,1:2] = 2.0*sqrt(diag(cov(boot.out\$t)))
}
else if (plot.errors.type == "quantiles") {
boot.err[,1:2] = t(sapply(as.data.frame(boot.out\$t),
function (y) {
quantile(y,probs = plot.errors.quantiles)
}))
boot.err[,1] = boot.out\$t0 - boot.err[,1]
boot.err[,2] = boot.err[,2] - boot.out\$t0
}
if (plot.errors.center == "bias-corrected")
boot.err[,3] <- 2*boot.out\$t0-colMeans(boot.out\$t)
list(boot.err = boot.err, bxp = all.bp)
}

compute.bootstrap.errors.bandwidth =
function(xdat,
exdat,
cdf,
slice.index,
plot.errors.boot.method,
plot.errors.boot.blocklen,
plot.errors.boot.num,
plot.errors.center,
plot.errors.type,
plot.errors.quantiles,
bws){
boot.err = matrix(data = NA, nrow = dim(exdat)[1], ncol = 3)

is.inid = plot.errors.boot.method=="inid"

strf = ifelse(is.inid, "function(data,indices){", "function(tsb){")
strt = ifelse(is.inid, "tdat = xdat[indices,],",
"tdat = tsb,")

boofun = eval(parse(text=paste(strf,
ifelse(cdf, "npudist(", "npudens("),
strt,
"edat = exdat, bws = bws)\$",
ifelse(cdf, "dist", "dens"),
"}", sep="")))

if (is.inid) {
boot.out = boot(data = data.frame(xdat), statistic = boofun,
R = plot.errors.boot.num)
} else {
boot.out = tsboot(tseries = data.frame(xdat), statistic = boofun,
R = plot.errors.boot.num,
l = plot.errors.boot.blocklen,
sim = plot.errors.boot.method)
}

all.bp <- list()

if (slice.index > 0 && (bws\$xdati\$iord | bws\$xdati\$iuno)[slice.index]){
boot.frame <- as.data.frame(boot.out\$t)
u.lev <- bws\$xdati\$all.ulev[[slice.index]]

## if we are bootstrapping a factor, there should be one
## set of replications for each level
stopifnot(length(u.lev)==ncol(boot.frame))

all.bp\$stats <- matrix(data = NA, nrow = 5, ncol = length(u.lev))
all.bp\$conf <- matrix(data = NA, nrow = 2, ncol = length(u.lev))

for (i in 1:length(u.lev)){
t.bp <- boxplot.stats(boot.frame[,i])
all.bp\$stats[,i] <- t.bp\$stats
all.bp\$conf[,i] <- t.bp\$conf
all.bp\$out <- c(all.bp\$out,t.bp\$out)
all.bp\$group <- c(all.bp\$group, rep.int(i,length(t.bp\$out)))
}
all.bp\$n <- rep.int(plot.errors.boot.num, length(u.lev))
all.bp\$names <- bws\$xdati\$all.lev[[slice.index]]
rm(boot.frame)
}

if (plot.errors.type == "standard") {
boot.err[,1:2] = 2.0*sqrt(diag(cov(boot.out\$t)))
}
else if (plot.errors.type == "quantiles") {
boot.err[,1:2] = t(sapply(as.data.frame(boot.out\$t),
function (y) {
quantile(y,probs = plot.errors.quantiles)
}))
boot.err[,1] = boot.out\$t0 - boot.err[,1]
boot.err[,2] = boot.err[,2] - boot.out\$t0
}
if (plot.errors.center == "bias-corrected")
boot.err[,3] <- 2*boot.out\$t0-colMeans(boot.out\$t)
list(boot.err = boot.err, bxp = all.bp)
}

compute.bootstrap.errors.dbandwidth <-
function(xdat,
exdat,
slice.index,
plot.errors.boot.method,
plot.errors.boot.blocklen,
plot.errors.boot.num,
plot.errors.center,
plot.errors.type,
plot.errors.quantiles,
bws){
boot.err = matrix(data = NA, nrow = dim(exdat)[1], ncol = 3)

is.inid = plot.errors.boot.method=="inid"

strf = ifelse(is.inid, "function(data,indices){", "function(tsb){")
strt = ifelse(is.inid, "tdat = xdat[indices,],",
"tdat = tsb,")

boofun = eval(parse(text=paste(strf, "npudist(",
strt, "edat = exdat, bws = bws)\$dist}", sep="")))

if (is.inid) {
boot.out = boot(data = data.frame(xdat), statistic = boofun,
R = plot.errors.boot.num)
} else {
boot.out = tsboot(tseries = data.frame(xdat), statistic = boofun,
R = plot.errors.boot.num,
l = plot.errors.boot.blocklen,
sim = plot.errors.boot.method)
}

all.bp <- list()

if (slice.index > 0 && (bws\$xdati\$iord | bws\$xdati\$iuno)[slice.index]){
boot.frame <- as.data.frame(boot.out\$t)
u.lev <- bws\$xdati\$all.ulev[[slice.index]]

## if we are bootstrapping a factor, there should be one
## set of replications for each level
stopifnot(length(u.lev)==ncol(boot.frame))

all.bp\$stats <- matrix(data = NA, nrow = 5, ncol = length(u.lev))
all.bp\$conf <- matrix(data = NA, nrow = 2, ncol = length(u.lev))

for (i in 1:length(u.lev)){
t.bp <- boxplot.stats(boot.frame[,i])
all.bp\$stats[,i] <- t.bp\$stats
all.bp\$conf[,i] <- t.bp\$conf
all.bp\$out <- c(all.bp\$out,t.bp\$out)
all.bp\$group <- c(all.bp\$group, rep.int(i,length(t.bp\$out)))
}
all.bp\$n <- rep.int(plot.errors.boot.num, length(u.lev))
all.bp\$names <- bws\$xdati\$all.lev[[slice.index]]
rm(boot.frame)
}

if (plot.errors.type == "standard") {
boot.err[,1:2] = 2.0*sqrt(diag(cov(boot.out\$t)))
}
else if (plot.errors.type == "quantiles") {
boot.err[,1:2] = t(sapply(as.data.frame(boot.out\$t),
function (y) {
quantile(y,probs = plot.errors.quantiles)
}))
boot.err[,1] = boot.out\$t0 - boot.err[,1]
boot.err[,2] = boot.err[,2] - boot.out\$t0
}
if (plot.errors.center == "bias-corrected")
boot.err[,3] <- 2*boot.out\$t0-colMeans(boot.out\$t)
list(boot.err = boot.err, bxp = all.bp)
}

compute.bootstrap.errors.conbandwidth =
function(xdat, ydat,
exdat, eydat,
cdf,
quantreg,
tau,
slice.index,
plot.errors.boot.method,
plot.errors.boot.blocklen,
plot.errors.boot.num,
plot.errors.center,
plot.errors.type,
plot.errors.quantiles,
bws){
exdat = toFrame(exdat)
boot.err = matrix(data = NA, nrow = dim(exdat)[1], ncol = 3)

tboo =
if(quantreg) "quant"
else if (cdf) "dist"
else "dens"

is.inid = plot.errors.boot.method=="inid"

strf = ifelse(is.inid, "function(data,indices){", "function(tsb){")
strtx = ifelse(is.inid, "txdat = xdat[indices,],",
"txdat = tsb[,1:ncol(xdat),drop=FALSE],")
strty = ifelse(is.inid, "tydat = ydat[indices,],",
"tydat = tsb[,(ncol(xdat)+1):ncol(tsb), drop=FALSE],")

boofun = eval(parse(text=paste(strf,
switch(tboo,
"quant" = "npqreg(",
"dist" = "npcdist(",
"dens" = "npcdens("),
strtx, strty,
"exdat = exdat,",
ifelse(quantreg, "tau = tau", "eydat = eydat"),
switch(tboo,
"}", sep="")))
if (is.inid){
boot.out = boot(data = data.frame(xdat,ydat), statistic = boofun,
R = plot.errors.boot.num)
} else {
boot.out = tsboot(tseries = data.frame(xdat,ydat), statistic = boofun,
R = plot.errors.boot.num,
l = plot.errors.boot.blocklen,
sim = plot.errors.boot.method)
}

all.bp <- list()

if (slice.index <= bws\$xndim){
tdati <- bws\$xdati
ti <- slice.index
} else {
tdati <- bws\$ydati
ti <- slice.index - bws\$xndim
}

if (slice.index > 0 && (tdati\$iord | tdati\$iuno)[ti]){
boot.frame <- as.data.frame(boot.out\$t)
u.lev <- tdati\$all.ulev[[ti]]

## if we are bootstrapping a factor, there should be one
## set of replications for each level
stopifnot(length(u.lev)==ncol(boot.frame))

all.bp\$stats <- matrix(data = NA, nrow = 5, ncol = length(u.lev))
all.bp\$conf <- matrix(data = NA, nrow = 2, ncol = length(u.lev))

for (i in 1:length(u.lev)){
t.bp <- boxplot.stats(boot.frame[,i])
all.bp\$stats[,i] <- t.bp\$stats
all.bp\$conf[,i] <- t.bp\$conf
all.bp\$out <- c(all.bp\$out,t.bp\$out)
all.bp\$group <- c(all.bp\$group, rep.int(i,length(t.bp\$out)))
}
all.bp\$n <- rep.int(plot.errors.boot.num, length(u.lev))
all.bp\$names <- tdati\$all.lev[[ti]]
rm(boot.frame)
}

if (plot.errors.type == "standard") {
boot.err[,1:2] = 2.0*sqrt(diag(cov(boot.out\$t)))
}
else if (plot.errors.type == "quantiles") {
boot.err[,1:2] = t(sapply(as.data.frame(boot.out\$t),
function (y) {
quantile(y,probs = plot.errors.quantiles)
}))
boot.err[,1] = boot.out\$t0 - boot.err[,1]
boot.err[,2] = boot.err[,2] - boot.out\$t0
}
if (plot.errors.center == "bias-corrected")
boot.err[,3] <- 2*boot.out\$t0-colMeans(boot.out\$t)
list(boot.err = boot.err, bxp = all.bp)
}

compute.bootstrap.errors.condbandwidth =
function(xdat, ydat,
exdat, eydat,
cdf,
quantreg,
tau,
slice.index,
plot.errors.boot.method,
plot.errors.boot.blocklen,
plot.errors.boot.num,
plot.errors.center,
plot.errors.type,
plot.errors.quantiles,
bws){
exdat = toFrame(exdat)
boot.err = matrix(data = NA, nrow = dim(exdat)[1], ncol = 3)

tboo =
if(quantreg) "quant"
else if (cdf) "dist"
else "dens"

is.inid = plot.errors.boot.method=="inid"

strf = ifelse(is.inid, "function(data,indices){", "function(tsb){")
strtx = ifelse(is.inid, "txdat = xdat[indices,],",
"txdat = tsb[,1:ncol(xdat),drop=FALSE],")
strty = ifelse(is.inid, "tydat = ydat[indices,],",
"tydat = tsb[,(ncol(xdat)+1):ncol(tsb), drop=FALSE],")

boofun = eval(parse(text=paste(strf,
switch(tboo,
"quant" = "npqreg(",
"dist" = "npcdist(",
"dens" = "npcdens("),
strtx, strty,
"exdat = exdat,",
ifelse(quantreg, "tau = tau", "eydat = eydat"),
switch(tboo,
"}", sep="")))
if (is.inid){
boot.out = boot(data = data.frame(xdat,ydat), statistic = boofun,
R = plot.errors.boot.num)
} else {
boot.out = tsboot(tseries = data.frame(xdat,ydat), statistic = boofun,
R = plot.errors.boot.num,
l = plot.errors.boot.blocklen,
sim = plot.errors.boot.method)
}

all.bp <- list()

if (slice.index <= bws\$xndim){
tdati <- bws\$xdati
ti <- slice.index
} else {
tdati <- bws\$ydati
ti <- slice.index - bws\$xndim
}

if (slice.index > 0 && (tdati\$iord | tdati\$iuno)[ti]){
boot.frame <- as.data.frame(boot.out\$t)
u.lev <- tdati\$all.ulev[[ti]]

## if we are bootstrapping a factor, there should be one
## set of replications for each level
stopifnot(length(u.lev)==ncol(boot.frame))

all.bp\$stats <- matrix(data = NA, nrow = 5, ncol = length(u.lev))
all.bp\$conf <- matrix(data = NA, nrow = 2, ncol = length(u.lev))

for (i in 1:length(u.lev)){
t.bp <- boxplot.stats(boot.frame[,i])
all.bp\$stats[,i] <- t.bp\$stats
all.bp\$conf[,i] <- t.bp\$conf
all.bp\$out <- c(all.bp\$out,t.bp\$out)
all.bp\$group <- c(all.bp\$group, rep.int(i,length(t.bp\$out)))
}
all.bp\$n <- rep.int(plot.errors.boot.num, length(u.lev))
all.bp\$names <- tdati\$all.lev[[ti]]
rm(boot.frame)
}

if (plot.errors.type == "standard") {
boot.err[,1:2] = 2.0*sqrt(diag(cov(boot.out\$t)))
}
else if (plot.errors.type == "quantiles") {
boot.err[,1:2] = t(sapply(as.data.frame(boot.out\$t),
function (y) {
quantile(y,probs = plot.errors.quantiles)
}))
boot.err[,1] = boot.out\$t0 - boot.err[,1]
boot.err[,2] = boot.err[,2] - boot.out\$t0
}
if (plot.errors.center == "bias-corrected")
boot.err[,3] <- 2*boot.out\$t0-colMeans(boot.out\$t)
list(boot.err = boot.err, bxp = all.bp)
}

compute.bootstrap.errors.sibandwidth =
function(xdat, ydat,
plot.errors.boot.method,
plot.errors.boot.blocklen,
plot.errors.boot.num,
plot.errors.center,
plot.errors.type,
plot.errors.quantiles,
bws){

boot.err = matrix(data = NA, nrow = nrow(xdat), ncol = 3)

is.inid = plot.errors.boot.method=="inid"

strf = ifelse(is.inid, "function(data,indices){", "function(tsb){")
strtx = ifelse(is.inid, "txdat = xdat[indices,],",
"txdat = tsb[,1:(ncol(tsb)-1),drop=FALSE],")
strty = ifelse(is.inid, "tydat = ydat[indices],",
"tydat = tsb[,ncol(tsb)],")

## beta[1] is always 1.0, so use first column of gradients matrix ...

boofun = eval(parse(text=paste(strf, "npindex(", strtx, strty,
"exdat = xdat, bws = bws,",

if (is.inid){
boot.out = boot(data = data.frame(xdat,ydat), statistic = boofun,
R = plot.errors.boot.num)
} else {
boot.out = tsboot(tseries = data.frame(xdat,ydat), statistic = boofun,
R = plot.errors.boot.num,
l = plot.errors.boot.blocklen,
sim = plot.errors.boot.method)
}

if (plot.errors.type == "standard") {
boot.err[,1:2] = 2.0*sqrt(diag(cov(boot.out\$t)))
}
else if (plot.errors.type == "quantiles") {
boot.err[,1:2] = t(sapply(as.data.frame(boot.out\$t),
function (y) {
quantile(y,probs = plot.errors.quantiles)
}))
boot.err[,1] = boot.out\$t0 - boot.err[,1]
boot.err[,2] = boot.err[,2] - boot.out\$t0
}
if (plot.errors.center == "bias-corrected")
boot.err[,3] <- 2*boot.out\$t0-colMeans(boot.out\$t)
boot.err
}

uocquantile <- function(x, prob) {
if(any(prob < 0 | prob > 1)) stop("'prob' outside [0,1]")
if(any(is.na(x) | is.nan(x))) stop("missing values and NaN's not allowed")
if (is.ordered(x)){
tq = unclass(table(x))
tq = tq / sum(tq)
tq[length(tq)] <- 1.0
bscape <- sort(unique(x))
tq <- sapply(1:length(tq), function(y){ sum(tq[1:y]) })
j <- sapply(prob, function(p){ which(tq >= p)[1] })
bscape[j]
} else if (is.factor(x)) {
## just returns mode
tq = unclass(table(x))
j = which(tq == max(tq))[1]
sort(unique(x))[j]
} else {
quantile(x, probs = prob)
}
}

trim.quantiles <- function(dat, trim){
if (sign(trim) == sign(-1)){
trim <- abs(trim)
tq <- quantile(dat, probs = c(0.0, 0.0+trim, 1.0-trim,1.0))
tq <- c(2.0*tq[1]-tq[2], 2.0*tq[4]-tq[3])
}
else {
tq <- quantile(dat, probs = c(0.0+trim, 1.0-trim))
}
tq
}

npplot <- function(bws = stop("'bws' has not been set"), ..., random.seed = 42){
## Save seed prior to setting
if(exists(".Random.seed", .GlobalEnv)) {
save.seed <- get(".Random.seed", .GlobalEnv)
exists.seed = TRUE
} else {
exists.seed = FALSE
}

set.seed(random.seed)

## Restore seed
on.exit(if(exists.seed) assign(".Random.seed", save.seed, .GlobalEnv))

UseMethod("npplot",bws)
}

npplot.rbandwidth <-
function(bws,
xdat,
ydat,
data = NULL,
xq = 0.5,
xtrim = 0.0,
neval = 50,
common.scale = TRUE,
perspective = TRUE,
main = NULL,
type = NULL,
border = NULL,
cex.axis = NULL,
cex.lab = NULL,
cex.main = NULL,
cex.sub = NULL,
col = NULL,
ylab = NULL,
xlab = NULL,
zlab = NULL,
sub = NULL,
ylim = NULL,
xlim = NULL,
zlim = NULL,
lty = NULL,
lwd = NULL,
theta = 0.0,
phi = 10.0,
view = c("rotate","fixed"),
plot.behavior = c("plot","plot-data","data"),
plot.errors.method = c("none","bootstrap","asymptotic"),
plot.errors.boot.num = 399,
plot.errors.boot.method = c("inid", "fixed", "geom"),
plot.errors.boot.blocklen = NULL,
plot.errors.center = c("estimate","bias-corrected"),
plot.errors.type = c("standard","quantiles"),
plot.errors.quantiles = c(0.025,0.975),
plot.errors.style = c("band","bar"),
plot.errors.bar = c("|","I"),
plot.errors.bar.num = min(neval,25),
plot.bxp = FALSE,
plot.bxp.out = TRUE,
plot.par.mfrow = TRUE,
...,
random.seed){

if(!is.null(options('plot.par.mfrow')\$plot.par.mfrow))
plot.par.mfrow <- options('plot.par.mfrow')\$plot.par.mfrow

miss.xy = c(missing(xdat),missing(ydat))

if (any(miss.xy) && !all(miss.xy))
stop("one of, but not both, xdat and ydat was specified")
else if(all(miss.xy) && !is.null(bws\$formula)){
tt <- terms(bws)
m <- match(c("formula", "data", "subset", "na.action"),
names(bws\$call), nomatch = 0)
tmf <- bws\$call[c(1,m)]
tmf[[1]] <- as.name("model.frame")
tmf[["formula"]] <- tt
umf <- tmf <- eval(tmf, envir = environment(tt))

ydat <- model.response(tmf)
xdat <- tmf[, attr(attr(tmf, "terms"),"term.labels"), drop = FALSE]
} else {
if(all(miss.xy) && !is.null(bws\$call)){
xdat <- data.frame(eval(bws\$call[["xdat"]], environment(bws\$call)))
ydat = eval(bws\$call[["ydat"]], environment(bws\$call))
}

## catch and destroy NA's
xdat = toFrame(xdat)
goodrows = 1:dim(xdat)[1]
rows.omit = attr(na.omit(data.frame(xdat,ydat)), "na.action")
goodrows[rows.omit] = 0

if (all(goodrows==0))
stop("Data has no rows without NAs")

xdat = xdat[goodrows,,drop = FALSE]
ydat = ydat[goodrows]
}

## ydat = as.double(ydat)

xq = double(bws\$ndim)+xq
xtrim = double(bws\$ndim)+xtrim

if (missing(plot.errors.method) &
any(!missing(plot.errors.boot.num), !missing(plot.errors.boot.method),
!missing(plot.errors.boot.blocklen))){
warning(paste("plot.errors.method must be set to 'bootstrap' to use bootstrapping.",
"\nProceeding without bootstrapping."))
}

plot.behavior = match.arg(plot.behavior)
plot.errors.method = match.arg(plot.errors.method)
plot.errors.boot.method = match.arg(plot.errors.boot.method)
plot.errors.center = match.arg(plot.errors.center)
plot.errors.type = match.arg(plot.errors.type)
plot.errors.style = match.arg(plot.errors.style)
plot.errors.bar = match.arg(plot.errors.bar)

common.scale = common.scale | (!is.null(ylim))

if (plot.errors.method == "asymptotic") {
if (plot.errors.type == "quantiles"){
warning("quantiles cannot be calculated with asymptotics, calculating standard errors")
plot.errors.type = "standard"
}

if (plot.errors.center == "bias-corrected") {
warning("no bias corrections can be calculated with asymptotics, centering on estimate")
plot.errors.center = "estimate"
}
}

if (is.element(plot.errors.boot.method, c("fixed", "geom")) &&
is.null(plot.errors.boot.blocklen))
plot.errors.boot.blocklen = b.star(xdat,round=TRUE)[1,1]

plot.errors = (plot.errors.method != "none")

if ((bws\$ncon + bws\$nord == 2) & (bws\$nuno == 0) & perspective & !gradients &
!any(xor(bws\$xdati\$iord, bws\$xdati\$inumord))){

view = match.arg(view)
rotate = (view == "rotate")

if (is.ordered(xdat[,1])){
x1.eval = bws\$xdati\$all.ulev[[1]]
x1.neval = length(x1.eval)
} else {
x1.neval = neval
qi = trim.quantiles(xdat[,1], xtrim[1])
x1.eval = seq(qi[1], qi[2], length.out = x1.neval)
}

if (is.ordered(xdat[,2])){
x2.eval = bws\$xdati\$all.ulev[[2]]
x2.neval = length(x2.eval)
} else {
x2.neval = neval
qi = trim.quantiles(xdat[,2], xtrim[2])
x2.eval = seq(qi[1], qi[2], length.out = x2.neval)
}

x.eval <- expand.grid(x1.eval, x2.eval)

if (is.ordered(xdat[,1]))
x1.eval <- (bws\$xdati\$all.dlev[[1]])[as.integer(x1.eval)]

if (is.ordered(xdat[,2]))
x2.eval <- (bws\$xdati\$all.dlev[[2]])[as.integer(x2.eval)]

tobj = npreg(txdat = xdat, tydat = ydat,
exdat = x.eval, bws = bws)

terr = matrix(data = tobj\$merr, nrow = dim(x.eval)[1], ncol = 3)
terr[,3] = NA

treg = matrix(data = tobj\$mean,
nrow = x1.neval, ncol = x2.neval, byrow = FALSE)

if (plot.errors.method == "bootstrap"){
terr <- compute.bootstrap.errors(xdat = xdat, ydat = ydat,
exdat = x.eval,
slice.index = 0,
plot.errors.boot.method = plot.errors.boot.method,
plot.errors.boot.blocklen = plot.errors.boot.blocklen,
plot.errors.boot.num = plot.errors.boot.num,
plot.errors.center = plot.errors.center,
plot.errors.type = plot.errors.type,
plot.errors.quantiles = plot.errors.quantiles,
bws = bws)[["boot.err"]]

pc = (plot.errors.center == "bias-corrected")

lerr = matrix(data = if(pc) {terr[,3]} else {tobj\$mean}
-terr[,1],
nrow = x1.neval, ncol = x2.neval, byrow = FALSE)

herr = matrix(data = if(pc) {terr[,3]} else {tobj\$mean}
+terr[,2],
nrow = x1.neval, ncol = x2.neval, byrow = FALSE)

} else if (plot.errors.method == "asymptotic") {
lerr = matrix(data = tobj\$mean - 2.0*tobj\$merr,
nrow = x1.neval, ncol = x2.neval, byrow = FALSE)

herr = matrix(data = tobj\$mean + 2.0*tobj\$merr,
nrow = x1.neval, ncol = x2.neval, byrow = FALSE)

}

if(is.null(zlim)) {
zlim =
if (plot.errors)
c(min(lerr),max(herr))
else
c(min(tobj\$mean),max(tobj\$mean))
}

if (plot.behavior != "plot"){
r1 = npregression(bws = bws,
eval = x.eval,
mean = as.double(treg),
merr = terr[,1:2],
ntrain = dim(xdat)[1])
r1\$bias = NA

if (plot.errors.center == "bias-corrected")
r1\$bias = terr[,3] - treg

if (plot.behavior == "data")
return ( list(r1 = r1) )

}

dtheta = 5.0
dphi = 10.0

persp.col = ifelse(plot.errors, FALSE, ifelse(!is.null(col),col,"lightblue"))

for (i in 0:((360 %/% dtheta - 1)*rotate)*dtheta+theta){
if (plot.errors){
persp(x1.eval,
x2.eval,
lerr,
zlim = zlim,
cex.axis = ifelse(!is.null(cex.axis),cex.axis,par()\$cex.axis),
cex.lab = ifelse(!is.null(cex.lab),cex.lab,par()\$cex.lab),
cex.main = ifelse(!is.null(cex.main),cex.main,par()\$cex.main),
cex.sub = ifelse(!is.null(cex.sub),cex.sub,par()\$cex.sub),
col = persp.col,
border = ifelse(!is.null(border),border,"grey"),
ticktype = "detailed",
xlab = "",
ylab = "",
zlab = "",
theta = i,
phi = phi,
lwd = ifelse(!is.null(lwd),lwd,par()\$lwd))
par(new = TRUE)
}

persp(x1.eval,
x2.eval,
treg,
zlim = zlim,
col = persp.col,
border = ifelse(!is.null(border),border,"black"),
ticktype = "detailed",
cex.axis = ifelse(!is.null(cex.axis),cex.axis,par()\$cex.axis),
cex.lab =  ifelse(!is.null(cex.lab),cex.lab,par()\$cex.lab),
cex.main = ifelse(!is.null(cex.main),cex.main,par()\$cex.main),
cex.sub = ifelse(!is.null(cex.sub),cex.sub,par()\$cex.sub),
xlab = ifelse(!is.null(xlab),xlab,gen.label(bws\$xnames[1], "X1")),
ylab = ifelse(!is.null(ylab),ylab,gen.label(bws\$xnames[2], "X2")),
zlab = ifelse(!is.null(zlab),zlab,gen.label(bws\$ynames,"Conditional Mean")),
theta = i,
phi = phi,
main = gen.tflabel(!is.null(main), main, paste("[theta= ", i,", phi= ", phi,"]", sep="")))

if (plot.errors){
par(new = TRUE)
persp(x1.eval,
x2.eval,
herr,
zlim = zlim,
cex.axis = ifelse(!is.null(cex.axis),cex.axis,par()\$cex.axis),
cex.lab = ifelse(!is.null(cex.lab),cex.lab,par()\$cex.lab),
cex.main = ifelse(!is.null(cex.main),cex.main,par()\$cex.main),
cex.sub = ifelse(!is.null(cex.sub),cex.sub,par()\$cex.sub),
col = persp.col,
border = ifelse(!is.null(border),border,"grey"),
ticktype = "detailed",
xlab = "",
ylab = "",
zlab = "",
theta = i,
phi = phi,
lwd = ifelse(!is.null(lwd),lwd,par()\$lwd))
}

Sys.sleep(0.5)
}

if (plot.behavior == "plot-data")
return ( list(r1 = r1) )

} else {

if (plot.behavior != "data" && plot.par.mfrow)
par(mfrow=n2mfrow(bws\$ndim),cex=par()\$cex)

ev = xdat[1,,drop = FALSE]

for (i in 1:bws\$ndim)
ev[1,i] = uocquantile(xdat[,i], prob=xq[i])

maxneval = max(c(sapply(xdat,nlevels),neval))

exdat = as.data.frame(matrix(data = 0, nrow = maxneval, ncol = bws\$ndim))

for (i in 1:bws\$ndim)
exdat[,i] = ev[1,i]

if (common.scale){
data.eval = matrix(data = NA, nrow = maxneval, ncol = bws\$ndim)
data.err = matrix(data = NA, nrow = maxneval, ncol = 3*bws\$ndim)
allei = as.data.frame(matrix(data = NA, nrow = maxneval, ncol = bws\$ndim))
all.bxp = list()
}

plot.out = list()

temp.err = matrix(data = NA, nrow = maxneval, ncol = 3)
temp.mean = replicate(maxneval, NA)

## plotting expressions

plot.bootstrap = plot.errors.method == "bootstrap"

pfunE = expression(ifelse(xi.factor,
ifelse(plot.bootstrap & plot.bxp,"bxp","plotFactor"), "plot"))
pxE = expression(ifelse(common.scale,
ifelse(xi.factor,
ifelse(plot.bootstrap & plot.bxp, "z = all.bxp[[i]],", "f = allei[,i],"),
"x = allei[,i],"),
ifelse(xi.factor,
ifelse(plot.bootstrap & plot.bxp, "z = temp.boot,", "f = ei,"), "x = ei,")))

pyE = expression(ifelse(xi.factor & plot.bootstrap & plot.bxp, "",
ifelse(common.scale,"y = data.eval[,i],", "y = temp.mean,")))

pylimE = ifelse(common.scale, "ylim = c(y.min,y.max),",
ifelse(plot.errors, "ylim = c(min(na.omit(c(temp.mean - temp.err[,1], temp.err[,3] - temp.err[,1]))),
max(na.omit(c(temp.mean + temp.err[,2], temp.err[,3] + temp.err[,2])))),", ""))

pxlabE = "xlab = ifelse(!is.null(xlab),xlab,gen.label(bws\$xnames[i], paste('X', i, sep = ''))),"
paste('Gradient Component ', i, ' of', sep=''), ''),
gen.label(bws\$ynames, 'Conditional Mean'))),"

prestE = expression(ifelse(xi.factor,"", "type = ifelse(!is.null(type),type,'l'), lty = ifelse(!is.null(lty),lty,par()\$lty), col = ifelse(!is.null(col),col,par()\$col), lwd = ifelse(!is.null(lwd),lwd,par()\$lwd), cex.axis = ifelse(!is.null(cex.axis),cex.axis,par()\$cex.axis), cex.lab = ifelse(!is.null(cex.lab),cex.lab,par()\$cex.lab), cex.main = ifelse(!is.null(cex.main),cex.main,par()\$cex.main), cex.sub = ifelse(!is.null(cex.sub),cex.sub,par()\$cex.sub),"))

pmainE = "main = ifelse(!is.null(main),main,''), sub = ifelse(!is.null(sub),sub,''),"

## error plotting expressions
plotOnEstimate = (plot.errors.center == "estimate")

efunE = "draw.errors"
eexE = expression(ifelse(common.scale, "ex = as.numeric(na.omit(allei[,i])),",
"ex = as.numeric(na.omit(ei)),"))
eelyE = expression(ifelse(common.scale,
ifelse(plotOnEstimate, "ely = na.omit(data.eval[,i] - data.err[,3*i-2]),",
"ely = na.omit(data.err[,3*i] - data.err[,3*i-2]),"),
ifelse(plotOnEstimate, "ely = na.omit(temp.mean - temp.err[,1]),",
"ely = na.omit(temp.err[,3] - temp.err[,1]),")))
eehyE = expression(ifelse(common.scale,
ifelse(plotOnEstimate, "ehy = na.omit(data.eval[,i] + data.err[,3*i-1]),",
"ehy = na.omit(data.err[,3*i] + data.err[,3*i-1]),"),
ifelse(plotOnEstimate, "ehy = na.omit(temp.mean + temp.err[,2]),",
"ehy = na.omit(temp.err[,3] + temp.err[,2]),")))

erestE = "plot.errors.style = ifelse(xi.factor,'bar',plot.errors.style),
plot.errors.bar = ifelse(xi.factor,'I',plot.errors.bar),
plot.errors.bar.num = plot.errors.bar.num,
lty = ifelse(xi.factor,1,2)"

for (i in 1:bws\$ndim){
temp.err[,] = NA
temp.mean[] =  NA
temp.boot = list()

xi.factor = is.factor(xdat[,i])

if (xi.factor){
ei = bws\$xdati\$all.ulev[[i]]
xi.neval = length(ei)
} else {
xi.neval = neval
qi = trim.quantiles(xdat[,i], xtrim[i])
ei = seq(qi[1], qi[2], length.out = neval)
}

if (xi.neval < maxneval){
ei[(xi.neval+1):maxneval] = NA
}

tr = npreg(txdat = xdat, tydat = ydat,
exdat = subcol(exdat,ei,i)[1:xi.neval,, drop = FALSE], bws = bws,

if (plot.errors){
if (plot.errors.method == "asymptotic")
temp.err[1:xi.neval,1:2] = replicate(2,2.0*(if(gradients) tr\$gerr[,i] else tr\$merr))
else if (plot.errors.method == "bootstrap"){
temp.boot <- compute.bootstrap.errors(
xdat = xdat, ydat = ydat,
exdat = subcol(exdat,ei,i)[1:xi.neval,, drop = FALSE],
slice.index = i,
plot.errors.boot.method = plot.errors.boot.method,
plot.errors.boot.blocklen = plot.errors.boot.blocklen,
plot.errors.boot.num = plot.errors.boot.num,
plot.errors.center = plot.errors.center,
plot.errors.type = plot.errors.type,
plot.errors.quantiles = plot.errors.quantiles,
bws = bws)
temp.err[1:xi.neval,] = temp.boot[["boot.err"]]
temp.boot <- temp.boot[["bxp"]]
if (!plot.bxp.out){
temp.boot\$out <- numeric()
temp.boot\$group <- integer()
}
}
}

if (common.scale){
allei[,i] = ei
data.eval[,i] = temp.mean
if (plot.errors){
all.bxp[i] = NA
all.bxp[[i]] = temp.boot

data.err[,c(3*i-2,3*i-1,3*i)] = temp.err
}
} else if (plot.behavior != "data") {
## plot evaluation
eval(parse(text = paste(eval(pfunE), "(", eval(pxE), eval(pyE),
eval(pylimE), eval(pxlabE), eval(pylabE), eval(prestE),
eval(pmainE), ")")))

## error plotting evaluation
if (plot.errors && !(xi.factor & plot.bootstrap & plot.bxp)){
if (!xi.factor && !plotOnEstimate)
lines(na.omit(ei), na.omit(temp.err[,3]), lty = 3)

eval(parse(text = paste(eval(efunE), "(", eval(eexE), eval(eelyE),
eval(eehyE), eval(erestE), ")")))
}

}

if (plot.behavior != "plot") {
plot.out[i] = NA
plot.out[[i]] =
npregression(bws = bws,
eval = as.data.frame(subcol(exdat,ei,i)[1:xi.neval,]),
mean = tr\$mean,
merr = tr\$merr,
gerr = na.omit(cbind(-temp.err[,1],
temp.err[,2])),
ntrain = dim(xdat)[1])
plot.out[[i]]\$gbias = na.omit(temp.mean - temp.err[,3])
} else {
plot.out[[i]] =
npregression(bws = bws,
eval = as.data.frame(subcol(exdat,ei,i)[1:xi.neval,]),
mean = na.omit(temp.mean),
merr = na.omit(cbind(-temp.err[,1],
temp.err[,2])),
ntrain = dim(xdat)[1])
plot.out[[i]]\$bias = na.omit(temp.mean - temp.err[,3])
}
plot.out[[i]]\$bxp = temp.boot
}
}

if (common.scale & (plot.behavior != "data")){
jj = 1:bws\$ndim*3

if (plot.errors.center == "estimate" | !plot.errors) {
y.max = max(na.omit(as.double(data.eval)) +
if (plot.errors) na.omit(as.double(data.err[,jj-1]))
else 0
)
y.min = min(na.omit(as.double(data.eval)) -
if (plot.errors) na.omit(as.double(data.err[,jj-2]))
else 0
)
} else if (plot.errors.center == "bias-corrected") {
y.max = max(na.omit(as.double(data.err[,jj] + data.err[,jj-1])))
y.min = min(na.omit(as.double(data.err[,jj] - data.err[,jj-2])))
}

if(!is.null(ylim)){
y.min = ylim[1]
y.max = ylim[2]
}

for (i in 1:bws\$ndim){
xi.factor = is.factor(xdat[,i])

## plot evaluation
eval(parse(text = paste(eval(pfunE), "(", eval(pxE), eval(pyE),
eval(pylimE), eval(pxlabE), eval(pylabE), eval(prestE),
eval(pmainE), ")")))

## error plotting evaluation
if (plot.errors && !(xi.factor & plot.bootstrap & plot.bxp)){
if (!xi.factor && !plotOnEstimate)
lines(na.omit(allei[,i]), na.omit(data.err[,3*i]), lty = 3)

eval(parse(text = paste(eval(efunE), "(", eval(eexE), eval(eelyE),
eval(eehyE), eval(erestE), ")")))
}
}
}

if (plot.behavior != "data" && plot.par.mfrow)
par(mfrow=c(1,1),cex=par()\$cex)

if (plot.behavior != "plot"){
names(plot.out) =
paste("rg",1:bws\$ndim,sep="")
else
paste("r",1:bws\$ndim,sep="")

return (plot.out)
}

}
}

npplot.scbandwidth <-
function(bws,
xdat,
ydat,
zdat = NULL,
data = NULL,
xq = 0.5,
zq = 0.5,
xtrim = 0.0,
ztrim = 0.0,
neval = 50,
common.scale = TRUE,
perspective = TRUE,
main = NULL,
type = NULL,
border = NULL,
cex.axis = NULL,
cex.lab = NULL,
cex.main = NULL,
cex.sub = NULL,
col = NULL,
ylab = NULL,
xlab = NULL,
zlab = NULL,
sub = NULL,
ylim = NULL,
xlim = NULL,
zlim = NULL,
lty = NULL,
lwd = NULL,
theta = 0.0,
phi = 10.0,
view = c("rotate","fixed"),
plot.behavior = c("plot","plot-data","data"),
plot.errors.method = c("none","bootstrap","asymptotic"),
plot.errors.boot.num = 399,
plot.errors.boot.method = c("inid", "fixed", "geom"),
plot.errors.boot.blocklen = NULL,
plot.errors.center = c("estimate","bias-corrected"),
plot.errors.type = c("standard","quantiles"),
plot.errors.quantiles = c(0.025,0.975),
plot.errors.style = c("band","bar"),
plot.errors.bar = c("|","I"),
plot.errors.bar.num = min(neval,25),
plot.bxp = FALSE,
plot.bxp.out = TRUE,
plot.par.mfrow = TRUE,
...,
random.seed){

if(!is.null(options('plot.par.mfrow')\$plot.par.mfrow))
plot.par.mfrow <- options('plot.par.mfrow')\$plot.par.mfrow

stop("gradients not supported with smooth coefficient models.")

miss.xy = c(missing(xdat),missing(ydat))
miss.z = missing(zdat) & is.null(bws\$zdati)

if (any(miss.xy) && !all(miss.xy))
stop("one of, but not both, xdat and ydat was specified")
else if(all(miss.xy) && !is.null(bws\$formula)){
tt <- terms(bws)
m <- match(c("formula", "data", "subset", "na.action"),
names(bws\$call), nomatch = 0)
tmf <- bws\$call[c(1,m)]
tmf[[1]] <- as.name("model.frame")
tmf[["formula"]] <- tt
umf <- tmf <- eval(tmf, envir = environment(tt))

ydat <- model.response(tmf)
xdat <- tmf[, bws\$chromoly[[2]], drop = FALSE]
if (!miss.z)
zdat <- tmf[, bws\$chromoly[[3]], drop = FALSE]
} else {
if(all(miss.xy) && !is.null(bws\$call)){
xdat <- data.frame(eval(bws\$call[["xdat"]], environment(bws\$call)))
ydat <- eval(bws\$call[["ydat"]], environment(bws\$call))
if (!miss.z)
zdat <- data.frame(eval(bws\$call[["zdat"]], environment(bws\$call)))
}

## catch and destroy NA's
xdat = toFrame(xdat)
if(!miss.z)
zdat <- toFrame(zdat)

goodrows = 1:dim(xdat)[1]
rows.omit =
eval(parse(text = paste("attr(na.omit(data.frame(xdat, ydat",
ifelse(miss.z,'',',zdat'),')), "na.action")')))

attr(na.omit(data.frame(xdat,ydat,zdat)), "na.action")
goodrows[rows.omit] = 0

if (all(goodrows==0))
stop("Data has no rows without NAs")

xdat = xdat[goodrows,,drop = FALSE]

if(!miss.z)
zdat <- zdat[goodrows,,drop = FALSE]

ydat = ydat[goodrows]
}

## ydat = as.double(ydat)

xq = double(ncol(xdat))+xq
xtrim = double(ncol(xdat))+xtrim

if (!miss.z){
zq = double(ncol(zdat))+zq
ztrim = double(ncol(zdat))+ztrim
}

if (missing(plot.errors.method) &
any(!missing(plot.errors.boot.num), !missing(plot.errors.boot.method),
!missing(plot.errors.boot.blocklen))){
warning(paste("plot.errors.method must be set to 'bootstrap' to use bootstrapping.",
"\nProceeding without bootstrapping."))
}

plot.behavior = match.arg(plot.behavior)
plot.errors.method = match.arg(plot.errors.method)
plot.errors.boot.method = match.arg(plot.errors.boot.method)
plot.errors.center = match.arg(plot.errors.center)
plot.errors.type = match.arg(plot.errors.type)
plot.errors.style = match.arg(plot.errors.style)
plot.errors.bar = match.arg(plot.errors.bar)

common.scale = common.scale | (!is.null(ylim))

if (plot.errors.method == "asymptotic") {
if (plot.errors.type == "quantiles"){
warning("quantiles cannot be calculated with asymptotics, calculating standard errors")
plot.errors.type = "standard"
}

if (plot.errors.center == "bias-corrected") {
warning("no bias corrections can be calculated with asymptotics, centering on estimate")
plot.errors.center = "estimate"
}
}

if (is.element(plot.errors.boot.method, c("fixed", "geom")) &&
is.null(plot.errors.boot.blocklen))
plot.errors.boot.blocklen = b.star(xdat,round=TRUE)[1,1]

plot.errors = (plot.errors.method != "none")

if ((sum(c(bws\$xdati\$icon, bws\$xdati\$iord, bws\$zdati\$icon, bws\$zdati\$iord))== 2) & (sum(c(bws\$xdati\$iuno, bws\$zdati\$iuno)) == 0) & perspective & !gradients &
!any(xor(c(bws\$xdati\$iord, bws\$zdati\$iord), c(bws\$xdati\$inumord, bws\$zdati\$inumord)))){

view = match.arg(view)
rotate = (view == "rotate")

if (is.ordered(xdat[,1])){
x1.eval = bws\$xdati\$all.ulev[[1]]
x1.neval = length(x1.eval)
} else {
x1.neval = neval
qi = trim.quantiles(xdat[,1], xtrim[1])
x1.eval = seq(qi[1], qi[2], length.out = x1.neval)
}

if(!miss.z){
tdat <- zdat[,1]
ti <- 1
tdati <-  bws\$zdati
ttrim <- ztrim
x2.names <- bws\$znames
} else {
tdat <- xdat[,2]
ti <- 2
tdati <- bws\$xdati
ttrim <- xtrim
x2.names <- bws\$xnames
}

if (is.ordered(tdat)){
x2.eval = tdati\$all.ulev[[ti]]
x2.neval = length(x2.eval)
} else {
x2.neval = neval
qi = trim.quantiles(tdat, ttrim[ti])
x2.eval = seq(qi[1], qi[2], length.out = x2.neval)
}

x.eval <- expand.grid(x1.eval, x2.eval)

if (is.ordered(xdat[,1]))
x1.eval <- (bws\$xdati\$all.dlev[[1]])[as.integer(x1.eval)]

if (is.ordered(tdat))
x2.eval <- (tdati\$all.dlev[[ti]])[as.integer(x2.eval)]

tobj <- eval(parse(text = paste('npscoef(txdat = xdat, tydat = ydat,',
ifelse(miss.z,'','tzdat = zdat,'),
ifelse(miss.z,'exdat = x.eval,','exdat = x.eval[,1, drop = FALSE], ezdat = x.eval[,2, drop = FALSE],'),
'bws = bws, iterate = FALSE, errors = plot.errors)')))

terr = matrix(data = tobj\$merr, nrow = dim(x.eval)[1], ncol = 3)
terr[,3] = NA

treg = matrix(data = tobj\$mean,
nrow = x1.neval, ncol = x2.neval, byrow = FALSE)

if (plot.errors.method == "bootstrap"){
terr <-
eval(parse(text = paste(
'compute.bootstrap.errors(xdat = xdat, ydat = ydat,',
ifelse(miss.z,'','zdat = zdat,'),
ifelse(miss.z, 'exdat = x.eval,',
'exdat = x.eval[,1, drop = FALSE], ezdat = x.eval[,1, drop = FALSE],'),
' gradients = FALSE, slice.index = 0,',
'plot.errors.boot.method = plot.errors.boot.method,',
'plot.errors.boot.blocklen = plot.errors.boot.blocklen,',
'plot.errors.boot.num = plot.errors.boot.num,',
'plot.errors.center = plot.errors.center,',
'plot.errors.type = plot.errors.type,',
'plot.errors.quantiles = plot.errors.quantiles,',
'bws = bws)[["boot.err"]]')))

pc = (plot.errors.center == "bias-corrected")

lerr = matrix(data = if(pc) {terr[,3]} else {tobj\$mean}
-terr[,1],
nrow = x1.neval, ncol = x2.neval, byrow = FALSE)

herr = matrix(data = if(pc) {terr[,3]} else {tobj\$mean}
+terr[,2],
nrow = x1.neval, ncol = x2.neval, byrow = FALSE)

} else if (plot.errors.method == "asymptotic") {
lerr = matrix(data = tobj\$mean - 2.0*tobj\$merr,
nrow = x1.neval, ncol = x2.neval, byrow = FALSE)

herr = matrix(data = tobj\$mean + 2.0*tobj\$merr,
nrow = x1.neval, ncol = x2.neval, byrow = FALSE)

}

if(is.null(zlim)) {
zlim =
if (plot.errors)
c(min(lerr),max(herr))
else
c(min(tobj\$mean),max(tobj\$mean))
}

if (plot.behavior != "plot"){
r1 <-
eval(parse(text = paste("smoothcoefficient(bws = bws,",
ifelse(miss.z, "eval = x.eval,", "eval = list(exdat = x.eval[,1, drop = FALSE], ezdat = x.eval[,2, drop = FALSE])"),
"mean = as.double(treg),",
"merr = terr[,1:2],",
"ntrain = dim(xdat)[1])")))
r1\$bias = NA

if (plot.errors.center == "bias-corrected")
r1\$bias = terr[,3] - treg

if (plot.behavior == "data")
return ( list(r1 = r1) )

}

dtheta = 5.0
dphi = 10.0

persp.col = ifelse(plot.errors, FALSE, ifelse(!is.null(col),col,"lightblue"))

##for (j in 0:((50 %/% dphi - 1)*rotate)*dphi+phi){
for (i in 0:((360 %/% dtheta - 1)*rotate)*dtheta+theta){
if (plot.errors){
persp(x1.eval,
x2.eval,
lerr,
zlim = zlim,
cex.axis = ifelse(!is.null(cex.axis),cex.axis,par()\$cex.axis),
cex.lab = ifelse(!is.null(cex.lab),cex.lab,par()\$cex.lab),
cex.main = ifelse(!is.null(cex.main),cex.main,par()\$cex.main),
cex.sub = ifelse(!is.null(cex.sub),cex.sub,par()\$cex.sub),
col = persp.col,
border = ifelse(!is.null(border),border,"grey"),
ticktype = "detailed",
xlab = "",
ylab = "",
zlab = "",
theta = i,
phi = phi,
lwd = ifelse(!is.null(lwd),lwd,par()\$lwd))
par(new = TRUE)
}

persp(x1.eval,
x2.eval,
treg,
zlim = zlim,
col = persp.col,
border = ifelse(!is.null(border),border,"black"),
ticktype = "detailed",
cex.axis = ifelse(!is.null(cex.axis),cex.axis,par()\$cex.axis),
cex.lab =  ifelse(!is.null(cex.lab),cex.lab,par()\$cex.lab),
cex.main = ifelse(!is.null(cex.main),cex.main,par()\$cex.main),
cex.sub = ifelse(!is.null(cex.sub),cex.sub,par()\$cex.sub),
xlab = ifelse(!is.null(xlab),xlab,gen.label(bws\$xnames[1], "X1")),
ylab = ifelse(!is.null(ylab),ylab,gen.label(x2.names[1], "X2")),
zlab = ifelse(!is.null(zlab),zlab,gen.label(bws\$ynames,"Conditional Mean")),
theta = i,
phi = phi,
main = gen.tflabel(!is.null(main), main, paste("[theta= ", i,", phi= ", phi,"]", sep="")))

if (plot.errors){
par(new = TRUE)
persp(x1.eval,
x2.eval,
herr,
zlim = zlim,
cex.axis = ifelse(!is.null(cex.axis),cex.axis,par()\$cex.axis),
cex.lab = ifelse(!is.null(cex.lab),cex.lab,par()\$cex.lab),
cex.main = ifelse(!is.null(cex.main),cex.main,par()\$cex.main),
cex.sub = ifelse(!is.null(cex.sub),cex.sub,par()\$cex.sub),
col = persp.col,
border = ifelse(!is.null(border),border,"grey"),
ticktype = "detailed",
xlab = "",
ylab = "",
zlab = "",
theta = i,
phi = phi,
lwd = ifelse(!is.null(lwd),lwd,par()\$lwd))
}

Sys.sleep(0.5)
}
##}

if (plot.behavior == "plot-data")
return ( list(r1 = r1) )

} else {

tot.dim <- (bws\$xndim <- length(bws\$xdati\$icon)) + (bws\$zndim <- length(bws\$zdati\$icon))

if (plot.behavior != "data" && plot.par.mfrow)
par(mfrow=n2mfrow(tot.dim),cex=par()\$cex)

maxneval = max(c(sapply(xdat,nlevels), unlist(sapply(zdat,nlevels)), neval))
all.isFactor = c(sapply(xdat, is.factor), unlist(sapply(zdat, is.factor)))

x.ev = xdat[1,,drop = FALSE]

for (i in 1:bws\$xndim)
x.ev[1,i] = uocquantile(xdat[,i], prob=xq[i])

exdat = as.data.frame(matrix(data = 0, nrow = maxneval, ncol = bws\$xndim))

for (i in 1:bws\$xndim)
exdat[,i] = x.ev[1,i]

if (!miss.z){
z.ev = zdat[1,,drop = FALSE]

for (i in 1:bws\$zndim)
z.ev[1,i] = uocquantile(zdat[,i], prob=zq[i])

ezdat = as.data.frame(matrix(data = 0, nrow = maxneval, ncol = bws\$zndim))

for (i in 1:bws\$zndim)
ezdat[,i] = z.ev[1,i]

}

if (common.scale){
data.eval = matrix(data = NA, nrow = maxneval,
ncol = tot.dim)

data.err = matrix(data = NA, nrow = maxneval,
ncol = 3*tot.dim)

allei = as.data.frame(matrix(data = NA, nrow = maxneval,
ncol = tot.dim))

all.bxp = list()
}

plot.out = list()

temp.err = matrix(data = NA, nrow = maxneval, ncol = 3)
temp.mean = replicate(maxneval, NA)

## plotting expressions

plot.bootstrap = plot.errors.method == "bootstrap"

pfunE = expression(ifelse(xi.factor,
ifelse(plot.bootstrap & plot.bxp,"bxp","plotFactor"), "plot"))
pxE = expression(ifelse(common.scale,
ifelse(xi.factor,
ifelse(plot.bootstrap & plot.bxp, "all.bxp[[plot.index]],",
"f = allei[,plot.index],"),
"x = allei[,plot.index],"),
ifelse(xi.factor,
ifelse(plot.bootstrap & plot.bxp, "z = temp.boot,", "f = ei,"),
"x = ei,")))

pyE = expression(ifelse(xi.factor & plot.bootstrap & plot.bxp, "",
ifelse(common.scale,"y = data.eval[,plot.index],", "y = temp.mean,")))

pylimE = ifelse(common.scale, "ylim = c(y.min,y.max),",
ifelse(plot.errors, "ylim = c(min(na.omit(c(temp.mean - temp.err[,1], temp.err[,3] - temp.err[,1]))),
max(na.omit(c(temp.mean + temp.err[,2], temp.err[,3] + temp.err[,2])))),", ""))

pxlabE = expression(paste("xlab = gen.label(bws\$",
xOrZ, "names[i], paste('", toupper(xOrZ),"', i, sep = '')),",sep=''))

paste('Gradient Component ', i, ' of', sep=''), ''),
gen.label(bws\$ynames, 'Conditional Mean')),"

prestE = expression(ifelse(xi.factor,"", "type = ifelse(!is.null(type),type,'l'), lty = ifelse(!is.null(lty),lty,par()\$lty), col = ifelse(!is.null(col),col,par()\$col), lwd = ifelse(!is.null(lwd),lwd,par()\$lwd), cex.axis = ifelse(!is.null(cex.axis),cex.axis,par()\$cex.axis), cex.lab = ifelse(!is.null(cex.lab),cex.lab,par()\$cex.lab), cex.main = ifelse(!is.null(cex.main),cex.main,par()\$cex.main), cex.sub = ifelse(!is.null(cex.sub),cex.sub,par()\$cex.sub),"))

pmainE = "main = ifelse(!is.null(main),main,''), sub = ifelse(!is.null(sub),sub,''),"

txobjE <-
parse(text = paste("npscoef(txdat = xdat, tydat = ydat,",
ifelse(miss.z,"","tzdat = zdat,"),
"exdat = subcol(exdat,ei,i)[1:xi.neval,, drop = FALSE],",
ifelse(miss.z,"","ezdat = ezdat[1:xi.neval,, drop = FALSE],"),
"bws = bws, errors = plot.errors)"))

## error plotting expressions
plotOnEstimate = (plot.errors.center == "estimate")

efunE = "draw.errors"
eexE = expression(ifelse(common.scale, "ex = as.numeric(na.omit(allei[,plot.index])),",
"ex = as.numeric(na.omit(ei)),"))
eelyE = expression(ifelse(common.scale,
ifelse(plotOnEstimate, "ely = na.omit(data.eval[,plot.index] - data.err[,3*plot.index-2]),",
"ely = na.omit(data.err[,3*plot.index] - data.err[,3*plot.index-2]),"),
ifelse(plotOnEstimate, "ely = na.omit(temp.mean - temp.err[,1]),",
"ely = na.omit(temp.err[,3] - temp.err[,1]),")))
eehyE = expression(ifelse(common.scale,
ifelse(plotOnEstimate, "ehy = na.omit(data.eval[,plot.index] + data.err[,3*plot.index-1]),",
"ehy = na.omit(data.err[,3*plot.index] + data.err[,3*plot.index-1]),"),
ifelse(plotOnEstimate, "ehy = na.omit(temp.mean + temp.err[,2]),",
"ehy = na.omit(temp.err[,3] + temp.err[,2]),")))

erestE = "plot.errors.style = ifelse(xi.factor,'bar',plot.errors.style),
plot.errors.bar = ifelse(xi.factor,'I',plot.errors.bar),
plot.errors.bar.num = plot.errors.bar.num,
lty = ifelse(xi.factor,1,2)"

plot.index = 0
xOrZ = "x"

for (i in 1:bws\$xndim){
plot.index = plot.index + 1
temp.err[,] = NA
temp.mean[] =  NA
temp.boot = list()

xi.factor = all.isFactor[plot.index]

if (xi.factor){
ei = bws\$xdati\$all.ulev[[i]]
xi.neval = length(ei)
} else {
xi.neval = neval
qi = trim.quantiles(xdat[,i], xtrim[i])
ei = seq(qi[1], qi[2], length.out = neval)
}

if (xi.neval < maxneval){
ei[(xi.neval+1):maxneval] = NA
}

tobj <- eval(txobjE)

temp.mean[1:xi.neval] = tobj\$mean

if (plot.errors){
if (plot.errors.method == "asymptotic")
temp.err[1:xi.neval,1:2] = 2.0*tobj\$merr
else if (plot.errors.method == "bootstrap"){
temp.boot <- eval(parse(text = paste("compute.bootstrap.errors(",
"xdat = xdat, ydat = ydat,",
ifelse(miss.z, "", "zdat = zdat,"),
"exdat = subcol(exdat,ei,i)[1:xi.neval,, drop = FALSE],",
ifelse(miss.z,"","ezdat = ezdat[1:xi.neval,, drop = FALSE],"),
"slice.index = plot.index,",
"plot.errors.boot.method = plot.errors.boot.method,",
"plot.errors.boot.blocklen = plot.errors.boot.blocklen,",
"plot.errors.boot.num = plot.errors.boot.num,",
"plot.errors.center = plot.errors.center,",
"plot.errors.type = plot.errors.type,",
"plot.errors.quantiles = plot.errors.quantiles,",
"bws = bws)")))
temp.err[1:xi.neval,] <- temp.boot[["boot.err"]]
temp.boot <- temp.boot[["bxp"]]
if (!plot.bxp.out){
temp.boot\$out <- numeric()
temp.boot\$group <- integer()
}
}
}

if (common.scale){
allei[,plot.index] = ei
data.eval[, plot.index] = temp.mean
if (plot.errors){
all.bxp[plot.index] = NA
all.bxp[[plot.index]] = temp.boot

data.err[, c(3*plot.index-2,3*plot.index-1,3*plot.index)] = temp.err
}
} else if (plot.behavior != "data") {
## plot evaluation
eval(parse(text = paste(eval(pfunE), "(", eval(pxE), eval(pyE),
eval(pylimE), eval(pxlabE), eval(pylabE), eval(prestE),
eval(pmainE), ")")))

## error plotting evaluation
if (plot.errors && !(xi.factor & plot.bootstrap & plot.bxp)){
if (!xi.factor && !plotOnEstimate)
lines(na.omit(ei), na.omit(temp.err[,3]), lty = 3)

eval(parse(text = paste(eval(efunE), "(", eval(eexE), eval(eelyE),
eval(eehyE), eval(erestE), ")")))

}
}

if (plot.behavior != "plot") {
plot.out[plot.index] = NA
} else {
plot.out[[plot.index]] =
eval(parse(text = paste("smoothcoefficient(bws = bws,",
"eval = ", ifelse(miss.z, "subcol(exdat,ei,i)[1:xi.neval,, drop = FALSE],",
"list(exdat = subcol(exdat,ei,i)[1:xi.neval,, drop = FALSE], ezdat = ezdat[1:xi.neval,, drop = FALSE]),"),
"mean = na.omit(temp.mean),",
"ntrain = dim(xdat)[1],",
"trainiseval = FALSE,",
"xtra = c(0, 0, 0, 0, 0, 0))")))
plot.out[[plot.index]]\$merr = NA
plot.out[[plot.index]]\$bias = NA

if (plot.errors)
plot.out[[plot.index]]\$merr = temp.err[,1:2]

if (plot.errors.center == "bias-corrected")
plot.out[[plot.index]]\$bias = temp.err[,3] - temp.mean
plot.out[[plot.index]]\$bxp = temp.boot
}
}
}

if (!miss.z){
xOrZ = "z"
for (i in 1:bws\$zndim){
plot.index = plot.index + 1
temp.err[,] = NA
temp.mean[] =  NA
temp.boot = list()

xi.factor = all.isFactor[plot.index]

if (xi.factor){
ei = bws\$zdati\$all.ulev[[i]]
xi.neval = length(ei)
} else {
xi.neval = neval
qi = trim.quantiles(zdat[,i], ztrim[i])
ei = seq(qi[1], qi[2], length.out = neval)
}

if (xi.neval < maxneval){
ei[(xi.neval+1):maxneval] = NA
}

tobj <- npscoef(txdat = xdat, tydat = ydat, tzdat = zdat,
exdat = exdat[1:xi.neval,, drop = FALSE],
ezdat = subcol(ezdat,ei,i)[1:xi.neval,, drop = FALSE],
bws = bws)

temp.mean[1:xi.neval] = tobj\$mean

if (plot.errors){
if (plot.errors.method == "asymptotic")
temp.err[1:xi.neval,1:2] = 2.0*tobj\$merr
else if (plot.errors.method == "bootstrap"){
temp.boot <- compute.bootstrap.errors(
xdat = xdat,
ydat = ydat,
zdat = zdat,
exdat = exdat[1:xi.neval,, drop = FALSE],
ezdat = subcol(ezdat,ei,i)[1:xi.neval,, drop = FALSE],
slice.index = plot.index,
plot.errors.boot.method = plot.errors.boot.method,
plot.errors.boot.blocklen = plot.errors.boot.blocklen,
plot.errors.boot.num = plot.errors.boot.num,
plot.errors.center = plot.errors.center,
plot.errors.type = plot.errors.type,
plot.errors.quantiles = plot.errors.quantiles,
bws = bws)
temp.err[1:xi.neval,] <- temp.boot[["boot.err"]]
temp.boot <- temp.boot[["bxp"]]
if (!plot.bxp.out){
temp.boot\$out <- numeric()
temp.boot\$group <- integer()
}
}
}

if (common.scale){
allei[,plot.index] = ei
data.eval[, plot.index] = temp.mean
if (plot.errors){
all.bxp[plot.index] = NA
all.bxp[[plot.index]] = temp.boot

data.err[, c(3*plot.index-2,3*plot.index-1,3*plot.index)] = temp.err
}
} else if (plot.behavior != "data") {
## plot evaluation
eval(parse(text = paste(eval(pfunE), "(", eval(pxE), eval(pyE),
eval(pylimE), eval(pxlabE), eval(pylabE), eval(prestE),
eval(pmainE), ")")))

## error plotting evaluation
if (plot.errors && !(xi.factor & plot.bootstrap & plot.bxp)){
if (!xi.factor && !plotOnEstimate)
lines(na.omit(ei), na.omit(temp.err[,3]), lty = 3)

eval(parse(text = paste(eval(efunE), "(", eval(eexE), eval(eelyE),
eval(eehyE), eval(erestE), ")")))

}
}

if (plot.behavior != "plot") {
plot.out[plot.index] = NA
} else {
plot.out[[plot.index]] =
smoothcoefficient(bws = bws,
eval = list(exdat = exdat[1:xi.neval,, drop = FALSE],
ezdat = subcol(ezdat,ei,i)[1:xi.neval,, drop = FALSE]),
mean = na.omit(temp.mean),
ntrain = dim(zdat)[1],
trainiseval = FALSE)

plot.out[[plot.index]]\$merr = NA
plot.out[[plot.index]]\$bias = NA

if (plot.errors)
plot.out[[plot.index]]\$merr = temp.err[,1:2]

if (plot.errors.center == "bias-corrected")
plot.out[[plot.index]]\$bias = temp.err[,3] - temp.mean
plot.out[[plot.index]]\$bxp = temp.boot
}
}
}
}

if (common.scale & (plot.behavior != "data")){
jj = 1:(bws\$xndim + bws\$zndim)*3

if (plot.errors.center == "estimate" | !plot.errors) {
y.max = max(na.omit(as.double(data.eval)) +
if (plot.errors) na.omit(as.double(data.err[,jj-1]))
else 0)
y.min = min(na.omit(as.double(data.eval)) -
if (plot.errors) na.omit(as.double(data.err[,jj-2]))
else 0)
} else if (plot.errors.center == "bias-corrected") {
y.max = max(na.omit(as.double(data.err[,jj] + data.err[,jj-1])))
y.min = min(na.omit(as.double(data.err[,jj] - data.err[,jj-2])))
}

if(!is.null(ylim)){
y.min = ylim[1]
y.max = ylim[2]
}

xOrZ = "x"

for (plot.index in 1:(bws\$xndim + bws\$zndim)){
i = ifelse(plot.index <= bws\$xndim, plot.index, plot.index - bws\$xndim)

if (plot.index > bws\$xndim)
xOrZ <- "z"

xi.factor = all.isFactor[plot.index]

## plot evaluation
eval(parse(text = paste(eval(pfunE), "(", eval(pxE), eval(pyE),
eval(pylimE), eval(pxlabE), eval(pylabE), eval(prestE),
eval(pmainE), ")")))

## error plotting evaluation
if (plot.errors && !(xi.factor & plot.bootstrap & plot.bxp)){
if (!xi.factor && !plotOnEstimate)
lines(na.omit(ei), na.omit(temp.err[,3]), lty = 3)

eval(parse(text = paste(eval(efunE), "(", eval(eexE), eval(eelyE),
eval(eehyE), eval(erestE), ")")))

}
}
}

if (plot.behavior != "data" && plot.par.mfrow)
par(mfrow=c(1,1),cex=par()\$cex)

if (plot.behavior != "plot"){
names(plot.out) =
else
paste("sc",1:(bws\$xndim+bws\$zndim),sep="")

return (plot.out)
}
}

}

npplot.plbandwidth <-
function(bws,
xdat,
ydat,
zdat,
data = NULL,
xq = 0.5,
zq = 0.5,
xtrim = 0.0,
ztrim = 0.0,
neval = 50,
common.scale = TRUE,
perspective = TRUE,
main = NULL,
type = NULL,
border = NULL,
cex.axis = NULL,
cex.lab = NULL,
cex.main = NULL,
cex.sub = NULL,
col = NULL,
ylab = NULL,
xlab = NULL,
zlab = NULL,
sub = NULL,
ylim = NULL,
xlim = NULL,
zlim = NULL,
lty = NULL,
lwd = NULL,
theta = 0.0,
phi = 10.0,
view = c("rotate","fixed"),
plot.behavior = c("plot","plot-data","data"),
plot.errors.method = c("none","bootstrap","asymptotic"),
plot.errors.boot.method = c("inid", "fixed", "geom"),
plot.errors.boot.blocklen = NULL,
plot.errors.boot.num = 399,
plot.errors.center = c("estimate","bias-corrected"),
plot.errors.type = c("standard","quantiles"),
plot.errors.quantiles = c(0.025,0.975),
plot.errors.style = c("band","bar"),
plot.errors.bar = c("|","I"),
plot.errors.bar.num = min(neval,25),
plot.bxp = FALSE,
plot.bxp.out = TRUE,
plot.par.mfrow = TRUE,
...,
random.seed){

if(!is.null(options('plot.par.mfrow')\$plot.par.mfrow))
plot.par.mfrow <- options('plot.par.mfrow')\$plot.par.mfrow

stop("gradients not supported with partially linear models. Coefficients may be extracted with coef()")

miss.xyz = c(missing(xdat), missing(ydat), missing(zdat))

if (any(miss.xyz) && !all(miss.xyz))
stop("one of, but not both, xdat and ydat was specified")
else if(all(miss.xyz) && !is.null(bws\$formula)){
tt <- terms(bws)
m <- match(c("formula", "data", "subset", "na.action"),
names(bws\$call), nomatch = 0)

tmf.xf <- tmf.x <- tmf <- bws\$call[c(1,m)]
tmf.x[[1]] <- as.name("model.matrix")
tmf.xf[[1]] <- tmf[[1]] <- as.name("model.frame")
tmf[["formula"]] <- tt
umf <- tmf <- eval(tmf, envir = environment(tt))

bronze <- lapply(bws\$chromoly, paste, collapse = " + ")

tmf.x[["object"]] <- as.formula(paste(" ~ ", bronze[[2]]),
env = environment(formula))
tmf.x <- eval(tmf.x,parent.frame())

tmf.xf[["formula"]] <- as.formula(paste(" ~ ", bronze[[2]]),
env = environment(formula))
tmf.xf <- eval(tmf.xf,parent.frame())

ydat <- model.response(tmf)
xdat <- as.data.frame(tmf.x[,-1, drop = FALSE])

cc <- attr(tmf.x,'assign')[-1]

for(i in 1:length(cc))
xdat[,i] <- cast(xdat[,i], tmf.xf[,cc[i]], same.levels = FALSE)

zdat <- tmf[, bws\$chromoly[[3]], drop = FALSE]
} else {
if(all(miss.xyz) && !is.null(bws\$call)){
xdat <- data.frame(eval(bws\$call[["xdat"]], environment(bws\$call)))
ydat = eval(bws\$call[["ydat"]], environment(bws\$call))
zdat <- data.frame(eval(bws\$call[["zdat"]], environment(bws\$call)))
}
xdat = toFrame(xdat)
zdat = toFrame(zdat)

## catch and destroy NA's
goodrows = 1:dim(xdat)[1]
rows.omit = attr(na.omit(data.frame(xdat,ydat,zdat)), "na.action")
goodrows[rows.omit] = 0

if (all(goodrows==0))
stop("Training data has no rows without NAs")

xdat = xdat[goodrows,,drop = FALSE]
ydat = ydat[goodrows]
zdat = zdat[goodrows,,drop = FALSE]
}

nxcon = sum(bws\$xdati\$icon)
nxuno = sum(bws\$xdati\$iuno)
nxord = sum(bws\$xdati\$iord)

nzcon = sum(bws\$zdati\$icon)
nzuno = sum(bws\$zdati\$iuno)
nzord = sum(bws\$zdati\$iord)

xq = double(bws\$xndim)+xq
zq = double(bws\$zndim)+zq

xtrim = double(bws\$xndim)+xtrim
ztrim = double(bws\$zndim)+ztrim

if (missing(plot.errors.method) &
any(!missing(plot.errors.boot.num), !missing(plot.errors.boot.method),
!missing(plot.errors.boot.blocklen))){
warning(paste("plot.errors.method must be set to 'bootstrap' to use bootstrapping.",
"\nProceeding without bootstrapping."))
}

plot.behavior = match.arg(plot.behavior)
plot.errors.method = match.arg(plot.errors.method)
plot.errors.boot.method = match.arg(plot.errors.boot.method)
plot.errors.center = match.arg(plot.errors.center)
plot.errors.type = match.arg(plot.errors.type)
plot.errors.style = match.arg(plot.errors.style)
plot.errors.bar = match.arg(plot.errors.bar)

common.scale = common.scale | (!is.null(ylim))

if (plot.errors.method == "asymptotic") {
warning(paste("asymptotic errors are not supported with partially linear regression.\n",
"Proceeding without calculating errors"))
plot.errors.method = "none"
}

if (is.element(plot.errors.boot.method, c("fixed", "geom")) &&
is.null(plot.errors.boot.blocklen))
plot.errors.boot.blocklen = b.star(xdat,round=TRUE)[1,1]

plot.errors = (plot.errors.method != "none")

if ((nxcon + nxord == 1) & (nzcon + nzord == 1) & (nxuno + nzuno == 0) &
perspective & !gradients & !any(xor(bws\$xdati\$iord, bws\$xdati\$inumord)) &
!any(xor(bws\$zdati\$iord, bws\$zdati\$inumord))){

view = match.arg(view)
rotate = (view == "rotate")

if (is.ordered(xdat[,1])){
x1.eval = bws\$xdati\$all.ulev[[1]]
x1.neval = length(x1.eval)
} else {
x1.neval = neval
qi = trim.quantiles(xdat[,1], xtrim[1])
x1.eval = seq(qi[1], qi[2], length.out = x1.neval)
}

if (is.ordered(zdat[,1])){
z1.eval = bws\$zdati\$all.ulev[[1]]
z1.neval = length(z1.eval)
} else {
z1.neval = neval
qi = trim.quantiles(zdat[,1], ztrim[1])
z1.eval = seq(qi[1], qi[2], length.out = z1.neval)
}

x.eval <- expand.grid(x1.eval, z1.eval)

if (bws\$xdati\$iord[1])
x1.eval <- (bws\$xdati\$all.dlev[[1]])[as.integer(x1.eval)]

if (bws\$zdati\$iord[1])
z1.eval <- (bws\$zdati\$all.dlev[[1]])[as.integer(z1.eval)]

tobj = npplreg(txdat = xdat, tydat = ydat, tzdat = zdat,
exdat = x.eval[,1], ezdat = x.eval[,2], bws = bws)

terr = matrix(data = NA, nrow = nrow(x.eval), ncol = 3)

treg = matrix(data = tobj\$mean,
nrow = x1.neval, ncol = z1.neval, byrow = FALSE)

if (plot.errors.method == "bootstrap"){
terr <- compute.bootstrap.errors(
xdat = xdat, ydat = ydat, zdat = zdat,
exdat = x.eval[,1], ezdat = x.eval[,2],
slice.index = 0,
plot.errors.boot.method = plot.errors.boot.method,
plot.errors.boot.blocklen = plot.errors.boot.blocklen,
plot.errors.boot.num = plot.errors.boot.num,
plot.errors.center = plot.errors.center,
plot.errors.type = plot.errors.type,
plot.errors.quantiles = plot.errors.quantiles,
bws = bws)[["boot.err"]]

pc = (plot.errors.center == "bias-corrected")

lerr = matrix(data = if(pc) {terr[,3]} else {tobj\$mean}
-terr[,1],
nrow = x1.neval, ncol = z1.neval, byrow = FALSE)

herr = matrix(data = if(pc) {terr[,3]} else {tobj\$mean}
+terr[,2],
nrow = x1.neval, ncol = z1.neval, byrow = FALSE)

}

if(is.null(zlim)) {
zlim =
if (plot.errors)
c(min(lerr),max(herr))
else
c(min(tobj\$mean),max(tobj\$mean))
}

if (plot.behavior != "plot"){
r1 = plregression(bws = bws, xcoef = tobj\$xcoef,
xcoefvcov = vcov(tobj),
xcoeferr = tobj\$xcoeferr,
evalx =  x.eval[,1],
evalz =  x.eval[,2],
mean = tobj\$mean,
ntrain = dim(xdat)[1],
trainiseval = FALSE,
xtra=c(tobj\$RSQ,tobj\$MSE,0,0,0,0))

r1\$merr = NA
r1\$bias = NA

if (plot.errors)
r1\$merr = terr[,1:2]

if (plot.errors.center == "bias-corrected")
r1\$bias = terr[,3] - treg

if (plot.behavior == "data")
return ( list(r1 = r1) )

}

dtheta = 5.0
dphi = 10.0

persp.col = ifelse(plot.errors, FALSE, ifelse(!is.null(col),col,"lightblue"))

##for (j in 0:((50 %/% dphi - 1)*rotate)*dphi+phi){
for (i in 0:((360 %/% dtheta - 1)*rotate)*dtheta+theta){
if (plot.errors){
persp(x1.eval,
z1.eval,
lerr,
zlim = zlim,
cex.axis = ifelse(!is.null(cex.axis),cex.axis,par()\$cex.axis),
cex.lab = ifelse(!is.null(cex.lab),cex.lab,par()\$cex.lab),
cex.main = ifelse(!is.null(cex.main),cex.main,par()\$cex.main),
cex.sub = ifelse(!is.null(cex.sub),cex.sub,par()\$cex.sub),
col = persp.col,
border = ifelse(!is.null(border),border,"grey"),
ticktype = "detailed",
xlab = "",
ylab = "",
zlab = "",
theta = i,
phi = phi,
lwd = ifelse(!is.null(lwd),lwd,par()\$lwd))
par(new = TRUE)
}

persp(x1.eval,
z1.eval,
treg,
zlim = zlim,
col = persp.col,
border = ifelse(!is.null(border),border,"black"),
ticktype = "detailed",
cex.axis = ifelse(!is.null(cex.axis),cex.axis,par()\$cex.axis),
cex.lab =  ifelse(!is.null(cex.lab),cex.lab,par()\$cex.lab),
cex.main = ifelse(!is.null(cex.main),cex.main,par()\$cex.main),
cex.sub = ifelse(!is.null(cex.sub),cex.sub,par()\$cex.sub),
xlab = ifelse(!is.null(xlab),xlab,gen.label(names(xdat)[1], "X1")),
ylab = ifelse(!is.null(ylab),ylab,gen.label(names(xdat)[2], "Z1")),
zlab = ifelse(!is.null(zlab),zlab,gen.label(names(ydat),"Conditional Mean")),
theta = i,
phi = phi,
main = gen.tflabel(!is.null(main), main, paste("[theta= ", i,", phi= ", phi,"]", sep="")))

if (plot.errors){
par(new = TRUE)
persp(x1.eval,
z1.eval,
herr,
zlim = zlim,
cex.axis = ifelse(!is.null(cex.axis),cex.axis,par()\$cex.axis),
cex.lab = ifelse(!is.null(cex.lab),cex.lab,par()\$cex.lab),
cex.main = ifelse(!is.null(cex.main),cex.main,par()\$cex.main),
cex.sub = ifelse(!is.null(cex.sub),cex.sub,par()\$cex.sub),
col = persp.col,
border = ifelse(!is.null(border),border,"grey"),
ticktype = "detailed",
xlab = "",
ylab = "",
zlab = "",
theta = i,
phi = phi,
lwd = ifelse(!is.null(lwd),lwd,par()\$lwd))
}

Sys.sleep(0.5)
}
##}

if (plot.behavior == "plot-data")
return ( list(r1 = r1) )

} else {
##      stop("not yet supported!")

if (plot.behavior != "data" && plot.par.mfrow)
par(mfrow=n2mfrow(bws\$xndim + bws\$zndim),cex=par()\$cex)

x.ev = xdat[1,,drop = FALSE]
z.ev = zdat[1,,drop = FALSE]

for (i in 1:bws\$xndim)
x.ev[1,i] = uocquantile(xdat[,i], prob=xq[i])

for (i in 1:bws\$zndim)
z.ev[1,i] = uocquantile(zdat[,i], prob=zq[i])

maxneval = max(c(sapply(xdat,nlevels), sapply(zdat,nlevels), neval))

exdat = as.data.frame(matrix(data = 0, nrow = maxneval, ncol = bws\$xndim))
ezdat = as.data.frame(matrix(data = 0, nrow = maxneval, ncol = bws\$zndim))

for (i in 1:bws\$xndim)
exdat[,i] = x.ev[1,i]

for (i in 1:bws\$zndim)
ezdat[,i] = z.ev[1,i]

if (common.scale){
data.eval = matrix(data = NA, nrow = maxneval,
ncol = (bws\$xndim + bws\$zndim))

data.err = matrix(data = NA, nrow = maxneval,
ncol = 3*(bws\$xndim + bws\$zndim))

allei = as.data.frame(matrix(data = NA, nrow = maxneval,
ncol = bws\$xndim + bws\$zndim))

all.bxp = list()
}

all.isFactor = c(sapply(xdat, is.factor), sapply(zdat, is.factor))

plot.out = list()

temp.err = matrix(data = NA, nrow = maxneval, ncol = 3)
temp.mean = replicate(maxneval, NA)

## plotting expressions

plot.bootstrap = plot.errors.method == "bootstrap"

pfunE = expression(ifelse(xi.factor,
ifelse(plot.bootstrap & plot.bxp,"bxp","plotFactor"), "plot"))
pxE = expression(ifelse(common.scale,
ifelse(xi.factor,
ifelse(plot.bootstrap & plot.bxp, "all.bxp[[plot.index]],",
"f = allei[,plot.index],"),
"x = allei[,plot.index],"),
ifelse(xi.factor,
ifelse(plot.bootstrap & plot.bxp, "z = temp.boot,", "f = ei,"),
"x = ei,")))

pyE = expression(ifelse(xi.factor & plot.bootstrap & plot.bxp, "",
ifelse(common.scale,"y = data.eval[,plot.index],", "y = temp.mean,")))

pylimE = ifelse(common.scale, "ylim = c(y.min,y.max),",
ifelse(plot.errors, "ylim = c(min(na.omit(c(temp.mean - temp.err[,1], temp.err[,3] - temp.err[,1]))),
max(na.omit(c(temp.mean + temp.err[,2], temp.err[,3] + temp.err[,2])))),", ""))

pxlabE = expression(paste("xlab = gen.label(bws\$",
xOrZ, "names[i], paste('", toupper(xOrZ),"', i, sep = '')),",sep=''))

paste('Gradient Component ', i, ' of', sep=''), ''),
gen.label(bws\$ynames, 'Conditional Mean')),"

prestE = expression(ifelse(xi.factor,"", "type = ifelse(!is.null(type),type,'l'), lty = ifelse(!is.null(lty),lty,par()\$lty), col = ifelse(!is.null(col),col,par()\$col), lwd = ifelse(!is.null(lwd),lwd,par()\$lwd), cex.axis = ifelse(!is.null(cex.axis),cex.axis,par()\$cex.axis), cex.lab = ifelse(!is.null(cex.lab),cex.lab,par()\$cex.lab), cex.main = ifelse(!is.null(cex.main),cex.main,par()\$cex.main), cex.sub = ifelse(!is.null(cex.sub),cex.sub,par()\$cex.sub),"))

pmainE = "main = ifelse(!is.null(main),main,''), sub = ifelse(!is.null(sub),sub,''),"

## error plotting expressions
plotOnEstimate = (plot.errors.center == "estimate")

efunE = "draw.errors"
eexE = expression(ifelse(common.scale, "ex = as.numeric(na.omit(allei[,plot.index])),",
"ex = as.numeric(na.omit(ei)),"))
eelyE = expression(ifelse(common.scale,
ifelse(plotOnEstimate, "ely = na.omit(data.eval[,plot.index] - data.err[,3*plot.index-2]),",
"ely = na.omit(data.err[,3*plot.index] - data.err[,3*plot.index-2]),"),
ifelse(plotOnEstimate, "ely = na.omit(temp.mean - temp.err[,1]),",
"ely = na.omit(temp.err[,3] - temp.err[,1]),")))
eehyE = expression(ifelse(common.scale,
ifelse(plotOnEstimate, "ehy = na.omit(data.eval[,plot.index] + data.err[,3*plot.index-1]),",
"ehy = na.omit(data.err[,3*plot.index] + data.err[,3*plot.index-1]),"),
ifelse(plotOnEstimate, "ehy = na.omit(temp.mean + temp.err[,2]),",
"ehy = na.omit(temp.err[,3] + temp.err[,2]),")))

erestE = "plot.errors.style = ifelse(xi.factor,'bar',plot.errors.style),
plot.errors.bar = ifelse(xi.factor,'I',plot.errors.bar),
plot.errors.bar.num = plot.errors.bar.num,
lty = ifelse(xi.factor,1,2)"

plot.index = 0
xOrZ = "x"

for (i in 1:bws\$xndim){
plot.index = plot.index + 1
temp.err[,] = NA
temp.mean[] =  NA
temp.boot = list()

xi.factor = all.isFactor[plot.index]

if (xi.factor){
ei = bws\$xdati\$all.ulev[[i]]
xi.neval = length(ei)
} else {
xi.neval = neval
qi = trim.quantiles(xdat[,i], xtrim[i])
ei = seq(qi[1], qi[2], length.out = neval)
}

if (xi.neval < maxneval){
ei[(xi.neval+1):maxneval] = NA
}

tobj <- npplreg(txdat = xdat, tydat = ydat, tzdat = zdat,
exdat = subcol(exdat,ei,i)[1:xi.neval,, drop = FALSE],
ezdat = ezdat[1:xi.neval,, drop = FALSE],
bws = bws)

temp.mean[1:xi.neval] = tobj\$mean

if (plot.errors){
if (plot.errors.method == "bootstrap"){
temp.boot <- compute.bootstrap.errors(
xdat = xdat,
ydat = ydat,
zdat = zdat,
exdat = subcol(exdat,ei,i)[1:xi.neval,, drop = FALSE],
ezdat = ezdat[1:xi.neval,, drop = FALSE],
slice.index = plot.index,
plot.errors.boot.method = plot.errors.boot.method,
plot.errors.boot.blocklen = plot.errors.boot.blocklen,
plot.errors.boot.num = plot.errors.boot.num,
plot.errors.center = plot.errors.center,
plot.errors.type = plot.errors.type,
plot.errors.quantiles = plot.errors.quantiles,
bws = bws)
temp.err[1:xi.neval,] <- temp.boot[["boot.err"]]
temp.boot <- temp.boot[["bxp"]]
if (!plot.bxp.out){
temp.boot\$out <- numeric()
temp.boot\$group <- integer()
}
}
}

if (common.scale){
allei[,plot.index] = ei
data.eval[, plot.index] = temp.mean
if (plot.errors){
all.bxp[plot.index] = NA
all.bxp[[plot.index]] = temp.boot

data.err[, c(3*plot.index-2,3*plot.index-1,3*plot.index)] = temp.err
}
} else if (plot.behavior != "data") {
## plot evaluation
eval(parse(text = paste(eval(pfunE), "(", eval(pxE), eval(pyE),
eval(pylimE), eval(pxlabE), eval(pylabE), eval(prestE),
eval(pmainE), ")")))

## error plotting evaluation
if (plot.errors && !(xi.factor & plot.bootstrap & plot.bxp)){
if (!xi.factor && !plotOnEstimate)
lines(na.omit(ei), na.omit(temp.err[,3]), lty = 3)

eval(parse(text = paste(eval(efunE), "(", eval(eexE), eval(eelyE),
eval(eehyE), eval(erestE), ")")))

}
}

if (plot.behavior != "plot") {
plot.out[plot.index] = NA
} else {
plot.out[[plot.index]] =
plregression(bws = bws, xcoef = tobj\$xcoef, xcoefvcov = vcov(tobj),
xcoeferr = tobj\$xcoeferr,
evalx = subcol(exdat,ei,i)[1:xi.neval,, drop = FALSE],
evalz = ezdat[1:xi.neval,, drop = FALSE],
mean = na.omit(temp.mean),
ntrain = dim(xdat)[1],
trainiseval = FALSE,
xtra = c(tobj\$RSQ, tobj\$MSE, 0, 0, 0, 0))
plot.out[[plot.index]]\$merr = NA
plot.out[[plot.index]]\$bias = NA

if (plot.errors)
plot.out[[plot.index]]\$merr = temp.err[,1:2]

if (plot.errors.center == "bias-corrected")
plot.out[[plot.index]]\$bias = temp.err[,3] - temp.mean
plot.out[[plot.index]]\$bxp = temp.boot
}
}
}

xOrZ = "z"
for (i in 1:bws\$zndim){
plot.index = plot.index + 1
temp.err[,] = NA
temp.mean[] =  NA
temp.boot = list()

xi.factor = all.isFactor[plot.index]

if (xi.factor){
ei = bws\$zdati\$all.ulev[[i]]
xi.neval = length(ei)
} else {
xi.neval = neval
qi = trim.quantiles(zdat[,i], ztrim[i])
ei = seq(qi[1], qi[2], length.out = neval)
}

if (xi.neval < maxneval){
ei[(xi.neval+1):maxneval] = NA
}

tobj <- npplreg(txdat = xdat, tydat = ydat, tzdat = zdat,
exdat = exdat[1:xi.neval,, drop = FALSE],
ezdat = subcol(ezdat,ei,i)[1:xi.neval,, drop = FALSE],
bws = bws)

temp.mean[1:xi.neval] = tobj\$mean

if (plot.errors){
if (plot.errors.method == "bootstrap"){
temp.boot <- compute.bootstrap.errors(
xdat = xdat,
ydat = ydat,
zdat = zdat,
exdat = exdat[1:xi.neval,, drop = FALSE],
ezdat = subcol(ezdat,ei,i)[1:xi.neval,, drop = FALSE],
slice.index = plot.index,
plot.errors.boot.method = plot.errors.boot.method,
plot.errors.boot.blocklen = plot.errors.boot.blocklen,
plot.errors.boot.num = plot.errors.boot.num,
plot.errors.center = plot.errors.center,
plot.errors.type = plot.errors.type,
plot.errors.quantiles = plot.errors.quantiles,
bws = bws)
temp.err[1:xi.neval,] <- temp.boot[["boot.err"]]
temp.boot <- temp.boot[["bxp"]]
if (!plot.bxp.out){
temp.boot\$out <- numeric()
temp.boot\$group <- integer()
}
}
}

if (common.scale){
allei[,plot.index] = ei
data.eval[, plot.index] = temp.mean
if (plot.errors){
all.bxp[plot.index] = NA
all.bxp[[plot.index]] = temp.boot

data.err[, c(3*plot.index-2,3*plot.index-1,3*plot.index)] = temp.err
}
} else if (plot.behavior != "data") {
## plot evaluation
eval(parse(text = paste(eval(pfunE), "(", eval(pxE), eval(pyE),
eval(pylimE), eval(pxlabE), eval(pylabE), eval(prestE),
eval(pmainE), ")")))

## error plotting evaluation
if (plot.errors && !(xi.factor & plot.bootstrap & plot.bxp)){
if (!xi.factor && !plotOnEstimate)
lines(na.omit(ei), na.omit(temp.err[,3]), lty = 3)

eval(parse(text = paste(eval(efunE), "(", eval(eexE), eval(eelyE),
eval(eehyE), eval(erestE), ")")))

}
}

if (plot.behavior != "plot") {
plot.out[plot.index] = NA
} else {
plot.out[[plot.index]] =
plregression(bws = bws, xcoef = tobj\$xcoef,
xcoeferr = tobj\$xcoeferr,
xcoefvcov = vcov(tobj),
evalx = exdat[1:xi.neval,, drop = FALSE],
evalz = subcol(ezdat,ei,i)[1:xi.neval,, drop = FALSE],
mean = na.omit(temp.mean),
ntrain = dim(zdat)[1],
trainiseval = FALSE,
xtra = c(tobj\$RSQ, tobj\$MSE, 0, 0, 0, 0))
plot.out[[plot.index]]\$merr = NA
plot.out[[plot.index]]\$bias = NA

if (plot.errors)
plot.out[[plot.index]]\$merr = temp.err[,1:2]

if (plot.errors.center == "bias-corrected")
plot.out[[plot.index]]\$bias = temp.err[,3] - temp.mean
plot.out[[plot.index]]\$bxp = temp.boot
}
}
}

if (common.scale & (plot.behavior != "data")){
jj = 1:(bws\$xndim + bws\$zndim)*3

if (plot.errors.center == "estimate" | !plot.errors) {
y.max = max(na.omit(as.double(data.eval)) +
if (plot.errors) na.omit(as.double(data.err[,jj-1]))
else 0)
y.min = min(na.omit(as.double(data.eval)) -
if (plot.errors) na.omit(as.double(data.err[,jj-2]))
else 0)
} else if (plot.errors.center == "bias-corrected") {
y.max = max(na.omit(as.double(data.err[,jj] + data.err[,jj-1])))
y.min = min(na.omit(as.double(data.err[,jj] - data.err[,jj-2])))
}

if(!is.null(ylim)){
y.min = ylim[1]
y.max = ylim[2]
}

xOrZ = "x"

for (plot.index in 1:(bws\$xndim + bws\$zndim)){
i = ifelse(plot.index <= bws\$xndim, plot.index, plot.index - bws\$xndim)

if (plot.index > bws\$xndim)
xOrZ <- "z"

xi.factor = all.isFactor[plot.index]

## plot evaluation
eval(parse(text = paste(eval(pfunE), "(", eval(pxE), eval(pyE),
eval(pylimE), eval(pxlabE), eval(pylabE), eval(prestE),
eval(pmainE), ")")))

## error plotting evaluation
if (plot.errors && !(xi.factor & plot.bootstrap & plot.bxp)){
if (!xi.factor && !plotOnEstimate)
lines(na.omit(ei), na.omit(temp.err[,3]), lty = 3)

eval(parse(text = paste(eval(efunE), "(", eval(eexE), eval(eelyE),
eval(eehyE), eval(erestE), ")")))

}
}
}

if (plot.behavior != "data" && plot.par.mfrow)
par(mfrow=c(1,1),cex=par()\$cex)

if (plot.behavior != "plot"){
names(plot.out) =
else
paste("plr",1:(bws\$xndim+bws\$zndim),sep="")

return (plot.out)
}
}
}

npplot.bandwidth <-
function(bws,
xdat,
data = NULL,
xq = 0.5,
xtrim = 0.0,
neval = 50,
common.scale = TRUE,
perspective = TRUE,
main = NULL,
type = NULL,
border = NULL,
cex.axis = NULL,
cex.lab = NULL,
cex.main = NULL,
cex.sub = NULL,
col = NULL,
ylab = NULL,
xlab = NULL,
zlab = NULL,
sub = NULL,
ylim = NULL,
xlim = NULL,
zlim = NULL,
lty = NULL,
lwd = NULL,
theta = 0.0,
phi = 10.0,
view = c("rotate","fixed"),
plot.behavior = c("plot","plot-data","data"),
plot.errors.method = c("none","bootstrap","asymptotic"),
plot.errors.boot.method = c("inid", "fixed", "geom"),
plot.errors.boot.blocklen = NULL,
plot.errors.boot.num = 399,
plot.errors.center = c("estimate","bias-corrected"),
plot.errors.type = c("standard","quantiles"),
plot.errors.quantiles = c(0.025,0.975),
plot.errors.style = c("band","bar"),
plot.errors.bar = c("|","I"),
plot.errors.bar.num = min(neval,25),
plot.bxp = FALSE,
plot.bxp.out = TRUE,
plot.par.mfrow = TRUE,
...,
random.seed){

if(!is.null(options('plot.par.mfrow')\$plot.par.mfrow))
plot.par.mfrow <- options('plot.par.mfrow')\$plot.par.mfrow

miss.x <- missing(xdat)

if(miss.x && !is.null(bws\$formula)){
tt <- terms(bws)
m <- match(c("formula", "data", "subset", "na.action"),
names(bws\$call), nomatch = 0)
tmf <- bws\$call[c(1,m)]
tmf[[1]] <- as.name("model.frame")
tmf[["formula"]] <- tt
umf <- tmf <- eval(tmf, envir = environment(tt))

xdat <- tmf[, attr(attr(tmf, "terms"),"term.labels"), drop = FALSE]
} else {
if(miss.x && !is.null(bws\$call)){
xdat <- data.frame(eval(bws\$call[["dat"]], environment(bws\$call)))
}
xdat = toFrame(xdat)
xdat = na.omit(xdat)
}

xq = double(bws\$ndim)+xq
xtrim = double(bws\$ndim)+xtrim

if (missing(plot.errors.method) &
any(!missing(plot.errors.boot.num), !missing(plot.errors.boot.method),
!missing(plot.errors.boot.blocklen))){
warning(paste("plot.errors.method must be set to 'bootstrap' to use bootstrapping.",
"\nProceeding without bootstrapping."))
}

plot.behavior = match.arg(plot.behavior)
plot.errors.method = match.arg(plot.errors.method)
plot.errors.boot.method = match.arg(plot.errors.boot.method)
plot.errors.center = match.arg(plot.errors.center)
plot.errors.type = match.arg(plot.errors.type)
plot.errors.style = match.arg(plot.errors.style)
plot.errors.bar = match.arg(plot.errors.bar)

common.scale = common.scale | (!is.null(ylim))

if (plot.errors.method == "asymptotic") {
if (plot.errors.type == "quantiles"){
warning("quantiles cannot be calculated with asymptotics, calculating standard errors")
plot.errors.type = "standard"
}

if (plot.errors.center == "bias-corrected") {
warning("no bias corrections can be calculated with asymptotics, centering on estimate")
plot.errors.center = "estimate"
}
}

if (is.element(plot.errors.boot.method, c("fixed", "geom")) &&
is.null(plot.errors.boot.blocklen))
plot.errors.boot.blocklen = b.star(xdat,round=TRUE)[1,1]

plot.errors = (plot.errors.method != "none")

if ((bws\$ncon + bws\$nord == 2) & (bws\$nuno == 0) & perspective &
!any(xor(bws\$xdati\$iord, bws\$xdati\$inumord))){
view = match.arg(view)
rotate = (view == "rotate")

if (is.ordered(xdat[,1])){
x1.eval = bws\$xdati\$all.ulev[[1]]
x1.neval = length(x1.eval)
} else {
x1.neval = neval
qi = trim.quantiles(xdat[,1], xtrim[1])
x1.eval = seq(qi[1], qi[2], length.out = x1.neval)
}

if (is.ordered(xdat[,2])){
x2.eval = bws\$xdati\$all.ulev[[2]]
x2.neval = length(x2.eval)
} else {
x2.neval = neval
qi = trim.quantiles(xdat[,2], xtrim[2])
x2.eval = seq(qi[1], qi[2], length.out = x2.neval)
}

x.eval <- expand.grid(x1.eval, x2.eval)

##x.eval = as.data.frame(matrix(data = NA, ncol = 2, nrow = x1.neval*x2.neval))
##x.eval[,1] = rep(x1.eval, times = x2.neval)
##x.eval[,2] = rep(x2.eval, each = x1.neval)

if (is.ordered(xdat[,1]))
x1.eval <- (bws\$xdati\$all.dlev[[1]])[as.integer(x1.eval)]

if (is.ordered(xdat[,2]))
x2.eval <- (bws\$xdati\$all.dlev[[2]])[as.integer(x2.eval)]

tobj =  npudens(tdat = xdat, edat = x.eval, bws = bws)

tcomp = parse(text="tobj\$dens")

tdens = matrix(data = eval(tcomp),
nrow = x1.neval, ncol = x2.neval, byrow = FALSE)

terr = matrix(data = tobj\$derr, nrow = nrow(x.eval), ncol = 3)
terr[,3] = NA

if (plot.errors.method == "bootstrap"){
terr <- compute.bootstrap.errors(xdat = xdat,
exdat = x.eval,
cdf = FALSE,
slice.index = 0,
plot.errors.boot.method = plot.errors.boot.method,
plot.errors.boot.blocklen = plot.errors.boot.blocklen,
plot.errors.boot.num = plot.errors.boot.num,
plot.errors.center = plot.errors.center,
plot.errors.type = plot.errors.type,
plot.errors.quantiles = plot.errors.quantiles,
bws = bws)[["boot.err"]]

pc = (plot.errors.center == "bias-corrected")

lerr = matrix(data = if(pc) {terr[,3]} else {eval(tcomp)}
-terr[,1],
nrow = x1.neval, ncol = x2.neval, byrow = FALSE)

herr = matrix(data = if(pc) {terr[,3]} else {eval(tcomp)}
+terr[,2],
nrow = x1.neval, ncol = x2.neval, byrow = FALSE)

} else if (plot.errors.method == "asymptotic") {
lerr = matrix(data = eval(tcomp) - 2.0*tobj\$derr,
nrow = x1.neval, ncol = x2.neval, byrow = FALSE)

herr = matrix(data = eval(tcomp) + 2.0*tobj\$derr,
nrow = x1.neval, ncol = x2.neval, byrow = FALSE)

}

if(is.null(zlim)) {
zlim =
if (plot.errors)
c(min(lerr),max(herr))
else
c(min(eval(tcomp)),max(eval(tcomp)))
}

tret = parse(text=paste(
"npdensity",
"(bws = bws, eval = x.eval,",
"dens",
" = eval(tcomp), derr = terr[,1:2], ntrain = nrow(xdat))", sep=""))

if (plot.behavior != "plot"){
d1 = eval(tret)
d1\$bias = NA

if (plot.errors.center == "bias-corrected")
d1\$bias = terr[,3] - eval(tcomp)

if (plot.behavior == "data")
return ( list(d1 = d1) )
}

# rows = constant x2
# cols = constant x1

dtheta = 5.0
dphi = 10.0

persp.col = ifelse(plot.errors, FALSE, ifelse(!is.null(col),col,"lightblue"))

for (i in 0:((360 %/% dtheta - 1)*rotate)*dtheta+theta){
if (plot.errors){
persp(x1.eval,
x2.eval,
lerr,
zlim = zlim,
cex.axis = ifelse(!is.null(cex.axis),cex.axis,par()\$cex.axis),
cex.lab = ifelse(!is.null(cex.lab),cex.lab,par()\$cex.lab),
cex.main = ifelse(!is.null(cex.main),cex.main,par()\$cex.main),
cex.sub = ifelse(!is.null(cex.sub),cex.sub,par()\$cex.sub),
col = persp.col,
border = ifelse(!is.null(border),border,"grey"),
ticktype = "detailed",
xlab = "",
ylab = "",
zlab = "",
theta = i,
phi = phi,
lwd = ifelse(!is.null(lwd),lwd,par()\$lwd))
par(new = TRUE)
}

persp(x1.eval,
x2.eval,
tdens,
zlim = zlim,
col = persp.col,
border = ifelse(!is.null(border),border,"black"),
ticktype = "detailed",
cex.axis = ifelse(!is.null(cex.axis),cex.axis,par()\$cex.axis),
cex.lab =  ifelse(!is.null(cex.lab),cex.lab,par()\$cex.lab),
cex.main = ifelse(!is.null(cex.main),cex.main,par()\$cex.main),
cex.sub = ifelse(!is.null(cex.sub),cex.sub,par()\$cex.sub),
xlab = ifelse(!is.null(xlab),xlab,gen.label(names(xdat)[1], "X1")),
ylab = ifelse(!is.null(ylab),ylab,gen.label(names(xdat)[2], "X2")),
zlab = ifelse(!is.null(zlab),zlab,"Joint Density"),
theta = i,
phi = phi,
main = gen.tflabel(!is.null(main), main, paste("[theta= ", i,", phi= ", phi,"]", sep="")))

if (plot.errors){
par(new = TRUE)
persp(x1.eval,
x2.eval,
herr,
zlim = zlim,
cex.axis = ifelse(!is.null(cex.axis),cex.axis,par()\$cex.axis),
cex.lab = ifelse(!is.null(cex.lab),cex.lab,par()\$cex.lab),
cex.main = ifelse(!is.null(cex.main),cex.main,par()\$cex.main),
cex.sub = ifelse(!is.null(cex.sub),cex.sub,par()\$cex.sub),
col = persp.col,
border = ifelse(!is.null(border),border,"grey"),
ticktype = "detailed",
xlab = "",
ylab = "",
zlab = "",
theta = i,
phi = phi,
lwd = ifelse(!is.null(lwd),lwd,par()\$lwd))
}

Sys.sleep(0.5)
}

if (plot.behavior == "plot-data")
return ( list(d1 = d1) )

} else {

if (plot.behavior != "data" && plot.par.mfrow)
par(mfrow=n2mfrow(bws\$ndim),cex=par()\$cex)

ev = xdat[1,,drop = FALSE]

for (i in 1:bws\$ndim)
ev[1,i] = uocquantile(xdat[,i], prob=xq[i])

maxneval = max(c(sapply(xdat,nlevels),neval))

exdat = as.data.frame(matrix(data = 0, nrow = maxneval, ncol = bws\$ndim))

for (i in 1:bws\$ndim)
exdat[,i] = ev[1,i]

if (common.scale){
data.eval = matrix(data = NA, nrow = maxneval, ncol = bws\$ndim)
data.err = matrix(data = NA, nrow = maxneval, ncol = 3*bws\$ndim)
allei = as.data.frame(matrix(data = NA, nrow = maxneval, ncol = bws\$ndim))
all.bxp = list()
}

plot.out = list()

temp.err = matrix(data = NA, nrow = maxneval, ncol = 3)
temp.dens = replicate(maxneval, NA)

## plotting expressions

plot.bootstrap = plot.errors.method == "bootstrap"

pfunE = expression(ifelse(xi.factor,
ifelse(plot.bootstrap & plot.bxp,"bxp","plotFactor"), "plot"))

pxE = expression(ifelse(common.scale,
ifelse(xi.factor,
ifelse(plot.bootstrap & plot.bxp, "z = all.bxp[[i]],", "f = allei[,i],"),
"x = allei[,i],"),
ifelse(xi.factor,
ifelse(plot.bootstrap & plot.bxp, "z = temp.boot,", "f = ei,"), "x = ei,")))

pyE = expression(ifelse(xi.factor & plot.bootstrap & plot.bxp, "",
ifelse(common.scale,"y = data.eval[,i],", "y = temp.dens,")))

pylimE = ifelse(common.scale, "ylim = c(y.min,y.max),",
ifelse(plot.errors, "ylim = c(min(na.omit(c(temp.dens - temp.err[,1], temp.err[,3] - temp.err[,1]))),
max(na.omit(c(temp.dens + temp.err[,2], temp.err[,3] + temp.err[,2])))),", ""))

pxlabE = "xlab = ifelse(!is.null(xlab),xlab,gen.label(bws\$xnames[i], paste('X', i, sep = ''))),"

pylabE = paste("ylab = ", "ifelse(!is.null(ylab),ylab,'Density')", ",")

prestE = expression(ifelse(xi.factor,"", "type = ifelse(!is.null(type),type,'l'), lty = ifelse(!is.null(lty),lty,par()\$lty), col = ifelse(!is.null(col),col,par()\$col), lwd = ifelse(!is.null(lwd),lwd,par()\$lwd), cex.axis = ifelse(!is.null(cex.axis),cex.axis,par()\$cex.axis), cex.lab = ifelse(!is.null(cex.lab),cex.lab,par()\$cex.lab), cex.main = ifelse(!is.null(cex.main),cex.main,par()\$cex.main), cex.sub = ifelse(!is.null(cex.sub),cex.sub,par()\$cex.sub),"))

pmainE = "main = ifelse(!is.null(main),main,''), sub = ifelse(!is.null(sub),sub,''),"

## error plotting expressions
plotOnEstimate = (plot.errors.center == "estimate")

efunE = "draw.errors"
eexE = ifelse(common.scale, "ex = as.numeric(na.omit(allei[,i])),",
"ex = as.numeric(na.omit(ei)),")
eelyE = ifelse(common.scale,
ifelse(plotOnEstimate, "ely = na.omit(data.eval[,i] - data.err[,3*i-2]),",
"ely = na.omit(data.err[,3*i] - data.err[,3*i-2]),"),
ifelse(plotOnEstimate, "ely = na.omit(temp.dens - temp.err[,1]),",
"ely = na.omit(temp.err[,3] - temp.err[,1]),"))
eehyE = ifelse(common.scale,
ifelse(plotOnEstimate, "ehy = na.omit(data.eval[,i] + data.err[,3*i-1]),",
"ehy = na.omit(data.err[,3*i] + data.err[,3*i-1]),"),
ifelse(plotOnEstimate, "ehy = na.omit(temp.dens + temp.err[,2]),",
"ehy = na.omit(temp.err[,3] + temp.err[,2]),"))

erestE = "plot.errors.style = ifelse(xi.factor,'bar',plot.errors.style),
plot.errors.bar = ifelse(xi.factor,'I',plot.errors.bar),
plot.errors.bar.num = plot.errors.bar.num,
lty = ifelse(!is.null(lty),lty,ifelse(xi.factor,1,2))"

## density / distribution expressions

devalE = parse(text=paste("npudens",
"(tdat = xdat, edat = subcol(exdat,ei,i)[1:xi.neval,, drop = FALSE], bws = bws)",
sep=""))

dcompE = parse(text="tobj\$dens")

doutE = parse(text=paste("npdensity",
"(bws = bws, eval = subcol(exdat,ei,i)[1:xi.neval,, drop = FALSE],",
"dens",
"= na.omit(temp.dens), derr = na.omit(cbind(-temp.err[,1], temp.err[,2])), ntrain = bws\$nobs)",
sep=""))

for (i in 1:bws\$ndim){
temp.err[,] = NA
temp.dens[] =  NA
temp.boot = list()

xi.factor = is.factor(xdat[,i])

if (xi.factor){
ei = bws\$xdati\$all.ulev[[i]]
xi.neval = length(ei)
} else {
xi.neval = neval
qi = trim.quantiles(xdat[,i], xtrim[i])
ei = seq(qi[1], qi[2], length.out = neval)
}

if (xi.neval < maxneval){
ei[(xi.neval+1):maxneval] = NA
}

tobj = eval(devalE)

temp.dens[1:xi.neval] = eval(dcompE)

if (plot.errors){
if (plot.errors.method == "asymptotic")
temp.err[1:xi.neval,1:2] = replicate(2,2.0*tobj\$derr)
else if (plot.errors.method == "bootstrap"){
temp.boot <- compute.bootstrap.errors(
xdat = xdat,
exdat = subcol(exdat,ei,i)[1:xi.neval,, drop = FALSE],
cdf = FALSE,
slice.index = i,
plot.errors.boot.method = plot.errors.boot.method,
plot.errors.boot.blocklen = plot.errors.boot.blocklen,
plot.errors.boot.num = plot.errors.boot.num,
plot.errors.center = plot.errors.center,
plot.errors.type = plot.errors.type,
plot.errors.quantiles = plot.errors.quantiles,
bws = bws)
temp.err[1:xi.neval,] = temp.boot[["boot.err"]]
temp.boot <- temp.boot[["bxp"]]
if (!plot.bxp.out){
temp.boot\$out <- numeric()
temp.boot\$group <- integer()
}
}
}

if (common.scale){
allei[,i] = ei
data.eval[,i] = temp.dens
if (plot.errors){
all.bxp[i] = NA
all.bxp[[i]] = temp.boot

data.err[,c(3*i-2,3*i-1,3*i)] = temp.err
}
} else if (plot.behavior != "data") {
## plot evaluation
eval(parse(text = paste(eval(pfunE), "(", eval(pxE), eval(pyE),
eval(pylimE), eval(pxlabE), eval(pylabE), eval(prestE),
eval(pmainE), ")")))

## error plotting evaluation
if (plot.errors && !(xi.factor & plot.bootstrap & plot.bxp)){
if (!xi.factor && !plotOnEstimate)
lines(na.omit(ei), na.omit(temp.err[,3]), lty = 3)

eval(parse(text = paste(eval(efunE), "(", eval(eexE), eval(eelyE),
eval(eehyE), eval(erestE), ")")))
}
}

if (plot.behavior != "plot") {
plot.out[i] = NA
plot.out[[i]] = eval(doutE)
plot.out[[i]]\$bias = na.omit(temp.dens - temp.err[,3])
plot.out[[i]]\$bxp = temp.boot
}
}

if (common.scale & (plot.behavior != "data")){
jj = 1:bws\$ndim*3

if (plot.errors.center == "estimate" | !plot.errors) {
y.max = max(na.omit(as.double(data.eval)) +
if (plot.errors) na.omit(as.double(data.err[,jj-1]))
else 0
)
y.min = min(na.omit(as.double(data.eval)) -
if (plot.errors) na.omit(as.double(data.err[,jj-2]))
else 0
)
} else if (plot.errors.center == "bias-corrected") {
y.max = max(na.omit(as.double(data.err[,jj] + data.err[,jj-1])))
y.min = min(na.omit(as.double(data.err[,jj] - data.err[,jj-2])))
}

if(!is.null(ylim)){
y.min = ylim[1]
y.max = ylim[2]
}

for (i in 1:bws\$ndim){
xi.factor = is.factor(xdat[,i])

## plot evaluation
eval(parse(text = paste(eval(pfunE), "(", eval(pxE), eval(pyE),
eval(pylimE), eval(pxlabE), eval(pylabE), eval(prestE),
eval(pmainE), ")")))

## error plotting evaluation
if (plot.errors && !(xi.factor & plot.bootstrap & plot.bxp)){
if (!xi.factor && !plotOnEstimate)
lines(na.omit(ei), na.omit(temp.err[,3]), lty = 3)

eval(parse(text = paste(eval(efunE), "(", eval(eexE), eval(eelyE),
eval(eehyE), eval(erestE), ")")))
}
}
}

if (plot.behavior != "data" && plot.par.mfrow)
par(mfrow=c(1,1),cex=par()\$cex)

if (plot.behavior != "plot"){
names(plot.out) = paste("d",1:bws\$ndim,sep="")
return (plot.out)
}
}
}

npplot.dbandwidth <-
function(bws,
xdat,
data = NULL,
xq = 0.5,
xtrim = 0.0,
neval = 50,
common.scale = TRUE,
perspective = TRUE,
main = NULL,
type = NULL,
```