"fig1.1" = 
function(ps = F)
{
  if(ps)
    ps.start("spen1.ps", c(4.5, 2.8))
  if(is.R())
    data(spencer)
  plot(spencer$age, spencer$mortality, xlab = "Age (Years)", ylab = 
    "Mortality Rate")
  fit <- lm(mortality ~ age, data = spencer)
  lines(spencer$age, fitted(fit))
  if(ps)
    graphics.off()
}

"fig1.2" = 
function(ps = F, trel = !is.R())
{
  if(ps)
    ps.start("spen2.ps", c(4.5, 5), trellis = T)
  if(is.R())
    data(spencer)
  xl <- "Age (Years)"
  yl <- "Mortality Rate"
  if(trel) {
    lv <- factor(1:2, labels = c("21-point rule", "15-point rule"))
    print(xyplot(range(spencer$mortality) ~ I(1:2) | lv, xlim = c(20, 45), xlab
       = xl, ylab = yl, panel = function(x, y, ...)
    {
      panel.xyplot(spencer$age, spencer$mortality)
      if(x == 2)
        grad <- spence.15(spencer$mortality)
      else grad <- spence.21(spencer$mortality)
      panel.xyplot(spencer$age, grad, type = "l")
    }
    , strip = strip.loc, layout = c(1, 2)))
  }
  else {
    plot(spencer$age, spencer$mortality, xlab = xl, ylab = yl)
    lines(spencer$age, spence.15(spencer$mortality))
    plot(spencer$age, spencer$mortality, xlab = xl, ylab = yl)
    lines(spencer$age, spence.21(spencer$mortality))
  }
  if(ps)
    graphics.off()
}

"fig10.1" = 
function(ps = F)
{
  if(ps)
    ps.start("olfs4.ps", c(4.5, 2.5), trellis = T)
  if(is.R())
    data(geyser)
  a <- c(0.5, 0.75, 1.)
  show(xyplot(a ~ c(1, 4, 6) | paste("k =", a), ylim = c(0, 60), xlab = 
    "Eruption Duration", ylab = "(f''(x))^2", panel = function(x, y, ...)
  {
    fit <- locfit.raw(geyser, deg = 2, link = "ident", kern = "gauss", ev = 
      grid(501, ll = 1, ur = 6), alpha = c(0, y), deriv = c(1, 1))
    z <- knots(fit, what = c("x", "coef"))
    panel.xyplot(z[, 1], z[, 2]^2, type = "l")
    itgl <- (sum(c(1, rep(2, 499), 1) * z[, 2]^2) * 5)/1000
    cat("integral: ", itgl, "\n")
  }
  , strip = strip.loc))
  if(ps)
    graphics.off()
  invisible(NULL)
}

"fig10.2" = 
function(ps = F)
{
  if(ps)
    ps.start("olfs3.ps", c(4.5, 2.5))
  if(is.R())
    data(geyser)
  gf <- 2.5
  a <- seq(0.05, 0.7, length.out = 100)
  z <- sjpi(geyser, seq(0.05, 0.7, length.out = 100))
  print(gf * z)
  plot(gf * z[, 2], gf * z[, 1], type = "l", xlab = "Pilot Bandwidth k", ylab
     = "Bandwidth h")
  lines(gf * z[, 3], gf * z[, 1], lty = 2)
  lines(gf * a, (gf * a)/107^(1/10), lty = 3)
  legend(gf * 0.05, gf * 0.4, lty = c(1, 2, 3), legend = c("Plug-in", 
    "SJ assumed", "GKK assumed"))
  if(ps)
    graphics.off()
  invisible(NULL)
}

"fig10.3" = 
function(ps = F)
{
  if(ps)
    ps.start("olfs1.ps", c(4.5, 3), trellis = T)
  if(is.R())
    data(geyser)
  h <- c(0.162, 0.315, 0.249, 0.705, 0.516, 0.625)
  mm <- c("AIC", "LCV", "LSCV", "BCV", "SJPI", "Visual")
  z <- matrix(ncol = 3, nrow = 0)
  k <- length(mm)
  for(i in 1:k) {
    fit <- locfit( ~ geyser, alpha = c(0, h[i]), kern = "gauss", deg = 0, ev = 
      grid(200, ll = 1, ur = 6))
    z <- rbind(z, cbind(knots(fit, what = c("x", "coef")), i))
  }
  me <- factor(z[, 3], labels = paste(mm, ": h=", h, sep = "")[1:k])
  show(xyplot(z[, 2] ~ z[, 1] | me, panel = function(x, y)
  {
    panel.xyplot(x, y, type = "l")
    rug(geyser)
  }
  , xlab = "Eruption Duration (minutes)", ylab = "Density", strip = strip.loc))
  if(ps)
    graphics.off()
  invisible(NULL)
}

"fig10.4" = 
function(ps = F, nsim = 0)
{
  if(ps)
    ps.start("olfs2.ps", c(4.5, 3), trellis = T)
  if(is.R())
    data(geyser)
  nmmt <- c("LSCV", "BCV", "SJPI")
  if(nsim > 0) {
    bb <- matrix(nrow = nsim, ncol = 6)
    for(i in 1:nsim) {
      print(i)
      x <- sample(geyser, 107, replace = T) + rnorm(107, sd = 0.219)
      bb[i, 1:3] <- kdeb(x, 0.01, 1., meth = nmmt)
      x <- sample(geyser, 107, replace = T) + rnorm(107, sd = 0.069)
      bb[i, 4:6] <- kdeb(x, 0.01, 1., meth = nmmt)
    }
    assign("ofsim", bb, where = 1)
  }
  nsim <- nrow(ofsim)
  bw <- as.numeric(ofsim)
  nm <- length(nmmt)
  meth <- factor(rep(rep(1:nm, rep(nsim, nm)), 2), levels = 1:nm, labels = nmmt
    )
  sig <- factor(rep(c("sigma=0.219", "sigma=0.070"), rep(nsim * nm, 2)))
  y <- rep(c(0.219, 0.07), c(nsim * nm, nsim * nm))
  show(xyplot(y ~ bw | (meth * sig), xlim = c(0.1, 0.8), ylim = c(0, 10.2), 
    panel = function(x, y)
  {
    fit <- locfit( ~ x, ev = rbox(ll = 0.1, ur = 0.8))
    lines(fit)
    rug(x)
    if(y[1] == 0.219)
      abline(v = 0.516, lty = 2)
    else abline(v = 0.249, lty = 2)
  }
  , xlab = "Bandwidth", ylab = "Density", strip = strip.loc))
  if(ps)
    graphics.off()
  invisible(NULL)
}

"fig10.5" = 
function(ps = F)
{
  if(ps)
    ps.start("claw2.ps", c(4.5, 2.5), c(1, 2))
  if(is.R())
    data(claw54)
  fit1 <- locfit( ~ claw54, deg = 0, kern = "gauss", alpha = c(0, 0.315), ev = 
    grid(100, ll = -3.5, ur = 2.7))
  fit2 <- locfit( ~ claw54, deg = 0, kern = "gauss", alpha = c(0, 0.985), ev = 
    grid(100, ll = -3.5, ur = 2.7))
  x <- seq(-3.5, 2.7, length.out = 200)
  y <- dnorm(x, -1., 0.1) + dnorm(x, -0.5, 0.1) + dnorm(x, 0, 0.1) + dnorm(x, 
    0.5, 0.1) + dnorm(x, 1., 0.1)
  y <- (y + 5 * dnorm(x))/10
  plot(fit1, get.data = T, main = "h=0.315", ylim = c(0, max(y)))
  lines(x, y, lty = 2)
  plot(fit2, get.data = T, main = "h=0.985", ylim = c(0, max(y)))
  lines(x, y, lty = 2)
  if(ps)
    graphics.off()
  invisible(NULL)
}

"fig10.6" = 
function(ps = F, nsim = 0, pl = T)
{
  meth <- c("LSCV", "BCV", "SJPI")
  if(nsim > 0) {
    rclaw <- function(n)
    {
      mu <- c(0, seq(-2, 2, by = 1))/2
      sd <- c(1, rep(0.1, 5))
      pr <- c(0.5, rep(0.1, 5))
      z <- sample(1:6, n, replace = T, prob = pr)
      mu[z] + sd[z] * rnorm(n)
    }
    claws <- matrix(nrow = nsim, ncol = 9)
    for(i in 1:nsim) {
      print(i)
      x <- rclaw(54)
      print(x)
      claws[i, 1:3] <- kdeb(x, meth = meth)
      x <- rclaw(193)
      claws[i, 4:6] <- kdeb(x, meth = meth)
      x <- rclaw(400)
      claws[i, 7:9] <- kdeb(x, meth = meth)
      print(claws[i,  ])
    }
    assign("claws", claws, where = 1)
  }
  if(pl) {
    if(ps)
      ps.start("claw1.ps", c(4.5, 4.5), trellis = T)
    nsim <- nrow(claws)
    z <- rep(0, 900)
    m <- 1:100
    x <- seq(0.05, 1.75, length.out = nsim)
    clawsim <- as.numeric(claws)
    ssize <- factor(rep(1:3, 3 * rep(nsim, 3)), labels = c("n=54", "n=193", 
      "n=400"))
    meth <- factor(rep(rep(1:3, rep(nsim, 3)), 3), labels = c("LSCV", "BCV", 
      "SJPI"))
    yy <- rep(c(54, 193, 400), rep(3 * nsim, 3))
    claw.pan <- function(x, y, subscripts)
    {
      fit <- locfit( ~ x, data = data.frame(x = x), alpha = 0.3)
      xev <- seq(0.05, 1.75, length = 200)
      z <- as(predict(fit, xev), "numeric")
      panel.xyplot(xev, z/max(z), type = "l")
      if(y[1] == 54) {
        abline(v = c(0.315, 0.985), lty = 2)
      }
      else if(y[1] == 193) {
        abline(v = 0.195, lty = 2)
      }
      else abline(v = 0.16, lty = 2)
      rug(x)
    }
    print(xyplot(yy ~ clawsim | meth * ssize, xlim = c(0.05, 1.75), data = 
      data.frame(clawsim = clawsim, yy = yy, ssize = ssize, meth = meth), ylim
       = c(0, 1), xlab = "Bandwidth", ylab = "Density", panel = claw.pan, strip
       = strip.loc))
    if(ps)
      graphics.off()
  }
  invisible(NULL)
}

"fig10.7" = 
function(ps = F)
{
  if(ps)
    ps.start("aisp1.ps", trellis = T, c(4.5, 3))
  if(is.R())
    data(ais)
  print(dotplot(sport ~ LBM | sex, data = ais, strip = strip.loc))
  if(ps)
    graphics.off()
  invisible(NULL)
}

"fig10.8" = 
function(ps = F)
{
  if(ps)
    ps.start("aisp2.ps", trellis = T, c(4.5, 3))
  if(is.R())
    data(ais)
  print(xyplot(c(0, 0.05) ~ c(2.658367, 9.386367) | as.factor(c("LSCV h=2.658", 
    "SJPI h=9.386")), xlim = c(30, 110), panel = function(x, y, ...)
  {
    fit <- locfit( ~ LBM, data = ais, kern = "gauss", deg = 0, alpha = c(0, x), 
      ev = grid(100, ll = 30, ur = 110))
    plot(fit, add = T)
    rug(ais$LBM)
  }
  , xlab = "LBM", ylab = "Density", strip = strip.loc))
  if(ps)
    graphics.off()
  invisible(NULL)
}

"fig10.9" = 
function(ps = F)
{
  if(ps)
    ps.start("aisp3.ps", c(4.5, 3), c(1, 2))
  if(is.R())
    data(ais)
  gf <- 2.5
  h <- c(1, 1.05, 1.1, 1.15, 1.2, 1.3, 1.4, 1.6, 1.8, 2, 2.25, 2.5, 3, 3.5, 4, 
    5, 6, 7, 8, 10, 12.5, 15, 17.5, 20, 23, 26, 30, 35, 40, 45)
  plot(lscvplot(ais$LBM, alpha = h, exact = T), type = "l")
  z <- sjpi(ais$LBM, h/4)
  plot(gf * z[, 2], gf * z[, 1], type = "l", xlab = "Pilot Bandwidth", ylab = 
    "Selected Bandwidth")
  lines(gf * z[, 3], gf * z[, 1], lty = 2)
  legend(10, 3, lty = 1:2, cex = 0.8, legend = c("plug-in", "SJ ass'd"))
  if(ps)
    graphics.off()
  invisible(NULL)
}

"fig11.1" = 
function(ps = F)
{
  if(ps)
    ps.start("ethm2.ps", c(4.5, 2.3), mfrow = c(1, 2))
  if(is.R())
    data(ethanol)
  fit <- locfit(NOx ~ E, data = ethanol, alpha = 0.3, deg = c(0, 3))
  plot(fit, get.data = T)
  plot(preplot(fit, where = "fitp", what = "deg"), type = "p", ylab = "Degree")
  if(ps)
    graphics.off()
  invisible(NULL)
}

"fig11.2" = 
function(ps = F)
{
  if(ps)
    ps.start(file = "mcyc4.ps", c(4, 3))
  tmp <- matrix(c(1.45, 1.68732, 1.885, 1.40911, 2.4505, 1.29482, 3.18565, 
    1.34986, 4.14135, 1.34611, 5.38375, 1.33071, 6.99887, 1.29307, 9.09853, 
    1.23902, 11.8281, 1.31084, 15.37652, 1.36923, 19.98948, 2.15223, 10.64529, 
    1.27582, 9.58076, 1.24452, 8.62268, 1.23951, 7.76041, 1.25705), ncol = 2, 
    byrow = T)
  plot(tmp[, 1], tmp[, 2], xlab = "Bandwidth", ylab = "Local CP", type = "n")
  text(tmp[, 1], tmp[, 2], 1:nrow(tmp), cex = 0.8)
  if(ps)
    graphics.off()
}

"fig11.3" = 
function(ps = F, col = 1:2)
{
  if(ps)
    ps.start(file = "mcyc5.ps", c(4.5, 5), mfrow = c(2, 1))
  if(is.R())
    data(mcyc)
  fit <- locfit(accel ~ time, alpha = c(0, 0, 2), weights = 1/vp, data = mcyc.n,
    acri = "cp")
  plot(mcyc$time, mcyc$accel, type = "n", xlab = "Acceleration", ylab = "Time")
  x <- knots(fit, what = "x")
  coef <- predict(fit, where = "fitp")
  h <- predict(fit, where = "fitp", what = "band")
  points(x, coef)
  segments(x - h, coef, x + h, coef, col = col[2])
  plot(fit, get.data = T, xlab = "Acceleration", ylab = "Time")
  if(ps)
    graphics.off()
}

"fig11.4" = 
function(ps = F)
{
  if(ps)
    ps.start("dopl1.ps", c(4.5, 6), c(3, 1), mar = c(3, 3, 1, 1) + 0.1)
  .Random.seed <- c(13, 5, 27, 35, 18, 2, 63, 16, 55, 35, 50, 0)
  x <- seq(0, 1, length.out = 2048)
  mean <- 20 * sqrt(x * (1 - x)) * sin((2 * pi * 1.05)/(x + 0.05))
  y <- mean + rnorm(2048)
  plot(x, y, pch = ".")
  fit <- locfit.raw(x, y, alpha = c(0, 0, 4), acri = "cp", maxk = 500)
  print(summary(fit))
  print(sum((fitted(fit) - mean)^2))
  plot(fit, m = 2048)
  plot(preplot(fit, what = "band", where = "fitp"), type = "p", log = "y", ylab
     = "Bandwidth")
  if(ps)
    graphics.off()
}

"fig11.5" = 
function(ps = F)
{
  if(ps)
    ps.start("dopl3.ps", c(4.5, 4), c(2, 1), mar = c(3, 3, 1, 1) + 0.1)
  .Random.seed <- c(13, 5, 27, 35, 18, 2, 63, 16, 55, 35, 50, 0)
  x <- seq(0, 1, length.out = 2048)
  mean <- 20 * sqrt(x * (1 - x)) * sin((2 * pi * 1.05)/(x + 0.05))
  y <- mean + rnorm(2048)
  fit <- locfit.raw(x, y, alpha = c(0, 0, 1.1), acri = "ici", maxk = 500)
  print(summary(fit))
  print(sum((fitted(fit) - mean)^2))
  plot(fit, m = 2048)
  plot(preplot(fit, what = "band", where = "fitp"), type = "p", log = "y", ylab
     = "bandwidth")
  if(ps)
    graphics.off()
}

"fig11.6" = 
function(ps = F)
{
  if(ps)
    ps.start("oldf4.ps", c(4.5, 2.8))
  if(is.R())
    data(geyser)
  geyser.round <<- data.frame(duration = seq(1.05, 5.95, by = 0.05), count = 
    as.numeric(table(cut(geyser, breaks = seq(1.025, 5.975, length.out = 100)))
    ))
  fit <- locfit(count ~ duration, weights = rep(107 * 0.05, 99), data = 
    geyser.round, alpha = c(0, 0, 2), acri = "cp", family = "poisson")
  plot(fit, get.data = T, m = 200, xlab = "Eruption Duration", ylab = "Density"
    )
  if(ps)
    graphics.off()
  fit
}

"fig12.1" = 
function(ps = F)
{
  if(ps)
    ps.start(file = "evst2.ps", c(4, 2.8))
  plot(0, 0, type = "l", xlim = c(-0.3, 2.), ylim = c(-0.5, 1.5), xlab = "", 
    ylab = "")
  lines(c(0, 1.5, 1.5, 0, 0), c(0, 0, 1, 1, 0))
  segments(c(0.75, 0.75), c(0, 0.5), c(0.75, 1.5), c(1., 0.5), lty = 2)
  text(c(0, 1.5, 0, 1.5, 0.75, 0.75, 0.55, 1.7), c(-0.1, -0.1, 1.1, 1.1, -0.1, 
    1.1, 0.5, 0.5), c("v0 h=1.8", "v1 h=1.0", "v2 h=1.6", "v3 h=1.1", 
    "v4 h=1.3", "v5 h=1.5", "v6 h=1.4", "v7 h=1.0"))
  if(ps)
    graphics.off()
  invisible(NULL)
}

"fig12.2" = 
function(ps = F)
{
  if(ps)
    ps.start(file = "evst5.ps", c(4.5, 3), c(1, 2), mar = c(2, 2, 1, 0) + 0.1)
  fit <- locfit(formula =  ~ x0 + x1, data = trimod, deg = 1, alpha = 0.35, ev
     = rbox(type = "kdtree"))
  plot.eval(fit)
  points(trimod$x0, trimod$x1, cex = 0.3)
  fit <- locfit(formula =  ~ x0 + x1, data = trimod, deg = 1, alpha = 0.35, ev
     = rbox(cut = 0.85))
  plot.eval(fit)
  points(trimod$x0, trimod$x1, cex = 0.3)
  if(ps)
    graphics.off()
  invisible(NULL)
}

"fig12.3" = 
function(ps = F)
{
  if(ps)
    ps.start("code1.ps", c(2.5, 2.5))
  v <- rbind(c(0, 0), c(1, 0), c(0, 1), c(1, 1))
  x <- c(0.4, 0.6)
  plot(v[c(1, 2, 4, 3, 1), 1], v[c(1, 2, 4, 3, 1), 2], type = "b", xlab = "x0", 
    ylab = "x1")
  lines(c(0, 1), c(x[2], x[2]), lty = 2)
  points(x[1], x[2], pch = "X")
  text(rep(c(0.08, 0.92), 3), rep(c(0.08, 0.92, 0.68), c(2, 2, 2)), paste("v", 
    0:5, sep = ""), adj = 0.5)
  if(ps)
    graphics.off()
  invisible(NULL)
}

"fig13.1" = 
function(ps = F)
{
  if(ps)
    ps.start("minm1.ps", c(4.5, 3))
  if(is.R())
    data(mmsamp)
  fit <- locfit(y ~ x, data = mmsamp, deg = 1, kern = "minmax", alpha = 4000, 
    ev = "grid", mg = 200, flim = c(0, 1))
  print(fit@dp)
  plot(preplot(fit, where = "fitp", get.data = T))
  fit <- locfit(y ~ x, data = mmsamp, deg = 1, alpha = c(0, 0.05), ev = "grid", 
    mg = 100, flim = c(0, 1))
  print(fit@dp)
  lines(preplot(fit, where = "fitp"), lty = 2)
  xx <- seq(0, 1, length.out = 100)
  yy <- 2 - 5 * xx + 5 * exp( - (20 * xx - 10)^2)
  lines(xx, yy, lty = 3)
  legend(0, -3, legend = c("Minimax", "Constant h", "True mean"), lty = 1:3)
  if(ps)
    graphics.off()
  invisible(NULL)
}

"fig13.2" = 
function(ps = F)
{
  if(ps)
    ps.start("minm2.ps", c(3.5, 3.5))
  m <- c(100, 200, 300, 500, 750, 1000, 1500, 2000, 3000, 5000, 7500, 10000)
  plot(m, apply(minmax, 2, mean), ylim = c(0, 120), xlab = "M", ylab = "S.S.E", 
    type = "b", log = "x")
  points(m, apply(mindex2, 2, mean), pch = "+", type = "b")
  points(m[1:4], apply(mband, 2, mean)[1:4], pch = "*", type = "b")
  par(xaxt = "s")
  legend(3., 120, c("Minimax", "M-Index", "h constant"), pch = c("o", "+", "*")
    )
  if(ps)
    graphics.off()
  invisible(NULL)
}

"fig2.1" = 
function()
{
  if(is.R())
    data(ethanol)
  demo.point <- function(x0, x, y, h, col)
  {
    W <- function(u)
    {
      u <- abs(u)
      v <- pmax(1 - u * u * u, 0)
      v * v * v
    }
    rx <- range(x)
    rx[1] <- max(rx[1], x0 - h)
    rx[2] <- min(rx[2], x0 + h)
    xx <- seq(rx[1], rx[2], length.out = 100)
    yy <- W((xx - x0)/h) - 1
    lines(xx, yy, col = col[2])
    my.dat <- data.frame(x = x, y = y, w = W((x0 - x)/0.2))
    fit <- lm(y ~ x + I(x * x), data = my.dat, weights = w)
    fv <- predict(fit, data.frame(x = xx))
    lines(xx, fv, col = col[3])
    points(x0, predict(fit, data.frame(x = x0)), cex = 2)
  }
  col <- 1:3
  x <- ethanol$E
  y <- ethanol$NOx
  plot(NOx ~ E, data = ethanol, ylim = c(-1, 4.2), pch = "+")
  demo.point(min(ethanol$E), ethanol$E, ethanol$NOx, 0.2, col)
  demo.point(0.95, ethanol$E, ethanol$NOx, 0.15, col)
}

"fig2.2" = 
function()
{
  if(is.R())
    data(ethanol)
  fit <- locfit(NOx ~ E, alpha = 0.5, data = ethanol)
  col <- 1:10
  plot(fit, get.data = T, col = col[2])
}

"fig2.3" = 
function(ps = F, trel = !is.R())
{
  if(ps)
    ps.start("baet1.ps", c(4.5, 4), trellis = T)
  if(is.R())
    data(ethanol)
  a <- c(0.2, 0.4, 0.6, 0.8)
  xl <- "Equivalence Ratio"
  yl <- "NOx"
  if(trel) {
    show(xyplot(a ~ I(a) | as.factor(paste("a =", a)), xlim = range(ethanol$E), 
      ylim = range(ethanol$NOx), panel = function(x, ...)
    {
      xd <- ethanol$E
      yd <- ethanol$NOx
      fit <- locfit.raw(xd, yd, alpha = x)
      panel.xyplot(xd, yd, cex = 0.7)
      lines.locfit(fit)
    }
    , strip = strip.loc, xlab = xl, ylab = yl))
  }
  else {
    for(i in 1:4) {
      fit <- locfit(NOx ~ E, data = ethanol, alpha = a[i])
      plot(fit, get.data = T, main = paste("a =", a[i]), xlab = xl, ylab = yl)
    }
  }
  if(ps)
    graphics.off()
}

"fig2.4" = 
function(ps = F)
{
  if(ps)
    ps.start("ethn5.ps", c(4.5, 4), trellis = T)
  if(is.R())
    data(ethanol)
  d <- 0:3
  a <- c(0.25, 0.3, 0.49, 0.59)
  f <- factor(1:4, labels = c("Local Constant", "Local Linear", 
    "Local Quadratic", "Local Cubic"))
  show(xyplot(a ~ d | f, xlim = range(ethanol$E), ylim = range(ethanol$NOx), 
    panel = function(x, y, ...)
  {
    xd <- ethanol$E
    yd <- ethanol$NOx
    fit <- locfit.raw(xd, yd, alpha = y, deg = x, ev = grid(100))
    panel.xyplot(xd, yd, cex = 0.7)
    lines.locfit(fit)
  }
  , strip = strip.loc, xlab = "Equivalence Ratio", ylab = "NOx"))
  if(ps)
    graphics.off()
}

"fig2.5" = 
function(ps = F)
{
  if(ps)
    ps.start("baet2.ps", c(4.5, 4), trellis = T)
  if(is.R())
    data(ethanol)
  a <- c(0.2, 0.4, 0.6, 0.8)
  show(xyplot(a ~ I(a) | as.factor(paste("a =", a)), xlim = range(ethanol$E), 
    ylim = c(-1, 1), panel = function(x, ...)
  {
    xd <- ethanol$E
    fit <- locfit.raw(xd, ethanol$NOx, alpha = x)
    print(fit)
    res <- residuals(fit)
    fit2 <- locfit.raw(xd, res, alpha = 0.2)
    panel.xyplot(xd, res, cex = 0.7)
    lines.locfit(fit2)
    panel.abline(h = 0, lty = 2)
  }
  , strip = strip.loc, xlab = "Equivalence Ratio", ylab = "Residual"))
  if(ps)
    graphics.off()
}

"fig2.6" = 
function(ps = F)
{
  if(ps)
    ps.start("baet4.ps", c(4.5, 2.5), trellis = T)
  if(is.R())
    data(ethanol)
  d <- 2:3
  a <- c(0.49, 0.59)
  show(xyplot(a ~ d | factor(1:2, labels = c("Local Quadratic", "Local Cubic")),
    xlim = range(ethanol$E), ylim = c(0, 0.5), panel = function(x, y, ...)
  {
    xd <- ethanol$E
    yd <- ethanol$NOx
    fit <- locfit.raw(xd, yd, alpha = y, deg = x, ev = grid(100))
    print(fit@dp)
    plot(preplot(fit, where = "fitp", what = "infl"), add = T)
    plot(preplot(fit, where = "fitp", what = "vari"), add = T, lty = 2)
    rug(xd)
  }
  , strip = strip.loc, xlab = "Equivalence Ratio", ylab = "Influence"))
  if(ps)
    graphics.off()
}

"fig2.7" = 
function(ps = F)
{
  if(ps)
    ps.start("demo3.ps", c(4.5, 3), c(1, 2))
  if(is.R())
    data(ethanol)
  alpha <- seq(0.2, 0.8, by = 0.05)
  plot(gcvplot(NOx ~ E, data = ethanol, alpha = alpha, df = 3), ylim = c(0, 0.2
    ))
  plot(cpplot(NOx ~ E, data = ethanol, alpha = alpha), ylim = c(0, 85))
  abline(0, 1)
  if(ps)
    graphics.off()
}

"fig3.1" = 
function(ps = F)
{
  if(ps)
    ps.start("ethn1.ps", c(4.5, 7), mfrow = c(2, 1))
  if(is.R())
    data(ethanol)
  fit <- locfit(NOx ~ E + C, data = ethanol, alpha = 0.5, scale = 0)
  plot(fit, get.data = T)
  plot(fit, type = "persp")
  if(ps)
    graphics.off()
}

"fig3.2" = 
function(ps = F)
{
  if(ps)
    ps.start("ethn4.ps", c(4.5, 4.5), mfrow = c(2, 1))
  if(is.R())
    data(ethanol)
  fit <- locfit(NOx ~ E + C, data = ethanol, alpha = 0.5, scale = 0)
  show(plot(fit, pv = "C", tv = "E", mtv = 9, get.data = T))
  if(ps)
    graphics.off()
}

"fig3.3" = 
function()
{
  if(is.R())
    data(ethanol)
  fit.ad <- gam(NOx ~ lf(E, alpha = 0.5) + lf(C, deg = 1), data = ethanol)
  plot(fit.ad)
}

"fig4.1" = 
function(ps = F)
{
  if(ps)
    ps.start("minf1.ps", c(4, 2.8))
  if(is.R())
    data(mine)
  fit <- locfit(frac ~ extrp, data = mine, family = "poisson", deg = 1, alpha
     = 0.6)
  plot(fit, band = "g", get.data = T, xlab = "Extraction Percentage", ylab = 
    "Fractures")
  if(ps)
    graphics.off()
}

"fig4.2" = 
function(ps = F)
{
  if(ps)
    ps.start("hsmo1.ps", c(4.5, 2.8))
  if(is.R())
    data(morths)
  fit <- locfit(deaths ~ age, weights = n, family = "binomial", data = morths, 
    alpha = 0.5)
  plot(fit, get.data = T, band = "g", xlab = "Age (years)", ylab = 
    "Death Proportion")
  if(ps)
    graphics.off()
}

"fig4.3" = 
function(ps = F)
{
  if(ps)
    ps.start("hsmo2.ps", c(4.5, 3.8), c(2, 2))
  if(is.R())
    data(morths)
  fit <- locfit(deaths ~ age, weights = n, family = "binomial", data = morths, 
    alpha = 0.5)
  for(ty in c("deviance", "pearson", "response", "ldot")) {
    res <- residuals(fit, type = ty)
    plot(morths$age, res, main = ty, type = "b", xlab = "Age", ylab = 
      "Residual")
    abline(h = 0, lty = 2)
  }
  if(ps)
    graphics.off()
}

"fig4.4" = 
function(ps = F)
{
  if(ps)
    ps.start("minf2.ps", c(3.5, 2.8))
  if(is.R())
    data(mine)
  a <- seq(0.4, 1, by = 0.05)
  show(aicplot(frac ~ extrp, data = mine, family = "poisson", deg = 1, alpha = 
    a))
  if(ps)
    graphics.off()
}

"fig5.1" = 
function(ps = F)
{
  if(ps)
    ps.start("oldf1.ps", c(4.5, 2.5))
  if(is.R())
    data(geyser)
  fit <- locfit( ~ geyser, alpha = c(0.1, 0.8), ev = rbox(ll = 1, ur = 6))
  plot(fit, m = 200, xlab = "Old Faithful Eruption Duration", ylab = "Density", 
    get.data = T)
  if(ps)
    graphics.off()
}

"fig5.2" = 
function(ps = F)
{
  if(ps)
    ps.start("oldf2.ps", c(4.5, 2.5))
  if(is.R())
    data(geyser)
  fit <- locfit( ~ geyser, alpha = c(0.1, 0.6), ev = rbox(ll = 1, ur = 6), link
     = "ident")
  plot(fit, m = 200, xlab = "Old Faithful Eruption Duration", ylab = "Density", 
    get.data = T)
  if(ps)
    graphics.off()
}

"fig5.3" = 
function(ps = F)
{
  if(ps)
    ps.start("stam1.ps", c(4.5, 3.2), colors = ps.colors)
  if(is.R())
    data(stamp)
  n <- sum(stamp$count)
  loc.stamp <<- cbind(stamp, w = rep(0.001 * n, 76))
  fit <- locfit(count ~ thick, weights = w, data = loc.stamp, family = 
    "poisson", alpha = c(0, 0.004))
  plot(fit, m = 200, get.data = T, ylab = "Density", xlab = "Thickness (m.m.)")
  if(ps)
    graphics.off()
}

"fig5.4" = 
function(ps = F)
{
  if(ps)
    ps.start("trim1.ps", c(4, 3), ...)
  if(is.R())
    data(trimod)
  fit.trim <- locfit( ~ x0 + x1, data = trimod, alpha = 0.35)
  plot(fit.trim, type = "persp")
  if(ps)
    graphics.off()
}

"fig5.5" = 
function(ps = F)
{
  if(ps)
    ps.start("trim2.ps", c(4, 3.5))
  if(is.R())
    data(trimod)
  fit.trim <- locfit( ~ x0 + x1, data = trimod, alpha = 0.35)
  emp <- sort(fitted(fit.trim))
  plot(fit.trim, v = emp[floor(c(0.05, 0.5) * 225)])
  points(trimod$x0, trimod$x1, cex = 0.5)
  if(ps)
    graphics.off()
}

"fig5.6" = 
function(ps = F)
{
  if(ps)
    ps.start("akde2.ps", c(4.5, 2.8))
  if(is.R())
    data(geyser)
  fit <- locfit( ~ geyser, alpha = c(0.1, 1.2), ev = rbox(ll = 1, ur = 6), 
    renorm = T)
  x <- seq(1, 6, by = 0.01)
  z <- predict(fit, x)
  plot(x, 0.01 * cumsum(z), type = "l", xlab = "Eruption Duration (Minutes)", 
    ylab = "Cumulative Distribution Function")
  lines(sort(geyser), (1:107)/107, type = "s")
  if(ps)
    graphics.off()
}

"fig5.7" = 
function(ps = F)
{
  if(ps)
    ps.start("akde3.ps", c(4.5, 7), c(3, 2), mar = c(2, 2, 0.5, 0) + 0.1)
  if(is.R())
    data(geyser.round)
  for(a in c(1.2, 0.8, 0.5)) {
    fit <- locfit(count ~ duration, data = geyser.round, weights = rep(107 * 
      0.05, 99), alpha = c(0.1, a), family = "poisson")
    plot(fit, get.data = T, xlab = "", ylab = "")
    res <- residuals(fit)
    fitr <- locfit.raw(geyser.round$duration, res, alpha = 0.1)
    plot(fitr, get.data = T, xlab = "", ylab = "")
    abline(h = 0, lty = 2)
  }
  if(ps)
    graphics.off()
  invisible(NULL)
}

"fig6.1" = 
function(ps = F)
{
  if(ps)
    ps.start("oldf5.ps", c(4.5, 2.5))
  if(is.R())
    data(geyser)
  fit1 <- locfit( ~ geyser, alpha = c(0.1, 0.6), ev = rbox(ll = 1, ur = 6))
  fit2 <- locfit( ~ geyser, alpha = c(0.1, 0.6), ev = rbox(ll = 1, ur = 6), 
    deriv = 1)
  z <- lfmarg(fit1, 200)
  plot(preplot(fit1, z) * preplot(fit2, z), xlab = 
    "Eruption Duration (minutes)", ylab = "Density Derivative")
  if(ps)
    graphics.off()
}

"fig6.2" = 
function(ps = F)
{
  if(ps)
    ps.start("codi1.ps", c(4.5, 2.5))
  if(is.R())
    data(co2)
  fit1 <- locfit(co2 ~ I(year + month/12), data = co2, alpha = 0.5, deg = 1)
  plot(fit1, ylim = range(co2$co2), xlab = "date", ylab = "co2")
  points(co2$year + co2$month/12, co2$co2, cex = 0.5)
  if(ps)
    graphics.off()
}

"fig6.3" = 
function(ps = F)
{
  if(ps)
    ps.start("codi2.ps", c(3.5, 2.5))
  if(is.R())
    data(co2)
  fit1 <- locfit(co2 ~ I(year + month/12), data = co2, alpha = 0.5, deg = 1)
  loc.co2 <<- cbind(co2, res = residuals(fit1))
  fit2 <- locfit(res ~ ang(year + month/12), data = loc.co2, scale = 1/(2 * pi),
    alpha = c(0, 2))
  plot(fit2, xlim = c(0, 1))
  if(ps)
    graphics.off()
}

"fig6.4" = 
function(ps = F)
{
  if(ps)
    ps.start("codi3.ps", c(4.5, 2.5))
  if(is.R())
    data(co2)
  fit1 <- locfit(co2 ~ I(year + month/12), data = co2, alpha = 0.5, deg = 1)
  loc.co2 <<- cbind(co2, res = residuals(fit1))
  fit2 <- locfit(res ~ ang(year + month/12), data = loc.co2, scale = 1/(2 * pi),
    alpha = c(0, 2))
  plot(co2$year + co2$month/12, fitted(fit1) + fitted(fit2), type = "l")
  if(ps)
    graphics.off()
}

"fig6.5" = 
function(ps = F)
{
  if(ps)
    ps.start("codi4.ps", c(4, 4))
  if(is.R())
    data(co2)
  fit <- locfit(co2 ~ ang(month) + I(year + month/12), data = co2, scale = c(12/
    (2 * pi), 10), alpha = 0.2)
  plot(fit, xlab = "Month", ylab = "Year")
  if(ps)
    graphics.off()
}

"fig6.6" = 
function(ps = F)
{
  if(ps)
    ps.start("penn1.ps", c(4.5, 2.8))
  if(is.R())
    data(penny)
  midp <- (1945:1988) + 0.5
  fitl <- locfit(thickness ~ left(year), data = penny, alpha = c(0, 10), deg = 
    1, ev = midp)
  fitr <- locfit(thickness ~ right(year), data = penny, alpha = c(0, 10), deg
     = 1, ev = midp)
  plot((preplot(fitr) - preplot(fitl))^2, type = "b", xlab = "Year", ylab = 
    "Delta^2(t)")
  if(ps)
    graphics.off()
}

"fig6.7" = 
function(ps = F)
{
  if(ps)
    ps.start("penn2.ps", c(4.5, 2.8))
  if(is.R())
    data(penny)
  fit0 <- locfit(thickness ~ year, data = penny, alpha = c(0, 10), deg = 1, 
    subset = (year <= 1958))
  fit1 <- locfit(thickness ~ year, data = penny, alpha = c(0, 10), deg = 1, 
    subset = (year >= 1959) & (year <= 1974))
  fit2 <- locfit(thickness ~ year, data = penny, alpha = c(0, 10), deg = 1, 
    subset = (year >= 1975))
  plot(penny$year, penny$thickness, xlab = "Year", ylab = 
    "Thickness (inches*1000)")
  lines(fit0)
  lines(fit1)
  lines(fit2)
  if(ps)
    graphics.off()
}

"fig7.1" = 
function(ps = F)
{
  if(ps)
    ps.start("hear1.ps", c(4.5, 2.8))
  if(is.R())
    data(heart)
  fit <- locfit( ~ surv, data = heart, cens = cens, family = "hazard", alpha = 
    0.4, xlim = c(0, 100000))
  plot(fit, m = 300, ylim = c(0, 0.004), xlab = "Survival Time", ylab = 
    "Hazard Rate", get.data = T)
  if(ps)
    graphics.off()
}

"fig7.2" = 
function(ps = F)
{
  if(ps)
    ps.start("livm1.ps", c(4.5, 4))
  if(is.R())
    data(livmet)
  fit <- locfit( ~ t + dm, data = livmet, cens = 1 - z, scale = 0, deg = 1, 
    family = "hazard", alpha = 0.5, xlim = list(t = c(0, 10000)))
  plot(fit, ylab = "Diameter (c.m.)", xlab = "Survival Time (months)", get.data
     = T)
  if(ps)
    graphics.off()
}

"fig7.3" = 
function(ps = F)
{
  if(ps)
    ps.start("hear2.ps", c(4.5, 4))
  if(is.R())
    data(heart)
  col <- 1:10
  plotbyfactor(heart$age, 0.5 + heart$surv, heart$cens, pch = c("O", "+"), col
     = col, xlab = "Age at Transplant (years)", ylab = 
    "0.5+Survival Time (Days)", ylim = c(0.5, 16000), log = "y")
  fit <- locfit(log10(0.5 + surv) ~ age, data = heart)
  lines.locfit(fit, tr = function(x)
  10^x)
  fit <- locfit.censor(log10(0.5 + surv) ~ age, cens = cens, data = heart)
  lines.locfit(fit, lty = 2, tr = function(x)
  10^x)
  fit <- locfit.censor(log10(0.5 + surv) ~ age, cens = cens, data = heart, km
     = T)
  lines.locfit(fit, lty = 3, tr = function(x)
  10^x)
  par(yaxt = "s")
  legend(12, log10(16000), legend = c("Raw data", "Normal model", 
    "Kaplan-Meier"), lty = 1:3)
  if(ps)
    graphics.off()
}

"fig7.4" = 
function(ps = F)
{
  if(ps)
    ps.start("cric1.ps", c(4.5, 2.8), colors = colors)
  if(is.R())
    data(border)
  fit0 <- locfit(runs ~ day, cens = no, data = border, family = "geom", alpha
     = 0.7)
  plot(fit0, get.data = T)
  if(ps)
    graphics.off()
}

"fig7.5" = 
function(ps = F)
{
  if(ps)
    ps.start("cric3.ps", c(4.5, 2.8))
  if(is.R())
    data(border)
  fit0 <- locfit(runs ~ day, cens = no, data = border, family = "geom", alpha
     = 0.7)
  fit1 <- locfit(runs ~ day, weights = rep(0.8, 265), cens = no, data = border, 
    family = "geom", alpha = 0.7)
  plot(0.8 * preplot(fit1) - preplot(fit0))
  if(ps)
    graphics.off()
}

"fig7.6" = 
function(ps = F)
{
  if(ps)
    ps.start("hear3.ps", c(4.5, 4))
  if(is.R())
    data(heart)
  col <- 1:10
  fit <- locfit(I((surv + 0.5)^0.625) ~ age, cens = cens, data = heart, alpha
     = 0.8, family = "gamma")
  y <- log(heart$surv + 0.5)
  plotbyfactor(heart$age, heart$surv + 0.5, heart$cens, pch = c("O", "+"), col
     = col, log = "y", xlab = "Age at Transplant (years)", ylab = 
    "0.5+Survival Time (Days)")
  plot(fit, tr = function(x)
  exp(x/0.625) * gamma(1 + 1/0.625), add = T)
  if(ps)
    graphics.off()
}

"fig8.1" = 
function(ps = F)
{
  if(ps)
    ps.start("disc1.ps", c(3.5, 3.5))
  if(is.R()) {
    data(cltrain)
    data(cltest)
  }
  fit <- locfit(I(y == 1) ~ x1 + x2, data = cltrain, scale = 0)
  plot(fit, v = 0.5)
  text(cltrain$x1, cltrain$x2, cltrain$y, cex = 0.7)
  t1 <- table(fitted(fit) > 0.5, cltrain$y)
  t2 <- table(predict(fit, cltest) > 0.5, cltest$y)
  if(ps)
    graphics.off()
  cbind(t1, t2)
}

"fig8.2" = 
function(ps = F)
{
  if(ps)
    ps.start("disc2.ps", c(3.5, 3.5))
  if(is.R())
    data(cltrain)
  if(is.R())
    data(cltest)
  fit1 <- locfit( ~ x1 + x2, data = cltrain, subset = y == 0, family = "rate", 
    scale = 0)
  fit2 <- locfit( ~ x1 + x2, data = cltrain, subset = y == 1, family = "rate", 
    scale = 0)
  pr <- lfmarg(c(-3, -2.2, 3, 2), c(50, 50))
  id <- function(x)
  x
  plot(preplot(fit2, pr, tr = id) - preplot(fit1, pr, tr = id), v = 0)
  text(cltest$x1, cltest$x2, cltest$y, cex = 0.7)
  if(ps)
    graphics.off()
  fiy1 <- predict(fit1, cltrain, tr = id)
  fiy2 <- predict(fit2, cltrain, tr = id)
  t3 <- table(fiy2 - fiy1 > 0, cltrain$y)
  pry1 <- predict(fit1, cltest, tr = id)
  pry2 <- predict(fit2, cltest, tr = id)
  t4 <- table(pry2 - pry1 > 0, cltest$y)
  cbind(t3, t4)
}

"fig8.3" = 
function(ps = F)
{
  if(ps)
    ps.start("irid2.ps", c(4, 3), colors = colors)
  if(is.R())
    data(iris)
  a <- (2:9)/10
  z <- numeric(length(a))
  for(i in 1:length(a)) {
    fit <- locfit(I(species == "virginica") ~ petal.wid + petal.len, data = 
      iris, scale = 0, deg = 1, alpha = a[i], ev = dat(cv = T))
    tab <- table(fitted(fit) < 0.5, iris$species == "virginica")
    z[i] <- sum(diag(tab))
  }
  fit <- locfit(I(species == "virginica") ~ petal.wid + petal.len, data = iris, 
    scale = 0, deg = 1)
  plotbyfactor(petal.wid, petal.len, species, data = iris, pch = c("o", "+"), 
    lg = c(1, 6.8))
  plot(fit, v = 0.5, add = T)
  if(ps)
    graphics.off()
  cbind(a, z)
}

"fig8.4" = 
function(ps = F, col = rep(1, 3))
{
  if(ps)
    ps.start("chdb1.ps", c(4, 4))
  if(is.R())
    data(chemdiab)
  plotbyfactor(fpg, ga, cc, data = chemdiab, col = col, pch = c("X", "O", "+"), 
    xlab = "Fasting Plasma Glucose", ylab = "Glucose Area", lg = c(70, 1550))
  fit1 <- locfit(I(cc == "Overt_Diabetic") ~ fpg + ga, data = chemdiab, scale
     = 0, deg = 1, ev = dat(cv = T))
  fit2 <- locfit(I(cc == "Chemical_Diabetic") ~ fpg + ga, data = chemdiab, 
    scale = 0, deg = 1, ev = dat(cv = T))
  fit3 <- locfit(I(cc == "Normal") ~ fpg + ga, data = chemdiab, scale = 0, deg
     = 1, ev = dat(cv = T))
  p1 <- fitted(fit1)
  p2 <- fitted(fit2)
  p3 <- fitted(fit3)
  z <- (p1 > pmax(p2, p3)) + 2 * (p2 > pmax(p1, p3)) + 3 * (p3 > pmax(p1, p2))
  print(table(chemdiab$cc, z))
  if(ps)
    graphics.off()
}

"fig8.5" = 
function(ps = F, col = rep(1, 4))
{
  if(ps)
    ps.start("kang1.ps", c(4, 4))
  if(is.R())
    data(kangaroo)
  plotbyfactor(cr.w, asc.h, spec, data = kangaroo, col = col[2:4], pch = c("X", 
    "O", "+"), lg = c(10, 620), xlab = "Crest Width", ylab = 
    "Ascending Ramus Height")
  fit1 <- locfit(I(spec == "giganteus") ~ cr.w + asc.h, data = kangaroo, scale
     = 0, deg = 1, ev = dat(T))
  fit2 <- locfit(I(spec == "melanops") ~ cr.w + asc.h, data = kangaroo, scale
     = 0, deg = 1, ev = dat(T))
  fit3 <- locfit(I(spec == "fuliginosus") ~ cr.w + asc.h, data = kangaroo, 
    scale = 0, deg = 1, ev = dat(T))
  p1 <- fitted(fit1)
  p2 <- fitted(fit2)
  p3 <- fitted(fit3)
  z <- (p1 > pmax(p2, p3)) + 2 * (p2 > pmax(p1, p3)) + 3 * (p3 > pmax(p1, p2))
  if(ps)
    graphics.off()
  table(kangaroo$spec, z)
}

"fig8.6" = 
function(ps = F, col = rep(1, 4))
{
  if(ps)
    ps.start("kang2.ps", c(4, 4))
  if(is.R())
    data(kangaroo)
  Z <- ZZ <- as.matrix(kangaroo[, 3:20])
  sp <- kangaroo$spec
  for(i in 1:18)
    ZZ[, i] <- residuals(aov(Z[, i] ~ sp))
  Sigma <- var(ZZ)
  eig <- eigen(Sigma)
  Z <- Z %*% eig$vectors/rep(sqrt(eig$values), rep(101, 18))
  m1 <- apply(Z[sp == "giganteus",  ], 2, mean)
  m2 <- apply(Z[sp == "melanops",  ], 2, mean)
  m3 <- apply(Z[sp == "fuliginosus",  ], 2, mean)
  x1 <- Z %*% (m1 - m2)
  x2 <- Z %*% (m2 - m3)
  plotbyfactor(x1, x2, sp, col = col[2:4], pch = c("X", "O", "+"), lg = c(2, 
    -50))
  fit1 <- locfit.raw(cbind(x1, x2), sp == "giganteus", scale = 0, deg = 1, ev
     = dat(T))
  fit2 <- locfit.raw(cbind(x1, x2), sp == "melanops", scale = 0, deg = 1, ev = 
    dat(T))
  fit3 <- locfit.raw(cbind(x1, x2), sp == "fuliginosus", scale = 0, deg = 1, ev
     = dat(T))
  p1 <- fitted(fit1)
  p2 <- fitted(fit2)
  p3 <- fitted(fit3)
  z <- (p1 > pmax(p2, p3)) + 2 * (p2 > pmax(p1, p3)) + 3 * (p3 > pmax(p1, p2))
  fit1 <- locfit.raw(cbind(x1, x2), sp == "giganteus", scale = 0, deg = 1)
  fit2 <- locfit.raw(cbind(x1, x2), sp == "melanops", scale = 0, deg = 1)
  fit3 <- locfit.raw(cbind(x1, x2), sp == "fuliginosus", scale = 0, deg = 1)
  newx <- lfmarg(c(min(x1), min(x2), max(x1), max(x2)), c(50, 50))
  p1 <- matrix(predict(fit1, newx), ncol = 50)
  p2 <- matrix(predict(fit2, newx), ncol = 50)
  p3 <- matrix(predict(fit3, newx), ncol = 50)
  contour(newx[[1]], newx[[2]], p1 - pmax(p2, p3), v = 0, add = T, labex = 0)
  contour(newx[[1]], newx[[2]], p2 - pmax(p1, p3), v = 0, add = T, labex = 0)
  contour(newx[[1]], newx[[2]], p3 - pmax(p1, p2), v = 0, add = T, labex = 0)
  if(ps)
    graphics.off()
  table(sp, z)
}

"fig8.7" = 
function(ps = F)
{
  if(ps)
    ps.start("clde1.ps", c(4.5, 2.8))
  if(is.R())
    data(cldem)
  fit <- locfit(y ~ x, data = cldem)
  z <- summary(fit)
  print(z)
  print(fit)
  plot(fit, get.data = T)
  x <- seq(min(cldem$x), max(cldem$x), length.out = 100)
  lines(x, exp(x)/(1 + exp(x)), lty = 2)
  abline(h = 0.5, lty = 3)
  polygon(c(-0.1086, 0, 0, -0.1086), c(-0.03, -0.03, 1.03, 1.03), density = 30, 
    lty = 2)
  if(ps)
    graphics.off()
  invisible(NULL)
}

"fig9.1" = 
function(ps = F)
{
  if(ps)
    ps.start("mcyc1.ps", c(4.5, 2.7))
  if(is.R())
    data(mcyc)
  plot(mcyc$time, mcyc$accel, xlab = "Time", ylab = "Acceleration")
  if(ps)
    graphics.off()
}

"fig9.2" = 
function(ps = F, pl = T)
{
  if(ps)
    ps.start("mcyc2.ps", c(4.5, 2.7))
  if(is.R())
    data(mcyc)
  fit <- locfit(accel ~ time, data = mcyc, alpha = 0.1)
  x <- knots(fit, what = "x")
  y <- -2 * predict(fit, what = "lik", where = "fitp")
  w <- predict(fit, what = "rdf", where = "fitp")
  fitv <- locfit(y ~ x, weights = w, family = "gamma", alpha = 0.4)
  if(pl)
    plot(fitv, get.data = T, xlab = "Time", ylab = "Variance")
  if(ps)
    graphics.off()
  invisible(predict(fitv, mcyc$time))
}

"fig9.3" = 
function(ps = F, pl = T)
{
  if(ps)
    ps.start("hear4.ps", c(4.5, 3))
  if(is.R())
    data(heart)
  fit <- locfit(I((surv + 0.5)^0.625) ~ age, cens = cens, data = heart, alpha
     = 0.8, family = "gamma", deriv = 1)
  crit(fit) <- kappa0(fit)
  plot(fit, band = "global", xlab = "Age at Transplant", ylab = "Local Slope")
  abline(h = 0, lty = 3)
  if(ps)
    graphics.off()
  invisible(fit)
}

"strip.loc" = 
function(which.given, which.panel, var.name, factor.levels, shingle.intervals, 
  par.strip.text = trellis.par.get("add.text"), strip.names = c(F, T), style = 
  3)
{
  if(!missing(strip.names) && length(strip.names) == 1)
    strip.names <- rep(strip.names, 2)
  which <- which.panel[which.given]
  var.name <- var.name[which.given]
  if(!is.null(factor.levels)) {
    strip.names <- if(strip.names[1]) paste(var.name, factor.levels[which], sep
         = ": ") else factor.levels[which]
    n <- length(factor.levels)
    do.call("text", c(list(x = 0.5, y = 0.5, labels = strip.names), 
      par.strip.text))
  }
  else if(!is.null(shingle.intervals)) {
    strip.shingle <- trellis.par.get("strip.shingle")
    strip.background <- trellis.par.get("strip.background")
    bgdcol <- strip.background$col[(which.given - 1) %% length(strip.background$
      col) + 1]
    polygon(c(0, 1, 1, 0), c(0, 0, 1, 1), border = F, col = bgdcol)
    xminmax <- shingle.intervals[which,  ]
    e <- par("cxy")[1]/5
    if(xminmax[1] == xminmax[2]) {
      xminmax[1] <- xminmax[1] - e
      xminmax[2] <- xminmax[2] + e
    }
    par(usr = par("usr") + c( - e, e, 0, 0))
    polygon(xminmax[c(1:2, 2:1)], c(0, 0, 1, 1), border = F, col = 
      strip.shingle$col[(which.given - 1) %% length(strip.shingle$col) + 1])
    if(strip.names[2])
      do.call("text", c(list(x = 0.5, y = 0.5, labels = var.name), 
        par.strip.text))
  }
  else stop("strip must be either a factor or a shingle")
}

"is.R" = 
function()
exists("version") && !is.null(vl <- version$language) && vl == "R"

"store" = 
function()
{
  dump(c("fig1.1", "fig1.2", "fig10.1", "fig10.2", "fig10.3", "fig10.4", 
    "fig10.5", "fig10.6", "fig10.7", "fig10.8", "fig10.9", "fig11.1", "fig11.2",
    "fig11.3", "fig11.4", "fig11.5", "fig11.6", "fig12.1", "fig12.2", "fig12.3",
    "fig13.1", "fig13.2", "fig2.1", "fig2.2", "fig2.3", "fig2.4", "fig2.5", 
    "fig2.6", "fig2.7", "fig3.1", "fig3.2", "fig3.3", "fig4.1", "fig4.2", 
    "fig4.3", "fig4.4", "fig5.1", "fig5.2", "fig5.3", "fig5.4", "fig5.5", 
    "fig5.6", "fig5.7", "fig6.1", "fig6.2", "fig6.3", "fig6.4", "fig6.5", 
    "fig6.6", "fig6.7", "fig7.1", "fig7.2", "fig7.3", "fig7.4", "fig7.5", 
    "fig7.6", "fig8.1", "fig8.2", "fig8.3", "fig8.4", "fig8.5", "fig8.6", 
    "fig8.7", "fig9.1", "fig9.2", "fig9.3", "strip.loc", "is.R", "store"), 
    "../S/lffigs.s")
}