kda.r

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

##############################################################################
# Kernel discriminant analysis
###############################################################################


###############################################################################
# Find bandwidths for each class in training set, for 2- to 6-dim 
#
# Parameters
# x - data values
# group - group variable
# bw - type of bandwidth selector
# nstage, pilot, pre - parameters for plugin bandwidths
# diag - FALSE - use full b/w matrices
#      - TRUE - use diag b/w matrices
#
# Returns
# Matrix of bandwidths for each group in training set
###############################################################################

hkda <- function(x, x.group, bw="plugin", nstage=2, binned=TRUE, bgridsize)
{
  grlab <- sort(unique(x.group))
  m <- length(grlab)
  bw <- substr(tolower(bw),1,1)
  hs <- numeric(0)

  if (missing(bgridsize)) bgridsize <- 401
  
  for (i in 1:m)
  {
    y <- x[x.group==grlab[i]]
    if (bw=="p")
      h <- hpi(y, nstage=nstage, binned=TRUE, bgridsize=bgridsize) 
    hs <- c(hs, h)
  }

  return(hs)
}
   
Hkda <- function(x, x.group, Hstart, bw="plugin", nstage=2, pilot="samse",
                 pre="sphere", binned=FALSE, bgridsize)
{
  d <- ncol(x)
  grlab <- sort(unique(x.group))
  m <- length(grlab)
  bw <- substr(tolower(bw),1,1)
  Hs <- numeric(0)

  if (missing(bgridsize) & binned) bgridsize <- default.gridsize(d)
  
  for (i in 1:m)
  {
    y <- x[x.group==grlab[i],]
    if (!missing(Hstart)) 
    {
      Hstarty <- Hstart[((i-1)*d+1) : (i*d),]
      if (bw=="l")        
        H <- Hlscv(y, Hstart=Hstarty)
      else if (bw=="s")
        H <- Hscv(y, pre=pre, Hstart=Hstarty, binned=binned, bgridsize=bgridsize)
      else if (bw=="p")
        H <- Hpi(y, nstage=nstage, pilot=pilot, pre=pre, Hstart=Hstarty,
                 binned=binned, bgridsize=bgridsize) 
    }
    else
    {
      if (bw=="l")
        H <- Hlscv(y)
      else if (bw=="s")
        H <- Hscv(y, pre=pre, binned=binned, bgridsize=bgridsize)
      else if (bw=="p")
        H <- Hpi(y, nstage=nstage, pilot=pilot, pre=pre, binned=binned, bgridsize=bgridsize)
    }
    Hs <- rbind(Hs, H)
  }

  return(Hs)   
}

Hkda.diag <- function(x, x.group, bw="plugin", nstage=2, pilot="samse",
                 pre="sphere", binned=FALSE, bgridsize)
{
  d <- ncol(x)
  grlab <- sort(unique(x.group))
  m <- length(grlab)
  bw <- substr(tolower(bw),1,1)
  Hs <- numeric(0)

  if (missing(bgridsize) & binned) bgridsize <- default.gridsize(d)

  for (i in 1:m)
  {
    y <- x[x.group==grlab[i],]
    if (bw=="l")
      H <- Hlscv.diag(y, binned=binned, bgridsize=bgridsize)
    else if (bw=="p") 
      H <- Hpi.diag(y, nstage=nstage, pilot=pilot, pre=pre, binned=binned, bgridsize=bgridsize)
    else if (bw=="s")
      H <- Hscv.diag(y, pre=pre, binned=binned, bgridsize=bgridsize)
    Hs <- rbind(Hs, H)
  }

  return(Hs)   
}



###############################################################################
# Classify data set according to discriminant analysis based on training data
# for 1- to 6-dim
#
# Parameter
# x - training data
# x.group - group variable for x
# y - data values to be classified
# Hs - bandwidth matrices
# prior.prob - prior probabilities
#
# Returns
# Group classification of data set y
###############################################################################

kda <- function(x, x.group, Hs, hs, y, prior.prob=NULL)
{
  if (is.vector(x))
  {
    disc.gr <- kda.1d(x=x, x.group=x.group, hs=hs, y=y, prior.prob=prior.prob)
  }
  else
  {  
    if (is.data.frame(x)) x <- as.matrix(x)
    if (is.data.frame(y)) y <- as.matrix(y)
    gr <- sort(unique(x.group))

    ##d <- ncol(x)
   
    ## if prior.prob is NULL then use sample proportions
    if (is.null(prior.prob))
    {
      prior.prob <- rep(0, length(gr))
      for (j in 1:length(gr))
        prior.prob[j] <- length(which(x.group==gr[j]))
      prior.prob <- prior.prob/nrow(x)
    }
    
    if (!(identical(all.equal(sum(prior.prob), 1), TRUE)))  
      stop("Sum of prior weights not equal to 1")
    
    ## Compute KDE and weighted KDE 
    m <- length(gr)
    fhat <- kda.kde(x, x.group, Hs=Hs, eval.points=y)
    fhat.wt <- matrix(0, ncol=m, nrow=nrow(y))  
    
    for (j in 1:m)
      fhat.wt[,j] <- fhat$est[[j]]* prior.prob[j]
    
    ## Assign y according largest weighted density value 
    disc.gr.temp <- apply(fhat.wt, 1, which.max)
    
    disc.gr <- gr
    for (j in 1:m)
    {
      ind <- which(disc.gr.temp==j)
      disc.gr[ind] <- gr[j]
    }
  }
  
  return(disc.gr) 
}

kda.1d <- function(x, x.group, hs, y, prior.prob=NULL)
{ 
  gr <- sort(unique(x.group))

  # if prior.prob is NULL then use sample proportions
  if (is.null(prior.prob))
  {
    prior.prob <- rep(0, length(gr))
    for (j in 1:length(gr))
      prior.prob[j] <- length(which(x.group==gr[j]))
    prior.prob <- prior.prob/length(x)
  }
  
  if (!(identical(all.equal(sum(prior.prob), 1), TRUE)))  
    stop("Sum of prior weights not equal to 1")

  ## Compute KDE and weighted KDE 
  m <- length(gr)
  fhat <- kda.kde(x, x.group, hs=hs, eval.points=y)
  fhat.wt <- matrix(0, ncol=m, nrow=length(y))  
  
  for (j in 1:m)
    fhat.wt[,j] <- fhat$estimate[[j]]* prior.prob[j]

  ## Assign y according largest weighted density value 
  disc.gr.temp <- apply(fhat.wt, 1, which.max)

  disc.gr <- gr
  for (j in 1:m)
  {
    ind <- which(disc.gr.temp==j)
    disc.gr[ind] <- gr[j]
  }
 
  return(disc.gr) 
}


###############################################################################
# Compares true group classification with an estimated one
#
# Parameters
# group - true group variable
# est.group - estimated group variable
#
# Returns
# List with components
# comp - cross-classification table of groupings - true groups are the rows,
#        estimated groups are the columns
# error - total mis-classification rate
###############################################################################

compare <- function(x.group, est.group, by.group=FALSE)
{
  if (length(x.group)!=length(est.group))
    stop("Group label vectors not the same length")
  
  grlab <- sort(unique(x.group))
  m <- length(grlab)
  comp <- matrix(0, nr=m, nc=m)
  
  for (i in 1:m)
    for (j in 1:m)
      comp[i,j] <- sum((x.group==grlab[i]) & (est.group==grlab[j]))  

  if (by.group)
  {
    er <- vector()
    for (i in 1:m)
      er[i] <- 1-comp[i,i]/rowSums(comp)[i]
    er <- matrix(er, nc=1)
    er <- rbind(er, 1 - sum(diag(comp))/sum(comp)) 
    rownames(er) <- c(as.character(paste(grlab, "(true)")), "Total")
    colnames(er) <- "error"
    
  }
  else 
    er <- 1 - sum(diag(comp))/sum(comp)
  
  comp <- cbind(comp, rowSums(comp))
  comp <- rbind(comp, colSums(comp))

  colnames(comp) <- c(as.character(paste(grlab, "(est.)")), "Total")
  rownames(comp) <- c(as.character(paste(grlab, "(true)")), "Total")

  return(list(cross=comp, error=er))
 
}

###############################################################################
# Computes cross-validated misclassification rates (for use when test data is
# not independent of training data) for KDA
#
# Parameters
# x - training data
# x.group - group variable for x
# y - data values to be classified
# Hs - bandwidth matrices
# prior.prob - prior probabilities
#
# Returns
# List with components
# comp - cross-classification table of groupings - true groups are the rows,
#        estimated groups are the columns
# error - total mis-classification rate
###############################################################################

compare.kda.cv <- function(x, x.group, bw="plugin",
    prior.prob=NULL, Hstart, by.group=FALSE, trace=FALSE, binned=FALSE,
    bgridsize, recompute=FALSE, ...)
{
  ## 1-d
  if (is.vector(x))
  {
    n <- length(x)
    h <- hkda(x, x.group, bw=bw, binned=binned, bgridsize=bgridsize, ...)
    gr <- sort(unique(x.group)) 
    kda.cv.gr <- x.group

    for (i in 1:n)
    {
      h.mod <- h
      ## find group that x[i] belongs to 
      ind <- which(x.group[i]==gr)
      indx <- x.group==gr[ind]
      indx[i] <- FALSE

      if (substr(bw,1,1)=="p")
        h.temp <- hpi(x[indx], binned=binned, bgridsize=bgridsize, ...)

      h.mod[ind] <- h.temp
    
      ## recompute KDA estimate of groups with x[i] excluded
      
      if (trace)
        cat(paste("Processing data item:", i, "\n"))
      
      kda.cv.gr[i] <- kda(x[-i], x.group[-i], hs=h.mod, y=x, prior.prob=prior.prob)[i]
    }

    return(compare(x.group, kda.cv.gr, by.group=by.group)) 
  }


  ## multi-dimensional   
  n <- nrow(x)
  d <- ncol(x)
  
  if (!missing(Hstart))
    H <- Hkda(x, x.group, bw=bw, Hstart=Hstart, binned=binned, bgridsize=bgridsize, ...)
  else
    H <- Hkda(x, x.group, bw=bw, binned=binned, bgridsize=bgridsize, ...)

  ### classify data x using KDA rules based on x itself
  ##kda.group <- kda(x, x.group, Hs=H, y=x, prior.prob=prior.prob)
  ##comp <- compare(x.group, kda.group)
 
  gr <- sort(unique(x.group)) 
  kda.cv.gr <- x.group
  
  for (i in 1:n)
  {
    H.mod <- H
    ### find group that x[i] belongs to 
    ind <- which(x.group[i]==gr)
    indx <- x.group==gr[ind]
    indx[i] <- FALSE

    if (recompute)
    {
      ## compute b/w matrix for that group with x[i] excluded
      if (!missing(Hstart))
      {  
        Hstart.temp <- Hstart[((ind-1)*d+1):(ind*d),]
        
        if (substr(bw,1,1)=="p")
          H.temp <- Hpi(x[indx,], Hstart=Hstart.temp, binned=binned, bgridsize=bgridsize,...)
        else if (substr(bw,1,1)=="s")
          H.temp <- Hscv(x[indx,],  Hstart=Hstart.temp, binned=binned, bgridsize=bgridsize,...)
        else if (substr(bw,1,1)=="l") 
          H.temp <- Hlscv(x[indx,],  Hstart=Hstart.temp, ...)
      }
      else
      {
        if (substr(bw,1,1)=="p")
          H.temp <- Hpi(x[indx,], binned=binned, bgridsize=bgridsize, ...)
        else if (substr(bw,1,1)=="s")
          H.temp <- Hscv(x[indx,], binned=binned, bgridsize=bgridsize, ...)
        else if (substr(bw,1,1)=="l")
          H.temp <- Hlscv(x[indx,], ...) 
      }
      
      H.mod[((ind-1)*d+1):(ind*d),] <- H.temp
    }
    ## recompute KDA estimate of groups with x[i] excluded
      
    if (trace)
      cat(paste("Processing data item:", i, "\n"))

    kda.cv.gr[i] <- kda(x[-i,], x.group[-i], Hs=H.mod, y=x, prior.prob=prior.prob)[i]
  }
  
  return(compare(x.group, kda.cv.gr, by.group=by.group)) 
}

###############################################################################
### Same as compare.kda.cv except uses diagonal b/w matrices
###############################################################################

compare.kda.diag.cv <- function(x, x.group, bw="plugin", prior.prob=NULL,
   by.group=FALSE, trace=FALSE, binned=FALSE, bgridsize, recompute=FALSE, ...)
{
  n <- nrow(x)
  d <- ncol(x)

  H <- Hkda.diag(x, x.group, bw=bw, binned=binned, bgridsize=bgridsize, ...)
  ##kda.group <- kda(x, x.group, Hs=H, y=x, prior.prob=prior.prob)
  ##comp <- compare(x.group, kda.group)
 
  gr <- sort(unique(x.group)) 
  kda.cv.gr <- x.group
  
  for (i in 1:n)
  {
    H.mod <- H

    if (recompute)
    {
      ind <- which(x.group[i]==gr)
      indx <- x.group==gr[ind]
      indx[i] <- FALSE
      if (substr(bw,1,1)=="p")
        H.temp <- Hpi.diag(x[indx,],  binned=binned, bgridsize=bgridsize, ...)
      else if (substr(bw,1,1)=="l")
        H.temp <- Hlscv.diag(x[indx,], binned=binned, bgridsize=bgridsize, ...)
      
      H.mod[((ind-1)*d+1):(ind*d),] <- H.temp
    }
    
    if (trace)
      cat(paste("Processing data item:", i, "\n"))

    kda.cv.gr[i] <- kda(x[-i,], x.group[-i], Hs=H.mod, y=x, prior.prob=prior.prob)[i]  
  }
  
  return(compare(x.group, kda.cv.gr, by.group=by.group)) 
}



###############################################################################
# KDEs of individual densities for KDA - 1- to 3-dim
#
# Parameters
# x - data values
# group - group variable
# Hs - bandwidth matrices
#
# Returns
# List with components (class dade)
# x - list of data values
# eval.points - evaluation points of dnesity estimate
# estimate - list of density estimate
# H - list of bandwidth matrices
##############################################################################

kda.kde <- function(x, x.group, Hs, hs, prior.prob=NULL, gridsize, xmin, xmax, supp=3.7, eval.points=NULL, binned=FALSE, bgridsize)
{
  if (is.vector(x))
  {
    if (missing(gridsize))  gridsize <- 101
    if (missing(bgridsize)) bgridsize <- 401
    fhat.list <- kda.kde.1d(x=x, x.group=x.group, hs=hs, prior.prob=prior.prob, gridsize=gridsize, eval.points=eval.points, supp=supp, binned=binned, bgridsize=bgridsize, xmin=xmin, xmax=xmax)
  }
  else
  {
    if (is.data.frame(x)) x <- as.matrix(x)
    grlab <- sort(unique(x.group))
    m <- length(grlab)
    d <- ncol(x)

    ## find largest bandwidth matrix to initialise grid
    detH <- vector() 
    for (j in 1:m)
      detH[j] <- det(Hs[((j-1)*d+1) : (j*d),])  
    Hmax.ind <- which.max(detH)
    Hmax <- Hs[((Hmax.ind-1)*d+1) : (Hmax.ind*d),]