if (missing(xmin)) xmin <- apply(x, 2, min) - supp*det(Hmax)    if (missing(xmax)) xmax <- apply(x, 2, max) + supp*det(Hmax)        if (d > 4)      stop("Binning only available for 1- to 4-dim data")        if (missing(bgridsize)) bgridsize <- default.gridsize(d)    if (missing(gridsize)) gridsize <- default.gridsize(d)           fhat.list <- list()    for (j in 1:m)    {      xx <- x[x.group==grlab[j],]           H <- Hs[((j-1)*d+1) : (j*d),]                 ## compute individual density estimate      if (binned)        fhat.temp <- kde.binned(x=xx, bgridsize=bgridsize, H=H, xmin=xmin, xmax=xmax)      else if (is.null(eval.points))        fhat.temp <- kde(x=xx, H=H, supp=supp, xmin=xmin, xmax=xmax, gridsize=gridsize)      else        fhat.temp <- kde(x=xx, H=H, eval.points=eval.points)            fhat.list$estimate <- c(fhat.list$estimate, list(fhat.temp$estimate))      fhat.list$eval.points <- fhat.temp$eval.points      fhat.list$x <- c(fhat.list$x, list(xx))      fhat.list$H <- c(fhat.list$H, list(H))    }        fhat.list$x.group <- x.group    pr <- rep(0, length(grlab))    for (j in 1:length(grlab))      pr[j] <- length(which(x.group==grlab[j]))    pr <- pr/nrow(x)    fhat.list$prior.prob <- pr    fhat.list$binned <- binned        class(fhat.list) <- "kda.kde"  }    return(fhat.list)}kda.kde.1d <- function(x, x.group, hs, prior.prob, gridsize, supp, eval.points, binned, bgridsize, xmin, xmax){  grlab <- sort(unique(x.group))  m <- length(grlab)  hmax <- max(hs)  if (missing(xmin)) xmin <- min(x) - supp*hmax  if (missing(xmax)) xmax <- max(x) + supp*hmax    fhat.list <- list()  for (j in 1:m)  {    xx <- x[x.group==grlab[j]]    h <- hs[j]        ## compute individual density estimate    if (binned)      fhat.temp <- kde.binned(x=xx, h=h, xmin=xmin, xmax=xmax, bgridsize=bgridsize)    else if (is.null(eval.points))      fhat.temp <- kde(x=xx, h=h, supp=supp, xmin=xmin, xmax=xmax, gridsize=gridsize)    else      fhat.temp <- kde(x=xx, h=h, eval.points=eval.points)        fhat.list$estimate <- c(fhat.list$estimate, list(fhat.temp$estimate))    fhat.list$eval.points <- fhat.temp$eval.points    fhat.list$x <- c(fhat.list$x, list(xx))    fhat.list$h <- c(fhat.list$h, h)  }      fhat.list$H <- fhat.list$h^2  fhat.list$binned <- binned  fhat.list$x.group <- x.group    if (is.null(prior.prob))  {    pr <- rep(0, length(grlab))    for (j in 1:length(grlab))      pr[j] <- length(which(x.group==grlab[j]))    pr <- pr/length(x)    fhat.list$prior.prob <- pr  }  else    fhat.list$prior.prob <- prior.prob   class(fhat.list) <- "kda.kde"    return(fhat.list)  }################################################################################ Contour method for kda.kde cobjects################################################################################contourLevels.kda.kde <- function(x, prob, cont, nlevels=5, ...) {  fhat <- x  m <- length(fhat$x)  hts <- list()    for (j in 1:m)  {    fhatj <- list(x=fhat$x[[j]], eval.points=fhat$eval.points,                  estimate=fhat$estimate[[j]], H=fhat$H[[j]], binned=fhat$binned)    class(fhatj) <- "kde"    hts[[j]] <- contourLevels(x=fhatj, prob=prob, cont=cont, nlevels=nlevels, ...)  }     return(hts) }############################################################################### Plot KDE of individual densities and partition - only for 2-dim## Parameters# fhat - output from `kda.kde'# y - data points (separate from training data inside fhat)# y.group - data group labels# prior.prob - vector of prior probabilities# disp - "part" - plot partition#      - "" - don't plot partition##############################################################################plot.kda.kde <- function(x, y, y.group, drawpoints=FALSE, ...) {      if (is.vector(x$x[[1]]))    plotkda.kde.1d(x=x, y=y, y.group=y.group, drawpoints=drawpoints, ...)  else  {      d <- ncol(x$x[[1]])        if (d==2)      plotkda.kde.2d(x=x, y=y, y.group=y.group, drawpoints=drawpoints, ...)     else if (d==3)         plotkda.kde.3d(x=x, y=y, y.group=y.group, drawpoints=drawpoints, ...)   }}plotkda.kde.1d <- function(x, y, y.group, prior.prob=NULL, xlim, ylim, xlab="x", ylab="Weighted density function", drawpoints=FALSE, col, partcol, ptcol, lty, jitter=TRUE, ...){   fhat <- x    m <- length(fhat$x)  ##eval1 <- fhat$eval.points    if (is.null(prior.prob))    prior.prob <- fhat$prior.prob    if (m != length(prior.prob))    stop("Prior prob. vector not same length as number of components in fhat")  if (!(identical(all.equal(sum(prior.prob), 1), TRUE)))      stop("Sum of prior weights not equal to 1")  weighted.fhat <- matrix(0, nrow=length(fhat$eval.points), ncol=m)   for (j in 1:m)    weighted.fhat[,j] <- fhat$estimate[[j]]*fhat$prior.prob[j]    if (missing(xlim)) xlim <- range(fhat$eval.points)  if (missing(ylim)) ylim <- range(weighted.fhat)  if (missing(lty)) lty <- rep(1, m)  if (length(lty) < m) lty <- rep(lty, m)  if (missing(col)) col <- 1:m  if (length(col) < m) col <- rep(col, m)  if (missing(ptcol)) ptcol <- col  if (length(ptcol) < m) ptcol <- rep(ptcol, m)  if (missing(partcol)) partcol <- col  if (length(partcol) < m) partcol <- rep(partcol, m)    ## plot each training group's KDE in separate colour and line type   plot(fhat$eval.points, weighted.fhat[,1], type="l", xlab=xlab, ylab=ylab, xlim=xlim, ylim=ylim, lty=lty[1], col=col[1], ...)    if (m > 1)    for (j in 2:m)      lines(fhat$eval.points, weighted.fhat[,j], lty=lty[j], col=col[j], ...)  ##eval.points.gr <- apply(weighted.fhat, 1, which.max)  ydata <- seq(min(fhat$eval.points), max(fhat$eval.points), length=401)  ydata.gr <- kda(unlist(fhat$x), x.group=fhat$x.group, hs=fhat$h, y=ydata, prior.prob=fhat$prior.prob)  ## draw partition class as rug plot with ticks facing inwards    for (j in 1:m)    rug(ydata[ydata.gr==levels(fhat$x.group)[j]], col=partcol[j])    for (j in 1:m)  {      ## draw data points    if (drawpoints)    {      if (missing(y))        if (jitter)          rug(jitter(fhat$x[[j]]), col=ptcol[1], ticksize=-0.03)        else          rug(fhat$x[[j]], col=ptcol[1], ticksize=-0.03)      else       {        if (missing(y.group))          if (jitter)            rug(jitter(y), col=ptcol[1], ticksize=-0.03)          else            rug(y, col=ptcol[1], ticksize=-0.03)        else          if (jitter)            rug(jitter(y[y.group==levels(y.group)[j]]), col=ptcol[j], ticksize=-0.03)          else            rug(y[y.group==levels(y.group)[j]], col=ptcol[j], ticksize=-0.03)       }    }  }   }plotkda.kde.2d <- function(x, y, y.group, prior.prob=NULL,     cont=c(25,50,75), abs.cont, xlim, ylim, xlab, ylab,    drawpoints=FALSE, drawlabels=TRUE, cex=1, pch, lty, col, partcol, ptcol, ...){   fhat <- x    ##d <- 2  m <- length(fhat$x)  ##eval1 <- fhat$eval.points[[1]]  ##eval2 <- fhat$eval.points[[2]]    xtemp <- numeric()  for (j in 1:m)     xtemp <- rbind(xtemp, fhat$x[[j]])   if (missing(xlim)) xlim <- range(xtemp[,1])  if (missing(ylim)) ylim <- range(xtemp[,2])   if (missing(pch)) pch <- 1:m  if (missing(lty)) lty <- rep(1, m)  if (length(lty) < m) lty <- rep(lty, m)  if (missing(col)) col <- 1:m  if (length(col) < m) col <- rep(col, m)  if (missing(partcol)) partcol <- grey.colors(m, start=0.7, end=1)   if (missing(ptcol))    if (missing(y.group))      ptcol <- rep("blue", m)    else      ptcol <- 1:m  if (length(ptcol)==1) ptcol <- rep(ptcol, m)                x.names <- colnames(fhat$x[[1]])   if (!is.null(x.names))  {    if (missing(xlab)) xlab <- x.names[1]    if (missing(ylab)) ylab <- x.names[2]  }  else  {    xlab="x"    ylab="y"  }  if (is.null(prior.prob))    prior.prob <- fhat$prior.prob  if (m != length(prior.prob))    stop("Prior prob. vector not same length as number of components in fhat")  if (!(identical(all.equal(sum(prior.prob), 1), TRUE)))      stop("Sum of prior weights not equal to 1")  ## set up plot  if (missing(y))     plot(fhat$x[[1]], type="n", xlab=xlab, ylab=ylab, xlim=xlim, ylim=ylim, ...)  else    plot(y, type="n", xlab=xlab, ylab=ylab, xlim=xlim, ylim=ylim, ...)    ## set up common grid for all densities   class.grid <- array(0, dim=dim(fhat$est[[1]]))  temp <- matrix(0, ncol=length(fhat$est), nrow=nrow(fhat$est[[1]]))  for (j in 1:ncol(fhat$est[[1]]))  {    for (k in 1:length(fhat$est))      temp[,k] <- fhat$est[[k]][,j]* prior.prob[k]    class.grid[,j] <- max.col(temp)      }    ## draw partition  image(fhat$eval[[1]], fhat$eval[[2]], class.grid, col=partcol, xlim=xlim,        ylim=ylim, add=TRUE, ...)  box()  ## common contour levels removed from >= v1.5.3   if (missing(abs.cont))  {    hts <- contourLevels(fhat, prob=(100-cont)/100)    nhts <- length(hts[[1]])  }  else  {    hts <- abs.cont    nhts <- length(hts)  }   ## draw contours  for (j in 1:m)  {    for (i in 1:nhts)     {      if (missing(abs.cont))      {        scale <- cont[i]/hts[[j]][i]        contour(fhat$eval.points[[1]], fhat$eval.points[[2]],                 fhat$estimate[[j]]*scale, level=hts[[j]][i]*scale, add=TRUE,                 drawlabels=drawlabels, lty=lty[j], col=col[j],                ...)      }      else      {        contour(fhat$eval.points[[1]], fhat$eval.points[[2]],                 fhat$estimate[[j]], level=hts[i], add=TRUE,                 drawlabels=drawlabels, lty=lty[j], col=col[j],                ...)      }    }  }    for (j in 1:m)  {      ## draw data points    if (drawpoints)    {      if (missing(y))        points(fhat$x[[j]], pch=pch[j], col=ptcol[1], cex=cex)      else       {        if (missing(y.group))          points(y, col=ptcol[1], cex=cex)        else          points(y[y.group==levels(y.group)[j],], pch=pch[j], col=ptcol[j], cex=cex)       }    }  }   }plotkda.kde.3d <- function(x, y, y.group, prior.prob=NULL,    cont=c(25,50,75), abs.cont, colors, alphavec, xlab, ylab, zlab,    drawpoints=FALSE, size=3, ptcol="blue", ...){  require(rgl)  require(misc3d)    fhat <- x     ##d <- 3  m <- length(fhat$x)   ##eval1 <- fhat$eval.points[[1]]  ##eval2 <- fhat$eval.points[[2]]  ##eval3 <- fhat$eval.points[[3]]  if (is.null(prior.prob))    prior.prob <- fhat$prior.prob  if (m != length(prior.prob))    stop("Prior prob. vector not same length as number of components in fhat")  if (!(identical(all.equal(sum(prior.prob), 1), TRUE)))      stop("Sum of prior weights not equal to 1")  x.names <- colnames(fhat$x[[1]])  if (missing(xlab))    if (is.null(x.names)) xlab <- "x" else xlab <- x.names[1]  if (missing(ylab))    if (is.null(x.names)) ylab <- "y" else ylab <- x.names[2]  if (missing(zlab))    if (is.null(x.names)) zlab <- "z" else zlab <- x.names[3]               ##dobs <- numeric(0)  xx <- numeric(0)  for (j in 1:m)    xx <- rbind(xx, fhat$x[[j]])  ##x.gr <- sort(unique(fhat$x.group))  ##if (fhat$binned)  ##  bin.par.xx <- dfltCounts.ks(xx, gridsize=dim(fhat$est[[j]]), sqrt(diag(fhat$H[[j]])), supp=3.7)  ## common contour levels removed from >= v1.5.3   if (missing(abs.cont))  {    hts <- contourLevels(fhat, prob=(100-cont)/100)    nhts <- length(hts[[1]])  }  else  {    hts <- abs.cont    nhts <- length(hts)  }    if (missing(alphavec)) alphavec <- seq(0.1,0.3,length=nhts)  if (missing(colors)) colors <- rainbow(m)  if (missing(ptcol))    if (missing(y.group))      ptcol <- rep("blue", m)    else      ptcol <- 1:m  if (length(ptcol)==1) ptcol <- rep(ptcol, m)    clear3d()  ##bg3d(color="white")  plot3d(x=fhat$eval.points[[1]], y=fhat$eval.points[[2]],         z=fhat$eval.points[[3]], type="n", xlab=xlab, ylab=ylab, zlab=zlab,         ...)    for (j in 1:m)  {    for (i in 1:nhts)       contour3d(x=fhat$eval.points[[1]], y=fhat$eval.points[[2]],                z=fhat$eval.points[[3]], f=fhat$estimate[[j]],                level=hts[[j]][nhts-i+1],                add=TRUE, alpha=alphavec[i], color=colors[j],...)    if (drawpoints)   ## plot points    {      if (missing(y))        points3d(fhat$x[[j]][,1], fhat$x[[j]][,2], fhat$x[[j]][,3],                    color=ptcol[j], size=size, alpha=1)      else      {        if (missing(y.group))          points3d(y[,1], y[,2], y[,3], color=ptcol, size=size, alpha=1)        else        {          y.temp <- y[y.group==levels(y.group)[j],]          if (nrow(y.temp)>0)            points3d(y.temp[,1], y.temp[,2], y.temp[,3], color=ptcol[j], size=size, alpha=1)        }      }    }  }}