Hi,
A. In a nutshell:
The training error, obtained as "error (ret)", from the return value
of a ksvm () call for a eps-svr model is (likely) being computed
wrongly. "nu-svr" and "eps-bsvr" suffer from this as well.
I am attaching three files: (1) ksvm.R from the the kernlab package,
un-edited, (2) ksvm_eps-svr.txt: (for easier reading) containing only
eps-svr pertinent blocks of code from ksvm.R with line numbers
prefixed to each line (I hope no licenses are being violated!) and (3)
trace_output_concise.txt: relevant output from a trace () call for the
example in C (my commands there are marked by "NOTE THIS" strings).
B. In detail:
(the line numbers refer to either ksvm.R or the "prefixed" numbers in
ksvm_eps-svr.txt)
(you may also refer to the trace_output in parallel)
1. the given response vector, y, is standardized at line 143:
142 if
(is.numeric(y)&&(type(ret)!="C-svc"&&type(ret)!="nu-svc"&&type(ret)!="C-bsvc"&&type(ret)!="spoc-svc"&&type(ret)!="kbb-svc"))
{
143 y <- scale(y)
144 y.scale <- attributes(y)[c("scaled:center","scaled:scale")]
145 y <- as.vector(y)
146 }
2. fitted response, fitted(ret), is obtained at lines 826,827:
826 fitted(ret) <- if (fit)
827 predict(ret, x) else NULL
(note: during the fitting process the return value of predict(ret,x)
is _NOT_ scaled BACK to the original units, see lines 2832-2838)
3. fitted response is now scaled BACK to the original units of "y" at
lines 839,840. So, when assigning "error(ret) <-" at line 844,
fitted(ret) is in the original units but y is in standardized units:
828
829 if(any(scaled))
830 scaling(ret) <- list(scaled = scaled, x.scale = x.scale,
y.scale = y.scale)
831
832 if (fit){
...
837
if(type(ret)=="eps-svr"||type(ret)=="nu-svr"||type(ret)=="eps-bsvr"){
838 if (!is.null(scaling(ret)$y.scale)){
839 scal <- scaling(ret)$y.scale$"scaled:scale"
840 fitted(ret) <- fitted(ret) *
scaling(ret)$y.scale$"scaled:scale" +
scaling(ret)$y.scale$"scaled:center"
841 }
842 else
843 scal <- 1
844 error(ret) <- drop((scal^2)*crossprod(fitted(ret) - y)/m)
845 }
846 }
C. Finally, an example (taken from ?ksvm):
require (kernlab)
seed (1234)
x <- seq(-20,20,0.1); x <- x[x != 0]
y <- sin(x)/x + rnorm(400,sd=0.03)
regm <- ksvm(x,y,epsilon=0.01,kpar=list(sigma=16),cross=3)
te <- crossprod (fitted(regm)-y)/400
s <- (scaling(regm)$y.scale[["scaled:scale"]])^2
error (regm) # 0.03891344
te # 0.0008958718
te * s # 6.37252e-05
te / s # 0.01259449
These numbers can also be seen in the trace_output_concise.txt. In
particular, compare the output marked by "** NOTE THIS", "?? NOTE
THIS", and "++ NOTE THIS" there.
Basically, compute the error either before scaling back "fitted (ret)"
or scale back y as well!!
Can anyone confirm / refute / agree / disagree?
Thanks,
--
Prasenjit
1 ## Support Vector Machines
2 ## author : alexandros karatzoglou
3 ## updated : 08.02.06
4
5 setGeneric("ksvm", function(x, ...) standardGeneric("ksvm"))
6 setMethod("ksvm",signature(x="formula"),
7 function (x, data=NULL, ..., subset, na.action = na.omit, scaled =
TRUE){
...
38 })
39
...
47 setMethod("ksvm",signature(x="matrix"),
48 function (x,
49 y = NULL,
50 scaled = TRUE,
51 type = NULL,
52 kernel = "rbfdot",
53 kpar = "automatic",
54 C = 1,
55 nu = 0.2,
56 epsilon = 0.1,
57 prob.model = FALSE,
58 class.weights = NULL,
59 cross = 0,
60 fit = TRUE,
61 cache = 40,
62 tol = 0.001,
63 shrinking = TRUE,
64 ...
65 ,subset
66 ,na.action = na.omit)
67 {
...
74 sparse <- FALSE
75
76 if(is.character(kernel)){
77 kernel <-
match.arg(kernel,c("rbfdot","polydot","tanhdot","vanilladot","laplacedot","besseldot","anovadot","splinedot","matrix"))
78
79 if(kernel == "matrix")
80 if(dim(x)[1]==dim(x)[2])
81 return(ksvm(as.kernelMatrix(x), y = y, type = type, C = C, nu =
nu, epsilon = epsilon, prob.model = prob.model, class.weights = class.weights,
cross = cross, fit = fit, cache = cache, tol = tol, shrinking = shrinking, ...))
82 else
83 stop(" kernel matrix not square!")
84
85 if(is.character(kpar))
86 if((kernel == "tanhdot" || kernel == "vanilladot" || kernel ==
"polydot"|| kernel == "besseldot" || kernel== "anovadot"|| kernel=="splinedot")
&& kpar=="automatic" )
87 {
88 cat (" Setting default kernel parameters ","\n")
89 kpar <- list()
90 }
91 }
92
93 ## subsetting and na-handling for matrices
94 ret <- new("ksvm")
95 if (!missing(subset)) x <- x[subset,]
96 if (is.null(y))
97 x <- na.action(x)
98 else {
99 df <- na.action(data.frame(y, x))
100 y <- df[,1]
101 x <- as.matrix(df[,-1])
102 }
103 n.action(ret) <- na.action
104
105 if (is.null(type)) type(ret) <- if (is.null(y)) "one-svc" else if
(is.factor(y)) "C-svc" else "eps-svr"
106
107 if(!is.null(type))
108 type(ret) <- match.arg(type,c("C-svc",
109 "nu-svc",
110 "kbb-svc",
111 "spoc-svc",
112 "C-bsvc",
113 "one-svc",
114 "eps-svr",
115 "eps-bsvr",
116 "nu-svr"))
124
125 x.scale <- y.scale <- NULL
126 ## scaling
127 if (length(scaled) == 1)
128 scaled <- rep(scaled, ncol(x))
129 if (any(scaled)) {
130 co <- !apply(x[,scaled, drop = FALSE], 2, var)
131 if (any(co)) {
132 scaled <- rep(FALSE, ncol(x))
133 warning(paste("Variable(s)",
134 paste("`",colnames(x[,scaled, drop = FALSE])[co],
135 "'", sep="", collapse=" and "),
136 "constant. Cannot scale data.")
137 )
138 } else {
139 xtmp <- scale(x[,scaled])
140 x[,scaled] <- xtmp
141 x.scale <- attributes(xtmp)[c("scaled:center","scaled:scale")]
142 if
(is.numeric(y)&&(type(ret)!="C-svc"&&type(ret)!="nu-svc"&&type(ret)!="C-bsvc"&&type(ret)!="spoc-svc"&&type(ret)!="kbb-svc"))
{
143 y <- scale(y)
144 y.scale <- attributes(y)[c("scaled:center","scaled:scale")]
145 y <- as.vector(y)
146 }
147 }
148 }
149 ncols <- ncol(x)
150 m <- nrows <- nrow(x)
151
152 if (!is.function(kernel))
153 if (!is.list(kpar)&&is.character(kpar)&&(class(kernel)=="rbfkernel"
|| class(kernel) =="laplacedot" || kernel == "laplacedot"|| kernel=="rbfdot")){
154 kp <- match.arg(kpar,"automatic")
155 if(kp=="automatic")
156 kpar <- list(sigma=mean(sigest(x,scaled=FALSE)[c(1,3)]))
157 cat("Using automatic sigma estimation (sigest) for RBF or laplace
kernel","\n")
158
159 }
160 if(!is(kernel,"kernel"))
161 {
162 if(is(kernel,"function")) kernel <- deparse(substitute(kernel))
163 kernel <- do.call(kernel, kpar)
164 }
165
166 if(!is(kernel,"kernel")) stop("kernel must inherit from class
`kernel'")
167
168 if (!is(y,"vector") && !is.factor (y) & is(y,"matrix") &
!(type(ret)=="one-svc")) stop("y must be a vector or a factor.")
169
170 if(!(type(ret)=="one-svc"))
171 if(is(y,"vector") | is(y,"factor") ) ym <- length(y) else
if(is(y,"matrix")) ym <- dim(y)[1] else stop("y must be a matrix or a vector")
172
173 if ((type(ret) != "one-svc") && ym != m) stop("x and y don't
match.")
174
175 if(nu > 1|| nu <0) stop("nu must be between 0 an 1.")
176
177 weightlabels <- NULL
178 nweights <- 0
179 weight <- 0
180 wl <- 0
181 ## in case of classification: transform factors into integers
182 if (type(ret) == "one-svc") # one class classification --> set dummy
183 y <- 1
184 else
185 if (is.factor(y)) {
186 lev(ret) <- levels (y)
187 y <- as.integer (y)
188 if (!is.null(class.weights)) {
189 weightlabels <- match (names(class.weights),lev(ret))
190 if (any(is.na(weightlabels)))
191 stop ("At least one level name is missing or misspelled.")
192 }
193 }
194 else {
...
198 if (type(ret) != "eps-svr" || type(ret) != "nu-svr"||
type(ret)!="eps-bsvr")
199 lev(ret) <- sort(unique (y))
200 }
201 ## initialize
202 nclass(ret) <- length (unique(y))
203 p <- 0
204 K <- 0
205 svindex <- problem <- NULL
206 sigma <- 0.1
207 degree <- offset <- scale <- 1
208
209 switch(is(kernel)[1],
210 "rbfkernel" =
211 {
212 sigma <- kpar(kernel)$sigma
213 ktype <- 2
214 },
215 "tanhkernel" =
216 {
217 sigma <- kpar(kernel)$scale
218 offset <- kpar(kernel)$offset
219 ktype <- 3
220 },
221 "polykernel" =
222 {
223 degree <- kpar(kernel)$degree
224 sigma <- kpar(kernel)$scale
225 offset <- kpar(kernel)$offset
226 ktype <- 1
227 },
228 "vanillakernel" =
229 {
230 ktype <- 0
231 },
232 "laplacekernel" =
233 {
234 ktype <- 5
235 sigma <- kpar(kernel)$sigma
236 },
237 "besselkernel" =
238 {
239 ktype <- 6
240 sigma <- kpar(kernel)$sigma
241 degree <- kpar(kernel)$order
242 offset <- kpar(kernel)$degree
243 },
244 "anovakernel" =
245 {
246 ktype <- 7
247 sigma <- kpar(kernel)$sigma
248 degree <- kpar(kernel)$degree
249 },
250 "splinekernel" =
251 {
252 ktype <- 8
253 },
254 {
255 ktype <- 4
256 }
257 )
258 prior(ret) <- list(NULL)
259
260 ## C classification
261 if(type(ret) == "C-svc"){
...
353 }
354
355 ## nu classification
356 if(type(ret) == "nu-svc"){
...
441 }
442
443 ## Bound constraint C classification
444 if(type(ret) == "C-bsvc"){
...
538 }
539
540 ## SPOC multiclass classification
541 if(type(ret) =="spoc-svc")
542 {
...
595 }
596
597 ## KBB multiclass classification
598 if(type(ret) =="kbb-svc")
599 {
...
648 }
649
650 ## Novelty detection
651 if(type(ret) =="one-svc")
652 {
...
689 }
690
691 ## epsilon regression
692 if(type(ret) =="eps-svr")
693 {
694 if(ktype==4)
695 K <- kernelMatrix(kernel,x)
696
697 resv <- .Call("smo_optim",
698 as.double(t(x)),
699 as.integer(nrow(x)),
700 as.integer(ncol(x)),
701 as.double(y),
702 as.double(K),
703 as.integer(if (sparse) x...@ia else 0),
704 as.integer(if (sparse) x...@ja else 0),
705 as.integer(sparse),
706 as.double(matrix(rep(-1,m))),
707 as.integer(ktype),
708 as.integer(3),
709 as.double(C),
710 as.double(nu),
711 as.double(epsilon),
712 as.double(sigma),
713 as.integer(degree),
714 as.double(offset),
715 as.integer(0), #weightlabl.
716 as.double(0),
717 as.integer(0),
718 as.double(cache),
719 as.double(tol),
720 as.integer(shrinking),
721 PACKAGE="kernlab")
722 tmpres <- resv[c(-(m+1),-(m+2))]
723 alpha(ret) <- coef(ret) <- tmpres[tmpres != 0]
724 svindex <- alphaindex(ret) <- which(tmpres != 0)
725 xmatrix(ret) <- x[svindex, ,drop=FALSE]
726 b(ret) <- resv[(m+1)]
727 obj(ret) <- resv[(m+2)]
728 param(ret)$epsilon <- epsilon
729 param(ret)$C <- C
730 }
731
732 ## nu regression
733 if(type(ret) =="nu-svr")
734 {
735 if(ktype==4)
736 K <- kernelMatrix(kernel,x)
737
738 resv <- .Call("smo_optim",
739 as.double(t(x)),
740 as.integer(nrow(x)),
741 as.integer(ncol(x)),
742 as.double(y),
743 as.double(K),
744 as.integer(if (sparse) x...@ia else 0),
745 as.integer(if (sparse) x...@ja else 0),
746 as.integer(sparse),
747 as.double(matrix(rep(-1,m))),
748 as.integer(ktype),
749 as.integer(4),
750 as.double(C),
751 as.double(nu),
752 as.double(epsilon),
753 as.double(sigma),
754 as.integer(degree),
755 as.double(offset),
756 as.integer(0),
757 as.double(0),
758 as.integer(0),
759 as.double(cache),
760 as.double(tol),
761 as.integer(shrinking),
762 PACKAGE="kernlab")
763 tmpres <- resv[c(-(m+1),-(m+2))]
764 alpha(ret) <- coef(ret) <- tmpres[tmpres!=0]
765 svindex <- alphaindex(ret) <- which(tmpres != 0)
766 xmatrix(ret) <- x[svindex,,drop=FALSE]
767 b(ret) <- resv[(m+1)]
768 obj(ret) <- resv[(m+2)]
769 param(ret)$epsilon <- epsilon
770 param(ret)$nu <- nu
771 }
772
773 ## bound constraint eps regression
774 if(type(ret) =="eps-bsvr")
775 {
776 if(ktype==4)
777 K <- kernelMatrix(kernel,x)
778
779 resv <- .Call("tron_optim",
780 as.double(t(x)),
781 as.integer(nrow(x)),
782 as.integer(ncol(x)),
783 as.double(y),
784 as.double(K),
785 as.integer(if (sparse) x...@ia else 0),
786 as.integer(if (sparse) x...@ja else 0),
787 as.integer(sparse),
788 as.integer(2),
789 as.integer(0),
790 as.integer(ktype),
791 as.integer(6),
792 as.double(C),
793 as.double(epsilon),
794 as.double(sigma),
795 as.integer(degree),
796 as.double(offset),
797 as.double(1), #Cbegin
798 as.double(2), #Cstep
799 as.integer(0), #weightlabl.
800 as.double(0),
801 as.integer(0),
802 as.double(0),
803 as.double(cache),
804 as.double(tol),
805 as.integer(10), #qpsize
806 as.integer(shrinking),
807 PACKAGE="kernlab")
808 tmpres <- resv[-(m + 1)]
809 alpha(ret) <- coef(ret) <- tmpres[tmpres!=0]
810 svindex <- alphaindex(ret) <- which(tmpres != 0)
811 xmatrix(ret) <- x[svindex,,drop=FALSE]
812 b(ret) <- -sum(alpha(ret))
813 obj(ret) <- resv[(m + 1)]
814 param(ret)$epsilon <- epsilon
815 param(ret)$C <- C
816 }
817
818
819 kcall(ret) <- match.call()
820 kernelf(ret) <- kernel
821 ymatrix(ret) <- y
822 SVindex(ret) <- sort(unique(svindex),method="quick")
823 nSV(ret) <- length(unique(svindex))
824 if(nSV(ret)==0)
825 stop("No Support Vectors found. You may want to change your
parameters")
826 fitted(ret) <- if (fit)
827 predict(ret, x) else NULL
828
829 if(any(scaled))
830 scaling(ret) <- list(scaled = scaled, x.scale = x.scale, y.scale =
y.scale)
831
832 if (fit){
...
837
if(type(ret)=="eps-svr"||type(ret)=="nu-svr"||type(ret)=="eps-bsvr"){
838 if (!is.null(scaling(ret)$y.scale)){
839 scal <- scaling(ret)$y.scale$"scaled:scale"
840 fitted(ret) <- fitted(ret) *
scaling(ret)$y.scale$"scaled:scale" + scaling(ret)$y.scale$"scaled:center"
841 }
842 else
843 scal <- 1
844 error(ret) <- drop((scal^2)*crossprod(fitted(ret) - y)/m)
845 }
846 }
847
848 cross(ret) <- -1
849 if(cross == 1)
850 cat("\n","cross should be >1 no cross-validation done!","\n","\n")
851 else if (cross > 1)
852 {
853
854 cerror <- 0
855 suppressWarnings(vgr<-split(sample(1:m,m),1:cross))
856 for(i in 1:cross)
857 {
858
859 cind <-
unsplit(vgr[-i],factor(rep((1:cross)[-i],unlist(lapply(vgr[-i],length)))))
860
if(type(ret)=="C-svc"||type(ret)=="nu-svc"||type(ret)=="spoc-svc"||type(ret)=="kbb-svc"||type(ret)=="C-bsvc")
861 {
...
868 }
869 if(type(ret)=="one-svc")
870 {
...
874 }
875
876
if(type(ret)=="eps-svr"||type(ret)=="nu-svr"||type(ret)=="eps-bsvr")
877 {
878 cret <-
ksvm(x[cind,],y[cind],type=type(ret),kernel=kernel,kpar =
NULL,C=C,nu=nu,epsilon=epsilon,tol=tol,scaled=FALSE, cross = 0, fit = FALSE,
cache = cache, prob.model = FALSE)
879 cres <- predict(cret, x[vgr[[i]],,drop=FALSE])
880 if (!is.null(scaling(ret)$y.scale))
881 scal <- scaling(ret)$y.scale$"scaled:scale"
882 else
883 scal <- 1
884 cerror <- drop((scal^2)*crossprod(cres - y[vgr[[i]]])/m)
+ cerror
885 }
886 }
887 cross(ret) <- cerror
888 }
889
890 prob.model(ret) <- list(NULL)
891
892 if(prob.model)
893 {
894 if(type(ret)=="C-svc"||type(ret)=="nu-svc"||type(ret)=="C-bsvc")
895 {
...
940 }
941 if(type(ret) == "eps-svr"||type(ret) ==
"nu-svr"||type(ret)=="eps-bsvr"){
942 suppressWarnings(vgr<-split(sample(1:m,m),1:3))
943 pres <- NULL
944 for(i in 1:3)
945 {
946 cind <-
unsplit(vgr[-i],factor(rep((1:3)[-i],unlist(lapply(vgr[-i],length)))))
947
948 cret <-
ksvm(x[cind,],y[cind],type=type(ret),kernel=kernel,kpar =
NULL,C=C,nu=nu,epsilon=epsilon,tol=tol,scaled=FALSE, cross = 0, fit = FALSE,
cache = cache, prob.model = FALSE)
949 cres <- predict(cret, x[vgr[[i]],])
950 if (!is.null(scaling(ret)$y.scale))
951 cres <- cres * scaling(ret)$y.scale$"scaled:scale" +
scaling(ret)$y.scale$"scaled:center"
952 pres <- rbind(pres, cres)
953 }
954 pres[abs(pres) > (5*sd(pres))] <- 0
955 prob.model(ret) <- list(sum(abs(pres))/dim(pres)[1])
956 }
957 }
958
959 return(ret)
960 })
961
962
963
964
965 ## kernelmatrix interface
966
967 setMethod("ksvm",signature(x="kernelMatrix"),
968 function (x,
969 y = NULL,
970 type = NULL,
971 C = 1,
972 nu = 0.2,
973 epsilon = 0.1,
974 prob.model = FALSE,
975 class.weights = NULL,
976 cross = 0,
977 fit = TRUE,
978 cache = 40,
979 tol = 0.001,
980 shrinking = TRUE,
981 ...)
982 {
...
1723 })
1724
1725
1726
1727 .classAgreement <- function (tab) {
...
1737 }
1738
1739 ## List Interface
1740
1741
1742 setMethod("ksvm",signature(x="list"),
1743 function (x,
1744 y = NULL,
1745 type = NULL,
1746 kernel = "stringdot",
1747 kpar = list(length = 4, lambda = 0.5),
1748 C = 1,
1749 nu = 0.2,
1750 epsilon = 0.1,
1751 prob.model = FALSE,
1752 class.weights = NULL,
1753 cross = 0,
1754 fit = TRUE,
1755 cache = 40,
1756 tol = 0.001,
1757 shrinking = TRUE,
1758 ...
1759 ,na.action = na.omit)
1760 {
...
2630 })
2631
2632 ##**************************************************************#
2633 ## predict for matrix, data.frame input
2634
2635 setMethod("predict", signature(object = "ksvm"),
2636 function (object, newdata, type = "response", coupler = "minpair")
2637 {
2638 type <-
match.arg(type,c("response","probabilities","votes","decision"))
2639 if (missing(newdata) && type=="response" & !is.null(fitted(object)))
2640 return(fitted(object))
2641 else if(missing(newdata))
2642 stop("Missing data !")
2643
2644 if(!is(newdata,"list")){
2645 if (!is.null(terms(object)) & !is(newdata,"kernelMatrix"))
2646 {
2647 if(!is.matrix(newdata))
2648 newdata <- model.matrix(delete.response(terms(object)),
as.data.frame(newdata), na.action = n.action(object))
2649 }
2650 else
2651 newdata <- if (is.vector(newdata)) t(t(newdata)) else
as.matrix(newdata)
2652
2653
2654 newnrows <- nrow(newdata)
2655 newncols <- ncol(newdata)
2656 if(!is(newdata,"kernelMatrix") && !is.null(xmatrix(object))){
2657 if(is(xmatrix(object),"list") &&
is(xmatrix(object)[[1]],"matrix")) oldco <- ncol(xmatrix(object)[[1]])
2658 if(is(xmatrix(object),"matrix")) oldco <- ncol(xmatrix(object))
2659 if (oldco != newncols) stop ("test vector does not match model !")
2660 }
2661 }
2662 else
2663 newnrows <- length(newdata)
2664
2665 p <- 0
2666
2667 if (is.list(scaling(object)))
2668 newdata[,scaling(object)$scaled] <-
2669 scale(newdata[,scaling(object)$scaled, drop = FALSE],
2670 center = scaling(object)$x.scale$"scaled:center", scale =
scaling(object)$x.scale$"scaled:scale")
2671
2672 if(type == "response" || type =="decision" || type=="votes")
2673 {
2674
if(type(object)=="C-svc"||type(object)=="nu-svc"||type(object)=="C-bsvc")
2675 {
...
2706 }
2707
2708 if(type(object) == "spoc-svc")
2709 {
...
2721 }
2722
2723 if(type(object) == "kbb-svc")
2724 {
...
2756 }
2757 }
2758
2759 if(type == "probabilities")
2760 {
2761 if(is.null(prob.model(object)[[1]]))
2762 stop("ksvm object contains no probability model. Make sure you
set the paramater prob.model in ksvm during training.")
2763
2764
if(type(object)=="C-svc"||type(object)=="nu-svc"||type(object)=="C-bsvc")
2765 {
...
2780 }
2781 else
2782 stop("probability estimates only supported for C-svc, C-bsvc
and nu-svc")
2783 }
2784
2785 if(type(object) == "one-svc")
2786 {
...
2799 }
2800 else {
2801
if(type(object)=="eps-svr"||type(object)=="nu-svr"||type(object)=="eps-bsvr")
2802 {
2803 if(is(newdata,"kernelMatrix"))
2804 predres <- newdata %*% coef(object) - b(object)
2805 else
2806 predres <-
kernelMult(kernelf(object),newdata,xmatrix(object),coef(object)) - b(object)
2807 }
2808 else {
...
2829 }
2830 }
2831
2832 if (!is.null(scaling(object)$y.scale) & !is(newdata,"kernelMatrix") &
!is(newdata,"list"))
2833 ## return raw values, possibly scaled back
2834 return(predres * scaling(object)$y.scale$"scaled:scale" +
scaling(object)$y.scale$"scaled:center")
2835 else
2836 ##else: return raw values
2837 return(predres)
2838 })
2839
2840
2841
2842
#****************************************************************************************#
2843
2844 setMethod("show","ksvm",
2845 function(object){
...
2893 })
2894
2895
2896 setMethod("plot", signature(x = "ksvm", y = "missing"),
2897 function(x, data = NULL, grid = 50, slice = list(), ...) {
...
2985 })
2986
2987
2988 setGeneric(".probPlatt", function(deci, yres)
standardGeneric(".probPlatt"))
2989 setMethod(".probPlatt",signature(deci="ANY"),
2990 function(deci,yres)
2991 {
...
3098 })
3099
3100 ## Sigmoid predict function
3101
3102 .SigmoidPredict <- function(deci, A, B)
3103 {
...
3111 }
3112
Browse[2]> summary (as.vector (y)) ## << NOTE THIS
Min. 1st Qu. Median Mean 3rd Qu. Max.
-0.246000 -0.059480 0.006395 0.075260 0.086470 1.014000
Browse[2]> c(mean=mean(y), sd=sd(y)) ## << NOTE THIS
mean sd
0.07525931 0.26670593
Browse[2]>
debug: {
xtmp <- scale(x[, scaled])
x[, scaled] <- xtmp
x.scale <- attributes(xtmp)[c("scaled:center", "scaled:scale")]
if (is.numeric(y) && (type(ret) != "C-svc" && type(ret) !=
"nu-svc" && type(ret) != "C-bsvc" && type(ret) != "spoc-svc" &&
type(ret) != "kbb-svc")) {
y <- scale(y)
y.scale <- attributes(y)[c("scaled:center", "scaled:scale")]
y <- as.vector(y)
}
}
Browse[2]>
debug: xtmp <- scale(x[, scaled])
Browse[2]>
debug: x[, scaled] <- xtmp
Browse[2]>
debug: x.scale <- attributes(xtmp)[c("scaled:center", "scaled:scale")]
Browse[2]>
debug: if (is.numeric(y) && (type(ret) != "C-svc" && type(ret) != "nu-svc" &&
type(ret) != "C-bsvc" && type(ret) != "spoc-svc" && type(ret) !=
"kbb-svc")) {
y <- scale(y)
y.scale <- attributes(y)[c("scaled:center", "scaled:scale")]
y <- as.vector(y)
}
Browse[2]>
debug: {
y <- scale(y)
y.scale <- attributes(y)[c("scaled:center", "scaled:scale")]
y <- as.vector(y)
}
Browse[2]>
debug: y <- scale(y)
Browse[2]>
debug: y.scale <- attributes(y)[c("scaled:center", "scaled:scale")]
Browse[2]>
debug: y <- as.vector(y)
Browse[2]>
debug: ncols <- ncol(x)
Browse[2]> summary (y) ## << NOTE THIS
Min. 1st Qu. Median Mean 3rd Qu. Max.
-1.205e+00 -5.052e-01 -2.582e-01 -8.429e-18 4.202e-02 3.521e+00
Browse[2]> sd(y) ## << NOTE THIS
[1] 1
Browse[2]>
debug: if (type(ret) == "eps-svr") {
if (ktype == 4)
K <- kernelMatrix(kernel, x)
resv <- .Call("smo_optim", as.double(t(x)), as.integer(nrow(x)),
as.integer(ncol(x)), as.double(y), as.double(K), as.integer(if (sparse)
x...@ia else 0),
as.integer(if (sparse) x...@ja else 0), as.integer(sparse),
as.double(matrix(rep(-1, m))), as.integer(ktype), as.integer(3),
as.double(C), as.double(nu), as.double(epsilon), as.double(sigma),
as.integer(degree), as.double(offset), as.integer(0),
as.double(0), as.integer(0), as.double(cache), as.double(tol),
as.integer(shrinking), PACKAGE = "kernlab")
tmpres <- resv[c(-(m + 1), -(m + 2))]
alpha(ret) <- coef(ret) <- tmpres[tmpres != 0]
svindex <- alphaindex(ret) <- which(tmpres != 0)
xmatrix(ret) <- x[svindex, , drop = FALSE]
b(ret) <- resv[(m + 1)]
obj(ret) <- resv[(m + 2)]
param(ret)$epsilon <- epsilon
param(ret)$C <- C
}
debug: {
if (ktype == 4)
K <- kernelMatrix(kernel, x)
resv <- .Call("smo_optim", as.double(t(x)), as.integer(nrow(x)),
as.integer(ncol(x)), as.double(y), as.double(K), as.integer(if (sparse)
x...@ia else 0),
as.integer(if (sparse) x...@ja else 0), as.integer(sparse),
as.double(matrix(rep(-1, m))), as.integer(ktype), as.integer(3),
as.double(C), as.double(nu), as.double(epsilon), as.double(sigma),
as.integer(degree), as.double(offset), as.integer(0),
as.double(0), as.integer(0), as.double(cache), as.double(tol),
as.integer(shrinking), PACKAGE = "kernlab")
tmpres <- resv[c(-(m + 1), -(m + 2))]
alpha(ret) <- coef(ret) <- tmpres[tmpres != 0]
svindex <- alphaindex(ret) <- which(tmpres != 0)
xmatrix(ret) <- x[svindex, , drop = FALSE]
b(ret) <- resv[(m + 1)]
obj(ret) <- resv[(m + 2)]
param(ret)$epsilon <- epsilon
param(ret)$C <- C
}
Browse[2]>
debug: if (ktype == 4) K <- kernelMatrix(kernel, x)
Browse[2]>
debug: NULL
Browse[2]>
debug: resv <- .Call("smo_optim", as.double(t(x)), as.integer(nrow(x)),
as.integer(ncol(x)), as.double(y), as.double(K), as.integer(if (sparse)
x...@ia else 0),
as.integer(if (sparse) x...@ja else 0), as.integer(sparse),
as.double(matrix(rep(-1, m))), as.integer(ktype), as.integer(3),
as.double(C), as.double(nu), as.double(epsilon), as.double(sigma),
as.integer(degree), as.double(offset), as.integer(0), as.double(0),
as.integer(0), as.double(cache), as.double(tol), as.integer(shrinking),
PACKAGE = "kernlab")
Browse[2]>
debug: [1] 0
Browse[2]>
debug: [1] 0
Browse[2]>
debug: tmpres <- resv[c(-(m + 1), -(m + 2))]
Browse[2]>
debug: alpha(ret) <- coef(ret) <- tmpres[tmpres != 0]
Browse[2]>
debug: svindex <- alphaindex(ret) <- which(tmpres != 0)
Browse[2]>
debug: xmatrix(ret) <- x[svindex, , drop = FALSE]
Browse[2]>
debug: b(ret) <- resv[(m + 1)]
Browse[2]>
debug: obj(ret) <- resv[(m + 2)]
Browse[2]>
debug: param(ret)$epsilon <- epsilon
Browse[2]>
debug: param(ret)$C <- C
Browse[2]>
debug: kcall(ret) <- match.call()
Browse[2]>
debug: kernelf(ret) <- kernel
Browse[2]>
debug: ymatrix(ret) <- y
Browse[2]>
debug: SVindex(ret) <- sort(unique(svindex), method = "quick")
Browse[2]>
debug: nSV(ret) <- length(unique(svindex))
Browse[2]>
debug: if (nSV(ret) == 0) stop("No Support Vectors found. You may want to
change your parameters")
Browse[2]>
debug: NULL
Browse[2]>
debug: fitted(ret) <- if (fit) predict(ret, x) else NULL
Browse[2]> scaling(ret) ## << NOTE THIS
NULL
Browse[2]> fitted(ret) ## << NOTE THIS
NULL
Browse[2]>
debug: predict(ret, x)
Browse[2]>
debug: if (any(scaled)) scaling(ret) <- list(scaled = scaled, x.scale = x.scale,
y.scale = y.scale)
Browse[2]>
debug: scaling(ret) <- list(scaled = scaled, x.scale = x.scale, y.scale =
y.scale)
Browse[2]>
debug: if (fit) {
if (type(ret) == "C-svc" || type(ret) == "nu-svc" || type(ret) ==
"spoc-svc" || type(ret) == "kbb-svc" || type(ret) ==
"C-bsvc")
error(ret) <- 1 - .classAgreement(table(y, as.integer(fitted(ret))))
if (type(ret) == "one-svc")
error(ret) <- sum(!fitted(ret))/m
if (type(ret) == "eps-svr" || type(ret) == "nu-svr" || type(ret) ==
"eps-bsvr") {
if (!is.null(scaling(ret)$y.scale)) {
scal <- scaling(ret)$y.scale$"scaled:scale"
fitted(ret) <- fitted(ret) * scaling(ret)$y.scale$"scaled:scale" +
scaling(ret)$y.scale$"scaled:center"
}
else scal <- 1
error(ret) <- drop((scal^2) * crossprod(fitted(ret) -
y)/m)
}
}
Browse[2]> scaling(ret)$y.scale ## << NOTE THIS
$`scaled:center`
[1] 0.07525931
$`scaled:scale`
[1] 0.2667059
Browse[2]>
debug: {
if (type(ret) == "C-svc" || type(ret) == "nu-svc" || type(ret) ==
"spoc-svc" || type(ret) == "kbb-svc" || type(ret) ==
"C-bsvc")
error(ret) <- 1 - .classAgreement(table(y, as.integer(fitted(ret))))
if (type(ret) == "one-svc")
error(ret) <- sum(!fitted(ret))/m
if (type(ret) == "eps-svr" || type(ret) == "nu-svr" || type(ret) ==
"eps-bsvr") {
if (!is.null(scaling(ret)$y.scale)) {
scal <- scaling(ret)$y.scale$"scaled:scale"
fitted(ret) <- fitted(ret) * scaling(ret)$y.scale$"scaled:scale" +
scaling(ret)$y.scale$"scaled:center"
}
else scal <- 1
error(ret) <- drop((scal^2) * crossprod(fitted(ret) -
y)/m)
}
}
Browse[2]>
debug: if (type(ret) == "eps-svr" || type(ret) == "nu-svr" || type(ret) ==
"eps-bsvr") {
if (!is.null(scaling(ret)$y.scale)) {
scal <- scaling(ret)$y.scale$"scaled:scale"
fitted(ret) <- fitted(ret) * scaling(ret)$y.scale$"scaled:scale" +
scaling(ret)$y.scale$"scaled:center"
}
else scal <- 1
error(ret) <- drop((scal^2) * crossprod(fitted(ret) - y)/m)
}
Browse[2]>
debug: {
if (!is.null(scaling(ret)$y.scale)) {
scal <- scaling(ret)$y.scale$"scaled:scale"
fitted(ret) <- fitted(ret) * scaling(ret)$y.scale$"scaled:scale" +
scaling(ret)$y.scale$"scaled:center"
}
else scal <- 1
error(ret) <- drop((scal^2) * crossprod(fitted(ret) - y)/m)
}
Browse[2]>
debug: if (!is.null(scaling(ret)$y.scale)) {
scal <- scaling(ret)$y.scale$"scaled:scale"
fitted(ret) <- fitted(ret) * scaling(ret)$y.scale$"scaled:scale" +
scaling(ret)$y.scale$"scaled:center"
} else scal <- 1
Browse[2]>
debug: {
scal <- scaling(ret)$y.scale$"scaled:scale"
fitted(ret) <- fitted(ret) * scaling(ret)$y.scale$"scaled:scale" +
scaling(ret)$y.scale$"scaled:center"
}
Browse[2]>
debug: scal <- scaling(ret)$y.scale$"scaled:scale"
Browse[2]>
debug: fitted(ret) <- fitted(ret) * scaling(ret)$y.scale$"scaled:scale" +
scaling(ret)$y.scale$"scaled:center"
Browse[2]> drop((scal^2) * crossprod(fitted(ret) - y)/m) ## ** NOTE THIS, and
compare to ?? and ++ below
[1] 0.0008958718
Browse[2]> summary(as.vector(fitted(ret))) ## << NOTE THIS
Min. 1st Qu. Median Mean 3rd Qu. Max.
-1.091000 -0.510000 -0.263900 -0.001303 -0.022230 3.491000
Browse[2]> sd(as.vector(fitted(ret))) ## << NOTE THIS
[1] 0.9956
Browse[2]> c(Mean=mean(y), SD=sd(y)) ## << NOTE THIS
Mean SD
-8.428655e-18 1.000000e+00
Browse[2]>
debug: error(ret) <- drop((scal^2) * crossprod(fitted(ret) - y)/m)
Browse[2]> summary(as.vector(fitted(ret))) ## << NOTE THIS
Min. 1st Qu. Median Mean 3rd Qu. Max.
-0.215600 -0.060770 0.004874 0.074910 0.069330 1.006000
Browse[2]> sd(as.vector(fitted(ret))) ## << NOTE THIS
[1] 0.2655324
Browse[2]> summary(y) ## << NOTE THIS
Min. 1st Qu. Median Mean 3rd Qu. Max.
-1.205e+00 -5.052e-01 -2.582e-01 -8.429e-18 4.202e-02 3.521e+00
Browse[2]> sd(y) ## << NOTE THIS
[1] 1
Browse[2]>
debug: cross(ret) <- -1
Browse[2]> error (ret) ## ?? NOTE THIS
[1] 0.03891344
Browse[2]> drop (crossprod (fitted(ret) - (y * scal +
scaling(ret)$y.scale$"scaled:center"))/m) ## ++ NOTE THIS, and compare to **
and ?? above
[1] 0.0008958718
Browse[2]> drop (crossprod (fitted(ret) - (y * scal +
scaling(ret)$y.scale$"scaled:center"))/m) / (scal^2) ## << NOTE THIS
[1] 0.01259449
Browse[2]> drop (crossprod (fitted(ret) - (y * scal +
scaling(ret)$y.scale$"scaled:center"))/m) * (scal^2) ## << NOTE THIS
[1] 6.37252e-05
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
