SUPPORT VECTOR MACHINE. EXAMPLES WITH MATLAB
J. SMITH
CONTENTS SUPPORT VECTOR MACHINE 1.1 SUPPORT VECTOR MACHINE INTRODUCTION 1.1.1 Applications 1.1.2 Multiclass SVM 1.1.3 SVM regression 1.2 IMPLEMENTATION CLASSIFICATION WITH SUPPORT VECTOR MACHINE 2.1 SUPPORT VECTOR MACHINE 2.1.1 Separable Data 2.1.2 Nonseparable Data 2.1.3 Nonlinear Transformation with Kernels 2.2 USING SUPPORT VECTOR MACHINES 2.2.1 Training an SVM Classifier 2.2.2 Classifying New Data with an SVM Classifier 2.2.3 Tuning an SVM Classifier 2.2.4 Train SVM Classifiers Using a Gaussian Kernel 2.2.5 Train SVM Classifier Using Custom Kernel
2.2.6 Optimize a Cross-Validated SVM Classifier Using 2.2.7 Plot Posterior Probability Regions for SVM Classification Models 2.2.8 Analyze Images Using Linear Support Vector Machines 2.3 FUNCTIONS FOR SUPPORT VECTOR MACHINE CLASSIFICATION 2.3.1 fitcsvm 2.3.2 fitSVMPosterior 2.3.3 predict 2.3.4 templateSVM 2.3.5 fitclinear 2.3.6 templateLinear 2.3.7 fitcecoc 2.3.8 templateECOC 2.4 TRAIN SUPPORT VECTOR MACHINES USING CLASSIFICATION LEARNER APP 2.5 TRAIN CLASSIFICATION MODELS IN CLASSIFICATION LEARNER APP 2.5.1 What Is Supervised Machine Learning? 2.5.2 Automated Classifier Training 2.5.3 Manual Classifier Training 2.5.4 Parallel Classifier Training 2.5.5 Compare and Improve Classification Models 2.6 CHOOSE CLASSIFIER OPTIONS 2.6.1 Choose a Classifier Type 2.6.2 Decision Trees 2.6.3 Discriminant Analysis 2.6.4 Logistic Regression 2.6.5 Support Vector Machines 2.6.6 Nearest Neighbor Classifiers 2.6.7 Ensemble Classifiers 2.7 ASSESS CLASSIFIER PERFORMANCE IN CLASSIFICATION LEARNER
2.7.1 Check Performance in the History List 2.7.2 Plot Classifier Results 2.7.3 Check Performance Per Class in the Confusion Matrix 2.7.4 Check the ROC Curve FIT MULTICLASS MODELS FOR SUPPORT VECTOR MACHINES OR OTHER CLASSIFIERS. CLASSIFICATION LEARNER APP 3.1 MULTICLASS ECOC FUNCTIONS 3.1.1 fitcecoc 3.1.2 CompactClassificationECOC class 3.2 TRAIN ENSEMBLE CLASSIFIERS USING CLASSIFICATION LEARNER APP PREDICTIVE ANALYTICS. SUPPORT VECTOR MACHINE REGRESSION (SVR) AND GAUSSIAN PROCCES REGRESSION (GPR) 4.1 SUPPORT VECTOR MACHINE REGRESSION 4.1.1 fitrsvm 4.1.2 predict 4.1.3 fitrlinear 4.1.4 predict 4.2 GAUSSIAN PROCESS REGRESSION MODELS 4.2.1 fitrgp 4.2.2 predict 4.2.3 postFitStatistics
1Chapter 1
SUPPORT VECTOR MACHINE
1.1 SUPPORT VECTOR MACHINE INTRODUCTION In machine learning, support vector machines (SVMs, also support vector networks) are supervised learning models with associated learning algorithms that analyze data used for classification and regression analysis. Given a set of training examples, each marked as belonging to one or the other of two categories, an SVM training algorithm builds a model that assigns new examples to one category or the other, making it a nonprobabilistic binary linear classifier. An SVM model is a representation of the examples as points in space, mapped so that the examples of the separate categories are divided by a clear gap that is as wide as possible. New examples are then mapped into that same space and predicted to belong to a category based on which side of the gap they fall. In addition to performing linear classification, SVMs can efficiently perform a non-linear classification using what is called the kernel trick, implicitly mapping their inputs into high-dimensional feature spaces. When data are not labeled, supervised learning is not possible, and an unsupervised learning approach is required, which attempts to find natural clustering of the data to groups, and then map new data to these formed groups. The clustering algorithm which provides an improvement to the support vector machines is called support vector clustering and is often used in industrial applications either when data are not labeled or when only some data are labeled as a preprocessing for a classification pass. More formally, a support vector machine constructs a hyperplane or set of hyperplanes in a high- or infinite-dimensional space, which can be used for classification, regression, or other tasks. Intuitively, a good separation is achieved by the hyperplane that has the largest distance to the nearest trainingdata point of any class (so-called functional margin), since in general the larger the margin the lower the generalization error of the classifier. Whereas the original problem may be stated in a finite dimensional space, it often happens that the sets to discriminate are not linearly separable in that space. For this reason, it was proposed that the original finite-dimensional space be mapped into a much higher-dimensional space, presumably making the separation easier in that space. To keep the computational load reasonable,
the mappings used by SVM schemes are designed to ensure that dot products may be computed easily in terms of the variables in the original space, by defining them in terms of a kernel function selected to suit the problem. The hyperplanes in the higher-dimensional space are defined as the set of points whose dot product with a vector in that space is constant. The vectors defining the hyperplanes can be chosen to be linear combinations with parameters of images of feature vectors that occur in the data base. With this choice of a hyperplane, the points x in the feature space that are mapped into the hyperplane are defined by the relation: Note that if becomes small as grows further away from x, each term in the sum measures the degree of closeness of the test point x to the corresponding data base point . In this way, the sum of kernels above can be used to measure the relative nearness of each test point to the data points originating in one or the other of the sets to be discriminated. Note the fact that the set of points x mapped into any hyperplane can be quite convoluted as a result, allowing much more complex discrimination between sets which are not convex at all in the original space.
1.1.1 Applications SVMs can be used to solve various real world problems: ·
SVMs are helpful in text and hypertext categorization as their application can significantly reduce the need for labeled training instances in both the standard inductive and transductive settings.
·
Classification of images can also be performed using SVMs. Experimental results show that SVMs achieve significantly higher search accuracy than traditional query refinement schemes after just three to four rounds of relevance feedback. This is also true of image segmentation systems, including those using a modified version SVM that uses the privileged approach as suggested by Vapnik.
· ·
Hand-written characters can be recognized using SVM. The SVM algorithm has been widely applied in the biological and other sciences. They have been used to classify proteins with up to 90% of the compounds classified correctly. Permutation tests based on SVM weights have been suggested as a mechanism for interpretation of SVM models. Support vector machine weights have also been used to interpret SVM models in the past. Posthoc interpretation of support vector machine models in order to identify features used by the model to make predictions is a relatively new area of research with special significance in the biological sciences.
1.1.2 Multiclass SVM Multiclass SVM aims to assign labels to instances by using support vector machines, where the labels are drawn from a finite set of several elements. The dominant approach for doing so is to reduce the single multiclass problem into multiple binary classification problems. Common methods for such reduction include: ·
· ·
Building binary classifiers which distinguish (i) between one of the labels and the rest (one-versus-all) or (ii) between every pair of classes (one-versus-one). Classification of new instances for the one-versus-all case is done by a winner-takes-all strategy, in which the classifier with the highest output function assigns the class (it is important that the output functions be calibrated to produce comparable scores). For the one-versusone approach, classification is done by a max-wins voting strategy, in which every classifier assigns the instance to one of the two classes, then the vote for the assigned class is increased by one vote, and finally the class with the most votes determines the instance classification. Directed acyclic graph SVM (DAGSVM) Error-correcting output codesCrammer and Singer proposed a multiclass SVM method which casts the multiclass classification problem into a single optimization problem, rather than decomposing it into multiple binary classification problems.
1.1.3 SVM regression A version of SVM for regression was proposed in 1996 by Vladimir N. Vapnik, Harris Drucker, Christopher J. C. Burges, Linda Kaufman and Alexander J. Smola. This method is called support vector regression (SVR). The model produced by support vector classification (as described above) depends only on a subset of the training data, because the cost function for building the model does not care about training points that lie beyond the margin. Analogously, the model produced by SVR depends only on a subset of the training data, because the cost function for building the model ignores any training data close to the model prediction. Another SVM version known as least squares support vector machine (LS-SVM) has been proposed by Suykens and Vandewalle. Training the original SVR means solving:
where is a training sample with target value . The inner product plus intercept is the prediction for that sample, and is a free parameter that serves as a threshold: all predictions have to be within an range of the true predictions. Slack variables are usually added into the above to allow for errors and to allow approximation in the case the above problem is infeasible.
1.2 IMPLEMENTATION The parameters of the maximum-margin hyperplane are derived by solving the optimization. There exist several specialized algorithms for quickly solving the QP problem that arises from SVMs, mostly relying on heuristics for breaking the problem down into smaller, more-manageable chunks. Another approach is to use an interior point method that uses Newton-like iterations to find a solution of the Karush–Kuhn–Tucker conditions of the primal and dual problems. Instead of solving a sequence of broken down problems, this approach directly solves the problem altogether. To avoid solving a linear system involving the large kernel matrix, a low rank approximation to the matrix is often used in the kernel trick. Another common method is Platt's sequential minimal optimization (SMO) algorithm, which breaks the problem down into 2-dimensional sub-problems that are solved analytically, eliminating the need for a numerical optimization algorithm and matrix storage. This algorithm is conceptually simple, easy to implement, generally faster, and has better scaling properties for difficult SVM problems. The special case of linear support vector machines can be solved more efficiently by the same kind of algorithms used to optimize its close cousin, logistic regression; this class of algorithms includes sub-gradient descent (e.g., PEGASOS) and coordinate descent (e.g., LIBLINEAR). LIBLINEAR has some attractive training time properties. Each convergence iteration takes time linear in the time taken to read the train data and the iterations also have a Q-Linear Convergence property, making the algorithm extremely fast. The general kernel SVMs can also be solved more efficiently using subgradient descent (e.g. P-packSVM, especially when parallelization is allowed. Kernel SVMs are available in many machine learning toolkits, including LIBSVM, MATLAB, SAS, SVMlight, kernlab, scikitlearn, Shogun, Weka, Shark, JKernelMachines, OpenCV and others.
2Chapter 2
CLASSIFICATION WITH SUPPORT VECTOR MACHINE
2.1 SUPPORT VECTOR MACHINE 2.1.1 Separable Data You can use a support vector machine (SVM) when your data has exactly two classes. An SVM classifies data by finding the best hyperplane that separates all data points of one class from those of the other class. The best hyperplane for an SVM means the one with the largest margin between the two classes. Margin means the maximal width of the slab parallel to the hyperplane that has no interior data points. The support vectors are the data points that are closest to the separating hyperplane; these points are on the boundary of the slab. The following figure illustrates these definitions, with + indicating data points of type 1, and – indicating data points of type –1.
Mathematical Formulation: Primal. This discussion follows Hastie, Tibshirani, and Friedman and Christianini and Shawe-Taylor . The data for training is a set of points (vectors) xj along with their categories yj . For some dimension d, the xj ∊ Rd, and the yj = ±1. The equation of a hyperplane is f(x)=x′β+b=0
where β ∊ Rd and b is a real number. The following problem defines the best separating hyperplane (i.e., the decision boundary). Find β and b that minimize ||β|| such that for all data points (xj ,yj ),
The support vectors are the xj on the boundary, those for which . For mathematical convenience, the problem is usually given as the equivalent problem of minimizing . This is a quadratic programming problem. The optimal solution follows:
enables classification of a vector z as
is the classification score and represents the distance z is from the decision boundary. Mathematical Formulation: Dual. It is computationally simpler to solve the dual quadratic programming problem. To obtain the dual, take positive Lagrange multipliers αj multiplied by each constraint, and subtract from the objective function:
where you look for a stationary point of LP over β and b. Setting the gradient of LP to 0, you get
Substituting into LP, you get the dual LD:
which you maximize over αj ≥ 0. In general, many αj are 0 at the maximum. The nonzero αj in the solution to the dual problem define the hyperplane, as seen in this page, which gives β as the sum of αjyj xj . The data points xj corresponding to nonzero αj are the support vectors. The derivative of LD with respect to a nonzero αj is 0 at an optimum. This gives
In particular, this gives the value of b at the solution, by taking any j with nonzero αj. The dual is a standard quadratic programming problem. For example, the Optimization Toolbox™ quadprog solver solves this type of problem.
2.1.2 Nonseparable Data Your data might not allow for a separating hyperplane. In that case, SVM can use a soft margin, meaning a hyperplane that separates many, but not all data points. There are two standard formulations of soft margins. Both involve adding slack variables ξj and a penalty parameter C.
The L1-norm refers to using ξj as slack variables instead of their squares. The three solver options SMO, ISDA, and L1QPof fitcsvm minimize the L1norm problem.
subject to the same constraints. In these formulations, you can see that increasing C places more weight on the slack variables ξj, meaning the optimization attempts to make a stricter separation between classes. Equivalently, reducing C towards 0 makes misclassification less important. Mathematical Formulation: Dual. For easier calculations, consider the L1 dual problem to this soft-margin formulation. Using Lagrange multipliers μj, the function to minimize for the L1-norm problem is:
where you look for a stationary point of LP over β, b, and positive ξj. Setting the gradient of LP to 0, you get
The final set of inequalities, 0 ≤ αj ≤ C, shows why C is sometimes called a box constraint. C keeps the allowable values of the Lagrange multipliers αj in a "box", a bounded region. The gradient equation for b gives the solution b in terms of the set of nonzero αj, which correspond to the support vectors. You can write and solve the dual of the L2-norm problem in an analogous manner. For details, see Christianini and Shawe-Taylor, Chapter 6. fitcsvm Implementation. Both dual soft-margin problems are quadratic programming problems. Internally, fitcsvm has several different algorithms for solving the problems. ·
For one-class or binary classification, if you do not set a fraction of expected outliers in the data, then the default solver is Sequential Minimal Optimization (SMO). SMO minimizes the one-norm problem
by a series of two-point minimizations. During optimization, SMO respects the linear constraint ,and explicitly includes the bias term in the model. SMO is relatively fast. For more details on SMO, see. · For binary classification, if you set a fraction of expected outliers in the data, then the default solver is the Iterative Single Data Algorithm. Like SMO, ISDA solves the one-norm problem. Unlike SMO, ISDA minimizes by a series on one-point minimizations, does not respect the linear constraint, and does not explicitly include the bias term in the model. For more details on ISDA, see . · For one-class or binary classification, and if you have an Optimization Toolbox license, you can choose to use quadprog to solve the one-norm problem. quadprog uses a good deal of memory, but solves quadratic programs to a high degree of precision.
2.1.3 Nonlinear Transformation with Kernels Some binary classification problems do not have a simple hyperplane as a useful separating criterion. For those problems, there is a variant of the mathematical approach that retains nearly all the simplicity of an SVM separating hyperplane. This approach uses these results from the theory of reproducing kernels: ·
There is a class of functions G(x1,x2) with the following property. There is a linear space S and a function φ mapping x to S such that G(x1,x2) = <φ(x1),φ(x2)>.
The dot product takes place in the space S. · This class of functions includes: · Polynomials: For some positive integer p, G(x1,x2) = (1 + x1′x2)p. ·
Radial basis function (Gaussian): G(x1,x2) = exp(–∥x1–x2)∥2). · Multilayer perceptron or sigmoid (neural network): For a positive number p1 and a negative number p2, G(x1,x2) = tanh(p1x1′x2 + p2). The mathematical approach using kernels relies on the computational method of hyperplanes. All the calculations for hyperplane classification use nothing more than dot products. Therefore, nonlinear kernels can use identical calculations and solution algorithms, and obtain classifiers that are nonlinear. The resulting classifiers are hypersurfaces in some space S, but the space S does not have to be identified or examined.
2.2 USING SUPPORT VECTOR MACHINES As with any supervised learning model, you first train a support vector machine, and then cross validate the classifier. Use the trained machine to classify (predict) new data. In addition, to obtain satisfactory predictive accuracy, you can use various SVM kernel functions, and you must tune the parameters of the kernel functions.
2.2.1 Training an SVM Classifier Train, and optionally cross validate, an SVM classifier using fitcsvm. The most common syntax is: SVMModel = fitcsvm(X,Y,'KernelFunction','rbf','Standardize',true,'ClassNames', {'negClass','posClass'}); The inputs are: ·
X — Matrix of predictor data, where each row is one observation, and each column is one predictor. · Y — Array of class labels with each row corresponding to the value of the corresponding row in X. Y can be a character array, categorical, logical or numeric vector, or cell vector of character vectors. Column vector with each row corresponding to the value of the corresponding row in X. Y can be a categorical or character array, logical or numeric vector, or cell array of character vectors. ·
KernelFunction — The default value is 'linear' for two-class learning, which separates the data by a hyperplane. The value 'gaussian' (or 'rbf') is the default for one-class learning, and specifies to use the Gaussian (or radial basis function) kernel. An important step to successfully train an SVM classifier is to choose an appropriate kernel function.
·
— Flag indicating whether the software should standardize the predictors before training the classifier. ClassNames — Distinguishes between the negative and positive
·
Standardize
classes, or specifies which classes to include in the data. The negative class is the first element (or row of a character array), e.g., 'negClass', and the positive class is the second element (or row of a character array),
e.g., 'posClass'. ClassNames must be the same data type as Y. It is good practice to specify the class names, especially if you are comparing the performance of different classifiers.
The resulting, trained model (SVMModel) contains the optimized parameters from the SVM algorithm, enabling you to classify new data. · For more name-value pairs you can use to control the training, see the fitcsvm reference page. ·
2.2.2 Classifying New Data with an SVM Classifier Classify new data using predict. The syntax for classifying new data using a trained SVM classifier (SVMModel)
is:
[label,score] = predict(SVMModel,newX); The resulting vector, label, represents the classification of each row in X. score is an n-by-2 matrix of soft scores. Each row corresponds to a row in X, which is a new observation. The first column contains the scores for the observations being classified in the negative class, and the second column contains the scores observations being classified in the positive class. To estimate posterior probabilities rather than scores, first pass the trained SVM classifier (SVMModel) to fitPosterior, which fits a score-toposterior-probability transformation function to the scores. The syntax is: ScoreSVMModel = fitPosterior(SVMModel,X,Y); The property ScoreTransform of the classifier ScoreSVMModel contains the optimal transformation function. Pass ScoreSVMModel to predict. Rather than returning the scores, the output argument score contains the posterior probabilities of an observation being classified in the negative (column 1 of score) or positive (column 2 of score) class.
2.2.3 Tuning an SVM Classifier Try tuning parameters of your classifier according to this scheme: 1. Pass the data to fitcsvm, and set the name-value pair arguments 'KernelScale','auto'. Suppose that the trained SVM model is called SVMModel. The software uses a heuristic procedure to select the kernel scale. The heuristic procedure uses subsampling. Therefore, to reproduce results, set a random number seed using rng before training the classifier. 2. Cross validate the classifier by passing it to crossval. By default, the software conducts 10-fold cross validation. 3. Pass the cross-validated SVM model to kFoldLoss to estimate and retain the classification error. 4. Retrain the SVM classifier, but adjust the 'KernelScale' and 'BoxConstraint' name-value pair arguments. ·
BoxConstraint — One strategy is to try a geometric sequence of the box constraint parameter. For example, take 11 values, from 1e5 to 1e5 by a factor of 10. Increasing BoxConstraint might decrease the number of support vectors, but also might increase training time.
·
KernelScale — One strategy is to try a geometric sequence of the RBF sigma parameter scaled at the original kernel scale. Do this by:
a.
Retrieving the original kernel scale, e.g., ks, using dot
notation: ks = SVMModel.KernelParameters.Scale. b. Use as new kernel scales factors of the original. For example, multiply ks by the 11 values 1e-5 to 1e5, increasing by a factor of 10. Choose the model that yields the lowest classification error. You might want to further refine your parameters to obtain better accuracy. Start with your initial parameters and perform another cross-validation step, this time using a factor of 1.2. Alternatively, optimize your parameters with fminsearch.
2.2.4 Train SVM Classifiers Using a Gaussian Kernel This example shows how to generate a nonlinear classifier with Gaussian kernel function. First, generate one class of points inside the unit disk in two dimensions, and another class of points in the annulus from radius 1 to radius 2. Then, generates a classifier based on the data with the Gaussian radial basis function kernel. The default linear classifier is obviously unsuitable for this problem, since the model is circularly symmetric. Set the box constraint parameter to Inf to make a strict classification, meaning no misclassified training points. Other kernel functions might not work with this strict box constraint, since they might be unable to provide a strict classification. Even though the rbf classifier can separate the classes, the result can be overtrained. Generate 100 points uniformly distributed in the unit disk. To do so, generate a radius r as the square root of a uniform random variable, generate an angle t uniformly in (0, ), and put the point at (r cos( t ), r sin( t )). rng(1); % For reproducibility r = sqrt(rand(100,1)); % Radius t = 2*pi*rand(100,1); % Angle data1 = [r.*cos(t), r.*sin(t)]; % Points Generate 100 points uniformly distributed in the annulus. The radius is again proportional to a square root, this time a square root of the uniform distribution from 1 through 4. r2 = sqrt(3*rand(100,1)+1); % Radius t2 = 2*pi*rand(100,1); % Angle data2 = [r2.*cos(t2), r2.*sin(t2)]; % points Plot the points, and plot circles of radii 1 and 2 for comparison. figure; plot(data1(:,1),data1(:,2),'r.','MarkerSize',15) hold on plot(data2(:,1),data2(:,2),'b.','MarkerSize',15) ezpolar(@(x)1);ezpolar(@(x)2);
axis equal hold off
Put the data in one matrix, and make a vector of classifications.
data3 = [data1;data2]; theclass = ones(200,1); theclass(1:100) = -1;
to 'rbf' and BoxConstraint set to Inf. Plot the decision boundary and flag the support vectors. Train an SVM classifier with KernelFunction set
%Train the SVM Classifier cl = fitcsvm(data3,theclass,'KernelFunction','rbf',... 'BoxConstraint',Inf,'ClassNames',[-1,1]); % Predict scores over the grid d = 0.02; [x1Grid,x2Grid] = meshgrid(min(data3(:,1)):d:max(data3(:,1)),... min(data3(:,2)):d:max(data3(:,2))); xGrid = [x1Grid(:),x2Grid(:)]; [~,scores] = predict(cl,xGrid);
% Plot the data and the decision boundary figure; h(1:2) = gscatter(data3(:,1),data3(:,2),theclass,'rb','.'); hold on ezpolar(@(x)1); h(3) = plot(data3(cl.IsSupportVector,1),data3(cl.IsSupportVector,2),'ko'); contour(x1Grid,x2Grid,reshape(scores(:,2),size(x1Grid)),[0 0],'k'); legend(h,{'-1','+1','Support Vectors'}); axis equal hold off
fitcsvm generates a classifier that is close to a circle of radius 1. The difference is due to the random training data. Training with the default parameters makes a more nearly circular classification boundary, but one that misclassifies some training data. Also, the default value of BoxConstraint is 1, and, therefore, there are more support vectors. cl2 = fitcsvm(data3,theclass,'KernelFunction','rbf');
[~,scores2] = predict(cl2,xGrid); figure; h(1:2) = gscatter(data3(:,1),data3(:,2),theclass,'rb','.'); hold on ezpolar(@(x)1); h(3) = plot(data3(cl2.IsSupportVector,1),data3(cl2.IsSupportVector,2),'ko'); contour(x1Grid,x2Grid,reshape(scores2(:,2),size(x1Grid)),[0 0],'k'); legend(h,{'-1','+1','Support Vectors'}); axis equal hold off
2.2.5 Train SVM Classifier Using Custom Kernel This example shows how to use a custom kernel function, such as the sigmoid kernel, to train SVM classifiers, and adjust custom kernel function parameters. Generate a random set of points within the unit circle. Label points in the first and third quadrants as belonging to the positive class, and those in the second and fourth quadrants in the negative class. rng(1); % For reproducibility n = 100; % Number of points per quadrant r1 = sqrt(rand(2*n,1)); % Random radii t1 = [pi/2*rand(n,1); (pi/2*rand(n,1)+pi)]; % Random angles for Q1 and Q3 X1 = [r1.*cos(t1) r1.*sin(t1)]; % Polar-to-Cartesian conversion r2 = sqrt(rand(2*n,1)); t2 = [pi/2*rand(n,1)+pi/2; (pi/2*rand(n,1)-pi/2)]; % Random angles for Q2 and Q4 X2 = [r2.*cos(t2) r2.*sin(t2)]; X = [X1; X2]; % Predictors Y = ones(4*n,1); Y(2*n + 1:end) = -1; % Labels Plot the data. figure; gscatter(X(:,1),X(:,2),Y); title('Scatter Diagram of Simulated Data')
Write a function that accepts two matrices in the feature space as inputs, and transforms them into a Gram matrix using the sigmoid kernel. function G = mysigmoid(U,V) % Sigmoid kernel function with slope gamma and intercept c gamma = 1; c = -1; G = tanh(gamma*U*V' + c); end
Save this code as a file named mysigmoid on your MATLAB® path. Train an SVM classifier using the sigmoid kernel function. It is good practice to standardize the data. Mdl1 = fitcsvm(X,Y,'KernelFunction','mysigmoid','Standardize',true); Mdl1 is a ClassificationSVM classifier containing the estimated parameters.
Plot the data, and identify the support vectors and the decision boundary. % Compute the scores over a grid d = 0.02; % Step size of the grid [x1Grid,x2Grid] = meshgrid(min(X(:,1)):d:max(X(:,1)),... min(X(:,2)):d:max(X(:,2))); xGrid = [x1Grid(:),x2Grid(:)]; % The grid [~,scores1] = predict(Mdl1,xGrid); % The scores figure; h(1:2) = gscatter(X(:,1),X(:,2),Y); hold on h(3) = plot(X(Mdl1.IsSupportVector,1),... X(Mdl1.IsSupportVector,2),'ko','MarkerSize',10); % Support vectors contour(x1Grid,x2Grid,reshape(scores1(:,2),size(x1Grid)),[0 0],'k'); % Decision boundary title('Scatter Diagram with the Decision Boundary') legend({'-1','1','Support Vectors'},'Location','Best'); hold off
You can adjust the kernel parameters in an attempt to improve the shape of the decision boundary. This might also decrease the within-sample misclassification rate, but, you should first determine the out-of-sample misclassification rate. Determine the out-of-sample misclassification rate by using 10-fold cross validation. CVMdl1 = crossval(Mdl1); misclass1 = kfoldLoss(CVMdl1); misclass1 misclass1 = 0.1350 The out-of-sample misclassification rate is 13.5%. Write another sigmoid function, but Set gamma = 0.5;. function G = mysigmoid2(U,V) % Sigmoid kernel function with slope gamma and intercept c gamma = 0.5; c = -1; G = tanh(gamma*U*V' + c); end Save this code as a file named mysigmoid2 on your MATLAB® path. Train another SVM classifier using the adjusted sigmoid kernel. Plot the data and the decision region, and determine the out-of-sample misclassification rate. Mdl2 = fitcsvm(X,Y,'KernelFunction','mysigmoid2','Standardize',true); [~,scores2] = predict(Mdl2,xGrid); figure; h(1:2) = gscatter(X(:,1),X(:,2),Y);
hold on h(3) = plot(X(Mdl2.IsSupportVector,1),... X(Mdl2.IsSupportVector,2),'ko','MarkerSize',10); title('Scatter Diagram with the Decision Boundary') contour(x1Grid,x2Grid,reshape(scores2(:,2),size(x1Grid)),[0 0],'k'); legend({'-1','1','Support Vectors'},'Location','Best'); hold off CVMdl2 = crossval(Mdl2); misclass2 = kfoldLoss(CVMdl2); misclass2 misclass2 = 0.0450
After the sigmoid slope adjustment, the new decision boundary seems to provide a better within-sample fit, and the cross-validation rate contracts by more than 66%.
2.2.6 Optimize a Cross-Validated SVM Classifier Using This example shows how to optimize an SVM classification. The classification works on locations of points from a Gaussian mixture model. In The Elements of Statistical Learning, Hastie, Tibshirani, and Friedman (2009), page 17 describes the model. The model begins with generating 10 base points for a "green" class, distributed as 2-D independent normals with mean (1,0) and unit variance. It also generates 10 base points for a "red" class, distributed as 2-D independent normals with mean (0,1) and unit variance. For each class (green and red), generate 100 random points as follows: Choose a base point m of the appropriate color uniformly at random. 2. Generate an independent random point with 2-D normal distribution with mean m and variance I/5, where I is the 2-by-2 identity matrix. In this example, use a variance I/50 to show the advantage of optimization more clearly. After generating 100 green and 100 red points, classify them using fitcsvm. Then use bayesopt to optimize the parameters of the resulting SVM model with respect to cross validation. 1.
Generate the Points and Classifier Generate the 10 base points for each class. rng default grnpop = mvnrnd([1,0],eye(2),10); redpop = mvnrnd([0,1],eye(2),10);
View the base points. plot(grnpop(:,1),grnpop(:,2),'go') hold on plot(redpop(:,1),redpop(:,2),'ro') hold off
Since some red base points are close to green base points, it can be difficult to classify the data points based on location alone. Generate the 100 data points of each class. redpts = zeros(100,2);grnpts = redpts; for i = 1:100 grnpts(i,:) = mvnrnd(grnpop(randi(10),:),eye(2)*0.02); redpts(i,:) = mvnrnd(redpop(randi(10),:),eye(2)*0.02); end View the data points. figure plot(grnpts(:,1),grnpts(:,2),'go') hold on plot(redpts(:,1),redpts(:,2),'ro') hold off
Prepare Data For Classification Put the data into one matrix, and make a vector grp that labels the class of each point. cdata = [grnpts;redpts]; grp = ones(200,1); % Green label 1, red label -1 grp(101:200) = -1;
Prepare Cross-Validation Set up a partition for cross-validation. This step fixes the train and test sets that the optimization uses at each step. c = cvpartition(200,'KFold',10);
Prepare Variables for Bayesian Optimization Set up a function that takes an input z = [rbf_sigma,boxconstraint] and returns the cross-validation loss value of z. Take the components of z as positive, log-transformed variables between 1e-5 and 1e5. Choose a wide range, because you don't know which values are likely to be good.
sigma = optimizableVariable('sigma',[1e-5,1e5],'Transform','log'); box = optimizableVariable('box',[1e-5,1e5],'Transform','log');
Objective Function This function handle computes the cross-validation loss at parameters [sigma,box]. For details,
see ClassificationPartitionedModel.kfoldLoss. bayesopt passes the variable z to the objective function as a one-row table. minfn = @(z)kfoldLoss(fitcsvm(cdata,grp,'CVPartition',c,... 'KernelFunction','rbf','BoxConstraint',z.box,... 'KernelScale',z.sigma));
Optimize Classifier Search for the best parameters [sigma,box] using bayesopt. For reproducibility, choose the 'expected-improvement-plus' acquisition function. The default acquisition function depends on run time, and so can give varying results. results = bayesopt(minfn,[sigma,box],'IsObjectiveDeterministic',true,... 'AcquisitionFunctionName','expected-improvement-plus') |=================================================================================================| | Iter | Eval | Objective | Objective | BestSoFar | BestSoFar | sigma | box | | | result | | runtime | (observed) | (estim.) | | | |=================================================================================================| | 1 | Best | 0.61 | 0.30782 | 0.61 | 0.61 | 1.5413e-05 | 2.1437e-05 | | 2 | Best | 0.135 | 0.23405 | 0.135 | 0.135 | 0.10189 | 94776 | | 3 | Accept | 0.565 | 0.18192 | 0.135 | 0.135 | 0.0071441 | 0.13909 | | 4 | Accept | 0.345 | 0.21426 | 0.135 | 0.135 | 67574 | 0.0033331 | | 5 | Accept | 0.24 | 1.7557 | 0.135 | 0.135 | 20.917 | 94267 | | 6 | Accept | 0.61 | 0.2171 | 0.135 | 0.135 | 0.00051047 | 93559 | | 7 | Best | 0.085 | 0.58957 | 0.085 | 0.085 | 0.59941 | 86673 | | 8 | Accept | 0.1 | 0.24869 | 0.085 | 0.085 | 1.0415 | 41.262 | | 9 | Accept | 0.095 | 0.15796 | 0.085 | 0.085 | 0.29333 | 0.0050124 | | 10 | Accept | 0.31 | 0.15137 | 0.085 | 0.085 | 1071.6 | 22723 | | 11 | Accept | 0.105 | 0.27237 | 0.085 | 0.085 | 0.40285 | 98592 |
| 12 | Accept | 0.155 | 0.17668 | 0.085 | 0.085 | 0.66026 | 1.2988e-05 | | 13 | Accept | 0.115 | 0.8829 | 0.085 | 0.085 | 1.3592 | 98804 | | 14 | Best | 0.07 | 0.54174 | 0.07 | 0.07 | 0.55194 | 26624 | | 15 | Accept | 0.08 | 0.14238 | 0.07 | 0.07 | 0.32927 | 33.727 | | 16 | Accept | 0.1 | 0.16678 | 0.07 | 0.07 | 0.6138 | 1513.8 | | 17 | Accept | 0.24 | 0.13882 | 0.07 | 0.07 | 1.0342 | 0.05635 | | 18 | Accept | 0.15 | 0.13611 | 0.07 | 0.07 | 0.15297 | 8.1938e-05 | | 19 | Accept | 0.095 | 0.14209 | 0.07 | 0.07 | 0.20544 | 0.58666 | | 20 | Accept | 0.13 | 0.14041 | 0.07 | 0.07 | 0.11718 | 471.22 | |=================================================================================================| | Iter | Eval | Objective | Objective | BestSoFar | BestSoFar | sigma | box | | | result | | runtime | (observed) | (estim.) | | | |=================================================================================================| | 21 | Accept | 0.21 | 0.18427 | 0.07 | 0.07 | 4.8378 | 277.63 | | 22 | Accept | 0.085 | 0.11853 | 0.07 | 0.07 | 0.41308 | 4.9196 | | 23 | Accept | 0.345 | 0.13407 | 0.07 | 0.07 | 93941 | 89740 | | 24 | Accept | 0.145 | 0.14859 | 0.07 | 0.07 | 0.14794 | 0.011919 | | 25 | Accept | 0.345 | 0.14352 | 0.07 | 0.07 | 606.29 | 0.0088404 | | 26 | Accept | 0.115 | 0.25864 | 0.07 | 0.07 | 0.7624 | 19123 | | 27 | Accept | 0.11 | 0.20794 | 0.07 | 0.07 | 1.6927 | 365.26 | | 28 | Accept | 0.11 | 0.17592 | 0.07 | 0.07 | 0.38303 | 1.1387e-05 | | 29 | Accept | 0.3 | 0.16991 | 0.07 | 0.07 | 0.051596 | 3870.2 | | 30 | Accept | 0.115 | 0.18234 | 0.07 | 0.07 | 0.17599 | 15024 |
__________________________________________________________ Optimization completed. MaxObjectiveEvaluations of 30 reached. Total function evaluations: 30 Total elapsed time: 49.8954 seconds. Total objective function evaluation time: 8.5224 Best observed feasible point: sigma box _______ _____ 0.55194
26624
Observed objective function value = 0.07 Estimated objective function value = 0.07 Function evaluation time = 0.54174 Best estimated feasible point (according to models): sigma box _______ _____
0.55194
26624
Estimated objective function value = 0.07 Estimated function evaluation time = 0.35211
results = BayesianOptimization with properties: ObjectiveFcn: [function_handle] VariableDescriptions: [1×2 optimizableVariable] Options: [1×1 struct] MinObjective: 0.0700 XAtMinObjective: [1×2 table] MinEstimatedObjective: 0.0700 XAtMinEstimatedObjective: [1×2 table] NumObjectiveEvaluations: 30 TotalElapsedTime: 49.8954 NextPoint: [1×2 table] XTrace: [30×2 table] ObjectiveTrace: [30×1 double] ConstraintsTrace: [] UserDataTrace: {30×1 cell} ObjectiveEvaluationTimeTrace: [30×1 double] IterationTimeTrace: [30×1 double] ErrorTrace: [30×1 double] FeasibilityTrace: [30×1 logical] FeasibilityProbabilityTrace: [30×1 double] IndexOfMinimumTrace: [30×1 double] ObjectiveMinimumTrace: [30×1 double] EstimatedObjectiveMinimumTrace: [30×1 double]
Use the results to train a new, optimized SVM classifier. z(1) = results.XAtMinObjective.sigma; z(2) = results.XAtMinObjective.box; SVMModel = fitcsvm(cdata,grp,'KernelFunction','rbf',... 'KernelScale',z(1),'BoxConstraint',z(2));
Plot the classification boundaries. To visualize the support vector classifier, predict scores over a grid. d = 0.02; [x1Grid,x2Grid] = meshgrid(min(cdata(:,1)):d:max(cdata(:,1)),... min(cdata(:,2)):d:max(cdata(:,2))); xGrid = [x1Grid(:),x2Grid(:)]; [~,scores] = predict(SVMModel,xGrid); h = nan(3,1); % Preallocation figure; h(1:2) = gscatter(cdata(:,1),cdata(:,2),grp,'rg','+*'); hold on h(3) = plot(cdata(SVMModel.IsSupportVector,1),... cdata(SVMModel.IsSupportVector,2),'ko'); contour(x1Grid,x2Grid,reshape(scores(:,2),size(x1Grid)),[0 0],'k'); legend(h,{'-1','+1','Support Vectors'},'Location','Southeast'); axis equal hold off
Evaluate Accuracy on New Data
Generate and classify some new data points. grnobj = gmdistribution(grnpop,.2*eye(2)); redobj = gmdistribution(redpop,.2*eye(2)); newData = random(grnobj,10); newData = [newData;random(redobj,10)]; grpData = ones(20,1); grpData(11:20) = -1; % red = -1 v = predict(SVMModel,newData); g = nan(7,1); figure; h(1:2) = gscatter(cdata(:,1),cdata(:,2),grp,'rg','+*'); hold on h(3:4) = gscatter(newData(:,1),newData(:,2),v,'mc','**'); h(5) = plot(cdata(SVMModel.IsSupportVector,1),... cdata(SVMModel.IsSupportVector,2),'ko'); contour(x1Grid,x2Grid,reshape(scores(:,2),size(x1Grid)),[0 0],'k'); legend(h(1:5),{'-1 (training)','+1 (training)','-1 (classified)',... '+1 (classified)','Support Vectors'},'Location','Southeast'); axis equal hold off
See which new data points are correctly classified. Circle the correctly classified points in red, and the incorrectly classified points in black. mydiff = (v == grpData); % Classified correctly figure; h(1:2) = gscatter(cdata(:,1),cdata(:,2),grp,'rg','+*'); hold on h(3:4) = gscatter(newData(:,1),newData(:,2),v,'mc','**'); h(5) = plot(cdata(SVMModel.IsSupportVector,1),... cdata(SVMModel.IsSupportVector,2),'ko'); contour(x1Grid,x2Grid,reshape(scores(:,2),size(x1Grid)),[0 0],'k'); for ii = mydiff % Plot red squares around correct pts h(6) = plot(newData(ii,1),newData(ii,2),'rs','MarkerSize',12); end for ii = not(mydiff) % Plot black squares around incorrect pts h(7) = plot(newData(ii,1),newData(ii,2),'ks','MarkerSize',12); end legend(h,{'-1 (training)','+1 (training)','-1 (classified)',... '+1 (classified)','Support Vectors','Correctly Classified',...
'Misclassified'},'Location','Southeast'); hold off
2.2.7 Plot Posterior Probability Regions for SVM Classification Models This example shows how to predict posterior probabilities of SVM models over a grid of observations, and then plot the posterior probabilities over the grid. Plotting posterior probabilities exposes decision boundaries. Load Fisher's iris data set. Train the classifier using the petal lengths and widths, and remove the virginica species from the data. load fisheriris classKeep = ~strcmp(species,'virginica'); X = meas(classKeep,3:4); y = species(classKeep); Train an SVM classifier using the data. It is good practice to specify the order of the classes. SVMModel = fitcsvm(X,y,'ClassNames',{'setosa','versicolor'}); Estimate the optimal score transformation function. rng(1); % For reproducibility [SVMModel,ScoreParameters] = fitPosterior(SVMModel); ScoreParameters Warning: Classes are perfectly separated. The optimal score-to-posterior transformation is a step function. ScoreParameters = struct with fields: Type: 'step' LowerBound: -0.8431 UpperBound: 0.6897 PositiveClassProbability: 0.5000 The optimal score transformation function is the step function because the
classes are separable. The fields LowerBoundand UpperBound of ScoreParameters indicate the lower and upper end points of the interval of scores corresponding to observations within the class-separating hyperplanes (the margin). No training observation falls within the margin. If a new score is in the interval, then the software assigns the corresonding observation a positive class posterior probability, i.e., the value in the PositiveClassProbability field of ScoreParameters. Define a grid of values in the observed predictor space. Predict the posterior probabilities for each instance in the grid. xMax = max(X); xMin = min(X); d = 0.01; [x1Grid,x2Grid] = meshgrid(xMin(1):d:xMax(1),xMin(2):d:xMax(2)); [~,PosteriorRegion] = predict(SVMModel,[x1Grid(:),x2Grid(:)]); Plot the positive class posterior probability region and the training data. figure; contourf(x1Grid,x2Grid,... reshape(PosteriorRegion(:,2),size(x1Grid,1),size(x1Grid,2))); h = colorbar; h.Label.String = 'P({\it{versicolor}})'; h.YLabel.FontSize = 16; caxis([0 1]); colormap jet; hold on gscatter(X(:,1),X(:,2),y,'mc','.x',[15,10]); sv = X(SVMModel.IsSupportVector,:); plot(sv(:,1),sv(:,2),'yo','MarkerSize',15,'LineWidth',2); axis tight hold off
In two-class learning, if the classes are separable, then there are three regions: one where observations have positive class posterior probability 0, one where it is 1, and the other where it is the postiive class prior probability.
2.2.8 Analyze Images Using Linear Support Vector Machines This example shows how to determine which quadrant of an image a shape occupies by training an error-correcting output codes (ECOC) model comprised of linear SVM binary learners. This example also illustrates the disk-space consumption of ECOC models that store support vectors, their labels, and the estimated coefficients.
Create the Data Set Randomly place a circle with radius five in a 50-by-50 image. Make 5000 images. Create a label for each image indicating the quadrant that the circle occupies. Quadrant 1 is in the upper right, quadrant 2 is in the upper left, quadrant 3 is in the lower left, and quadrant 4 is in the lower right. The predictors are the intensities of each pixel. d = 50; % Height and width of the images in pixels n = 5e4; % Sample size X = zeros(n,d^2); % Predictor matrix preallocation Y = zeros(n,1); % Label preallocation theta = 0:(1/d):(2*pi); r = 5; % Circle radius rng(1); % For reproducibility for j = 1:n; figmat = zeros(d); % Empty image c = datasample((r + 1):(d - r - 1),2); % Random circle center x = r*cos(theta) + c(1); % Make the circle y = r*sin(theta) + c(2); idx = sub2ind([d d],round(y),round(x)); % Convert to linear indexing figmat(idx) = 1; % Draw the circle X(j,:) = figmat(:); % Store the data Y(j) = (c(2) >= floor(d/2)) + 2*(c(2) < floor(d/2)) + ... (c(1) < floor(d/2)) + ... 2*((c(1) >= floor(d/2)) & (c(2) < floor(d/2))); % Determine the
quadrant end
Plot an observation. figure; imagesc(figmat); h = gca; h.YDir = 'normal'; title(sprintf('Quadrant %d',Y(end)));
Train the ECOC Model Use a 25% holdout sample and specify the training and holdout sample indices. p = 0.25; CVP = cvpartition(Y,'Holdout',p); % Cross-validation data partition isIdx = training(CVP); % Training sample indices oosIdx = test(CVP); % Test sample indices Create an SVM template that specifies storing the support vectors of the binary learners. Pass it and the training data to fitcecoc to train the model.
Determine the training sample classification error. t = templateSVM('SaveSupportVectors',true); MdlSV = fitcecoc(X(isIdx,:),Y(isIdx),'Learners',t); isLoss = resubLoss(MdlSV) isLoss = 0 MdlSV is a trained ClassificationECOC multiclass model. It stores the training data and the support vectors of each binary learner. For large data sets, such as those in image analysis, the model can consume a lot of memory. Determine the amount of disk space that the ECOC model consumes. infoMdlSV = whos('MdlSV'); mbMdlSV = infoMdlSV.bytes/1.049e6 mbMdlSV = 763.5919 The model consumes
1477.5 MB.
Improve Model Efficiency You can assess out-of-sample performance. You can also assess whether the model has been overfit with a compacted model that does not contain the support vectors, their related parameters, and the training data. Discard the support vectors and related parameters from the trained ECOC model. Then, discard the training data from the resulting model by using compact. Mdl = discardSupportVectors(MdlSV); CMdl = compact(Mdl); info = whos('Mdl','CMdl');
[bytesCMdl,bytesMdl] = info.bytes; memReduction = 1 - [bytesMdl bytesCMdl]/infoMdlSV.bytes memReduction = 0.0626
0.9996
In this case, discarding the support vectors reduces the memory consumption by about 3%. Compacting and discarding support vectors reduces the size by about 99.99%. An alternative way to manage support vectors is to reduce their numbers during training by specifying a larger box constraint, such as 100. Though SVM models that use fewer support vectors are more desirable and consume less memory, increasing the value of the box constraint tends to increase the training time. Remove MdlSV and Mdl from the workspace. clear Mdl MdlSV;
Assess Holdout Sample Performance Calculate the classification error of the holdout sample. Plot a sample of the holdout sample predictions. oosLoss = loss(CMdl,X(oosIdx,:),Y(oosIdx)) yHat = predict(CMdl,X(oosIdx,:)); nVec = 1:size(X,1); oosIdx = nVec(oosIdx); figure; for j = 1:9; subplot(3,3,j) imagesc(reshape(X(oosIdx(j),:),[d d])); h = gca; h.YDir = 'normal'; title(sprintf('Quadrant: %d',yHat(j))) end
text(-1.33*d,4.5*d + 1,'Predictions','FontSize',17) oosLoss = 0
The model does not misclassify any holdout sample observations.
2.3 FUNCTIONS FOR SUPPORT VECTOR MACHINE CLASSIFICATION For greater accuracy and kernel-function choices on low- through mediumdimensional data sets, train a binary SVM model or a multiclass errorcorrecting output codes (ECOC) model containing SVM binary learners using the Classification Learner app. For greater flexibility, use the command-line interface to train a binary SVM model using fitcsvm or train a multiclass ECOC model composed of binary SVM learners using fitcecoc. For reduced computation time on high-dimensional data sets that fit in the MATLAB® Workspace, efficiently train a binary, linear classification model, such as a linear SVM model, using fitclinear or train a multiclass ECOC model composed of SVM models using fitcecoc.
Functions: Binary SVM fitcsvm fitSVM Posterior predict templateSVM
Train binary support vector machine classifier Fit posterior probabilities Predict labels using support vector machine classification model Support vector machine template
Binary Linear Classification fitclinear predict templateLinear
Fit linear classification model to high-dimensional data Predict labels for linear classification models Linear classification learner template
Multiclass ECOC fitcecoc predict templateECOC
Fit multiclass models for support vector machines or other classifiers Predict labels using multiclass, error-correcting output codes model Error-correcting output codes learner template
2.3.1 fitcsvm Train binary support vector machine classifier fitcsvm trains or cross-validates a support vector machine (SVM) model for two-class (binary) classification on a low- through moderatedimensional predictor data set. fitcsvmsupports mapping the predictor data using kernel functions, and supports SMO, ISDA, or L1 soft-margin minimization via quadratic programming for objective-function minimization. To train a linear SVM model for binary classification on a high-dimensional data set, that is, data sets that include many predictor variables, use fitclinear instead. For multiclass learning by combining binary SVM models, use errorcorrecting output codes (ECOC). For more details, see fitcecoc. To train an SVM regression model, see fitrsvm for low- through moderatedimensional predictor data sets, or fitrlinear for high-dimensional data sets.
Syntax Mdl = fitcsvm(Tbl,ResponseVarName) Mdl = fitcsvm(Tbl,formula) Mdl = fitcsvm(Tbl,Y) Mdl = fitcsvm(X,Y) Mdl = fitcsvm(___,Name,Value)
Description Mdl = fitcsvm(Tbl,ResponseVarName) returns a support vector machine classifier Mdl trained using the sample data contained in a table (Tbl). ResponseVarName is the name of the variable in Tbl that contains
the class labels for one- or two-class classification. Mdl = fitcsvm(Tbl,formula) returns an SVM classifer trained using the sample data contained in a table (Tbl). formula is an explanatory model of the response and a subset of predictor variables in Tbl used to fit Mdl. Mdl = fitcsvm(Tbl,Y) returns an SVM classifer trained using the predictor variables in table Tbl and class labels in vector Y. Mdl = fitcsvm(X,Y) returns an SVM classifier trained using the predictors in the matrix X and class labels in vector Y for one- or two-class classification. Mdl = fitcsvm(___,Name,Value) returns a support vector machine classifier with additional options specified by one or more Name,Value pair arguments, using any of the previous syntaxes. For example, you can specify the type of cross-validation, the cost for misclassification, or the type of score transformation function.
Examples: Train a Support Vector Machine Classifier Load Fisher's iris data set. Remove the sepal lengths and widths, and all observed setosa irises. load fisheriris inds = ~strcmp(species,'setosa'); X = meas(inds,3:4); y = species(inds); Train an SVM classifier using the processed data set. SVMModel = fitcsvm(X,y) SVMModel = ClassificationSVM
ResponseName: 'Y' CategoricalPredictors: [] ClassNames: {'versicolor' 'virginica'} ScoreTransform: 'none' NumObservations: 100 Alpha: [24×1 double] Bias: -14.4149 KernelParameters: [1×1 struct] BoxConstraints: [100×1 double] ConvergenceInfo: [1×1 struct] IsSupportVector: [100×1 logical] Solver: 'SMO' The Command Window shows that SVMModel is a trained ClassificationSVM classifier and a property list. Display the properties of SVMModel, for example, to determine the class order, by using dot notation. classOrder = SVMModel.ClassNames classOrder = 2×1 cell array 'versicolor' 'virginica' The first class ('versicolor') is the negative class, and the second ('virginica') is the positive class. You can change the class order during training by using the 'ClassNames' name-value pair argument. Plot a scatter diagram of the data and circle the support vectors. sv = SVMModel.SupportVectors; figure gscatter(X(:,1),X(:,2),y) hold on
plot(sv(:,1),sv(:,2),'ko','MarkerSize',10) legend('versicolor','virginica','Support Vector') hold off
The support vectors are observations that occur on or beyond their estimated class boundaries. You can adjust the boundaries (and therefore the number of support vectors) by setting a box constraint during training using the 'BoxConstraint' namevalue pair argument.
Train and Cross Validate an SVM Classifier Load the ionosphere data set. load ionosphere rng(1); % For reproducibility Train an SVM classifier using the radial basis kernel. Let the software find a scale value for the kernel function. It is good practice to standardize the predictors. SVMModel = fitcsvm(X,Y,'Standardize',true,'KernelFunction','RBF',...
'KernelScale','auto'); SVMModel is a trained ClassificationSVM classifier. Cross validate the SVM classifier. By default, the software uses 10-fold cross validation. CVSVMModel = crossval(SVMModel); CVSVMModel is a ClassificationPartitionedModel cross-validated classifier. Estimate the out-of-sample misclassification rate. classLoss = kfoldLoss(CVSVMModel) classLoss = 0.0484 The generalization rate is approximately 5%.
Detect Outliers Using SVM and One-Class Learning Load Fisher's iris data set. Remove the petal lengths and widths. Treat all irises as coming from the same class. load fisheriris X = meas(:,1:2); y = ones(size(X,1),1); Train an SVM classifier using the processed data set. Assume that 5% of the observations are outliers. It is good practice to standardize the predictors. rng(1); SVMModel = fitcsvm(X,y,'KernelScale','auto','Standardize',true,... 'OutlierFraction',0.05); SVMModel is a trained ClassificationSVM classifier. By default, the software uses the Gaussian kernel for one-class learning. Plot the observations and the decision boundary. Flag the support vectors
and potential outliers. svInd = SVMModel.IsSupportVector; h = 0.02; % Mesh grid step size [X1,X2] = meshgrid(min(X(:,1)):h:max(X(:,1)),... min(X(:,2)):h:max(X(:,2))); [~,score] = predict(SVMModel,[X1(:),X2(:)]); scoreGrid = reshape(score,size(X1,1),size(X2,2)); figure plot(X(:,1),X(:,2),'k.') hold on plot(X(svInd,1),X(svInd,2),'ro','MarkerSize',10) contour(X1,X2,scoreGrid) colorbar; title('{\bf Iris Outlier Detection via One-Class SVM}') xlabel('Sepal Length (cm)') ylabel('Sepal Width (cm)') legend('Observation','Support Vector') hold off
The boundary separating the outliers from the rest of the data occurs where
the contour value is 0. Verify that the fraction of observations with negative scores in the crossvalidated data is close to 5%. CVSVMModel = crossval(SVMModel); [~,scorePred] = kfoldPredict(CVSVMModel); outlierRate = mean(scorePred<0) outlierRate = 0.0467
Find Multiple Class Boundaries Using Binary SVM Load Fisher's iris data set. Use the petal lengths and widths. load fisheriris X = meas(:,3:4); Y = species; Examine a scatter plot of the data. figure gscatter(X(:,1),X(:,2),Y); h = gca; lims = [h.XLim h.YLim]; % Extract the x and y axis limits title('{\bf Scatter Diagram of Iris Measurements}'); xlabel('Petal Length (cm)'); ylabel('Petal Width (cm)'); legend('Location','Northwest');
There are three classes, one of which is linearly separable from the others. For each class: 1. Create a logical vector (indx) indicating whether an observation is a member of the class. 2. Train an SVM classifier using the predictor data and indx. 3. Store the classifier in a cell of a cell array. It is good practice to define the class order. SVMModels = cell(3,1); classes = unique(Y); rng(1); % For reproducibility for j = 1:numel(classes); indx = strcmp(Y,classes(j)); % Create binary classes for each classifier SVMModels{j} = fitcsvm(X,indx,'ClassNames',[false true],'Standardize',true,... 'KernelFunction','rbf','BoxConstraint',1); end SVMModels is a 3-by-1 cell array, with each cell containing
a ClassificationSVM classifier. For each cell, the positive class is setosa, versicolor, and virginica, respectively. Define a fine grid within the plot, and treat the coordinates as new observations from the distribution of the training data. Estimate the score of the new observations using each classifier. d = 0.02; [x1Grid,x2Grid] = meshgrid(min(X(:,1)):d:max(X(:,1)),... min(X(:,2)):d:max(X(:,2))); xGrid = [x1Grid(:),x2Grid(:)]; N = size(xGrid,1); Scores = zeros(N,numel(classes)); for j = 1:numel(classes); [~,score] = predict(SVMModels{j},xGrid); Scores(:,j) = score(:,2); % Second column contains positive-class scores end Each row of Scores contains three scores. The index of the element with the largest score is the index of the class to which the new class observation most likely belongs. Associate each new observation with the classifier that gives it the maximum score. [~,maxScore] = max(Scores,[],2); Color in the regions of the plot based on which class the corresponding new observation belongs. figure h(1:3) = gscatter(xGrid(:,1),xGrid(:,2),maxScore,... [0.1 0.5 0.5; 0.5 0.1 0.5; 0.5 0.5 0.1]); hold on h(4:6) = gscatter(X(:,1),X(:,2),Y); title('{\bf Iris Classification Regions}'); xlabel('Petal Length (cm)'); ylabel('Petal Width (cm)'); legend(h,{'setosa region','versicolor region','virginica region',...
'observed setosa','observed versicolor','observed virginica'},... 'Location','Northwest'); axis tight hold off
Optimize SVM Classifier This example shows how to optimize hyperparameters automatically using fitcsvm. The example uses the ionosphere data. Load the data. load ionosphere Find hyperparameters that minimize five-fold cross-validation loss by using automatic hyperparameter optimization. For reproducibility, set the random seed and use the 'expectedimprovement-plus' acquisition function. rng default Mdl = fitcsvm(X,Y,'OptimizeHyperparameters','auto',... 'HyperparameterOptimizationOptions',struct('AcquisitionFunctionName',...
'expected-improvement-plus')) |=================================================================================================| | Iter | Eval | Objective | Objective | BestSoFar | BestSoFar | BoxConstrain-| KernelScale | | | result | | runtime | (observed) | (estim.) | t | | |=================================================================================================| | 1 | Best | 0.12536 | 25.826 | 0.12536 | 0.12536 | 0.42371 | 0.006703 | | 2 | Accept | 0.13675 | 0.50478 | 0.12536 | 0.1261 | 0.71122 | 4.9578 | | 3 | Accept | 0.21083 | 25.42 | 0.12536 | 0.12557 | 462.63 | 0.018515 | | 4 | Accept | 0.35897 | 0.32799 | 0.12536 | 0.12539 | 0.0016958 | 247.22 | | 5 | Accept | 0.35897 | 0.27284 | 0.12536 | 0.12538 | 0.13395 | 926.61 | | 6 | Best | 0.12251 | 13.403 | 0.12251 | 0.12251 | 0.0010949 | 0.0021493 | | 7 | Accept | 0.24501 | 26.737 | 0.12251 | 0.12251 | 248.37 | 0.0014376 | | 8 | Accept | 0.12536 | 0.28826 | 0.12251 | 0.12252 | 0.0055388 | 0.10366 | | 9 | Accept | 0.13675 | 0.52324 | 0.12251 | 0.12251 | 1.1173 | 0.33611 | | 10 | Accept | 0.1339 | 0.76113 | 0.12251 | 0.12573 | 0.0034957 | 0.012172 | | 11 | Accept | 0.12251 | 23.311 | 0.12251 | 0.12242 | 0.011113 | 0.0010276 | | 12 | Accept | 0.13675 | 0.22352 | 0.12251 | 0.12252 | 0.0093098 | 0.84291 | | 13 | Accept | 0.1339 | 0.83004 | 0.12251 | 0.12256 | 977.4 | 6.4103 | | 14 | Accept | 0.12251 | 16.514 | 0.12251 | 0.12249 | 0.019781 | 0.0030908 | | 15 | Accept | 0.12821 | 20.518 | 0.12251 | 0.12112 | 0.0010033 | 0.0010166 | | 16 | Accept | 0.12821 | 0.35842 | 0.12251 | 0.12108 | 30.029 | 2.6682 | | 17 | Accept | 0.12536 | 2.185 | 0.12251 | 0.12113 | 0.089673 | 0.027124 | | 18 | Accept | 0.12536 | 20.008 | 0.12251 | 0.12276 | 0.023151 | 0.0027814 | | 19 | Accept | 0.12536 | 20.818 | 0.12251 | 0.12336 | 0.0092371 | 0.0018045 | | 20 | Best | 0.11966 | 18.16 | 0.11966 | 0.12191 | 0.078369 | 0.0066965 | |=================================================================================================| | Iter | Eval | Objective | Objective | BestSoFar | BestSoFar | BoxConstrain-| KernelScale | | | result | | runtime | (observed) | (estim.) | t | | |=================================================================================================| | 21 | Accept | 0.12251 | 13.601 | 0.11966 | 0.12201 | 0.078379 | 0.0073195 | | 22 | Accept | 0.12536 | 0.3259 | 0.11966 | 0.1221 | 61.648 | 10.882 | | 23 | Accept | 0.13105 | 0.28265 | 0.11966 | 0.12213 | 943.7 | 40.587 | | 24 | Accept | 0.12536 | 20.634 | 0.11966 | 0.12283 | 0.085431 | 0.0055761 | | 25 | Accept | 0.12251 | 13.001 | 0.11966 | 0.12275 | 0.087864 | 0.0087906 | | 26 | Accept | 0.35897 | 0.31639 | 0.11966 | 0.12237 | 984.11 | 993.79 | | 27 | Accept | 0.12821 | 0.40254 | 0.11966 | 0.12241 | 956.16 | 19.107 | | 28 | Accept | 0.12251 | 0.2791 | 0.11966 | 0.1224 | 0.028887 | 0.17111 | | 29 | Accept | 0.13105 | 0.38401 | 0.11966 | 0.12231 | 0.046715 | 0.097528 | | 30 | Accept | 0.1396 | 0.35153 | 0.11966 | 0.12239 | 0.0010417 | 0.19914 |
__________________________________________________________ Optimization completed. MaxObjectiveEvaluations of 30 reached. Total function evaluations: 30 Total elapsed time: 324.2057 seconds. Total objective function evaluation time: 266.5659 Best observed feasible point: BoxConstraint KernelScale
_____________ 0.078369
___________
0.0066965
Observed objective function value = 0.11966 Estimated objective function value = 0.12239 Function evaluation time = 18.1597 Best estimated feasible point (according to models): BoxConstraint KernelScale _____________ ___________ 0.078379
0.0073195
Estimated objective function value = 0.12239 Estimated function evaluation time = 14.8562
Mdl = ClassificationSVM ResponseName: 'Y' CategoricalPredictors: [] ClassNames: {'b' 'g'} ScoreTransform: 'none' NumObservations: 351 HyperparameterOptimizationResults: [1×1 BayesianOptimization] Alpha: [70×1 double] Bias: -20.9179 KernelParameters: [1×1 struct] BoxConstraints: [351×1 double] ConvergenceInfo: [1×1 struct] IsSupportVector: [351×1 logical] Solver: 'SMO'
2.3.2 fitSVMPosterior Fit posterior probabilities
Syntax ScoreSVMModel = fitSVMPosterior(SVMModel) ScoreSVMModel = fitSVMPosterior(SVMModel,TBL,ResponseVarName) ScoreSVMModel = fitSVMPosterior(SVMModel,TBL,Y) ScoreSVMModel = fitSVMPosterior(SVMModel,X,Y) ScoreSVMModel = fitSVMPosterior(___,Name,Value) [ScoreSVMModel,ScoreTransform] = fitSVMPosterior(___)
Description ScoreSVMModel = fitSVMPosterior(SVMModel) returns ScoreSVMModel, which is a trained, support vector machine (SVM) classifier containing the optimal score-toposterior-probability transformation function for two-class learning. The software fits the appropriate score-to-posterior-probability transformation function using the SVM classifierSVMModel, and by cross validation using the stored predictor data (SVMModel.X) and the class labels (SVMModel.Y). The transformation function computes the posterior probability that an observation is classified into the positive class (SVMModel.Classnames(2)). · If the classes are inseparable, then the transformation function is the sigmoid function. · If the classes are perfectly separable, the transformation function is the step function. · In two-class learning, if one of the two classes has a relative frequency of 0, then the transformation function is the constant function. fitSVMPosterior is not appropriate for one-class learning. · If SVMModel is a ClassificationSVM classifier, then the software
estimates the optimal transformation function by 10-fold cross validation as outlined in [1]. Otherwise, SVMModel must be a ClassificationPartitionedModel classifier. SVMModel specifies the cross-validation method. · The software stores the optimal transformation function in ScoreSVMModel.ScoreTransform. ScoreSVMModel = fitSVMPosterior(SVMModel,TBL,ResponseVarName) returns a trained support vector classifier containing the transformation function from the trained, compact SVM classifier SVMModel. The software estimates the score transformation function using predictor data in the table TBL and class labels TBL.ResponseVarName. ScoreSVMModel = fitSVMPosterior(SVMModel,TBL,Y) returns a trained support vector classifier containing the transformation function from the trained, compact SVM classifier SVMModel. The software estimates the score transformation function using predictor data in the table TBL and class labels Y. ScoreSVMModel = fitSVMPosterior(SVMModel,X,Y) returns a trained support vector classifier containing the transformation function from the trained, compact SVM classifier SVMModel. The software estimates the score transformation function using predictor data X and class labels Y. ScoreSVMModel = fitSVMPosterior(___,Name,Value) uses additional options specified by one or more Name,Value pair arguments provided SVMModel is a ClassificationSVM classifier. For example, you can specify the number of folds to use in k-fold cross validation. [ScoreSVMModel,ScoreTransform] = fitSVMPosterior(___) additionally returns the transformation function parameters (ScoreTransform) using any of the input arguments in the previous syntaxes.
Examples: Fit the Score-to-Posterior Probability Function for Separable Classes
Load Fisher's iris data set. Train the classifier using the petal lengths and widths, and remove the virginica species from the data. load fisheriris classKeep = ~strcmp(species,'virginica'); X = meas(classKeep,3:4); y = species(classKeep); gscatter(X(:,1),X(:,2),y); title('Scatter Diagram of Iris Measurements') xlabel('Petal length') ylabel('Petal width') legend('Setosa','Versicolor')
The classes are perfectly separable. Therefore, the score transformation function is a step function. Train an SVM classifier using the data. Cross validate the classifer using 10-fold cross validation (the default). rng(1); CVSVMModel = fitcsvm(X,y,'CrossVal','on');
CVSVMModel is a trained ClassificationPartitionedModel SVM classifier. Estimate the step function that transforms scores to posterior probabilities. [ScoreCVSVMModel,ScoreParameters] = fitSVMPosterior(CVSVMModel); Warning: Classes are perfectly separated. The optimal score-to-posterior transformation is a step function. fitSVMPosterior does the following: ·
Uses the data that the software stored in CVSVMModel to fit the transformation function · Warns whenever the classes are separable · Stores the step function in ScoreCSVMModel.ScoreTransform Display the score function type and its parameter values. ScoreParameters ScoreParameters = struct with fields: Type: 'step' LowerBound: -0.8431 UpperBound: 0.6897 PositiveClassProbability: 0.5000 ScoreParameters is a structure array with four fields: · ·
The score transformation function type (Type) The score corresponding to the negative class boundary (LowerBound) · The score corresponding to the positive class boundary (UpperBound) · The positive class probability (PositiveClassProbability) Since the classes are separable, the step function transforms the score to
either 0 or 1, which is the posterior probability that an observation is a versicolor iris.
Fit the Score-to-Posterior Probability Function for Inseparable Classes Load the ionosphere data set. load ionosphere The classes of this data set are not separable. Train an SVM classifier. Cross validate using 10-fold cross validation (the default). It is good practice to standardize the predictors and specify the class order. rng(1) % For reproducibility CVSVMModel = fitcsvm(X,Y,'ClassNames',{'b','g'},'Standardize',true,... 'CrossVal','on'); ScoreTransform = CVSVMModel.ScoreTransform ScoreTransform = 'none' CVSVMModel is a trained ClassificationPartitionedModel SVM classifier. The positive class is 'g'. The ScoreTransform property is none. Estimate the optimal score function for mapping observation scores to posterior probabilities of an observation being classified as 'g'. [ScoreCVSVMModel,ScoreParameters] = fitSVMPosterior(CVSVMModel); ScoreTransform = ScoreCVSVMModel.ScoreTransform ScoreParameters ScoreTransform = '@(S)sigmoid(S,-9.481576e-01,-1.218300e-01)' ScoreParameters =
struct with fields: Type: 'sigmoid' Slope: -0.9482 Intercept: -0.1218 ScoreTransform is the optimal score transform function. ScoreParameters contains the score transformation function, slope estimate, and the intercept estimate. You can estimate test-sample, posterior probabilities by passing ScoreCVSVMModel to kfoldPredict.
Estimate Posterior Probabilities for Test Samples Estimate positive class posterior probabilities for the test set of an SVM algorithm. Load the ionosphere data set. load ionosphere Train an SVM classifier. Specify a 20% holdout sample. It is good practice to standardize the predictors and specify the class order. rng(1) % For reproducibility CVSVMModel = fitcsvm(X,Y,'Holdout',0.2,'Standardize',true,... 'ClassNames',{'b','g'}); CVSVMModel is a trained ClassificationPartitionedModel cross-validated classifier. Estimate the optimal score function for mapping observation scores to posterior probabilities of an observation being classified as 'g'. ScoreCVSVMModel = fitSVMPosterior(CVSVMModel); ScoreSVMModel is a trained ClassificationPartitionedModel crossvalidated classifier containing the optimal score transformation function estimated from the training data. Estimate the out-of-sample positive class posterior probabilities. Display
the results for the first 10 out-of-sample observations. [~,OOSPostProbs] = kfoldPredict(ScoreCVSVMModel); indx = ~isnan(OOSPostProbs(:,2)); hoObs = find(indx); % Holdout observation numbers OOSPostProbs = [hoObs, OOSPostProbs(indx,2)]; table(OOSPostProbs(1:10,1),OOSPostProbs(1:10,2),... 'VariableNames',{'ObservationIndex','PosteriorProbability'}) ans = 10×2 table ObservationIndex PosteriorProbability ________________ ____________________ 6 7 8 9 16 22 23 24 38 41
0.17378 0.89637 0.0076609 0.91602 0.026718 4.6081e-06 0.9024 2.4129e-06 0.00042697 0.86427
2.3.3 predict Predict labels using support vector machine classification model
Syntax label = predict(SVMModel,X) [label,score] = predict(SVMModel,X)
Description label = predict(SVMModel,X) returns a vector of predicted class labels for the predictor data in the table or matrix X, based on the full or compact, trained SVM classification model SVMModel. [label,score] = predict(SVMModel,X) also returns a matrix of scores (score), indicating the likelihood that a label comes from a particular class. For SVM, likelihood measures are either classification scores or class posterior probabilities. For each observation in X, the predicted class label corresponds to the maximum score among all classes.
Examples: Label Test Sample Observations of SVM Classifiers Load the ionosphere data set. load ionosphere rng(1); % For reproducibility Train an SVM classifier. Specify a 15% holdout sample for testing. It is good practice to specify the class order and standardize the data. CVSVMModel = fitcsvm(X,Y,'Holdout',0.15,'ClassNames',{'b','g'},... 'Standardize',true); CompactSVMModel = CVSVMModel.Trained{1}; % Extract trained, compact classifier testInds = test(CVSVMModel.Partition); % Extract the test indices XTest = X(testInds,:); YTest = Y(testInds,:);
CVSVMModel is a ClassificationPartitionedModel classifier. It contains the property Trained, which is a 1-by-1 cell array holding a CompactClassificationSVMclassifier that the software trained using the training set. Label the test sample observations. Display the results for the first 10 observations in the test sample. [label,score] = predict(CompactSVMModel,XTest); table(YTest(1:10),label(1:10),score(1:10,2),'VariableNames',... {'TrueLabel','PredictedLabel','Score'}) ans = 10×3 table TrueLabel PredictedLabel Score _________ ______________ ________ 'b' 'g' 'b' 'g' 'b' 'g' 'b' 'b' 'g' 'g'
'b' 'g' 'b' 'g' 'b' 'g' 'b' 'b' 'g' 'g'
-1.7177 2.0003 -9.6841 2.5618 -1.548 2.0984 -2.7018 -0.66291 1.6046 1.7731
Predict Labels and Posterior Probabilities of SVM Classifiers A goal of classification is to predict labels of new observations using a trained algorithm. Many applications train algorithms on large data sets, which can use resources that are better used elsewhere. This example shows how to efficiently label new observations using an SVM classifier.
Load the ionosphere data set. Suppose that the last 10 observations become available after training the SVM classifier. load ionosphere n = size(X,1); % Training sample size isInds = 1:(n-10); % In-sample indices oosInds = (n-9):n; % Out-of-sample indices Train an SVM classifier. It is good practice to standardize the predictors and specify the order of the classes. Conserve memory by reducing the size of the trained SVM classifier. SVMModel = fitcsvm(X(isInds,:),Y(isInds),'Standardize',true,... 'ClassNames',{'b','g'}); CompactSVMModel = compact(SVMModel); whos('SVMModel','CompactSVMModel') Name
Size Attributes
Bytes Class
CompactSVMModel 1x1 29936 classreg.learning.classif.CompactClassificationSVM SVMModel 1x1 137268 ClassificationSVM The positive class is 'g'. The CompactClassificationSVM classifier (CompactSVMModel) uses less space than the ClassificationSVM classifier (SVMModel) because the latter stores the data. Estimate the optimal score-to-posterior-probability-transformation function. CompactSVMModel = fitPosterior(CompactSVMModel,... X(isInds,:),Y(isInds)) CompactSVMModel = classreg.learning.classif.CompactClassificationSVM ResponseName: 'Y'
CategoricalPredictors: [] ClassNames: {'b' 'g'} ScoreTransform: '@(S)sigmoid(S,-1.968452e+00,3.121267e-01)' Alpha: [88×1 double] Bias: -0.2143 KernelParameters: [1×1 struct] Mu: [1×34 double] Sigma: [1×34 double] SupportVectors: [88×34 double] SupportVectorLabels: [88×1 double] The optimal score transformation function (CompactSVMModel.ScoreTransform) is the sigmoid function because the classes are inseparable. Predict the out-of-sample labels and positive class posterior probabilities. Since true labels are available, compare them with the predicted labels. [labels,PostProbs] = predict(CompactSVMModel,X(oosInds,:)); table(Y(oosInds),labels,PostProbs(:,2),'VariableNames',... {'TrueLabels','PredictedLabels','PosClassPosterior'}) ans = 10×3 table TrueLabels PredictedLabels PosClassPosterior __________ _______________ _________________ 'g' 'g' 'g' 'g' 'g' 'g' 'g' 'g' 'g'
'g' 'g' 'g' 'g' 'g' 'g' 'g' 'g' 'g'
0.98419 0.95545 0.67792 0.94448 0.98744 0.92482 0.97111 0.96986 0.97803
'g'
'g'
0.94361
PostProbs is a 10-by-2 matrix, where the first column is the negative class posterior probabilities, and the second column is the positive class posterior probabilities corresponding to the new observations.
2.3.4 templateSVM Support vector machine template
Syntax t = templateSVM() t = templateSVM(Name,Value)
Description t = templateSVM() returns a support vector machine (SVM) learner template suitable for training error-correcting output code (ECOC) multiclass models. If you specify a default template, then the software uses default values for all input arguments during training. Specify t as a binary learner, or one in a set of binary learners, in fitcecoc to train an ECOC multiclass classifer. t = templateSVM(Name,Value) returns a template with additional options specified by one or more name-value pair arguments. For example, you can specify the box constraint, the kernel function, or whether to standardize the predictors. If you display t in the Command Window, then all options appear empty ([]), except those that you specify using name-value pair arguments. During training, the software uses default values for empty options.
Examples: Create a Default Support Vector Machine Template Use templateSVM to specify a default SVM template. t = templateSVM()
t= Fit template for classification SVM. Alpha: [0×1 double] BoxConstraint: [] CacheSize: [] CachingMethod: '' ClipAlphas: [] DeltaGradientTolerance: [] Epsilon: [] GapTolerance: [] KKTTolerance: [] IterationLimit: [] KernelFunction: '' KernelScale: [] KernelOffset: [] KernelPolynomialOrder: [] NumPrint: [] Nu: [] OutlierFraction: [] RemoveDuplicates: [] ShrinkagePeriod: [] Solver: '' StandardizeData: [] SaveSupportVectors: [] VerbosityLevel: [] Version: 2 Method: 'SVM' Type: 'classification' All properties of the template object are empty except for Method and Type. When you pass t to the training function, the software fills in the empty properties with their respective default values. For example, the software fills the KernelFunction property with 'linear'. For details on other default
values, see fitcsvm. t is a plan for an SVM learner, and no computation occurs when you specify it. You can pass t to fitcecoc to specify SVM binary learners for ECOC multiclass learning. However, by default, fitcecoc uses default SVM binary learners.
Create an SVM Template for ECOC Multiclass Learning Create a nondefault SVM template for use in fitcecoc. Load Fisher's iris data set. load fisheriris Create a template for SVM binary classifiers, and specify to use a Gaussian kernel function. t = templateSVM('KernelFunction','gaussian') t= Fit template for classification SVM. Alpha: [0×1 double] BoxConstraint: [] CacheSize: [] CachingMethod: '' ClipAlphas: [] DeltaGradientTolerance: [] Epsilon: [] GapTolerance: [] KKTTolerance: [] IterationLimit: [] KernelFunction: 'gaussian' KernelScale: [] KernelOffset: [] KernelPolynomialOrder: [] NumPrint: []
Nu: [] OutlierFraction: [] RemoveDuplicates: [] ShrinkagePeriod: [] Solver: '' StandardizeData: [] SaveSupportVectors: [] VerbosityLevel: [] Version: 2 Method: 'SVM' Type: 'classification' All properties of the template object are empty except for DistributionNames, Method, and Type. When trained on, the software fills in the empty properties with their respective default values. Specify t as a binary learner for an ECOC multiclass model. Mdl = fitcecoc(meas,species,'Learners',t); Mdl is a ClassificationECOC multiclass classifier. By default, the software trains Mdl using the one-versus-one coding design. Display the in-sample (resubstitution) misclassification error. L = resubLoss(Mdl,'LossFun','classiferror') L= 0.0200
Retain and Discard Support Vectors of SVM Binary Learners By default, fitcecoc empties the Alpha, SupportVectorLabels, and SupportVectors properties of the linear, SVM binary learners stored in the BinaryLearners property of the trained ECOC model. You can retain the support vectors and related values, and then discard them from the model. Load Fisher's iris data set.
load fisheriris rng(1); % For reproducibility Train an ECOC model using the entire data set. Specify retaining the support vectors by passing in the appropriate SVM template. t = templateSVM('SaveSupportVectors',true); MdlSV = fitcecoc(meas,species,'Learners',t); Mdl is a trained ClassificationECOC model. By default, fitcecoc uses linear, SVM binary learners. It implements a one-versus-one coding design, which requires three binary learners for three-class learning. Access the estimated values using dot notation. alpha = cell(3,1); alpha{1} = MdlSV.BinaryLearners{1}.Alpha; alpha{2} = MdlSV.BinaryLearners{2}.Alpha; alpha{3} = MdlSV.BinaryLearners{3}.Alpha; alpha alpha = 3×1 cell array [ 3×1 double] [ 3×1 double] [23×1 double] alpha is a 3-by-1 cell array that stores the estimated values of . Discard the support vectors and related values from the ECOC model. Mdl = discardSupportVectors(MdlSV); Mdl is similar to MdlSV, except that the Alpha, SupportVectorLabels, and SupportVectors of all linear SVM binary learners are empty ([]). areAllEmpty = @(x)isempty([x.Alpha x.SupportVectors x.SupportVectorLabels]); cellfun(areAllEmpty,Mdl.BinaryLearners) ans =
3×1 logical array 1 1 1 Compare the sizes of the two ECOC models. vars = whos('MdlSV','Mdl'); 100*(1 - vars(1).bytes/vars(2).bytes) ans = 5.1575 Mdl is about 5% smaller than MdlSV. Reduce your memory footprint by compacting Mdl, and then clearing MdlSV and Mdl from the workspace. CMdl = compact(Mdl); clear MdlSV Mdl; Predict the label for a random row of the training data using the more efficient SVM model. idx = randsample(size(meas,1),1) predictedLabel = predict(CMdl,meas(idx,:)) trueLabel = species(idx) idx = 63 predictedLabel = cell 'versicolor'
trueLabel = cell 'versicolor'
2.3.5 fitclinear Fit linear classification model to high-dimensional data fitclinear trains linear classification models for two-class (binary) learning with high-dimensional, full or sparse predictor data. Available linear classification models include regularized support vector machines (SVM) and logistic regression models. fitclinear minimizes the objective function using techniques that reduce computing time (e.g., stochastic gradient descent). A high-dimensional data set includes many predictor variables. Although such a data set can consume a significant fraction of memory, it must fit in the MATLAB® Workspace. To train a linear classification model for multiclass learning by combining SVM or logistic regression binary classifiers using error-correcting output codes, see fitcecoc.
Syntax Mdl = fitclinear(X,Y) Mdl = fitclinear(X,Y,Name,Value) [Mdl,FitInfo] = fitclinear(___) [Mdl,FitInfo,HyperparameterOptimizationResults] = fitclinear(___)
Description Mdl = fitclinear(X,Y) returns a trained linear classification model object that contains the results of fitting a binary support vector machine to the predictors Xand class labels Y. Mdl = fitclinear(X,Y,Name,Value) returns a trained linear classification model with additional options specified by one or more Name,Value pair arguments. For example, you can specify that the columns of the predictor matrix correspond to observations, implement logistic regression, or specify to cross-validate. It is good practice to cross-validate using
the Kfold Name,Value pair argument. The cross-validation results determine how well the model generalizes. [Mdl,FitInfo] = fitclinear(___) also returns optimization details using any of the previous syntaxes. You cannot request FitInfo for cross-validated models. [Mdl,FitInfo,HyperparameterOptimizationResults] = fitclinear(___) also returns hyperparameter optimization details when you pass an OptimizeHyperparameters name-value pair.
Examples Train Linear Classification Model Train a binary, linear classification model using support vector machines, dual SGD, and ridge regularization. Load the NLP data set. load nlpdata X is a sparse matrix of predictor data, and Y is a categorical vector of class labels. There are more than two classes in the data. Identify the labels that correspond to the Statistics and Machine Learning Toolbox™ documentation web pages. Ystats = Y == 'stats'; Train a binary, linear classification model that can identify whether the word counts in a documentation web page are from the Statistics and Machine Learning Toolbox™ documentation. Train the algorithm using the entire data set. Determine how well the optimization algorithm fit the model to the data by extracting a fit summary. rng(1); % For reproducibility [Mdl,FitInfo] = fitclinear(X,Ystats) Mdl =
ClassificationLinear ResponseName: 'Y' ClassNames: [0 1] ScoreTransform: 'none' Beta: [34023×1 double] Bias: -1.0059 Lambda: 3.1674e-05 Learner: 'svm' FitInfo = struct with fields: Lambda: 3.1674e-05 Objective: 5.3783e-04 PassLimit: 10 NumPasses: 10 BatchLimit: [] NumIterations: 238561 GradientNorm: NaN GradientTolerance: 0 RelativeChangeInBeta: 0.0562 BetaTolerance: 1.0000e-04 DeltaGradient: 1.4582 DeltaGradientTolerance: 1 TerminationCode: 0 TerminationStatus: {'Maximal number of iterations has been reached.'} Alpha: [31572×1 double] History: [] FitTime: 0.3981 Solver: {'dual'} Mdl is a ClassificationLinear model. You can pass Mdl and the training or
new data to loss to inspect the in-sample classification error. Or, you can pass Mdl and new predictor data to predict to predict class labels for new observations. FitInfo is a structure array containing, among other things, the termination status (TerminationStatus) and how long the solver took to fit the model to the data (FitTime). It is good practice to use FitInfo to determine whether optimization-termination measurements are satisfactory. Because training time is small, you can try to retrain the model, but increase the number of passes through the data. This can improve measures like DeltaGradient.
Find Good Lasso Penalty Using Cross-Validation To determine a good lasso-penalty strength for a linear classification model that uses a logistic regression learner, implement 5-fold cross-validation. Load the NLP data set. load nlpdata X is a sparse matrix of predictor data, and Y is a categorical vector of class labels. There are more than two classes in the data. The models should identify whether the word counts in a web page are from the Statistics and Machine Learning Toolbox™ documentation. So, identify the labels that correspond to the Statistics and Machine Learning Toolbox™ documentation web pages. Ystats = Y == 'stats'; Create a set of 11 logarithmically-spaced regularization strengths from through . Lambda = logspace(-6,-0.5,11); Cross-validate the models. To increase execution speed, transpose the predictor data and specify that the observations are in columns. Estimate the coefficients using SpaRSA. Lower the tolerance on the gradient of the objective function to 1e-8.
X = X'; rng(10); % For reproducibility CVMdl = fitclinear(X,Ystats,'ObservationsIn','columns','KFold',5,... 'Learner','logistic','Solver','sparsa','Regularization','lasso',... 'Lambda',Lambda,'GradientTolerance',1e-8) numCLModels = numel(CVMdl.Trained) CVMdl = classreg.learning.partition.ClassificationPartitionedLinear CrossValidatedModel: 'Linear' ResponseName: 'Y' NumObservations: 31572 KFold: 5 Partition: [1×1 cvpartition] ClassNames: [0 1] ScoreTransform: 'none'
numCLModels = 5 CVMdl is a ClassificationPartitionedLinear model. Because fitclinear implements 5-fold cross-validation, CVMdl contains 5 ClassificationLinear models that the software trains on each fold. Display the first trained linear classification model. Mdl1 = CVMdl.Trained{1} Mdl1 = ClassificationLinear ResponseName: 'Y'
ClassNames: [0 1] ScoreTransform: 'logit' Beta: [34023×11 double] Bias: [1×11 double] Lambda: [1×11 double] Learner: 'logistic' Mdl1 is a ClassificationLinear model object. fitclinear constructed Mdl1 by training on the first four folds. Because Lambda is a sequence of regularization strengths, you can think of Mdl1 as 11 models, one for each regularization strength in Lambda. Estimate the cross-validated classification error. ce = kfoldLoss(CVMdl); Because there are 11 regularization strengths, ce is a 1-by-11 vector of classification error rates. Higher values of Lambda lead to predictor variable sparsity, which is a good quality of a classifier. For each regularization strength, train a linear classification model using the entire data set and the same options as when you cross-validated the models. Determine the number of nonzero coefficients per model. Mdl = fitclinear(X,Ystats,'ObservationsIn','columns',... 'Learner','logistic','Solver','sparsa','Regularization','lasso',... 'Lambda',Lambda,'GradientTolerance',1e-8); numNZCoeff = sum(Mdl.Beta~=0); In the same figure, plot the cross-validated, classification error rates and frequency of nonzero coefficients for each regularization strength. Plot all variables on the log scale. figure; [h,hL1,hL2] = plotyy(log10(Lambda),log10(ce),... log10(Lambda),log10(numNZCoeff)); hL1.Marker = 'o'; hL2.Marker = 'o';
ylabel(h(1),'log_{10} classification error') ylabel(h(2),'log_{10} nonzero-coefficient frequency') xlabel('log_{10} Lambda') title('Test-Sample Statistics') hold off
Choose the index of the regularization strength that balances predictor variable sparsity and low classification error. In this case, a value between to should suffice. idxFinal = 7; Select the model from Mdl with the chosen regularization strength. MdlFinal = selectModels(Mdl,idxFinal); MdlFinal is a ClassificationLinear model containing one regularization strength. To estimate labels for new observations, pass MdlFinal and the new data to predict.
Optimize Linear Classifier
This example shows how to minimize the cross-validation error in a linear classifier using fitclinear. The example uses the NLP data set. Load the NLP data set. load nlpdata X is a sparse matrix of predictor data, and Y is a categorical vector of class labels. There are more than two classes in the data. The models should identify whether the word counts in a web page are from the Statistics and Machine Learning Toolbox™ documentation. Identify the relevant labels. X = X'; Ystats = Y == 'stats'; Optimize the classification using the 'auto' parameters. For reproducibility, set the random seed and use the 'expectedimprovement-plus' acquisition function. rng default Mdl = fitclinear(X,Ystats,'ObservationsIn','columns','Solver','sparsa',... 'OptimizeHyperparameters','auto','HyperparameterOptimizationOptions',... struct('AcquisitionFunctionName','expected-improvement-plus')) |=================================================================================================| | Iter | Eval | Objective | Objective | BestSoFar | BestSoFar | Lambda | Learner | | | result | | runtime | (observed) | (estim.) | | | |=================================================================================================| | 1 | Best | 0.00082351 | 48.413 | 0.00082351 | 0.00082351 | 7.5709e-06 | svm | | 2 | Accept | 0.0010452 | 103.77 | 0.00082351 | 0.00093425 | 1.7949e-05 | logistic | | 3 | Accept | 0.049221 | 28.411 | 0.00082351 | 0.00087814 | 0.87651 | svm | | 4 | Best | 0.00069682 | 30.61 | 0.00069682 | 0.00076227 | 7.6379e-10 | logistic | | 5 | Accept | 0.00069682 | 38.516 | 0.00069682 | 0.00069538 | 9.4937e-09 | logistic | | 6 | Accept | 0.00082351 | 36.471 | 0.00069682 | 0.00069536 | 1.0303e-09 | svm | | 7 | Accept | 0.049221 | 2.9968 | 0.00069682 | 0.00069303 | 3.162 | logistic | | 8 | Accept | 0.00076017 | 42.406 | 0.00069682 | 0.00069167 | 1.1517e-07 | svm | | 9 | Accept | 0.00091854 | 128.34 | 0.00069682 | 0.00069748 | 8.5525e-07 | logistic | | 10 | Accept | 0.00076017 | 27.156 | 0.00069682 | 0.00069692 | 1.2319e-06 | svm | | 11 | Accept | 0.00079184 | 31.88 | 0.00069682 | 0.00070577 | 7.7482e-09 | svm | | 12 | Accept | 0.00085519 | 31.09 | 0.00069682 | 0.0007018 | 3.187e-10 | svm | | 13 | Best | 0.00066515 | 22.217 | 0.00066515 | 0.00066668 | 3.1741e-10 | logistic | | 14 | Accept | 0.00091854 | 93.362 | 0.00066515 | 0.00067434 | 4.3736e-06 | logistic | | 15 | Accept | 0.00069682 | 25.07 | 0.00066515 | 0.00068162 | 2.1114e-09 | logistic | | 16 | Accept | 0.00066515 | 17.759 | 0.00066515 | 0.00068773 | 4.998e-08 | logistic | | 17 | Accept | 0.00082351 | 30.693 | 0.00066515 | 0.00068794 | 2.3401e-06 | svm |
| 18 | Accept | 0.00069682 | 15.348 | 0.00066515 | 0.00068453 | 3.1745e-10 | logistic | | 19 | Accept | 0.00069682 | 14.783 | 0.00066515 | 0.00068214 | 2.0723e-08 | logistic | | 20 | Accept | 0.00069682 | 13.158 | 0.00066515 | 0.00068825 | 3.1958e-10 | logistic | |=================================================================================================| | Iter | Eval | Objective | Objective | BestSoFar | BestSoFar | Lambda | Learner | | | result | | runtime | (observed) | (estim.) | | | |=================================================================================================| | 21 | Accept | 0.00069682 | 13.542 | 0.00066515 | 0.0006881 | 6.3913e-09 | logistic | | 22 | Accept | 0.00069682 | 16.579 | 0.00066515 | 0.00068209 | 4.7906e-08 | logistic | | 23 | Accept | 0.00069682 | 17.758 | 0.00066515 | 0.000682 | 9.8311e-10 | logistic | | 24 | Accept | 0.00069682 | 18.664 | 0.00066515 | 0.00068545 | 1.8854e-08 | logistic | | 25 | Accept | 0.00085519 | 16.133 | 0.00066515 | 0.00068545 | 3.1392e-08 | svm | | 26 | Accept | 0.00069682 | 10.988 | 0.00066515 | 0.00068548 | 3.2059e-10 | logistic | | 27 | Accept | 0.00069682 | 12.791 | 0.00066515 | 0.00068545 | 9.4312e-10 | logistic | | 28 | Accept | 0.022172 | 25.117 | 0.00066515 | 0.00068504 | 0.014274 | logistic | | 29 | Accept | 0.00069682 | 8.4165 | 0.00066515 | 0.00068401 | 7.0427e-08 | logistic | | 30 | Accept | 0.0017104 | 15.095 | 0.00066515 | 0.00068978 | 0.0048555 | svm |
__________________________________________________________ Optimization completed. MaxObjectiveEvaluations of 30 reached. Total function evaluations: 30 Total elapsed time: 986.7883 seconds. Total objective function evaluation time: 937.5269 Best observed feasible point: Lambda Learner __________ ________ 3.1741e-10
logistic
Observed objective function value = 0.00066515 Estimated objective function value = 0.00068978 Function evaluation time = 22.2169 Best estimated feasible point (according to models): Lambda Learner __________ ________ 3.2059e-10
logistic
Estimated objective function value = 0.00068978
Estimated function evaluation time = 16.0241
Mdl = ClassificationLinear ResponseName: 'Y' ClassNames: [0 1] ScoreTransform: 'logit' Beta: [34023×1 double] Bias: -9.7737 Lambda: 3.2059e-10 Learner: 'logistic'
2.3.6 templateLinear Linear classification learner template templateLinear creates a template suitable for fitting a linear classification model to high-dimensional data for multiclass problems. The template specifies the binary learner model, regularization type and strength, and solver, among other things. After creating the template, train the model by passing the template and data to fitcecoc.
Syntax t = templateLinear() t = templateLinear(Name,Value)
Description t = templateLinear() returns a linear classification learner template. If you specify a default template, then the software uses default values for all input arguments during training. t = templateLinear(Name,Value) returns a template with additional options specified by one or more name-value pair arguments. For example, you can specify to implement logistic regression, specify the regularization type or strength, or specify the solver to use for objective-function minimization. If you display t in the Command Window, then all options appear empty ([]) except options that you specify using name-value pair arguments. During training, the software uses default values for empty options.
Examples: Train Multiclass Linear Classification Model Train an ECOC model composed of multiple binary, linear classification
models. Load the NLP data set. load nlpdata X is a sparse matrix of predictor data, and Y is a categorical vector of class labels. There are more than two classes in the data. Create a default linear-classification-model template. t = templateLinear(); To adjust the default values, see templateLinear. Train an ECOC model composed of multiple binary, linear classification models that can identify the product given the frequency distribution of words on a documentation web page. For faster training time, transpose the predictor data, and specify that observations correspond to columns. X = X'; rng(1); % For reproducibility Mdl = fitcecoc(X,Y,'Learners',t,'ObservationsIn','columns') Mdl = classreg.learning.classif.CompactClassificationECOC ResponseName: 'Y' ClassNames: [1×13 categorical] ScoreTransform: 'none' BinaryLearners: {78×1 cell} CodingMatrix: [13×78 double] Alternatively, you can train an ECOC model composed of default linear classification models using 'Learners','Linear'. To conserve memory, fitcecoc returns trained ECOC models composed of linear classification learners in CompactClassificationECOC model objects.
2.3.7 fitcecoc Fit multiclass models for support vector machines or other classifiers
Syntax Mdl = fitcecoc(Tbl,ResponseVarName) Mdl = fitcecoc(Tbl,formula) Mdl = fitcecoc(Tbl,Y) Mdl = fitcecoc(X,Y) Mdl = fitcecoc(___,Name,Value) [Mdl,HyperparameterOptimizationResults] = fitcecoc(___,Name,Value)
Description Mdl = fitcecoc(Tbl,ResponseVarName) returns a full, trained, multiclass, error-correcting output codes (ECOC) model using the predictors in table Tbl and the class labels in Tbl.ResponseVarName. fitcecocuses K(K – 1)/2 binary support vector machine (SVM) models using the one-versus-one coding design, where K is the number of unique class labels (levels). Mdl is a ClassificationECOC model. Mdl = fitcecoc(Tbl,formula) returns an ECOC model using the predictors in table Tbl and the class labels. formula is an explanatory model of the response and a subset of predictor variables in Tbl used for training. Mdl = fitcecoc(Tbl,Y) returns an ECOC model using the predictors in table Tbl and the class labels in vector Y. Mdl = fitcecoc(X,Y) returns a full, trained ECOC model using the predictors X and the class labels Y. Mdl = fitcecoc(___,Name,Value) returns an ECOC model with additional options specified by one or more Name,Value pair arguments, using any of the previous syntaxes. For example, specify different binary learners, a different coding design, or to cross-validate. It is good practice to cross-validate using the Kfold Name,Value pair argument. The cross-validation results
determine how well the model generalizes. [Mdl,HyperparameterOptimizationResults] = fitcecoc(___,Name,Value) also returns hyperparameter optimization details when you pass an OptimizeHyperparameters name-value pair with Learners = 'linear'. For other Learners, the HyperparameterOptimizationResults property of Mdlcontains the results.
Examples: Train Multiclass Model Using SVM Learners Train an error-correcting output codes (ECOC) multiclass model using support vector machine (SVM) binary learners. Load Fisher's iris data set. load fisheriris X = meas; Y = species; Train an ECOC multiclass model using the default options. Mdl = fitcecoc(X,Y) Mdl = ClassificationECOC ResponseName: 'Y' CategoricalPredictors: [] ClassNames: {'setosa' 'versicolor' 'virginica'} ScoreTransform: 'none' BinaryLearners: {3×1 cell} CodingName: 'onevsone' Mdl is a ClassificationECOC model. By default, fitcecoc uses SVM binary learners, and uses a one-versus-one coding design. You can access Mdl properties using dot notation. Display the coding design matrix.
Mdl.ClassNames CodingMat = Mdl.CodingMatrix ans = 3×1 cell array 'setosa' 'versicolor' 'virginica' CodingMat = 1 1 -1 0 0 -1
0 1 -1
A one-versus-one coding design on three classes yields three binary learners. Columns of CodingMat correspond to learners, and rows correspond to classes. The class order corresponds to the order in Mdl.ClassNames. For example, CodingMat(:,1) is [1; -1; 0] and indicates that the software trains the first SVM binary learner using all observations classified as 'setosa' and 'versicolor'. Since 'setosa' corresponds to 1, it is the positive class, and since 'versicolor' corresponds to -1, it is the negative class. You can access each binary learner using cell indexing and dot notation. Mdl.BinaryLearners{1}
% The first binary learner
ans = classreg.learning.classif.CompactClassificationSVM ResponseName: 'Y' CategoricalPredictors: [] ClassNames: [-1 1] ScoreTransform: 'none' Beta: [4×1 double]
Bias: 1.4505 KernelParameters: [1×1 struct] Compute the in-sample classification error. isLoss = resubLoss(Mdl) isLoss = 0.0067 The classification error is small, but the classifier might have been overfit. You can cross-validate the classifier using crossval.
Cross-Validate ECOC Classifier Load Fisher's iris data set. load fisheriris X = meas; Y = species; rng(1); % For reproducibility Create an SVM template. It is good practice to standardize the predictors. t = templateSVM('Standardize',1) t= Fit template for classification SVM. Alpha: [0×1 double] BoxConstraint: [] CacheSize: [] CachingMethod: '' ClipAlphas: [] DeltaGradientTolerance: [] Epsilon: [] GapTolerance: []
KKTTolerance: [] IterationLimit: [] KernelFunction: '' KernelScale: [] KernelOffset: [] KernelPolynomialOrder: [] NumPrint: [] Nu: [] OutlierFraction: [] RemoveDuplicates: [] ShrinkagePeriod: [] Solver: '' StandardizeData: 1 SaveSupportVectors: [] VerbosityLevel: [] Version: 2 Method: 'SVM' Type: 'classification' t is an SVM template. All of its properties are empty, except for StandardizeData, Method, and Type. When the software trains the ECOC classifier, it sets the applicable properties to their default values. Train the ECOC classifier. It is good practice to specify the class order. Mdl = fitcecoc(X,Y,'Learners',t,... 'ClassNames',{'setosa','versicolor','virginica'}); Mdl is a ClassificationECOC classifier. You can access its properties using dot notation. Cross-validate Mdl using 10-fold cross-validation. CVMdl = crossval(Mdl); CVMdl is a ClassificationPartitionedECOC cross-validated ECOC classifier. Estimate the generalization error. oosLoss = kfoldLoss(CVMdl)
oosLoss = 0.0400 The out-of-sample classification error is 4%, which indicates that the ECOC classifier generalizes fairly well.
Estimate Posterior Classifiers
Probabilities
Using
ECOC
Load Fisher's iris data set. Train the classifier using the petal dimensions as predictors. load fisheriris X = meas(:,3:4); Y = species; rng(1); % For reproducibility Create an SVM template, and specify the Gaussian kernel. It is good practice to standardize the predictors. t = templateSVM('Standardize',1,'KernelFunction','gaussian'); t is an SVM template. Most of its properties are empty. When the software trains the ECOC classifier, it sets the applicable properties to their default values. Train the ECOC classifier using the SVM template. Transform classification scores to class posterior probabilities (which are returned by predict or resubPredict) using the 'FitPosterior' name-value pair argument. Display diagnostic messages during the training using the 'Verbose' name-value pair argument. It is good practice to specify the class order. Mdl = fitcecoc(X,Y,'Learners',t,'FitPosterior',1,... 'ClassNames',{'setosa','versicolor','virginica'},... 'Verbose',2); Training binary learner 1 (SVM) out of 3 with 50 negative and 50 positive
observations. Negative class indices: 2 Positive class indices: 1 Fitting posterior probabilities for learner 1 (SVM). Training binary learner 2 (SVM) out of 3 with 50 negative and 50 positive observations. Negative class indices: 3 Positive class indices: 1 Fitting posterior probabilities for learner 2 (SVM). Training binary learner 3 (SVM) out of 3 with 50 negative and 50 positive observations. Negative class indices: 3 Positive class indices: 2 Fitting posterior probabilities for learner 3 (SVM). Mdl is a ClassificationECOC model. The same SVM template applies to each binary learner, but you can adjust options for each binary learner by passing in a cell vector of templates. Predict the in-sample labels and class posterior probabilities. Display diagnostic messages during the computation of labels and class posterior probabilities using the 'Verbose' name-value pair argument. [label,~,~,Posterior] = resubPredict(Mdl,'Verbose',1); Mdl.BinaryLoss Predictions from all learners have been computed. Loss for all observations has been computed. Computing posterior probabilities... ans = 'quadratic' The software assigns an observation to the class that yields the smallest
average binary loss. Since all binary learners are computing posterior probabilities, the binary loss function is quadratic. Display a random set of results. idx = randsample(size(X,1),10,1); Mdl.ClassNames table(Y(idx),label(idx),Posterior(idx,:),... 'VariableNames',{'TrueLabel','PredLabel','Posterior'}) ans = 3×1 cell array 'setosa' 'versicolor' 'virginica' ans = 10×3 table TrueLabel PredLabel Posterior ____________ ____________ ______________________________________ 'virginica' 'virginica' 'virginica' 'virginica' 'virginica' 'virginica' 'versicolor' 'versicolor' 'setosa' 'setosa' 'versicolor' 'virginica' 'versicolor' 'versicolor' 'setosa' 'setosa' 'versicolor' 'versicolor' 'setosa' 'setosa'
0.0039321 0.0039869 0.99208 0.017067 0.018263 0.96467 0.014948 0.015856 0.9692 2.2197e-14 0.87317 0.12683 0.999 0.00025091 0.00074639 2.2195e-14 0.059429 0.94057 2.2194e-14 0.97001 0.029986 0.999 0.0002499 0.00074741 0.0085646 0.98259 0.008849 0.999 0.00025013 0.00074718
The columns of Posterior correspond to the class order of Mdl.ClassNames. Define a grid of values in the observed predictor space. Predict the
posterior probabilities for each instance in the grid. xMax = max(X); xMin = min(X); x1Pts = linspace(xMin(1),xMax(1)); x2Pts = linspace(xMin(2),xMax(2)); [x1Grid,x2Grid] = meshgrid(x1Pts,x2Pts); [~,~,~,PosteriorRegion] = predict(Mdl,[x1Grid(:),x2Grid(:)]); For each coordinate on the grid, plot the maximum class posterior probability among all classes. figure; contourf(x1Grid,x2Grid,... reshape(max(PosteriorRegion,[],2),size(x1Grid,1),size(x1Grid,2))); h = colorbar; h.YLabel.String = 'Maximum posterior'; h.YLabel.FontSize = 15; hold on gh = gscatter(X(:,1),X(:,2),Y,'krk','*xd',8); gh(2).LineWidth = 2; gh(3).LineWidth = 2; title 'Iris Petal Measurements and Maximum Posterior'; xlabel 'Petal length (cm)'; ylabel 'Petal width (cm)'; axis tight legend(gh,'Location','NorthWest') hold off
Train Multiclass Linear Classification Model Train an ECOC model composed of multiple binary, linear classification models. Load the NLP data set. load nlpdata X is a sparse matrix of predictor data, and Y is a categorical vector of class labels. There are more than two classes in the data. Create a default linear-classification-model template. t = templateLinear(); To adjust the default values, see templateLinear. Train an ECOC model composed of multiple binary, linear classification models that can identify the product given the frequency distribution of words on a documentation web page. For faster training time, transpose the predictor data, and specify that observations correspond to columns. X = X'; rng(1); % For reproducibility
Mdl = fitcecoc(X,Y,'Learners',t,'ObservationsIn','columns') Mdl = classreg.learning.classif.CompactClassificationECOC ResponseName: 'Y' ClassNames: [1×13 categorical] ScoreTransform: 'none' BinaryLearners: {78×1 cell} CodingMatrix: [13×78 double] Alternatively, you can train an ECOC model composed of default linear classification models using 'Learners','Linear'. To conserve memory, fitcecoc returns trained ECOC models composed of linear classification learners in CompactClassificationECOC model objects.
Train ECOC Classifiers Using Ensembles and Parallel Computing Train a one-versus-all ECOC classifier using a GentleBoost ensemble of decision trees with surrogate splits. Estimate the classification error using 10-fold cross-validation. Load and inspect the arrhythmia data set. load arrhythmia [n,p] = size(X) isLabels = unique(Y); nLabels = numel(isLabels) tabulate(categorical(Y)) n= 452 p= 279
nLabels = 13 Value 1 2 3 4 5 6 7 8 9 10 14 15 16
Count Percent 245 54.20% 44 9.73% 15 3.32% 15 3.32% 13 2.88% 25 5.53% 3 0.66% 2 0.44% 9 1.99% 50 11.06% 4 0.88% 5 1.11% 22 4.87%
There are 279 predictors, and a relatively small sample size of 452. There are 16 distinct labels, but only 13 are represented in the response (Y), and each label describes various degrees of arrhythmia. 54.20% of the observations are in class 1. Create an ensemble template. You must specify at least three arguments: a method, a number of learners, and the type of learner. For this example, specify 'GentleBoost' for the method, 100 for the number of learners, and a decision tree template that uses surrogate splits since there are missing observations. tTree = templateTree('surrogate','on'); tEnsemble = templateEnsemble('GentleBoost',100,tTree); tEnsemble is a template object. Most of its properties are empty, but the software fills them with their default values during training. Train a one-versus-all ECOC classifier using the ensembles of decision
trees as binary learners. If you have a Parallel Computing Toolbox™ license, then you can speed up the computation by specifying to use parallel computing. This sends each binary learner to a worker in the pool (the number of workers depends on your system configuration). Also, specify that the prior probabilities are 1/K, where K = 13, which is the number of distinct classes. pool = parpool; % Invoke workers options = statset('UseParallel',1); Mdl = fitcecoc(X,Y,'Coding','onevsall','Learners',tEnsemble,... 'Prior','uniform','Options',options); Starting parallel pool (parpool) using the 'local' profile ... connected to 2 workers. Mdl is a ClassificationECOC model. Cross-validate the ECOC classifier using 10-fold cross-validation. CVMdl = crossval(Mdl,'Options',options); Warning: One or more folds do not contain points from all the groups. CVMdl is a ClassificationPartitionedECOC model. The warning indicates that some classes are not represented while the software trains at least one fold. Therefore, those folds cannot predict labels for the missing classes. You can inspect the results of a fold using cell indexing and dot notation, e.g., access the results of the first fold by entering CVMdl.Trained{1}. Your results might vary. Use the cross-validated ECOC classifier to predict out-of-fold labels. You can compute the confusion matrix using confusionmat. However, if you have a Neural Network Toolbox™ license, you can plot the confusion matrix using plotconfusion. The input arguments of plotconfusion are not vectors of the true and predicted labels like confusionmat, but indicator matrices of the true and predicted labels. Both start as K-by-n matrices of 0s. If observation j has label index k (or has predicted label k), then element (k,j) of the true label indicator matrix (or predicted label indicator matrix) is 1. You can convert label indices returned by predict, resubPredict, or kfoldPredict to label indicator matrices using linear indexing. For details
on linear indexing, see sub2ind and ind2sub. oofLabel = kfoldPredict(CVMdl,'Options',options); ConfMat = confusionmat(Y,oofLabel); % Convert the integer label vector to a class-identifier matrix. [~,grpOOF] = ismember(oofLabel,isLabels); oofLabelMat = zeros(nLabels,n); idxLinear = sub2ind([nLabels n],grpOOF,(1:n)'); oofLabelMat(idxLinear) = 1; % Flags the row corresponding to the class [~,grpY] = ismember(Y,isLabels); YMat = zeros(nLabels,n); idxLinearY = sub2ind([nLabels n],grpY,(1:n)'); YMat(idxLinearY) = 1; figure; plotconfusion(YMat,oofLabelMat); h = gca; h.XTickLabel = [num2cell(isLabels); {''}]; h.YTickLabel = [num2cell(isLabels); {''}];
Optimize ECOC Classifier This example shows how to optimize hyperparameters automatically using fitcecoc. The example uses Fisher's iris data. Load the data.
load fisheriris X = meas; Y = species;
Find hyperparameters that minimize five-fold cross-validation loss by using automatic hyperparameter optimization. For reproducibility, set the random seed and use the 'expectedimprovement-plus' acquisition function.
rng default
Mdl = fitcecoc(X,Y,'OptimizeHyperparameters','auto',
'HyperparameterOptimizationOptions',struct('Ac 'expected-improvement-plus')) |================================================================================================================|
| Iter | Eval KernelScale |
|
| Objective
| result | |
| Objective
| BestSoFar
| BestSoFar
| runtime
| (observed) | (estim.)
|
|
Coding | BoxConstrain-|
| t
|
|================================================================================================================|
| |
1 | Best | 310.42 |
0.33333 |
3.3156 |
0.33333 |
0.33333 |
onevsall |
0.0089562
| 2 | Best 0.0053534 |
|
0.053333 |
5.2976 |
0.053333 |
0.073138 |
onevsone |
67.106 |
| 3 | Best 0.017254 |
|
0.04 |
0.57165 |
0.04 |
0.040121 |
onevsall |
0.0027084 |
| 4 | Accept | 0.14755 |
0.06 |
13.097 |
0.04 |
0.040005 |
onevsall |
259.06 |
| |
5 | Accept | 1.0303 |
0.046667 |
0.66074 |
0.04 |
0.039999 |
onevsone |
127.23
| |
6 | Accept | 996.84 |
0.11333 |
0.41791 |
0.04 |
0.04 |
onevsone |
0.37412
| 7 | Accept | 0.0010002 |
0.48 |
20.956 |
0.04 |
0.040038 |
onevsall |
182.3 |
| 8 | Accept | 0.042615 |
0.04 |
0.76888 |
0.04 |
0.039982 |
onevsall |
0.090529 |
| 9 | Best 0.09393 |
|
0.02 |
0.54564 |
0.02 |
0.020009 |
onevsone |
0.0030779 |
|
10 | Accept | 0.02549 |
0.046667 |
0.41283 |
0.02 |
0.020015 |
onevsone |
1.9001 |
| |
11 | Accept | 24.057 |
0.11333 |
0.45093 |
0.02 |
0.020016 |
onevsone |
0.082855
| 12 | Accept | 0.26902 |
0.02 |
0.55456 |
0.02 |
0.019995 |
onevsone |
0.028965 |
| 13 | Accept | 0.0010042 |
0.046667 |
5.7777 |
0.02 |
0.020004 |
onevsone |
74.47 |
14 | Accept | 0.2665 |
0.33333 |
0.37761 |
0.02 |
0.020001 |
onevsall |
| 15 | Accept | 0.04193 |
0.26 |
19.695 |
0.02 |
0.020003 |
onevsall |
971.54 |
| 16 | Accept | 0.82502 |
0.04 |
0.82597 |
0.02 |
0.020003 |
onevsall |
83.197 |
| 17 | Accept | 0.0024331 |
0.046667 |
0.5434 |
0.02 |
0.020004 |
onevsone |
0.028079 |
| 18 | Accept | 0.0019763 |
0.04 |
5.9089 |
0.02 |
0.020207 |
onevsall |
0.0011615 |
| 19 | Accept | 0.0073833 |
0.04 |
1.8263 |
0.02 |
0.020004 |
onevsall |
0.020594 |
| |
0.0010257
| 20 | Accept | 0.19642 |
0.04 |
1.3375 |
0.02 |
0.020004 |
onevsall |
9.6782 |
|================================================================================================================|
| Iter | Eval KernelScale |
| |
| Objective
| result | |
| Objective
| BestSoFar
| BestSoFar
| runtime
| (observed) | (estim.)
|
|
Coding | BoxConstrain-|
| t
|================================================================================================================|
| |
21 | Accept | 39.212 |
0.02 |
0.3982 |
0.02 |
0.020005 |
onevsone |
965.24
| |
22 | Accept | 643.21 |
0.12 |
0.47837 |
0.02 |
0.020005 |
onevsone |
980.8
| |
23 | Accept | 8.1701 |
0.04 |
0.41142 |
0.02 |
0.020037 |
onevsone |
973.21
| 24 | Accept | 0.69618 |
0.11333 |
0.4301 |
0.02 |
0.020026 |
onevsone |
0.0010507 |
| 25 | Accept | 0.012414 |
0.033333 |
0.46536 |
0.02 |
0.020011 |
onevsone |
0.0011494 |
| 26 | Accept | 0.058611 |
0.033333 |
0.47349 |
0.02 |
0.020506 |
onevsone |
0.035003 |
| 27 | Accept | 0.34837 |
0.033333 |
0.3971 |
0.02 |
0.020059 |
onevsone |
1.2763 |
| 28 | Accept | 0.096199 |
0.046667 |
2.4195 |
0.02 |
0.020367 |
onevsone |
978.94 |
| |
29 | Accept | 34.743 |
0.046667 |
0.36933 |
0.02 |
0.020471 |
onevsone |
52.349
| 30 | Accept | 0.65742 |
0.046667 |
7.0458 |
0.02 |
0.022071 |
onevsall |
984.41 |
______________________________________________ Optimization completed. MaxObjectiveEvaluations of 30 reached. Total function evaluations: 30 Total elapsed time: 160.6519 seconds. Total objective function evaluation time: 96.2299
Best observed feasible point: Coding KernelScale
BoxConstraint
________ ___________
_____________
onevsone
0.0030779
0.09393
Observed objective function value = 0.02 Estimated objective function value = 0.022071 Function evaluation time = 0.54564
Best estimated feasible point (according to models): Coding KernelScale
BoxConstraint
________ ___________
_____________
onevsone 0.09393
0.0030779
Estimated objective function value = 0.022071 Estimated function evaluation time = 0.49997
Mdl =
ClassificationECOC ResponseName: 'Y' CategoricalPredictors: []
{'setosa'
'versicolor'
ClassNames: 'virginica'}
ScoreTransform: 'none' BinaryLearners: {3×1 cell} CodingName: 'onevsone' HyperparameterOptimizationResults: [1×1 BayesianOptimization]
2.3.8 templateECOC Error-correcting output codes learner template
Syntax t = templateECOC() t = templateECOC(Name,Value)
Description t = templateECOC() returns an error-correcting output codes (ECOC) classification learner template. If you specify a default template, then the software uses default values for all input arguments during training. t = templateECOC(Name,Value) returns a template with additional options specified by one or more name-value pair arguments. For example, you can specify a coding design, whether to fit posterior probabilities, or the types of binary learners. If you display t in the Command Window, then all options appear empty ([]), except those that you specify using name-value pair arguments. During training, the software uses default values for empty options.
Examples Create a Default ECOC Classification Learner Template Use templateECOC to create a default ECOC template. t = templateECOC() t= Fit template for classification ECOC. BinaryLearners: '' Coding: ''
FitPosterior: [] Options: [] VerbosityLevel: [] Version: 1 Method: 'ECOC' Type: 'classification' All properties of the template object are empty except for Method and Type. When you pass t to testckfold, the software fills in the empty properties with their respective default values. For example, the software fills the BinaryLearners property with 'SVM'. For details on other default values, see fitcecoc. t is a plan for an ECOC learner. When you create it, no computation occurs. You can pass t to testckfold to specify a plan for an ECOC classification model to statistically compare with another model.
Statistically Compare Performance of Two ECOC Classification Models One way to select predictors or features is to train two models where one that uses a subset of the predictors that trained the other. Statistically compare the predictive performances of the models. If there is sufficient evidence that model trained on fewer predictors performs better than the model trained using more of the predictors, then you can proceed with a more efficient model. Load Fisher's iris data set. Plot all 2-dimensional combinations of predictors. load fisheriris d = size(meas,2); % Number of predictors pairs = combnk(1:d,2); figure; for j = 1:size(pairs,1); subplot(3,2,j); gscatter(meas(:,pairs(j,1)),meas(:,pairs(j,2)),species);
xlabel(sprintf('meas(:,%d)',pairs(j,1))); ylabel(sprintf('meas(:,%d)',pairs(j,2))); legend off; end
Based on the scatterplot, meas(:,3) and meas(:,4) seem like they separate the groups well. Create an ECOC template. Specify to use a one-versus-all coding design. t = templateECOC('Coding','onevsall'); By default, the ECOC model uses linear SVM binary learners. You can choose other, supported algorithms by specifying them using the 'Learners' name-value pair argument. Test whether an ECOC model that is just trained using predictors 3 and 4 performs at most as well as an ECOC model that is trained using all predictors. Rejecting this null hypothesis means that the ECOC model trained using predictors 3 and 4 performs better than the ECOC model trained using all predictors. Suppose represents the classification error of the ECOC model trained using predictors 3 and 4 and represents the classification error of the
ECOC model trained using all predictors, then the test is:
By default, testckfold conducts a 5-by-2 k-fold F test, which is not appropriate as a one-tailed test. Specify to conduct a 5-by-2 k-fold t test. rng(1); % For reproducibility [h,pValue] = testckfold(t,t,meas(:,pairs(1,:)),meas,species,... 'Alternative','greater','Test','5x2t') h= logical 0 pValue = 0.8940 The h = 0 indicates that there is not enough evidence to suggest that the model trained using predictors 3 and 4 is more accurate than the model trained using all predictors.
2.4 TRAIN SUPPORT VECTOR MACHINES USING CLASSIFICATION LEARNER APP This example shows how to construct support vector machine (SVM) classifiers in the Classification Learner app, using the ionosphere data set that contains two classes. You can use a support vector machine (SVM) with two or more classes in Classification Learner. An SVM classifies data by finding the best hyperplane that separates all data points of one class from those of another class. In the ionosphere data, the response variable is categorical with two levels: g represents good radar returns, and b represents bad radar returns. 1. In MATLAB®, load the ionosphere data set and define some variables from the data set to use for a classification. 2. load ionosphere 3. ionosphere = array2table(X); ionosphere.Group = Y; 4. On the Apps tab, in the Math, Statistics and Optimization group, click Classification Learner. 5. On the Classification Learner tab, in the File section, click New Session. In the New Session dialog box, observe that the app has selected response and predictor variables based on their data type. The response variable Group has two levels. All the other variables are predictors. 6. Click Start Session. Classification Learner creates a scatter plot of the data. 7. Use the scatter plot to visualize which variables are useful for predicting the response. Select different variables in the X- and Y-axis controls. Observe which variables separate the class colors most clearly. 8. To create a selection of SVM models, on the Classification Learner tab, in the Model Type section, click the down arrow to expand the list of classifiers, and under Support Vector Machines, click All SVMs.
Then click Train. Tip: If you have Parallel Computing Toolbox™ then the first time you click Train you see a dialog while the app opens a parallel pool of workers. After the pool opens, you can train multiple classifiers at once and continue working.
Classification Learner trains one of each SVM classification option in the gallery, and highlights the best score. The app outlines in a box the Accuracy score of the best model.
9. Select a model in the History list to view the results. Examine the scatter plot for the trained model and try plotting different predictors. Misclassified points are shown as an X. 10. To inspect the accuracy of the predictions in each class, on the Classification Learner tab, in the Plots section, click Confusion
Matrix. View the matrix of true class and predicted class results. 11. Select the other models in the list to compare. 12. Choose the best model in the History list (the best score is highlighted in a box). To improve the model, try including different features in the model. See if you can improve the model by removing features with low predictive power. On the Classification Learner tab, in the Data Transformation section, click Features. In the Feature Selection dialog box, specify predictors to remove from the model, and click Train to train a new model using the new options. Compare results among the classifiers in the History list. 13. To investigate features to include or exclude, use the parallel coordinates plot. On the Classification Learner tab, in the Plots section, select Parallel Coordinates Plot. 14. Choose the best model in the History list. To try to improve the model further, try changing SVM settings. On the Classification Learner tab, in the Model Type section, click Advanced. Try changing a setting, then train the new model by clicking Train. 15. To export the trained model to the workspace, select the Classification Learner tab and click Export model. 16. To examine the code for training this classifier, select Export Model > Generate Code. Use the same workflow to evaluate and compare the other classifier types you can train in Classification Learner. To try all the classifier model presets available for your data set: 1. Click the arrow on the far right of the Model Type section to expand the list of classifiers. 2. Click All, then click Train.
2.5 TRAIN CLASSIFICATION MODELS IN CLASSIFICATION LEARNER APP
2.5.1 What Is Supervised Machine Learning? You can use Classification Learner to train models to classify data. Using this app, you can explore supervised machine learning using various classifiers. You can explore your data, select features, specify validation schemes, train models, and assess results. You can perform automated training to search for the best classification model type, including decision trees, discriminant analysis, support vector machines, logistic regression, nearest neighbors, and ensemble classification. Perform supervised machine learning by supplying a known set of input data (observations or examples) and known responses to the data (i.e., labels or classes). Use the data to train a model that generates predictions for the response to new data. To use the model with new data, or to learn about programmatic classification, you can export the model to the workspace or generate MATLAB® code to recreate the trained model.
Get started by training a selection of model types.
2.5.2 Automated Classifier Training You can use Classification Learner to automatically train a selection of different classification models on your data. ·
· 1. 2. 3.
Get started by automatically training multiple models at once. You can quickly try a selection of models, then explore promising models interactively. If you already know what classifier type you want, train individual classifiers instead. On the Apps tab, in the Math, Statistics and Optimization group, click Classification Learner. Click New Session and select data from the workspace or from file. Specify a response variable and variables to use as predictors. On the Classification Learner tab, in the Model Type section, click All Quick-To-Train. This option will train all the model presets available for your data set that are fast to fit.
4. Click Train . 5. A selection of model types appears in the History list. When they finish training, the best percentage Accuracy score is highlighted in a box.
6. Click models in the history list to explore results in the plots. 7. To try all the classifier model presets available for your data set, click All, then click Train.
2.5.3 Manual Classifier Training If you want to explore individual model types, or if you already know what classifier type you want, you can train classifiers one at a time, or a train a group of the same type. 1. Choose a classifier. On the Classification Learner tab, in the Model Type section, click a classifier type. To see all available classifier options, click the arrow on the far right of the Model Type section to expand the list of classifiers. The options in the Model Type gallery are preset starting points with different settings, suitable for a range of different classification problems. To read a description of each classifier, switch to the details view.
For more information on each option, see Choose Classifier Options.
2. After selecting a classifier, click Train. Repeat to try different classifiers. 3. If you decide you want to try all model types or train a group of the same type, try one of the All options in the Model Type gallery.
2.5.4 Parallel Classifier Training You can train models in parallel using Classification Learner if you have Parallel Computing Toolbox. When you train classifiers, the app automatically starts a parallel pool of workers, unless you turn off the default parallel preference Automatically create a parallel pool. If a pool is already open, the app uses it for training. Parallel training allows you to train multiple classifiers at once and continue working. 1. The first time you click Train, you see a dialog while the app opens a parallel pool of workers. After the pool opens, you can train multiple classifiers at once. 2. When classifiers are training in parallel, you see progress indicators on each training and queued model in the history list, and you can cancel individual models if you want. During training, you can examine results and plots from models, and initiate training of more classifiers.
To control parallel training, toggle the Use Parallel button in the app toolstrip. The Use Parallel button is only available if you have Parallel Computing Toolbox.
If you have Parallel Computing Toolbox, then parallel training is available in Classification Learner, and you do not need to set the UseParallel option of the statset function. If you turn off the parallel preference to Automatically create a parallel pool, then the app will not start a pool for you without asking first.
2.5.5 Compare and Improve Classification Models 1. Click models in the history list to explore the results in the plots. Compare model performance by inspecting results in the scatter plot and confusion matrix. Examine the percentage accuracy reported in the history list for each model. 2. Select the best model in the history list and then try including and excluding different features in the model. Click Feature Selection.
Try the parallel coordinates plot to help you identify features to remove. See if you can improve the model by removing features with low predictive power. Specify predictors to include in the model, and train new models using the new options. Compare results among the models in the history list. You can also try transforming features with PCA to reduce
dimensionality. 3. To improve the model further, you can try changing classifier parameter settings in the Advanced dialog box, and then train using the new options. 4. If feature selection, PCA, or new parameter settings improve your model, try training All model types with the new settings. See if another model type does better with the new settings. The figure shows the app with a history list containing various classifier types.
2.6 CHOOSE CLASSIFIER OPTIONS 2.6.1 Choose a Classifier Type You can use Classification Learner to automatically train a selection of different classification models on your data. Use automated training to quickly try a selection of model types, then explore promising models interactively. To get started, try these options first:
Get Started Classifier Buttons All Q uick-ToTrain All Linear All
Description Try this first. The app will train all the model types available for your data set that are typically fast to fit. Try this if you expect linear boundaries between the classes in your data. This option fits only Linear SVM and Linear Discriminant. Use this to train all available model types. Trains every type regardless of any prior trained models. Can be time-consuming.
If you want to explore classifiers one at a time, or you already know what classifier type you want, you can select individual models or train a group of the same type. To see all available classifier options, on the Classification Learner tab, click the arrow on the far right of the Model Type section to expand the list of classifiers. The options in the Model Type gallery are preset starting points with different settings, suitable for a range of different classification problems. For help choosing the best classifier type for your problem, see the table showing typical characteristics of different supervised learning algorithms. Use the table as a guide for your final choice of algorithms. Decide on the tradeoff you want in speed, memory usage, flexibility, and interpretability. The best classifier type depends on your data. Characteristics of Classifier Types
Classifier Decision Trees
Prediction Speed
Memory Usage
Interpretability
Fast
Small
Easy
Fast
Small for linear, large for quadratic
Easy
Fast
Medium
Easy
Medium for linear. Slow for others.
Medium for linear. All others: medium for multiclass, large for binary.
Easy for Linear SVM. Hard for all other kernel types.
Slow for cubic. Medium for others.
Medium
Hard
Fast to medium depending on choice of algorithm
Low to high depending on choice of algorithm
Hard
Discriminant Analysis
Logistic Regression
Support Vector Machines
Note: Supports C code generation for prediction. Nearest Neighbor Classifiers
Ensemble Classifiers
The tables on this page describe general characteristics of speed and memory usage for all the preset classifiers. The classifiers were tested with various data sets (up to 7000 observations, 80 predictors, and 50 classes), and the results define the following groups: Speed · Fast 0.01 second · Medium 1 second · Slow 100 seconds Memory · Small 1MB · Medium 4MB · Large 100MB These tables provide a general guide. Your results depend on your data and the speed of your machine.
To read a description of each classifier in Classification Learner, switch to the details view.
Categorical Predictor Support In Classification Learner, the Model Type gallery only shows classifier types that support your selected data. Classifier Decision Trees Discriminant Analysis Logistic Regression SVM Nearest Neighbor Ensembles
All predictors numeric Yes
Yes
Some categorical, some numeric Yes
Yes
No
No
Yes Yes Euclidean distance only
Yes Yes
Yes Yes
Hamming distance only
No
Yes, except Subspace Discriminant
Yes, except any Subspace
Yes
All predictors categorical
2.6.2 Decision Trees Decision trees are easy to interpret, fast for fitting and prediction, and low on memory usage, but they can have low predictive accuracy. Try to grow simpler trees to prevent overfitting. Control the depth with the Maximum number of splits setting.
Classifier Type Simple Tree
Prediction Speed
Memory Usage
Interpretability
Fast
Small
Easy
Fast
Small
Easy
Fast
Small
Easy
Model Flexibility Low.
Medium Tree
Few leaves to make coarse distinctions between classes (maximum number of splits is 4).
Medium
Complex Tree
Medium number of leaves for finer distinctions between classes (maximum number of splits is 20).
High Many leaves to make many fine distinctions between classes (maximum number of splits is 100).
You train classification trees to predict responses to data. To predict a response, follow the decisions in the tree from the root (beginning) node down to a leaf node. The leaf node contains the response. Statistics and Machine Learning Toolbox™ trees are binary. Each step in a prediction involves checking the value of one predictor (variable). For example, here is a simple classification tree:
This tree predicts classifications based on two predictors, x1 and x2. To
predict, start at the top node. At each decision, check the values of the predictors to decide which branch to follow. When the branches reach a leaf node, the data is classified either as type 0 or 1. You can visualize your decision tree model by exporting the model from the app, and then entering: view(trainedModel.ClassificationTree,'Mode','graph') The figure shows an example complex tree trained with the fisheriris data.
Advanced Tree Options Classification trees in Classification Learner use the fitctree function. You can set these options: Maximum number of splits Specify the maximum number of splits or branch points to control the depth of your tree. When you grow a decision tree, consider its simplicity and predictive power. To change the number of splits, click the buttons or enter
a positive integer value in the Maximum number of splits box. · A complex tree with many leaves is usually highly accurate on the training data. However, the tree might not show comparable accuracy on an independent test set. A leafy tree tends to overtrain, and its validation accuracy is often far lower than its training (or resubstitution) accuracy. · In contrast, a simple tree does not attain high training accuracy. But a simple tree can be more robust in that its training accuracy can approach that of a representative test set. Also, a simple tree is easy to interpret. Split criterion Specify the split criterion measure for deciding when to split nodes. Try each of the three settings to see if they improve the model with your data. Split criterion options are Gini's diversity index, Twoing rule, or Maximum deviance reduction (also known as cross entropy). The classification tree tries to optimize to pure nodes containing only one class. Gini's diversity index (the default) and the deviance criterion measure node impurity. The twoing rule is a different measure for deciding how to split a node, where maximizing the twoing rule expression increases node purity. For details of these split criteria, see ClassificationTree Definitions. Surrogate decision splits — Only for missing data. Specify surrogate use for decision splits. If you have data with missing values, use surrogate splits to improve the accuracy of predictions. When you set Surrogate decision splits to On, the classification tree finds at most 10 surrogate splits at each branch node. To change the number, click the buttons or enter a positive integer value in the Maximum surrogates per node box. When you set Surrogate decision splits to Find All, the classification tree finds all surrogate splits at each branch node. The Find All setting can use considerable time and memory.
2.6.3 Discriminant Analysis Discriminant analysis is a popular first classification algorithm to try because it is fast, accurate and easy to interpret. Discriminant analysis is good for wide datasets. Discriminant analysis assumes that different classes generate data based on different Gaussian distributions. To train a classifier, the fitting function estimates the parameters of a Gaussian distribution for each class. Classifier Type Linear Discriminant
Prediction Speed
Memory Usage
Interpretability
Model Flexibility
Fast
Small
Easy
Low. Creates linear boundaries between classes.
Fast
Large
Easy
Low. Creates nonlinear boundaries between classes (ellipse, parabola or hyperbola).
Quadratic Discriminant
Advanced Discriminant Options Discriminant analysis in Classification Learner uses the fitcdiscr function. For both linear and quadratic discriminants, you can change the Covariance structure option. If you have predictors with zero variance or if any of the covariance matrices of your predictors are singular, training can fail using the default, Full covariance structure. If training fails, select the Diagonal covariance structure instead.
2.6.4 Logistic Regression If you have 2 classes, logistic regression is a popular simple classification algorithm to try because it is easy to interpret. The classifier models the class probabilities as a function of the linear combination of predictors. Classifier Type Logistic Regression
Prediction Speed
Memory Usage
Interpretability
Model Flex ibility
Fast
Medium
Easy
Low. You cannot change any parameters to control model flexibility.
Logistic regression in Classification Learner uses the fitglm function. You cannot set any options for this classifier in the app.
2.6.5 Support Vector Machines In Classification Learner, you can train SVMs when your data has two or more classes. Classifier Type Linear SVM
Quadratic SVM
Cubic SVM
Fine Gaussian SVM
Medium Gaussian SVM
Coarse Gaussian SVM
Prediction Speed
Memory Usage
Interpretability
Model Flexibility
Binary: Fast Multiclass: Medium
Medium
Easy
Low Makes a simple linear separation between classes.
Binary: Fast Multiclass: Slow
Binary: Medium Multiclass: Large
Hard
Medium
Binary: Fast Multiclass: Slow
Binary: Medium Multiclass: Large
Hard
Medium
Binary: Fast Multiclass: Slow
Binary: Medium Multiclass: Large
Hard
High — decreases with kernel scale setting. Makes finely detailed distinctions between classes, with kernel scale set to sqrt(P)/4.
Binary: Fast Multiclass: Slow
Binary: Medium Multiclass: Large
Hard
Medium Medium distinctions, with kernel scale set to sqrt(P).
Binary: Fast Multiclass: Slow
Binary: Medium Multiclass: Large
Hard
Low Makes coarse distinctions between classes, with kernel scale set to sqrt(P)*4, where P is the number of predictors.
An SVM classifies data by finding the best hyperplane that separates data points of one class from those of the other class. The best hyperplane for an SVM means the one with the largest margin between the two classes. Margin means the maximal width of the slab parallel to the hyperplane that has no interior data points. The support vectors are the data points that are closest to the separating hyperplane; these points are on the boundary of the slab. The following figure illustrates these definitions, with + indicating data points of type 1, and – indicating data points of type –1.
SVMs can also use a soft margin, meaning a hyperplane that separates many, but not all data points. Advanced SVM Options If you have exactly two classes, Classification Learner uses the fitcsvm function to train the classifier. If you have more than two classes, the app uses the fitcecoc function to reduce the multiclass classification problem to a set of binary classification subproblems, with one SVM learner for each subproblem. To examine the code for the binary and multiclass classifier types, you can generate code from your trained classifiers in the app. You can set these options in the app: Kernel function Specify the Kernel function to compute the classifier. · Linear kernel, easiest to interpret · Gaussian or Radial Basis Function (RBF) kernel · Quadratic · Cubic Box constraint level Specify the box constraint to keep the allowable values of the Lagrange multipliers in a box, a bounded region. To tune your SVM classifier, try increasing the box constraint level. Click the buttons or enter a positive scalar value in the Box constraint level box.
Increasing the box constraint level can decrease the number of support vectors, but also can increase training time. The Box Constraint parameter is the soft-margin penalty known as C in the primal equations, and is a hard "box" constraint in the dual equations. Kernel scale mode Specify manual kernel scaling if desired. When you set Kernel scale mode to Auto, then the software uses a heuristic procedure to select the scale value. The heuristic procedure uses subsampling. Therefore, to reproduce results, set a random number seed using rng before training the classifier. When you set Kernel scale mode to Manual, you can specify a value. Click the buttons or enter a positive scalar value in the Manual kernel scale box. The software divides all elements of the predictor matrix by the value of the kernel scale. Then, the software applies the appropriate kernel norm to compute the Gram matrix. Multiclass method Only for data with 3 or more classes. This method reduces the multiclass classification problem to a set of binary classification subproblems, with one SVM learner for each subproblem. One-vs-One trains one learner for each pair of classes. It learns to distinguish one class from the other. Onevs-All trains one learner for each class. It learns to distinguish one class from all others. Standardize data Specify whether to scale each coordinate distance. If predictors have widely different scales, standardizing can improve the fit.
2.6.6 Nearest Neighbor Classifiers Nearest neighbor classifiers typically have good predictive accuracy in low dimensions, but might not in high dimensions. They have high memory usage, and are not easy to interpret.
Classifier Type Fine KNN
Prediction Speed
Memory Usage
Interpretability
Model Flexibility
Medium
Medium
Hard
Finely detailed distinctions between classes. The number of neighbors is set to 1.
Medium
Medium
Hard
Medium distinctions between classes. The number of neighbors is set to 10.
Medium
Medium
Hard
Coarse distinctions between classes. The number of neighbors is set to 100.
Medium
Medium
Hard
Medium distinctions between classes, using a Cosine distance metric. The number of neighbors is set to 10.
Slow
Medium
Hard
Medium distinctions between classes, using a cubic distance metric. The number of neighbors is set to 10.
Medium
Medium
Hard
Medium distinctions between classes, using a distance weight. The number of neighbors is set to 10.
Medium KNN
Coarse KNN
Cosine KNN
Cubic KNN
Weighted KNN
What is k-Nearest Neighbor classification? Categorizing query points based on their distance to points (or neighbors) in a training dataset can be a simple yet effective way of classifying new points. You can use various metrics to determine the distance. Given a set X of n points and a distance function, k-nearest neighbor (kNN) search lets you find the k closest points in X to a query point or set of points. kNN-based algorithms are widely used as benchmark machine learning rules.
Advanced KNN Options Nearest Neighbor classifiers in Classification Learner use the fitcknn function. You can set these options: Number of neighbors Specify the number of nearest neighbors to find for classifying each point when predicting. Specify a fine (low number) or coarse classifier (high number) by changing the number of neighbors. For example, a fine KNN uses one neighbor, and a coarse KNN uses 100. Many neighbors can be time consuming to fit. Distance metric You can use various metrics to determine the distance to points. For definitions, see the class ClassificationKNN. Distance weight Specify the distance weighting function. You can choose Equal (no weights), Inverse (weight is 1/distance), or Squared Inverse (weight is 1/distance2). Standardize data Specify whether to scale each coordinate distance. If predictors have widely different scales, standardizing can improve the fit.
2.6.7 Ensemble Classifiers Ensemble classifiers meld results from many weak learners into one highquality ensemble model. Qualities depend on the choice of algorithm.
Classifier Type
Prediction Speed
Memory Usage
Inter-pretability
Ensemble Method
Model Flexibility
Fast
Low
Hard
AdaBoost , with Decision Tree learners
Medium to high — increases with Number of learners or Maximum number of splits setting. Tip: Boosted trees can usually do better than bagged, but might require parameter tuning and more learners
Medium
High
Hard
Random forest Bag, with Decision Treelearners
High — increases with Number of learners setting. Tip: Try this classifier first.
Medium
Low
Hard
Subspace, with Discriminant learners
Medium — increases with Number of learners setting. Good for many predictors
Medium
Medium
Hard
Subspace, with Nearest Neighbor learners
Medium — increases with Number of learners setting. Good for many predictors
Hard
RUSBoost , with Decision Tree learners
Medium — increases with Number of learners or Maximum number of splits setting. Good for skewed data (with many more observations of 1 class)
Hard
GentleBoost or LogitBoost , with Decision Tree learners Choose Boosted Trees and change to GentleBoost method.
Medium — increases with Number of learners or Maximum number of splits setting. For binary classification only
Boosted Trees
Bagged Trees
Subspace Discriminant
Subspace KNN
RUSBoost Trees Fast
GentleBoost or LogitBoost — not available in the Model Type gallery. If you have 2 class data, select manually.
Fast
Low
Low
Bagged trees use Breiman's 'random forest' algorithm. For reference, see Breiman, L. Random Forests. Machine Learning 45, pp. 5–32, 2001. Tips ·
Try bagged trees first. Boosted trees can usually do better but might
require searching many parameter values, which is time-consuming. · Try training each of the ensemble classifier options in the Model Type gallery. Train them all to see which settings produce the best model with your data. Select the best model in the History list. To try to improve your model, try feature selection, PCA, and then (optionally) try changing some advanced options. · For boosting ensemble methods, you can get fine detail with either deeper trees or larger numbers of shallow trees. As with single tree classifiers, deep trees can cause overfitting. You need to experiment to choose the best tree depth for the trees in the ensemble, in order to tradeoff data fit with tree complexity. Use the Number of learners and Maximum number of splits settings. Advanced Ensemble Options Ensemble classifiers in Classification Learner use the fitcensemble function. You can set these options: ·
For help choosing Ensemble method and Learner type, see the Ensemble table. Try the presets first. · Maximum number of splits For boosting ensemble methods, specify the maximum number of splits or branch points to control the depth of your tree learners. Many branches tend to overfit, and simpler trees can be more robust and easy to interpret. Experiment to choose the best tree depth for the trees in the ensemble. Number of learners Try changing the number of learners to see if you can improve the model. Many learners can produce high accuracy, but can be time consuming to fit. Start with a few dozen learners, and then inspect the performance. An ensemble with good predictive power can need a few hundred learners. Learning rate Specify the learning rate for shrinkage. If you set the learning rate to less than 1, the ensemble requires more learning iterations but often achieves better accuracy. 0.1 is a popular choice Subspace dimension
For subspace ensembles, specify the number of predictors to sample in each learner. The app chooses a random subset of the predictors for each learner. The subsets chosen by different learners are independent.
2.7 ASSESS CLASSIFIER PERFORMANCE IN CLASSIFICATION LEARNER 2.7.1 Check Performance in the History List After training a model in Classification Learner, check the History list to see which model has the best overall accuracy in percent. The best Accuracy score is highlighted in a box. This score is the validation accuracy (unless you opted for no validation scheme). The validation accuracy score estimates a model's performance on new data compared to the training data. Use the score to help you choose the best model. ·
For cross-validation, the score is the accuracy on all observations, counting each observation when it was in a held-out fold. · For holdout validation, the score is the accuracy on the held-out observations. · For no validation, the score is the resubstitution accuracy against all the training data observations. The best overall score might not be the best model for your goal. A model with a slightly lower overall accuracy might be the best classifier for your goal. For example, false positives in a particular class might be important to you. You might want to exclude some predictors where data collection is expensive or difficult. To find out how the classifier performed in each class, examine the confusion matrix.
2.7.2 Plot Classifier Results In the scatter plot, view the classifier results. After you train a classifier, the scatter plot switches from displaying the data to showing model predictions. To investigate your results, use the controls on the right. You can: · ·
Choose whether to plot model predictions or the data alone Show or hide correct or incorrect results using the check boxes under Model predictions · Choose features to plot using the X and Y lists under Predictors · Visualize results by class by showing or hiding specific classes using the check boxes under Classes
2.7.3 Check Performance Per Class in the Confusion Matrix Use the confusion matrix plot to understand how the currently selected classifier performed in each class. To view the confusion matrix after training a model, on the Classification Learner tab, in the Plots section, click Confusion Matrix. The confusion matrix helps you identify the areas where the classifier has performed poorly. When you open the plot, the rows show the true class, and the columns show the predicted class. The diagonal cells show where the true class and predicted class match. If these cells are green, the classifier has performed well and classified observations of this true class correctly. The default view shows number of observations in each cell. To see how the classifier performed per class, under Plot, select the True Positive Rates, False Negative Rates option. The plot shows summaries per true class in the last two columns on the right.
In this example, using the carsmall data set, the top row shows all cars with true class France. The columns show the predicted classes. In the top row, 25% of the cars from France are correctly classified, so 25% is the true
positive rate for correctly classified points in this class, shown in the green cell in the True Positive Rate column. The other cars in the France row are misclassified: 50% of the cars are incorrectly classified as from Japan, and 25% are classified as from Sweden. 75% is the false negative rate for incorrectly classified points in this class, shown in the red cell in the False Negative Rate column. If you want to see numbers of observations (cars, in this example) instead of percentages, under Plot, select Number of observations. If false positives are important in your classification problem, plot results per predicted class (instead of true class) to investigate false discovery rates. To see results per predicted class, under Plot, select the Positive Predictive ValuesFalse Discovery Rates option. The confusion matrix now shows summary rows underneath the table. Positive predictive values are shown in green for the correctly predicted points in each class, and false discovery rates are shown below it in red for the incorrectly predicted points in each class. If you decide there are too many misclassified points in the classes of interest, try changing classifier settings or feature selection to search for a better model.
2.7.4 Check the ROC Curve To view the ROC curve after training a model, on the Classification Learner tab, in the Plots section, click ROC Curve. View the receiver operating characteristic (ROC) curve showing true and false positive rates. The ROC curve shows true positive rate versus false positive rate for the currently selected trained classifier. You can select different classes to plot. The marker on the plot shows the performance of the currently selected classifier. The marker shows the values of the false positive rate (FPR) and the true positive rate (TPR) for the currently selected classifier. For example, a false positive rate (FPR) of 0.2 indicates that the current classifier assigns 20% of the observations incorrectly to the positive class. A true positive rate of 0.9 indicates that the current classifier assigns 90% of the observations correctly to the positive class. A perfect result with no misclassified points is a right angle to the top left of the plot. A poor result that is no better than random is a line at 45 degrees. The Area Under Curve number is a measure of the overall quality of the classifier. Larger Area Under Curve values indicate better classifier performance. Compare classes and trained models to see if they perform differently in the ROC curve.
3Chapter 3
FIT MULTICLASS MODELS FOR SUPPORT VECTOR MACHINES OR OTHER CLASSIFIERS. CLASSIFICATION LEARNER APP
3.1 MULTICLASS ECOC FUNCTIONS fitcecoc templateSVM predict
Fit multiclass models for support vector machines or other classifiers Support vector machine template Predict labels using multiclass, error-correcting output codes model
3.1.1 fitcecoc Fit multiclass models for support vector machines or other classifiers Syntax Mdl = fitcecoc(Tbl,ResponseVarName) Mdl = fitcecoc(Tbl,formula) Mdl = fitcecoc(Tbl,Y) Mdl = fitcecoc(X,Y) Mdl = fitcecoc(___,Name,Value) [Mdl,HyperparameterOptimizationResults] = fitcecoc(___,Name,Value) Description Mdl = fitcecoc(Tbl,ResponseVarName) returns a full, trained, multiclass, error-correcting output codes (ECOC) model using the predictors in table Tbl and the class labels in Tbl.ResponseVarName. fitcecoc uses K(K – 1)/2 binary support vector machine (SVM) models using the one-versusone coding design, where K is the number of unique class labels (levels). Mdl is a ClassificationECOC model. Mdl = fitcecoc(Tbl,formula) returns an ECOC model using the predictors in table Tbl and the class labels. formula is an explanatory model of the response and a subset of predictor variables in Tbl used for training. Mdl = fitcecoc(Tbl,Y) returns an ECOC model using the predictors in table Tbl and the class labels in vector Y. Mdl = fitcecoc(X,Y) returns a full, trained ECOC model using the predictors X and the class labels Y. Mdl = fitcecoc(___,Name,Value) returns an ECOC model with additional options specified by one or more Name,Value pair arguments, using any of the previous syntaxes. For example, specify different binary learners, a different coding design, or to cross-validate. It is good practice to cross-validate using the Kfold Name,Value pair argument. The cross-validation results determine
how well the model generalizes. [Mdl,HyperparameterOptimizationResults] = fitcecoc(___,Name,Value) also returns hyperparameter optimization details when you pass an OptimizeHyperparameters name-value pair with Learners = 'linear'. For other Learners, the HyperparameterOptimizationResults property of Mdlcontains the results. Examples: Train Multiclass Model Using SVM Learners Train an error-correcting output codes (ECOC) multiclass model using support vector machine (SVM) binary learners. Load Fisher's iris data set. load fisheriris X = meas; Y = species; Train an ECOC multiclass model using the default options. Mdl = fitcecoc(X,Y) Mdl = ClassificationECOC ResponseName: 'Y' CategoricalPredictors: [] ClassNames: {'setosa' 'versicolor' 'virginica'} ScoreTransform: 'none' BinaryLearners: {3×1 cell} CodingName: 'onevsone' Mdl is a ClassificationECOC model. By default, fitcecoc uses SVM binary learners, and uses a one-versus-one coding design. You can access Mdl properties using dot notation. Display the coding design matrix.
Mdl.ClassNames CodingMat = Mdl.CodingMatrix ans = 3×1 cell array 'setosa' 'versicolor' 'virginica' CodingMat = 1 1 -1 0 0 -1
0 1 -1
A one-versus-one coding design on three classes yields three binary learners. Columns of CodingMat correspond to learners, and rows correspond to classes. The class order corresponds to the order in Mdl.ClassNames. For example, CodingMat(:,1) is [1; -1; 0] and indicates that the software trains the first SVM binary learner using all observations classified as 'setosa' and 'versicolor'. Since 'setosa' corresponds to 1, it is the positive class, and since 'versicolor' corresponds to -1, it is the negative class. You can access each binary learner using cell indexing and dot notation. Mdl.BinaryLearners{1}
% The first binary learner
ans = classreg.learning.classif.CompactClassificationSVM ResponseName: 'Y' CategoricalPredictors: [] ClassNames: [-1 1] ScoreTransform: 'none' Beta: [4×1 double] Bias: 1.4505
KernelParameters: [1×1 struct] Compute the in-sample classification error. isLoss = resubLoss(Mdl) isLoss = 0.0067 The classification error is small, but the classifier might have been overfit. You can cross-validate the classifier using crossval. Cross-Validate ECOC Classifier Load Fisher's iris data set. load fisheriris X = meas; Y = species; rng(1); % For reproducibility Create an SVM template. It is good practice to standardize the predictors. t = templateSVM('Standardize',1) t= Fit template for classification SVM. Alpha: [0×1 double] BoxConstraint: [] CacheSize: [] CachingMethod: '' ClipAlphas: [] DeltaGradientTolerance: [] Epsilon: [] GapTolerance: [] KKTTolerance: [] IterationLimit: []
KernelFunction: '' KernelScale: [] KernelOffset: [] KernelPolynomialOrder: [] NumPrint: [] Nu: [] OutlierFraction: [] RemoveDuplicates: [] ShrinkagePeriod: [] Solver: '' StandardizeData: 1 SaveSupportVectors: [] VerbosityLevel: [] Version: 2 Method: 'SVM' Type: 'classification' t is an SVM template. All of its properties are empty, except for StandardizeData, Method, and Type. When the software trains the ECOC classifier, it sets the applicable properties to their default values. Train the ECOC classifier. It is good practice to specify the class order. Mdl = fitcecoc(X,Y,'Learners',t,... 'ClassNames',{'setosa','versicolor','virginica'}); Mdl is a ClassificationECOC classifier. You can access its properties using dot notation. Cross-validate Mdl using 10-fold cross-validation. CVMdl = crossval(Mdl); CVMdl is a ClassificationPartitionedECOC cross-validated ECOC classifier. Estimate the generalization error. oosLoss = kfoldLoss(CVMdl) oosLoss = 0.0400
The out-of-sample classification error is 4%, which indicates that the ECOC classifier generalizes fairly well. Estimate Posterior Probabilities Using ECOC Classifiers Load Fisher's iris data set. Train the classifier using the petal dimensions as predictors. load fisheriris X = meas(:,3:4); Y = species; rng(1); % For reproducibility Create an SVM template, and specify the Gaussian kernel. It is good practice to standardize the predictors. t = templateSVM('Standardize',1,'KernelFunction','gaussian'); t is an SVM template. Most of its properties are empty. When the software trains the ECOC classifier, it sets the applicable properties to their default values. Train the ECOC classifier using the SVM template. Transform classification scores to class posterior probabilities (which are returned by predict or resubPredict) using the 'FitPosterior' name-value pair argument. Display diagnostic messages during the training using the 'Verbose' name-value pair argument. It is good practice to specify the class order. Mdl = fitcecoc(X,Y,'Learners',t,'FitPosterior',1,... 'ClassNames',{'setosa','versicolor','virginica'},... 'Verbose',2); Training binary learner 1 (SVM) out of 3 with 50 negative and 50 positive observations. Negative class indices: 2 Positive class indices: 1 Fitting posterior probabilities for learner 1 (SVM). Training binary learner 2 (SVM) out of 3 with 50 negative and 50 positive observations.
Negative class indices: 3 Positive class indices: 1 Fitting posterior probabilities for learner 2 (SVM). Training binary learner 3 (SVM) out of 3 with 50 negative and 50 positive observations. Negative class indices: 3 Positive class indices: 2 Fitting posterior probabilities for learner 3 (SVM). Mdl is a ClassificationECOC model. The same SVM template applies to each binary learner, but you can adjust options for each binary learner by passing in a cell vector of templates. Predict the in-sample labels and class posterior probabilities. Display diagnostic messages during the computation of labels and class posterior probabilities using the 'Verbose' name-value pair argument. [label,~,~,Posterior] = resubPredict(Mdl,'Verbose',1); Mdl.BinaryLoss Predictions from all learners have been computed. Loss for all observations has been computed. Computing posterior probabilities... ans = 'quadratic' The software assigns an observation to the class that yields the smallest average binary loss. Since all binary learners are computing posterior probabilities, the binary loss function is quadratic. Display a random set of results. idx = randsample(size(X,1),10,1); Mdl.ClassNames table(Y(idx),label(idx),Posterior(idx,:),... 'VariableNames',{'TrueLabel','PredLabel','Posterior'})
ans = 3×1 cell array 'setosa' 'versicolor' 'virginica' ans = 10×3 table TrueLabel PredLabel Posterior ____________ ____________ ______________________________________ 'virginica' 'virginica' 'virginica' 'virginica' 'virginica' 'virginica' 'versicolor' 'versicolor' 'setosa' 'setosa' 'versicolor' 'virginica' 'versicolor' 'versicolor' 'setosa' 'setosa' 'versicolor' 'versicolor' 'setosa' 'setosa'
0.0039321 0.0039869 0.99208 0.017067 0.018263 0.96467 0.014948 0.015856 0.9692 2.2197e-14 0.87317 0.12683 0.999 0.00025091 0.00074639 2.2195e-14 0.059429 0.94057 2.2194e-14 0.97001 0.029986 0.999 0.0002499 0.00074741 0.0085646 0.98259 0.008849 0.999 0.00025013 0.00074718
The columns of Posterior correspond to the class order of Mdl.ClassNames. Define a grid of values in the observed predictor space. Predict the posterior probabilities for each instance in the grid. xMax = max(X); xMin = min(X); x1Pts = linspace(xMin(1),xMax(1)); x2Pts = linspace(xMin(2),xMax(2)); [x1Grid,x2Grid] = meshgrid(x1Pts,x2Pts); [~,~,~,PosteriorRegion] = predict(Mdl,[x1Grid(:),x2Grid(:)]);
For each coordinate on the grid, plot the maximum class posterior probability among all classes. figure; contourf(x1Grid,x2Grid,... reshape(max(PosteriorRegion,[],2),size(x1Grid,1),size(x1Grid,2))); h = colorbar; h.YLabel.String = 'Maximum posterior'; h.YLabel.FontSize = 15; hold on gh = gscatter(X(:,1),X(:,2),Y,'krk','*xd',8); gh(2).LineWidth = 2; gh(3).LineWidth = 2; title 'Iris Petal Measurements and Maximum Posterior'; xlabel 'Petal length (cm)'; ylabel 'Petal width (cm)'; axis tight legend(gh,'Location','NorthWest') hold off
Train Multiclass Linear Classification Model Train an ECOC model composed of multiple binary, linear classification models. Load the NLP data set. load nlpdata X is a sparse matrix of predictor data, and Y is a categorical vector of class labels. There are more than two classes in the data. Create a default linear-classification-model template. t = templateLinear(); To adjust the default values, see templateLinear. Train an ECOC model composed of multiple binary, linear classification models that can identify the product given the frequency distribution of words on a documentation web page. For faster training time, transpose the predictor data, and specify that observations correspond to columns. X = X'; rng(1); % For reproducibility Mdl = fitcecoc(X,Y,'Learners',t,'ObservationsIn','columns') Mdl = classreg.learning.classif.CompactClassificationECOC ResponseName: 'Y' ClassNames: [1×13 categorical] ScoreTransform: 'none' BinaryLearners: {78×1 cell} CodingMatrix: [13×78 double] Alternatively, you can train an ECOC model composed of default linear classification models using 'Learners','Linear'. To conserve memory, fitcecoc returns trained ECOC models composed of linear classification learners in CompactClassificationECOC model objects.
Train ECOC Classifiers Using Ensembles and Parallel Computing Train a one-versus-all ECOC classifier using a GentleBoost ensemble of decision trees with surrogate splits. Estimate the classification error using 10fold cross-validation. Load and inspect the arrhythmia data set. load arrhythmia [n,p] = size(X) isLabels = unique(Y); nLabels = numel(isLabels) tabulate(categorical(Y)) n= 452 p= 279 nLabels = 13 Value 1 2 3 4 5 6 7 8 9 10 14
Count Percent 245 54.20% 44 9.73% 15 3.32% 15 3.32% 13 2.88% 25 5.53% 3 0.66% 2 0.44% 9 1.99% 50 11.06% 4 0.88%
15 16
5 22
1.11% 4.87%
There are 279 predictors, and a relatively small sample size of 452. There are 16 distinct labels, but only 13 are represented in the response (Y), and each label describes various degrees of arrhythmia. 54.20% of the observations are in class 1. Create an ensemble template. You must specify at least three arguments: a method, a number of learners, and the type of learner. For this example, specify 'GentleBoost'for the method, 100 for the number of learners, and a decision tree template that uses surrogate splits since there are missing observations. tTree = templateTree('surrogate','on'); tEnsemble = templateEnsemble('GentleBoost',100,tTree); tEnsemble is a template object. Most of its properties are empty, but the software fills them with their default values during training. Train a one-versus-all ECOC classifier using the ensembles of decision trees as binary learners. If you have a Parallel Computing Toolbox™ license, then you can speed up the computation by specifying to use parallel computing. This sends each binary learner to a worker in the pool (the number of workers depends on your system configuration). Also, specify that the prior probabilities are 1/K, where K = 13, which is the number of distinct classes. pool = parpool; % Invoke workers options = statset('UseParallel',1); Mdl = fitcecoc(X,Y,'Coding','onevsall','Learners',tEnsemble,... 'Prior','uniform','Options',options); Starting parallel pool (parpool) using the 'local' profile ... connected to 2 workers. Mdl is a ClassificationECOC model. Cross-validate the ECOC classifier using 10-fold cross-validation. CVMdl = crossval(Mdl,'Options',options); Warning: One or more folds do not contain points from all the groups. CVMdl is a ClassificationPartitionedECOC model. The warning indicates that some classes are not represented while the software trains at least one fold.
Therefore, those folds cannot predict labels for the missing classes. You can inspect the results of a fold using cell indexing and dot notation, e.g., access the results of the first fold by entering CVMdl.Trained{1}. Your results might vary. Use the cross-validated ECOC classifier to predict out-of-fold labels. You can compute the confusion matrix using confusionmat. However, if you have a Neural Network Toolbox™ license, you can plot the confusion matrix using plotconfusion. The input arguments of plotconfusion are not vectors of the true and predicted labels like confusionmat, but indicator matrices of the true and predicted labels. Both start as K-by-n matrices of 0s. If observation j has label index k (or has predicted label k), then element (k,j) of the true label indicator matrix (or predicted label indicator matrix) is 1. You can convert label indices returned by predict, resubPredict, or kfoldPredictto label indicator matrices using linear indexing. For details on linear indexing, see sub2ind and ind2sub. oofLabel = kfoldPredict(CVMdl,'Options',options); ConfMat = confusionmat(Y,oofLabel); % Convert the integer label vector to a class-identifier matrix. [~,grpOOF] = ismember(oofLabel,isLabels); oofLabelMat = zeros(nLabels,n); idxLinear = sub2ind([nLabels n],grpOOF,(1:n)'); oofLabelMat(idxLinear) = 1; % Flags the row corresponding to the class [~,grpY] = ismember(Y,isLabels); YMat = zeros(nLabels,n); idxLinearY = sub2ind([nLabels n],grpY,(1:n)'); YMat(idxLinearY) = 1; figure; plotconfusion(YMat,oofLabelMat); h = gca; h.XTickLabel = [num2cell(isLabels); {''}]; h.YTickLabel = [num2cell(isLabels); {''}];
Optimize ECOC Classifier This example shows how to optimize hyperparameters automatically using fitcecoc. The example uses Fisher's iris data. Load the data. load fisheriris X = meas; Y = species; Find hyperparameters that minimize five-fold cross-validation loss by using automatic hyperparameter optimization. For reproducibility, set the random seed and use the 'expected-improvementplus' acquisition function. rng default Mdl = fitcecoc(X,Y,'OptimizeHyperparameters','auto',... 'HyperparameterOptimizationOptions',struct('AcquisitionFunctionName',... 'expected-improvement-plus'))
|======================================================================================================= | Iter | Eval | Objective | Objective | BestSoFar | BestSoFar | Coding | BoxConstrain-| KernelScale | | | result | | runtime | (observed) | (estim.) | |t | | |======================================================================================================= | 1 | Best | 0.33333 | 3.3156 | 0.33333 | 0.33333 | onevsall | 0.0089562 | 310.42 | | 2 | Best | 0.053333 | 5.2976 | 0.053333 | 0.073138 | onevsone | 67.106 | 0.0053534 | | 3 | Best | 0.04 | 0.57165 | 0.04 | 0.040121 | onevsall | 0.0027084 | 0.017254 | | 4 | Accept | 0.06 | 13.097 | 0.04 | 0.040005 | onevsall | 259.06 | 0.14755 | | 5 | Accept | 0.046667 | 0.66074 | 0.04 | 0.039999 | onevsone | 127.23 | 1.0303 | | 6 | Accept | 0.11333 | 0.41791 | 0.04 | 0.04 | onevsone | 0.37412 | 996.84 | | 7 | Accept | 0.48 | 20.956 | 0.04 | 0.040038 | onevsall | 182.3 | 0.0010002 | | 8 | Accept | 0.04 | 0.76888 | 0.04 | 0.039982 | onevsall | 0.090529 | 0.042615 | | 9 | Best | 0.02 | 0.54564 | 0.02 | 0.020009 | onevsone | 0.0030779 | 0.09393 | | 10 | Accept | 0.046667 | 0.41283 | 0.02 | 0.020015 | onevsone | 1.9001 | 0.02549 | | 11 | Accept | 0.11333 | 0.45093 | 0.02 | 0.020016 | onevsone | 0.082855 | 24.057 | | 12 | Accept | 0.02 | 0.55456 | 0.02 | 0.019995 | onevsone | 0.028965 | 0.26902 | | 13 | Accept | 0.046667 | 5.7777 | 0.02 | 0.020004 | onevsone | 74.47 | 0.0010042 | | 14 | Accept | 0.33333 | 0.37761 | 0.02 | 0.020001 | onevsall | 0.0010257 | 0.2665 | | 15 | Accept | 0.26 | 19.695 | 0.02 | 0.020003 | onevsall | 971.54 | 0.04193 | | 16 | Accept | 0.04 | 0.82597 | 0.02 | 0.020003 | onevsall | 83.197 | 0.82502 | | 17 | Accept | 0.046667 | 0.5434 | 0.02 | 0.020004 | onevsone | 0.028079 | 0.0024331 | | 18 | Accept | 0.04 | 5.9089 | 0.02 | 0.020207 | onevsall | 0.0011615 | 0.0019763 | | 19 | Accept | 0.04 | 1.8263 | 0.02 | 0.020004 | onevsall | 0.020594 | 0.0073833 | | 20 | Accept | 0.04 | 1.3375 | 0.02 | 0.020004 | onevsall | 9.6782 | 0.19642 | |======================================================================================================= | Iter | Eval | Objective | Objective | BestSoFar | BestSoFar | Coding | BoxConstrain-| KernelScale | | | result | | runtime | (observed) | (estim.) | |t | | |======================================================================================================= | 21 | Accept | 0.02 | 0.3982 | 0.02 | 0.020005 | onevsone | 965.24 | 39.212 | | 22 | Accept | 0.12 | 0.47837 | 0.02 | 0.020005 | onevsone | 980.8 | 643.21 | | 23 | Accept | 0.04 | 0.41142 | 0.02 | 0.020037 | onevsone | 973.21 | 8.1701 | | 24 | Accept | 0.11333 | 0.4301 | 0.02 | 0.020026 | onevsone | 0.0010507 | 0.69618 | | 25 | Accept | 0.033333 | 0.46536 | 0.02 | 0.020011 | onevsone | 0.0011494 | 0.012414 | | 26 | Accept | 0.033333 | 0.47349 | 0.02 | 0.020506 | onevsone | 0.035003 | 0.058611 | | 27 | Accept | 0.033333 | 0.3971 | 0.02 | 0.020059 | onevsone | 1.2763 | 0.34837 | | 28 | Accept | 0.046667 | 2.4195 | 0.02 | 0.020367 | onevsone | 978.94 | 0.096199 | | 29 | Accept | 0.046667 | 0.36933 | 0.02 | 0.020471 | onevsone | 52.349 | 34.743 | | 30 | Accept | 0.046667 | 7.0458 | 0.02 | 0.022071 | onevsall | 984.41 | 0.65742 |
__________________________________________________________ Optimization completed. MaxObjectiveEvaluations of 30 reached. Total function evaluations: 30 Total elapsed time: 160.6519 seconds. Total objective function evaluation time: 96.2299 Best observed feasible point: Coding BoxConstraint KernelScale ________ _____________ ___________
onevsone
0.0030779
0.09393
Observed objective function value = 0.02 Estimated objective function value = 0.022071 Function evaluation time = 0.54564 Best estimated feasible point (according to models): Coding BoxConstraint KernelScale ________ _____________ ___________ onevsone
0.0030779
0.09393
Estimated objective function value = 0.022071 Estimated function evaluation time = 0.49997 Mdl = ClassificationECOC ResponseName: 'Y' CategoricalPredictors: [] ClassNames: {'setosa' 'versicolor' 'virginica'} ScoreTransform: 'none' BinaryLearners: {3×1 cell} CodingName: 'onevsone' HyperparameterOptimizationResults: [1×1 BayesianOptimization]
3.1.2 CompactClassificationECOC class Compact multiclass model for support vector machines or other classifiers
Description CompactClassificationECOC is a compact, multiclass, error-correcting output codes (ECOC) model. The compact classifier does not include the data used for training the ECOC multiclass model. Therefore, you cannot perform tasks, such as crossvalidation, using the compact classifier. Use a compact ECOC multiclass model for labeling new data (in other words, predicting the labels of new data).
Construction CompactMdl = compact(Mdl) returns a compact, multiclass ECOC model (CompactModel) from a full, trained, multiclass, ECOC model (Mdl). CompactMdl = fitcecoc(X,Y,'Learners',t) returns a compact ECOC model composed of linear classification models (ClassificationLinear model objects). t must be the character vector 'linear' or a linear-classification-model template object returned by templateLinear. For more details, see fitcecoc
Examples: Reduce the Size of Full ECOC Models Full ECOC models (i.e., ClassificationECOC classifiers) hold the training data. For efficiency, you might not want to predict new labels using a large classifier. Load Fisher's iris data set. load fisheriris X = meas;
Y = species; Train an ECOC model using default SVMs as binary learners. Mdl = fitcecoc(X,Y) Mdl = ClassificationECOC ResponseName: 'Y' CategoricalPredictors: [] ClassNames: {'setosa' 'versicolor' 'virginica'} ScoreTransform: 'none' BinaryLearners: {3×1 cell} CodingName: 'onevsone' Mdl is a ClassificationECOC model. Reduce the size of the ECOC model. CMdl = compact(Mdl) CMdl = classreg.learning.classif.CompactClassificationECOC ResponseName: 'Y' CategoricalPredictors: [] ClassNames: {'setosa' 'versicolor' 'virginica'} ScoreTransform: 'none' BinaryLearners: {3×1 cell} CodingMatrix: [3×3 double] CMdl is a CompactClassificationECOC model. Display how much memory each classifier consumes. whos('Mdl','CMdl') Name
Size
Bytes Class
CMdl 1x1 12764 classreg.learning.classif.CompactClassificationECOC
Attributes
Mdl
1x1
25337 ClassificationECOC
The full ECOC model (Mdl) is almost twice the size of the compact ECOC model (CMdl). You can remove Mdl from the MATLAB® Workspace, and pass CMdl and new predictor values to predict to efficiently label new observations. Train and Cross-Validate ECOC Classifiers Train an ECOC classifier using different binary learners and the one-versus-all coding design. Then, cross validate the classifier. Load Fisher's iris data set. load fisheriris X = meas; Y = species; classNames = unique(species(~strcmp(species,''))) % Remove empty classes K = numel(classNames) % Number of classes rng(1); % For reproducibility classNames = 3×1 cell array 'setosa' 'versicolor' 'virginica' K= 3 classNames are the unique classes in the data set, and K is the number of classes. You can use classNames to specify the order of the classes during training.
For a one-versus-all coding design, there are K = 3 binary learners. Specify templates for the binary learners such that: Binary learner 1 and 2 are naive Bayes classifiers. By default, each predictor is conditionally, normally distributed given its label. Binary learner 3 is an SVM classifier. Specify to use the Gaussian kernel. tNB = templateNaiveBayes(); tSVM = templateSVM('KernelFunction','gaussian'); tLearners = {tNB tNB tSVM}; tNB and tSVM are template objects for naive Bayes and SVM learning, respectively. They indicate what options to use during training. Most of their properties are empty, except for those specified using name-value pair arguments. The software fills in the empty properties with their default values during training. Train and cross-validate an ECOC classifier using the binary learner templates and the one-versus-all coding design. Specify the order of the classes. By default, naive Bayes classifiers use posterior probabilities as scores, whereas SVM classifiers use distance from the decision boundary. Therefore, to aggregate the binary learners, you must specify fitting posterior probabilities. CVMdl = fitcecoc(X,Y,'ClassNames',classNames,'CrossVal','on',... 'Learners',tLearners,'FitPosterior',1); CVMdl is not a ClassificationECOC model, but a ClassificationPartitionedECOC cross-validated model. By default, the software implements 10-fold cross-validation. The scores across the binary learners are in the same form (that is, they are posterior probabilities), so the software can aggregate the results of the binary classifications properly. Inspect one of the trained folds using dot notation. CVMdl.Trained{1} ans = classreg.learning.classif.CompactClassificationECOC ResponseName: 'Y'
CategoricalPredictors: [] ClassNames: {'setosa' 'versicolor' 'virginica'} ScoreTransform: 'none' BinaryLearners: {3×1 cell} CodingMatrix: [3×3 double] Each fold is a CompactClassificationECOC model trained on 90% of the data. You can access the results of the binary learners using dot notation and cell indexing. Display the trained SVM classifier (the third binary learner) in the first fold. CVMdl.Trained{1}.BinaryLearners{3} ans = classreg.learning.classif.CompactClassificationSVM ResponseName: 'Y' CategoricalPredictors: [] ClassNames: [-1 1] ScoreTransform: '@(S)sigmoid(S,-4.016268e+00,-3.243751e-01)' Alpha: [33×1 double] Bias: -0.1345 KernelParameters: [1×1 struct] SupportVectors: [33×4 double] SupportVectorLabels: [33×1 double] Estimate the generalization error. genError = kfoldLoss(CVMdl) genError = 0.0333 On average, the generalization error is approximately 3%.
3.2 TRAIN ENSEMBLE CLASSIFIERS USING CLASSIFICATION LEARNER APP This example shows how to construct ensembles of classifiers in the Classification Learner app. Ensemble classifiers meld results from many weak learners into one high-quality ensemble predictor. Qualities depend on the choice of algorithm, but ensemble classifiers tend to be slow to fit because they often need many learners. 1.
In MATLAB®, load the fisheriris data set and define some variables from the data set to use for a classification. fishertable = readtable('fisheriris.csv');
2.
On the Apps tab, in the Math, Statistics and Optimization group, click Classification Learner. On the Classification Learner tab, in the File section, click New Session. In the New Session dialog box, observe that the app has selected response and predictor variables based on their data type. Petal and sepal length and width are predictors. Species is the response that you want to classify. For this example, do not change the selections. Click Start Session. Classification Learner creates a scatter plot of the data. Use the scatter plot to investigate which variables are useful for predicting the response. Select different variables in the X- and Y-axis controls to visualize the distribution of species and measurements. Observe which variables separate the species colors most clearly. To create a selection of ensemble models, on the Classification Learner tab, in the Model Type section, click the down arrow to expand the list of classifiers, then under Ensemble Classifiers, click All Ensembles.
3.
4. 5.
6.
7.
In the Training section, click Train.
Tip: If you have Parallel Computing Toolbox™ then the first time you click Train you see a dialog while the app opens a parallel pool of workers. After the pool opens, you can train multiple classifiers at once and continue working.
8.
Classification Learner trains one of each ensemble classification option in the gallery, and highlights the best score. The app outlines in a box the Accuracy score of the best model. 9. Select a model in the History list to view the results. Examine the scatter plot for the trained model. Misclassified points are shown as an X.
10. To inspect the accuracy of the predictions in each class, on the Classification Learner tab, in the Plots section, click Confusion Matrix. View the matrix of true class and predicted class results. 11. Select the other models in the list to compare. 12. Choose the best model in the History list (the best score is highlighted in the Accuracy box). To improve the model, try including different features in the model. See if you can improve the model by removing features with low
predictive power. On the Classification Learner tab, in the Data Transformation section, click Features. In the Feature Selection dialog box, specify predictors to remove from the model, and click Train to train a new model using the new options. Compare results among the classifiers in the History list. 13. To investigate features to include or exclude, use the scatter and parallel coordinates plots. On the Classification Learner tab, in the Plots section, select Parallel Coordinates Plot. 14. Choose the best model in the History list. To try to improve the model further, try changing settings. On the Classification Learner tab, in the Model Type section, click Advanced. Try changing a setting, then train the new model by clicking Train. 15. To export the trained model to the workspace, select the Classification Learner tab and click Export model. 16. To examine the code for training this classifier, select Export Model > Generate Code. Use the same workflow to evaluate and compare the other classifier types you can train in Classification Learner. To try all the classifier model presets available for your data set: 1.
Click the arrow on the far right of the Model Type section to expand the list of classifiers. 2. Click All, then click Train.
4Chapter 4
PREDICTIVE ANALYTICS. SUPPORT VECTOR MACHINE REGRESSION (SVR) AND GAUSSIAN PROCCES REGRESSION (GPR)
4.1 SUPPORT VECTOR MACHINE REGRESSION MATLAB use functions for support vector machine regression. For greater accuracy on low- through medium-dimensional data sets, train a support vector machine (SVM) model using fitrsvm. For reduced computation time on high-dimensional data sets that fit in the MATLAB® Workspace, efficiently train a linear regression model, such as a linear SVM model, using fitrlinear.
4.1.1 fitrsvm Fit a support vector machine regression model fitrsvm trains or cross-validates a support vector machine (SVM) regression model on a low- through moderate-dimensional predictor data set. fitrsvm supports mapping the predictor data using kernel functions, and supports SMO, ISDA, or L1 soft-margin minimization via quadratic programming for objective-function minimization. To train a linear SVM regression model on a high-dimensional data set, that is, data sets that include many predictor variables, use fitrlinear instead. To train an SVM model for binary classification, see fitcsvm for low- through moderate-dimensional predictor data sets, or fitclinear for high-dimensional data sets.
Syntax Mdl = fitrsvm(Tbl,ResponseVarName) Mdl = fitrsvm(Tbl,formula) Mdl = fitrsvm(Tbl,Y) Mdl = fitrsvm(X,Y) Mdl = fitrsvm(___,Name,Value)
Description Mdl = fitrsvm(Tbl,ResponseVarName) returns a full, trained support vector machine (SVM) regression model Mdl trained using the predictors values in the table Tbl and the response values in Tbl.ResponseVarName. Mdl = fitrsvm(Tbl,formula) returns a full SVM regression model trained using the predictors values in the table Tbl. formula is an explanatory model of the response and a subset of predictor variables in Tblused to fit Mdl.
Mdl = fitrsvm(Tbl,Y) returns a full, trained SVM regression model trained using the predictors values in the table Tbl and the response values in the vector Y. Mdl = fitrsvm(X,Y) returns a full, trained SVM regression model trained using the predictors values in the matrix X and the response values in the vector Y. Mdl = fitrsvm(___,Name,Value) returns an SVM regression model with additional options specified by one or more name-value pair arguments, using any of the previous syntaxes. For example, you can specify the kernel function or train a cross-validated model.
Examples: Train Linear Support Vector Machine Regression Model Train a support vector machine (SVM) regression model using sample data stored in matrices. Load the carsmall data set. load carsmall rng 'default' % For reproducibility Specify Horsepower and Weight as the predictor variables (X) and MPG as the response variable (Y). X = [Horsepower,Weight]; Y = MPG; Train a default SVM regression model. Mdl = fitrsvm(X,Y) Mdl = RegressionSVM ResponseName: 'Y'
CategoricalPredictors: [] ResponseTransform: 'none' Alpha: [76×1 double] Bias: 43.1920 KernelParameters: [1×1 struct] NumObservations: 93 BoxConstraints: [93×1 double] ConvergenceInfo: [1×1 struct] IsSupportVector: [93×1 logical] Solver: 'SMO' Mdl is a trained RegressionSVM model. Check the model for convergence. Mdl.ConvergenceInfo.Converged ans = logical 0 0 indicates that the model did not converge. Retrain the model using standardized data. MdlStd = fitrsvm(X,Y,'Standardize',true) MdlStd = RegressionSVM ResponseName: 'Y' CategoricalPredictors: [] ResponseTransform: 'none' Alpha: [77×1 double] Bias: 22.9131 KernelParameters: [1×1 struct] Mu: [109.3441 2.9625e+03]
Sigma: [45.3545 805.9668] NumObservations: 93 BoxConstraints: [93×1 double] ConvergenceInfo: [1×1 struct] IsSupportVector: [93×1 logical] Solver: 'SMO' Check the model for convergence. MdlStd.ConvergenceInfo.Converged ans = logical 1 1 indicates that the model did converge. Compute the resubstitution (in-sample) mean-squared error for the new model. lStd = resubLoss(MdlStd) lStd = 17.0256
Train Support Vector Machine Regression Model Train a support vector machine regression model using the abalone data from the UCI Machine Learning Repository. Download the data and save it in your current folder with the name 'abalone.csv'. url = 'https://archive.ics.uci.edu/ml/machine-learningdatabases/abalone/abalone.data'; websave('abalone.csv',url);
Read the data into a table. Specify the variable names. varnames = {'Sex'; 'Length'; 'Diameter'; 'Height'; 'Whole_weight';... 'Shucked_weight'; 'Viscera_weight'; 'Shell_weight'; 'Rings'}; Tbl = readtable('abalone.csv','Filetype','text','ReadVariableNames',false); Tbl.Properties.VariableNames = varnames; The sample data contains 4177 observations. All the predictor variables are continuous except for Sex, which is a categorical variable with possible values 'M' (for males), 'F' (for females), and 'I' (for infants). The goal is to predict the number of rings (stored in Rings) on the abalone and determine its age using physical measurements. Train an SVM regression model, using a Gaussian kernel function with an automatic kernel scale. Standardize the data. rng default % For reproducibility Mdl = fitrsvm(Tbl,'Rings','KernelFunction','gaussian','KernelScale','auto',... 'Standardize',true) Mdl = RegressionSVM PredictorNames: {1×8 cell} ResponseName: 'Rings' CategoricalPredictors: 1 ResponseTransform: 'none' Alpha: [3635×1 double] Bias: 10.8144 KernelParameters: [1×1 struct] Mu: [1×10 double] Sigma: [1×10 double] NumObservations: 4177 BoxConstraints: [4177×1 double] ConvergenceInfo: [1×1 struct] IsSupportVector: [4177×1 logical] Solver: 'SMO'
The Command Window shows that Mdl is a trained RegressionSVM model and displays a property list. Display the properties of Mdl using dot notation. For example, check to confirm whether the model converged and how many iterations it completed. conv = Mdl.ConvergenceInfo.Converged iter = Mdl.NumIterations conv = logical 1 iter = 2759 The returned results indicate that the model converged after 2759 iterations.
Cross-Validate SVM Regression Model Load the carsmall data set. load carsmall rng 'default' % For reproducibility Specify Horsepower and Weight as the predictor variables (X) and MPG as the response variable (Y). X = [Horsepower Weight]; Y = MPG; Cross-validate two SVM regression models using 5-fold cross-validation. For both models, specify to standardize the predictors. For one of the models, specify to train using the default linear kernel, and the Gaussian kernel for the other model. MdlLin = fitrsvm(X,Y,'Standardize',true,'KFold',5) MdlGau = fitrsvm(X,Y,'Standardize',true,'KFold',5,'KernelFunction','gaussian')
MdlLin.Trained MdlLin = classreg.learning.partition.RegressionPartitionedSVM CrossValidatedModel: 'SVM' PredictorNames: {'x1' 'x2'} ResponseName: 'Y' NumObservations: 94 KFold: 5 Partition: [1×1 cvpartition] ResponseTransform: 'none' MdlGau = classreg.learning.partition.RegressionPartitionedSVM CrossValidatedModel: 'SVM' PredictorNames: {'x1' 'x2'} ResponseName: 'Y' NumObservations: 94 KFold: 5 Partition: [1×1 cvpartition] ResponseTransform: 'none' ans = 5×1 cell array [1×1 classreg.learning.regr.CompactRegressionSVM] [1×1 classreg.learning.regr.CompactRegressionSVM] [1×1 classreg.learning.regr.CompactRegressionSVM] [1×1 classreg.learning.regr.CompactRegressionSVM] [1×1 classreg.learning.regr.CompactRegressionSVM] MdlLin and MdlGau are RegressionPartitionedSVM cross-validated models.
The Trained property of each model is a 5-by-1 cell array of CompactRegressionSVM models. The models in the cell store the results of training on 4 folds of observations, and leaving one fold of observations out. Compare the generalization error of the models. In this case, the generalization error is the out-of-sample mean-squared error. mseLin = kfoldLoss(MdlLin) mseGau = kfoldLoss(MdlGau) mseLin = 17.4417 mseGau = 16.7397 The SVM regression model using the Gaussian kernel performs better than the one using the linear kernel. Create a model suitable for making predictions by passing the entire data set to fitrsvm, and specify all name-value pair arguments that yielded the betterperforming model. However, do not specify any cross-validation options. MdlGau = fitrsvm(X,Y,'Standardize',true,'KernelFunction','gaussian'); To predict the MPG of a set of cars, pass Mdl and a table containing the horsepower and weight measurements of the cars to predict.
Optimize SVM Regression This example shows how to optimize hyperparameters automatically using fitrsvm. The example uses the carsmall data. Load the carsmall data set.
load carsmall Specify Horsepower and Weight as the predictor variables (X) and MPG as the response variable (Y). X = [Horsepower Weight]; Y = MPG; Find hyperparameters that minimize five-fold cross-validation loss by using automatic hyperparameter optimization. For reproducibility, set the random seed and use the 'expected-improvementplus' acquisition function. rng default Mdl = fitrsvm(X,Y,'OptimizeHyperparameters','auto',... 'HyperparameterOptimizationOptions',struct('AcquisitionFunctionName',... 'expected-improvement-plus'))
|======================================================================================================= | Iter | Eval | Objective | Objective | BestSoFar | BestSoFar | BoxConstrain-| KernelScale | Epsilon | | | result | | runtime | (observed) | (estim.) | t | | | |======================================================================================================= | 1 | Best | 4.1768 | 1.9579 | 4.1768 | 4.1768 | 0.0023098 | 0.0089562 | 349.56 | | 2 | Best | 2.9448 | 1.7426 | 2.9448 | 3.0114 | 168.32 | 67.106 | 0.037505 | | 3 | Accept | 14.892 | 11.94 | 2.9448 | 3.4392 | 0.0082243 | 0.0027084 | 0.099459 | | 4 | Accept | 3.2603 | 0.11411 | 2.9448 | 2.9452 | 0.011994 | 259.06 | 0.59478 | | 5 | Best | 2.9188 | 0.33923 | 2.9188 | 2.9189 | 49.411 | 994.96 | 0.015002 | | 6 | Accept | 4.1768 | 0.076366 | 2.9188 | 2.9188 | 0.10907 | 4.5556 | 926.55 | | 7 | Accept | 4.1768 | 0.069372 | 2.9188 | 2.9187 | 15.476 | 998.99 | 186.25 | | 8 | Accept | 4.146 | 0.072151 | 2.9188 | 2.9187 | 0.015651 | 13.914 | 15.354 | | 9 | Accept | 3.8433 | 0.067662 | 2.9188 | 2.9191 | 0.0019192 | 210.7 | 0.0093026 | | 10 | Accept | 2.9398 | 0.13751 | 2.9188 | 2.9209 | 997.32 | 838.82 | 0.20365 | | 11 | Accept | 4.1768 | 0.071144 | 2.9188 | 2.9207 | 985.11 | 0.04371 | 919.64 | | 12 | Accept | 2.9471 | 0.63568 | 2.9188 | 2.8196 | 990.29 | 260 | 0.070438 | | 13 | Accept | 4.1768 | 0.079277 | 2.9188 | 2.8069 | 0.0033522 | 0.0010129 | 917.85 | | 14 | Accept | 2.9338 | 0.13213 | 2.9188 | 2.9234 | 982.04 | 976.47 | 0.056459 | | 15 | Accept | 3.0346 | 10.343 | 2.9188 | 2.9173 | 977.86 | 45.317 | 0.29352 | | 16 | Accept | 2.9385 | 0.1003 | 2.9188 | 2.9252 | 794.26 | 992.75 | 0.0097194 | | 17 | Accept | 4.1768 | 0.066994 | 2.9188 | 2.8937 | 0.0011152 | 0.18074 | 902.85 | | 18 | Best | 2.8849 | 0.10203 | 2.8849 | 2.9041 | 761.43 | 988.48 | 2.9861 | | 19 | Accept | 2.9152 | 0.42196 | 2.8849 | 2.8867 | 993.79 | 281.46 | 2.6184 | | 20 | Accept | 2.9531 | 0.066503 | 2.8849 | 2.8782 | 2.8244 | 996.96 | 1.7726 | |======================================================================================================= | Iter | Eval | Objective | Objective | BestSoFar | BestSoFar | BoxConstrain-| KernelScale | Epsilon | | | result | | runtime | (observed) | (estim.) | t | | | |======================================================================================================= | 21 | Accept | 2.941 | 0.15447 | 2.8849 | 2.8892 | 967.04 | 956.88 | 1.0803 | | 22 | Accept | 2.9527 | 5.7576 | 2.8849 | 2.8957 | 961.09 | 80.75 | 0.092167 | | 23 | Accept | 2.9427 | 2.4948 | 2.8849 | 2.8963 | 950.56 | 136.59 | 0.0094631 | | 24 | Accept | 2.9281 | 0.081381 | 2.8849 | 2.8971 | 6.2426 | 981.69 | 0.12581 |
| | | | | |
25 | Accept | 26 | Accept | 27 | Accept | 28 | Accept | 29 | Accept | 30 | Accept |
2.9153 | 3.0003 | 2.9514 | 2.9326 | 2.9464 | 2.9185 |
0.096715 | 2.8849 | 2.8965 | 79.343 | 507.59 | 2.258 | 0.093436 | 2.8849 | 2.9028 | 880.99 | 436.41 | 6.9026 | 0.27316 | 2.8849 | 2.9026 | 101.47 | 189.42 | 0.59609 | 0.21596 | 2.8849 | 2.9031 | 286.51 | 291 | 0.020006 | 4.097 | 2.8849 | 2.9037 | 988.08 | 103.96 | 0.024607 | 0.085007 | 2.8849 | 2.9046 | 262.95 | 913.07 | 1.9138 |
__________________________________________________________ Optimization completed. MaxObjectiveEvaluations of 30 reached. Total function evaluations: 30 Total elapsed time: 114.1812 seconds. Total objective function evaluation time: 41.8852 Best observed feasible point: BoxConstraint KernelScale Epsilon _____________ ___________ _______ 761.43
988.48
2.9861
Observed objective function value = 2.8849 Estimated objective function value = 2.9046 Function evaluation time = 0.10203 Best estimated feasible point (according to models): BoxConstraint KernelScale Epsilon _____________ ___________ _______ 262.95
913.07
1.9138
Estimated objective function value = 2.9046 Estimated function evaluation time = 0.084523 Mdl = RegressionSVM ResponseName: 'Y' CategoricalPredictors: [] ResponseTransform: 'none'
Alpha: [58×1 double] Bias: 45.7794 KernelParameters: [1×1 struct] NumObservations: 93 HyperparameterOptimizationResults: [1×1 BayesianOptimization] BoxConstraints: [93×1 double] ConvergenceInfo: [1×1 struct] IsSupportVector: [93×1 logical] Solver: 'SMO'
The optimization searched over BoxConstraint, KernelScale, and Epsilon. The output is the regression with the minimum estimated cross-validation loss.
4.1.2 predict Predict responses using support vector machine regression model
Syntax yfit = predict(Mdl,X)
Description yfit = predict(Mdl,X) returns a vector of predicted responses for the predictor data in the table or matrix X, based on the full or compact, trained support vector machine (SVM) regression model Mdl.
Input Arguments: Mdl — SVM regression model RegressionSVM object | CompactRegressionSVM object SVM regression model, specified as a RegressionSVM model or a CompactRegressionSVM model, returned by fitrsvm or compact, respectively. X — Predictor data used to generate responses numeric matrix | table Predictor data used to generate responses, specified as a numeric matrix or table. Each row of X corresponds to one observation, and each column corresponds to one variable. ·
For a numeric matrix: o The variables making up the columns of X must have the same order as the predictor variables that trained Mdl. o If you trained Mdl using a table (for example, Tbl), then X can be a numeric matrix if Tbl contains all numeric predictor variables. To treat numeric predictors in Tbl as categorical during training, identify categorical predictors using the CategoricalPredictors name-value pair argument of fitrsvm. If Tbl contains heterogeneous predictor variables (for example, numeric and categorical data types) and X is
·
a numeric matrix, then predict throws an error. For a table: o predict does not support multi-column variables and cell arrays other than cell arrays of character vectors. o If you trained Mdl using a table (for example, Tbl), then all predictor variables in X must have the same variable names and data types as those that trained Mdl (stored in Mdl.PredictorNames). However, the column order of X does not need to correspond to the column order of Tbl. Tbl and X can contain additional variables (response variables, observation weights, etc.), but predict ignores them. o If you trained Mdl using a numeric matrix, then the predictor names in Mdl.PredictorNames and corresponding predictor variable names in X must be the same. To specify predictor names during training, see the PredictorNames name-value pair argument of fitrsvm. All predictor variables in X must be numeric vectors. X can contain additional variables (response variables, observation weights, etc.), but predictignores them. If you set 'Standardize',true in fitrsvm to train Mdl, then the software standardizes the columns of X using the corresponding means in Mdl.Mu and standard deviations in Mdl.Sigma. Data Types: table | double | single
Output Arguments yfit — Predicted responses vector Predicted responses, returned as a vector of length n. n is the number of observations in the training data.
Examples: Predict Test Sample Response for SVM Regression
Model Load the carsmall data set. Consider a model that predicts a car's fuel efficiency given its horsepower and weight. Determine the sample size. load carsmall tbl = table(Horsepower,Weight,MPG); N = size(tbl,1); Partition the data into training and test sets. Hold out 10% of the data for testing. rng(10); % For reproducibility cvp = cvpartition(N,'Holdout',0.1); idxTrn = training(cvp); % Training set indices idxTest = test(cvp); % Test set indices Train a linear SVM regression model. Standardize the data. Mdl = fitrsvm(tbl(idxTrn,:),'MPG','Standardize',true); Mdl is a RegressionSVM model. Predict responses for the test set. YFit = predict(Mdl,tbl(idxTest,:)); Create a table containing the observed response values and the predicted response values side by side. table(tbl.MPG(idxTest),YFit,'VariableNames',... {'ObservedValue','PredictedValue'}) ans = 10×2 table ObservedValue PredictedValue _____________ ______________ 14 9.4833 27 28.938 10 7.765
28 22 29 24.5 18.5 32 28
27.155 21.054 31.484 30.306 19.12 28.225 26.632
4.1.3 fitrlinear Fit linear regression model to high-dimensional data fitrlinear efficiently trains linear regression models with high-dimensional, full or sparse predictor data. Available linear regression models include regularized support vector machines (SVM) and least-squares regression methods. fitrlinear minimizes the objective function using techniques that reduce computing time (e.g., stochastic gradient descent). A high-dimensional data set includes many predictor variables. Although such a data set can consume a significant fraction of memory, it must fit in the MATLAB® Workspace. For low- through medium-dimensional predictor data sets,
Syntax Mdl = fitrlinear(X,Y) Mdl = fitrlinear(X,Y,Name,Value) [Mdl,FitInfo] = fitrlinear(___) [Mdl,FitInfo,HyperparameterOptimizationResults] = fitrlinear(___)
Description Mdl = fitrlinear(X,Y) returns a trained regression model object that contains the results of fitting a support vector machine regression model to the predictors X and response Y. Mdl = fitrlinear(X,Y,Name,Value) returns a trained linear regression model with additional options specified by one or more Name,Value pair arguments. For example, you can specify implement least-squares regression, specify to cross-validate, or specify the type of regularization. It is good practice to cross-validate using the Kfold Name,Value pair argument. The crossvalidation results determine how well the model generalizes. [Mdl,FitInfo] = fitrlinear(___) also returns optimization details using any of the previous syntaxes. You cannot request FitInfo for cross-validated models.
[Mdl,FitInfo,HyperparameterOptimizationResults] = fitrlinear(___) also returns hyperparameter optimization details when you pass an OptimizeHyperparameters name-value pair.
Examples Train Linear Regression Model Train a linear regression model using SVM, dual SGD, and ridge regularization. Simulate 10000 observations from this model ·
is a 10000-by-1000 sparse matrix with 10% nonzero standard normal elements. · e is random normal error with mean 0 and standard deviation 0.3. rng(1) % For reproducibility n = 1e4; d = 1e3; nz = 0.1; X = sprandn(n,d,nz); Y = X(:,100) + 2*X(:,200) + 0.3*randn(n,1); Train a linear regression model. By default, fitrlinear uses support vector machines with a ridge penalty, and optimizes using dual SGD for SVM. Determine how well the optimization algorithm fit the model to the data by extracting a fit summary. [Mdl,FitInfo] = fitrlinear(X,Y) Mdl = RegressionLinear ResponseName: 'Y' ResponseTransform: 'none'
Beta: [1000×1 double] Bias: -0.0056 Lambda: 1.0000e-04 Learner: 'svm' FitInfo = struct with fields: Lambda: 1.0000e-04 Objective: 0.2726 PassLimit: 10 NumPasses: 10 BatchLimit: [] NumIterations: 100000 GradientNorm: NaN GradientTolerance: 0 RelativeChangeInBeta: 0.4908 BetaTolerance: 1.0000e-04 DeltaGradient: 1.5818 DeltaGradientTolerance: 0.1000 TerminationCode: 0 TerminationStatus: {'Maximal number of iterations has been reached.'} Alpha: [10000×1 double] History: [] FitTime: 0.0879 Solver: {'dual'} Mdl is a RegressionLinear model. You can pass Mdl and the training or new data to loss to inspect the in-sample mean-squared error. Or, you can pass Mdl and new predictor data to predict to predict responses for new observations. FitInfo is a structure array containing, among other things, the termination status (TerminationStatus) and how long the solver took to fit the model to the data (FitTime). It is good practice to use FitInfo to determine whether
optimization-termination measurements are satisfactory. In this case, fitrlinear reached the maximum number of iterations. Because training time is fast, you can retrain the model, but increase the number of passes through the data. Or, try another solver, such as LBFGS.
Find Good Lasso Penalty Using Cross-Validation To determine a good lasso-penalty strength for a linear regression model that uses least squares, implement 5-fold cross-validation. Simulate 10000 observations from this model ·
is a 10000-by-1000 sparse matrix with 10% nonzero standard normal elements. · e is random normal error with mean 0 and standard deviation 0.3. rng(1) % For reproducibility n = 1e4; d = 1e3; nz = 0.1; X = sprandn(n,d,nz); Y = X(:,100) + 2*X(:,200) + 0.3*randn(n,1); Create a set of 15 logarithmically-spaced regularization strengths from through . Lambda = logspace(-5,-1,15); Cross-validate the models. To increase execution speed, transpose the predictor data and specify that the observations are in columns. Optimize the objective function using SpaRSA. X = X'; CVMdl = fitrlinear(X,Y,'ObservationsIn','columns','KFold',5,'Lambda',Lambda,... 'Learner','leastsquares','Solver','sparsa','Regularization','lasso'); numCLModels = numel(CVMdl.Trained)
numCLModels = 5 CVMdl is a RegressionPartitionedLinear model. Because fitrlinear implements 5-fold cross-validation, CVMdlcontains 5 RegressionLinear models that the software trains on each fold. Display the first trained linear regression model. Mdl1 = CVMdl.Trained{1} Mdl1 = RegressionLinear ResponseName: 'Y' ResponseTransform: 'none' Beta: [1000×15 double] Bias: [1×15 double] Lambda: [1×15 double] Learner: 'leastsquares' Mdl1 is a RegressionLinear model object. fitrlinear constructed Mdl1 by training on the first four folds. Because Lambda is a sequence of regularization strengths, you can think of Mdl1 as 11 models, one for each regularization strength in Lambda. Estimate the cross-validated MSE. mse = kfoldLoss(CVMdl); Higher values of Lambda lead to predictor variable sparsity, which is a good quality of a regression model. For each regularization strength, train a linear regression model using the entire data set and the same options as when you cross-validated the models. Determine the number of nonzero coefficients per model. Mdl = fitrlinear(X,Y,'ObservationsIn','columns','Lambda',Lambda,... 'Learner','leastsquares','Solver','sparsa','Regularization','lasso');
numNZCoeff = sum(Mdl.Beta~=0); In the same figure, plot the cross-validated MSE and frequency of nonzero coefficients for each regularization strength. Plot all variables on the log scale. figure; [h,hL1,hL2] = plotyy(log10(Lambda),log10(mse),... log10(Lambda),log10(numNZCoeff)); hL1.Marker = 'o'; hL2.Marker = 'o'; ylabel(h(1),'log_{10} MSE') ylabel(h(2),'log_{10} nonzero-coefficient frequency') xlabel('log_{10} Lambda') hold off
Choose the index of the regularization strength that balances predictor variable sparsity and low MSE (for example, Lambda(10)). idxFinal = 10; Extract the model with corresponding to the minimal MSE. MdlFinal = selectModels(Mdl,idxFinal) idxNZCoeff = find(MdlFinal.Beta~=0) EstCoeff = Mdl.Beta(idxNZCoeff)
MdlFinal = RegressionLinear ResponseName: 'Y' ResponseTransform: 'none' Beta: [1000×1 double] Bias: -0.0050 Lambda: 0.0037 Learner: 'leastsquares' idxNZCoeff = 100 200
EstCoeff = 1.0051 1.9965 MdlFinal is a RegressionLinear model with one regularization strength. The nonzero coefficients EstCoeff are close to the coefficients that simulated the data.
Optimize a Linear Regression This example shows how to optimize hyperparameters automatically using fitrlinear. The example uses artificial (simulated) data for the model ·
is a 10000-by-1000 sparse matrix with 10% nonzero standard normal elements. · e is random normal error with mean 0 and standard deviation 0.3. rng(1) % For reproducibility n = 1e4;
d = 1e3; nz = 0.1; X = sprandn(n,d,nz); Y = X(:,100) + 2*X(:,200) + 0.3*randn(n,1); Find hyperparameters that minimize five-fold cross validation loss by using automatic hyperparameter optimization. For reproducibility, use the 'expected-improvement-plus' acquisition function. hyperopts = struct('AcquisitionFunctionName','expected-improvement-plus'); [Mdl,FitInfo,HyperparameterOptimizationResults] = fitrlinear(X,Y,... 'OptimizeHyperparameters','auto',... 'HyperparameterOptimizationOptions',hyperopts)
|=================================================================================================| | Iter | Eval | Objective | Objective | BestSoFar | BestSoFar | Lambda | Learner | | | result | | runtime | (observed) | (estim.) | | | |=================================================================================================| | 1 | Best | 0.10681 | 6.5103 | 0.10681 | 0.10681 | 0.016588 | svm | | 2 | Accept | 0.16238 | 0.82723 | 0.10681 | 0.11045 | 1.246e-08 | svm | | 3 | Accept | 0.13601 | 0.7158 | 0.10681 | 0.10831 | 1.5348e-07 | leastsquares | | 4 | Accept | 0.43611 | 0.53262 | 0.10681 | 0.21022 | 2.8871 | leastsquares | | 5 | Accept | 0.13088 | 0.88313 | 0.10681 | 0.10686 | 0.0027653 | svm | | 6 | Accept | 0.48011 | 0.53856 | 0.10681 | 0.10684 | 9.9895 | svm | | 7 | Accept | 0.13383 | 0.50301 | 0.10681 | 0.10684 | 7.6954e-09 | leastsquares | | 8 | Accept | 0.16061 | 0.70661 | 0.10681 | 0.10685 | 1.2282e-06 | svm | | 9 | Accept | 0.13436 | 0.49813 | 0.10681 | 0.10685 | 1.4623e-05 | leastsquares | | 10 | Accept | 0.13268 | 0.46243 | 0.10681 | 0.10685 | 0.00043155 | leastsquares | | 11 | Accept | 0.13393 | 0.57175 | 0.10681 | 0.10687 | 1.0066e-09 | leastsquares | | 12 | Accept | 0.16057 | 0.64287 | 0.10681 | 0.10687 | 1.0163e-09 | svm | | 13 | Accept | 0.16105 | 0.65268 | 0.10681 | 0.1069 | 4.105e-05 | svm | | 14 | Accept | 0.10682 | 0.79268 | 0.10681 | 0.10681 | 0.011705 | svm | | 15 | Accept | 0.13505 | 0.53991 | 0.10681 | 0.10678 | 9.3502e-05 | leastsquares | | 16 | Best | 0.10664 | 0.62679 | 0.10664 | 0.10654 | 0.01453 | svm | | 17 | Accept | 0.10684 | 0.69147 | 0.10664 | 0.10661 | 0.014328 | svm | | 18 | Accept | 0.10675 | 0.63836 | 0.10664 | 0.10664 | 0.014207 | svm | | 19 | Accept | 0.14223 | 0.51168 | 0.10664 | 0.10663 | 0.043736 | leastsquares | | 20 | Accept | 0.11133 | 0.54519 | 0.10664 | 0.10664 | 0.006649 | leastsquares | |=================================================================================================| | Iter | Eval | Objective | Objective | BestSoFar | BestSoFar | Lambda | Learner | | | result | | runtime | (observed) | (estim.) | | | |=================================================================================================| | 21 | Accept | 0.1067 | 0.63698 | 0.10664 | 0.10663 | 0.0158 | svm | | 22 | Accept | 0.41058 | 0.64622 | 0.10664 | 0.10648 | 0.54655 | svm | | 23 | Accept | 0.13203 | 0.57725 | 0.10664 | 0.10647 | 1.3948e-06 | leastsquares | | 24 | Accept | 0.12096 | 0.55958 | 0.10664 | 0.10647 | 0.0025687 | leastsquares | | 25 | Accept | 0.16124 | 0.64974 | 0.10664 | 0.10646 | 1.3205e-07 | svm | | 26 | Accept | 0.11057 | 0.48293 | 0.10664 | 0.10646 | 0.012941 | leastsquares |
| | | |
27 | Best | 0.10645 | 0.54893 | 0.10645 | 0.10646 | 0.012594 | svm | 28 | Accept | 0.10655 | 0.636 | 0.10645 | 0.10647 | 0.012839 | svm | 29 | Accept | 0.45006 | 0.48319 | 0.10645 | 0.10647 | 9.9931 | leastsquares | 30 | Accept | 0.15988 | 0.61594 | 0.10645 | 0.10648 | 0.00027628 | svm |
__________________________________________________________ Optimization completed. MaxObjectiveEvaluations of 30 reached. Total function evaluations: 30 Total elapsed time: 81.5182 seconds. Total objective function evaluation time: 24.228 Best observed feasible point: Lambda Learner ________ _______ 0.012594
svm
Observed objective function value = 0.10645 Estimated objective function value = 0.10648 Function evaluation time = 0.54893 Best estimated feasible point (according to models): Lambda Learner ________ _______ 0.014207 svm Estimated objective function value = 0.10648 Estimated function evaluation time = 0.71414 Mdl = RegressionLinear ResponseName: 'Y' ResponseTransform: 'none' Beta: [1000×1 double] Bias: -0.0029 Lambda: 0.0142
Learner: 'svm' FitInfo = struct with fields: Lambda: 0.0142 Objective: 0.2277 PassLimit: 10 NumPasses: 10 BatchLimit: [] NumIterations: 99993 GradientNorm: NaN GradientTolerance: 0 RelativeChangeInBeta: 0.0730 BetaTolerance: 1.0000e-04 DeltaGradient: 1.1375 DeltaGradientTolerance: 0.1000 TerminationCode: 0 TerminationStatus: {'Maximal number of iterations has been reached.'} Alpha: [10000×1 double] History: [] FitTime: 0.0881 Solver: {'dual'} HyperparameterOptimizationResults = BayesianOptimization with properties: ObjectiveFcn: @createObjFcn/theObjFcn VariableDescriptions: [3×1 optimizableVariable] Options: [1×1 struct] MinObjective: 0.1065 XAtMinObjective: [1×2 table] MinEstimatedObjective: 0.1065 XAtMinEstimatedObjective: [1×2 table] NumObjectiveEvaluations: 30
TotalElapsedTime: 81.5182 NextPoint: [1×2 table] XTrace: [30×2 table] ObjectiveTrace: [30×1 double] ConstraintsTrace: [] UserDataTrace: {30×1 cell} ObjectiveEvaluationTimeTrace: [30×1 double] IterationTimeTrace: [30×1 double] ErrorTrace: [30×1 double] FeasibilityTrace: [30×1 logical] FeasibilityProbabilityTrace: [30×1 double] IndexOfMinimumTrace: [30×1 double] ObjectiveMinimumTrace: [30×1 double] EstimatedObjectiveMinimumTrace: [30×1 double]
4.1.4 predict Predict response of linear regression model
Syntax YHat = predict(Mdl,X) YHat = predict(Mdl,X,Name,Value)
Description YHat = predict(Mdl,X) returns predicted responses for each observation in the predictor data Xbased on the trained linear regression model Mdl. YHat contains responses for each regularization strength in Mdl. YHat = predict(Mdl,X,Name,Value) returns predicted responses with additional options specified by one or more Name,Value pair arguments. For example, specify that columns in the predictor data correspond to observations
Examples: Predict Test-Sample Responses Simulate 10000 observations from this model ·
is a 10000-by-1000 sparse matrix with 10% nonzero standard normal elements. · e is random normal error with mean 0 and standard deviation 0.3. rng(1) % For reproducibility n = 1e4; d = 1e3; nz = 0.1; X = sprandn(n,d,nz); Y = X(:,100) + 2*X(:,200) + 0.3*randn(n,1); Train a linear regression model. Reserve 30% of the observations as a holdout sample. CVMdl = fitrlinear(X,Y,'Holdout',0.3);
Mdl = CVMdl.Trained{1} Mdl = RegressionLinear ResponseName: 'Y' ResponseTransform: 'none' Beta: [1000×1 double] Bias: -0.0066 Lambda: 1.4286e-04 Learner: 'svm' CVMdl is a RegressionPartitionedLinear model. It contains the property Trained, which is a 1-by-1 cell array holding a RegressionLinear model that the software trained using the training set. Extract the training and test data from the partition definition. trainIdx = training(CVMdl.Partition); testIdx = test(CVMdl.Partition); Predict the training- and test-sample responses. yHatTrain = predict(Mdl,X(trainIdx,:)); yHatTest = predict(Mdl,X(testIdx,:)); Because there is one regularization strength in Mdl, yHatTrain and yHatTest are numeric vectors.
Predict from Best-Performing Model Predict responses from the best-performing, linear regression model that uses a lasso-penalty and least squares. Simulate 10000 observations as in Predict Test-Sample Responses. rng(1) % For reproducibility n = 1e4; d = 1e3; nz = 0.1; X = sprandn(n,d,nz);
Y = X(:,100) + 2*X(:,200) + 0.3*randn(n,1); Create a set of 15 logarithmically-spaced regularization strengths from through . Lambda = logspace(-5,-1,15); Cross-validate the models. To increase execution speed, transpose the predictor data and specify that the observations are in columns. Optimizing the objective function using SpaRSA. X = X'; CVMdl = fitrlinear(X,Y,'ObservationsIn','columns','KFold',5,'Lambda',Lambda,... 'Learner','leastsquares','Solver','sparsa','Regularization','lasso'); numCLModels = numel(CVMdl.Trained) numCLModels = 5 CVMdl is a RegressionPartitionedLinear model. Because fitrlinear implements 5-fold cross-validation, CVMdlcontains 5 RegressionLinear models that the software trains on each fold. Display the first trained linear regression model. Mdl1 = CVMdl.Trained{1} Mdl1 = RegressionLinear ResponseName: 'Y' ResponseTransform: 'none' Beta: [1000×15 double] Bias: [1×15 double] Lambda: [1×15 double] Learner: 'leastsquares'
Mdl1 is a RegressionLinear model object. fitrlinear constructed Mdl1 by training on the first four folds. Because Lambda is a sequence of regularization strengths, you can think of Mdl1 as 11 models, one for each regularization strength in Lambda. Estimate the cross-validated MSE. mse = kfoldLoss(CVMdl); Higher values of Lambda lead to predictor variable sparsity, which is a good quality of a regression model. For each regularization strength, train a linear regression model using the entire data set and the same options as when you cross-validated the models. Determine the number of nonzero coefficients per model. Mdl = fitrlinear(X,Y,'ObservationsIn','columns','Lambda',Lambda,... 'Learner','leastsquares','Solver','sparsa','Regularization','lasso'); numNZCoeff = sum(Mdl.Beta~=0); In the same figure, plot the cross-validated MSE and frequency of nonzero coefficients for each regularization strength. Plot all variables on the log scale. figure; [h,hL1,hL2] = plotyy(log10(Lambda),log10(mse),... log10(Lambda),log10(numNZCoeff)); hL1.Marker = 'o'; hL2.Marker = 'o'; ylabel(h(1),'log_{10} MSE') ylabel(h(2),'log_{10} nonzero-coefficient frequency') xlabel('log_{10} Lambda') hold off
Choose the index of the regularization strength that balances predictor variable sparsity and low MSE (for example, Lambda(10)). idxFinal = 10; Extract the model with corresponding to the minimal MSE. MdlFinal = selectModels(Mdl,idxFinal) idxNZCoeff = find(MdlFinal.Beta~=0) EstCoeff = Mdl.Beta(idxNZCoeff) MdlFinal = RegressionLinear ResponseName: 'Y' ResponseTransform: 'none' Beta: [1000×1 double] Bias: -0.0050 Lambda: 0.0037 Learner: 'leastsquares' idxNZCoeff = 100
200 EstCoeff = 1.0051 1.9965 MdlFinal is a RegressionLinear model with one regularization strength. The nonzero coefficients EstCoeff are close to the coefficients that simulated the data. Simulate 10 new observations, and predict corresponding responses using the best-performing model. XNew = sprandn(d,10,nz); YHat = predict(MdlFinal,XNew,'ObservationsIn','columns');
4.2 GAUSSIAN PROCESS REGRESSION MODELS Gaussian process regression (GPR) models are nonparametric kernel-based probabilistic models. You can train a GPR model using the fitrgp function.
4.2.1 fitrgp Fit a Gaussian process regression (GPR) model
Syntax gprMdl = fitrgp(tbl,ResponseVarName) gprMdl = fitrgp(tbl,formula) gprMdl = fitrgp(tbl,y) gprMdl = fitrgp(X,y) gprMdl = fitrgp(___,Name,Value)
Description gprMdl = fitrgp(tbl,ResponseVarName) returns a Gaussian process regression (GPR) model trained using the sample data in tbl, where ResponseVarName is the name of the response variable in tbl. gprMdl = fitrgp(tbl,formula) returns a Gaussian process regression (GPR) model, trained using the sample data in tbl, for the predictor variables and response variables identified by formula. gprMdl = fitrgp(tbl,y) returns a GPR model for the predictors in table tbl and continuous response vector y. gprMdl = fitrgp(X,y) returns a GPR model for predictors X and continuous response vector y. gprMdl = fitrgp(___,Name,Value) returns a GPR model for any of the input arguments in the previous syntaxes, with additional options specified by one or more Name,Value pair arguments. For example, you can specify the fitting method, the prediction method, the covariance function, or the active set selection method. You can also train a cross-validated model. gprMdl is a RegressionGP object. For methods and properties of this class, see RegressionGP class page. If you train a cross-validated model, then gprMdl is a RegressionPartitionedModel object. For further analysis on the crossvalidated object, use the methods of RegressionPartitionedModel class. For
the methods of this class, see the RegressionPartitionedModel class page.
Examples: Train GPR Model Using Data in Table This example uses the abalone data , from the UCI Machine Learning Repository . Download the data and save it in your current folder with the name ‘abalone.data'. Store the data into a table. Display the first seven rows.
tbl = readtable('abalone.data','Filetype','text','ReadVariableNames',false);tbl.Properties.VariableN = {'Sex','Length','Diameter','Height','WWeight','SWeight','VWeight','ShWeight','NoShellRings'} tbl(1:7,:) ans = Sex Length Diameter Height WWeight SWeight VWeight ShWeight NoShellRings ___ ______ ________ ______ _______ _______ _______ ________ ____________ 'M' 0.455 0.365 0.095 0.514 0.2245 0.101 0.15 15 'M' 0.35 0.265 0.09 0.2255 0.0995 0.0485 0.07 7 'F' 0.53 0.42 0.135 0.677 0.2565 0.1415 0.21 9 'M' 0.44 0.365 0.125 0.516 0.2155 0.114 0.155 10 'I' 0.33 0.255 0.08 0.205 0.0895 0.0395 0.055 7 'I' 0.425 0.3 0.095 0.3515 0.141 0.0775 0.12 8 'F' 0.53 0.415 0.15 0.7775 0.237 0.1415 0.33 20
The dataset has 4177 observations. The goal is to predict the age of abalone from eight physical measurements. The last variable, number of shell rings shows the age of the abalone. The first predictor is a categorical variable. The last variable in the table is the response variable. Fit a GPR model using the subset of regressors method for parameter estimation and fully independent conditional method for prediction. Standardize the predictors. gprMdl = fitrgp(tbl,'NoShellRings','KernelFunction','ardsquaredexponential',... 'FitMethod','sr','PredictMethod','fic','Standardize',1) rMdl =
RegressionGP PredictorNames: {1x8 cell} ResponseName: 'Var9' ResponseTransform: 'none' NumObservations: 4177 KernelFunction: 'ARDSquaredExponential' KernelInformation: [1x1 struct] BasisFunction: 'Constant' Beta: 10.9148 Sigma: 2.0243 PredictorLocation: [10x1 double] PredictorScale: [10x1 double] Alpha: [1000x1 double] ActiveSetVectors: [1000x10 double] PredictMethod: 'FIC' ActiveSetSize: 1000 FitMethod: 'SR' ActiveSetMethod: 'Random' IsActiveSetVector: [4177x1 logical] LogLikelihood: -9.0013e+03 ActiveSetHistory: [1x1 struct] BCDInformation: [] Predict the responses using the trained model. ypred = resubPredict(gprMdl); Plot the true response and the predicted responses. figure(); plot(tbl.NoShellRings,'r.'); hold on plot(ypred,'b'); xlabel('x'); ylabel('y'); legend({'data','predictions'},'Location','Best');
axis([0 4300 0 30]); hold off;
Compute the regression loss on the training data (resubstitution loss) for the trained model. L = resubLoss(gprMdl) L= 4.0064
Train GPR Model and Plot Predictions Generate sample data. rng(0,'twister'); % For reproducibility n = 1000; x = linspace(-10,10,n)'; y = 1 + x*5e-2 + sin(x)./x + 0.2*randn(n,1); Fit a GPR model using a linear basis function and the exact fitting method to estimate the parameters. Also use the exact prediction method. gprMdl = fitrgp(x,y,'Basis','linear',... 'FitMethod','exact','PredictMethod','exact'); Predict the response corresponding to the rows of x (resubstitution predictions) using the trained model.
ypred = resubPredict(gprMdl); Plot the true response with the predicted values. plot(x,y,'b.'); hold on; plot(x,ypred,'r','LineWidth',1.5); xlabel('x'); ylabel('y'); legend('Data','GPR predictions'); hold off
Impact of Specifying Initial Kernel Parameter Values Load the sample data. load(fullfile(matlabroot,'examples','stats','gprdata2.mat')) The data has one predictor variable and continuous response. This is simulated data. Fit a GPR model using the squared exponential kernel function with default kernel parameters. gprMdl1 = fitrgp(x,y,'KernelFunction','squaredexponential');
Now, fit a second model, where you specify the initial values for the kernel parameters. sigma0 = 0.2; kparams0 = [3.5, 6.2]; gprMdl2 = fitrgp(x,y,'KernelFunction','squaredexponential',... 'KernelParameters',kparams0,'Sigma',sigma0); Compute the resubstitution predictions from both models. ypred1 = resubPredict(gprMdl1); ypred2 = resubPredict(gprMdl2); Plot the response predictions from both models and the responses in training data. figure(); plot(x,y,'r.'); hold on plot(x,ypred1,'b'); plot(x,ypred2,'g'); xlabel('x'); ylabel('y'); legend({'data','default kernel parameters',... 'kparams0 = [3.5,6.2], sigma0 = 0.2'},... 'Location','Best'); title('Impact of initial kernel parameter values'); hold off
The marginal log likelihood that fitrgp maximizes to estimate GPR parameters has multiple local solutions; the solution that it converges to depends on the initial point. Each local solution corresponds to a particular interpretation of the data. In this example, the solution with the default initial kernel parameters corresponds to a low frequency signal with high noise whereas the second solution with custom initial kernel parameters corresponds to a high frequency signal with low noise.
Use Separate Length Scales for Predictors Load the sample data. load(fullfile(matlabroot,'examples','stats','gprdata.mat')) There are six continuous predictor variables. There are 500 observations in the training data set and 100 observations in the test data set. This is simulated data. Fit a GPR model using the squared exponential kernel function with a separate length scale for each predictor. This covariance function is defined as:
where represents the length scale for predictor , = 1, 2, ..., and the signal standard deviation. The unconstrained parametrization is
is
Initialize length scales of the kernel function at 10 and signal and noise standard deviations at the standard deviation of the response. sigma0 = std(ytrain); sigmaF0 = sigma0; d = size(Xtrain,2); sigmaM0 = 10*ones(d,1); Fit the GPR model using the initial kernel parameter values. Standardize the predictors in the training data. Use the exact fitting and prediction methods. gprMdl = fitrgp(Xtrain,ytrain,'Basis','constant','FitMethod','exact',... 'PredictMethod','exact','KernelFunction','ardsquaredexponential',... 'KernelParameters',[sigmaM0;sigmaF0],'Sigma',sigma0,'Standardize',1); Compute the regression loss on the test data. L = loss(gprMdl,Xtest,ytest) L= 0.6919 Access the kernel information. gprMdl.KernelInformation ans = struct with fields: Name: 'ARDSquaredExponential' KernelParameters: [7×1 double] KernelParameterNames: {7×1 cell} Display the kernel parameter names.
gprMdl.KernelInformation.KernelParameterNames ans = 7×1 cell array 'LengthScale1' 'LengthScale2' 'LengthScale3' 'LengthScale4' 'LengthScale5' 'LengthScale6' 'SigmaF' Display the kernel parameters. sigmaM = gprMdl.KernelInformation.KernelParameters(1:end-1,1) sigmaF = gprMdl.KernelInformation.KernelParameters(end) sigma = gprMdl.Sigma sigmaM = 1.0e+04 * 0.0004 0.0007 0.0004 4.1729 0.1018 0.0056 sigmaF = 28.1718 sigma = 0.8162
Plot the log of learned length scales. figure() plot((1:d)',log(sigmaM),'ro-'); xlabel('Length scale number'); ylabel('Log of length scale');
The log of length scale for the 4th and 5th predictor variables are high relative to the others. These predictor variables do not seem to be as influential on the response as the other predictor variables. Fit the GPR model without using the 4th and 5th variables as the predictor variables. X = [Xtrain(:,1:3) Xtrain(:,6)]; sigma0 = std(ytrain); sigmaF0 = sigma0; d = size(X,2); sigmaM0 = 10*ones(d,1); gprMdl = fitrgp(X,ytrain,'Basis','constant','FitMethod','exact',... 'PredictMethod','exact','KernelFunction','ardsquaredexponential',...
'KernelParameters',[sigmaM0;sigmaF0],'Sigma',sigma0,'Standardize',1); Compute the regression error on the test data. xtest = [Xtest(:,1:3) Xtest(:,6)]; L = loss(gprMdl,xtest,ytest) L= 0.6928 The loss is similar to the one when all variables are used as predictor variables. Compute the predicted response for the test data. ypred = predict(gprMdl,xtest); Plot the original response along with the fitted values. figure; plot(ytest,'r'); hold on; plot(ypred,'b'); legend('True response','GPR predicted values','Location','Best'); hold off
Optimize GPR Regression This example shows how to optimize hyperparameters automatically using fitrgp. The example uses the gprdata2 data that ships with your software. Load the data. load(fullfile(matlabroot,'examples','stats','gprdata2.mat')) The data has one predictor variable and continuous response. This is simulated data. Fit a GPR model using the squared exponential kernel function with default kernel parameters. gprMdl1 = fitrgp(x,y,'KernelFunction','squaredexponential'); Find hyperparameters that minimize five-fold cross-validation loss by using automatic hyperparameter optimization. For reproducibility, set the random seed and use the 'expected-improvementplus' acquisition function. rng default gprMdl2 = fitrgp(x,y,'KernelFunction','squaredexponential',...
'OptimizeHyperparameters','auto','HyperparameterOptimizationOptions',... struct('AcquisitionFunctionName','expected-improvement-plus')); |==================================================================================| | Iter | Eval | Objective | Objective | BestSoFar | BestSoFar | Sigma | | | result | | runtime | (observed) | (estim.) | | |==================================================================================| | 1 | Best | 0.44511 | 1.5781 | 0.44511 | 0.44511 | 8.5486 | | 2 | Best | 0.31949 | 1.9102 | 0.31949 | 0.32719 | 0.009463 | | 3 | Accept | 0.31949 | 1.998 | 0.31949 | 0.31949 | 0.00027981 | | 4 | Best | 0.13146 | 1.4224 | 0.13146 | 0.14737 | 0.39244 | | 5 | Accept | 0.21452 | 1.5708 | 0.13146 | 0.17013 | 0.44668 | | 6 | Best | 0.12927 | 1.4485 | 0.12927 | 0.12927 | 0.35388 | | 7 | Best | 0.037962 | 1.5693 | 0.037962 | 0.050194 | 0.27207 | | 8 | Best | 0.037893 | 1.7106 | 0.037893 | 0.038108 | 0.20304 | | 9 | Accept | 0.037906 | 1.4331 | 0.037893 | 0.034856 | 0.22437 | | 10 | Accept | 0.037938 | 1.4814 | 0.037893 | 0.030237 | 0.12373 | | 11 | Best | 0.03789 | 1.5218 | 0.03789 | 0.033148 | 0.16967 | | 12 | Accept | 0.2237 | 1.5856 | 0.03789 | 0.033458 | 0.065552 | | 13 | Accept | 0.037893 | 1.5801 | 0.03789 | 0.034842 | 0.16041 | | 14 | Accept | 0.037897 | 1.5758 | 0.03789 | 0.036168 | 0.15318 | | 15 | Accept | 0.0379 | 1.6014 | 0.03789 | 0.03611 | 0.21631 | | 16 | Accept | 0.037904 | 1.3969 | 0.03789 | 0.036677 | 0.14487 | | 17 | Accept | 0.037907 | 1.5767 | 0.03789 | 0.036582 | 0.22534 | | 18 | Accept | 0.037906 | 1.6101 | 0.03789 | 0.037004 | 0.14245 | | 19 | Accept | 0.03791 | 1.7306 | 0.03789 | 0.036729 | 0.22943 | | 20 | Accept | 0.037907 | 1.5169 | 0.03789 | 0.036716 | 0.1414 | |==================================================================================| | Iter | Eval | Objective | Objective | BestSoFar | BestSoFar | Sigma | | | result | | runtime | (observed) | (estim.) | | |==================================================================================| | 21 | Accept | 0.31949 | 1.9515 | 0.03789 | 0.036727 | 0.0016322 | | 22 | Accept | 2.0999 | 1.322 | 0.03789 | 0.036665 | 29.262 | | 23 | Accept | 0.42096 | 1.3101 | 0.03789 | 0.036678 | 3.3239 | | 24 | Accept | 0.31949 | 1.9251 | 0.03789 | 0.03668 | 0.00010003 | | 25 | Accept | 0.31949 | 2.1799 | 0.03789 | 0.036698 | 0.0039898 | | 26 | Accept | 0.31949 | 2.0027 | 0.03789 | 0.036723 | 0.00066626 | | 27 | Accept | 0.31949 | 2.0363 | 0.03789 | 0.036687 | 0.023583 | | 28 | Accept | 0.42073 | 1.4276 | 0.03789 | 0.036663 | 1.4756 | | 29 | Accept | 0.31949 | 2.0018 | 0.03789 | 0.036562 | 0.00015967 | | 30 | Accept | 0.63932 | 1.0931 | 0.03789 | 0.037714 | 15.25 | __________________________________________________________
Optimization completed. MaxObjectiveEvaluations of 30 reached. Total function evaluations: 30
Total elapsed time: 80.2692 seconds. Total objective function evaluation time: 49.0683 Best observed feasible point: Sigma _______ 0.16967 Observed objective function value = 0.03789 Estimated objective function value = 0.037714 Function evaluation time = 1.5218 Best estimated feasible point (according to models): Sigma _______ 0.22534 Estimated objective function value = 0.037714 Estimated function evaluation time = 1.5446
Compare the pre- and post-optimization fits. ypred1 = resubPredict(gprMdl1); ypred2 = resubPredict(gprMdl2);
figure(); plot(x,y,'r.'); hold on plot(x,ypred1,'b'); plot(x,ypred2,'k','LineWidth',2); xlabel('x'); ylabel('y'); legend({'data','Initial Fit','Optimized Fit'},'Location','Best'); title('Impact of Optimization'); hold off
Train GPR Model Using Cross-Validation This example uses the abalone data [1], [2], from the UCI Machine Learning Repository [3]. Download the data and save it in your current folder with the name ‘abalone.data'. Store the data into a table. Display the first seven rows.
tbl = readtable('abalone.data','Filetype','text','ReadVariableNames',false);tbl.Properties.VariableN =
{'Sex','Length','Diameter','Height','WWeight','SWeight','VWeight','ShWeight','NoShellRings'} tbl(1:7,:) ans = Sex Length Diameter Height WWeight SWeight VWeight ShWeight NoShellRings ___ ______ ________ ______ _______ _______ _______ ________ ____________ 'M' 0.455 0.365 0.095 0.514 0.2245 0.101 0.15 15 'M' 0.35 0.265 0.09 0.2255 0.0995 0.0485 0.07 7 'F' 0.53 0.42 0.135 0.677 0.2565 0.1415 0.21 9 'M' 0.44 0.365 0.125 0.516 0.2155 0.114 0.155 10 'I' 0.33 0.255 0.08 0.205 0.0895 0.0395 0.055 7 'I' 0.425 0.3 0.095 0.3515 0.141 0.0775 0.12 8 'F' 0.53 0.415 0.15 0.7775 0.237 0.1415 0.33 20
The dataset has 4177 observations. The goal is to predict the age of abalone from eight physical measurements. The last variable, number of shell rings shows the age of the abalone. The first predictor is a categorical variable. The last variable in the table is the response variable. Train a cross-validated GPR model using the 25% of the data for validation. rng('default') % For reproducibility cvgprMdl = fitrgp(tbl,'NoShellRings','Standardize',1,'Holdout',0.25); Compute the average loss on folds using models trained on out-of-fold observations. kfoldLoss(cvgprMdl) ans = 4.6409 Predict the responses for out-of-fold data. ypred = kfoldPredict(cvgprMdl); Plot the true responses used for testing and the predictions. figure(); plot(ypred(cvgprMdl.Partition.test)); hold on; y = table2array(tbl(:,end)); plot(y(cvgprMdl.Partition.test),'r.');
axis([0 1050 0 30]); xlabel('x') ylabel('y') hold off;
Fit GPR Model Using Custom Kernel Function Generate the sample data. rng(0,'twister'); % For reproducibility n = 1000; x = linspace(-10,10,n)'; y = 1 + x*5e-2 + sin(x)./x + 0.2*randn(n,1); Define the squared exponential kernel function as a custom kernel function. You can compute the squared exponential kernel function as
where is the signal standard deviation, is the length scale. Both and must be greater than zero. This condition can be enforced by the unconstrained
parametrization, and , for some unconstrained parametrization vector . Hence, you can define the squared exponential kernel function as a custom kernel function as follows: kfcn = @(XN,XM,theta) (exp(theta(2))^2)*exp((pdist2(XN,XM).^2)/(2*exp(theta(1))^2)); Here pdist2(XN,XM).^2 computes the distance matrix. Fit a GPR model using the custom kernel function, kfcn. Specify the initial values of the kernel parameters (Because you use a custom kernel function, you must provide initial values for the unconstrained parametrization vector, theta). theta0 = [1.5,0.2]; gprMdl = fitrgp(x,y,'KernelFunction',kfcn,'KernelParameters',theta0); fitrgp uses analytical derivatives to estimate parameters when using a built-in kernel function, whereas when using a custom kernel function it uses numerical derivatives. Compute the resubstitution loss for this model. L = resubLoss(gprMdl) L= 0.0391 Fit the GPR model using the built-in squared exponential kernel function option. Specify the initial values of the kernel parameters (Because you use the built-in custom kernel function and specifying initial parameter values, you must provide the initial values for the signal standard deviation and length scale(s) directly). sigmaL0 = exp(1.5); sigmaF0 = exp(0.2); gprMdl2 = fitrgp(x,y,'KernelFunction','squaredexponential','KernelParameters', [sigmaL0,sigmaF0]); Compute the resubstitution loss for this model. L2 = resubLoss(gprMdl2) L2 =
0.0391 The two loss values are the same as expected.
4.2.2 predict Predict response of Gaussian process regression model
Syntax ypred = predict(gprMdl,Xnew) [ypred,ysd] = predict(gprMdl,Xnew) [ypred,ysd,yint] = predict(gprMdl,Xnew) [ypred,ysd,yint] = predict(gprMdl,Xnew,Name,Value)
Description ypred = predict(gprMdl,Xnew) returns the predicted responses ypred for the full or compact Gaussian process regression (GPR) model, gprMdl, and the predictor values in Xnew. [ypred,ysd] = predict(gprMdl,Xnew) also returns the estimated standard deviations for the new responses at the predictor values in Xnew from a trained GPR model. [ypred,ysd,yint] = predict(gprMdl,Xnew) also returns the 95% prediction intervals, yint, for the true responses corresponding to each row of Xnew. [ypred,ysd,yint] = predict(gprMdl,Xnew,Name,Value) also returns the prediction intervals with additional options specified by one or more Name,Value pair arguments. For example, you can specify the confidence level of the prediction interval.
Examples Compute Predicted Responses Generate the sample data. n = 10000; rng(1) % For reproducibility x = linspace(0.5,2.5,n)'; y = sin(10*pi.*x) ./ (2.*x)+(x-1).^4 + 1.5*rand(n,1); Fit a GPR model using the Matern 3/2 kernel function with separate length
scale for each predictor and an active set size of 100. Use the subset of regressors approximation method for parameter estimation and fully independent conditional method for prediction. gprMdl = fitrgp(x,y,'KernelFunction','ardmatern32',... 'ActiveSetSize',100,'FitMethod','sr','PredictMethod','fic'); Compute the predictions. [ypred,~,yci] = predict(gprMdl,x); Plot the data along with the predictions and prediction intervals. plot(x,y,'r.'); hold on plot(x,ypred); plot(x,yci(:,1),'k--'); plot(x,yci(:,2),'k--'); xlabel('x'); ylabel('y');
Compute Predictions When Data in Table Load the sample data and store in a table.
load fisheriris tbl = table(meas(:,1),meas(:,2),meas(:,3),meas(:,4),species,... 'VariableNames',{'meas1','meas2','meas3','meas4','species'}); Fit a GPR model using the first measurement as the response and the other variables as the predictors. mdl = fitrgp(tbl,'meas1'); Compute the predictions and the 99% confidence intervals. [ypred,~,yci] = predict(mdl,tbl,'Alpha',0.01); Plot the true response and the predictions along with the prediction intervals. figure(); plot(mdl.Y,'r.'); hold on; plot(ypred); plot(yci(:,1),'k:'); plot(yci(:,2),'k:'); legend('True response','GPR predictions',... 'Lower prediction limit','Upper prediction limit',... 'Location','Best');
Plot Predicted Response for Test Data Navigate to a folder containing sample data. cd(matlabroot) cd('help/toolbox/stats/examples') Load the sample data. load gprdata The data contains training and test data. There are 500 observations in training data and 100 observations in test data. The data has 8 predictor variables. This is simulated data. Fit a GPR model using the squared exponential kernel function with a separate length scale for each predictor. Standardize predictors in the training data. Use the exact fitting and prediction methods. gprMdl = fitrgp(Xtrain,ytrain,'Basis','constant',... 'FitMethod','exact','PredictMethod','exact',... 'KernelFunction','ardsquaredexponential','Standardize',1); Predict the responses for test data. [ytestpred,~,ytestci] = predict(gprMdl,Xtest); Plot the test response along with the predictions. figure; plot(ytest,'r'); hold on; plot(ytestpred,'b'); plot(ytestci(:,1),'k:'); plot(ytestci(:,2),'k:'); legend('Actual response','GPR predictions',... '95% lower','95% upper','Location','Best'); hold off
Tips ·
You can choose the prediction method while training the GPR model using the PredictMethod name-value pair argument in fitrgp. The default prediction method is 'exact' for n ≤ 10000, where n is the number of observations in the training data, and 'bcd' (block coordinate descent), otherwise. · Computation of standard deviations, ysd, and prediction intervals, yint, is not supported when PredictMethod is 'bcd'. · If gprMdl is a CompactRegressionGP object, you cannot compute standard deviations, ysd, or prediction intervals, yint, for PredictMethod equal to 'sr' or 'fic'. To compute ysd and yint for PredictMethod equal to 'sr' or 'fic', use the full regression (RegressionGP) object.
4.2.3 postFitStatistics Compute post-fit statistics for the exact Gaussian process regression model
Syntax loores = postFitStatistics(gprMdl) [loores,neff] = postFitStatistics(gprMdl)
Description loores = postFitStatistics(gprMdl) returns the leave-one-out residuals, loores, for the trained Gaussian process regression (GPR) model. [loores,neff] = postFitStatistics(gprMdl) also returns the number of effective parameters, neff.
Examples Compute Post-Fit Statistics Generate sample data. rng(0,'twister'); % For reproducibility n = 1500; x = linspace(-10,10,n)'; y = sin(3*x).*cos(3*x) + sin(2*x).*cos(2*x) + sin(x) + cos(x) + 0.2*randn(n,1); Fit a GPR model using the exact method for fitting and prediction. gprMdl = fitrgp(x,y,'Basis','linear','FitMethod','exact',... 'PredictMethod','exact','KernelFunction','matern52'); Compute the leave-one-out residuals and the effective number of parameters in the trained model. [loores,neff] = postFitStatistics(gprMdl); Predict the responses using the trained model. ypred = resubPredict(gprMdl);
Plot the true and predicted responses, and display effective number of parameters in the fit. figure() plot(x,y,'r.'); hold on; plot(x,ypred,'b'); xlabel('x'); ylabel('y'); legend('Data','GPR prediction','Location','Best'); title(['Effective number of parameters = ',num2str(neff)]); hold off
Plot leave-one-out residuals. figure() plot(x,loores,'r.-'); xlabel('x'); ylabel('leave-one-out residuals');
Tips ·
You can only compute the post-fit statistics when PredictMethod is 'exact'. · If FitMethod is 'exact', then postFitStatistics accounts for the fact that the fixed basis function coefficients are estimated from the data. · If FitMethod is different than 'exact', then postFitStatistics treats the fixed basis function coefficients as known. · For all PredictMethod and FitMethod options, postFitStatistics treats the estimated kernel parameters and noise standard deviation as known.
