Note:
The example fitting will throw a warning in R2009a but will not in later releases due to changes in how POLYFIT throws warnings. The numerical properties of the problem have not changed, it is still ill-conditioned.
The warning message occurs is due to the scaling of the data. It produces an ill-conditioned least squares problem which POLYFIT flags. One way to address this is to scale the data, such that each "X" in the original input is scaled into "X_new" as follows:
X_new = (X – mean(X))/std(X)
This scales the input values so that x' is contained in [-1 1]. The third calling syntax of POLYFIT does this scaling before fitting the polynomial and returns a third output:
mu = [mean(X) std(X)]
This will return a different set of coefficients than as if the data were to be fitted with no scaling.
you will receive:
%scaled polynomial
p*(X_new) = A_new*(X_new)^2 +B_new*(X_new) + C_new
%For notational purposes, the unscaled polynomial is
p(X) = A*(X)^2 +B*(X) + C
These new coefficients are correct in the sense that p(X_new) will produce the same Y value as its corresponding X and the fitting can be used to estimate any query point if it is scaled properly. These scaled coefficients can be transformed back to their non-scaled form with some caveats which I will address later. Since these two equations must be equal for corresponding X and X_new values, there is a relationship between the coefficients for a scaled (A_new,B_new,C_new) and un-scaled (A,B,C).
A_new*(X_new)^2 +B_new *(X_new) + C_new = A*(X)^2 + B*(X) + C
By writing X_new in terms of X as above, this relationship can be manipulated to generate closed form expressions for A,B, and C:
[p,S,mu] = polyfit(X,y,2)
A_new = p(1); B_new = p(2); C_new = p(3); mu1 = mu(1); mu2 = mu(2);
A = A_new/(mu2)^2
B = (B_new*mu2)/mu2 – (2*A_new*mu1)/(mu2)^2
C=A_new*(mu1)^2/(mu2)^2-(B_new*(mu1/mu2))+C_new
However, this approach can be generalized for higher order fittings using the binomial theorem. (https://en.wikipedia.org/wiki/Binomial_theorem ). Consider the closed form equations above, but in matrix form, with A, b, and C being the sums of the columns. \n\n
M = [nchoosek(2,0) nchoosek(2,1)*(-mu1) nchoosek(2,2)*(-mu1)^2;
0, nchoosek(1,0), nchoosek(1,1)*(-mu1);
0, 0, choosek(0,0)];
T = ( p.' ./ [mu2^2 mu2 1].' ) .* M
This should hopefully make it a bit clearer that the rows of the matrix M are the coefficients of the expansion of:
(X - mu1)^k
The general form of which is given by the Binomial Theorem.
Lastly, a couple of caveats; This is not a recommended workflow to manually convert from scaled to non-scaled coefficients. Manually transforming the data is largely is redundant as you can produce the un-scaled coefficients by simply calling POLYFIT with the 2 output syntax in addition to the 3 output syntax. However, please note that these results can differ slightly from the manually transformed coefficients. See the well-conditioned example from the documentation below using the generalized matrix form of the transformation to demonstrate.
>> load census
>> p = polyfit(cdate,pop,2);
>> [phat,S,mu] = polyfit(cdate,pop,2);
>> T = ( phat.' ./ [mu(2)^2 mu(2) 1].' ) .* [nchoosek(2,0) nchoosek(2,1)*(-mu(1)) nchoosek(2,2)*(-mu(1))^2; 0 nchoosek(1,0) nchoosek(1,1)*(-mu(1)); 0 0 nchoosek(0,0)]
T =
1.0e+04 *
0.0000 -0.0025 2.3366
0 0.0001 -0.2298
0 0 0.0062
>> pT = sum(T) % sum over cols
pT =
1.0e+04 *
0.0000 -0.0024 2.1130
>> p
p =
1.0e+04 *
0.0000 -0.0024 2.1130
>> p - pT
ans =
1.0e-09 *
-0.0000 0.0008 -0.7785
Note that even in this relatively nice example there are some numerical differences since the coefficients have been calculated differently. Now let's turn our attention back to the original example
>> x = -100000:10000:200000;
f = @(x) 1000*x.^2+ 3*x + 0.023;
y = f(x);
>> p = polyfit(x,y,2);
[phat,S,mu] = polyfit(x,y,2);
T = ( phat.' ./ [mu(2)^2 mu(2) 1].' ) .* [nchoosek(2,0) nchoosek(2,1)*(-mu(1)) nchoosek(2,2)*(-mu(1))^2; 0 nchoosek(1,0) nchoosek(1,1)*(-mu(1)); 0 0 nchoosek(0,0)]
T =
1.0e+12 *
0.0000 -0.0001 2.5000
0 0.0001 -5.0000
0 0 2.5000
>> pT = sum(T)
pT =
1.0e+03 *
1.0000 0.0030 0.0000
>> p
p =
1.0e+03 *
1.0000 0.0030 0.0000
>> p - pT
ans =
0.0000 -0.0000 0.0025
In this case you will see that the numerical difference has increased by nearly 10^7. For similar situations, it is likely that the newly calculated non-scaled coefficients will not be an accurate fit for your data, and all you have done is essentially circumvented the warnings for an ill-conditioned problem.
