I'd like to introduce a new guest blogger - John D'Errico - an applied mathematician, now retired from Eastman Kodak, where he used MATLAB for over 20 years. Since then, MATLAB is still in his blood, so you will often find him answering questions on the newsgroup and writing new utilities to add to MATLAB Central.
Contents
I'll assume you have some data points through which you wish to pass a curve, interpolating your data. (Initially, I will only talk about problems with one independent variable.) In these coming blogs, I'll try to show some ways to do exactly this, i.e., find a curve that passes through your data. Along the way I'll try to give some pointers on curve fitting, interpolation, modeling, approximation, etc.
Polynomial regression
A valid question for some to ask is why start out with a discussion about polynomial regression , when we really wanted to talk about interpolation. Many people mistake the ideas of interpolation with the approximation produced by a regression model, calling both of these things interpolation. So I'm starting out with some discussion about what interpolation is not. Plus, I want to assure an understanding of polynomials, since many of the tools for interpolation are polynomial based in some way.
Let us start by creating some data. An exponential is a good place to start, a simple curve shape that is easy to fit.
x = -1:.1:1; y = exp(x);
Plot your data
It is always a good idea to plot your data. In fact, I'll suggest that you should plot everything. Plots are useful, since your eye and brain are splendid at things like pattern recognition. Only you know your data, as the scientist, engineer, or analyst. You will always benefit if you can employ your knowledge of a system as part of the modeling process.
Next, always think in advance about your goals for any model.
- Will you use this model purely for interpolation, i.e., for predictive purposes only?
- Do you need to derive some understanding about the process from your model? Perhaps you need to estimate an upper asymptote of your process. If so, then you may want a model that has such an upper asymptote built into it.
- Will you need to include the model coefficients in some paper that you expect to write? Will you need to use this model in MATLAB, or in some other tool?
- Must the model be simple to evaluate?
- Must the model be efficient to evaluate? You should know whether you will need to evaluate this model millions of times, or just once.
- What do you know about your data? Is there noise in the data? May you ignore that noise, or must you smooth the noise away?
- What do you know about the underlying functional relationship? Is it monotone? Increasing? Decreasing? Positively/negatively curved? Must it pass through a specific point?
- Must the interpolant have specific requirements in terms of continuity or differentiability?
For example, I once had a problem where I knew I had some significant noise in our process, but I chose not to do any smoothing anyway. Any such smoothing would also have smoothed out some potentially important features of our process. Since I could survive with the noise in my interpolant, I chose the lesser evil.
Always know your goals for any such task.
plot(x,y,'bo') xlabel X ylabel Y grid on title 'Exponential data'
This is a nice, well-behaved function. It is of the form y = f(x). I'll pretend for the moment that I have no idea what was in the underlying functional relationship.
One thing I learned in some long past calculus course was that a Taylor series will provide an approximation for many functions. Polynomials are useful things. They are simple to use, simple to build, simple to work with. And a truncated Taylor series is basically a nice, simple polynomial.
So we will start with a linear polynomial approximation for this curve, built using polyfit . This is a utility provided in MATLAB to estimate a polynomial model using linear regression techniques. We could also use many other tools to build our polynomial model, but polyfit is a useful one, and easy to use.
A linear, or first degree polynomial (many use the words "order" and "degree" interchangeably), might be written mathematically as y(x) = a1*x + a2. In MATLAB we will merely store the coefficients, as a vector [a1,a0]. Note that a polynomial in MATLAB has it's coefficients stored with the highest order term first.
P1 = polyfit(x,y,1)
P1 =
1.1140 1.1937
We can evaluate the polynomial with polyval.
yhat = polyval(P1,x); plot(x,y,'bo',x,yhat,'r-') xlabel X ylabel Y grid on title 'Linear polynomial fit'
Look at the residuals
In fact, I'll claim the relationship we are modeling is not terribly well represented by a linear model. Depending on your needs for this model, you might have decided differently.
When you build a regression model, look at the residuals.
ALWAYS PLOT EVERYTHING!
Plot your residuals. In general, some good ways to plot the residuals are versus
- the independent variable. Look for patterns. Patterns here in the residuals are often a hint that you should look more deeply.
- the measurement order, just in case there was a problem with your equipment. (I've seen many cases where measurement apparatus was re-calibrated at the end of every day, but an experiment spanned more than one day.)
- the dependent variable. This might help pick out cases of non-uniform variance.
Look at your plots. Think about what you see there. Compare it to your expectations.
Since this data was very simply generated, I'll dispense with some of those plots for brevity. Note that the residuals for this linear fit look vaguely like a quadratic polynomial. This is often the case when there is lack of fit in a polynomial. That lack of fit often looks like the first term we truncated from the Taylor series.
res = y - yhat; plot(x,res,'bo') xlabel X ylabel Residuals grid on title 'Residuals for the linear fit'
If the residuals looked vaguely parabolic in shape, then it might make sense to use a second order (quadratic) polynomial for our fit.
P2 = polyfit(x,y,2) yhat = polyval(P2,x); plot(x,y,'bo',x,yhat,'r-') xlabel X ylabel Y grid on title 'Quadratic polynomial fit'
P2 =
0.5402 1.1140 0.9956
Note that the residuals here look vaguely like a cubic polynomial, although they are much smaller in magnitude than the previous fit. Again, at each step as we increase the order of the model, the residuals will often tend to look much like a polynomial of the next higher order.
res = y - yhat; plot(x,res,'bo') xlabel X ylabel Residuals grid on title 'Residuals for the quadratic fit'
Higher order polynomials - are more terms always better?
Polynomial modeling with polyfit is indeed simple and easy to do. In fact, we might decide to just use higher and higher order polynomials, always chasing the term we truncated in the previous model. After all, if a quadratic model is better than a linear one, then why not go to a cubic? Ten terms must be better than nine. Look at a tenth order model.
P10 = polyfit(x,y,10); disp(P10(1:6)) disp(P10(7:11))
0.0000 0.0000 0.0000 0.0002 0.0014 0.0083
0.0417 0.1667 0.5000 1.0000 1.0000
We can write the Maclaurin series representation for the exponential function as
We can compare P10 to the coefficients of the known series. How well did we do in the fit?
series = 1./factorial(10:-1:0); disp(series(1:6)) disp(series(7:11))
0.0000 0.0000 0.0000 0.0002 0.0014 0.0083
0.0417 0.1667 0.5000 1.0000 1.0000
Note that the higher order coefficients deviate somewhat from the known series, although the lower order terms appear to be quite accurate.
yhat = polyval(P10,x); plot(x,y,'bo',x,yhat,'r-') xlabel X ylabel Y grid on title 'Tenth order polynomial fit'
The residuals oscillate tightly around zero. In fact, they are so small that this last polynomial begins to approach a true interpolating polynomial. Perhaps if we just add a few more terms we may get there. The numerical issues of floating point arithmetic will often preclude true interpolation down to the least significant bit anyway.
res = y - yhat; plot(x,res,'bo') xlabel X ylabel Residuals grid on title 'Residuals for the tenth order fit'
What do you do with interpolation?
I'll start talking about true interpolation in my next blog. But remember that interpolation is different from the approximations provided by polyfit or any other regression modeling tool.
Until that time, please give me your comments on this blog, or ideas for future blog topics on interpolation or modeling in general. Do you have some interesting applications of interpolation? Some interesting problems?
Get
the MATLAB code
Published with MATLAB® 7.6
John,
Do you have any thoughts in the way of rules of thumb for polynomial fitting with relatively sparse data. I will often have a situation where I’ll have ten samples at each of 10 discrete values of the independent variable. Can one issue a general statement such as “Make sure your polynomial has at most N-3 terms”? (Where N would be the number of discrete values.) Likewise, are there any great rules of thumb if one has to consider the possibility of extrapolation?
Thanks,
Dan
PS. At least as it is being displayed, there is no difference between the series coefficients and the polyfit values.
Dan,
Sorry. I looks like a format difference when the blog was published. The two sets of coefficients truly are a little different, mainly in the higher order terms. Here they are as I saw them:
P10 =
Columns 1 through 6
2.8157e-07 2.8227e-06 2.4795e-05 0.00019835 0.0013889 0.0083334
Columns 7 through 11
0.041667 0.16667 0.5 1 1
ans =
Columns 1 through 6
2.7557e-07 2.7557e-06 2.4802e-05 0.00019841 0.0013889 0.0083333
Columns 7 through 11
0.041667 0.16667 0.5 1 1
Good question about the choice of model order. I’m not sure that a rule of N-k terms, or only 50% as many, or ay such arbitrary rule will work, since I am also sure that I could cook up an example where such a chosen rule will fail.
There are some things that you can try however. The extra data that you have is a boon. For example, try this example:
x = -5:5;
x = repmat(-5:5,10,1);
xfit = x(1:5,:);
xval = x(6:10,:);
y = sin(x/4) + randn(size(x))/20;
plot(x,y,’.')
yfit = y(1:5,:);
yval = y(6:10,:);
valerror = zeros(1,10);
for i = 1:10
P = polyfit(xfit,yfit,i);
valres = yval – polyval(P,xval);
valerror(i) = std(valres(:));
end
plot(1:10,valerror)
xlabel ‘Polynomial order’
ylabel ‘Validation error’
It is fairly logical that a first or second order polynomial would be inadequate for this data, since there is a clear inflection points in the curve. What I did was to split the data in half, using some of the points for the fit, the rest purely for validation purposes. Look at the curve of validation error as a function of order. As I added the higher order terms beyond cubic, I start to chase the noise. Those higher order terms are not justified here, not in terms of the signal to noise ratio in this data. The minimum validation error came from a cubic model.
My very last step would be to use the entire data set to fit the chosen order model, a cubic in this case.
P = polyfit(x,y,3);
Finally, I’ll be honest and state that my own choice would not be to fit a polynomial model at all for many such problems. I’d use a least squares spline, constrained if necessary to have the behavior I’d expect, here for example to be monotone.
John
Dan,
There are often domain-specific rules of thumb for the model order. For example, when fitting speech waveforms to an autoregressive model, a good approximation for AR model order is 2 * the number of spectral peaks [poles] one expects to find (the idea is that it takes two AR coefficients to capture information about one spectral peak).
A more general approach is to use an information-theoretic criterion: Akaike Information Criterion (AIC) or Minimum Description Length (MDL). See
http://en.wikipedia.org/wiki/Model_selection
These are heuristics for estimating the marginal “bang for the buck” – assume that each additional parameter comes with a certain cost, and assess when this cost outweighs the decrease in variance. The different criteria (AIC, BIC, MDL, etc.) differ in how they assign the cost, their expectation of noise in the parameter estimation, and so on. Since they are heuristics rather than proven optimal methods, knowing which one to apply can be a bit of an art.
Gautam
Wonderful column that introduced me to polyfit. I figured I’d learn it some day, but haven’t had the need.
Hi.
In the blog post I see a question stated as to wheather higher orders are better, and looking at the examples one might be led to believe that is so.
At this time one might add a word of caution against extremely high polynomial orders are they tend to become very moody close to the edges of the measured data and are completely undefined outside the domain of the measured data.
I have to say I am looking forward to the rest of this series as interpolation is something I allways found to be somewhat of a hassle but have often seen the benefit of if I was more familiar with the tools.
Daniel,
High order polynomials are almost never the interpolation tool of choice. But I wanted to start out this way, increasing the order until you run out of room. I’ll give some examples in my next blog that show off polynomial failures.
While I have plans for a few more blogs, I’d love to hear if anyone has specific areas I should cover.
John
John and Gautam,
Thank you both for your responses to my query. I look forward to reading more on the interpolation aspects of Matlab (especially when it regards rapid computation).
Regards,
Dan
I use the thin-plate spline for interpolating velocity fields for a particle image velocimetry (PIV) algorithm. The velocity vectors are calculated at discrete locations – and so the spline is needed to evaluate the velocity field at any arbitrary set of image coordinates. A spline can be fit to each component of the velocity vector (vx, vy). I like to think of this procedure as “surface” interpolation.
Hi John,
Thank you for the blog post. I have developed a image processing program to track the shape of fluctuating protein filaments. The program outputs a bunch of coordinates representing the ‘skeleton’ in sub-pixel resolution.
At the end, I need to fit a function/model to the coordinate so I can analyze the shape; thus the function/model will be treated as the filament itself.
I naturally went to polynomial fit of the data, generally using 10th order. I do notice the fit becomes ‘moody’ at the edge and introduces error. May I know if there’d be a better model, or a proper implementation of polyfit so I can fit the coordinates (I need them fitted very tightly) while not chasing noise? S/N ratio should be pretty high in my case.
For n data points, polyfit (x,y,n-1) is an interpolant. So polyfit can be used for interpolation.
For the concept that high order interpolation is a bad idea, check out
http://autarkaw.wordpress.com/2008/06/14/higher-order-interpolation-is-a-bad-idea/
http://autarkaw.wordpress.com/2008/06/16/a-simple-matlab-program-to-show-that-high-order-interpolation-is-a-bad-idea/
In regards to finding the optimum order of polynomial to use for regression, one finds where
sr(m)/(n-m-1) becomes either a minimum or does not vary much by increasing m any further.
sr(m)= sum of square of residuals for mth order polynomial
n= number of data points
m= order of polynomial (m<n)
Yes, you can use the overall sum of squares of residuals to choose an order. Given enough data, the validation trick I suggested should be as useful though, and it will often have an easily identified minimum to pick out.
The next step (my next blog) will point out that polyfit(x,y,n-1) is indeed an interpolant, although it does more work than is necessary. In fact, a quick test for a 14th order polynomial shows that it is roughly 15 times faster to use a simple backslash. (I’ve warmed up vander and polyfit in the test below.)
x = linspace(-1,1,15)’;
y = rand(size(x));
tic,polyfit(x,y,14);toc
Elapsed time is 0.009136 seconds.
tic,vander(x)\y;toc
Elapsed time is 0.000605 seconds.
Campion,
I’d suggest avoiding the use of high order polynomials here. While they are indeed simple and easy to use, they can be difficult too as you have found.
Can you use a least squares spline instead? They are often much more well behaved, and they can be controlled nicely in terms of monotonicity, curvature, etc., at least if you have knowledge of their expected behavior.
John
I remember two problems. One concerning interpolation the other one approximation. The interpolation problem was the question how to connect points in the plane by clothoidal splines. I found an article or maybe diploma theses on this at the internet. But this is years ago and the paper work has gone. For this think of the problem to build a track for a road or a railway. The spline should fulfill the condition that the curvature is linear. Is there always or sometimes a solution for cubic splines ? I don´t remember the mathematics and I´m not a mathematician.
Hello John,
I have enjoyed gleaning information from you on this topic for years. I look forward to future blogs, but I would like to ask if you can recommend a good book on this subject.
To let you know what kind of book I have in mind, I did double major in Physics/Mathematics, so an applied book at upper-division level or graduate level would be great. The only books I have found are Springer publishing books (the yellow ones) that were too specialized. I would like a book that briefly covered the basics then went on to advanced applied topics. Thanks.
Matt,
The classic is de Boor, “A Practical Guide to Splines”. It is quite readable, and covers some very useful material. It was the first book I read on the subject, and still the first I’d recommend.
A book that I always thought quite readable, even for an advanced undergrad, was Lancaster and Salkauskas. “Curve and Surface Fitting, an Introduction”.
john
Sorry. I missed the question from Campion. 10th order polynomial interpolation/approximation will indeed be “moody” at the ends. This is why I strongly recommend least squares splines instead. A least squares spline is much easier to control than is a high order polynomial, and a cubic spline tends to be often the best choice, not TOO flexible, but just capable enough to fit almost any curve shape of interest. Splines are also nice things to handle for modeling.
For example, the second derivative of a cubic spline can be no more complex than piecewise linear, so it is quite easy to constrain a spline to be positively or negatively curved over a domain.
Not quite so easy is to constrain a cubic spline to be monotonic, but that too is possible in a linear regression context. Other constraints too are possible. In total, splines can form a splendid way to approximate data where you do not have a mechanistic or physical model for the process, yet you need something that behaves itself.
The next blogs I do will begin to discuss splines, starting with baby steps as piecewise linear interpolants, then next into the realm of cubic Hermite functions and cubic splines. Perhaps in a future blog I could even spend some time on least squares approximation by splines.
John
That’s truly wonderful and helpful.
I am working with some real time series data and plan to use interpolation to take account of missing datas from a radar scaterring of space observation. However, I am afraid that importantly significant frequency componetns will be lost if I use interpolation. Also I would like to know what the effect of this interpolation will be when I take the fft of the interpolated data.
I am not sure, if you could have delt with the interpolation before I send you this comment. Even then, i really love to read it, given the site.
Thx a lot
Abiyu
Hi Abiyu,
I’m afraid that I can’t help you too much here. You can look at the MTF for an interpolation method to learn what it will do to the frequencies of interest. See pages 73+ in Vollmerhausen & Driggers, “Analysis of Sampled Imaging Systems” for example. (You can find a large part of this book on Google Books.)
John
I’m trying to use polyfit to fit a nearly vertical straight line and a very poor fit results.
This is not mentioned in the polyfit documentation as a known issue.
Can you comment on this?
Thanks,
Jennifer
Jennifer,
I’d not expect this to be mentioned as an “issue”. A vertical line has an infinite slope. It has no finite value for the slope. Since polyfit estimates a model of the form y = a1*x + a2, what coefficients would you wish polyfit to return? For example,
x = 2 + randn(5,1)*1e-8; y = randn(5,1); [x,y] ans = 1.99999999813291 0.11393131352081 2.00000000725791 1.06676821135919 1.99999999411683 0.0592814605236053 2.00000002183186 -0.095648405483669 1.99999999863604 -0.832349463650022 plot(x,y,'o') axis equalThis is a vertical line, to within any rational limits. The equation of the line is essentially x=2. That is, y is not a function of x at all. So when you try to fit a line using polyfit,
P = polyfit(x,y,1) P = 10177980.6056492 -20355961.1895639The coefficients are garbage, purely so. But we should expect that! Garbage in, garbage out. When you ask polyfit to fit the model it is designed to fit on this set of data, why would you expect differently?
Having said that, there are several alternatives. One can swap the axes, fitting a model of the form x = b1*y + b2, still using polyfit.
P = polyfit(y,x,1) P = 2.71191329645078e-09 2.0000000038259Again, as expected. The model suggests that x=2, is independent of the value of y.
You might also choose to use an orthogonal regression code. There should be several such tools to be found by a quick search on the file exchange.
John
Hello Jon,
Is there any easy way to interpolate a general function f from R^3 into R^3 given a list of N points that the function maps? I.e. I have a list of N “fixed points” where the behavior of the function f is known: f([x_1,y_1,z_1])=[u1,v1,w1], f([x_2,y_2,z_2])=[u2,v2,w2], …, f([x_N,y_N,z_N])=[uN,vN,wN] and I want to estimate f([x,y,z]) for any point [x,y,z] in R^3.
Thanks,
Geoff
There are several schemes for doing this type of interpolation. Griddata3 (or griddatan) provides one way. Here we treat each output independently as a 3-1 variable mapping of the form f(x,y,z). In that case, griddata forms a tessellation of the scattered domain (x,y,z). Then any point must lie inside one of these simplexes, or along a possibly shared facet thereof. Linear interpolation is then done within the simplex. Think of this as building three independent higher dimensional surfaces to form the mapping.
The second scheme is to consider this as a deformation of a domain in three dimensions. To do this, we might perhaps use a finite element code to solve for the deformation of an elastic body in three dimensions. Any single point in the domain then maps to a corresponding point in the range space.
Both of these schemes have their associated merits. A problem with the elastic body solver is what I like to call unwanted baggage in the metaphorical model. Think of it like this. Suppose I were to compress a piece of rubber along one axis. When I do so, my expectation is that the rubber will bulge out into the other dimensions. This is expected behavior for an elastic body. (A parameter in a physical model is Poisson’s ratio, which helps us to quantify how a given material behaves under deformation.) But suppose I were to do a similar operation on a set of data points in some three dimensional domain? Here a simple scaling along one axis need not be accompanied by a bulge in the other directions. It turns out that if you do use an elastic body solver here, then I’ll argue that it probably makes sense to use a Poisson’s ration of zero.
My point is that use of a metaphorical model such as the deformations of an elastic body to interpolate data may make some sense. A cubic spline is my favorite member of the family of metaphorical models, wherein we use a mathematical model for one physical system to predict the behavior or a completely, unrelated physical system. But at the same time any such metaphorical model carries along baggage. That baggage is the set of sometimes strange predictions that we will see if we are less than careful in our usage of a metaphor.
Climbing down now from my metaphorical soapbox…
John
Hi John, thanks for these blogs on interpolation- they’ve been very helpful so far.
However, I am having trouble working out how to find the polynomial fit to a 3D surface, when instead of having 3 vectors for x,y and z I simply have a single matrix, Z, of size [m x n]. The Z values represent the elevations/surface points, obviuosly I can view this with surf(Z) and see my real surface with no problems. I want to find the function Z = f(x,y) that best fits my surface.
I have seen your polyfitn function but I cannot work out how to prepare my matrix Z for input. Do you have any suggestions?
Many thanks,
Gareth
Hi Gareth,
There are two issues here. One is to build the independent variables for the fit itself. (This is the easy part, but even there we will trip over issues to deal with.) The second issue is what model one will choose.
I’ll use the MATLAB standard here to define x and y. So if we have an mxn matrix Z, then we have m points in y, and n in x. We can simply choose to define x and y as living on the integers 1:n and 1:m respectively. Don’t forget that we have a n entire matrix though! The MESHGRID function will allow us to build arrays of he right size and shape for x and y. For example, given a 3×5 array for Z, we might do this:
m = 3; n = 5; [ x, y ] = meshgrid(1:n,1:m) x = 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 y = 1 1 1 1 1 2 2 2 2 2 3 3 3 3 3Here x and y are now both mxn matrices, of the same size as Z. How might we build the polynomial model itself? For this, x, suppose we simply wanted a linear model in both x and y? We can use either backslash to do the work, or we might use regression codes like polyfitn or regress from the statistics toolbox. Or, if you have the curvefitting toolbox, it now allows you to do fits in two independent variables. I’ll show you how to do it using some of those alternatives.
Before I do any actual examples, I’ll need some actual data. Ok, made up data.
Z = -3 + x + 2.5*y + randn(m,n)/3 Z = 0.53798 1.4681 2.0546 3.2694 3.9688 3.3556 3.7226 5.2381 6.286 6.5197 5.5198 6.5981 8.0412 8.918 9.6904 surf(x,y,Z)My model will be of the general form
First, solve this with backslash. For a model that is linear in x and y, this is easy.
coef = [ones(numel(x),1),x(:),y(:)] \ Z(:) coef = -3.3027 0.94052 2.7469Not too bad. I had dumped some noise into the surface when I built it, so you should not expect perfect estimates of the numbers.
My own POLYFITN (from the file exchange) also does this nicely enough of course, along with some useful statistics on the model and the parameters.
model = polyfitn([x(:),y(:)],Z(:),'constant x y') model = ModelTerms: [3x2 double] Coefficients: [-3.3027 0.94052 2.7469] ParameterVar: [0.057808 0.0025134 0.0075402] ParameterStd: [0.24043 0.050134 0.086834] R2: 0.99121 RMSE: 0.2456 VarNames: {'x' 'y'}Or if you have the statistics toolbox, REGRESS is a very good choice. Here we see 95% confidence intervals returned on the parameter estimates.
[B,BINT] = regress(Z(:),[ones(numel(x),1),x(:),y(:)]) B = -3.3027 0.94052 2.7469 BINT = -3.8266 -2.7789 0.83129 1.0498 2.5577 2.9361Having estimated a linear model so easily, one might choose to expand to higher order models. At some point, you will see numerical conditioning problems. You might be forced to work with orthogonal polynomials, or at the very least, to consider centering and scaling the variables to be lie in the interval [-1,1]. This will greatly reduce the potentiality for singular or nearly singular matrices.
Even so, what happens when we try to estimate a higher order model here? Admittedly, we don’t have much data to estimate that model.
model = polyfitn([x(:),y(:)],Z(:),'constant x y x^2 x*y y^2') model = ModelTerms: [6x2 double] Coefficients: [-3.0534 0.99008 2.5181 -0.041565 0.099913 -0.017752] ParameterVar: [0.49589 0.076522 0.37543 0.0016708 0.0035087 0.021052] ParameterStd: [0.7042 0.27663 0.61272 0.040875 0.059234 0.14509] R2: 0.99386 RMSE: 0.20519 VarNames: {'x' 'y'}Yes, we can estimate a fully quadratic model on this data. In fact though, we know the true underlying relationship to have been linear. So we might hope that to see the higher order terms to be statistically insignificant. Look at the ratio of the parameter coefficients to the parameter standard deviations. The absolute value of that ratio should be compared to a student’s t statistic, here with 15-6=9 degrees of freedom.
model.Coefficients./model.ParameterStd ans = -4.336 3.5791 4.1097 -1.0169 1.6867 -0.12235The critical level for the Student’s t is given here by TINV, from the stats toolbox, or from a set of tables.
tinv(0.975,9) ans = 2.2622So it would appear that the x^2, x*y, and y^2 terms are not significantly different from zero in this model, based on this very limited amount of data. I won’t go any further along these lines. Better is to refer you to any competent book on regression analysis. I like Draper and Smith, as that is what I used many years ago.
Higher order polynomial models are not in general something I recommend though. Linear or quadratic models often are useful. They provide information about the parameters, about the relative importance of the variables in a model. But when people start throwing fifth or tenth order models in several dimensions at their data, I’ll claim this is silly. Many dozens of coefficients will mean nothing to them anyway. Worse, high order polynomial models do terrible things when they are used to extrapolate. They can even do strange and nasty things between your data points. (I call this “intrapolation”.)
In this event, I’ll always recommend a spline type model for your system. For example, my own gridfit from the file exchange is one such member of the general class of spline-like models. Or use the splines toolbox here. These tools do not give you a model that you can simply write down. But why does that matter, since an adequate polynomial model might have had dozens of nonsensical coefficients for the fit?
Finally, it seems that I often get requests from people wanting to know how to find the nonlinear “function” that fits their surface. The nonlinear modeling question is an entirely different one, worth several books worth of writing. Having written virtually a small book here, I’ll stop now.
HTH,
John
Hi John,
Wow! thanks so much for such a detailed reply- Your examples really help- and I can see where I was going wrong, the biggest problem was using meshgrid correctly for my scenario and getting the syntax right for polyfitn and the backslash.
polyfitn will definitely be the most useful to me to begin with I think- as I want to experiment with fitting to a surface that looks like a peak/impulse, so I’m looking at at least order 2 to get a parabola shape function.
Thanks again.
Gareth
Dear Jun,
This is doable, but not anything that I’d truly recommend doing. In fact, it is something that I wrote code for a year or so ago, but never posted on the File Exchange. Let me see if I can explain why I did not.
Your problem is related to numerical differentiation, in that you wish to find the function such that when integrated, gives you the resulting set of points along a curve. Thus, between times defined by the vector t, we have a given set of volumes:
I’ll do it the simple way here, and define a piecewise linear function for our result. As such, I will parameterize it in a linear Hermite form initially, so that I need only work with the value of the function at each of the 13 points in t. Think of it as the connect-the-dots function, where the values at each location in t are unknowns.
Now, can we find those unknown values? Write it as a linear system of equations. If the unknown value, y(t) is known to be linear over each interval, but integrate to the given volumes over that interval, then we know that
(y(1) + y(2))/2 = volume(1)
(y(2) + y(3))/2 = volume(2)
(y(3) + y(4))/2 = volume(3)
…
This just effectively uses the trapezoidal rule to do the integration. Remember, the values in the vector y are unknowns! Write it out using linear algebra. I’ve done it here as a bi-diagonal matrix, A.
We can solve for the unknown vector y, such that this system of equations is satisfied, but you must realize that the solution is not unique because a is rank deficient. The standard solutions in MATLAB to the problem A*y = rhs, are y = A\rhs, or y = pinv(A)*rhs. In fact, it turns out that we could choose ANY of an infinite number of solutions. I’ll plot two of them. (Sorry that I can’t include the actual figure in this response, so I’ll also print them out.)
plot(t,[pinv(A)*rhs,A\rhs]) [pinv(A)*rhs,A\rhs] ans = 692.31 -1.0797e-12 1307.7 2000 4692.3 4000 3307.7 4000 1692.3 1000 2307.7 3000 4692.3 4000 4307.7 5000 5692.3 5000 6307.7 7000 2692.3 2000 3307.7 4000 692.31 0In fact, the infinite family of solutions that results can be written in the form:
where k is some arbitrary scalar constant.
A little trick that I might use here is to try to find the smoothest possible such solution. I can get it using a constrained linear least squares (I could have done it using quadratic programming too.) Here I’ll use my own lse from the file exchange in case you don’t have the optimization toolbox. I also show the corresponding lsqlin call. Here is a link to lse if you want it.
http://www.mathworks.com/matlabcentral/fileexchange/13835
reg = full(spdiags(repmat([1 -2 1],11,1),[0 1 2],11,13)); % y_opt = lsqlin(reg,zeros(11,1),[],[],A,rhs); y_opt = lse(reg,zeros(11,1),A,rhs,'pinv') y_opt = 681.82 1318.2 4681.8 3318.2 1681.8 2318.2 4681.8 4318.2 5681.8 6318.2 2681.8 3318.2 681.82Above, reg is a matrix that essentially forms the second difference of y. The regression solver then finds the vector that minimizes the overall curvature of y, while at the same time satisfying the integral constraint for the “interpolation” problem.
We can plot the solution of course, although, really it is not that much smoother than the alternatives I’ve offered.
If you like what you see, this is easy enough to encode into a pp form, so you could use the MATLAB splines tools to evaluate the approximated function. Use mkpp, as I do here:
This may be acceptable to you. Personally, I’m more picky, so I won’t post this on the file exchange, or the cubic spline version that I once wrote. Perhaps it does not seem that poor, and really, this may be as good as you can ask for. So I’ll give an example from a known smooth function, and see if I can reconstruct the function.
x = randn(1000000,1); min(x) ans = -4.8521 max(x) ans = 4.6365Now, build a histogram. It looks pretty smooth to me from the plot. After all, I had 1e6 data points to build it from.
Can I reconstruct the original Gaussian shape from only the counts? My code would look like this:
Again, since the plot won’t come through this response, here is the result. See what happens in the tails.
y_opt = 596.79 -592.79 594.79 -572.79 612.79 -514.79 818.79 -170.79 1730.8 1533.2 4938.8 6951.2 13993 20167 33075 44527 61465 75207 88643 97129 1.0023e+05 96623 88613 75299 60327 46391 31485 22363 12363 8511.2 3382.8 3051.2 204.79 1425.2 -681.21 977.21 -849.21 883.21 -877.21 883.21 -883.21While it does look roughly like a Gaussian, it has an annoying habit of oscillating in the tails. Those oscillations reflect the problem of implicit numerical differentiation from raw data.
John
Dirga-
You might consider looking at the function interp2.
–loren
Dear John,
I need to find the intersection of an 3D object with its orthogonal plane(to the centerline of a non-straight cylindrical object). For this reason I need to do following step:
1-finding the intersection of the orthogonal plane and 3D vessel (I have the normal vector and a point of plan)
2-Mapping the interpolated value of each voxel(pixel in 3d) to related area on plane(cross-sectional area)
each voxel only have value on top and button surface
3-Transform voxels on the mapped plane such that plane normal coincides with the z axis
for mapping and interpolation I need to do some thing like interpolation on the face of a triangle, like the GPU does on the face of a triangle for the texture coordinates. I am not interested in the perspective correction aspect as this is only needed for the CPU and not for visualization. Do you have any suggestion for me?
I use matlab for implementation.
best regards,
mehdi
Hi john
I need to approximate Ackleys function of order 2 and 10 with a Polynomial Regression method, so I have a M×2 or M×10 Matrix as an input variable (M is the number of sample points). But I can’t use “polyfit” structure since the size of x and y should be the same.
I would like to know how can I approximate a multi variable function with Polynomial Regression in Matlab?
Best regards,
Golnaz Mir
Hi Golnaz,
There are many ways to do this in matlab. If you have the statistics toolbox, then the regress function is available. Or, if that is not an option, then the optimization toolbox (use lsqlin) or the curve fitting toolboxes can help. In matlab itself, there is lscov, or just the backslash operator.
Any of those tools CAN solve the problem, but they will all require you to build the matrix that defines your model. If you don’t know much about regression modeling, this is not a bad thing to learn, at least if you will solve this class of problem again in the future.
This design matrix is a matrix with one column for each term in your model. For example, a linear model in x and y (where x and y are column vectors) would have you build the matrix A:
And then solve the problem as
The column of ones in A corresponds to the constant term in the model.
You can build the design matrix A for more complex polynomial models too of course. However, rather than force people to know how to do all of this, as well as go through the effort of building these matrices, I posted a function several years ago on the file exchange to help them with this task. It is called polyfitn, which suggests that it is just an extension of polyfit to problems with several independent variables. You can download it from this link:
http://www.mathworks.com/matlabcentral/fileexchange/10065
Polyfitn will allow you to estimate any simple polynomial model as you need to build. Of course, polyfitn will not help you to decide what order polynomial mode to use. That is a problem that usually requires knowledge of your system, in terms of how much noise is present in the data. It also requires some appreciation of your goals in the process, and what you will do with the model. How accurately do you need to fit the data? Will you need to use it for extrapolation?
A few extra thoughts after my last response as a cautionary tale…
Too often I see people decide to use a polynomial mode without understanding what they are doing, who have no idea what order model they should use. So they fit a low order model, and it does not give them as good a fit as they want/hope to see. So they try a higher order model. It fits a bit better, so then they decide to fit a really high order model. Of course at some point they end up with such a high order model that it fits the data points perfectly, yet what they don’t realize is this high order model is usually completely meaningless. They have managed to over-fit their data so horrendously that while the polynomial passes through the data points, it also has so much flexibility that it does all sorts of crazy things between the data points.
Polynomial modeling is both easy to do and difficult, since it is so easy to over-fit your data. Which terms belong in the model? How high an order do you use?
For these reasons I rarely ever recommend fitting a high order polynomial model. Stick with the low order polynomials, linear, quadratics, cubics where necessary. In several dimensions even a full cubic model has many terms in it. If you need more flexibility that that, then it is often a good idea to look at spline models, the shape of which can be controlled using the right tools. Alternatively, look into nonlinear model that have the desired shape built right into the model.
Hello Dear John
Several months ago I have asked you a question about approximation of Ackley’s function; actually I would like to use the polynomial approximated form of Ackley as a meta-model in an optimization problem.
As you recommended I used the “polyfitn” tool for approximation, the second order polynomial was my choice in both 2 and 10 dimensional problem, but even for 2 dimensional Ackley I could not get good results.
The code I used for 2 dim approximations was this:
inrnd =200; % Number of initial random data generation for the first meta-modelling
p(:,:) = xmin*ones(N,inrnd)+(xmax-xmin)*rand(N,inrnd);
t(:,:) = 20-20*exp(-0.2*sqrt(1/N*sum(p(:,1:inrnd).^2)))-exp(1/N*sum(cos(2*pi*p(:,1:inrnd))))+exp(1);
inrnallb = 10;
p(:,inrnd+1:inrnd+inrnallb) = xmin*ones(N,inrnallb)+(xmax-xmin)*randint(N,inrnallb,[0,1]);
t(:,inrnd+1:inrnd+inrnallb) = 20-20*exp(-0.2*sqrt(1/N*sum(p(:,inrnd+1:inrnd+inrnallb).^2)))-exp(1/N*sum(cos(2*pi*p(:,inrnd+1:inrnd+inrnallb))))+exp(1);
inrnoneb = 20; % Each variable has its boundary
for i = (inrnallb+inrnd+1) : (inrnallb+inrnd+N*inrnoneb)
st = ceil((i-(inrnallb+inrnd))/inrnoneb);
if st == 1
p(1,i) = xmin*ones(1,1)+(xmax-xmin)*randint(1,1,[0,1]);
else
p(1,i) = xmin*ones(1,1)+(xmax-xmin)*rand(1,1);
end
if st == 2
p(2,i) = xmin*ones(1,1)+(xmax-xmin)*randint(1,1,[0,1]);
else
p(2,i) = xmin*ones(1,1)+(xmax-xmin)*rand(1,1);
end
t(:,i) = 20-20*exp(-0.2*sqrt(1/N*sum(p(:,i).^2)))-exp(1/N*sum(cos(2*pi*p(:,i))))+exp(1);
end
w=size(t);
[ptn,ps] = mapstd(p);
[tn,ts] = mapstd(t);
[R,Q(1)] = size(ptn);
numexactfiteval = Q(1);
iitr = [1:5:Q(1) 2:5:Q(1) 3:5:Q(1) 4:5:Q(1) 5:5:Q(1)];
ptr = ptn(:,iitr); ttr = tn(:,iitr);
%……………..function approximation…………….
ptr=ptr’;
modelterms = [0 0 ;1 0;0 1;1 1;2 0;1 2];
polymodel = polyfitn(ptr,ttr,modelterms);
zn=polyvaln(polymodel,ptr);
z=mapstd(‘reverse’,zn,ts);
z=z’;
for i=1:1
figure(i)
[m(i),b(i),r(i)] = postreg(z(i,:),t(i,:)); % regression khati beine natayeje shabake asabi va maghadire asli(target)
pause(1);
end
I would like to know if I have made any mistakes in using “polyfitn” tools?
Thanks,
Best Regards,
Golnaz mir.
This is a hard one for me to answer, since I don’t know the values of some of your parameters. What are N, xmin, and xmax? A wild guess is that N = 2, since you are talking about 2 dimensions in this example. But what are xmin, xmax?
I also don’t have the neural net toolbox, so I can’t test your code. A check online tells me that mapstd shifts and scale the columns of your array to have mean 0, standard deviation 1. And since I don’t have the communications toolbox, I don’t have randint either, but that too I was able to replace.
Finally, I tried N = 2, xmin = -1, xmax = 1. These seemed like good, generic numbers, although I have no idea if they are reasonable or correct. When I did that, and wrote a simple version of mapstd of my own, the very first thing I tried was to plot what you are trying to fit.
NEVER, EVER just blindly fit a function without plotting it if you can do so, without knowing what you are trying to fit, and without thinking about whether the model you pose is meaningful for that data.
When I do so, I see something that is wildly NOT a polynomial. And worse, it is surely not even close to being something that a low order (quadratic) polynomial could ever fit. It has what appears to be possibly a singularity at (0,0). No polynomial model exists that has a singularity in it.
So, no, you cannot use polyfitn for this problem, with ANY order of polynomial model.
If I have tried the wrong parameters and so am wrong about what I’ve seen, please ask again.
John
Hello Dear John
Well first thanks a lot for answering me,
Actually I forgot to specify the parameters for you, Xmin=-5 and Xmax=5 and N=2, I used the polyfitn toolbox from matlab exchange files in Mathwork.com.
Best Regards,
Golnaz Mir
Hi Golnaz,
Given those values of xmin and xmax, try changing inrnd to 10000. Then add this line near the end, but before you pother to try fitting any models.
Rotate this plot around, and look at it carefully. First of all, does this look like something that a quadratic model would generate? Recall that a quadratic model as you will form it will take a parabolic shape, with two independent variables. This data is not remotely parabolic. It has a singularity at (0,0), something that no polynomial model will ever yield.
Worse that that, the surface is actually rather bumpy, not at all smooth like a low order polynomial will yield. In fact, it is rather complex.
I’m sorry to say that no simple model will allow you to easily approximate your surface, although you might find some other modeling technique to be of utility. As an offhand idea, perhaps wavelets might be of value.
John
hello ..a have a question..how do i write Ferguson and Hermite-Coons interpolator in matlab,for conditions: r0:=[1,2,3]; r1:=[2,5,1]; r`0:=[0,-1,0]; r`1:=[1,0,1]? thanks
Hello John,
It seems the link to polyfitn is invalid. Could you please post a valid link?
Thanks,
Marc
Marc-
John elected to remove some of this files from the file exchange.
–Loren
Hello,
I am not sure if this blog is still active but I hope it is.
I am looking at interpolating a rather simple set of data. I have a matrix of weekly data (dates and values) for one year in the past. For the last week (say 24th Jul to 30th July), I need to interpolate values for each day, based on available past data.
How can I do this in matlab?
Regards.