Code covered by the BSD License  

Highlights from
polyfitZero

4.66667

4.7 | 6 ratings Rate this file 100 Downloads (last 30 days) File Size: 47.1 KB File ID: #35401
image thumbnail

polyfitZero

by

 

01 Mar 2012 (Updated )

Fit polynomial to data, forcing y-intercept to zero or arb. value and slope to zero or arb. value

| Watch this File

File Information
Description

[update 2014-02-19: fix error in polyfitBM, add root and select-terms tools]
This submission contains four convenience polynomial fitting functions similar to POLYFIT.
1. POLYFITZERO - fit polynomial to data, forcing y-intercept to zero.
2. POLYFITB - force y-intercept to "b".
3. POLYFITB0 - force y-intercept to "b" and slope at (0,b) to zero.
4. POLYFITBM - force y-intercept to "b" and slope at x=0 to "m", e.g.: dy/dx = m.
5. POLYFITBROOT - force intercept and root
6. POLYFITBMROOT - force intercept, slope and root
7. POLYFITBMROOTTERMS - force intercept, slope, root and terms

If you get to POLYFITBMROOTTERMS please just use POLYFITN by John D'Errico:
http://www.mathworks.com/matlabcentral/fileexchange/34765-polyfitn

Forcing the y-intercept to zero is accomplished by noting for polynomial p = [aN, ..., a3, a2, a1, a0],

IE: y = aN*y^N + ... + a3*y^3 + a2*y^2 + a1*y + a0

when x is zero, then y is the constant term, "a0". Therefore a0 = 0, or in the general case when forcing the y-intercept to an arbitrary value, a0 = b.

Forcing the slope at x=0, is accomplished by noticing that the derivative of p

IE: dy/dx = N*aN*y^(N-1) + ... + 3*a3*y^2 + 2*a2*y + a1

evaluated at x = 0 yields "a1". Therefore a1 = 0 for zero slope, or in the general case when forcing the slope to an arbitrary value, a1 = m.

Required Products MATLAB
MATLAB release MATLAB 7.13 (R2011b)
Tags for This File   Please login to tag files.
Please login to add a comment or rating.
Comments and Ratings (12)
12 Mar 2014 Lucia Morales-Rivas  
02 Dec 2013 Subhashree Mishra  
17 Oct 2013 Zoe  
23 Feb 2013 giga  
30 May 2012 Mark Mikofski

John, Sorry for my slow response. You can accomplish this the same way I zeroed out the last coefficient.

E.g. to force the 1st order coefficient to zero do the following:

% ...
x = x(:);
y = y(:);

% since you want might want to include the constant term,
% change "z" to be a matrix of ones, and make it one column wider
% to include the 0th-order term

% z = zeros(dim,degree); % comment out this line which doesn't include constant
z = ones(dim,degree+1); % add this line which now includes 0th-order term

for n = 1:degree
z(:,n) = x.^(degree-n+1);
end

% you don't have to calculate the constant term,
% because it will be a column of ones (i.e. x^0 = 1)
% now remove any columns that you don't want to fit
% e.g. we want to force the 1st order coefficient to zero
orderOfZeroCoeffs = 1; % a vector of terms to omit from regression, e.g. [1,0],
% you could make "orderOfZeroCoeffs" an input to the function

% now, remove the terms whose coeffs you want to zero
orderOfNonZeroCoeffs= ~ismember(degree:-1:0, orderOfZeroCoeffs);
z = z(:,orderOfNonZeroCoeffs); % replace z with the columns removed
% there are a 1000 other ways to write this code!

% comment out the following lines
% because now it's more complicated
% pZero = z\y;
% pZero = [pZero;0]';

% instead create a vector of zeros
% any terms not calculated are already set to zero
pZero = zeros(1,degree+1);
% substitute the coefficients directly into "pZero" by index assignment
pZero(orderOfNonZeroCoeffs) = z\y;

% ...

Try that out. I'm curious what you're using this for. Also check out the other polyfit0 (fex272) function that Jennifer Sanders found above. You can modify it in exactly the same way I show here, but you will also have to modify the other output parameters.

24 May 2012 John

I wonder if you could modify the routine so that higher powers could also be forced to be zero. For example, the current routine set the constant coefficient to zero. But, can the 1st order coefficient also be forced to be zero. And so on...

09 Apr 2012 Jennifer Sanders

Gahh! It didn't post my initial message which was asking how to generate the structure S as you can ordinarily so with polyfit. Thus allowing you to calculate the error on the gradient. And then I found a file that does just that on the file exchange (See above post: which it posted twice?!).

09 Apr 2012 Jennifer Sanders

Oh, I just found there is an m-file that does what I want; polyfit0.m, found at:

http://www.mathworks.com/matlabcentral/fileexchange/272

Thanks anyway,
Jen

09 Apr 2012 Jennifer Sanders

Oh, I just found there is an m-file that does what I want; polyfit0.m, found at:

http://www.mathworks.com/matlabcentral/fileexchange/272

Thanks anyway,
Jen

02 Mar 2012 Mark Mikofski

Thanks John - You are absolutely right! I originally wrote polyfitZero for a coworker who wanted to constrain the y-intercept to zero as in Excel, but later when I realized that I could use conv(), I liked it for the silly reasons that I rarely get to divide polynomials or use convolution. But you were right to point out the dangers of dividing by zero, and I should have been more explicit in the comment above that polyfitZero_v2 is not suited for functions that cross the y-axis, for all the reasons you pointed out. I wish I could delete that comment. Thanks again.

02 Mar 2012 John D'Errico

Sigh. I don't like to do this, but Mark's comment is simply a terribly bad idea. It might be correct in the absence of error. It might make sense if you have no points where x is actually near zero. But it is a terrible idea in any other case. Please, don't follow his advice!!!!!

What happens when you do have an exact zero in x? Yes, you get a divide by zero. But even for a near zero in any element of x, what happens is you are essentially doing a weighted regression, where the point with a near zero gets a huge weight. This will be the only point that matters in the regression analysis.

The fact is, a polynomial fit with no constant term is trivial to accomplish by a variety of means. This code uses a reasonable method, that treats the data properly and avoids risk of divide by zeros for normal data.

Again, please avoid following Mark's advice.

02 Mar 2012 Mark Mikofski

There is also an even easier way to do this. Just divide by x first and then call polyfit, then use conv to multiply the 2 polynomials

function pZero2 = polyfitZero_v2(x,y,n)
q = y./x;
qnozero = q(x~=0);xnozero = x(x~=0);
pZero2 = polyfit(xnozero,qnozero,n-1);
pZero2 = conv([1 0],pZero2);
end

Updates
16 Oct 2013

output delta's for error estimation using polyval. add example

19 Feb 2014

add fits that force b to arb value and force slope at (0,b) to arb value.

19 Feb 2014

fix error in polyfitBM, add root and select-terms tools

Contact us