How do I rescale the polyfit coefficients back to normal coordinate system?

I fit some data with 1) [p,S,mu] = polyfit(x,y,1); These polynomial coefficients are different from those you would get from 2) [p,S] = polyfit(x,y,1); I need the coeffecients that correspond to a fit on my original, unscaled coordinate system. However, I want to use the 1st equation to benefit from the centering and scaling. How do I use the 1st equation (with mu) and obtain coefficients that correspond to my unscaled system?

Answer 1

[Frank Fournel]

Hello, i've found a way to do this : [pp(i,:),~,mu(i,:)]=polyfit(x1,y(i,:),order);%Fit par renormalisation normx(i,:)=[1/mu(i,2) -mu(i,1)/mu(i,2)]; r(i,:) = sym2poly(subs(poly2sym(pp(i,:)),poly2sym(normx(i,:)))); yy2(i,:)=polyval(r(i,:),x1);

but I would like to have a stand alone application and it is not possible to compile symbolic function. Have you found another way to do this?

[Neil Campbell]

I guess-and-checked and found a dumbed-down version of what you did. Thanks for the answer!

[James Thompson]

What did you end up doing? I am having trouble doing exactly this and also don't want to use symbolic functions.

Answer 2

[Jannik M.]

If phat are the output coefficients of polyfit in standardized x units, and mu = [mean(x) std(x)] the standardization values, you can get the coefficients p in units of x like this:

[phat, ~, mu] = polyfit(x, y, n)

phat = flip(phat); % Flip order to [p0, ..., pn]
p = zeros(size(q));

for i = 0:n
    for k = i:n
        p(i+1) = p(i+1) + nchoosek(k, k-i) * phat(k+1)/mu(2)^k * (-mu(1))^(k-i);
    end
end

p = flip(p); % Back to original order [pn, ..., p0]

Background

The standardization polyfit applies to the x-values is

,

where (mu(1)) and (mu(2)) are the mean and standard deviation of the original x-values. Express the polynomial function in terms of and the standardized coefficients and use the binomial formula, such that

If you write this out explicitly (e.g., for ) and regroup the prefactors by the power of x, you will see that you can express in terms of p,

using the transformation:

which is the formula applied in the code.

[wonsoek lee]

what is the q ? in the 3rd line

[Al]

I used p and it worked for me. Note you will need to plot anything taken from poly2sym since it buggy and randomly gives bad results. If using the equation in a SPICE code or other program that needs the equation double check, it will save you a lot of headaches. If it is wrong try changing the polyfit to a different order fit and re-check.

[Steven Lord]

Note you will need to plot anything taken from poly2sym since it buggy and randomly gives bad results.

Do you have a concrete small example of a random "bad result" from poly2sym?

Answer 3

[Faruk Telibecirovic]

I also stumbled on this problem and find relief in @Jannik M.code. Precision can be improved by changing this line:

p(i+1) = p(i+1) + nchoosek(k, k-i) * phat(k+1)/mu(2)^k * (-mu(1))^(k-i);

to fuse division of large powers of 'mu':

p(i+1) = p(i+1) + nchoosek(k, k-i) * phat(k+1) * ((-mu(1))^(k-i) / mu(2)^k);

If further precision is needed, without use of symbolic functions, I implemented some Dekker's double-double math:

function  [p] = rescalePhatDekker(phat,mu)
n = numel(phat) - 1;
phat = flip(phat); % Flip order to [p0, ..., pn]
p = zeros(size(phat));
for i = 0:n
    po_s = [0.0,0.0];
    for k = i:n
        nck_s = split(nchoosek(k, k-i));
        ph_s = split(phat(k+1));
        nckph_s = dd_mul(nck_s, ph_s);
        m1_s = split(-mu(1));
        m1p_s = dd_pow(m1_s,(k-i));
        m2_s = split(mu(2));
        m2p_s = dd_pow(m2_s, k);
        m1d2p_s = dd_div(m1p_s,m2p_s);
        rowadd_s = dd_mul(nckph_s, m1d2p_s);
        po_s = dd_add(po_s,rowadd_s);
    end
    p(i+1) = po_s(1) + po_s(2);
end
p = flip(p); % Back to original order [pn, ..., p0]
end

function xs = split(x)
% Splits double x into two non-overlapping components
t = 134217729 * x;  % 2^27 + 1
a = t - (t - x);
b = x - a;
xs = [a, b];
end

function r = mul12(a, b)
% Multiplies double a and b, returning the result as two non-overlapping components
as = split(a);
a1 = as(1);
a2 = as(2);
bs = split(b);
b1 = bs(1);
b2 = bs(2);
r1 = a * b;
r2 = a2 * b2 - (((r1 - a1 * b1) - a2 * b1) - a1 * b2);
r = [r1, r2];
end

function [r] = dd_div(a, b)
u = a(1) / b(1);
t = mul12(u, b(1));
L = (a(1) - t(1) - t(2) + a(2) - u * b(2) ) / b(1);
r1 = u + L;
r2 = u - r1 +L;
r = [r1, r2];
end

function z = dd_add(x, y)
A = x(1);
a = x(2);
B = y(1);
b = y(2);
R = A+B;
r = 0.0;
 if (abs(A) > abs(B))
    r = A-R+B+b+a;
 else
    r = B-R+A+a+b;
 end
    C = R + r;
    c = R - C + r;
    z = [C, c];
end
    
function z = dd_sub(x, y)
A = x(1);
a = x(2);
B = y(1);
b = y(2);
R = A-B;
r = 0.0;
if (abs(A) > abs(B))
    r = A-R-B-b+a;
else
    r = -B-R+A+a-b;
end
    C = R + r;
    c = R - C + r;
    z = [C, c];
end
    
function p = dd_mul(a, b)
% Double-double precision multiplication
a1 = a(1);
a2 = a(2);
b1 = b(1);
b2 = b(2);
s = mul12(a1, b1);
s1 = s(1);
s2 = s(2);
s2 = s2 + a1 * b2;
s2 = s2 + a2 * b1;
p = fast_two_sum(s1, s2);
end

function r = fast_two_sum(a, b)
% Helper function to sum two floating-point numbers with non-overlapping bits
r1 = a + b;
v = r1 - a;
r2 = b - v;
r = [r1, r2];
end

function r = dd_pow(a, n)
% Double-double exponentiation by squaring
% a is a double-double number, and n is an integer exponent
r = [1, 0]; % Initialize result as double-double 1
x = a;      % Initialize base
while n > 0
    if mod(n, 2) == 1
        r = dd_mul(r, x); % Multiply result by current base if bit is 1
    end
    x = dd_mul(x, x);    % Square the base
    n = floor(n / 2);    % Move to the next bit
end
end

For my problem (interpolation), error is even smaler than with 100 digit symbolic solution. Strange, but true. More than here used, helper functions, are there for convinience.

BR,

FT