www.gusucode.com > stats 源码程序 matlab案例代码 > stats/TuneRegularizationParameterForNCAForRegressionExample.m

    %% Tune Regularization Parameter in NCA for Regression
% Load the sample data.
load(fullfile(matlabroot,'examples','stats','robotarm.mat'))

%%
% The robotarm (pumadyn32nm) dataset is created using a robot arm simulator
% with 7168 training and 1024 test observations with 32 features [1], [2].
% This is a preprocessed version of the original data set. Data are
% preprocessed by subtracting off a linear regression fit followed by
% normalization of all features to unit variance.
%%
% Perform neighborhood component analysis feature selection for regression with
% default $\lambda$ (regularization parameter) value.

nca = fsrnca(Xtrain,ytrain,'FitMethod','exact',...
    'Solver','lbfgs');

%%
% Plot the selected values.
figure;
plot(nca.FeatureWeights,'ro');
xlabel('Feature index');
ylabel('Feature weight');
grid on;

%%
% Most of the feature weights are nonzero. Compute the loss using the test
% set as a measure of the performance using the selected features.
L = loss(nca,Xtest,ytest)

%%
% Let's see if the performance can be improved. Tune the regularization
% parameter $\lambda$ for feature selection using five-fold
% cross-validation. Tuning $\lambda$ means finding the $\lambda$ value that
% will produce the minimum regression loss. Here are the steps for
% tuning $\lambda$ using cross-validation:
%
% 1. First partition the data into five folds. For each fold, |cvpartition|
% assigns 4/5th of the data as a training set, and 1/5th of the data as
% a test set. 
rng(1) % For reproducibility 
n           = length(ytrain);
cvp         = cvpartition(length(ytrain),'kfold',5);
numvalidsets = cvp.NumTestSets;

%%
% Assign the lambda values for the search. Create an array to store the
% loss values.
lambdavals  = linspace(0,40,20)*std(ytrain)/n;
lossvals    = zeros(length(lambdavals),numvalidsets);

%%
% 2. Train the neighborhood component analysis (nca) model for each
% $\lambda$ value using the training set in each fold.
%
% 3. Compute the regression loss for the corresponding test
% set in the fold using the nca model. Record the loss value. 
%%
% 4. Repeat this for each $\lambda$ value and each fold.
for i = 1:length(lambdavals)
    for k = 1:numvalidsets
        X = Xtrain(cvp.training(k),:);
        y = ytrain(cvp.training(k),:);
        Xvalid  = Xtrain(cvp.test(k),:);
        yvalid  = ytrain(cvp.test(k),:);

        nca = fsrnca(X,y,'FitMethod','exact',...
             'Solver','lbfgs','Lambda',lambdavals(i));
         
        lossvals(i,k) = loss(nca,Xvalid,yvalid,'LossFunction','mse');
    end
end

%%
% Compute the average loss obtained from the folds for each $\lambda$ value.
meanloss = mean(lossvals,2);

%%
% Plot the mean loss versus the $\lambda$ values.
figure;
plot(lambdavals,meanloss,'ro-');
xlabel('Lambda');
ylabel('Loss (MSE)');
grid on;

%%
% Find the $\lambda$ value that gives the minimum loss value.
[~,idx] = min(meanloss)

bestlambda = lambdavals(idx)
bestloss   = meanloss(idx)
bestrmse   = sqrt(bestloss)

%%
% Fit the neighborhood component feature selection for regression using the
% best $\lambda$ value.
nca = fsrnca(Xtrain,ytrain,'FitMethod','exact',...
    'Solver','lbfgs','Lambda',bestlambda);

%%
% Plot the selected features.
figure;
plot(nca.FeatureWeights,'ro');
xlabel('Feature index');
ylabel('Feature weight');
grid on;

%%
% Most of the feature weights are zero. |fsrnca| identifies the most
% relevant four features.
%%
% Compute the loss for the test set.
L = loss(nca,Xtest,ytest)

%%
% Tuning the regularization parameter $\lambda$ eliminated more of the irrelevant
% features and improved the performance.