The following script is a tutorial on the methodology of what I call "Shape Prescriptive Modeling". Its a tool for modeling a function of a single variable
y = f(x) + error
The actual model will be in the form of a least squares spline. Shape prescriptive modeling is really a way of thinking about a model - it has a lot of very Bayesian underpinnings. I'll get to that in a bit. First, lets talk about modeling, who builds models, and why we choose a specific model form.
Why do we, as scientists and engineers, fit a curve to data? There are really two basic reasons why, if we ignore the common "My boss told me to do this." The first such reason is for purely predictive value. We may wish to reduce some set of data to a simply accessible functional relationship for future use. The process of modeling may allow us to smooth the data, producing a smooth curve where only distinct data points fell before. As well, this smooth curve is now a platform for interpolation. Note that if we have smoothed the data, removing unwanted noise that was present in the data, then the interpolant is more an approximant. I.e., one will not interpolate the data, but a smoother form of it.
The second common goal of curve fitting is often to predict/estimate some parameter of the relationship, such as a minimum or maximum value where the curve approaches an asymptote, or perhaps a maximum slope or a rate parameter.
If we are to fit a model to data, what kind of models do we use, and how does the model class reflect our goals in the modeling process? To answer this question, I'll first look at the common modeling tools we might see. First and most primeval is polyfit. Polyfit allows us to fit a polynomial model to our data. Such a model tells us little about our process beyond perhaps a simple slope. These models are really of value only for very simple relationships, where we are willing to invest little in the modeling process. Many of the individuals who use polynomial modeling tools do so for the wrong reasons. They use higher order polynomial fits to get a better approximation to their data, not realizing the problems inherent with those higher order polynomials.
At the other end of the spectrum, one sees individuals using nonlinear regression to fit a large variety of curves to data. Exponential terms, gaussian-like modes, sigmoid shapes, as well as many more. Did these model forms result from valid mechanistic arguments? I.e., has the analyst developed a real model of the process under study, where all that is left is to estimate the parameters of that model from data?
Sometimes it is so, but far more often the analyst individual knows something about the basic underlying process, and the model chosen has the correct basic functional shape. I'll argue that truly valid mechanistic models are nearly as rare as hen's teeth.
There is a subtle variation on the mechanistic model. I call it the metaphoric or metaphorical model. Here we use a mathematical model of one process that we do understand, used as a model for a process that we do not really understand. My favorite examples of such metaphorical models are cubic splines, where a mathematical model of a thin flexible beam is used to predict many other processes. A second good example is the use of epidemiological models as used to predict sales of some product. There are many other examples of course.
An advantage of a metaphorical model is it may help us to understand/ predict/even extrapolate future behavior based on our knowledge of the metaphor. Many of the nonlinear regression models that we see built in Matlab are of the shape variety. I.e., we see a bell shaped curve, so we fit a gaussian to it. We see a function with a lower and upper constant asymptote, and we fit some variety of logistic function, or an erf variant to it. If the fit is inadequate for our purposes, we try another functional form. I'll argue that this group of individuals is often using nonlinear regression for the wrong reason. They simply want some curve that behaves the way they know it should. Examples of this might be a photographic film modelor, who knows the film darkens as more photons land upon it, or a biologist, who knows that the count of bacteria in a petri dish will grow with time. (All of these examples have limits, but a good scientist/engineer will know those limits.)
The fact is, we often understand something about a process we are trying to model. What we may not have is the ability to build that knowledge and understanding into our model. Even if we know the underlying functional relationship is monotone increasing/decreasing, how do we estimate a model that has that behavior built into it? This is where shape prescriptive modeling shines. It is really a very Bayesian concept. It is the process of quantifying our basic knowledge of a process in mathematical terms, and then estimating a predictive model of the process that has that knowledge built into it. Our knowledge becomes a prescription for the final shape of the process. Sometimes we know much about a process, in which case we have a detailed prescription. Other times we have a very sketchy prescription.
The slm tools in this toolkit are designed to help the user to embed their knowledge into a model using a langauage of shape. One builds a prescription for the desired shape, then combines it with their data to build a model. These steps can be done distinctly, where an explicit prescription is created, then applied with a modeling engine to data, or the prescription can be built in a fluid manner using a gui. Here a point and click interface can provide the information necessary.
Some final notes before we begin actual experimentation: This toolkit assumes possession of the optimization toolbox.
Also, all examples will use the commandline engine for modeling, instead of the gui interface. One good reason is that at the time of writing this tutorial, I had not finished writing the gui. A better reason is that anything that can be done via the gui is also possible to do using the computational engine, and cell mode favors this approach.
Finally, ...
Contents
- Some data
- Polyfit?
- Less data, more noise, with polyfit
- A piecewise linear SLM fit
- Shape prescriptions
- Use of slmset
- Evaluating a model
- Plotting a model using plotslm
- slmengine can plot too
- slmpar: Computations on a spline model, min max, integral, etc.
- slmsolve: find all roots of a spline, f(x) == y, for any given value of y
- Model fitting using shape prescriptors
- Curvature constraints
- Both curvature and monotonicity constraints
- Minimum slope or maximum slope constraints
- Extrapolation in the fit, non-uniform knot spacings
- Piecewise constant models
- Fixed specific values
- Flawed information, "Garbage in, garbage out"
- Too many knots?
- Use in practice
- Returning a set of points, evaluated through the fitted curve.
- Combining various constraints to produce a better model
- Extrapolation of a spline model
- You can get information about the fit itself
- Knot placement optimization
Some data
Lets start out by creating some data. A simple functional form that we understand the shape of.
n = 100;
x1 = sort(rand(n,1));
y1 = exp(x1) + randn(n,1)/20;
plot(x1,y1,'ro')
Polyfit?
First, lets build a simple model using polyfit.
P = polyfit(x1,y1,6); xev = 0:.001:1; yev = polyval(P,xev);
and plot the results on top of the data
plot(xev,yev,'g-',x1,y1,'ro')
Are you happy with the resulting curve? Perhaps. Much depends upon our goals in building a model for this process. How good of a fit do we really need, how much noise is there in the data? Does the model behave as we believe it should? Does this model have lack of fit? Actually, it looks pretty good to me, at least given the amount of noise in the data. But see what happens when I use a high order polynomial model.
This one was just too much, with condition problems.
P = polyfit(x1,y1,15); xev = 0:.001:1; yev = polyval(P,xev);
Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT.
and plot the results on top of the data
plot(xev,yev,'g-',x1,y1,'ro')
Less data, more noise, with polyfit
Lets try again, only this time we'll use less data, with a much higher level of noise.
close n = 100; x2 = sort(4*rand(n,1)-3); y2 = exp(x2) + randn(n,1); plot(x2,y2,'ro') P = polyfit(x2,y2,6); xev = linspace(min(x2),max(x2),100); yhat = polyval(P,xev); hold on plot(xev,yhat,'g-') hold off
Its a clear example of overfitting this noisy data.
A piecewise linear SLM fit
Now lets try using the slm tools. After all, its why we are working through this tutorial.
slm = slmengine(x2,y2,'degree',1,'plot','on')
slm =
form: 'slm'
degree: 1
knots: [6x1 double]
coef: [6x1 double]
stats: [1x1 struct]
prescription: [1x1 struct]
x: [100x1 double]
y: [100x1 double]
Extrapolation: 'constant'
note that the result from slmengine is a structure. It contains 7 fields, and has much the same information as the 'pp' form does.
'form' is (in this case) 'slm'. This means that the function this structure contains is a slm model. We could have had it return a pp form instead.
'degree' is the polynomial degree of the piecewise segments. 1 indicates a piecewise linear function, continuous across the breakpoints (or knots.) Note that the splines toolbox tends to use the term 'order', one more than the degree.
'knots' contains the set of knots or breakpoints. The default is to use 6 equally spaced knots, from min to max of the data.
'coef' is the set of coefficients of the spline model. I've used a Hermite form here. Its an easy one to understand, as well as to look at the coefficients and visualize the shape of the curve. The field coef contains the value of the function at the corresponding knot.
'prescription' is the default prescription. A bayesian might call it an uninformative prior. Its the list of all settings for every optional parameter one can set for a model. Since we specified nothing but the data in our call to slmengine, we get all defaults.
Finally, there are fields for x and y. The model thus contains the data used to build it. This serves a documentation purpose but it also allows plotslm to work well later.
Shape prescriptions
We can look at the parameters in slm.prescription, but now is as good a tiem as any to see what slmset does for us. Here are the default set of properties.
prescription = slmset
prescription =
C2: 'on'
ConcaveDown: 'off'
ConcaveUp: 'off'
ConstantRegion: []
Decreasing: 'off'
Degree: 3
EndConditions: 'estimate'
Envelope: 'off'
ErrorBar: []
Extrapolation: 'constant'
Increasing: 'off'
Integral: []
InteriorKnots: 'fixed'
Jerk: ''
Knots: 6
LeftMaxSlope: []
LeftMaxValue: []
LeftMinSlope: []
LeftMinValue: []
LeftSlope: []
LeftValue: []
LinearRegion: []
MaxSlope: []
MaxValue: []
MinSlope: []
MinValue: []
NegativeInflection: []
Order: []
Plot: 'off'
PositiveInflection: []
Predictions: 1001
Regularization: 0.0001
Result: 'slm'
RightMaxSlope: []
RightMaxValue: []
RightMinSlope: []
RightMinValue: []
RightSlope: []
RightValue: []
Robust: 'off'
Scaling: 'on'
SegmentConstant: []
SegmentLinear: []
SimplePeak: []
SimpleValley: []
SumResiduals: []
Verbosity: 0
Weights: []
XY: []
XYP: []
XYPP: []
XYPPP: []
YScale: []
YShift: []
slmset is the tool that allows us to specify a variety of pieces of information to the engine.
Use of slmset
Help slmset will list all the options. It works on a property/value pair interface. There are a wide variety of properties that one can set, enough to satisfy almost any type of information a user can bring to a modeling problem. Here are the options one can set:
help slmset
slmset: defines the shape prescription for a model
usage 1: prescription=slmset(prop1,val1,pro2,val2,...)
usage 2: prescription=slmset(prescription,prop1,val1,pro2,val2,...)
A set of property/value pairs are parsed into a structure
to use as the prescription for slmengine and slmfit.
arguments:
prop, val - property/value pairs (see below for the complete list)
There is no limit on the upper number of property/value
pairs allowed.
Property names may be shortened, as long as they are
unambiguous. Case is ignored.
prescription - shape prescription structure, defining that
which the user is willing to say about the model.
The fields of the prescription structure reflect
the defaults, plus any modifications specified.
This structure will control the behaviours built
a curve fit by slmengine or the defaults for slmfit.
Property/value pairs:
Properties are character strings, chosen to be mnemonic of their
purpose. Any property name may be shortened, as long as the
shortened string is unambiguous. Thus since no other property
starts with the letter k, any of these alternatives are acceptable
shortenings for the 'knots' property: 'knot', 'kno', 'kn' or 'k'.
In the event that a given property is assigned more than once in the
list of property/value pairs, only the last value in the list is
assigned.
Property names & admissable values:
'C2 - Most users want their cubic spline curves to be as
smooth and differentiable as possible. For a cubic
spline, this is called a twice continuously differentiable
function, here simply 'C2'. Of course, if your spline is not
cubic, this specification is blithely ignored.
= 'on' --> Causes the result to be twice differentiable
across the knots.
= 'off' --> Allows a cubic spline to be only C1. The second
derivative of the spline MAY have discontinuities across
each knot of the spline.
DEFAULT VALUE: 'on'
Comment: I would rarely recommend changing this option
from the default.
'concavedown' - controls curvature of the function
= 'off' --> No part of the spline is constrained to be a
concave down function (i.e., a negative second derivative.)
= 'on' --> f''(x) >= 0 over the entire domain of the spline.
= vector of length 2 denoting the start and end points of a
region of the spline over which the second derivative is
negative.
= array of size nx2, each row of which denotes the start and
end points of a region of the spline over which the second
derivative is negative.
DEFAULT VALUE: 'off'
Comment: curvature properties do not apply to piecewise
constant functions.
Comment: This constraint is equivalent to a monotone
decreasing slope over the specified region.
'concaveup' - controls curvature of the function
= 'off' --> No part of the spline is constrained to be a
concave up function (i.e., a positive second derivative.)
= 'on' --> f''(x) >= 0 over the entire domain of the spline.
= vector of length 2 denoting the start and end points of a
region of the spline over which the second derivative is
positive.
= array of size nx2, each row of which denotes the start and
end points of a region of the spline over which the second
derivative is positive.
DEFAULT VALUE: 'off'
Comment: curvature properties do not apply to piecewise
constant functions.
Comment: This constraint is equivalent to a monotone
increasing slope over the specified region.
'constantregion' - defines a region over which the curve is forced
be constant, although the exact level is not defined.
= [] --> No region of the spline is forced to be a constant
function.
= vector of length 2 denoting the start and end points of a
region of the spline over which it is a constant function.
= array of size nx2, each row of which denotes the start and
end points of a region of the spline over which it is a
constant function.
DEFAULT VALUE: []
Comments: A segment which is forced to be constant over
only part of a knot interval must necessarily be constant
over that entire interval, since the curve is composed of
a single polynomial segment in a knot interval.
'decreasing' - controls monotonicity of the function
= 'off' --> No part of the spline is constrained to be a
decreasing function.
= 'on' --> the function will be decreasing over its entire domain.
= vector of length 2 --> denotes the start and end points of a
region of the curve over which it is monotone decreasing.
= array of size nx2 --> each row of which denotes the start
and end points of a region of the curve over which it is
monotone decreasing.
DEFAULT VALUE: 'off'
Comments: in actuality this property should be named
'non-increasing', since a constant function is admissible.
In addition, it is a sufficient constraint for monotonicity.
It is not a necessary constraint. There may exist another
spline which has a slightly lower sum of squares and is also
monotone.
'degree' - controls the degree of the piecewise Hermite function
= 'constant' --> Use a piecewise constant "Hermite"
= 'linear' --> Use a piecewise linear Hermite
= 'cubic' --> Use a piecewise cubic Hermite
As a concession to memory, valid synonyms for 'constant'
'linear', and 'cubic' are respectively the integers 0, 1, & 3
DEFAULT: 'cubic'
Comment: Some properties are inappropriate for all degrees
of function. E.g., it would be silly to specify a specific
value for the left hand end point slope of a piecewise
constant function. All information supplied will be used
to whatever extent possible.
The "order" of a form, as used by the spline toolbox, will
be one more than the degree. 'order' is also a valid property
in these tools, but it results only in setting the degree
field, where degree = order - 1.
'endconditions' - controls the end conditions applied to the spline
= 'natural' --> The "natural" spline conditions will be applied.
I.e., f''(x) = 0 at end end of the spline.
= 'notaknot' --> Not-a-knot end conditions applied.
= 'periodic' --> Periodic end conditions applied.
= 'estimate' --> end conditions are estimated from the data.
DEFAULT VALUE: 'estimate
Comment: Except for periodicity, end conditions are not
relevant to any degree model below cubic.
Periodic end conditions mean that the function values are
the same at each end of the curve for linear and cubic fits.
For cubic fits, periodicity means that the function has
also first and second derivative continuity across the
wrapped boundary.
For piecewise constant fits, end conditions do not apply,
and are ignored.
'envelope' - allows the user to solve for an envelope of the data
= 'off' --> the curve will be a simple least squares spline
= 'supremum' --> comute a model such that all residuals
(yhat - y) are positive. In effect, the curve will be a
"supremum" (least upper bound) function.
= 'infimum' --> comute a model such that all residuals
(yhat - y) are negative. In effect, the curve will be a
"infimum" (greatest lower bound) function.
DEFAULT VALUE: 'off'
Comment: this option should be rarely used, but its a cute
idea when it does come up.
'errorbar' - defines a set of lower and upper bounds for the function
value of the spline at each data point. If there are n data
points, then if the corresponding value is:
= [] --> No errorbar constraints will be imposed
= a scalar E --> error bars will be set at [Y-E,Y+E]
= a vector E --> error bars will be set at [Y-E,Y+E]
= an nx2 array --> error bars will be set at [Y-E(:,1),Y+E(:,2)]
DEFAULT VALUE: []
Comment: It is possible that depending on the choice of
knots, it will be impossible to satisfy some sets of
error bars. It may be best to simply use every single data
point as a knot. Of course, replicate data points with
non-overlappng error bars will always cause a failure.
'extrapolation' - tells SLMEVAL how to extrapolate the function outside
the defined knot points of the spline. This parameter has
NOTHING to do with how the curve will be fit, but only
with how it will be evaluated after the fit is done.
= 'error' --> No extrapolation will be allowed. An error
will be returned if the user tries to evaluate the
spline at any location outside of the knots.
= 'warning' --> A warning message will be generated when
extrapolation is attempted, but extrapolation will still
done as a constant value outside of the end point knots.
= 'constant' --> Constant extrapolation will be performed,
using the value of the spline at the associated end point
knot. This is the default, to be consistent with the
behavior of this tool in previous releases.
= 'linear' --> Linear extrapolation will be performed,
using the slope of the spline at the associated end point
knot.
= 'cubic' --> Cubic extrapolation will be performed,
using the shape of the spline at the associated end point
knot.
= 'NaN' --> NaNs will be inserted for any attempted extrapolation.
although no error or warning message will be generated.
DEFAULT value: 'constant'
'increasing' - controls monotonicity of the function
= 'off' --> No part of the spline is constrained to be an
increasing function.
= 'on' --> the function will be increasing over its entire domain.
= vector of length 2 --> denotes the start and end points of a
region of the curve over which it is monotone increasing.
= array of size nx2 --> each row of which denotes the start
and end points of a region of the curve over which it is
monotone increasing.
DEFAULT VALUE: 'off'
Comments: in actuality this property should be named
'non-decreasing', since a constant function is admissible.
In addition, it is a sufficient constraint for monotonicity.
It is not a necessary constraint. There may exist another
spline which has a slightly lower sum of squares and is also
monotone.
'integral' - known aim value for the integral of the curve over
its domain.
DEFAULT VALUE: []
'interiorknots' - allows for free knot placement (of the interior knots)
= 'fixed'
spaced knots. In this case, the first data point will
be the first knot.
= 'free' --> uses fmincon to optimize the interior knot
placement to minimize the overall rmse.
DEFAULT VALUE: 'fixed'
Comment: The initial values for the knot placement are
taken from the knots property.
Comment: Since the free knot placement is done by an
optimizer (fmincon), it is not allowed to both choose
a set of free knots and set the rmse of the fit.
Comment: The first and last knots are not adjusted by
the optimization, so there must be at least 3 knots.
Comment: Piecewise constant functions sometimes have
difficulty with free knots.
'jerk' - Controls the sign of the jerk function over the support
of the spline, to be either 'positive' or 'negative' (or no
constraint at all applied.) This a global sign
constraint on the third derivative of the function.
= 'positive' --> causes the third derivative to be
everywhere positive over the support of the function
= 'negative' --> causes the third derivative to be
everywhere negative over the support of the function
= '' --> No constraint applied
DEFAULT VALUE: ''
Comment: This property can be used to constrain the
curvature of the function to be increasing (or decreasing)
over the support, in combination with the concaveup or
concavedown properties as appropriate.
Comment: The third derivative is sometimes known as the
'jolt', but from my experience, 'jerk' seems to be the
most common name. I will add that the third derivative
quite difficult to estimate for a spline model, especially
if the data is at all noisy.t
'knots' - controls the number of knots used, or the number of
equally spaced knots.
= A scalar integer which denotes the number of equally
spaced knots. In this case, the first and last data
points will be the first and last knots.
= A vector containing the list of knots themselves.
The knots must be distinct and MUST wholly contain
the data or an error will be generated. No replicate
knots are allowed.
= A negative value K (no larger in absolute value than
the number of data points - 1) causes every K'th data
point to be used as a knot.
DEFAULT VALUE: 6
'leftmaxslope' - controls the maximum slope allowed at the left hand
end of the curve.
= [] --> No explicit value provided for the maximum slope
of the spline at its left hand end point.
= A numeric scalar --> the slope of the function will be
constrained to not rise above this value at its left
hand end point (i.e., the first knot.)
DEFAULT VALUE: []
'leftmaxvalue' - controls the maximum valued allowed at the left hand
end of the curve.
= [] --> No explicit value provided for the maximum value
of the spline at its left hand end point.
= A numeric scalar --> the function will be constrained
to not rise above this value at its left hand end point
(i.e., the first knot.)
DEFAULT VALUE: []
'leftminslope' - controls the minimum slope allowed at the left hand
end of the curve.
= [] --> No explicit value provided for the minimum slope
of the spline at its left hand end point.
= A numeric scalar --> the slope of the function will be
constrained to not fall below this value at its left
hand end point (i.e., the first knot.)
DEFAULT VALUE: []
'leftminvalue' - controls the minimum valued allowed at the left hand
end of the curve.
= [] --> No explicit value provided for the minimum value
of the spline at its left hand end point.
= A numeric scalar --> the function will be constrained
to not fall below this value at its left hand end point
(i.e., the first knot.)
DEFAULT VALUE: []
'leftslope' - controls the function slope at the left hand endpoint.
= [] --> No explicit value provided for the slope of the
curve at its left hand end point.
= A numeric scalar --> the function will be assigned this
slope at its left hand end point (i.e., the first knot.)
DEFAULT VALUE: []
'leftvalue' - controls the function value at its left hand endpoint.
= [] --> No explicit value provided for the value of the
curve at its left hand end point.
= A numeric scalar --> the function will be assigned this
value at its left hand end point (i.e., the first knot.)
DEFAULT VALUE: []
'linearregion' - defines a region over which the curve is forced
be linear, although the exact level is not defined.
= [] --> No region of the spline is forced to be a linear
function.
= vector of length 2 denoting the start and end points of a
region of the spline over which it is a linear function.
= array of size nx2, each row of which denotes the start and
end points of a region of the spline over which it is a
linear function.
DEFAULT VALUE: []
Comment 1: A segment which is forced to be linear over
only part of a knot interval must necessarily be constant
over that entire interval, since the curve is composed of
a single polynomial segment in a knot interval.
Comment 2: A linear region may extend across knots, in
which case the slope will take on the same value across
knot boundaries.
'maxslope' - controls the globally maximum slope of the function
= [] --> No explicit value provided for the globally maximum
slope of the curve.
= A numeric scalar --> the globally maximum slope of the spline
DEFAULT VALUE: []
Comment: This is a sufficient constraint for the maximum
slope of the spline. It is not a necessary constraint.
There may exist another spline which has a slightly lower
sum of squares and also has the same maximum slope.
'maxvalue' - controls the globally maximum value of the curve
= [] --> No explicit value provided for the globally
maximum value of the spline.
= A numeric scalar --> the function must lie no higher
than this maximum value.
DEFAULT VALUE: []
Comment 1: This constraint is only a necessary constraint.
It is not sufficient. In some circumstances the spline may
pass slightly above this maximum value
Comment 2: The location of the global minimizer is unspecified.
'minslope' - controls the globally minimum slope of the function
= [] --> No explicit value provided for the globally minimum
slope of the curve.
= A numeric scalar --> the globally minimum slope of the spline
DEFAULT VALUE: []
Comment: This is a sufficient constraint for the minimum
slope of the spline. It is not a necessary constraint.
There may exist another spline which has a slightly lower
sum of squares and also has the same minimum slope.
'minvalue' - controls the globally minimum value of the curve
= [] --> No explicit value provided for the globally
minimum value of the spline.
= A numeric scalar --> the function will lie no lower
than this minimum value.
DEFAULT VALUE: []
Comment 1: This constraint is only a necessary constraint.
It is not sufficient. In some circumstances the spline may
pass slightly below this minimum value
Comment 2: The location of the global minimizer is unspecified.
'negativeinflection' - controls the existence and placement of a
inflection point of the final curve. The second derivative
of the function will pass through zero at this point.
= [] --> no inflection point constraint employed
= a numeric scalar --> the function will have a point of
inflection at the supplied x location
DEFAULT VALUE: []
Comment: NegativeInflection is really just a composite property
of a curve. A negative inflection point at x == a is equivalent
to a ConcaveUP function for x<=a, and a ConcaveDown function
for x>=a. As such, NegativeInflection will override any other
curvature properties one has previously set.
Comment: Inflection points only apply to linear or cubic models
'order' - An implicit synonym for 'degree', controls the degree of
the piecewise Hermite function. Order MUST be an numeric
integer, from the set [1, 2, 4].
= 1 --> Use a piecewise constant "Hermite"
= 2 --> Use a piecewise linear Hermite
= 4 --> Use a piecewise cubic Hermite
DEFAULT: 4
Setting the order to some value has the effect of setting
the degree of the Hermite spline fit to one less than the
order.
Comment: Some properties are inappropriate for all degrees
of function. E.g., it would be silly to specify a specific
value for the left hand end point slope of a piecewise
constant function. All information supplied will be used
to whatever extent possible.
'plot' - controls whether a final plot is generated of the curve
= 'off' --> No plot
= 'on' --> plot the curve and data using slmplot
DEFAULT VALUE: 'off'
'positiveinflection' - controls the existence and placement of a
inflection point of the final curve. The second derivative
of the function will pass through zero at this point.
= [] --> no inflection point constraint employed
= a numeric scalar --> the function will have a point of
inflection at the supplied x location
DEFAULT VALUE: []
Comment: PositiveInflection is really just a composite property
of a curve. An inflection point at x == a is equivalent to a
ConcaveDown function for x<=a, and a ConcaveUp function for
x>=a. As such, PositiveInflection will override any other
curvature properties one has previously set.
Comment: Inflection points only apply to linear or cubic models
'predictions' - The number of points to evaluate the curve at. Assumed
to be equally spaced. The supplied value must be a positive
scalar, integer value, >= 2, or empty. If empty, no predictions
will be generated.
DEFAULT VALUE: []
Comment: plotslm generates 1001 points along the curve.
'regularization' -
= [] --> Uses the default regularization parameter of 0.0001.
= A Non-negative scalar value --> explicitly defines the weight
given to smoothness of the resulting curve.
= 'crossvalidation' --> Use cross validation method to choose
the regularization parameter.
= A NEGATIVE scalar value --> attempts to choose a
regularization parameter which has as its rmse the absolute
value of the supplied value.
= 'smoothest' --> Finds the smoothest curve that satisfies
the supplied prescription. This is not really a least
squares model, since the goal is purely to maximize the
smoothness. This option would very often be used in
conjunction with errorbars.
DEFAULT VALUE: 0.0001
Comment: Smaller values will yield less smoothing, larger
values more smoothing. In most cases this parameter should
be left alone. It is used to prevent numerical singularities
in the linear algebra, as well as help in the case of
extrapolation and intrapolation. Smoothness of the resulting
spline can be far better controlled by changing the number
of knots and their placement. Specifically, the regularization
parameter is a scale factor applied to the integral of the
squared second derivative of the spline.
Comment: It is possible that no value for the regularization
parameter will yield the given rmse. In this case the sky will
fall down.
Comment: Since the cross validation option and matching a
given rmse are both optimizations, it is not allowed to use
these options together with the interiorknots estimation.
Comment: Both the cross validation and rmse options can
both be quite slow. Remember that they are optimizations.
'result' - controls the output structure style
= 'pp' --> Returns a pp struct, use ppval to evaluate
= 'slm' --> Returns a slm struct in a Hermite form. Evaluate
using slmeval.
DEFAULT VALUE: 'slm'
'rightmaxslope' - controls the maximum slope allowed at the right
hand end of the curve.
= [] --> No explicit value provided for the maximum slope
of the spline at its right hand end point.
= A numeric scalar --> the slope of the function will be
constrained to not rise above this value at its right
hand end point (i.e., the last knot.)
DEFAULT VALUE: []
'rightmaxvalue' - controls the maximum valued allowed at the right
hand end of the curve.
= [] --> No explicit value provided for the maximum value
of the spline at its right hand end point.
= A numeric scalar --> the function will be constrained
to not rise above this value at its right hand end point
(i.e., the first knot.)
DEFAULT VALUE: []
'rightminslope' - controls the minimum slope allowed at the right
hand end of the curve.
= [] --> No explicit value provided for the minimum slope
of the spline at its right hand end point.
= A numeric scalar --> the slope of the function will be
constrained to not fall below this value at its right
hand end point (i.e., the last knot.)
DEFAULT VALUE: []
'rightminvalue' - controls the minimum valued allowed at the right
hand end of the curve.
= [] --> No explicit value provided for the minimum value
of the spline at its right hand end point.
= A numeric scalar --> the function will be constrained
to not fall below this value at its right hand end point
(i.e., the last knot.)
DEFAULT VALUE: []
'rightslope' - controls the function slope at the right hand endpoint.
= [] --> No explicit value provided for the slope of the
curve at its right hand end point.
= A numeric scalar --> the function will be assigned this
slope at its right hand end point (i.e., the last knot.)
DEFAULT VALUE: []
'rightvalue' - controls the function value at its right hand endpoint.
= [] --> No explicit value provided for the value of the
curve at its right hand end point.
= A numeric scalar --> the function will be assigned this
value at its right hand end point (i.e., the last knot.)
DEFAULT VALUE: []
'robust' - Controls the use of a simple robust solver, employing
iteratively reweighted least squares.
= 'off' --> Simple ordinary least squares is employed
= 'on' --> Performs iteratively re-weighted least squares
DEFAULT VALUE: 'off'
Comment: when weights are also supplied, iterative
re-weighting is applied on top of those weights.
Comment: Beware use of the robust fitting option combined
with other options such as:
free interior knots
an rmse goal
cross-validation as a choice of a regularization parameter
In combination with the robust fitting option, any of those
options may introduce instabilities in the result and will at
best be slow.
'scaling' - controls data scaling to avoid numerical problems
= 'on' --> data is shifted and scaled so as to minimize
any numerical problems that may result in the solution.
= 'off' --> No scaling is done.
DEFAULT VALUE: 'on'
Comments: There is no transformation that will positively
eliminate all problems, but this transformation will
somtimes make the solution a bit kess prone to problems.
All scalings are undone in the final spline coefficients,
so these parameters are completely transparent to the user.
In most cases you would see no difference in the solution,
with or without scaling.
Any scaling is automatically applied to both the data as
well as all provided prescription parameters.
The actual shift and scale parameters applied will be reported
as part of the stats.
'segmentconstant' - Allows the user to specify the index of a constant
segment or segments. This is for the rare case where the knots
may be free, yet there is a known internal constant region, so
perhaps one always needs the second polynomial segment to be
exactly constant even though the knots vary in their placement.
The index of the knot interval which is to be held fixed is
supplied. More than one interval may be so held fixed. If there
are N knots, then the supplied index must be an integer between 1
and N-1.
Comment: This property is meaningless for piecewise constant
splines.
Comment: Cubic splines with constant segments will still have
the overall degree of continuity. This may result in strange
looking curves if there are several segments held constant.
DEFAULT Value: []
'segmentlinear' - Allows the user to specify the index of a linear
segment or segments. This is for the rare case where the knots
may be free, yet there is a known internal linear region.
The index of the knot interval which is to be held fixed is
supplied. More than one interval may be so held fixed. If there
are N knots, then the supplied index must be an integer between 1
and N-1.
This property is meaningless for piecewise constant or linear
splines.
DEFAULT Value: []
'simplepeak' - controls the existence and placement of a single
maximizer of the final curve
= [] --> no peak placement constraint employed
= a numeric scalar --> the function will attain its maximum
at the supplied x location
DEFAULT VALUE: []
Comment: SimplePeak is really just a composite property of
a curve. A peak at x == a is equivalent to a monotone increasing
function for x<=a, and a monotone decreasing function for
x>=a. As such, simplepeak will override any other monotonicity
properties one has previously set.
'simplevalley' - controls the existence and placement of a single
minimizer of the final curve
= [] --> no valley placement constraint employed
= a numeric scalar --> the function will attain its minimum
at the supplied x location
DEFAULT VALUE: []
Comment: SimpleValley is really just a composite property
of a curve. A valley at x == a is equivalent to a monotone
decreasing function for x<=a, and a monotone increasing
function for x>=a. As such, simplevalley will override any
other monotonicity properties one has previously set.
'sumresiduals' - Allows the user to specify an explicit goal for the
sum of the residuals. This is a property that virtually
NOBODY should ever have a need for, since normally you get
a zero sum implicitly, as a freebie. This is due to the
presence of an effective constant term in the regression
model. Under some circumstances however, that property is
circumvented, for example by forcing the model to explicitly
pass through a given point. Use of non-unit weights in the
model could also cause a non-zero sum of residuals.
In practice, any real scalar value for the sum of the
residuals can be specified, but a value of 0 is the only
value that makes any sense that I can see. Weights on the
data points are disregarded when the residual sum is
computed.
The computation of residual here is defined as [yhat - y].
DEFAULT VALUE: []
'verbosity' - controls commandline feedback to the user
= 0 --> No response at the commandline
= 1 --> Basic fit statistics reported at the command line
= 2 --> Debug level output
DEFAULT VALUE: 0
'weights' - defines a weight vector
= [] --> all data points are assigned equal (unit) weight.
= vector of the same length as length(x), denotes relative
weights for each data point. If supplied, the length of
this vector of weights must be the same as the number of
data points.
DEFAULT VALUE: []
'xy' - Forces the curve through an individual point or set of points
= [] --> no points are forced
= a 1x2 vector --> an x-y pair that the curve must pass
through with no error.
= an nx2 array --> each row corresponds to a single point
that the curve passes through.
DEFAULT VALUE: []
Comment 1: The curve will pass through the desired point
to within computational error, IF it is possible to do so.
Comment 2: Multiple points that are inconsistent with each
other, or inconsistent with other parameters that are set
will cause failure of the least squares.
'xyp' - Forces the curve to have a specified slope at an individual
point or set of points
= [] --> no points have their slope enforced
= a 1x2 vector --> an x-yprime pair that the curve must
satisfy through with no error.
= an nx2 array --> each row corresponds to an x-yprime
pair that the curve must satisfy.
DEFAULT VALUE: []
Comment 1: The curve will satisfy the desired slope to
within computational error, IF it is possible to do so.
Comment 2: Multiple points that are inconsistent with each
other, or inconsistent with other parameters that are set
will cause failure of the least squares.
Comment 3: Setting the slope at some point of a zero'th
degree function will be ignored.
'xypp' - Forces the curve to have a specified second derivative
at an individual point or set of points
= [] --> no points have their 2nd derivative enforced
= a 1x2 vector --> an x-y'' pair that the curve must
satisfy through with no error.
= an nx2 array --> each row corresponds to an x-y''
pair that the curve must satisfy.
DEFAULT VALUE: []
Comment: The curve will satisfy the desired 2nd derivative
to within computational error, IF it is possible to do so.
Comment: Multiple points that are inconsistent with each
other, or inconsistent with other parameters that are set
will cause failure of the least squares.
Comment: Setting the 2nd derivative at some point of a
zero'th degree or linear function will be ignored.
Comment: This property only applies to cubic models.
'xyppp' - Forces the curve to have a specified third derivative
at an individual point or set of points.
Really, the only meaningful use of this property that I
see under normal circumstances is to force a segment
to be only quadratic, rather than cubic. Thus, by forcing
the third derivative to zero at some point in the segment,
since the third derivative of a cubic polynomial is a
constant, will reduce the cubic to a quadratic.
= [] --> no points have their 3rd derivative enforced
= a 1x2 vector --> an x-y''' pair that the curve must
satisfy through with no error.
= an nx2 array --> each row corresponds to an x-y'''
pair that the curve must satisfy.
DEFAULT VALUE: []
Comment: Setting the 3rd derivative at some point of a
zero'th degree or linear function will be ignored.
Comment: This property only applies to cubic models.
Evaluating a model
One more thing before we return to fitting data. We can use the models that we build in a predictive mode.
slmeval(.5,slm)
ans =
1.5204
slmeval can also differentiate the function (although piecewise constant models would have a zero derivative everywhere, with a singularity at every knot.) The third argument controls the order ofthe derivative to be generated.
slmeval(.5,slm,1)
ans =
2.0316
and even compute an inverse at any point. Note that slmeval computes only the "left-most" inverse. So if there is more than one solution that you wish to find, then you should use slmsolve.
slmeval(1.50,slm,-1)
ans =
0.48996
We could plot the function as fit uing slmeval,
yhat = slmeval(xev,slm,0); plot(x2,y2,'ro',xev,yhat,'b-')
Plotting a model using plotslm
Or we could have plotted it using a nice utility that I've provided. Plotslm plots the curve, the data that it was built from, and the knots/breaks of the spline model. Plotslm also allows you to plot the residuals of the model to the data, or the derivatives of the curve. Finally, you can control many features of this plot using the plot menu.
plotslm(slm)
slmengine can plot too
Or we could have just told our engine to plot after it finished estimating its curve fit. This time we will fit it as a cubic spline.
close all slm = slmengine(x2,y2,'degree',3,'knots',4,'plot','on');
slmpar: Computations on a spline model, min max, integral, etc.
plotslm can plot the function derivatives or integral, and slmeval can evaluate the function, its derivatives or the inverse function at specific points. However, it is occasionally useful to be able to find the minimum (or maximum) value of a spline over some interval, or find the minimum or maximum value of the slope, or compute the integral of the curve between a pair of points. slmpar also can convert a spline model into a symbolic form, for those who wish to see what the segments are in a functional form.
The slmpar function does all of these things, both for a slm form as well as a pp form.
x = -5:.1:5; y = exp(-x.^2/2)/sqrt(2*pi) + randn(size(x))/100; slm = slmengine(x,y,'plot','on','knots',15,'integral',1,'deg',3,'minval',0)
slm =
form: 'slm'
degree: 3
knots: [15x1 double]
coef: [15x2 double]
stats: [1x1 struct]
prescription: [1x1 struct]
x: [101x1 double]
y: [101x1 double]
Extrapolation: 'constant'
Of course, here the integral was forced to 1, so the integral should be correct.
res = slmpar(slm,'in')
res =
1
The maximum value should occur very near zero.
[res,loc] = slmpar(slm,'maxfun') % The next few examples show slmpar at work on a pp form, as created by spline: x = linspace(-pi,pi); y = sin(x); pp = spline(x,y);
res =
0.40277
loc =
-0.0054187
The overall integral should be zero, or as close as floating point computations can give us for the integral of a spline.
res = slmpar(pp,'integral')
res = -1.9993e-16
The integral from 0 to pi should be 2.
res = slmpar(pp,'integral',[0,pi])
res =
2
The maximum of a sine wave over its period is 1, and should occur at x = pi/2.
[res,loc] = slmpar(pp,'maxfun')
res =
1
loc =
1.5708
The minimum of a sine wave over [-1,1] should occur at the lower end point.
[res,loc] = slmpar(pp,'minfun',[-1,1])
res =
-0.84147
loc =
-1
A nice thing about plotslm is it overlays the knot positions on top of the curve. The knots are at the vertical green dotted lines.
slmsolve: find all roots of a spline, f(x) == y, for any given value of y
slmsolve will find all solutions for the this general problem, as opposed to slmeval, which returns only the left-most (closest to -inf) solution.
slmsolve works on any SLM toolbox produced spline model, as well as any pp form model, produced by pp, spline, etc.
x = linspace(0,1,200); y = x + sin(20*x); pp = spline(x,y); fnplt(pp)
Find the solutions of f(x) = 0.5 over the support of the spline.
xloc = slmsolve(pp,0.5) hold on plot(xloc',0.5,'ro',[min(xloc);max(xloc)],[0.5 0.5],'r-')
xloc =
Columns 1 through 6
0.024762 0.13859 0.32305 0.46972 0.62219 0.80067
Column 7
0.92076
Model fitting using shape prescriptors
Having seen the exponential data fit above, one might wonder why one would use a shape prescriptive model at all. It appears to be little better than that 6'th order polynomial model we fit before using polyfit.
What do we know about the underlying functional relationship? We want to think in terms of fundamental shape primitives.
I'd suggest that the function is known to be a monotone increasing function. It was a simple exponential after all. The first derivatve should never be negative.
close all slm = slmengine(x2,y2,'plot','on','increasing','on');
Curvature constraints
I'd also suggest that the function is known to be positively curved. The second derivatve should never be negative.
close all slm = slmengine(x2,y2,'plot','on','concaveup','on');
Both curvature and monotonicity constraints
We can also use both a curvature and monotonicity constraint, although on many data sets, one will be sufficient.
close all slm = slmengine(x2,y2,'plot','on','concaveup','on','increasing','on');
Minimum slope or maximum slope constraints
Try an erf function, with just a bit of noise. But then, show that the curve can be constrained to have an overall minimum and maximum slope.
n = 100; x = 4*rand(n,1) - 2; y = erf(x) + randn(size(x))/100;
First, with no slope constraints at all
slm = slmengine(x,y,'plot','on','knots',15);
Plot the first derivative. See that it should strictly be bounded by in the interval [0,1].
xev = linspace(-2,2,1000); d = slmeval(xev,slm,1); plot(xev,d,'-') grid on
Redo the fit, but with just a touch of a constraint on the slopes.
slm = slmengine(x,y,'plot','on','knots',15,'minslope',0.1,'maxslope',0.9); d = slmeval(xev,slm,1); figure plot(xev,d,'-') grid on
Extrapolation in the fit, non-uniform knot spacings
We might also decide to fit the curve over a wider knot range, extrapolating here beyond the data. Note that it is BETTER to build a spline to have the shape you want over a desired region then to try to extrapolate a spline that was built over a smaller region.
In fact, I try to force the user to do no extrpolation at all, unless they insist on doing so. The extrapolation options will be shown in some detail later on.
close all slm = slmengine(x2,y2,'plot','on','concaveup','on','increasing','on','knots',[-3, -2:.5:1 , 2]);
Piecewise constant models
Or even to use a lower order model. This version of my toolkit only allows piecewise constant or piecewise linear besides the cubic models.
Also, for brevity, we can always shorten any property names, as long as the name remains unambiguous. Capitalization is ignored.
slm = slmengine(x1,y1,'plot','on','kn',0:.1:1,'deg',0);
Note that I also used more knots in this last fit. I also returned to our original set of data with 100 points.
Fixed specific values
There are other shape parameters we could have set in our prescription. Kowing that the underlying relationship is an exponential, we might remember that at x == 0, that y == 1. Specific values like this are often known. If our data was from an MTF curve, we might choose to assume a modulation of 100% at a frequency of 0. Likewise, population growth data might have a known population at time == 0.
slm = slmengine(x1,y1,'plot','on','concaveup','on','knots',0:.2:1,'leftvalue',1);
Flawed information, "Garbage in, garbage out"
What if our information is flawed? The old saying always applies.
slm = slmengine(x1,y1,'plot','on','concaveup','on','knots',0:.2:1,'leftvalue',2);
Or suppose we specify a bit too small a value for the maximum slope. The result is lack of fit.
slm = slmengine(x1,y1,'plot','on','maxslope',1.5,'knots',0:.1:1);
Garbage in, garbage out. If we specify a monotone decreasing function on increasing data, then we have wasted a lot of machine cycles to compute the mean of our data.
slm = slmengine(x1,y1,'plot','on','decreasing','on','knots',0:.2:1,'degree',1);
We could have done the same with somewhat less work...
slm = slmengine(x1,y1,'plot','on','degree',0,'knots',2);
The really lazy would just use mean(y1). Perhaps efficient is a better word.
lets look at a function with slightly more interesting behaviour. still keep it simple enough that we know what it should look like.
close all n = 500; x = linspace(0,1,n)'; y = sin(x*pi) + randn(n,1)/5; slm = slmengine(x,y,'plot','on','knots',10);
If I remember correctly, sin(x) will have a negative second derivative over that interval. If I choose to apply this information to our curve, the fit gets quite good.
slm = slmengine(x,y,'plot','on','knots',10,'concavedown','on');
A simple way of specifying a curve shape with a single maximum value at some point is just this:
slm = slmengine(x,y,'plot','on','knots',0:.1:1,'simplepeak',.5);
all it does is to automatically specify increasing and decreasing regions of the curve.
lets change our data again now. Make it a fairly noisy sine wave.
clear
close all
n = 250;
x = linspace(0,1,n)';
y = sin(2*x*pi) + randn(n,1)/2;
this time we might choose to assume periodicity.
slm = slmengine(x,y,'plot','on','knots',0:.1:1,'endconditions','periodic');
That last curve had a lot of noise in it, I might not be surpised if it was a little bumpy. We also know that somewhere in the middle of the curve, this function goes through an inflection point. I'll just pick somehere to put that inflection point. Note that the curvature would be negative to the left, and positive to the right, so set a positive inflection point.
slm = slmengine(x,y,'plot','on','kn',0:.1:1,'endc','per','positiveinflection',.5);
There are certainly many ways to define a prescription for the shape of a curve. We have only covered a few of them here.
Lets return to slmset though. Note that the two lines below are the same as the direct call to slmengine.
prescription = slmset('plot','on','kn',0:.1:1,'endc','per','positiveinflection',.5); slm = slmengine(x,y,prescription);
The model as it is returned includes a field with the shape prescription that defined that model. It also includes the original data in named fields. Each of these fields is useful for documentation purposes.
slm
slm =
form: 'slm'
degree: 3
knots: [11x1 double]
coef: [11x2 double]
stats: [1x1 struct]
prescription: [1x1 struct]
x: [250x1 double]
y: [250x1 double]
Extrapolation: 'constant'
We can also modify an existing prescription, adding new pieces of information.
slm = slmengine(x,y,prescription,'maxvalue',0.5);
Too many knots?
Sometimes we have little data, but use too many knots. A traditional spline model will often fail in this case, with a singular matrix warning. Slmengine will not fail, although you may not always be happy with the result. Remember, when there is little information brought to a problem by the data, a reasonable fit may still be obtained if you can add enough information from your own kowledge of the process. In this particular example, I'm tryin only to show that the use of too many kots did not result in a singularity.
clear close all n = 10; x = rand(n,1); y = sin(2*x*pi) + randn(n,1)/10; slm = slmengine(x,y,'plot','on','knots',0:.05:1);
Fewer knots may be better, even though we did survive having way too many knots.
slm = slmengine(x,y,'plot','on','knots',[0 1]);
And while it may be a little slow, one can use cross validation to help smooth your curve. Don't try it with too many knots or too many points.
slm = slmengine(x,y,'plot','on','knots',0:.05:1,'regularization','cross');
Another option, if you had a decent estimate of the error standard deviation, is to supply that value. slmengine will build a fit with the desired residual error.
slm = slmengine(x,y,'plot','on','knots',0:.05:1,'regularization',-0.10);
Use in practice
for some more realistic data, try this one out. The function is known to be increasing, and everywhere non-negative. There is a little curvature at the bottom that is probably real. At the top end, we expect the curve to be moderately well behaved so extrapolation to x = 255 is linear.
close all load scann slm = slmengine(R,L,'plot','on','knots', ... [0 5 10 15 20 30 50 80 120 160 200 255],'leftminvalue',0,'increasing','on','deg',3,'linearregion',[200,255])
slm =
form: 'slm'
degree: 3
knots: [12x1 double]
coef: [12x2 double]
stats: [1x1 struct]
prescription: [1x1 struct]
x: [42x1 double]
y: [42x1 double]
Extrapolation: 'constant'
I often see people using splines wrongly. For example, given a relationship with a derivative singularity, how best to build a spline model?
x = linspace(0,1,20); truefun = @(x) nthroot(x,5); y = truefun(x); xfine = linspace(0,1,1000); yfine = truefun(xfine); plot(x,y,'-o') title 'A singularity at x = 0' xlabel X ylabel Y
Obviously, fitting a spline to this data, in the form y(x), we should expect a problem. Even with no noise present in the system.
slm = slmengine(x,y,'plot','on','knots',10) hold on plot(xfine,yfine,'b--') hold off
slm =
form: 'slm'
degree: 3
knots: [10x1 double]
coef: [10x2 double]
stats: [1x1 struct]
prescription: [1x1 struct]
x: [20x1 double]
y: [20x1 double]
Extrapolation: 'constant'
We might force the spline to be monotone, or even choose a non-uniform knot spacing, but the derivative singularity at x == 0 is too much for polynomial segments to handle well.
slm = slmengine(x,y,'plot','on','knots',[0 .01 .02 .05 .1 .3 .6 1],'incr','on','concavedown','on') hold on plot(xfine,yfine,'b--') hold off
slm =
form: 'slm'
degree: 3
knots: [8x1 double]
coef: [8x2 double]
stats: [1x1 struct]
prescription: [1x1 struct]
x: [20x1 double]
y: [20x1 double]
Extrapolation: 'constant'
A better solution is to recognize the essential singularity in the problem. Here, swap the dependent and independent variables. What was once a singularity is no longer so. Se that the spline fit and the original function now overlay each other almost perfectly here.
slm = slmengine(y,x,'plot','on','knots',10,'incr','on','concaveup','on','deg',3) hold on plot(yfine,xfine,'b--') hold off
slm =
form: 'slm'
degree: 3
knots: [10x1 double]
coef: [10x2 double]
stats: [1x1 struct]
prescription: [1x1 struct]
x: [20x1 double]
y: [20x1 double]
Extrapolation: 'constant'
Of course, prediction will be more difficult, since the resulting model must be interpolated in an inverse form. slmevel can do so. The result here would be 0.7, IF the spline were an exact representation of the fifth root of its input. Since we live in a floating point world, the result is as close as I might expect.
format long g slmeval(0.7.^5,slm,-1) format short g
ans =
0.700006517591832
Returning a set of points, evaluated through the fitted curve.
SLM has the ability to generate a list of points, then evaluate them through the resulting model. While this is easily enough done using slmeval, this is a friendly option that will save a set for the user. By default, 1001 points are generated along the curve, however, any number of equally spaced points can be generated. For the user who wishes for an unequal (user specified) spacing, there is always slmeval to fall back upon. In this example, only 11 equally spaced points are generated along the curve to be returned. The user can then plot these points as they choose. Of course, plotslm has significantly more functionality than a simple plot, since you can view more than just the simple function values.
close all format long g [slm,xp,yp] = slmengine(y,x,'knots',10,'incr','on','concaveup','on','rightvalue',1,'predictions',11)
slm =
form: 'slm'
degree: 3
knots: [10x1 double]
coef: [10x2 double]
stats: [1x1 struct]
prescription: [1x1 struct]
x: [20x1 double]
y: [20x1 double]
Extrapolation: 'constant'
xp =
Columns 1 through 3
0 0.1 0.2
Columns 4 through 6
0.3 0.4 0.5
Columns 7 through 9
0.6 0.7 0.8
Columns 10 through 11
0.9 1
yp =
Columns 1 through 3
-7.40407735122517e-13 -7.39976019063549e-13 1.92543951631985e-05
Columns 4 through 6
0.000916087696866076 0.00912838378939623 0.0311222564765782
Columns 7 through 9
0.0777495706060542 0.168062444052809 0.327684525849983
Columns 10 through 11
0.590505873560403 1
Combining various constraints to produce a better model
On this curve, we know the actual form. It should be a hyperbolic shape, that approaches a linear asymptote.
x = rand(40,1);
y = x - 1./(10*(x+.1)) + randn(size(x))/25;
plot(x,y,'o')
We can model the curve as simply a monotone increasing function. I'll use more knots than I really need to exacerbate the problem. (There really are way too many knots in this spline, possibly causing numerical problems.)
slm = slmengine(x,y,'incr','on','knots',35,'plot','on')
slm =
form: 'slm'
degree: 3
knots: [35x1 double]
coef: [35x2 double]
stats: [1x1 struct]
prescription: [1x1 struct]
x: [40x1 double]
y: [40x1 double]
Extrapolation: 'constant'
Now, try enforcing both monotonicity, and a concave down function, (a negative second derivative.)
slm = slmengine(x,y,'incr','on','concavedown','on','knots',35,'plot','on')
slm =
form: 'slm'
degree: 3
knots: [35x1 double]
coef: [35x2 double]
stats: [1x1 struct]
prescription: [1x1 struct]
x: [40x1 double]
y: [40x1 double]
Extrapolation: 'constant'
Our second try is better, but still not well enough behaved. See that by requiring that the jerk be positive everywhere, we get a nice, smooth curve, even with far too many knots for this simple curve. In fact, it even appears to be approaching a linear asymptote.
slm = slmengine(x,y,'incr','on','jerk','positive','concavedown','on','knots',35,'plot','on')
slm =
form: 'slm'
degree: 3
knots: [35x1 double]
coef: [35x2 double]
stats: [1x1 struct]
prescription: [1x1 struct]
x: [40x1 double]
y: [40x1 double]
Extrapolation: 'constant'
Extrapolation of a spline model
Normally, I try to push a user to NOT extrapolate a spline model. Cubic polynomial segments do squirrely things away from where they were used in a fit. So if you seriously tried to extrapolate a cubic spline, then don't be surprised to see it do silly things. For this reason the default for SLMEVAL is to only extrapolate as a constant function. Thus below the first knot or above the top knot, the spline prediction will take on the values of the function at the corresponding knot, but no more. Essentially only constant extrapolation is done.
SLMSET does provide for other options in extrapolation, but only if you ask for it. The 'extrapolation' prescription may take on the values:
{'error', 'warning', 'constant', 'linear', 'cubic', 'NaN'}The default is 'constant'. For this case, constant extrapolation is done with with no warning generated for all points outside of the knots. The other options seem pretty stratight forward.
'error' generates an error if any point falls outside the knots.
'warning' generates a warning message, then does constant extrapolation.
'constant' does constant extrapolation. (The default)
'linear' does linear extrapolation.
'cubic' does cubic extrapolation.
'nan' inserts NaN values for all points outside the knots.
close all
x = 0:.1:1;
y = sin(2*pi*x) + randn(size(x))/100;
See that the Extrapolation field has been set to 'constant' by default.
slm = slmengine(x,y,'knots',x,'plot','on')
slm =
form: 'slm'
degree: 3
knots: [11x1 double]
coef: [11x2 double]
stats: [1x1 struct]
prescription: [1x1 struct]
x: [11x1 double]
y: [11x1 double]
Extrapolation: 'constant'
See what happens when various extrapolation modes are specified though. Note that the extrapolation type is indicated when you create the spline itself. I won't show what happens when 'error' is chosen, since that will screw up this published document. You can guess though.
xev = -1:.01:2;
Warning throws a warning if extrapolation is tried, then does constant extrapolation, like the default.
slm = slmengine(x,y,'knots',x,'extrapolation','warning'); yev = slmeval(xev,slm,0); plot(x,y,'o',xev,yev,'r-') title 'Constant extrapolation (warning generated)'
Warning: One or more points fell outside of the knots. constant extrapolation performed
Linear simply does linear extrapolation, with no warnings or errors.
slm = slmengine(x,y,'knots',x,'extrapolation','linear'); yev = slmeval(xev,slm,0); plot(x,y,'o',xev,yev,'r-') title 'Linear extrapolation'
Cubic simply does cubic extrapolation, with no warnings or errors. I would rarely ever recommend 'cubic' extrpolation, certainly not if you extraapolate too far. See what happens here.
slm = slmengine(x,y,'knots',x,'extrapolation','cubic'); yev = slmeval(xev,slm,0); plot(x,y,'o',xev,yev,'r-') title 'Cubic extrapolation'
You can get information about the fit itself
slmengine always returns this information in the stats field of the returned model. However, the verbosity propert has three levels {0,1,2}. 0 tells slmengine to be quiet, with no output to the command line. 1 tells it to report fit statistics at the command line,
close all slm = slmengine(x,y,'verbosity',1);
========================================= lse ========================================= ========================================= MODEL STATISTICS REPORT Number of data points: 11 Scale factor applied to y 0.5212 Shift applied to y 1.1233 Total degrees of freedom: 12 Net degrees of freedom: 8 R-squared: 0.99998 Adjusted R-squared: 0.99994 RMSE: 0.0056764 Range of prediction errors: -0.0056066 0.0049007 Error quartiles (25%, 75%): -0.0022108 0.0021211 =========================================
and verbosity of 2 tells slmengine to be more expressive yet.
slm = slmengine(x,y,'verbosity',1);
========================================= lse ========================================= ========================================= MODEL STATISTICS REPORT Number of data points: 11 Scale factor applied to y 0.5212 Shift applied to y 1.1233 Total degrees of freedom: 12 Net degrees of freedom: 8 R-squared: 0.99998 Adjusted R-squared: 0.99994 RMSE: 0.0056764 Range of prediction errors: -0.0056066 0.0049007 Error quartiles (25%, 75%): -0.0022108 0.0021211 =========================================
Knot placement optimization
You can let an optimizer choose the knots. It may give you better results than a blind choice of knots. But, the fact is all optimizers can get lost at times.
close all x = -5:.1:5; y = erf(x) + randn(size(x))/100; slm = slmengine(x,y,'increasing','on','plot','on','degree',1, ... 'interiorknots','free','knots',10)
slm =
form: 'slm'
degree: 1
knots: [10x1 double]
coef: [10x1 double]
stats: [1x1 struct]
prescription: [1x1 struct]
x: [101x1 double]
y: [101x1 double]
Extrapolation: 'constant'
If however, you have some idea where the knots might belong based on the shape of your function, then it is a good idea to provide some input. Here for example, there is little happening with this function on both ends of the curve. So why waste the knots out there?
slm = slmengine(x,y,'increasing','on','plot','on','degree',1,... 'interiorknots','free','knots',[-5 -1.25 -1 -.75 -.5 .5 .75 1 1.25 5])
slm =
form: 'slm'
degree: 1
knots: [10x1 double]
coef: [10x1 double]
stats: [1x1 struct]
prescription: [1x1 struct]
x: [101x1 double]
y: [101x1 double]
Extrapolation: 'constant'
slmpar has some other options, that are not terribly useful for computation, but are nice to understand how splines work, as well as how the coefficients are stored in these spline containers. These last options in slmpar require the symbolic toolbox.
First, I'll create a piecewise linear spline with 4 segments.
x = linspace(-pi,pi,100); y = sin(x); slm = slmengine(x,y,'knots',5,'result','slm','degree',1)
slm =
form: 'slm'
degree: 1
knots: [5x1 double]
coef: [5x1 double]
stats: [1x1 struct]
prescription: [1x1 struct]
x: [100x1 double]
y: [100x1 double]
Extrapolation: 'constant'
The SLM storage form can be called a Hermite form. It contains a list of knots, and the function values of the curve at those knots. So we see below the first column isthe list of knots, while the second column the values of the sine function approximation at those knots.
[slm.knots,slm.coef]
ans =
-3.14159265358979 -0.0917040748632527
-1.5707963267949 -1.19301540478917
0 -1.33209993392464e-15
1.5707963267949 1.19301540478917
3.14159265358979 0.0917040748632509
However, we can see the actual linear polynoial segments using slmpar. This option tells slmpar to create a cell array that lists the knot intervals over which each segment is defined, as well as the actual linear polynomial used over that interval, as a symbolic polynomial in the variable x.
p = slmpar(slm,'symabs')
p{:}
p =
[1x2 double] [1x2 double] [1x2 double] [1x2 double]
[1x1 sym ] [1x1 sym ] [1x1 sym ] [1x1 sym ]
ans =
-3.14159265358979 -1.5707963267949
ans =
- 0.70111656816327772023811348844902*x - 2.2943267347150934757856043826599
ans =
-1.5707963267949 0
ans =
0.75949719542790039561452886118786*x - 0.0000000000000013139637232895128895315501627765
ans =
0 1.5707963267949
ans =
0.75949719542790206094906579892267*x - 0.0000000000000013320999339246355654677894198042
ans =
1.5707963267949 3.14159265358979
ans =
2.2943267347150951509722696831747 - 0.70111656816327883046113811360556*x
We can use slmpar to understand the pp form of a spline, and how the spline coefficients are stored. They are stored as the coefficients of a polynomial in standard matlab poly form, so highest order coefficient first. However, the pp form uses what I'll call a relative polynomial.
pp = slm2pp(slm)
pp =
form: 'pp'
breaks: [1x5 double]
coefs: [4x2 double]
pieces: 4
order: 2
dim: 1
prescription: [1x1 struct]
We can list out the breaks (or knots if you choose) of the spline. between each pair of breaks, there is here a linear polynomial segment defined. However, the i'th such polynomial will be evaluated as if it lives in the interval [0,pp.breaks(i+1) - pp.breaks(i)]
pp.breaks pp.coefs
ans =
Columns 1 through 3
-3.14159265358979 -1.5707963267949 0
Columns 4 through 5
1.5707963267949 3.14159265358979
ans =
-0.701116568163278 -0.0917040748632527
0.7594971954279 -1.19301540478917
0.759497195427902 -1.33209993392464e-15
-0.701116568163279 1.19301540478917
So the first linear polynomial segment actually lives on the interval
[pp.breaks(1),pp.breaks(2)]
ans =
-3.14159265358979 -1.5707963267949
However, it is evaluated on the translated interval
[0,pp.breaks(2)-pp.breaks(1)]
ans =
0 1.5707963267949
This translation of the interval is not terribly important when one is working with linear spline segments, but when cubic (or even higher order) segments are employed, there may well be numerical issues if MATLAB is forced to take powers of very large numbers. Suppose for example, our spline was forced to be evaluated for numbers on the order of 1e10? Now forming the cube of such large numbers will cause massive loss of accuracy in the floating point numbers.
Similarly, we might create a cubic spline. Here I'll use only 4 knots, thus three cubic segments, but I'll generate data that lies in a different domain.
x = linspace(1e6,1e6 + 3,100); y = sin(x); slm = slmengine(x,y,'knots',4,'result','slm','degree',3,'plot','on')
slm =
form: 'slm'
degree: 3
knots: [4x1 double]
coef: [4x2 double]
stats: [1x1 struct]
prescription: [1x1 struct]
x: [100x1 double]
y: [100x1 double]
Extrapolation: 'constant'
First, lets understand the Hermite form here. The first column below is the list of knots. The second column the function values at those knots, and the third column contains the first derivatives of the spline function at those knots. Personally, I've always had a preference for this form because one can easily visualize much about the shape of the cubic segments directly from those coefficients.
[slm.knots,slm.coef]
ans =
1000000 -0.348818178910341 0.923373690093335
1000001 0.600371047473741 0.797449369267934
1000002 0.998889373328152 -0.0700485104495925
1000003 0.478882305961742 -0.881721552005563
Converting the spline to its pp form, we get:
pp = slm2pp(slm)
pp =
form: 'pp'
breaks: [1000000 1000001 1000002 1000003]
coefs: [3x4 double]
pieces: 3
order: 4
dim: 1
prescription: [1x1 struct]
There are 4 breaks, so 3 segments. These coefficients presume the relative form for evaluation
pp.coefs
crel = slmpar(pp,'symrel')
crel{:}
ans =
Columns 1 through 3
-0.177555393406894 0.203370929697641 0.923373690093335
-0.0696357928904805 -0.329295250523043 0.797449369267934
0.0882440722776654 -0.538202629194483 -0.0700485104495925
Column 4
-0.348818178910341
0.600371047473741
0.998889373328152
crel =
[1x2 double] [1x2 double] [1x2 double]
[1x1 sym ] [1x1 sym ] [1x1 sym ]
ans =
1000000 1000001
ans =
- 0.17755539340689430183317654154962*u^3 + 0.20337092969764092309503666911041*u^2 + 0.92337369009333514213722082786262*u - 0.34881817891034055945098657502967
ans =
1000001 1000002
ans =
- 0.069635792890480510686224135952216*u^3 - 0.32929525052304264853830773063237*u^2 + 0.79744936926793408282776454143459*u + 0.60037104747374125945924561165157
ans =
1000002 1000003
ans =
0.088244072277665397407986347388942*u^3 - 0.53820262919448347282980193995172*u^2 - 0.070048510449592496507342787026573*u + 0.99888937332815230796256855683168
See however, what happens to the coefficients of the polynomials in the absolute form. Those coefficients are now a bit nasty for a cubic polynomial. If x is on the order of 1e6, I would expect to see a serious loss of accuracy.
cabs = slmpar(pp,'symabs')
cabs{:}
cabs =
[1x2 double] [1x2 double] [1x2 double]
[1x1 sym ] [1x1 sym ] [1x1 sym ]
ans =
1000000 1000001
ans =
- 0.17755539340689430183317654154962*x^3 + 532666.38359161260314045271968553*x^2 - 532666586961.61892709128247958006*x + 177555596776900625.43518812253381
ans =
1000001 1000002
ans =
- 0.069635792890480510686224135952216*x^3 + 208907.25828356968045755592822132*x^2 - 208907137894.45069636479687073816*x + 69635672501361526.136727299042292
ans =
1000002 1000003
ans =
0.088244072277665397407986347388942*x^3 - 264733.28450005905269981631988685*x^2 + 264734352170.26360404933392377058*x - 88245139947869950.477280652044594
To compare the accuracy problems one would expect, what is the sin of a moderately large number?
sin(sym('1e6 + 0.5'))
sin(1e6 + 0.5)
slmeval(1e6 + 0.5,slm)
ans =
0.14195469900074400352584942400447
ans =
0.141954699000744
ans =
0.141516974384875
So the 4 knot spline fit itself was only accurate to about 3 significant digits. Lets now use polyval to evaluate those polynomials at that same effective location. The relative form returns a value that matches slmeval exactly, as well as being a viable value for the sine function. However, the absolute form just returns numerical junk.
polyval(fliplr(double(coeffs(crel{2,1}))),0.5)
polyval(fliplr(double(coeffs(cabs{2,1}))),1e6 + 0.5)
close all
ans =
0.141516974384875
ans =
-32