Home: http://manmustbecool.github.io/MyWiki/

svm

0.1 Support Vector Machine (SVM)

• e1071

tune http://rgm2.lab.nig.ac.jp/RGM2/func.php?rd_id=e1071:tune

svm http://rgm2.lab.nig.ac.jp/RGM2/func.php?rd_id=e1071:svm

http://stackoverflow.com/questions/7390173/svm-equations-from-e1071-r-package

Recursive Feature Extraction (SVM-RFE)

http://stats.stackexchange.com/questions/10676/using-an-svm-for-feature-selection

0.2 SVM example 1

library(e1071)

day = c(0,1,2,3,4,5,6)
weather = c(1,0,0,0,0,0,0)
happy = factor(c(T,F,F,F,F,F,F))
d = data.frame(day=day, weather=weather, happy=happy)
d <- rbind(d, d)

model = svm(happy ~ day + weather, data = d)
plot(model, d)

# tune the paramaters 
tuned <- tune.svm( happy ~ ., data = d, 
                    gamma = 10^(-2:0), cost = 10^(0:2), kernel = 'radial')
summary(tuned)

## 
## Parameter tuning of 'svm':
## 
## - sampling method: 10-fold cross validation 
## 
## - best parameters:
##  gamma cost
##    0.1    1
## 
## - best performance: 0 
## 
## - Detailed performance results:
##   gamma cost error dispersion
## 1  0.01    1  0.15  0.3374743
## 2  0.10    1  0.00  0.0000000
## 3  1.00    1  0.00  0.0000000
## 4  0.01   10  0.00  0.0000000
## 5  0.10   10  0.00  0.0000000
## 6  1.00   10  0.00  0.0000000
## 7  0.01  100  0.00  0.0000000
## 8  0.10  100  0.00  0.0000000
## 9  1.00  100  0.00  0.0000000

model <- svm(happy ~ ., data = d, kernel= 'radial', 
                   gamma = tuned$best.parameters$gamma, 
                   cost= tuned$best.parameters$cost,
                   scale=FALSE, cross = 10)
print(model)

## 
## Call:
## svm(formula = happy ~ ., data = d, kernel = "radial", gamma = tuned$best.parameters$gamma, 
##     cost = tuned$best.parameters$cost, cross = 10, scale = FALSE)
## 
## 
## Parameters:
##    SVM-Type:  C-classification 
##  SVM-Kernel:  radial 
##        cost:  1 
##       gamma:  0.1 
## 
## Number of Support Vectors:  5

plot(model, d)

# support vector
SV <- d[c(model$index),]
SV

##    day weather happy
## 1    0       1  TRUE
## 8    0       1  TRUE
## 2    1       0 FALSE
## 9    1       0 FALSE
## 14   6       0 FALSE

0.3 SVM example 2

# install.packages('kernlab')
# install.packages('ROCR')

## Generate data ##
##---------------##

# number of data points
n <- 150
# dimension
p <- 2

sigma <- 2  # variance of the distribution
meanpos <- 0 # centre of the distribution of positive examples
meanneg <- 3 # centre of the distribution of negative examples
npos <- round(n/2) # number of positive examples
nneg <- n-npos # number of negative examples

# Generate the positive and negative examples
xpos <- matrix(rnorm(npos*p,mean=meanpos,sd=sigma),npos,p)
xneg <- matrix(rnorm(nneg*p,mean=meanneg,sd=sigma),npos,p)
x <- rbind(xpos,xneg)

# Generate the labels
y <- matrix(c(rep(1,npos),rep(-1,nneg)))

# Visualize the data
plot(x,col=ifelse(y>0,1,2))
legend("topleft",c('Positive','Negative'),col=seq(2),pch=1,text.col=seq(2))

## Prepare a training and a test set ##
ntrain <- round(n*0.8) # number of training examples
tindex <- sample(n,ntrain) # indices of training samples
xtrain <- x[tindex,]
xtest <- x[-tindex,]
ytrain <- y[tindex]
ytest <- y[-tindex]
istrain=rep(0,n)
istrain[tindex]=1

# Visualize
plot(x,col=ifelse(y>0,1,2),pch=ifelse(istrain==1,1,2))
legend("topleft",c('Positive Train','Positive Test','Negative Train','Negative Test'),col=c(1,1,2,2),pch=c(1,2,1,2),text.col=c(1,1,2,2))

## Train a linear SVM with C=1 ##
##-----------------------------##

# load the kernlab package

library(kernlab)

# train the SVM
svp <- ksvm(xtrain,ytrain,type="C-svc",kernel='vanilladot',C=1,scaled=c())

##  Setting default kernel parameters

# General summary
svp

## Support Vector Machine object of class "ksvm" 
## 
## SV type: C-svc  (classification) 
##  parameter : cost C = 1 
## 
## Linear (vanilla) kernel function. 
## 
## Number of Support Vectors : 36 
## 
## Objective Function Value : -33.72 
## Training error : 0.091667

# Attributes that you can access
attributes(svp)

## $param
## $param$C
## [1] 1
## 
## 
## $scaling
## `\001NULL\001`
## 
## $coef
## $coef[[1]]
##  [1]  1.00000000  1.00000000 -1.00000000  1.00000000 -0.04499051
##  [6]  1.00000000  0.74970416 -1.00000000 -1.00000000 -1.00000000
## [11]  0.29528635  1.00000000 -1.00000000  1.00000000  1.00000000
## [16]  1.00000000 -1.00000000 -1.00000000 -1.00000000 -1.00000000
## [21]  1.00000000 -1.00000000  1.00000000 -1.00000000  1.00000000
## [26] -1.00000000 -1.00000000 -1.00000000  1.00000000 -1.00000000
## [31]  1.00000000 -1.00000000  1.00000000  1.00000000 -1.00000000
## [36]  1.00000000
## 
## 
## $alphaindex
## $alphaindex[[1]]
##  [1]   3   6  12  13  18  26  28  29  31  38  41  42  45  58  61  62  64
## [18]  65  67  68  69  73  76  77  78  82  85  93  96 106 107 109 111 113
## [35] 116 119
## 
## 
## $b
## [1] -1.509243
## 
## $obj
## [1] -33.71998
## 
## $SVindex
##  [1]   3   6  12  13  18  26  28  29  31  38  41  42  45  58  61  62  64
## [18]  65  67  68  69  73  76  77  78  82  85  93  96 106 107 109 111 113
## [35] 116 119
## 
## $nSV
## [1] 36
## 
## $prior
## $prior[[1]]
## $prior[[1]]$prior1
## [1] 59
## 
## $prior[[1]]$prior0
## [1] 61
## 
## 
## 
## $prob.model
## $prob.model[[1]]
## NULL
## 
## 
## $alpha
## $alpha[[1]]
##  [1] 1.00000000 1.00000000 1.00000000 1.00000000 0.04499051 1.00000000
##  [7] 0.74970416 1.00000000 1.00000000 1.00000000 0.29528635 1.00000000
## [13] 1.00000000 1.00000000 1.00000000 1.00000000 1.00000000 1.00000000
## [19] 1.00000000 1.00000000 1.00000000 1.00000000 1.00000000 1.00000000
## [25] 1.00000000 1.00000000 1.00000000 1.00000000 1.00000000 1.00000000
## [31] 1.00000000 1.00000000 1.00000000 1.00000000 1.00000000 1.00000000
## 
## 
## $type
## [1] "C-svc"
## 
## $kernelf
## function (x, y = NULL) 
## {
##     if (!is(x, "vector")) 
##         stop("x must be a vector")
##     if (!is(y, "vector") && !is.null(y)) 
##         stop("y must be a vector")
##     if (is(x, "vector") && is.null(y)) {
##         crossprod(x)
##     }
##     if (is(x, "vector") && is(y, "vector")) {
##         if (!length(x) == length(y)) 
##             stop("number of dimension must be the same on both data points")
##         crossprod(x, y)
##     }
## }
## <environment: 0x000000000bbd4af0>
## attr(,"kpar")
## list()
## attr(,"class")
## [1] "vanillakernel"
## attr(,"class")attr(,"package")
## [1] "kernlab"
## 
## $kpar
## list()
## 
## $xmatrix
## $xmatrix[[1]]
##             X1          X2
## 3    3.8059092  1.16077313
## 6    1.3854210  3.57695610
## 12   3.3275995  0.08856863
## 13   2.5667411 -0.89491832
## 18  -3.1786713  5.61838461
## 26   0.3568473  2.52519359
## 28  -2.1029615  2.11112018
## 29   1.4481414  2.05639528
## 31   2.3518444  0.95942526
## 38   1.1989913  1.87183726
## 41   2.7775676 -1.13749692
## 42  -0.2690097  0.91360820
## 45   0.4490090  2.38097913
## 58  -0.7820423  2.55831075
## 61   0.3224972  1.35262231
## 62   3.6006022 -0.52254510
## 64  -1.7522648  0.56921490
## 65  -0.8353704  2.37649451
## 67   1.5164193  1.19274595
## 68   0.5751521 -0.56200090
## 69   2.1286079 -0.39315328
## 73  -0.5640993  1.18659336
## 76  -0.5161680  1.18723482
## 77   1.0538383  2.58332550
## 78   2.0032955 -0.30237108
## 82   4.5301881 -0.15709788
## 85  -2.8157637  2.28825347
## 93  -0.9999351  3.46626875
## 96  -3.3350431  4.17067419
## 106  1.7227528 -0.24388748
## 107 -0.3344997  1.33275655
## 109 -0.3867338 -0.18492715
## 111  1.3831066  0.04381767
## 113 -0.6644324  3.10944928
## 116  2.1851913  1.78297066
## 119  1.4900020  0.12652469
## 
## 
## $ymatrix
##   [1] -1  1  1  1 -1  1 -1 -1  1  1 -1 -1  1  1 -1 -1  1 -1 -1  1  1 -1  1
##  [24]  1  1  1 -1  1 -1 -1 -1 -1 -1  1 -1 -1 -1 -1  1  1  1  1  1  1 -1 -1
##  [47] -1 -1 -1 -1 -1 -1  1 -1 -1  1  1  1  1 -1  1  1 -1 -1 -1  1 -1 -1  1
##  [70] -1  1 -1 -1  1  1  1 -1  1 -1  1  1 -1  1 -1 -1 -1 -1 -1  1 -1  1  1
##  [93] -1 -1 -1  1 -1  1 -1  1  1  1 -1  1  1 -1  1  1 -1  1  1 -1  1  1  1
## [116] -1 -1  1  1 -1
## 
## $fitted
##   [1] -1  1 -1  1 -1 -1 -1 -1  1  1 -1 -1  1  1 -1 -1  1 -1 -1  1  1 -1  1
##  [24]  1  1 -1 -1  1 -1 -1 -1 -1 -1  1 -1 -1 -1 -1  1  1  1  1  1  1 -1 -1
##  [47] -1 -1 -1 -1 -1 -1  1 -1 -1  1  1  1  1 -1  1  1 -1  1  1  1 -1  1  1
##  [70] -1  1 -1  1  1  1  1 -1  1 -1  1  1 -1  1 -1  1 -1 -1 -1  1 -1  1  1
##  [93] -1 -1 -1  1 -1  1 -1  1  1  1 -1  1  1  1  1  1  1  1  1 -1 -1  1  1
## [116] -1 -1  1  1 -1
## 
## $lev
## [1] -1  1
## 
## $nclass
## [1] 2
## 
## $error
## [1] 0.09166667
## 
## $cross
## [1] -1
## 
## $n.action
## function (object, ...) 
## UseMethod("na.omit")
## <bytecode: 0x00000000078cdf78>
## <environment: namespace:stats>
## 
## $terms
## `\001NULL\001`
## 
## $kcall
## .local(x = x, y = ..1, scaled = ..5, type = "C-svc", kernel = "vanilladot", 
##     C = 1)
## 
## $class
## [1] "ksvm"
## attr(,"package")
## [1] "kernlab"

# For example, the support vectors
alpha(svp)

## [[1]]
##  [1] 1.00000000 1.00000000 1.00000000 1.00000000 0.04499051 1.00000000
##  [7] 0.74970416 1.00000000 1.00000000 1.00000000 0.29528635 1.00000000
## [13] 1.00000000 1.00000000 1.00000000 1.00000000 1.00000000 1.00000000
## [19] 1.00000000 1.00000000 1.00000000 1.00000000 1.00000000 1.00000000
## [25] 1.00000000 1.00000000 1.00000000 1.00000000 1.00000000 1.00000000
## [31] 1.00000000 1.00000000 1.00000000 1.00000000 1.00000000 1.00000000

alphaindex(svp)

## [[1]]
##  [1]   3   6  12  13  18  26  28  29  31  38  41  42  45  58  61  62  64
## [18]  65  67  68  69  73  76  77  78  82  85  93  96 106 107 109 111 113
## [35] 116 119

b(svp)

## [1] -1.509243

# Use the built-in function to pretty-plot the classifier
plot(svp,data=xtrain)

# Predict labels on test
ypred = predict(svp,xtest)
# Confusion table
table(ytest,ypred)

##      ypred
## ytest -1  1
##    -1 12  2
##    1   4 12

# Compute accuracy
sum(ypred==ytest)/length(ytest)

## [1] 0.8

# Compute at the prediction scores
ypredscore = predict(svp,xtest,type="decision")

# Check that the predicted labels are the signs of the scores
table(ypredscore > 0,ypred)

##        ypred
##         -1  1
##   FALSE 16  0
##   TRUE   0 14

# Package to compute ROC curve, precision-recall etc...
library(ROCR)

## Warning: package 'ROCR' was built under R version 3.3.3

## Loading required package: gplots

## Warning: package 'gplots' was built under R version 3.3.2

## 
## Attaching package: 'gplots'

## The following object is masked from 'package:stats':
## 
##     lowess

pred <- prediction(ypredscore,ytest)

# Plot ROC curve
perf <- performance(pred, measure = "tpr", x.measure = "fpr") 
plot(perf)

# Plot precision/recall curve
perf <- performance(pred, measure = "prec", x.measure = "rec") 
plot(perf)

# Plot accuracy as function of threshold
perf <- performance(pred, measure = "acc") 
plot(perf)

## Use cross-validation ##

cv.folds <-
function(n, folds = 10)
{
    split(sample(1:n), rep(1:folds, length = n))
}

cvpred.ksvm <- function(x,y,folds=3,predtype="response",...)
## Return a vector of predictions by cross-validation
## 'predtype' should be one of response (by default), decision or probabilities, depending the prediction we want (SVM label, score or probability, see predict.ksvm())
## Additional parameters are passed to ksvm() to train the SVM
{
   n <- length(y)
    ypred <- numeric(n)
    s <- cv.folds(n,folds)
    for (i in seq(folds)) {
        m <- ksvm(x[-s[[i]],],y[-s[[i]]],...)
        ypred[s[[i]]] <- predict(m,x[s[[i]],],type=predtype)
        }
    invisible(ypred)
 }

##----------------------##

# We compute the prediction vector by cross-validation
k=5  
ypredscorecv <- cvpred.ksvm(x,y,folds=k, type="C-svc", kernel='vanilladot', C=1, scaled=c(), predtype="decision")

##  Setting default kernel parameters  
##  Setting default kernel parameters  
##  Setting default kernel parameters  
##  Setting default kernel parameters  
##  Setting default kernel parameters

# Check the performance
print(table(ypredscorecv > 0,y))

##        y
##         -1  1
##   FALSE 64  9
##   TRUE  11 66

pred <- prediction(ypredscorecv,y)
perf <- performance(pred, measure = "tpr", x.measure = "fpr") 
plot(perf)

# Estimate the CV error with ksvm directly, and compare
svp <- ksvm(x,y,type="C-svc",kernel='vanilladot',C=1,scaled=c(),cross=5)

##  Setting default kernel parameters

print(cross(svp))

## [1] 0.12

print(1-sum((ypredscorecv>0)==(y==1))/n)

## [1] 0.1333333

code is from WWW.

^Home Page^