Editor's Note: This file was selected as MATLAB Central Pick of the Week
Sometimes we want to fit a polynomial to some data, but want to force a perfect match to some known values at one or more points. This function returns such a polynomial.
Are Mjaavatten (2019). polyfix(x,y,n,xfix,yfix,xder,dydx) (https://www.mathworks.com/matlabcentral/fileexchange/54207polyfixxynxfixyfixxderdydx), MATLAB Central File Exchange. Retrieved .
1.3.1.0  Fixed the error in version 1.3.0.0 pointed out by Jacob Bresler. 

1.3.0.0  As requested by Jacob Bresler, I have added optional outputs S and mu, to make the output from polyfix conform with polyfit 

1.2.0.0  Fixed trivial errors in the help section 

1.1.0.0  Version 1.1 includes an option for specifying derivatives. p is now a row vector, for consistency with polyfit 
Create scripts with code, output, and formatted text in a single executable document.
GuillyT (view profile)
Jacob Bresler (view profile)
Are Mjaavatten (view profile)
Jacob: Thanks for pointing out the errors. I have implemented your corrections (and added scaling of dydx). Obviously, I should have tested more thoroughly.
Jacob Bresler (view profile)
I went through the code, and the only changes you should need to make is the addition of the following lines
xfix = (xfix  mu(1))/mu(2); %this should be added in the new if nargout > 2 loop
also the following lines need to be added in the if nargin > 5 loop
if nargout > 2
xder = (xder(:)mu(1))/mu(2);
else
xder = xder(:);
end
Jacob Bresler (view profile)
thanks for the update, however the mu function does not appear to be properly applied to the fixed x values
Jacob Bresler (view profile)
could you add a couple lines to deal with the other standard outputs from polyfit such as error and the centering of the x values (mu=[mean(x);std(x)], x_c=(xmu(1))/mu(2)). Otherwise this is a great script.
Are Mjaavatten (view profile)
Jesse:
R squared will be negative whenever y = y_mean is a better model than your polynomial fit. This may happen if you fix a point that is not representative for the observed values. Example:
x = linspace(0,1,100)';y = 1 + .1*randn(100,1);
p = polyfix(x,y,2,0,0);f = polyval(p,x);
plot(x,y,'.',x,f);
Rsquared = 1((yf)'*(yf))/((ymean(y))'*(ymean(y)))
Jesse Kadosh (view profile)
Are,
Many thanks for this useful and efficient code!
Reading about Rsquared values, I see that they can be greater than 1 even for linear regression with loworder polynomials in cases where a yintercept is fixed at 0 (i.e., regression through the origin, which occurs when a constant is not included among the regression functions). I realize that polyfix (which uses polyval) probably always includes a constant term in the regression functions, but I am wondering if the mathematical approach used in polyfix to force the fixed points in the fit has a similar effect to omitting the constant term. That is, if we use polyfix to force the curve to pass through a chosen point, is that equivalent mathematically to "forcing" a yintercept value of zero (passing through the origin), as far as computing Rsquared is concerned? If so, I assume it is then possible for Rsquared values to surpass 1 when using polyfix and forcing one or more fixed points.
Some reference posts regarding Rsquared surpassing 1 when the yintercept is omitted (fixed at zero) are:
https://www.mathworks.com/matlabcentral/answers/159369polyfitandr2value
https://stats.stackexchange.com/questions/219810/rsquaredandhigherorderpolynomialregression
https://stats.stackexchange.com/questions/334004/canr2begreaterthan1
Amal Mansoor (view profile)
This is so helpful! I spent hours trying to figure this out. I was just about to give up and use excel (which is horrendous for large data sets) when I found this. Thanks!
Are Mjaavatten (view profile)
Avichay: For most of my own use of polyfit or polyfix, visual inspection is more appropriate than formal measures of the goodness of fit. Most of those measures are easily calculated from the inputs and outputs anyway, so I prefer to to leave that to the user.
Here is an example of how to calculate R^2::
x = linspace(0,pi,25);
y = sin(x);
p = polyfix(x,y,2,0,0);
f = polyval(p,x);
Rsquared = 1((yf)*(yf)')/((ymean(y))*(ymean(y))');
Avichay Efraim (view profile)
Great function, thanks!
Is there a way to calculate R^2 of the output fit curve?
Are Mjaavatten (view profile)
HeadInTheClouds:
I had a feeling that I might have misunderstood your question. I believe that what you want will require a bounded, nonlinear optimization that is best done outside of polyfix.
In the example below, nrm is an inline function that calculates the norm of the differences between the y values and the corresponding values of the best fit polynomial passing through(x0,0). I use fminbnd from the optimization toolbox to find the value of x0 within the given interval that minimizes nrm(x0). Given the spread of the y values, a straight line seems most appropriate in this case.
x = linspace(0,2,100)';
y = exp(x)0.25+ 0.3*randn(100,1);
nrm = @(x0) norm((ypolyval(polyfix(x,y,1,x0,0),x)));
x0 = fminbnd(nrm,1.3,1.6);
disp(x0)
p = polyfix(x,y,1,x0,0);
plot(x,y,'.',x,polyval(p,x),[1.3,1.6],[0,0],'k');
HeadInTheClouds (view profile)
Are Mjaavatten,
Thank you for the prompt response. My question was poorly worded which resulted in a misunderstanding. It should have been worded, "Is it possible to use this function in order to fit a polynomial through a specified yvalue and an xvalue within a range of specified xvalues?" By this, I mean that I want to force a polynomial to cross the xaxis but do not want to specify a particular xvalue in which it crosses the xaxis, I'd like to specify a range of possible xvalues in which the polynomial is fit through one point within the specified range.
Hope this helps clear up my question.
Are Mjaavatten (view profile)
HeadInTheClouds: If a polynomial is constant throughout an interval it must be constant everywhere, and you do not need polyfix to tell you that. In your case: p(x) = 0. Or did I misunderstand your question?
HeadInTheClouds (view profile)
Hi Are Mjaavatten,
Is it possible to use this function in order to fit a polynomial through a range of values instead of through a discrete point (e.g. if I want to force a function through y=0 between x = x1 and x = x2)?
Thanks in advance!
Vogel (view profile)
Great function!
Does it exist for 3D?
Giuseppina Buttitta (view profile)
Are Mjaavatten (view profile)
Answer to Md. Nahidul Islam:
Sorry I did not see your comment before now. I assume that what you want is to stretch and shift the original curve to pass through the specified points. Let p be the original polynomial and xfix and yfix the coordinates of the specified end points.Then the new polynomial is given by
scale = diff(yfix)/diff(polyval(p,xfix));
pnew = p*scale; % stretch
pnew(4) = pnew(4) + yfix(1)polyval(pnew,xfix(1)); % shift
Md. Nahidul Islam (view profile)
Hi Are Mjaavatten,
Suppose I have an equation
ynew = p1.*xnew.^3 + p2.*xnew.^2 + p3.*xnew + p4; %equation 1
where,
p1 = 4.095e15;
p2 = 3.3606e10;
p3 = 1.164e05;
p4 = 2.9954;
If
xnew=[5:5:22385]; then we have a curve.
I want to get a new curve, which will follow exactly the same trend as equation 1, but will be the start and end point as [5, 22385], [15, 14.7]
How can i do that?
Thank you.
Victoria Duke (view profile)
Fantastic code, would be very helpful to have for 3D!
DanaAdriana Botesteanu (view profile)
burges (view profile)
The code works perfectly. How about extending this to 3 dimensions?... Many thanks!
Dominik Vít (view profile)
Are Mjaavatten (view profile)
Answer to chef13:
It is certainly A tool, and you have already used it, I see. Whether it is an appropriate tool will depend on the broader context of your task.
chef13 (view profile)
Hi I want to solve this problem:
https://www.mathworks.com/matlabcentral/answers/327332generatetrajectorieswithsametangentinaspecificpoint
Do you think that this is the right tool ?
Thanks
Are Mjaavatten (view profile)
Answer to Anelechi Ibekw:
The file has been tested on versions up to R2014b, which is the latest to which have access.
burges (view profile)
Please, what version of Matlab has this function? I actually thought it was R2016b, but it isn't.
SEBIN JOSE (view profile)
Thank you