www.pudn.com > WebsiteCode.zip > polynomial.m, change:2008-07-24,size:3072b

function [BestOrder,Coeff,MU] = polynomial(X,Y) % Fits a one-variable polynomial to one-dimensional data % % Inputs: % X,Y - training data vectors % % Outputs: % BestOrder - the order of the polynomial, estimated using % cross-validation % Coeff - the coefficients of the polynomial % MU - normalisation vector % % Copyright 2007 A Sobester % % This program is free software: you can redistribute it and/or modify it % under the terms of the GNU Lesser General Public License as published by % the Free Software Foundation, either version 3 of the License, or any % later version. % % This program is distributed in the hope that it will be useful, but % WITHOUT ANY WARRANTY; without even the implied warranty of % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser % General Public License for more details. % % You should have received a copy of the GNU General Public License and GNU % Lesser General Public License along with this program. If not, see % <http://www.gnu.org/licenses/>. PlotPoints = 100; OrderLabel =... {'Linear',... 'Quadratic',... 'Cubic',... 'Quartic',... 'Fifth order',... 'Sixth order',... 'Seventh order',... 'Eighth order',... 'Ninth order',... 'Tenth order',... 'Eleventh order',... 'Twelveth order',... 'Thirteenth order',... 'Fourteenth order',... 'Fifteenth order'}; % The crossvalidation will split the data into this many subsets % (this may be changed if required) q = 5; % This is the highest order polynomial that will be considered MaxOrder = 15; n = length(X); % X split into q randomly selected subsets XS = randperm(n); FullXS = XS; % The beginnings of the subsets... From = (1:round(n/q):n-1); To = zeros(size(From)); % ...and their ends for i=1:q-1 To(i) = From(i+1)-1; end To(q) = n; CrossVal = zeros(1,MaxOrder); % Cycling through the possible orders for Order = 1:MaxOrder fprintf('Now considering the polynomial of order %d...\n',Order); CrossVal(Order) = 0; % Model fitting to subsets of the data for j = 1:q Removed = XS(From(j):To(j)); XS(From(j):To(j)) = []; [P,S,MU] = polyfit(X(XS),Y(XS),Order); CrossVal(Order) = CrossVal(Order) +... sum((Y(Removed) - polyval(P,X(Removed),S,MU)).^2)... /length(Removed); XS = FullXS; end end [MinCV, BestOrder] = min(CrossVal); [Coeff,S,MU] = polyfit(X,Y,BestOrder); % **** VISUALISING THE TRAINING DATA AND THE PREDICTOR ***** XPrediction = (min(X):(max(X)-min(X))/(PlotPoints-1):max(X)); [YPrediction,DELTA] = polyval(Coeff,XPrediction,S,MU); plot(X,Y,'o'); hold on plot(XPrediction,YPrediction,'LineWidth',2); % Optional error plot % plot(x,Y+DELTA,'g'); % plot(x,Y-DELTA,'g'); title(strcat(OrderLabel(BestOrder),... ' polynomial fitted to the training data')); xlabel('x'); ylabel('y');