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
This submission has inspired the following:
gapolyfitn

• MATLAB release: MATLAB 7.0.1 (R14SP1)
• Published M Files: polyfitn_demo

Other Files

license.txt,
PolyfitnTools/demo/html/polyfitn_demo.png,
PolyfitnTools/demo/html/polyfitn_demo_01.png,
PolyfitnTools/demo/html/polyfitn_demo_02.png,
PolyfitnTools/demo/polyfitn_demo.m,
PolyfitnTools/polyfitn.m,
PolyfitnTools/polyn2sym.m,
PolyfitnTools/polyn2sympoly.m,
PolyfitnTools/polyvaln.m,
PolyfitnTools/ReadMe.rtf,
PolyfitnTools/test/test_main.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 ...

[Thomas Clark]

Fit for it's purpose ;)

Thanks John, I had to chop and modify polyvaln for my application, but it was really good to have a way of generating the polynomial coefficients for a given n and order.

[Chao Xu]

Hallo John,
Nice work!
I am a newbie in matlab and still got some questions about 3d fit.
I have a 3d dataset F(128,128,16) and would like to get the polynomial function F(x,y,z), how should I do it?

thanks!

[John D'Errico]

Assuming that you have an array, that contains only w(x,y,z) for each element of the array, use ndgrid to generate the combinations of (x,y,z). Then call polyfitn with the data.

[Andreas]

Great work !
Thank you very much. It is really easy to use, allows a lot of user specifications and is quite fast imho.

[James]

Hi John,

I am new to all this stuff. actually to MATLAB itself.
I have 2D data I want to fit (to a 2-variable function) with 3 parameters: z = A + B*x +C*y. But that is scarce matrix, i.e. with holes. I do not have data (or z if you prefer) for each and every combination of x and y, but only for a sall number of them.

Any way around this ?

I appreciate your help in that matter.

James

[John D'Errico]

James,

You don't need every possible combination of (x,y,z) to be known to estimate a model, certainly not a first order model like this. To estimate the 3 unknown parameters in this model, in theory, you need at least 3 data points at a minimum. In practice, you won't get very good estimates of those parameters with only 3 data points.

In the past, I often had clients ask me how much data they needed to obtain. After all, data is often scarce, expensive, hard to obtain. It costs money, and time. You can't have too much data, at least if it is good data.

In practice, if you end up with 10 or 20 points to estimate these three parameters, that will (probably) be entirely adequate. More data will give you higher quality estimates. Even 6 points, carefully chosen, if the noise is low enough, will often be adequate. Only you know the quality of your data, and how much you have, as well as how accurately you need the model to be estimated.

If you have further questions, I can only suggest you send me an e-mail, with your data, and any other information you have on the problem.

John