Polyfitn

N-d polynomial regression model

Description

Yes, I know. This is yet another in the long list of polynomial regression modeling tools. What does polyfitn do that the others don't? You can supply any list of terms, with any set of exponents: positive, negative, or real. It works in n-d as well as 1-d. And it has a simple, clean interface to specify your terms.  
 
Simple statistics are generated, R^2, rmse, standard deviations and variances of the coefficients. I've supplied an evaluation tool to use the model.  
 
x = -2:.1:2;  
y = cos(x);  
p = polyfitn(x,y,'constant x^2 x^4 x^6');  
 
% Evaluation  
polyvaln(p,[0 .5 1])  
ans =  
      0.99996  
       0.8776  
      0.54031  
 
% Conversion to a sympoly  
polyn2sympoly(p)  
A scalar sympoly object  
    0.99996 - 0.49968*x^2 + 0.041242*x^4 - 0.0012079*x^6  
 
This last operation employs my own sympoly toolbox also found on the file exchange, although it is not necessary for use of polyfitn and polyvaln.  
 
http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=9577&objectType=FILE  
 
As another example, fit a complete 2-d model with all terms up to a second order model:  
x = randn(100,2);  
y = randn(100,1);  
p = polyfitn(x,y,2);  
 
% Display it as a sympoly, just makes it easy to read:  
polyn2sympoly(p)  
 
A scalar sympoly object  
    0.19175*X1^2 + 0.052769*X1*X2 + 0.093648*X1 + 0.039723*X2^2 - 0.04943*X2 - 0.19188  
 
% Only want a linear model in 2-d, but with no constant term? Its easy to specify:  
uv = rand(100,2);  
w = sin(sum(uv,2));  
p = polyfitn(uv,w,'u, v');  
 
polyn2sympoly(p)  
A scalar sympoly object  
    0.74604*u + 0.68452*v

Acknowledgements

The author wishes to acknowledge the following in the creation of this submission:
Symbolic Polynomial Manipulation

• MATLAB release: MATLAB 7.0.1 (R14SP1)
• Published M Files: polyfitn_demo
• Other Files: PolyfitnTools/ReadMe.rtf, ; PolyfitnTools/test/test_main.m, ; PolyfitnTools/demo/html/polyfitn_demo.png, ; PolyfitnTools/demo/html/polyfitn_demo_01.png, ; PolyfitnTools/demo/html/polyfitn_demo_02.png, ; PolyfitnTools/polyvaln.m, ; PolyfitnTools/polyfitn.m, ; PolyfitnTools/polyn2sym.m, ; PolyfitnTools/polyn2sympoly.m, ; PolyfitnTools/demo/polyfitn_demo.m

[Duane Hanselman]

Powerful general tool for N-D curvefitting. Great extension to polyfit/polyval. Flexible input syntax. Scales data to minimize numerical problems.

[Vassili Pastushenko]

I wanted to compare the speed of my SURFIT with that of POLYFITN - POLYVALN. My problem is that POLYFITN apparently accepts only a vector as dependent variable, whereas I have to get a polynomial approximation of a matrix. How should I do that?  
 
An example from my work: imagine that I have a data matrix Z (512 by 511) which I want to approximate by a polynomial of 2-nd degree in vertical direction and 0-th degree in horizontal direction. This means I want to use a bipolynomial of degrees [2 0]. With SURFIT the problem is solved easily:  
 
B=surfit(Z,[2 0]); %B = my bipolynomial approximation of Z (regression)  
 
What should I do with POLYFITN-POLYVALN to get B?  
 

[Vassili Pastushenko]

Finally I have succeeded to compare the speed, which is surprisingly almost 40 times slower by POLYFITN-POLYVALN than by SURFIT. I cannot understand this because I know how low does John rating other programs if they are with some examples even 2% slower than possible. Here is my comparison made for a bipolynomial of 2-nd degree in both directions (I was not able to resolve my case of degrees [2 0]). Much more surprising is higher accuracy of SURFIT. I should vote for deletion of POLYFITN were I sure that my syntax of input to it indeed corresponds to bipolynomial of degrees [2 2].  
 
tic,[X,Y]=meshgrid(1:size(Z,2),1:size(Z,1));  
mod=polyfitn([X(:),Y(:)],Z(:),2);  
B=polyvaln(mod,[X(:),Y(:)]);  
B=reshape(B,size(Z,1),size(Z,2));toc  
Elapsed time is 1.203000 seconds.  
 
tic,B2=surfit(Z,2);toc  
Elapsed time is 0.031000 seconds.  
 
sum((Z(:)-B(:)).^2)  
ans = 76181  
 
sum((Z(:)-B2(:)).^2)  
ans = 75534

[Vassili Pastushenko]

Yes, the comparison was done correctly. I have noticed that the ratio of errors depends on data. Occasionally, I used real data for comparison, which allowed to notice an inacceptably big difference in accuracy. Below I have included a model of a piece of these data as default.  
   
 
[B,B2]=john;  
Elapsed time is 0.063000 seconds. %POLY...N  
Elapsed time is 0.000000 seconds. %SURFIT  
 
JV = 1.1931e+005 95907 %Errors  
 
 function [B,B2]=john(Z)  
%======================  
if nargin<1  
y=zeros(100,160);  
y(1:30,20:60)=-8;  
y(80:100,40:120)=-8;  
Z=y+randn(size(y));  
end  
 
[R,C]=size(Z);  
[X,Y]=meshgrid((1:C)/C,(1:R)/R);  
z=[X(:),Y(:)];  
ZZ=Z(:);  
tic,  
mod=polyfitn(z,ZZ,2);  
B=polyvaln(mod,z);  
toc  
B=reshape(B,R,C);  
 
tic,B2=surfit(Z,2);toc  
JV=sum([ZZ-B(:),ZZ-B2(:)].^2)

[Tony Flipper]

did what I needed it to do - find the same regression equation that Excel gave me

[g hunter]

does what i need it to.  
 
note that Vassili Pastushenko wrote SURFIT.

[ye zhang]

Thanks a lot!!  
it is really helpful commmand...

[Jean-Luc Dellis]

Hi all,  
 
Good job John as usually you do. I planed to use bootstrap approach for some educational experimental works because uncertainties can touch both x and y variables. In that goal i compared your program with my own bootstrap program using polyfit. That gave similar results except a) as expected bootstrap was slower, but that is not very important in my case, b) bootstrap failed when data points were under around 10 with a polyfit warning message stating an ill-conditioned matrix. This is due to replicated points in the bootstrap samples.  
 
Maybe i could manage bootstrap in a more correct way. Anyway your program succeed more quickly and with reasonable uncertainties even when there was noise in the independant variable x.  
 
So, thanks a lot John. I know that if i have some clever task to do, i before have to check if John has done it.

[Oh-ig Kwoun]

This is exactly what I needed.  
It works great.  
Thank you, John

[Simon Ren]

This is excellent work. Thanks John.  
I have a question. I have a series points (x, y,z), but no any equation yet. These data was obtained from test. I am not sure if I can get the function z=f(x,y) with polyfitn? Could you help me? Thank you very much.  
 

[guido spinola]

 

[Rajneesh Kumar]

It's very easy to use and works excellent. Has been very useful to my research. Thank you.

[Vincent Morio]

Very powerful and efficient tool. I have used it to approximate the nonlinear entries of my linearized model by multivariate polynomials. Then, it is rather easy to put it under LFT form in order to synthesize an LPV controller. Many thanks.

[Muhammad Ihsan Salim]

Well, I need to try it out first, so i gave it four star. If it works on fitting my DIC displacement distribution, I gave you the whole star you want. But after all, thank you any way..

[Rebecca Bendick]

When I try to use polyfitn, I get lots of errors and no answer:  
 
??? Error using ==> class  
The CLASS function must be called from a class constructor.  
 
Error in ==> sympoly at 92  
    P=class(P,'sympoly');  
 
Error in ==> horzcat at 12  
  sp=sympoly(sp);  
 
Error in ==> repmat at 45  
    siz = [M N];  
 
Error in ==> polyfitn>buildcompletemodel at 241  
    modelterms = [modelterms;[repmat(k,nt,1),t]];  
 
Error in ==> polyfitn at 148  
  modelterms = buildcompletemodel(modelterms,p);

[John D'Errico]

This error is incomplete, since it does not show the calling syntax. As such, I cannot say why the problem arose. Please send me a direct e-mail with any questions, and I will resolve them immediately.  

[Richard Crozier]

 

[Richard Crozier]

Hello, the software is excellent, but there is a very minor error, it does not check for that the modelterm matrix passed in is not empty, resulting in an mtimes error in this case. I discovered this because of my poor coding, so I suppose I shouldn't expect polyfitn to check this for me!

[John D'Errico]

Sorry. Polyfitn should check this for you. It will do so in my next revision.

[Mohsen Nosratinia]

Many thanks! just what I needed. saved me a lot of coding

[Gavrilo Bozovic]

Thanks a lot!

[Agustin Pimstein]

John hello, 
 
Even running this example of yours, 
 
n = 1000; 
x = rand(n,2); 
y = exp(sum(x,2)) + randn(n,1)/100; 
p = polyfitn(x,y,3) 
 
 
...I'm getting this error: 
??? Undefined function or method 'polyfitn' for input arguments of type 'double'. 
 
Any idea of what the problem might be??? 
 
Thanks,  
 
Agustin 

[John D'Errico]

Agustin found his problem - this is always a path problem. Something got put where MATLAB does not look.

[Michal Kvasnicka]

John, could you add to the polyfitn package description file mentioned in your ReadMe file or any relevant reference on multidimensional polynomial regression.  
 
"In additions to these m-files, there are several other files of interest. For those who have an interest in understanding how polyfitn works (quite simple really) you should look to the file understanding_polyfitn.rtf, in the doc subdirectory."  
 
The underlying theory will be very helpful for me. Thanks in advance ...