ksvm.r

这是核学习的一个基础软件包

#**************************************************************#


setMethod("predict", signature(object = "ksvm"),
function (object, newdata, type = "response", coupler = "minpair")
{
  if (missing(newdata))
    return(fit(object))
  ncols <- ncol(xmatrix(object))
  nrows <- nrow(xmatrix(object))
  oldco <- ncols

  if (!is.null(kterms(object)))
    {  
      newdata <- model.matrix(delete.response(kterms(object)), as.data.frame(newdata), na.action = n.action(object))
    }
  else
    newdata  <- if (is.vector(newdata)) t(t(newdata)) else as.matrix(newdata)

  newcols  <- 0
  newnrows <- nrow(newdata)
  newncols <- ncol(newdata)
  newco    <- newncols
    
  if (oldco != newco) stop ("test vector does not match model !")
  p<-0

  if (is.list(scaling(object)))
    newdata[,scaling(object)$scaled] <-
      scale(newdata[,scaling(object)$scaled, drop = FALSE],
            center = scaling(object)$x.scale$"scaled:center",
            scale  = scaling(object)$x.scale$"scaled:scale"
            )

  type <- match.arg(type,c("response","probabilities","votes"))
  
  if(type == "response")
    {
  if(type(object)=="C-classification"||type(object)=="nu-classification")
    {
      predres <- 1:newnrows
      votematrix <- matrix(0,nclass(object),newnrows)
      for(i in 1:(nclass(object)-1))
        {
        jj <- i+1
        for(j in jj:nclass(object))
          {
            p <- p+1
            ret <- kernelMult(kernelf(object),newdata,as.matrix(xmatrix(object)[alphaindex(object)[[p]],]),coeff(object)[[p]]) - b(object)[p]
            votematrix[i,ret>0] <- votematrix[i,ret>0] + 1
            votematrix[j,ret<0] <- votematrix[j,ret<0] + 1
          }
      }
      predres <- sapply(predres, function(x) which.max(votematrix[,x]))
    }

  if(type(object) == "spoc-classification")
    {
      predres <- 1:newnrows
      votematrix <- matrix(0,nclass(object),newnrows)
      for(i in 1:nclass(object))
        votematrix[i,] <- kernelMult(kernelf(object),newdata,xmatrix(object)[alphaindex(object)[[i]],],coeff(object)[[i]])
      predres <- sapply(predres, function(x) which.max(votematrix[,x]))
    }

  if(type(object) == "kbb-classification")
    {
      predres <- 1:newnrows
      votematrix <- matrix(0,nclass(object),newnrows)
      se <- 1:(nclass(object)-1)
      A <- rowSums(alpha(object))
      Aindex <- which(A!=0)
      for(i in 1:nclass(object))
        {
          for(k in se[se<i])
            votematrix[k,] <- votematrix[k,] - (kernelMatrix(kernelf(object),newdata,xmatrix(object)[alphaindex(object)[[k]],]) +1)%*%coeff(object)[[k]]

          votematrix[i,] <- votematrix[i,] + (kernelMatrix(kernelf(object),newdata,xmatrix(object)[Aindex,])+1)%*%A[Aindex]      

          for(k in se[!se<i])
            votematrix[k+1,] <- votematrix[k+1,] - (kernelMatrix(kernelf(object),newdata,xmatrix(object)[alphaindex(object)[[k]],]) +1)%*%coeff(object)[[k]]
        }
      
      predres <- sapply(predres, function(x) which.max(votematrix[,x]))
    }
}

  if(type == "probabilities")
    {
      if(type(object)=="C-classification"||type(object)=="nu-classification")
        {
          binprob <- matrix(0, newnrows, nclass(object)*(nclass(object) - 1)/2)
          for(i in 1:(nclass(object)-1))
            {
              jj <- i+1
              for(j in jj:nclass(object))
                {
                  p <- p+1
                  binprob[,p] <- .probPlatt(as.vector(kernelMult(kernelf(object),newdata,as.matrix(xmatrix(object)[alphaindex(object)[[p]],]),coeff(object)[[p]]) - b(object)[p]))
                }
            }
          multiprob <- couple(binprob, coupler = coupler)
        }
      else
        stop("probability estimates only supported for C-classification and nu-classification")
    }
  
  if(type(object) == "one-classification")
    {
      ret <- kernelMult(kernelf(object),newdata,as.matrix(xmatrix(object)[alphaindex(object),]),coeff(object)) - b(object)
      ret[ret>0]<-1
    }
  
  if(type(object)=="eps-regression"||type(object)=="nu-regression")
    {
      predres <- kernelMult(kernelf(object),newdata,as.matrix(xmatrix(object)[alphaindex(object),]),coeff(object)) - b(object)
    }
  
  if (is.character(lev(object)))
    {
      ##classification & probabilities : return probability matrix
      if(type == "probabilities")
        {
          colnames(multiprob) <- lev(object)
          return(multiprob)
        }
      ##classification & type response: return factors
      if(type == "response")
        return(factor (lev(object)[predres], levels = lev(object)))
      ##classification & votes : return votematrix
      if(type == "votes")
        return(votematrix)
    }
  else if (type(object) == "one-classification")
    ##one-class-classification: return TRUE/FALSE (probabilities ?)
    return(ret == 1)
  else if (!is.null(scaling(object)$y.scale))
    ## return raw values, possibly scaled back
    return(predres * scaling(object)$y.scale$"scaled:scale" + scaling(object)$y.scale$"scaled:center")
  else
    ##else: return raw values
    return(predres)
})

#****************************************************************************************#

setMethod("show","ksvm",
function(object){
  cat("Support Vector Machine object of class \"ksvm\"","\n")
  cat("\n")
   cat(paste("SV type:", type(object), "\n"))
  switch(type(object),
         "C-classification" = cat(paste(" parameter : cost C =",param(object)$C, "\n")),
         "nu-classification" = cat(paste(" parameter : nu =", param(object)$nu, "\n")),
         "one-classification" = cat(paste(" parameter : nu =", param(object)$nu, "\n")),
         "spoc-classification" = cat(paste(" parameter : cost C =",param(object)$C, "\n")),
         "eps-regression" = cat(paste(" parameter : epsilon =",param(object)$epsilon,"\n")),
         "nu-regression" = cat(paste(" parameter : epsilon =", param(object)$epsilon, "nu =", param(object)$nu,"\n"))
         )
  cat("\n")
 show(kernelf(object))
  cat(paste("\nNumber of Support Vectors :", nSV(object),"\n"))
  if(!is.null(fit(object)))
  cat(paste("Training error :", round(error(object),6),"\n"))
  if(cross(object)!= -1)
    cat("Cross validation error :",round(cross(object),6),"\n")
  ##train error & loss
})

setGeneric(".probPlatt", function(deci) standardGeneric(".probPlatt"))
setMethod(".probPlatt",signature(deci="ANY"),
function(deci)
  {
    if (is.matrix(deci))
    deci <- as.vector(deci)
    if (!is.vector(deci))
      stop("input should be matrix or vector")
    ## Create label and count priors
    boolabel <- deci >= 0
    prior1 <- sum(boolabel)
    m <- length(deci)
    prior0 <- m - prior1

    ## set parameters (should be on the interface I guess)
    maxiter <- 100
    minstep <- 1e-10
    sigma <- 1e-3
    eps <- 1e-5

    ## Sigmoid predict function
    .SigmoidPredict <- function(deci, A, B)
      {
        fApB <- deci*A +B
        k<-length(deci)
        ret <- rep(0,k)
        bindex <- fApB >= 0
        ret[bindex] <- exp(-fApB[bindex])/(1 + exp(-fApB[bindex]))
        ret[!bindex] <- 1/(1 + exp(fApB[!bindex]))

        ret
      }
    
    ## Construct target support
    hiTarget <- (prior1 + 1)/(prior1 + 2)
    loTarget <- 1/(prior0 + 2)
    length <- prior1 + prior0
    t <- rep(loTarget, m)
    t[boolabel] <- hiTarget

    ##Initial Point & Initial Fun Value
    A <- 0
    B <- log((prior0 + 1)/(prior1 + 1))
    fval <- 0

    fApB <- deci*A + B
    bindex <- fApB >= 0

    fval <- sum(t[bindex]*fApB[bindex] + log(1 + exp(-fApB[bindex])))
    fval <- fval + sum((t[!bindex] - 1)*fApB[!bindex] + log(1+exp(fApB[!bindex])))


    for (it in 1:maxiter)
      {
        h11 <- h22 <- sigma
        h21 <- g1 <- g2 <- 0
        fApB <- deci*A + B
        
        for (l in 1:2)
          {
          if (l == 1)
            {
              bindex <- fApB >= 0
              p <- exp(-fApB[bindex])/(1 + exp(-fApB[bindex]))
              q <- 1/(1+exp(-fApB[bindex]))
            }
          if (l == 2)
            {
              bindex <- fApB < 0
              p <- 1/(1 + exp(fApB[bindex]))
              q <- exp(fApB[bindex])/(1 + exp(fApB[bindex]))
            }
          
          d2 <- p*q
          h11 <- h11 + sum(d2*deci[bindex]^2)
          h22 <- h22 + sum(d2)
          h21 <- h21 + sum(deci[bindex]*d2)
          d1 <- t[bindex] - p
          g1 <- g1 + sum(deci[bindex]*d1)
          g2 <- g2 + sum(d1)
        }
        ## Stopping Criteria
        if (abs(g1) < eps && abs(g2) < eps)
          break

        ## Finding Newton Direction -inv(t(H))%*%g
        det <- h11*h22 - h21^2
        dA <- -(h22*g1 - h21*g2) / det
        dB <- -(-h21*g1 + h11*g2) / det
        gd <- g1*dA + g2*dB

        ## Line Search

        stepsize <- 1

        while(stepsize >= minstep)
          {
            newA <- A + stepsize * dA
            newB <- B + stepsize * dB

            ## New function value
            newf <- 0
            fApB <- deci * newA + newB
            bindex <- fApB >= 0 
            newf <- sum(t[bindex] * fApB[bindex] + log(1 + exp(-fApB[bindex])))
            newf <- newf + sum((t[!bindex] - 1)*fApB[!bindex] + log(1 + exp(fApB[!bindex])))

           ## Check decrease
            if (newf < (fval + 0.0001 * stepsize * gd))
              {
                A <- newA
                B <- newB
                fval <- newf
                break
              }
            else
              stepsize <- stepsize/2
          }
        if (stepsize < minstep)
          {
            cat("line search fails", A, B, g1, g2, dA, dB, gd)
            ret <- .SigmoidPredict(deci, A, B)
            return(ret)
          }
      }
    if(it >= maxiter -1)
      cat("maximum number of iterations reached",g1,g2)

    ret <- .SigmoidPredict(deci, A, B)
    return(ret)
  })