Polyfitn is an extension of polyfit, allowing the user to create models with more than one independent variable. It also allows the user to specify a general model, for example, a quadratic model, with constant and quadratic terms, but no linear term.
For example, to fit a polynomial model to points selected from a cosine curve, we will only need the even ordered terms.
x = 2:.1:2;
y = cos(x);
p = polyfitn(x,y,'constant x^2 x^4 x^6');
p.Coefficients
ans =
[0.99996 0.49968 0.041242 0.0012079]
The coefficients won't be exact of course, as I used only a finite number of terms for what is essentially a truncated Taylor series, and I had only a finite amount of points to build the model from. The exact first 4 coefficients for a cosine series would have been:
>> [1 1/2 1/24 1/720]
ans =
1 0.5 0.041667 0.0013889
so we got the expected result.
Of course, polyfitn works in higher dimensions, as it was this problem it was really designed to solve.
x = rand(100,1);
y = rand(100,1);
z = exp(x+y) + randn(100,1)/100;
p = polyfitn([x,y],z,3);
The result can be converted into a symbolic form to view the model more simply. Here I'll use my sympoly toolbox, but there is also a polyn2sym function provided.
polyn2sympoly(p)
ans =
0.43896*X1^3 + 1.4919*X1^2*X2 + 0.041084*X1^2 + 1.4615*X1*X2^2  0.095977*X1*X2 + 1.2799*X1 + 0.56912*X2^3  0.15306*X2^2 + 1.361*X2 + 0.94819
And of course, parameter error estimates are generated for those who wish to determine the significance of the terms generated.
I've also provided tools for evaluating these models, and to differentiate a model.
A caveat  beware the use of high order polynomials to fit your data. Just because a low order model works, a high order model is not necessarily better. High order polynomials often suffer from severe ringing between the data points. Always plot your data. Think about the model you will be building. Then plot the resulting model. Use your eyes to validate the result, more so than simply looking at an rsquared coefficient (although I do return that parameter too.)
If you do find that a high order polynomial mode is necessary because your curve is simply too complicated, consider using a regression or smoothing spline model instead.
1.3  Documentation added for p and DoF, also an example was added in the help. 

1.2  Remove NaNs in data for fit 

1.1  Added a pvalue on the coefficients. 
Inspired: polyVal2D and polyFit2D, gapolyfitn
Issac Gan (view profile)
Qin Xu (view profile)
Tobias Eberhorn (view profile)
Hi John,
Thanks for your polyfitn toolbox, it's a great help!
I'm working on a problem with 45 input parameters (imported from an ExcelFile, one column vectors) and one parameter on which the other parameters depend on.
I validated my function with the CF toolbox with two parameters.
When I tried it with a third input parameter I got a warning:
Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.794064e17.
In addition I checked the R2 (=0.9495) and the result is in my opinion too well adjusted.
Could you think of a possible way to fix this or am I missing something?
I thought the toolbox was especially made for working with multidimensional problems.
If you could think of anything which could help me, I would really appreciate it!
Thanks a lot,
Cheers Tobias
Here's my coding.
x = Durchmesser3;
y = Plattendicke3;
w = EModulNiet3;
z = Kraft3;
p = polyfit([x,y,w],z,2);
if exist('sympoly') == 2
polyn2sympoly(p)
end
if exist('sym') == 2
polyn2sym(p)
end
John D'Errico (view profile)
@Ramon:
a + b*sqrt(x):
polyfitn(x,y,[0;1/2])
a + b/x:
polyfitn(x,y,[0;1])
Jaime (view profile)
John!!! THANKS. No more words, I was going to program it. THANKS! ;)
Ramon Dalmau (view profile)
It is polyfitn able to fit functions such that :
a + b*sqrt(x) or a + b/x
?
Michele Girardi (view profile)
Moritz Kick (view profile)
Thank you so much for this Toolbox, great help!
Andrejs Fedjajevs (view profile)
Is it able to substitute stepwiseglm functionality from statistics toolbox?
Sebastian Bahamonde (view profile)
Jay (view profile)
Is there any specific conference paper that explains whole process of polyfitn?
It would really help to dig into the algorithm, rather than making a try to understand the coding.
Thanks
PerJohan (view profile)
Rose Ann Haft (view profile)
Does anyone have an example of how to use this in 3d?
Thanks!
Till Huelnhagen (view profile)
I was looking for a way to fit 3D fields using polynomials. This function solves the task perfectly and quickly.
Thank you for sharing this!
Daniel Cipriano (view profile)
Hi, John! I'm a starter in Matlab. I need to find a polynomial fit from a surface, having two vectors of values x and y.
x is 1x191 and y is 1x51.
Then I have the plane z, which dimensions are 51x191.
I can plot it with:
surface(x,y,z)
But I cannot quiet understand how can I use your file in order to find my equation.
Thanks for your answer!
AndBu (view profile)
Robert Wylie (view profile)
I have overcome the rookiness and managed to use it successfully now. complete wizardry! thank you.
Robert Wylie (view profile)
Hi John,
Can i start by saying im a severe matlab rookie. I have three sets of recorded data, i used scatteredInterpolant and meshgrid to create a surface plot. I now want find a formula for this surface.but when inputting the three data sets i get "Error using betainc
Z must be real and
nonnegative.
Error in polyfitn (line
266)
polymodel.p =
betainc(polymodel.DoF./(t.^2
+
polymodel.DoF),polymodel.DoF/2,1/2);"
the Z variables are all positive numbers form 25 to 26. I'm not sure if this is a problem with using three specific data sets as oppose to getting data out of a function.
Many Thanks.
Alessio Scanziani (view profile)
John D'Errico (view profile)
Mimo  I'm not sure what you are asking.
Are you asking if you have only a single variable x, but multiple y vectors? polyfitn is not set up to solve a problem of that form, but it would be just a loop inside polyfitn anyway. I suppose I could have written polyfitn to loop over the columns of y. An idea for the next time I post a release. :)
The question is if you are using polyfitn to simply compute regression coefficients, or if you want the whole nine yards, i.e., standard errors, R^2, etc.
For example, IF your problem is as simpe as you seem to describe...
Assume that x is an nx1 column vector of data, the independent variable vector, fixed over all problems. y is an nxp array, so there are p separate problems that you need to solve. Then the regression coefficients are easily found as:
ab = [ones(n,1),x(:)]\y;
the array ab will be a 2xp array, where each column of ab is one set of regression coefficients. So the first row of ab will be the constant terms for each problem.
a = ab(1,:);
b = ab(2,:);
Again, this does not offer any statistics, R^2, etc. But it is a simple solution to the problem that I think you have posed.
Mimo (view profile)
Hi John,
This function is amazing. I just came across with this question that if it is possible to also find a multilinear regression in form of “y(i)=a(i)+b(i).x”, having only a single variable “x”, by finding the vectors a(i) and b(i)?
John D'Errico (view profile)
Assuming that phi and theta are variables, then just do this:
mdl = polyfitn([phi(:),theta(:)],z(:),2);
You may need to use meshgrid or ndgrid (as appropriate) to generate phi and theta.
Héctor Esteban Cabezos (view profile)
Hi guys! I've been working in a 3D plot, having phi (points from pi/3 to pi/3) and teta the same and I have a value (z) for each point, so finally I get a matrix for each possible combination of phi and teta. I want to approximate polynomically the points of the matrix I have in function of phi and teta. I've looking on the polyfitn but the points you have to put (z) could not be a matrix. Polyfit doesn't work either because it is only for 2D (phi,z) or (teta,z)
George Varverakis (view profile)
Nice work!
John D'Errico (view profile)
Nothing beyond the basic here. (I admit that basic is in the eye of the beholder.) Just classical multivariate linear regression. See any text like Draper and Smith on regression analysis. So any tool from the stats toolbox will do the same, as long as you prepare the design matrix properly for your model.
Sebastian Bahamonde (view profile)
Dear John
Thanks for your nice tool. I am actively using it, and it works very well. I want to learn the method you are using. I wonder if you know where I can find the associated theory.
Cheers,
Sebastian
netorha (view profile)
Hello John,
thanks a lot for your great script. It works very very well!
But, I've a question concerning a problem that is a little bit more difficult ... Actually I perform a datafit by using matlab fit() with the following code (where arr_node_coords_y, _x and arr_temperaturefield are 33x5Matrices):
[xData, yData, zData] = prepareSurfaceData( arr_node_coords_y, arr_node_coords_x, arr_temperaturefield );
% Set up fittype and options.
ft = fittype( 'poly24' );
opts = fitoptions( ft );
% Fit model to data.
[fitresult, gof] = fit( [xData, yData], zData, ft, opts );
% plot
subplot( 2, 2, 2 );
h = plot( fitresult, [xData, yData], zData );
axis([0 2.5 0 8 20 60]);
xlabel( 'Width [m]' );
ylabel( 'Height [m]' );
zlabel( 'Temperature [°C]' );
title('Temperature from Fit');
grid on;
daspect([8 8 60]);
colorbar;
... how is it possible, to perform the same thing with polyfitn()? ...  I tried it in a few ways but unfortunately I did not find a working solution.
Maybe someone can help me ... ?
rihab (view profile)
Balasundaram Mohan (view profile)
hello john
I'm having data of 100x1,100x1 as a independent variable and dependent variable as a 100x100 and also it's complex number.could you please help me out with polynomial fitting for this?
while i have been trying to do polyfitn for diagonal elements of a dependent variable it's showing 'Error using betainc
Inputs must be real, full, and double or
single'.
what it means and how can i overcome this problem?
John D'Errico (view profile)
@Sagar  but that is not a polynomial model. So you cannot use a tool for polynomial models to fit ANY general nonlinear model. Use a nonlinear regression tool. The curve fitting toolbox is the simplest choice, but there are many others, in the stats or optimization toolbox for starters. Or use an optimizer like fminsearch or fminbod, but there you will need to minimize a sum of squares of the residuals.
Sagar (view profile)
Hi, excellent function. I know we can do: poly_fit = polyfitn(x, y, 3); for testing all combinations of input variables (say x1 and x2). But what should I do if I want to construct customized input exponents, I want to fit something like this:
poly_fit = 2.3x1^3 + 1.5e^x2
How do I do this?
Akash (view profile)
John D'Errico (view profile)
Fixed the NaN issue. Sorry about that.
Chad Greene (view profile)
One of many excellent submissions by John D'Errico.
If the dependent variable contains any NaNs, a confusing betainc error results. I recommend adding this error check early in the code:
assert(sum(isnan(depvar))==0,'Dependent variable cannot contain any NaNs.')
Lukas Theisen (view profile)
JinSun (view profile)
Incredible. Thank you for your contribution!
akshay (view profile)
works great and should be included in standard package
Christoph (view profile)
Muhammad Irfan Zafar (view profile)
Hi John,
I have used 'polyfitn' to find polynomials for varying parameters of nonlinear system based on three input variables. Inserting these polynomials in three equations, i have then tried to 'fsolve' them but i couldn't converge to one solution by inserting polynomials of different orders.
I have been going through your file exchange list to look for some other related code. I was wondering if i should try to use some spline method instead of 'polyfitn' but i am not sure how i can insert it in the equation and solve them. Also is their any other method to solve three nonlinear equations simultaneously.
Thanks
Johan (view profile)
John D'Errico (view profile)
I'm sorry, but I can only assume that Gautam, like Mitch, is using a quite old release of MATLAB.
If you really want to resolve the problem, you would provide sufficient information for me to know if this is the case. Since a test case as you describe works perfectly, as I show here, my suggestion is to use a reasonably recent release, and since Mitch said that he was using a 7 year old release and got that same error, I assume this is the case.
If you are indeed using a recent release, and are still having a problem, then I can only assume that you are either using it improperly, or that you have something else installed that is causing the problem. For example, when I perform this test, it runs perfectly, but I don't have a 7 year old MATLAB release, nor can I even run something that old.
indvar = rand(12,2);
depvar = rand(12,1);
m = [0 0; 1 0; 0 1; 2 0; 1 1; 0 2; 3 0; 2 1; 1 2; 0 3];
p = polyfitn(indvar,depvar,m)
p =
ModelTerms: [10x2 double]
Coefficients: [0.0610 4.6966 5.5355 27.4400 24.9658 0.4228 32.4825 42.3187 18.4250 5.9917]
ParameterVar: [12.5177 249.5518 129.5089 850.2142 511.0641 282.0906 504.6530 953.0042 477.2848 209.7347]
ParameterStd: [3.5380 15.7972 11.3802 29.1584 22.6067 16.7956 22.4645 30.8708 21.8468 14.4822]
DoF: 2
p: [0.9878 0.7943 0.6748 0.4460 0.3845 0.9822 0.2851 0.3040 0.4878 0.7192]
R2: 0.6497
AdjustedR2: 0.9265
RMSE: 0.1327
VarNames: {'' ''}
Without useful information from you, I cannot guess what is the problem. The information I refer to in this context would involve telling me what MATLAB version you have. It would involve telling me the complete error message, not just a fragment of that error message.
By the way, trying fit a data set with only 12 data points with a model that has 10 terms will yield a virtually meaningless result. I would call that a case of MASSIVELY overfitting your data.
Gautam (view profile)
Mr. D'Errico, Thanks very much for your submission. It looks quite helpful, except that when I give a 12 x 2 matrix as the independent variable matrix and a 12 x 1 column matrix for the dependent variable, it displays the following result:
Undefined function or variable "varlist".
Error in polyfitn (line 233)
polymodel.VarNames = varlist;
The model I used is this:
m = [0 0; 1 0; 0 1; 2 0; 1 1; 0 2; 3 0; 2 1; 1 2; 0 3];
Any advice on how to get your function to work for me would be very helpful and greatly appreciated. Thank you very much.
Gautam.
John D'Errico (view profile)
MItch,
It is not relevant that you think it was properly defined in your test. I simply don't know what it was defined as, because you gave me everything else but that variable. So unless you can show me a constructible example that fails, then I have no idea what you are doing. That I used a different variable name is completely irrelevant. I could as easily have named the variables I passed in as Fred and Barney and it would have run.
I also don't know if there is an issue with that old of a MATLAB release. The site tells me I wrote the tool using R2007b, so it SHOULD work for you, but it is also conceivable that one of the mods made to it over the years is now making it fail in an old release. This is something I cannot test since I cannot use that old of a MATLAB release.
So my first suggestion is to run the test I gave. Verify that it does run properly. If it does not run, then you need to see what is different about the input arguments in the test you have. OR you can send me the data you have if you don't wish to post it online.
If the issue is that your old release is indeed the problem, then you will need to learn to solve the problem using backslash (not that difficult a task) or perhaps get a copy of the stats toolbox.
Mitch Martelli (view profile)
Hi John,
Thanks for your quickly replay.
I'm sorry for the missing data in my question, but my variable called err_func is your outpar and it was defined.
I check all the variables dimensions and they are correct.
Then I try to copy and past your example and appears the same error :
??? Undefined function or variable "varlist".
Error in ==> polyfitn at 232
polymodel.VarNames = varlist;
Could be a problem reletad to my Matlab version (R2007b)?
Thanks in advance
Mitch
John D'Errico (view profile)
MItch,
I don't know what you are doing. My initial guess is you tried to "run" the function, from the editor, or something like that. You don't run functions like that, although many people try to do so.
I do see that you have not defined the variable err_func in your example though, so I cannot really even guess what is your error.
inpar = rand(20,2);
outpar = rand(20,1);
model_exp= [1 0;2 0;0 1;1 1;0 2];
p = polyfitn(inpar,outpar,model_exp);
This test case runs nicely.
Mitch Martelli (view profile)
Hi John,
Thanks for your great submission.
I try to use it on my own data. But I’m not able to set correctly modelterms argument:
i.e.
in_par=[zeta theta];
model_exp= [1 0;2 0;0 1;1 1;0 2];
p = polyfitn(in_par,err_func,model_exp)
after run appears the following error:
??? Undefined function or variable "varlist".
Error in ==> polyfitn at 232
polymodel.VarNames = varlist;
Error in ==> launch_regr at 14
p = polyfitn(in_par,err_func,model_exp)
Please can you help me?
francesco (view profile)
John D'Errico (view profile)
I suppose it still begs the question, does polyfitn need to have the ability to fit multiple right hand sides, i.e., the same model, but for multiple dependent variables?
Personally, I've found this needed only somewhat rarely, and when we did use such capability in my previous incarnation supporting a large group of color scientists, it was really only useful for building monitor models. Since then, monitors are more complex, thus less linear. I would conjecture that the newer monitors out there are better handled using more sophisticated models. This is certainly true for printers, etc, which are best handled using 3d lookup tables, fit by a tool like my own gridfit.
But the fact is, I'm often surprised to find others use these tools in ways I did not expect.
So if I see at least a few people requesting this capability, I would consider adding it. Until then, my preference is for simplicity of code in this respect. (Hey, I am supposed to be retired!)
Mark Mikofski (view profile)
I completely agree with John; I threw out that suggestion flippantly, without thinking it all the way through. Sorry if I led you astray. Thanks John for your sound advice and superb mathematical tools  they are invaluable!
John D'Errico (view profile)
I'm sorry, but I completely, whole heartedly disagree with Mark's advice to Anderson.
Systems with multiple output variables require no more than a simple loop around polyfitn. The cost of that loop is simply not significant enough to increase the complexity of the code and the complexity of the interface.
To suggest instead that one use an iterative scheme like GaussNewton, lsqnonlin, or lsqcurvefit, etc., these are all absurd approaches. You will still need to loop, OR you will need to generate a large problem with essentially a block diagonal Jacobian matrix. All of this to avoid a simple loop? Or you will write a loop, just to avoid a loop? Note that in Mark's suggestions, you will need to write a bit of code just to evaluate the model. And if you use a tool like lsqnonlin or lsqcurvefit, then you either need to provide the Jacobian, or you will need to let the tool call your objective function multiple times just to compute finite differences to generate that Jacobian.
Don't throw an iterative tool at a problem to avoid a simple loop, on a problem that is truly solvable by basic linear algebra. See that even while the iterative tools will converge quickly, those tools will actually try to take spare iterations, until it effectively realizes that it has indeed converged. This will require a spare Jacobian evaluation.
Mark Mikofski (view profile)
@Anderson, you could use newtonraphson (http://www.mathworks.com/matlabcentral/fileexchange/43097newtonraphsonsolver) or if you have the optimization toolbox either lsqcurvfit or lsqnonlin (http://www.mathworks.com/help/optim/nonlinearleastsquarescurvefitting.html)
Anderson (view profile)
This tool looks very useful. But do I understand correctly that it will not fit systems with more than one output variable? If so, does anyone know of a tool that does polynomial fitting for systems with more than one output variable?
John D'Errico (view profile)
Eran  It might be interesting to some, but I think it will encourage people to fit models that are too high of an order, something that I am constantly forced to caution people not to do.
Since your idea would require adding more options to the code, it would make the code more complex. As well, using the tool would be less simple, since if many models are generated and returned, then tools like polyvaln and polyn2sym would not work directly, or would also need to be modified.
As far as a significant gain in efficiency is concerned, remember that whenever there are more options in a code, it becomes less efficient for the basic user. Complex code is also more likely to have unintended bugs.
So I would not put in this enhancement myself. Could you do so yourself? Well, yes, you could in theory. To be most efficient, one would need to wrap a loop around the modeling process, probably using the qrupdate function to add or remove terms from the model. I'm not sure that is worth the effort compared to a simple loop around the entire code, just calling polyfitn repeatedly unless you were using the tool extensively.
Eran (view profile)
Hi again,
Thank you so much for your last answer!
Of course it was a dumb question :)
I have another question, that might be a bit more tricky 
lets say that I have a set of n points.
I want to check what is the polynomial fit for every k<n (of course k is big enough to hold enough data for the fit). meaning I want to have (almost) n polynomials fitings.
of course I can do a "for" loop and get my answer... But that is very inefficient.
Is there any way to modify your code slighlty so that this would work?
I think that the main problem is in the part where you define the "design matrix".
Thanks again for the quick help!
John D'Errico (view profile)
As a continuation, let me explain it like this, because this is a common problem in regression.
When you create a dataset where y = x/2, your data lives in a region that does not fill a region of the (x,y) plane. Instead, it lives entirely on the line y = x/2. In effect, x and y are completely dependent on each other, given x, we know the value of y.
The problem is, we have no information about what happens for z when your points fall off that line, so no estimation of the coefficients of a full model can be satisfactory.
To help convince you of what has happened, lets start with the original function z that you created.
z=4+x+x.^2+x.*y;
Substitute the known relation for y into z, to get it as a function only of x. Then we would have:
z = 4 + x + 1.5*x^2
See that we can model z as a function only of x.
out=polyfitn(x,z,2)
out =
ModelTerms: [3x1 double]
Coefficients: [1.5 1 4]
ParameterVar: [1.5362e32 1.6407e28 7.6808e26]
ParameterStd: [1.2394e16 1.2809e14 2.7714e13]
R2: 1
AdjustedR2: 1
RMSE: 9.327e13
VarNames: {'X1'}
No warning was made here, and the fit was essentially perfect, yielding the values for the coefficients we would have expected. I could as easily have done the fit for z as a function of y.
John D'Errico (view profile)
HI Eran,
But think about what you are doing.
At first glance I made the rapid assessment that x was moderately large, and thus this was just a scaling issue. Then I looked at your problem. This is the line that causes the problem!
y = x/2
Your model is quadratic in x and y, in the sense that it has terms x.^2, x.*y, and y.^2. Of course, you hope that polyfitn will tell you that the coefficient of x.^2 was zero, as it implicitly was when you created the data.
But you created y as a simple multiple of x. You should see that x.*y == 2*y.^2 == x.^2/2.
There is NO linear algebra scheme in the world that can intelligently recover the coefficients of these terms, and know which set of coefficients you started with. The solution simply is not unique. This is why you get a warning about singularity. You will get an answer, but not a unique one.
The story here is garbage in, garbage out. Your data does not adequately support estimating the model you wish to fit.
Eran (view profile)
Hi Jhon,
Tahnks you very much for this soneeded function. It's a shame that matlab doesn't support this...
I tried using this function, and encountered a small problem.
I used:
x=0:100;
y=x/2;
z=4+x+x.^2+x.y;
in=[x' y'];
out=polyfitn(in,z,2);
and what I get is a warnning about the singularety of the matrix "R" produced by the function "qr". That of course when I use "\" that practicly uses the inverse of "R".
This happans twice in the lines after using the function "qr".
Do you have any advise on how to solve this?
I'm asking this question of course, becuase I get that the coeficients produced by "polyfitn" are very different from those you would expect...
Thanks you very much!
John D'Errico (view profile)
Hi Alex,
I think you wanted to ezsurf, not ezplot. For example, this worked for me:
xy = rand(100,2);
z = 3 + 2*sum(xy,2) + sum(xy.^2,2) + 4*prod(xy,2) + randn(100,1)/10;
model = polyfitn(xy,z,2);
fun = polyn2sym(model);
ezsurf(fun)
ECOSat Team (view profile)
Hi John,
These functions are wonderful!
I'm using it to model the intersect of two perpendicular cones, which provides me with a parabola in 3D space. I was wondering, is there an inbuilt way to plot fitted curves in their proper dimensions? I tried using your polyn2sym tool and then using ezplot, but it fails to properly show the curve in 3D.
Thanks,
Alex
Ecosat Team
John D'Errico (view profile)
Hi Raymond,
Those parameters are derived to apply to models with noise in them. For example, suppose I was given 10 measurements of a given process at a single point. The least squares best estimate of the constant model, thus y=c, is given by the mean of the samples. At the same time, we could compute an estimate of our uncertainty in that estimate of the mean, using it to compute approximate confidence intervals around the parameters. We can do similar things with more complex regression models, deriving estimates of both the parameters in the model and our uncertainty around those parameters based on the perceived noise in the data. The noise is perceived to be essentially what is left over after the model is estimated. Again, this will give us estimates of standard deviations that we can put on the parameters themselves. (Multiply the provided standard deviation by a Student's tstatistic for an approximate confidence interval at your required level of confidence. Using approximately +/ 2 or 3 sigma is a common thing to do here to get the confidence interval.)
In fact, most of the time when you have more data points than parameters to fit a model, we can also form estimates of our uncertainty in those parameters. The more data you have, the better of course, as it allows us better estimates of the parameters. A nice thing is we don't even need pure replicates to form that estimate of our uncertainty in the parameters. The uncertainty is based on what the model cannot predict in the data.
Traditionally for polynomial models, these standard deviation estimates are one way to decide to drop terms in a model. If the confidence interval includes zero as a possibility, then we might conclude that we cannot distinguish the provided estimate of a given parameter from zero.
At least this is true in a traditional model where we are trying to fit the correct model to our data, and all that obstructs this process is the presence of random (Gaussian) noise. The above approximations also assume the parameters are independent, so depending on the data, x and x^2 will usually not be orthogonal variables. That lack of orthogonality introduces correlation in the regression estimates, which is unfortunately ignored in the approximation used for the reported parameter standard deviations.
There are additional problems with this methodology when there is lack of fit in the model. Here we are trying to fit one model to the data, but in reality, the process is not truly based on that model. Our choice of model may be only an approximation to reality, or it may simply be wrong. Suppose you have a curved surface, and you attempt to estimate a linear model. In this case, there will be lack of fit, because the data is not well represented by the model.
Of course, the linear algebra does not understand that lack of fit is not just pure noise. Mathematics is a bit dumb in this respect. It cannot see a pattern in the residuals as we can. It simply believes the residuals are pure noise, and uses that noise to produce estimates of our uncertainty in the parameters using standard statistical techniques.
In the end, what I'm trying to say is these standard deviations are rough measures of how well the software thinks the parameters are estimated, based on the data you provided, based on several approximations, and based on the validity of your model for this data.
Raymond Tsang (view profile)
Hi John,
This is a great tool! I'm using it to fit a quadratic function to a magnetic field over a 3D lattice.
I'm interested in the uncertainty in the fit parameters. Is it correct to interpret "ParameterStd" as this uncertainty? (If not, what is?) And it is a little puzzling to me as to how "ParameterStd" is determined given that I have not provided the measurement uncertainty to the code. Can you explain a little bit about how this "ParameterStd" is computed?
Thank you!
Raymond
John D'Errico (view profile)
yonatan  Bound constraints on the coefficients are not one of the options. You can instead use a tool like lsqnonneg for basic nonnegativity constraints on the coefficients. For example, if x and y are column vectors then to get nonnegative coefficients for a and b as a vector, you would do this:
ab = lsqnonneg([x, x.^2],y);
For more complicated constraints on the coefficients, you would generally want a tool like lsqlin from the optimization toolbox. For example, if you wanted to constrain only one of the coefficients or you wanted constraints on the sum of the coefficients, then lsqlin would be useful. Lacking that toolbox, there are other tools like my own lee from the file exchange that can do some things for you along these lines, and with some effort you can always convert many general problems into one that lsqnonneg can handle.
yonatan (view profile)
This is Great!
Is there any way to constraint the function to be with positive coefficients only?
for example, f(x) = a*x + b*x^2 for positive a and b?
thank in advance.
Andrew (view profile)
This is great!
Syeda (view profile)
Chenyun Pan (view profile)
John D'Errico (view profile)
Martin  You would need to list the terms in the model. E.g., to specify a model that has total order no higher than 3, but no cubic terms in y, you might specify the modelterms argument as
modelterms = [3 0
2 1;
2 0;
1 2;
1 1;
1 0;
0 2;
0 1;
0 0]
Alternatively, you may do it by giving a list of terms as a string for that same model:
modelterms = 'x^3, x^2*y, x^2, x*y^2, x*y, x, y^2, y, constant'
Martin (view profile)
Hi John, first of all  great work!
My question: is it possible to use two different orders for my polynomial? Because i managed to fit my function (3D) with an polynomial of the order 6 (n=6). But it would be great to fit my function with a polynomial with the order e.g. 5 in one direction, and with 6 in the other one (m and n) .
z(x,y) = a(n)*x^n + b(m)*y^m + a(n1)*x^(n1) + b(m1)*y^(m1)...
John D'Errico (view profile)
Jeff  The reason why polyfitn returns a struct is because the model may be somewhat complicated, with general terms specified by the user in several variables. A 1d polynomial is far easier to describe as a list of coefficients, since the 1d polyfit always assumes a very specific model form.
However, it is no problem at all to extract the coefficients from the struct as returned. For example, consider this simple model z(x,y).
% Some simple data...
x = rand(10,1);
y = rand(10,1);
z = x + 2*y + 3 + randn(size(x))/10;
% Generate a model
mdl = polyfitn([x,y],z,{'x','y','constant'})
mdl =
ModelTerms: [3x2 double]
Coefficients: [0.9706 2.0162 3.0302]
ParameterVar: [0.0214 0.0266 0.0093]
ParameterStd: [0.1462 0.1631 0.0966]
R2: 0.9787
AdjustedR2: 0.9726
RMSE: 0.1180
VarNames: {'x' 'y'}
Of course, the coefficients are just numbers in a field of the struct. So we can extract the coefficients as a vector:
format long
C = mdl.Coefficients
C =
0.970599545975740 2.016179608143575 3.030225230230983
They are not exactly what I put in, but that random noise on the data corrupted things a bit.
Or we could even have converted the struct to a symbolic polynomial, either a sym, or my own sympoly form.
polyn2sympoly(mdl)
ans =
0.97059954597574*x + 2.01617960814358*y + 3.03022523023098
And we can evaluate the polynomial on our own without need of polyvaln, as long as you know what the model is. (And since the exponents for each term are stored in the struct, that is easy too.)
C(1)*.5 + C(2)*.25 + C(3)
ans =
4.019569905254746
For comparison:
polyvaln(mdl,[.5,.25])
ans =
4.019569905254746
Jeff (view profile)
John  wonderful work. Question: I think I've read that the coefficients solved for by polyfitn can only be interpretted by polyvaln. Is that correct? If I want to use polyfitn to find coefficients for a data fitting exercise, and then apply them to future datasets outside of matlab (and thus without polyvaln), is that not possible? Or, did I misunderstand something?
Daniel Higginbottom (view profile)
Eric H. (view profile)
Dennis Wouters,
Thanks for providing this fix. It works!
Dennis Wouters (view profile)
John, nice work.
I encountered a minor glitch.
When specifying the modelterms as an array of exponents, the variable 'varlist' was not assigned a value.
I fixed this by inserting the following at line 160:
else varlist = repmat({''},1,p);
John D'Errico (view profile)
Despite my having written polyfitn, high order polynomials (especially in multiple dimensions) are rarely a good substitute for a nonlinear model that is chosen to have the proper behavior. That is, IF you have such a model. Those polynomials tend to do nasty things to you, especially in extrapolation. And unless you have sufficient data, then the model may behave strangely even inside your data. I call this the problem of interpolation.
Basically, you need to be very careful with high order polynomial fits.
By the way, nlinfit DOES handle multiple independent variables. You suggest it does not, but there is no constraint on that. The problem is choosing an intelligent model in multiple dimensions. This is EXTREMELY difficult in general. One good reason is that a functional form like w = f(x,y,z) is so difficult to visualize.
Personally, my preference has always been to move to a splinelike model for these problems. I am often unable to pull a model out of thin air, and I would prefer to avoid a high order polynomial. My own gridfit is one simple example, with 2 independent variables.
Karan (view profile)
Hello John,
Thank you for your inputs. I really appreciate the help. From what i have understood is that the data is not sufficient enough for using such a high order polynomial fit. I am actually looking at fitting the data to my own custom equation but haven't found a model for it yet similar to nlin fit but with multiple independent variables. If you do have any suggestions please do let me know.
Thank you.
Karan
John D'Errico (view profile)
Hi Karan,
You've fallen into a common fallacy in regression modeling. It goes like this: "If n is good, then n+1 is better, and n+2 must be really wonderful."
The fallacy starts out easily. We try a linear model on our data. It works, all right, but not very well. After all, must problems are only pretty locally linear. So we try a quadratic model. It fits better. The lack of fit is reduced. But there are still problems in the model. It is inadequate to fit the shape of our data.
Since we did so well with the quadratic, try out a third order model. Yep, it does better still. Wow. Keep pushing things, until we have a 6th order model, or 10th order. A polynomial model is just like a Taylor series approximation, so more terms must give a better result. Eventually though, it all just turns to numerical crap before our eyes. So what happened?
Here you have 3 variables. A first order model has a constant term, as well as x, y, and z terms. So 4 terms in that model. You have 74 data points, and surely had no problem estimating the corresponding coefficients.
A quadratic model adds in 6 more terms, x^2, y^2, z^2, but also the cross product terms x*y, x*z, and y*z, so a total of 10 coefficients to estimate.
The curse of our fallacy is that as we move up in order, the number of terms grows faster and faster. At 3rd order, we have 20 terms to estimate, 35 at 4th order, 56 at 5th order, and 84 at 6th order. A 10th order model would have 286 terms.
How about your data? You have only 74 data points. (In my early days when I did some statistical consulting, my clients would always ask how much data they needed. My answer was almost always "More than you have.")
Every distinct data point is at most one additional piece of information. (Replicate points do not count here. All they do is reduce the uncertainty in your estimates. So replicates are good, but not quite as valuable as distinct points in some respects.)
Worse, new data points do not always provide new information. For example, two points determine a line. A third point on that line adds no new information towards the slope or intercept of that line. Again, like a replicate, all it does is reduce the uncertainty in your estimate of the coefficients.
So what happened when you tried to estimate a 6th order model with only 74 data points? You had too little information in your data to support such a large model. The message that comes back is from MATLAB, telling you that your matrix is numerically singular (rank deficient.)
Having said all this, things can get worse yet, making the problem yet more difficult. Some sets of data are worse than others. For example, if you add 1e6 to all of your independent variables, a polynomial of some given order will still be estimable mathematically (in theory), yet it may result in a matrix that is numerically singular. This is because the linear algebra is written using double precision floating point arithmetic. The numbers may simply lie over too wide of a dynamic range for the algebra to succeed. (What happens when you raise numbers that are on the order of 1e6 to successively higher powers?)
So when you have more coefficients to estimate than data points, you will always have a rank deficient matrix, but depending on the data itself, you may still have rank deficiency problems.
Some polynomial models will become more tractable to estimate if you center and scale the variables. Thus make each variable so it has a mean of zero, and lie roughly in the interval [1,1]. See that raising numbers that are on the order of +/1 to higher powers does NOT create a dynamic range issue. This will force you to transform your problem, but it often allows you to gain viable estimates for an otherwise intractable problem. (The use of orthogonal polynomials is another approach that often helps here.)
I hope this helps. Really, this is a subject that I've seen taught as a semester long course, so I have only touched the surface.
John
Karan (view profile)
Thank you so much for this. I am getting this warning
Warning: Rank deficient, rank = 64, tol = 5.8216e010.
I have used 3 independent variables, 74 data points and used degree of 6. Could you please help me understand what this means? Thank you.
John Mahoney (view profile)
John,
Great code as usual!
It seems to fail when given only one data point.
n = 1;
x = rand(n,2);
y = exp(sum(x,2)) + randn(n,1)/100;
p = polyfitn(x,y,3);
While failing might be desirable in some cases, it should really throw an error instead. And your code does supply an nparam output for n1 inputs, so it would not be inconsistent to handle this case.
I think it is in the logic after line 120  checking the size of the inputs and transposing accordingly.
Thanks!
wei (view profile)
Scott DeWolf (view profile)
Criz Andre Tamatoa Thomas Almonte Torrez Manlangit (view profile)
John D'Errico (view profile)
John  There are two possible questions here. First, if you have a problem with what is known as multiple right hand sides, then yes it is very vectorizable. By this I mean a problem with the same sets of independent variables, but different sets of dependent variables. Then backslash is the way to go, allowing you to factorize the equations once for a huge time savings. You would need to build the design matrix of course, but this is not a difficult thing to do.
If however, you have multiple problems, each with different sets of independent variables, then vectorization would require you to build a sparse block diagonal design matrix. As such, you can gain a bit in this effort, but not as significant of a gain as the first case. You need to build the design matrix as a sparse matrix. The trick here is to supply the blocks to blkdiag as individually sparse matrices. Thus, create the individual (sparse) design matrices in a cell array, then pass it to blkdiag as a comma separated list. Can you do this efficiently enough to make the backslash solve more efficient than a loop around the calls? Yes, but as I said, it will take some effort.
You can probably glean from my response that it will not be a trivial action to vectorize the result from this code, but that you could do things more efficiently with some time invested.
John (view profile)
Great set of functions!
Can you suggest any ways to vectorize their use? I have a scalar function that is dependent on its location in space f(x,y,z), and I calculate this function surrounding N bodies. Currently, I run a loop for each body, but is there a way to modify the code so that the backslash operator works across all of them at once, returning N sets of coefficients?
John D'Errico (view profile)
This should do it...
Look for these two lines below:
[Q,R,E] = qr(M,0);
polymodel.Coefficients(E) = R\(Q'*depvar);
Assume a column vector of weights as W, one for each data point. Change the above lines to read:
[Q,R,E] = qr(bsxfun(@times,W(:),M),0);
polymodel.Coefficients(E) = R\(Q'*(W(:).*depvar));
This should give you weighted estimates of the parameters. I think the statistics on the parameters will work too, as well as R^2, etc. Of course, you will need to pass in the weight vector too as an argument.
The above change assumes you are using a version that has bsxfun available. Its been around for a while so that should not be a problem. Even if not, a call to repmat will replace it.
Slaven (view profile)
Great function! Thank you. Would there be an easy way to convert it to do weighted least squares. I was looking to try to do it myself, perhaps you have some insights or warnings?
John D'Errico (view profile)
A linear regression (which is polyfitn) uses backslash at its core. For example, suppose one wished to use polyfit to estimate a quadratic model in one variable. The model would be
y = a1*x^2 + a2*x + a3 (+noise)
polyfit is able to estimate the three unknown coefficients [a1,a2,a3] by solving a linear (but overdetermined) system of equations in a way that minimizes the sum of squares of the residuals. Again, backslash is the engine underneath here.
Now, suppose one wished to estimate a more complex model, with two independent variables, z(x,y). I'll write it as a full two variable quadratic model, so there will be six unknown coefficients.
z = a1*x^2 + a2*x*y + a3*y^2 + a4*x + a5*y + a6 + noise
As long as you have sufficient data to estimate the six unknown coefficients, the problem is essentially the same. See that they enter the problem linearly in the same way as for the one variable problem we looked at before. Backslash again can solve the overdetermined problem to estimate the coefficients.
In fact, both polyfit and polyfitn use a QR factorization of the design matrix to achieve their goals. This is because they are both designed to also generate covariance matrix estimates on the parameters. A QR factorization is the easiest way to serve both purposes, but its all just linear algebra.
Mr D (view profile)
Thanks..your function solve my problem, but honestly I'm curious about the concept of linear regression in 2d surface. Do we have to do polyfit (1st order) on each row to get the best fit plane of the surface?
Thanks.