kda.r

r软件 另一款可以计算核估计的软件包 需安装r软件

    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)
        }
      }
    }
  }
}