Polyfit with scaling

[Ben Abbott]

Using the same approach as applied in wpolyfit, it is possible to  
improve the numerical stability/accuracy of polyfit.

However, it is not clear how the second output should be handled.

[P, S] = polyfit (X, Y, N)

S.R: The Cholesky factor of the Vandermonde matrix used to compute the  
polynomial coefficients.
S.X: The Vandermonde matrix used to compute the polynomial coefficients.
S.df: The degrees of freedom.
S.normr: The norm of the residuals.
S.yf: The values of the polynomial for each value of X.

Essentially, when there is no normalization of X the outputs are

         L = numel (X);
	V = (X * ones (1, N+1)) .^ (ones (L, 1) * (N : -1 : 0));
	P = (V \ Y).';
	S.yf = polyval (P, X);
	S.df = numel (X) - N;
	S.X  = V;
	s.R  = chol (V'*V);

When X is normalized

	Xn = (X - mean (X)) / std (X - mean (X));
	V = (Xn * ones (1, n+1)) .^ (ones (L, 1) * (n : -1 : 0));
	P = (V \ Y).';
         P = polyscale (P, std (X - mean (X)));
	P = polyshift (P, -mean (X));

I'm not sure what should be done with the second output when  
normalization is in place. I think the proper treatment depends upon  
what the results are used for.

Does anyone have experience/insight into how/where S.S, and S.R are  
used?

Ben