function [slm,xp,yp] = slmengine(x,y,varargin)
% slmengine: estimates a spline function from data plus a fit prescription
% usage 1: slm = slmengine(x,y);
% usage 2: slm = slmengine(x,y,prescription);
% usage 3: slm = slmengine(x,y,prop1,val1,prop2,val2,...);
% usage 4: slm = slmengine(x,y,prescription,prop1,val1,prop2,val2,...);
% usage 5: [slm,xp,yp] = slmengine(x,y,prescription,prop1,val1,prop2,val2,...);
%
% Note: slmengine is the command driven tool for fitting a model
% to data. It uses either a prescription structure (as supplied by
% slmset) or sets of property/value pairs. Those pairs are defined
% in the help to slmset. slmengine is also used by slmfit, the gui
% tool for fitting a curve to data.
%
% Note: The optimization toolbox (lsqlin) is generally required for
% most fits. Some models (those where the break points are also
% estimated) require fmincon.
%
% Note: Some prescriptive parameters or combinations will be
% inappropriate for the lower degree models.
%
% arguments: (input)
% x,y - vectors of data used to fit the model. Any
% NaN or inf elements (in x or y) will be excluded
% from the fit.
%
% prescription - structure generated by slmset to control the
% fitting process.
%
% Alternatively, one can simply supply property/value
% pairs directly to slmengine.
%
% Arguments: (output)
% slm - model as a shape-language-model structure,
% normally constructed by slmfit or slmengine. slm
% may also be returned in other forms, either a pp
% structure tht ppval can use, or as a simple array
% of coefficients.
%
% Upon return, slm will contain a strucure field 'stats'
% that holds fields which describe the quality of fit:
%
% slm.stats.TotalDoF = Total degrees of freedom in the model
%
% slm.stats.NetDoF = Net degrees of freedom in the spline model.
% Thus, NetDoF reflects any equality constraints in the model
%
% slm.stats.R2 = Traditional R-squared coefficient
%
% slm.stats.R2Adj = Adjusted R^2, accounting for changing
% degrees of freedom. For a definition, look here:
%
% http://en.wikipedia.org/wiki/Coefficient_of_determination
%
% slm.stats.RMSE = Root-Mean-Squared-Error
%
% slm.stats.ErrorRange = 1x2 vector that contains the minimmum
% and maximum errors, defined as (Yhat - Y)
%
% slm.stats.Quartiles = 1x2 vector that contains the 25% and 75%
% error quartiles. This gives a trimmed measure of the
% error magnitudes, less any outliers.
%
% xp, yp - predicted points along the fitted curve. The number of these
% points is defined by the predictions property. (see slmset)
% xp will be a list of equally spaced points.
%
% Example:
% (See the SLM_tutorial demo for other examples.)
% x = rand(1,50);
% y = sin(x*2*pi) + randn(size(x))/10;
%
% slm = slmengine(x,y,'knots',0:.2:1,'plot','on','concavedown','on','minvalue',0)
% slm =
% form: 'slm'
% degree: 3
% knots: [6x1 double]
% coef: [6x2 double]
% prescription: [1x1 struct]
% x: [50x1 double]
% y: [50x1 double]
%
% Returned in a form that ppval or fnval can use:
%
% slm = slmengine(x,y,'knots',0:.2:1,'plot','on','concavedown','on','minvalue',0,'res','pp')
% slm =
% form: 'pp'
% breaks: [0 0.2 0.4 0.6 0.8 1]
% coefs: [5x4 double]
% pieces: 5
% order: 4
% dim: 1
% prescription: [1x1 struct]
%
%
% See also: spap2, ppval, fnval
%
% Author: John D'Errico
% e-mail: woodchips@rochester.rr.com
% Release: 1.0
% Release date: 1/2/08
% which form does the shape prescription take? Or do we
% just use our defaults for the prescription?
if nargin<2
error('SLMENGINE:improperdata','Must supply two vectors to fit a model: x & y')
elseif nargin==2
prescription = slmset;
elseif (nargin>=3)
% prescription supplied directly or as property/value pairs,
% or as a prescription, modified by a few set of pairs.
% slmset will resolve all of these cases.
prescription = slmset(varargin{:});
end
% check the data for size, turning it into column vectors
x = x(:);
y = y(:);
n = length(x);
if n~=length(y)
error('SLMENGINE:inconsistentdata','x and y must be the same size')
end
% were there any NaN or inf elements in the data?
k = isnan(x) | isnan(y) | isinf(x) | isinf(y);
if any(k)
% drop them from the analysis
x(k) = [];
y(k) = [];
% also drop corresponding weights if they were supplied
if ~isempty(prescription.Weights)
prescription.Weights(k) = [];
end
n = length(x);
end
% if weights or errorbars were set, verify the sizes of these
% parameters, compared to the number of data points.
if ~isempty(prescription.Weights)
prescription.Weights = prescription.Weights(:);
if n~=length(prescription.Weights)
error('SLMENGINE:inconsistentweights','Weights vector must be the same length as # of data points')
end
end
if ~isempty(prescription.ErrorBar)
EB = prescription.ErrorBar;
if length(EB) == 1
prescription.ErrorBar = repmat(EB,n,2);
elseif (size(EB,1) == 1)
prescription.ErrorBar = repmat(EB',1,2);
elseif (size(EB,2) == 1)
prescription.ErrorBar = repmat(EB,1,2);
end
end
if prescription.Verbosity > 1
disp('=========================================')
disp('Model Prescription')
disp(prescription)
end
% is this a free knot problem?
if strcmp(prescription.InteriorKnots,'free')
% the alternative is normally 'fixed', but its free now
% get the starting values for the knots
% knots vector, dx
if length(prescription.Knots)==1
[knots,nk] = chooseknots(prescription.Knots,n,x);
else
% we should check that the knots contain the data
knots = sort(prescription.Knots(:));
nk= length(knots);
if (knots(1)>min(x)) || (knots(end)<max(x))
error('SLMENGINE:inadequateknots',['Knots do not contain the data. Data range: ',num2str([min(x),max(x)])])
end
end
% there must be at least 3 knots
if nk<3
error('SLMENGINE:inadequateknots','Free knot estimation requires at least 3 knots')
end
% set the constraints for fmincon
% 1. no two knots may lie too close to each other
tol = min(0.001,0.1/(nk-1));
mindelta = tol*(knots(end) - knots(1));
A = zeros(2,nk-2);
A(1,1) = -1;
A(2,nk-2) = 1;
b = [-knots(1) - mindelta ; knots(end) - mindelta];
A = [A;full(spdiags(repmat([1 -1],nk-3,1),[0 1],nk-3,nk-2))];
b = [b;repmat(-mindelta,nk-3,1)];
% set up the optimization parameters (fmincon)
fminconoptions = optimset('fmincon');
fminconoptions.LargeScale = 'off';
fminconoptions.Algorithm = 'active-set';
if prescription.Verbosity>1
fminconoptions.Display = 'iter';
else
fminconoptions.Display = 'off';
end
% change the prescription for the subsequent internal calls
prescrip = prescription;
prescrip.Knots = knots;
prescrip.InteriorKnots = 'fixed';
prescrip.Plot = 'off';
% call the optimizer
intknots = knots(2:(end-1));
intknots = fmincon(@free_knot_obj,intknots,A,b, ...
[],[],[],[],[],fminconoptions,x,y,prescrip);
% fit the curve one last time with the final knots
prescrip.Knots(2:(end-1)) = intknots;
slm = slmengine(x,y,prescrip);
else
% the knots are fixed, just estimate the model
% do the appropriate fit. break them down into special
% cases, purely for simplicity of the code.
switch prescription.Degree
case {0 'constant'}
slm = slmengine_constant(x,y,prescription);
case {1 'linear'}
slm = slmengine_linear(x,y,prescription);
case {3 'cubic'}
slm = slmengine_cubic(x,y,prescription);
end
end
% do we need to report anything? I.e., what was the
% prescription.Verbosity setting?
if prescription.Verbosity>0
disp('=========================================')
disp('MODEL STATISTICS REPORT')
disp(['Number of data points: ',num2str(n)])
disp(['Total degrees of freedom: ',num2str(slm.stats.TotalDoF)])
disp(['Net degrees of freedom: ',num2str(slm.stats.NetDoF)])
disp(['R-squared: ',num2str(slm.stats.R2)])
disp(['Adjusted R-squared: ',num2str(slm.stats.R2Adj)])
disp(['RMSE: ',num2str(slm.stats.RMSE)])
disp(['Range of prediction errors: ',num2str(slm.stats.ErrorRange)])
disp(['Error quartiles (25%, 75%): ',num2str(slm.stats.Quartiles)])
disp('=========================================')
end
% append the shape prescription to the final structure
% this can be used as a template for fitting other functions,
% but also as documentation for the curve fit. This is why
% it is attached to the model structure, rather than as a
% second return argument.
slm.prescription = prescription;
% also attach the data used to build the model, purely
% for documentation purposes.
slm.x = x;
slm.y = y;
% do we need to plot the curve? I.e., what was the
% prescription.Plot setting?
if strcmp(prescription.Plot,'on')
% plot the curve
plotslm(slm)
end
% do we generate points on the curve?
if isempty(prescription.Predictions)
% none were asked for
xp = [];
yp = [];
else
% generate the required points on the curve
xp = linspace(min(x),max(x),prescription.Predictions);
yp = slmeval(xp,slm,0);
end
% do they want the result in a 'pp' form?
if strcmpi(prescription.Result,'pp') && ~strcmpi(slm.form,'pp')
% convert to pp form
stats = slm.stats;
slm = slm2pp(slm);
% make sure the stats field gets copied into the pp form
slm.stats = stats;
end
% ========================================================
% =========== slmengines for each model degree ============
% ========================================================
% ========================================================
% ============== piecewise constant model ================
% ========================================================
function slm = slmengine_constant(x,y,prescription)
% fits a piecewise constant shape prescriptive model
% check for inappropriate properties for a piecewise
% constant model
property_check(prescription,'constant')
% simple things about data first...
% slmengine has already made it a column vector,
% and ensured compatibility in length between
% x and y
nx = length(x);
% knots vector, dx
if length(prescription.Knots)==1
[knots,nk] = chooseknots(prescription.Knots,nx,x);
else
% we should check that the knots contain the data
knots = sort(prescription.Knots(:));
nk= length(knots);
if (knots(1)>min(x)) || (knots(end)<max(x))
error('SLMENGINE:inadequateknots',['Knots do not contain the data. Data range: ',num2str([min(x),max(x)])])
end
end
dx = diff(knots);
if any(dx==0)
error('SLMENGINE:indistinctknots','Knots must be distinct.')
end
% number of coefficients to estimate
% a piecewise constant function has one coefficient
% at each knot, but the last knot is "not".
nc = nk-1;
% create empty design, equality, inequality
% constraints arrays. also rhs vectors for each
Mdes = zeros(0,nc);
Meq = Mdes;
Mineq = Mdes;
Mreg = Mdes;
rhseq = [];
rhsineq = [];
rhsreg = [];
% -------------------------------------
% build design matrix -
% first, bin the data - histc wll do it
[junk,xbin] = histc(x,knots); %#ok
% any point which falls at the top end, is said to
% be in the last bin.
xbin(xbin==nk)=nk-1;
% the design matrix is easy to build - one call to sparse
Mdes = accumarray([(1:nx)',xbin],1,[nx,nc]);
rhs = y;
% have we been given a constraint on the sum of the residuals?
if ~isempty(prescription.SumResiduals)
Meq = [Meq;sum(Mdes,1)];
rhseq = [rhseq;prescription.SumResiduals + sum(rhs)];
end
% apply weights
W = prescription.Weights;
if ~isempty(W)
W = W(:);
if length(W)~=nx
error('SLMENGINE:inconsistentweights','Weight vector is not the same length as data')
end
% scale weight vector
W = sqrt(nx)*W./norm(W);
W = spdiags(W,0,nx,nx);
Mdes = W*Mdes;
rhs = W*rhs;
end
% -------------------------------------
% end conditions do not apply for a piecewise
% constant function
% -------------------------------------
% build regularizer
if nc>2
% use a simple second order finite difference
dx1 = dx(1:(nc-2));
dx2 = dx(2:(nc-1));
fda = [-2./(dx1.*(dx1+dx2)), 2./(dx1.*dx2), ...
-2./(dx2.*(dx1+dx2))];
Mreg = zeros(nc-2,nc);
for i=1:(nc-2);
Mreg(i,i+[0 1 2]) = fda(i,:);
end
rhsreg = zeros(nc-2,1);
end
% scale the regularizer before we apply the
% regularization parameter.
Mreg = Mreg/norm(Mreg,1);
% -------------------------------------
% single point equality constraints
% left hand side
if ~isempty(prescription.LeftValue)
M = zeros(1,nc);
M(1) = 1;
Meq = [Meq;M];
rhseq = [rhseq;prescription.LeftValue];
end
% right hand side
if ~isempty(prescription.RightValue)
M = zeros(1,nc);
M(end) = 1;
Meq = [Meq;M];
rhseq = [rhseq;prescription.RightValue];
end
% force curve through an x-y pair
xy = prescription.XY;
if ~isempty(xy)
n = size(xy,1);
if any(xy(:,1)<knots(1)) || any(xy(:,1)>knots(end))
error('SLMENGINE:improperconstraint','XY pairs to force the curve through must lie inside the knots')
end
[junk,ind] = histc(xy(:,1),knots); %#ok
ind(ind==(nc+1))=nc;
M = sparse((1:n)',ind,1,n,nc);
Meq = [Meq;M];
rhseq = [rhseq;xy(:,2)];
end
% -------------------------------------
% Integral equality constraint
if ~isempty(prescription.Integral)
% Rectangle rule. The last knot point
% adds no contribution to the area, but
% there are only nk-1 knots anyway.
M = dx(:)';
Meq = [Meq;M];
rhseq = [rhseq;prescription.Integral];
end
% -------------------------------------
% single point inequality constraints
% left hand side minimum
if ~isempty(prescription.LeftMinValue)
M = zeros(1,nc);
M(1) = -1;
Mineq = [Mineq;M];
rhsineq = [rhsineq;-prescription.LeftMinValue];
end
% left hand side maximum
if ~isempty(prescription.LeftMaxValue)
M = zeros(1,nc);
M(1) = 1;
Mineq = [Mineq;M];
rhsineq = [rhsineq;prescription.LeftMinValue];
end
% right hand side min
if ~isempty(prescription.RightMinValue)
M = zeros(1,nc);
M(end) = -1;
Mineq = [Mineq;M];
rhsineq = [rhsineq;-prescription.RightMinValue];
end
% right hand side max
if ~isempty(prescription.RightMaxValue)
M = zeros(1,nc);
M(end) = 1;
Mineq = [Mineq;M];
rhsineq = [rhsineq;prescription.RightMaxValue];
end
% -------------------------------------
% Error bar inequality constraints
if ~isempty(prescription.ErrorBar)
% lower bounds
Mineq = [Mineq;-Mdes];
rhsineq = [rhsineq;-(y-prescription.ErrorBar(:,1))];
% upper bounds
Mineq = [Mineq;Mdes];
rhsineq = [rhsineq;y+prescription.ErrorBar(:,2)];
end
% -------------------------------------
% constant region(s)
if ~isempty(prescription.ConstantRegion)
% there was at least one constant regions specified
cr = prescription.ConstantRegion;
M = zeros(0,nc);
n = 0;
for i=1:size(cr,1)
% enforce constancy between the given range limits
for j=1:(nc-1)
if (knots(j+1)>=cr(i,1)) && (knots(j+1)<cr(i,2))
n=n+1;
M(n,j+[0 1]) = [-1 1];
end
end
end
Meq = [Meq;M];
rhseq = [rhseq;zeros(size(M,1),1)];
end
% -------------------------------------
% overall inequalities
% global min value
if ~isempty(prescription.MinValue)
M = -eye(nc,nc);
Mineq = [Mineq;M];
rhsineq = [rhsineq;repmat(-prescription.MinValue,nc,1)];
end
% global max value
if ~isempty(prescription.MaxValue)
M = eye(nc,nc);
Mineq = [Mineq;M];
rhsineq = [rhsineq;repmat(prescription.MaxValue,nc,1)];
end
% SimplePeak and SimpleValley are really just composite
% properties. A peak at x == a is equivalent to a monotone
% increasing function for x<=a, and a monotone decreasing
% function for x>=a. Likewise a valley is just the opposite.
%
% So specifying a peak or a valley will cause any monotonicity
% specs to be ignored. slmset has already checked for this,
% and turned off any other monotonicity based properties.
sp = prescription.SimplePeak;
if isnumeric(sp) && ~isempty(sp)
prescription.Increasing = [knots(1),sp];
prescription.Decreasing = [sp,knots(end)];
end
sv = prescription.SimpleValley;
if isnumeric(sv) && ~isempty(sv)
prescription.Decreasing = [knots(1), sv];
prescription.Increasing = [sv,knots(end)];
end
% monotonicity?
% increasing regions
incR = prescription.Increasing;
L=0;
if ischar(incR)
if strcmp(incR,'on')
L=L+1;
mono(L).knotlist = 1:nc;
mono(L).direction = 1;
mono(L).range = [];
end
elseif ~isempty(incR)
for i=1:size(incR,1)
L=L+1;
mono(L).knotlist = []; %#ok
mono(L).direction = 1; %#ok
mono(L).range = sort(incR(i,:)); %#ok
end
end
% decreasing regions
decR = prescription.Decreasing;
if ischar(decR)
if strcmp(decR,'on')
L=L+1;
mono(L).knotlist = 1:nc;
mono(L).direction = -1;
mono(L).range = [];
end
elseif ~isempty(decR)
for i=1:size(decR,1)
L=L+1;
mono(L).knotlist = [];
mono(L).direction = -1;
mono(L).range = sort(decR(i,:));
end
end
if L>0
% there were at least some monotone regions specified
M = zeros(0,nc);
n = 0;
for i=1:L
if isempty(mono(L).range)
% the entire range was specified to be monotone
for j=1:(nc-1)
n=n+1;
M(n,j+[0 1]) = [1 -1]*mono(i).direction;
end
else
% only enforce monotonicity between the given range limits
for j=1:(nc-1)
if (knots(j)<mono(i).range(2)) && ...
(knots(j+1)>=mono(i).range(1))
n=n+1;
M(n,j+[0 1]) = [1 -1]*mono(i).direction;
end
end
end
end
Mineq = [Mineq;M];
rhsineq = [rhsineq;zeros(size(M,1),1)];
end
% -------------------------------------
% Use envelope inequalities?
switch prescription.Envelope
case 'off'
% nothing to do in this case
case 'supremum'
% yhat - y >= 0
Mineq = [Mineq;-Mdes];
rhsineq = [rhsineq;-rhs];
case 'infimum'
% yhat - y <= 0
Mineq = [Mineq;Mdes];
rhsineq = [rhsineq;rhs];
end
% -------------------------------------
% scale equalities for unit absolute row sum
if ~isempty(Meq)
rs = diag(1./sum(abs(Meq),2));
Meq = rs*Meq;
rhseq = rs*rhseq;
end
% scale inequalities for unit absolute row sum
if ~isempty(Mineq)
rs = diag(1./sum(abs(Mineq),2));
Mineq = rs*Mineq;
rhsineq = rs*rhsineq;
end
% -------------------------------------
% now worry about the regularization. There are three
% possible cases.
% 1. We have a given regularization parameter
% 2. We have a given rmse that we wish to match
% 3. We must use cross validation to choose the parameter
RP = prescription.Regularization;
if (isnumeric(RP) && (RP>=0)) || ((ischar(RP)) && (strcmpi(RP,'smoothest')))
% solve the problem using the given regularization parameter
coef = solve_slm_system(RP,Mdes,rhs,Mreg,rhsreg, ...
Meq,rhseq,Mineq,rhsineq,prescription);
elseif isnumeric(RP) && (RP<0)
% we must match abs(RP) as the rmse.
aim_rmse = abs(RP);
fminbndoptions = optimset('fminbnd');
RP = fminbnd(@match_rmse,-6,6,fminbndoptions, ...
Mdes,rhs,Mreg,rhsreg,Meq,rhseq,Mineq,rhsineq, ...
aim_rmse,prescription);
% we logged the parameter in the optimization. undo that for
% the final call
RP = 10^RP;
% do one final call to get the final coefficients
coef = solve_slm_system(RP,Mdes,rhs,Mreg,rhsreg, ...
Meq,rhseq,Mineq,rhsineq,prescription);
elseif ischar(RP)
% its a cross validation problem to solve. drop out each point
% in turn from the model, then use fminbnd to minimize the
% predicted sum of squares at the missing points.
fminbndoptions = optimset('fminbnd');
% fminbndoptions.Display = 'iter';
RP = fminbnd(@min_cv,-6,6,fminbndoptions, ...
Mdes,rhs,Mreg,rhsreg,Meq,rhseq,Mineq,rhsineq,prescription);
% we logged the parameter in the optimization. undo that for
% the final call
RP = 10^RP;
% do one final call to get the final coefficients
coef = solve_slm_system(RP,Mdes,rhs,Mreg,rhsreg, ...
Meq,rhseq,Mineq,rhsineq,prescription);
end
% -------------------------------------
% unpack coefficients into the result structure
slm.form = 'slm';
slm.degree = 0;
slm.knots = knots;
slm.coef = coef;
% degrees of freedom available
slmstats.TotalDoF = nk - 1;
slmstats.NetDoF = slmstats.TotalDoF - size(Meq,1);
% this function does all of the stats, stuffing into slmstats
slm.stats = modelstatistics(slmstats,Mdes,y,coef);
% ========================================================
% =============== piecewise linear model =================
% ========================================================
function slm = slmengine_linear(x,y,prescription)
% fits a piecewise linear shape prescriptive model
% check for inappropriate properties for a piecewise
% linear model
property_check(prescription,'linear')
% simple things about data first...
% slmengine has already made it a column vector,
% and ensured compatibility in length between
% x and y
nx = length(x);
% knots vector, dx
if length(prescription.Knots)==1
% just given a number of knots
[knots,nk] = chooseknots(prescription.Knots,nx,x);
else
% we should check that the knots contain the data
knots = sort(prescription.Knots(:));
nk= length(knots);
if (knots(1)>min(x)) || (knots(end)<max(x))
error('SLMENGINE:inadequateknots',['Knots do not contain the data. Data range: ',num2str([min(x),max(x)])])
end
end
dx = diff(knots);
if any(dx==0)
error('SLMENGINE:indistinctknots','Knots must be distinct.')
end
% number of coefficients to estimate.
% a piecewise linear Hermite has one coefficient
% at each knot.
nc = nk;
% create empty design, equality, inequality
% constraints arrays. also rhs vectors for each
Mdes = zeros(0,nc);
Meq = Mdes;
Mineq = Mdes;
Mreg = Mdes;
rhseq = [];
rhsineq = [];
rhsreg = [];
% -------------------------------------
% build design matrix -
% first, bin the data - histc wll do it
[junk,xbin] = histc(x,knots); %#ok
% any point which falls at the top end, is said to
% be in the last bin.
xbin(xbin==nk)=nk-1;
% build design matrix
t = (x - knots(xbin))./dx(xbin);
Mdes = accumarray([repmat((1:nx)',2,1), ...
[xbin;xbin+1]],[1-t;t],[nx,nc]);
rhs = y;
% have we been given a constraint on the sum of the residuals?
if ~isempty(prescription.SumResiduals)
Meq = [Meq;sum(Mdes,1)];
rhseq = [rhseq;prescription.SumResiduals + sum(rhs)];
end
% apply weights
W = prescription.Weights;
if ~isempty(W)
W = W(:);
if length(W)~=nx
error('SLMENGINE:inconsistentweights','Weight vector is not the same length as data')
end
% scale weight vector
W = sqrt(nx)*W./norm(W);
W = spdiags(W,0,nx,nx);
Mdes = W*Mdes;
rhs = W*rhs;
end
% -------------------------------------
% build regularizer, only bother if at least 3 knots
if nc>2
% use a simple second order finite difference
dx1 = dx(1:(nc-2));
dx2 = dx(2:(nc-1));
fda = [-2./(dx1.*(dx1+dx2)), 2./(dx1.*dx2), ...
-2./(dx2.*(dx1+dx2))];
Mreg = zeros(nc-2,nc);
for i=1:(nc-2);
Mreg(i,i+[0 1 2]) = fda(i,:);
end
rhsreg = zeros(nc-2,1);
end
% scale the regularizer before we apply the
% regularization parameter.
Mreg = Mreg/norm(Mreg,1);
% -------------------------------------
% end conditions only apply for a piecewise linear
% function, IF the conditions are periodicity
if strcmp(prescription.EndConditions,'periodic')
% set the first and last function values equal
M = zeros(1,nc);
M([1,end]) = [-1 1];
Meq = [Meq;M];
rhseq = [rhseq;0];
end
% -------------------------------------
% single point equality constraints
% left hand side
if ~isempty(prescription.LeftValue)
M = zeros(1,nc);
M(1) = 1;
Meq = [Meq;M];
rhseq = [rhseq;prescription.LeftValue];
end
% right hand side
if ~isempty(prescription.RightValue)
M = zeros(1,nc);
M(end) = 1;
Meq = [Meq;M];
rhseq = [rhseq;prescription.RightValue];
end
% left end slope
if ~isempty(prescription.LeftSlope)
M = zeros(1,nc);
M(1:2) = [-1 1]/dx(1);
Meq = [Meq;M];
rhseq = [rhseq;prescription.LeftSlope];
end
% Right end slope
if ~isempty(prescription.RightSlope)
M = zeros(1,nc);
M([nc-1, nc]) = [-1 1]/dx(end);
Meq = [Meq;M];
rhseq = [rhseq;prescription.RightSlope];
end
% force curve through an x-y pair
xy = prescription.XY;
if ~isempty(xy)
n = size(xy,1);
if any(xy(:,1)<knots(1)) || any(xy(:,1)>knots(end))
error('SLMENGINE:improperconstraint','XY pairs to force the curve through must lie inside the knots')
end
[junk,ind] = histc(xy(:,1),knots); %#ok
ind(ind==(nc+1))=nc;
t = (xy(:,1) - knots(ind))./dx(ind);
M = sparse(repmat((1:n)',1,2),[ind,ind+1],[1-t,t],n,nc);
Meq = [Meq;M];
rhseq = [rhseq;xy(:,2)];
end
% force slope at a point, or set of points
xyp = prescription.XYP;
if ~isempty(xyp)
n = size(xyp,1);
if any(xyp(:,1)<knots(1)) || any(xyp(:,1)>knots(end))
error('SLMENGINE:improperconstraint','X-Y'' pairs to enforce the slope at must lie inside the knots')
end
[junk,ind] = histc(xyp(:,1),knots); %#ok
ind(ind==(nc+1))=nc;
M = sparse(repmat((1:n)',1,2),[ind,ind+1],(1./dx(ind))*[-1,1],n,nc);
Meq = [Meq;M];
rhseq = [rhseq;xyp(:,2)];
end
% -------------------------------------
% constant region(s)
if ~isempty(prescription.ConstantRegion)
% there was at least one constant region specified
cr = prescription.ConstantRegion;
M = zeros(0,nc);
n = 0;
for i=1:size(cr,1)
% enforce constancy between the given range limits
for j=1:(nc-1)
if (knots(j+1)>=cr(i,1)) && (knots(j)<cr(i,2))
n=n+1;
M(n,j+[0 1]) = [-1 1];
end
end
end
Meq = [Meq;M];
rhseq = [rhseq;zeros(size(M,1),1)];
end
% -------------------------------------
% linear region(s)
% a linear region imply means that across any knot
% inside that region, the slopes must be constant
if ~isempty(prescription.LinearRegion) && (nc>2)
% there was at least one linear region specified
lr = prescription.LinearRegion;
M = zeros(0,nc);
n = 0;
for i=1:size(lr,1)
% enforce constancy of slope between the given range limits
for j=1:(nc-2)
if (knots(j+1)>=lr(i,1)) && (knots(j+1)<lr(i,2))
n=n+1;
M(n,j+[0 1 2]) = [-1 1 0]/dx(j) - [0 -1 1]/dx(j+1);
end
end
end
Meq = [Meq;M];
rhseq = [rhseq;zeros(size(M,1),1)];
end
% -------------------------------------
% Integral equality constraint
if ~isempty(prescription.Integral)
% Trapezoidal rule over the support.
M = ([dx(:)',0] + [0,dx(:)'])/2;
Meq = [Meq;M];
rhseq = [rhseq;prescription.Integral];
end
% -------------------------------------
% single point inequality constraints
% left hand side minimum
if ~isempty(prescription.LeftMinValue)
M = zeros(1,nc);
M(1) = -1;
Mineq = [Mineq;M];
rhsineq = [rhsineq;-prescription.LeftMinValue];
end
% left hand side maximum
if ~isempty(prescription.LeftMaxValue)
M = zeros(1,nc);
M(1) = 1;
Mineq = [Mineq;M];
rhsineq = [rhsineq;prescription.LeftMinValue];
end
% right hand side min
if ~isempty(prescription.RightMinValue)
M = zeros(1,nc);
M(end) = -1;
Mineq = [Mineq;M];
rhsineq = [rhsineq;-prescription.RightMinValue];
end
% right hand side max
if ~isempty(prescription.RightMaxValue)
M = zeros(1,nc);
M(end) = 1;
Mineq = [Mineq;M];
rhsineq = [rhsineq;prescription.RightMaxValue];
end
% left end slope min
if ~isempty(prescription.LeftMinSlope)
M = zeros(1,nc);
M(1:2) = [1 -1]/dx(1);
Mineq = [Mineq;M];
rhsineq = [rhsineq;-prescription.LeftMinSlope];
end
% left end slope max
if ~isempty(prescription.LeftMaxSlope)
M = zeros(1,nc);
M(1:2) = [-1 1]/dx(1);
Mineq = [Mineq;M];
rhsineq = [rhsineq;prescription.LeftMaxSlope];
end
% right end slope min
if ~isempty(prescription.RightMinSlope)
M = zeros(1,nc);
M([nc-1,nc]) = [1 -1]/dx(end);
Mineq = [Mineq;M];
rhsineq = [rhsineq;-prescription.RightMinSlope];
end
% right end slope max
if ~isempty(prescription.RightMaxSlope)
M = zeros(1,nc);
M([nc-1,nc]) = [-1 1]/dx(end);
Mineq = [Mineq;M];
rhsineq = [rhsineq;prescription.RightMaxSlope];
end
% -------------------------------------
% overall inequalities
% global min value
if ~isempty(prescription.MinValue)
M = -eye(nc,nc);
Mineq = [Mineq;M];
rhsineq = [rhsineq;repmat(-prescription.MinValue,nc,1)];
end
% global max value
if ~isempty(prescription.MaxValue)
M = eye(nc,nc);
Mineq = [Mineq;M];
rhsineq = [rhsineq;repmat(prescription.MaxValue,nc,1)];
end
% global min slope
if ~isempty(prescription.MinSlope)
M = full(spdiags((1./dx)*[1 -1],[0 1],nc-1,nc));
Mineq = [Mineq;M];
rhsineq = [rhsineq;repmat(-prescription.MinSlope,nc-1,1)];
end
% global max slope
if ~isempty(prescription.MaxSlope)
M = full(spdiags((1./dx)*[-1 1],[0 1],nc-1,nc));
Mineq = [Mineq;M];
rhsineq = [rhsineq;repmat(prescription.MaxSlope,nc-1,1)];
end
% -------------------------------------
% SimplePeak and SimpleValley are really just composite
% properties. A peak at x == a is equivalent to a monotone
% increasing function for x<=a, and a monotone decreasing
% function for x>=a. Likewise a valley is just the opposite.
%
% So specifying a peak or a valley will cause any monotonicity
% specs to be ignored. slmset has already checked for this,
% and turned off any other monotonicity based properties.
sp = prescription.SimplePeak;
if isnumeric(sp) && ~isempty(sp)
prescription.Increasing = [knots(1),sp];
prescription.Decreasing = [sp,knots(end)];
end
sv = prescription.SimpleValley;
if isnumeric(sv) && ~isempty(sv)
prescription.Decreasing = [knots(1), sv];
prescription.Increasing = [sv,knots(end)];
end
% -------------------------------------
% monotonicity?
% increasing regions
incR = prescription.Increasing;
L=0;
if ischar(incR)
if strcmp(incR,'on')
L=L+1;
mono(L).knotlist = 'all';
mono(L).direction = 1;
mono(L).range = [];
end
elseif ~isempty(incR)
for i=1:size(incR,1)
L=L+1;
mono(L).knotlist = []; %#ok
mono(L).direction = 1; %#ok
mono(L).range = sort(incR(i,:)); %#ok
end
end
% decreasing regions
decR = prescription.Decreasing;
if ischar(decR)
if strcmp(decR,'on')
L=L+1;
mono(L).knotlist = 'all';
mono(L).direction = -1;
mono(L).range = [];
end
elseif ~isempty(decR)
for i=1:size(decR,1)
L=L+1;
mono(L).knotlist = [];
mono(L).direction = -1;
mono(L).range = sort(decR(i,:));
end
end
if L>0
% there were at least some monotone regions specified
M = zeros(0,nc);
n = 0;
for i=1:L
if isempty(mono(L).range)
% the entire range was specified to be monotone
for j=1:(nc-1)
n=n+1;
M(n,j+[0 1]) = [1 -1]*mono(i).direction;
end
else
% only enforce monotonicity between the given range limits
for j=1:(nc-1)
if (knots(j)<mono(i).range(2)) && ...
(knots(j+1)>=mono(i).range(1))
n=n+1;
M(n,j+[0 1]) = [1 -1]*mono(i).direction;
end
end
end
end
Mineq = [Mineq;M];
rhsineq = [rhsineq;zeros(size(M,1),1)];
end
% -------------------------------------
% PositiveInflection and NegativeInflection are really just
% composite properties. A point of inflection at x == a is
% equivalent to a change in the sign of the curvature at
% x == a.
pinf = prescription.PositiveInflection;
if isnumeric(pinf) && ~isempty(pinf)
prescription.ConcaveDown = [knots(1),pinf];
prescription.ConcaveUp = [pinf,knots(end)];
end
ninf = prescription.NegativeInflection;
if isnumeric(ninf) && ~isempty(ninf)
prescription.ConcaveUp = [knots(1),ninf];
prescription.ConcaveDown = [ninf,knots(end)];
end
% -------------------------------------
% concave up or down?
% regions of postive curvature
cuR = prescription.ConcaveUp;
L=0;
if ischar(cuR)
if strcmp(cuR,'on')
L=L+1;
curv(L).knotlist = 'all';
curv(L).direction = 1;
curv(L).range = [];
end
elseif ~isempty(cuR)
for i=1:size(cuR,1)
L=L+1;
curv(L).knotlist = []; %#ok
curv(L).direction = 1; %#ok
curv(L).range = sort(cuR(i,:)); %#ok
end
end
% decreasing regions
cdR = prescription.ConcaveDown;
if ischar(cdR)
if strcmp(cdR,'on')
L=L+1;
curv(L).knotlist = 'all';
curv(L).direction = -1;
curv(L).range = [];
end
elseif ~isempty(cdR)
for i=1:size(cdR,1)
L=L+1;
curv(L).knotlist = [];
curv(L).direction = -1;
curv(L).range = sort(cdR(i,:));
end
end
if L>0
% there were at least some regions with specified curvature
M = zeros(0,nc);
n = 0;
for i=1:L
if isempty(curv(L).range)
% the entire range was specified to be
% curved in some direction
for j=2:(nc-1)
n=n+1;
M(n,j+[-1 0 1]) = ([-1 1 0]/dx(j-1) - [0 -1 1]/dx(j))*curv(i).direction;
end
else
% only enforce curvature between the given range limits
for j=2:(nc-1)
if (knots(j)<curv(i).range(2)) && ...
(knots(j)>=curv(i).range(1))
n=n+1;
M(n,j+[-1 0 1]) = ([-1 1 0]/dx(j-1) - [0 -1 1]/dx(j))*curv(i).direction;
end
end
end
end
Mineq = [Mineq;M];
rhsineq = [rhsineq;zeros(size(M,1),1)];
end
% -------------------------------------
% Error bar inequality constraints
if ~isempty(prescription.ErrorBar)
% lower bounds
Mineq = [Mineq;-Mdes];
rhsineq = [rhsineq;-(y-prescription.ErrorBar(:,1))];
% upper bounds
Mineq = [Mineq;Mdes];
rhsineq = [rhsineq;y+prescription.ErrorBar(:,2)];
end
% -------------------------------------
% Use envelope inequalities?
switch prescription.Envelope
case 'off'
% nothing to do in this case
case 'supremum'
% yhat - y >= 0
Mineq = [Mineq;-Mdes];
rhsineq = [rhsineq;-rhs];
case 'infimum'
% yhat - y <= 0
Mineq = [Mineq;Mdes];
rhsineq = [rhsineq;rhs];
end
% -------------------------------------
% scale equalities for unit absolute row sum
if ~isempty(Meq)
rs = diag(1./sum(abs(Meq),2));
Meq = rs*Meq;
rhseq = rs*rhseq;
end
% scale inequalities for unit absolute row sum
if ~isempty(Mineq)
rs = diag(1./sum(abs(Mineq),2));
Mineq = rs*Mineq;
rhsineq = rs*rhsineq;
end
% -------------------------------------
% now worry about the regularization. There are four
% possible cases.
% 1. We have a given regularization parameter
% 2. We have a given rmse that we wish to match
% 3. We must use cross validation to choose the parameter
RP = prescription.Regularization;
if (isnumeric(RP) && (RP>=0)) || ((ischar(RP)) && (strcmpi(RP,'smoothest')))
% solve the problem using the given regularization parameter
coef = solve_slm_system(RP,Mdes,rhs,Mreg,rhsreg, ...
Meq,rhseq,Mineq,rhsineq,prescription);
elseif isnumeric(RP) && (RP<0)
% we must match abs(RP) as the rmse.
aim_rmse = abs(RP);
fminbndoptions = optimset('fminbnd');
RP = fminbnd(@match_rmse,-6,6,fminbndoptions, ...
Mdes,rhs,Mreg,rhsreg,Meq,rhseq,Mineq,rhsineq, ...
aim_rmse,prescription);
% we logged the parameter in the optimization. undo that for
% the final call
RP = 10^RP;
% do one final call to get the final coefficients
coef = solve_slm_system(RP,Mdes,rhs,Mreg,rhsreg, ...
Meq,rhseq,Mineq,rhsineq,prescription);
elseif ischar(RP)
% its a cross validation problem to solve. drop out each point
% in turn from the model, then use fminbnd to minimize the
% predicted sum of squares at the missing points.
fminbndoptions = optimset('fminbnd');
% fminbndoptions.Display = 'iter';
RP = fminbnd(@min_cv,-6,6,fminbndoptions, ...
Mdes,rhs,Mreg,rhsreg,Meq,rhseq,Mineq,rhsineq,prescription);
% we logged the parameter in the optimization. undo that for
% the final call
RP = 10^RP;
% do one final call to get the final coefficients
coef = solve_slm_system(RP,Mdes,rhs,Mreg,rhsreg, ...
Meq,rhseq,Mineq,rhsineq,prescription);
end
% -------------------------------------
% unpack coefficients into the result structure
slm.form = 'slm';
slm.degree = 1;
slm.knots = knots;
slm.coef = coef;
% degrees of freedom available
slmstats.TotalDoF = nk;
slmstats.NetDoF = slmstats.TotalDoF - size(Meq,1);
% this function does all of the stats, stuffing into slmstats
slm.stats = modelstatistics(slmstats,Mdes,y,coef);
% ========================================================
% =============== piecewise cubic model ==================
% ========================================================
function slm = slmengine_cubic(x,y,prescription)
% fits a piecewise cubic shape prescriptive model
% check for inappropriate properties for a cubic model
property_check(prescription,'cubic')
% simple things about data first...
% slmengine has already made it a column vector,
% and ensured compatibility in length between
% x and y
nx = length(x);
% knots vector, dx
if length(prescription.Knots)==1
% just given a number of knots
[knots,nk] = chooseknots(prescription.Knots,nx,x);
else
% we should check that the knots contain the data
knots = sort(prescription.Knots(:));
nk= length(knots);
if (knots(1)>min(x)) || (knots(end)<max(x))
error('SLMENGINE:inadequateknots',['Knots do not contain the data. Data range: ',num2str([min(x),max(x)])])
end
end
dx = diff(knots);
if any(dx==0)
error('SLMENGINE:indistinctknots','Knots must be distinct.')
end
% number of coefficients to estimate.
% a piecewise cubic Hermite has two coefficients
% at each knot.
nc = 2*nk;
% create empty inequality constraints arrays. also
% a rhs vector
Mineq = zeros(0,nc);
rhsineq = [];
Meq = zeros(0,nc);
rhseq = [];
% -------------------------------------
% build design matrix -
% first, bin the data - histc wll do it best
[junk,xbin] = histc(x,knots); %#ok
% any point which falls at the top end, is said to
% be in the last bin.
xbin(xbin==nk)=nk-1;
% build design matrix
t = (x - knots(xbin))./dx(xbin);
t2 = t.^2;
t3 = t.^3;
s2 = (1-t).^2;
s3 = (1-t).^3;
vals = [3*s2-2*s3 ; 3*t2-2*t3 ; ...
-(s3-s2).*dx(xbin) ; (t3-t2).*dx(xbin)];
% the coefficients will be stored in two blocks,
% first nk function values, then nk derivatives.
Mdes = accumarray([repmat((1:nx)',4,1), ...
[xbin;xbin+1;nk+xbin;nk+xbin+1]],vals,[nx,nc]);
rhs = y;
% have we been given a constraint on the sum of the residuals?
if ~isempty(prescription.SumResiduals)
Meq = [Meq;sum(Mdes,1)];
rhseq = [rhseq;prescription.SumResiduals + sum(rhs)];
end
% apply weights
W = prescription.Weights;
if ~isempty(W)
W = W(:);
if length(W)~=nx
error('SLMENGINE:inconsistentweights','Weight vector is not the same length as data')
end
% scale weight vector
W = sqrt(nx)*W./norm(W);
W = spdiags(W,0,nx,nx);
Mdes = W*Mdes;
rhs = W*rhs;
end
% -------------------------------------
% Regularizer
% We are integrating the piecewise linear f''(x), as
% a quadratic form in terms of the (unknown) second
% derivatives at the knots.
Mreg=zeros(nk,nk);
Mreg(1,1:2)=[dx(1)/3 , dx(1)/6];
Mreg(nk,nk+[-1 0])=[dx(end)/6 , dx(end)/3];
for i=2:(nk-1)
Mreg(i,i+[-1 0 1])=[dx(i-1)/6 , (dx(i-1)+dx(i))/3 , dx(i)/6];
end
% do a matrix square root. cholesky is simplest & ok since regmat is
% positive definite. this way we can write the quadratic form as:
% s'*r'*r*s,
% where s is the vector of second derivatives at the knots.
Mreg=chol(Mreg);
% next, write the second derivatives as a function of the
% function values and first derivatives. s = [sf,sd]*[f;d]
sf = zeros(nk,nk);
sd = zeros(nk,nk);
for i = 1:(nk-1)
sf(i,i+[0 1]) = [-1 1].*(6/(dx(i)^2));
sd(i,i+[0 1]) = [-4 -2]./dx(i);
end
sf(nk,nk+[-1 0]) = [1 -1].*(6/(dx(end)^2));
sd(nk,nk+[-1 0]) = [2 4]./dx(end);
Mreg = Mreg*[sf,sd];
% scale the regularizer before we apply the
% regularization parameter.
Mreg = Mreg/norm(Mreg,1);
rhsreg = zeros(nk,1);
% -------------------------------------
% C2 continuity across knots
if strcmp(prescription.C2,'on')
MC2 = zeros(nk-2,nc);
for i = 1:(nk-2)
MC2(i,[i,i+1]) = [6 -6]./(dx(i).^2);
MC2(i,[i+1,i+2]) = MC2(i,[i+1,i+2]) + [6 -6]./(dx(i+1).^2);
MC2(i,nk+[i,i+1]) = [2 4]./dx(i);
MC2(i,nk+[i+1,i+2]) = MC2(i,nk+[i+1,i+2]) + [4 2]./dx(i+1);
end
Meq = [Meq;MC2];
rhseq = [rhseq;zeros(nk-2,1)];
end
% -------------------------------------
% end conditions (do nothing here if they are 'estimate')
if strcmp(prescription.EndConditions,'periodic')
% set the first and last function values equal
M = zeros(3,nc);
M(1,[1,nk]) = [-1 1];
% as well as the first derivatives
M(2,[nk+1,2*nk]) = [-1 1];
% and the second derivatives
M(3,[1 2]) = [-6 6]/(dx(1).^2);
M(3,[nk+1,nk+2]) = [-4 -2]/dx(1);
M(3,[nk-1,nk]) = [-6 6]/(dx(end).^2);
M(3,[2*nk-1,2*nk]) = [-2 -4]/dx(end);
Meq = [Meq;M];
rhseq = [rhseq;zeros(3,1)];
elseif strcmp(prescription.EndConditions,'not-a-knot') && nk==2
% third derivative is continuous across 2nd knot
M = zeros(1,nc);
M(1,[1 2]) = [12 -12]/(dx(1).^3);
M(1,[2 3]) = M(1,2:3) - [12 -12]/(dx(2).^3);
M(1,nk+[1 2]) = [6 6]/(dx(1).^2);
M(1,nk+[2 3]) = M(1,nk+[2 3]) - [6 6]/(dx(2).^2);
Meq = [Meq;M];
rhseq = [rhseq;0];
elseif strcmp(prescription.EndConditions,'notaknot') && nk>2
% third derivative is continuous at 2nd & penultimate knots
M = zeros(2,nc);
M(1,[1 2]) = [12 -12]/(dx(1).^3);
M(1,[2 3]) = M(1,2:3) - [12 -12]/(dx(2).^3);
M(1,nk+[1 2]) = [6 6]/(dx(1).^2);
M(1,nk+[2 3]) = M(1,nk+[2 3]) - [6 6]/(dx(2).^2);
M(2,nk + [-2 -1]) = [12 -12]/(dx(end - 1).^3);
M(2,nk + [-1 0]) = M(2,nk + [-1 0]) - [12 -12]/(dx(end).^3);
M(2,2*nk+[-2 -1]) = [6 6]/(dx(end-1).^2);
M(2,2*nk+[-1 0]) = M(2,2*nk+[-1 0]) - [6 6]/(dx(end).^2);
Meq = [Meq;M];
rhseq = [rhseq;0;0];
elseif strcmp(prescription.EndConditions,'natural')
% second derivative at each end is zero
M = zeros(2,nc);
M(1,[1 2]) = [-6 6]/(dx(1).^2);
M(1,nk+[1 2]) = [-4 -2]/dx(1);
M(2,nk+[-1 0]) = [6 -6]/(dx(end).^2);
M(2,2*nk+[-1 0]) = [2 4]/dx(end);
Meq = [Meq;M];
rhseq = [rhseq;0;0];
end
% -------------------------------------
% single point equality constraints at a knot
% left hand side
if ~isempty(prescription.LeftValue)
M = zeros(1,nc);
M(1,1) = 1;
Meq = [Meq;M];
rhseq = [rhseq;prescription.LeftValue];
end
% right hand side
if ~isempty(prescription.RightValue)
M = zeros(1,nc);
M(1,nk) = 1;
Meq = [Meq;M];
rhseq = [rhseq;prescription.RightValue];
end
% left end slope
if ~isempty(prescription.LeftSlope)
M = zeros(1,nc);
M(1,nk+1) = 1;
Meq = [Meq;M];
rhseq = [rhseq;prescription.LeftSlope];
end
% Right end slope
if ~isempty(prescription.RightSlope)
M = zeros(1,nc);
M(1,2*nk) = 1;
Meq = [Meq;M];
rhseq = [rhseq;prescription.RightSlope];
end
% -------------------------------------
% force curve through an x-y pair or pairs
xy = prescription.XY;
if ~isempty(xy)
n = size(xy,1);
if any(xy(:,1)<knots(1)) || any(xy(:,1)>knots(end))
error('SLMENGINE:improperconstraint','XY pairs to force the curve through must lie inside the knots')
end
[junk,ind] = histc(xy(:,1),knots); %#ok
ind(ind==(nk))=nk-1;
t = (xy(:,1) - knots(ind))./dx(ind);
t2 = t.^2;
t3 = t.^3;
s2 = (1-t).^2;
s3 = (1-t).^3;
vals = [3*s2-2*s3 ; 3*t2-2*t3 ; ...
-(s3-s2).*dx(ind) ; (t3-t2).*dx(ind)];
M = accumarray([repmat((1:n)',4,1), ...
[ind;ind+1;nk+ind;nk+ind+1]],vals,[n,nc]);
Meq = [Meq;M];
rhseq = [rhseq;xy(:,2)];
end
% -------------------------------------
% force slope at a point, or a set of points
xyp = prescription.XYP;
if ~isempty(xyp)
n = size(xyp,1);
if any(xyp(:,1)<knots(1)) || any(xyp(:,1)>knots(end))
error('SLMENGINE:improperconstraint','X-Y'' pairs to enforce the slope at must lie inside the knots')
end
[junk,ind] = histc(xyp(:,1),knots); %#ok
ind(ind==(nk))=nk-1;
t = (xyp(:,1) - knots(ind))./dx(ind);
t2 = t.^2;
s2 = (1-t).^2;
s = (1-t);
vals = [(-6*s+6*s2)./dx(ind) ; ...
(6*t-6*t2)./dx(ind) ; (3*s2-2*s) ; (3*t2-2*t)];
M = accumarray([repmat((1:n)',4,1), ...
[ind;ind+1;nk+ind;nk+ind+1]],vals,[n,nc]);
Meq = [Meq;M];
rhseq = [rhseq;xyp(:,2)];
end
% -------------------------------------
% force second derivative at a point, or set of points
xypp = prescription.XYPP;
if ~isempty(xypp)
n = size(xypp,1);
if any(xypp(:,1)<knots(1)) || any(xypp(:,1)>knots(end))
error('SLMENGINE:improperconstraint','X-Y'' pairs to enforce y" at must lie inside the knots')
end
[junk,ind] = histc(xypp(:,1),knots); %#ok
ind(ind==(nk))=nk-1;
t = (xypp(:,1) - knots(ind))./dx(ind);
s = (1-t);
vals = [(6 - 12*s)./dx(ind).^2; ...
(6 - 12*t)./dx(ind).^2 ; ...
-(6*s - 2)./dx(ind) ; ...
(6*t - 2)./dx(ind)];
M = accumarray([repmat((1:n)',4,1), ...
[ind;ind+1;nk+ind;nk+ind+1]],vals,[n,nc]);
Meq = [Meq;M];
rhseq = [rhseq;xypp(:,2)];
end
% -------------------------------------
% force third derivative at a point, or set of points
xyppp = prescription.XYPPP;
if ~isempty(xyppp)
n = size(xyppp,1);
if any(xyppp(:,1)<knots(1)) || any(xyppp(:,1)>knots(end))
error('SLMENGINE:improperconstraint','X-Y'''''' pairs to enforce y'''''' at must lie inside the knots')
end
[junk,ind] = histc(xyppp(:,1),knots); %#ok
ind(ind==(nk))=nk-1;
vals = [12./dx(ind).^3; -12./dx(ind).^3 ; ...
6./dx(ind).^2 ;6./dx(ind).^2];
M = accumarray([repmat((1:n)',4,1), ...
[ind;ind+1;nk+ind;nk+ind+1]],vals,[n,nc]);
Meq = [Meq;M];
rhseq = [rhseq;xyppp(:,2)];
end
% -------------------------------------
% constant region(s)
if ~isempty(prescription.ConstantRegion)
% there was at least one constant region specified
cr = prescription.ConstantRegion;
M = zeros(0,nc);
n = 0;
for i=1:size(cr,1)
% enforce constancy between the given range limits
flag = true;
for j=1:(nk-1)
if (knots(j+1)>=cr(i,1)) && (knots(j)<cr(i,2))
% f(j) == f(j+1)
n=n+1;
M(n,j+[0 1]) = [-1 1];
% also zero f'' on this interval
if flag
% only set this constraint for the first knot
% interval in a constant region
n=n+1;
M(n,j+[0 1 nk nk+1]) = [-6/(dx(j).^2) ; ...
6/(dx(j).^2) ; -4/dx(j) ; -2/dx(j)];
flag = false;
end
n=n+1;
M(n,j+[0 1 nk nk+1]) = [-6/(dx(j).^2) ; ...
6/(dx(j).^2) ; 2/dx(j) ; 4/dx(j)];
end
end
end
Meq = [Meq;M];
rhseq = [rhseq;zeros(size(M,1),1)];
end
% -------------------------------------
% linear region(s)
% a linear region simply means that across any knot
% inside that region, the slopes must be constant
if ~isempty(prescription.LinearRegion) && (nc>2)
% there was at least one linear region specified
lr = prescription.LinearRegion;
M = zeros(0,nc);
n = 0;
for i=1:size(lr,1)
% enforce constancy of slope between the given range limits
flag = true;
for j=1:(nk-1)
if (knots(j+1)>=lr(i,1)) && (knots(j)<lr(i,2))
% f'(j) == f'(j+1)
n=n+1;
M(n,j+nk+[0 1]) = [-1 1];
% also zero f'' at each end of this interval
if flag
% only set this constraint for the first knot
% interval in a linear region
n=n+1;
M(n,j+[0 1 nk nk+1]) = [-6/(dx(j).^2) ; ...
6/(dx(j).^2) ; -4/dx(j) ; -2/dx(j)];
flag = false;
end
n=n+1;
M(n,j+[0 1 nk nk+1]) = [6/(dx(j).^2) ; ...
-6/(dx(j).^2) ; 2/dx(j) ; 4/dx(j)];
end
end
end
Meq = [Meq;M];
rhseq = [rhseq;zeros(size(M,1),1)];
end
% -------------------------------------
% Integral equality constraint
if ~isempty(prescription.Integral)
% Simpson's rule over the support will be exact,
% evaluating at each midpoint. Don't forget that
% the effective "stepsize" is actually dx(i)/2
M = zeros(1,nc);
for i = 1:(nk-1)
M(1,i+[0 1]) = M(1,i+[0 1]) + dx(i)/6;
% the midpoint
M(1,i+[0 1 nk nk+1]) = M(1,i+[0 1 nk nk+1]) + ...
[1/2 , 1/2 , (1/8).*dx(i) , (-1/8).*dx(i)]*4*dx(i)/6;
end
Meq = [Meq;M];
rhseq = [rhseq;prescription.Integral];
end
% -------------------------------------
% single point inequality constraints
% left hand side minimum
if ~isempty(prescription.LeftMinValue)
M = zeros(1,nc);
M(1,1) = -1;
Mineq = [Mineq;M];
rhsineq = [rhsineq;-prescription.LeftMinValue];
end
% left hand side maximum
if ~isempty(prescription.LeftMaxValue)
M = zeros(1,nc);
M(1,1) = 1;
Mineq = [Mineq;M];
rhsineq = [rhsineq;prescription.LeftMaxValue];
end
% right hand side min
if ~isempty(prescription.RightMinValue)
M = zeros(1,nc);
M(1,nk) = -1;
Mineq = [Mineq;M];
rhsineq = [rhsineq;-prescription.RightMinValue];
end
% right hand side max
if ~isempty(prescription.RightMaxValue)
M = zeros(1,nc);
M(1,nk) = 1;
Mineq = [Mineq;M];
rhsineq = [rhsineq;prescription.RightMaxValue];
end
% left end slope min
if ~isempty(prescription.LeftMinSlope)
M = zeros(1,nc);
M(1,nk+1) = -1;
Mineq = [Mineq;M];
rhsineq = [rhsineq;-prescription.LeftMinSlope];
end
% left end slope max
if ~isempty(prescription.LeftMaxSlope)
M = zeros(1,nc);
M(1,nk+1) = 1;
Mineq = [Mineq;M];
rhsineq = [rhsineq;prescription.LeftMaxSlope];
end
% right end slope min
if ~isempty(prescription.RightMinSlope)
M = zeros(1,nc);
M(1,nc) = -1;
Mineq = [Mineq;M];
rhsineq = [rhsineq;-prescription.RightMinSlope];
end
% right end slope max
if ~isempty(prescription.RightMaxSlope)
M = zeros(1,nc);
M(1,nc) = 1;
Mineq = [Mineq;M];
rhsineq = [rhsineq;prescription.RightMaxSlope];
end
% -------------------------------------
% overall inequalities
% global min & max value
if ~isempty(prescription.MinValue) || ~isempty(prescription.MaxValue)
% sample points for necessary conditions are chebychev
% nodes. No good reason here, except that they are
% as good a set of points as any other.
tm = [0.017037; 0.066987; 0.14645; 0.25; 0.37059; ...
0.5; 0.62941; 0.75; 0.85355; 0.93301; 0.98296];
nsamp = length(tm);
ntot = nk+(nk-1)*nsamp;
Mmin = zeros(ntot,nc);
Mmax = Mmin;
% constrain values at knots
if ~isempty(prescription.MinValue)
Mmin(1:nk,1:nk) = -eye(nk,nk);
end
if ~isempty(prescription.MaxValue)
Mmax(1:nk,1:nk) = eye(nk,nk);
end
% and intermediate sample points
t2 = tm.^2;
t3 = tm.^3;
s2 = (1-tm).^2;
s3 = (1-tm).^3;
vals = [3*s2-2*s3 , 3*t2-2*t3 , -s3+s2 , t3-t2];
for j = 1:(nk-1)
if ~isempty(prescription.MinValue)
Mmin((1:nsamp) + (j-1)*nsamp + nk , j+[0 1 nk nk+1]) = ...
-vals*diag([1 1 dx(j) dx(j)]);
end
if ~isempty(prescription.MaxValue)
Mmax((1:nsamp) + (j-1)*nsamp + nk , j+[0 1 nk nk+1]) = ...
vals*diag([1 1 dx(j) dx(j)]);
end
end
if ~isempty(prescription.MinValue)
Mineq = [Mineq;Mmin];
rhsineq = [rhsineq;repmat(-prescription.MinValue,ntot,1)];
end
if ~isempty(prescription.MaxValue)
Mineq = [Mineq;Mmax];
rhsineq = [rhsineq;repmat(prescription.MaxValue,ntot,1)];
end
end
% global min or max slope
if ~isempty(prescription.MinSlope) || ~isempty(prescription.MaxSlope)
% its a transformed monotonicity problem. See the
% explanation provided below for sufficient conditions
% for monotonicity. The trick here is to build them
% for a minimum slope. Do that by effectively building
% an offset into the slope. I.e., assume that theta is the
% minimum slope allowed. Implicitly subtract the line
% theta*x from the curve, then use the standard linear
% system derived below for monotonicity. This will yield
% a system of linear constraints that are sufficient for
% a minimum (or maximum) slope of theta.
%
MS = zeros(5*(nk - 1) + nk,nc);
rhsMS = zeros(5*(nk-1) + nk,1);
% constrain the derivatives at each knot
for n = 1:nk
MS(n,nk + n) = -1;
end
rhsMS(1:nk) = -1;
n = nk + 1;
for j = 1:(nk - 1)
% The constraints applied for an increasing function
% are (in a form that lsqlin will like), i.e., A*x <= b
%
% -delta + theta <= 0
%
% (etc.)
MS(n,j + [0 1]) = [1 -1];
rhsMS(n) = -dx(j);
MS(n + 1,j + [0, 1, nk,nk + 1]) = [3, -3, [-1, 1]*dx(j)];
rhsMS(n + 1) = -3*dx(j);
MS(n + 2,j + [0, 1, nk,nk + 1]) = [3, -3, [1, -1]*dx(j)];
rhsMS(n + 2) = -3*dx(j);
MS(n + 3,j + [0, 1, nk,nk + 1]) = [9, -9, [1, 2]*dx(j)];
rhsMS(n + 3) = -6*dx(j);
MS(n + 4,j + [0, 1, nk,nk + 1]) = [9, -9, [2, 1]*dx(j)];
rhsMS(n + 4) = -6*dx(j);
n = n + 5;
end
% fold these inequalities into any others
if ~isempty(prescription.MinSlope)
theta = prescription.MinSlope;
Mineq = [Mineq;MS];
rhsineq = [rhsineq;theta*rhsMS];
end
if ~isempty(prescription.MaxSlope)
% the max slope case is as if we have negated
% the min slope case.
theta = prescription.MaxSlope;
Mineq = [Mineq;-MS];
rhsineq = [rhsineq;-theta*rhsMS];
end
end
% -------------------------------------
% SimplePeak and SimpleValley are really just composite
% properties. A peak at x == a is equivalent to a monotone
% increasing function for x<=a, and a monotone decreasing
% function for x>=a. Likewise a valley is just the opposite.
%
% It is best if a peak or valley occurs at a knot.
%
% So specifying a peak or a valley will cause any monotonicity
% specs to be ignored. slmset has already checked for this,
% and turned off any other monotonicity based properties.
sp = prescription.SimplePeak;
if isnumeric(sp) && ~isempty(sp)
prescription.Increasing = [knots(1),sp];
prescription.Decreasing = [sp,knots(end)];
end
sv = prescription.SimpleValley;
if isnumeric(sv) && ~isempty(sv)
prescription.Decreasing = [knots(1), sv];
prescription.Increasing = [sv,knots(end)];
end
% -------------------------------------
% monotonicity?
% increasing regions
totalmonotoneintervals = 0;
incR = prescription.Increasing;
L=0;
if ischar(incR)
if strcmp(incR,'on')
L=L+1;
mono(L).knotlist = [1,nk];
mono(L).direction = 1;
totalmonotoneintervals = totalmonotoneintervals + nk - 1;
end
elseif ~isempty(incR)
for i=1:size(incR,1)
L=L+1;
mono(L).knotlist = [max(1,sum(knots <= incR(i,1))), ...
min(nk,1 + sum(knots < incR(i,2)))]; %#ok
mono(L).direction = 1; %#ok
totalmonotoneintervals = totalmonotoneintervals + diff(mono(L).knotlist);
end
end
% decreasing regions
decR = prescription.Decreasing;
if ischar(decR)
if strcmp(decR,'on')
L=L+1;
mono(L).knotlist = [1,nk];
mono(L).direction = -1;
totalmonotoneintervals = totalmonotoneintervals + nk - 1;
end
elseif ~isempty(decR)
for i=1:size(decR,1)
L=L+1;
mono(L).knotlist = [max(1,sum(knots <= decR(i,1))), ...
min(nk,1 + sum(knots < decR(i,2)))];
mono(L).direction = -1;
totalmonotoneintervals = totalmonotoneintervals + diff(mono(L).knotlist);
end
end
if L>0
% there were at least some monotone regions specified
M = zeros(7*totalmonotoneintervals,nc);
n = 0;
for i=1:L
for j = mono(i).knotlist(1):(mono(i).knotlist(2) - 1)
% the function must be monotone between
% knots j and j + 1. The direction over
% that interval is specified. The constraint
% system used comes from Fritsch & Carlson, see here:
%
% http://en.wikipedia.org/wiki/Monotone_cubic_interpolation
%
% Define delta = (y(i+1) - y(i))/(x(i+1) - x(i))
% Thus delta is the secant slope of the curve across
% a knot interval. Further, define alpha and beta as
% the ratio of the derivative at each end of an
% interval to the secant slope.
%
% alpha = d(i)/delta
% beta = d(i+1)/delta
%
% Then we have an elliptically bounded region in the
% first quadrant that defines the set of monotone cubic
% segments. We cannot define that elliptical region
% using a set of linear constraints. However, by use
% of a system of 7 linear constraints, we can form a
% set of sufficient conditions such that the curve
% will be monotone. There will be some few cubic
% segments that are actually monotone, yet lie outside
% of the linear system formed. This is acceptable,
% as our linear approximation here is a sufficient
% one for monotonicity, although not a necessary one.
% It merely says that the spline may be slightly over
% constrained, i.e., slightly less flexible than is
% absolutely necessary. (So?)
%
% The 7 constraints applied for an increasing function
% are (in a form that lsqlin will like):
%
% -delta <= 0
% -alpha <= 0
% -beta <= 0
% -alpha + beta <= 3
% alpha - beta <= 3
% alpha + 2*beta <= 9
% 2*alpha + beta <= 9
%
% Multiply these inequalities by (y(i+1) - y(i)) to
% put them into a linear form.
M(n + 1,j+[0 1]) = [1 -1]*mono(i).direction;
M(n + 2,nk + j) = -mono(i).direction;
M(n + 3,nk + j + 1) = -mono(i).direction;
M(n + 4,j + [0, 1, nk,nk + 1]) = mono(i).direction*[3, -3, [-1, 1]*dx(j)];
M(n + 5,j + [0, 1, nk,nk + 1]) = mono(i).direction*[3, -3, [1, -1]*dx(j)];
M(n + 6,j + [0, 1, nk,nk + 1]) = mono(i).direction*[9, -9, [1, 2]*dx(j)];
M(n + 7,j + [0, 1, nk,nk + 1]) = mono(i).direction*[9, -9, [2, 1]*dx(j)];
n = n + 7;
end
end
Mineq = [Mineq;M];
rhsineq = [rhsineq;zeros(size(M,1),1)];
end
% -------------------------------------
% PositiveInflection and NegativeInflection are really just
% composite properties. A point of inflection at x == a is
% equivalent to a change in the sign of the curvature at
% x == a.
pinf = prescription.PositiveInflection;
if isnumeric(pinf) && ~isempty(pinf)
prescription.ConcaveDown = [knots(1),pinf];
prescription.ConcaveUp = [pinf,knots(end)];
end
ninf = prescription.NegativeInflection;
if isnumeric(ninf) && ~isempty(ninf)
prescription.ConcaveUp = [knots(1),ninf];
prescription.ConcaveDown = [ninf,knots(end)];
end
% -------------------------------------
% concave up or down?
% regions of positive curvature
cuR = prescription.ConcaveUp;
L=0;
if ischar(cuR)
if strcmp(cuR,'on')
L=L+1;
curv(L).knotlist = 'all';
curv(L).direction = 1;
curv(L).range = [];
end
elseif ~isempty(cuR)
for i=1:size(cuR,1)
L=L+1;
curv(L).knotlist = []; %#ok
curv(L).direction = 1; %#ok
curv(L).range = sort(cuR(i,:)); %#ok
end
end
% negative curvature regions
cdR = prescription.ConcaveDown;
if ischar(cdR)
if strcmp(cdR,'on')
L=L+1;
curv(L).knotlist = 'all';
curv(L).direction = -1;
curv(L).range = [];
end
elseif ~isempty(cdR)
for i=1:size(cdR,1)
L=L+1;
curv(L).knotlist = [];
curv(L).direction = -1;
curv(L).range = sort(cdR(i,:));
end
end
if L>0
% there were at least some regions with specified curvature
M = zeros(0,nc);
n = 0;
for i=1:L
if isempty(curv(L).range)
% the entire domain was specified to be
% curved in some direction
for j=1:(nk-1)
n=n+1;
M(n,j+[0 1]) = curv(i).direction*[6 -6]/(dx(j).^2);
M(n,nk+j+[0 1]) = curv(i).direction*[4 2]/dx(j);
end
n=n+1;
M(n,nk+[-1 0]) = curv(i).direction*[-6 6]/(dx(end).^2);
M(n,2*nk+[-1 0]) = curv(i).direction*[-2 -4]/dx(end);
else
% only enforce curvature between the given range limits
% do each knot first.
for j=1:(nk-1)
if (knots(j)<curv(i).range(2)) && ...
(knots(j)>=curv(i).range(1))
n=n+1;
M(n,j+[0 1]) = curv(i).direction*[6 -6]/(dx(j).^2);
M(n,nk+j+[0 1]) = curv(i).direction*[4 2]/dx(j);
end
end
% also constrain at the endpoints of the range
curv(i).range = max(min(curv(i).range(:),knots(end)),knots(1));
[junk,ind] = histc(curv(i).range,knots); %#ok
ind(ind==(nk))=nk-1;
t = (curv(i).range - knots(ind))./dx(ind);
s = 1-t;
for j = 1:numel(ind)
M(n+j,ind(j)+[0 1 nk nk+1]) = -curv(i).direction* ...
[(6 - 12*s(j))./(dx(ind(j)).^2), (6 - 12*t(j))./(dx(ind(j)).^2) , ...
(2 - 6*s(j))./dx(ind(j)), (6*t(j) - 2)./dx(ind(j))];
end
n = n + numel(ind);
end
end
Mineq = [Mineq;M];
rhsineq = [rhsineq;zeros(size(M,1),1)];
end
% -------------------------------------
% Jerk function increasing or decreasing?
% Just a constraint on the sign of the third derivative.
switch prescription.Jerk
case 'positive'
fpppsign = 1;
case 'negative'
fpppsign = -1;
otherwise
fpppsign = 0;
end
if fpppsign
% constrain the third derivative sign
M = zeros(nk-1,nc);
for i=1:(nk-1)
% one constraint on every interval
M(i,i+[0 1]) = fpppsign*[-12 12]/(dx(i).^3);
M(i,nk+i+[0 1]) = fpppsign*[-6 -6]/dx(i).^2;
end
Mineq = [Mineq;M];
rhsineq = [rhsineq;zeros(size(M,1),1)];
end
% -------------------------------------
% Error bar inequality constraints
if ~isempty(prescription.ErrorBar)
% lower bounds
Mineq = [Mineq;-Mdes];
rhsineq = [rhsineq;-(y-prescription.ErrorBar(:,1))];
% upper bounds
Mineq = [Mineq;Mdes];
rhsineq = [rhsineq;y+prescription.ErrorBar(:,2)];
end
% -------------------------------------
% Use envelope inequalities?
switch prescription.Envelope
case 'off'
% nothing to do in this case
case 'supremum'
% yhat - y >= 0
Mineq = [Mineq;-Mdes];
rhsineq = [rhsineq;-rhs];
case 'infimum'
% yhat - y <= 0
Mineq = [Mineq;Mdes];
rhsineq = [rhsineq;rhs];
end
% -------------------------------------
% scale equalities for unit absolute row sum
if ~isempty(Meq)
rs = diag(1./sum(abs(Meq),2));
Meq = rs*Meq;
rhseq = rs*rhseq;
end
% scale inequalities for unit absolute row sum
if ~isempty(Mineq)
rs = diag(1./sum(abs(Mineq),2));
Mineq = rs*Mineq;
rhsineq = rs*rhsineq;
end
% -------------------------------------
% now worry about the regularization. There are three
% possible cases.
% 1. We have a given regularization parameter
% 2. We have a given rmse that we wish to match
% 3. We must use cross validation to choose the parameter
RP = prescription.Regularization;
if (isnumeric(RP) && (RP>=0)) || ((ischar(RP)) && (strcmpi(RP,'smoothest')))
% solve the problem using the given regularization parameter
coef = solve_slm_system(RP,Mdes,rhs,Mreg,rhsreg, ...
Meq,rhseq,Mineq,rhsineq,prescription);
elseif isnumeric(RP) && (RP<0)
% we must match abs(RP) as the rmse.
aim_rmse = abs(RP);
fminbndoptions = optimset('fminbnd');
RP = fminbnd(@match_rmse,-6,6,fminbndoptions, ...
Mdes,rhs,Mreg,rhsreg,Meq,rhseq,Mineq,rhsineq, ...
aim_rmse,prescription);
% we logged the parameter in the optimization. undo that for
% the final call
RP = 10^RP;
% do one final call to get the final coefficients
coef = solve_slm_system(RP,Mdes,rhs,Mreg,rhsreg, ...
Meq,rhseq,Mineq,rhsineq,prescription);
elseif ischar(RP)
% its a cross validation problem to solve. drop out each point
% in turn from the model, then use fminbnd to minimize the
% predicted sum of squares at the missing points.
fminbndoptions = optimset('fminbnd');
% fminbndoptions.Display = 'iter';
RP = fminbnd(@min_cv,-6,6,fminbndoptions, ...
Mdes,rhs,Mreg,rhsreg,Meq,rhseq,Mineq,rhsineq,prescription);
% we logged the parameter in the optimization. undo that for
% the final call
RP = 10^RP;
% do one final call to get the final coefficients
coef = solve_slm_system(RP,Mdes,rhs,Mreg,rhsreg, ...
Meq,rhseq,Mineq,rhsineq,prescription);
end
% -------------------------------------
% unpack coefficients into the result structure
slm.form = 'slm';
slm.degree = 3;
slm.knots = knots;
slm.coef = reshape(coef,nk,2);
% generate model statistics
slmstats.TotalDoF = 2*nk;
slmstats.NetDoF = slmstats.TotalDoF - size(Meq,1);
% this function does all of the stats, stuffing into slmstats
slm.stats = modelstatistics(slmstats,Mdes,y,coef);
% ========================================================
% ========= basic linear system solve ====================
% ========================================================
function [coef,lambda] = solve_slm_system(RP,Mdes,rhs,Mreg,rhsreg,Meq,rhseq,Mineq,rhsineq,prescription)
% solves the final linear system of equations for the model coefficients
if prescription.Verbosity > 1
% Linear solve parameters
disp('=========================================')
disp('LINEAR SYSTEM SOLVER')
disp(['Design matrix shape: ',num2str(size(Mdes))])
disp(['Regularizer shape: ',num2str(size(Mreg))])
disp(['Equality constraints: ',num2str(size(Meq))])
disp(['Inequality constraints: ',num2str(size(Mineq))])
end
% combine design matrix and regularizer
if strcmpi(prescription.Regularization,'smoothest')
% find the smoothest possible solution
Mdes = [0.000001*Mdes;Mreg];
rhs = [0.000001*rhs;rhsreg];
else
Mdes = [Mdes;RP*Mreg];
rhs = [rhs;RP*rhsreg];
end
if prescription.Verbosity > 1
disp(' ')
disp(['Condition number of the regression: ',num2str(cond(Mdes))])
disp(' ')
end
% -------------------------------------
% solve the linear system
if isempty(Mineq) && isempty(Meq)
coef = Mdes\rhs;
lambda.eqlin=[];
lambda.ineqlin=[];
elseif isempty(Mineq)
% with no inequality constraints, lse is faster than
% is lsqlin. This also allows the use of slm when
% the optimization toolbox is not present if there
% are no inequality constraints.
coef = lse(Mdes,rhs,full(Meq),rhseq);
else
% use lsqlin. first, set the options
options = optimset('lsqlin');
if prescription.Verbosity > 1
options.Display = 'final';
else
options.Display = 'off';
end
% the Largescale solver will not allow general constraints,
% either equality or inequality
options.LargeScale='off';
% and solve
[coef,junk,junk,exitflag,junk,lambda] = ...
lsqlin(Mdes,rhs,Mineq,rhsineq,Meq,rhseq,[],[],[],options); %#ok
% was there a feasible solution?
if exitflag == -2
coef = nan(size(coef));
end
end
% ========================================================
% ========== chose RP to hit aim rmse ====================
% ========================================================
function delta_rmse = match_rmse(RP,Mdes,rhs,Mreg,rhsreg,Meq,rhseq,Mineq,rhsineq,aim_rmse,prescription)
% returns delta rmse for a given regularization parameter
% the regularization parameter has been logged.
% exponentiate for the solve
RP = 10^RP;
% solve the system at the given RP
coef = solve_slm_system(RP,Mdes,rhs,Mreg,rhsreg, ...
Meq,rhseq,Mineq,rhsineq,prescription);
% compute rmse
rmse = sqrt(mean((rhs - Mdes*coef).^2));
% delta rmse error from aim
delta_rmse = abs(aim_rmse - rmse);
% ========================================================
% ======= chose RP to minimize cv prediction error =======
% ========================================================
function cvpred = min_cv(RP,Mdes,rhs,Mreg,rhsreg,Meq,rhseq,Mineq,rhsineq,prescription)
% returns delta rmse for a given regularization parameter
% the regularization parameter has been logged.
% exponentiate for the solve
RP = 10^RP;
% solve the system at the given RP, dropping out one point
% at a time from the data.
n = size(Mdes,1);
uselist = true(n,1);
cvpred = 0;
for i=1:n
ul = uselist;
ul(i) = false;
coef = solve_slm_system(RP,Mdes(ul,:),rhs(ul),Mreg, ...
rhsreg,Meq,rhseq,Mineq,rhsineq,prescription);
pred = Mdes(i,:)*coef;
cvpred = cvpred + (pred - rhs(i)).^2;
end
% ========================================================
% ========= choose knot placement for free knots =========
% ========================================================
function rss = free_knot_obj(intknots,x,y,prescription)
% returns residual sum of squares for fmincon to work on
% insert the interior knots from the optimizer
prescription.Knots(2:(end-1)) = intknots;
slm = slmengine(x,y,prescription);
% get predictions, form residuals, square, and sum
switch prescription.Result
case 'slm'
rss = sum((slmeval(x,slm,0) - y).^2);
case 'pp'
rss = sum((ppval(slm,x) - y).^2);
end
% ========================================================
% ========= set scaling as necessary =========
% ========================================================
function [ytrans,prescription] = scaleproblem(x,y,prescription)
% chooses an appropriate scaling for the problem
if strcmp(prescription.Scaling,'on')
% scale y so that the minimum value is 1, and the maximum value 2.
else
% no scaling is done
prescription.yshift = 0;
prescription.yscale = 1;
ytrans = y;
end
% ========================================================
% ========= chooseknots =========
% ========================================================
function [knots,nk] = chooseknots(K,n,x)
% choose every K'th data point as a knot
if K > 0
% just given a number of knots, equally spaced
nk = K;
knots = linspace(min(x),max(x),nk)';
knots(end)=max(x); % just to make sure
else
% we need every abs(K)'th knot
if (K == 0) || (K >= n)
error('SLMENGINE:knots','Every K''th data point was specified, but K was too large or zero')
end
Kind = 1:abs(K):n;
% in case K was not an integer, or did not
% go to the very end.
if (n - Kind(end)) < (K/2)
% expand the last knot interval
Kind(end) = n;
else
% append one extra knot for the last data point
Kind(end + 1) = n;
end
% rounding makes them into indices
Kind = round(Kind);
% in case of replicate data points
x = sort(x);
knots = unique(x(Kind));
nk = numel(knots);
end
% ========================================================
% ========= model statistics =========
% ========================================================
function slmstats = modelstatistics(slmstats,Mdes,y,coef)
% generate model statistics, stuffing them into slmstats
% residuals, as yhat - y
resids = Mdes*coef - y;
% RMSE: Root Mean Squared Error
slmstats.RMSE = sqrt(mean(resids.^2));
% R-squared
slmstats.R2 = 1 - sum(resids.^2)./sum((y - mean(y)).^2);
% adjusted R^2
ndata = numel(y);
slmstats.R2Adj = 1 - (1-slmstats.R2)*(ndata - 1)./(ndata - slmstats.NetDoF);
% range of the errors, min to max, as yhat - y
slmstats.ErrorRange = [min(resids),max(resids)];
% compute the 25% and 75% points (quartiles) of the residuals
% (This is consistent with prctile, from the stats TB.)
resids = sort(resids.');
ind = 0.5 + ndata*[0.25 0.75];
f = ind - floor(ind);
ind = min(ndata - 1,max(1,floor(ind)));
slmstats.Quartiles = resids(ind).*(1-f) + resids(ind+1).*f;
% ========================================================
% ========== check for inappropriate properties ==========
% ========================================================
function property_check(prescription,model_degree)
% issues warning messages for inappropriate properties
% for the given model degree.
switch model_degree
case {3 'cubic'}
% no properties flagged for cubic models
case {1 'linear'}
% only (some) curvature properties are flagged
% for the linear model
if ~isempty(prescription.Jerk)
warning('SLMENGINE:ignoredconstraint', ...
'Property ignored for a linear model: Jerk')
end
if ~isempty(prescription.XYPP)
warning('SLMENGINE:ignoredconstraint', ...
'Property ignored for a linear model: XYPP')
end
if ~isempty(prescription.XYPPP)
warning('SLMENGINE:ignoredconstraint', ...
'Property ignored for a linear model: XYPPP')
end
case {0 'constant'}
% curvature properties are flagged for
% the constant model, as well as most slope
% related properties
if ~isempty(prescription.Jerk)
warning('SLMENGINE:ignoredconstraint', ...
'Property ignored for a constant model: Jerk')
end
if ~isempty(prescription.LeftMaxSlope)
warning('SLMENGINE:ignoredconstraint', ...
'Property ignored for a constant model: LeftMaxSlope')
end
if ~isempty(prescription.LeftMinSlope)
warning('SLMENGINE:ignoredconstraint', ...
'Property ignored for a constant model: LeftMinSlope')
end
if ~isempty(prescription.LeftSlope)
warning('SLMENGINE:ignoredconstraint', ...
'Property ignored for a constant model: LeftSlope')
end
if ~isempty(prescription.LinearRegion)
warning('SLMENGINE:ignoredconstraint', ...
'Property ignored for a constant model: LinearRegion')
end
if ~isempty(prescription.MaxSlope)
warning('SLMENGINE:ignoredconstraint', ...
'Property ignored for a constant model: MaxSlope')
end
if ~isempty(prescription.MinSlope)
warning('SLMENGINE:ignoredconstraint', ...
'Property ignored for a constant model: MinSlope')
end
if ~isempty(prescription.NegativeInflection)
warning('SLMENGINE:ignoredconstraint', ...
'Property ignored for a constant model: NegativeInflection')
end
if ~isempty(prescription.PositiveInflection)
warning('SLMENGINE:ignoredconstraint', ...
'Property ignored for a constant model: PositiveInflection')
end
if ~isempty(prescription.RightMaxSlope)
warning('SLMENGINE:ignoredconstraint', ...
'Property ignored for a constant model: RightMaxSlope')
end
if ~isempty(prescription.RightMinSlope)
warning('SLMENGINE:ignoredconstraint', ...
'Property ignored for a constant model: RightMinSlope')
end
if ~isempty(prescription.RightSlope)
warning('SLMENGINE:ignoredconstraint', ...
'Property ignored for a constant model: RightSlope')
end
if ~isempty(prescription.XYP)
warning('SLMENGINE:ignoredconstraint', ...
'Property ignored for a constant model: XYP')
end
if ~isempty(prescription.XYPP)
warning('SLMENGINE:ignoredconstraint', ...
'Property ignored for a constant model: XYPP')
end
if ~isempty(prescription.XYPPP)
warning('SLMENGINE:ignoredconstraint', ...
'Property ignored for a constant model: XYPPP')
end
end
