uncertainty in polyfit from measurements?
167 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
I have some measure point and ive fitted it w/ polyfit. The guesses alright but how can i find out the uncertainty in that coefficients? I used to use Origin for this but it crashes all te time so i decided to switch to matlab.
0 commentaires
Réponses (5)
John D'Errico
le 9 Fév 2011
There are several ways to do this, but I'm not certain what you are looking for. Are you looking for the uncertainty in a prediction of a fitted polynomial? Or are you just asking for the uncertainty in the fitted coefficients? Here is an example.
x = randn(100,1);
y = sin(x) + randn(size(x))/100;
P = polyfitn(x,y,3)
P =
ModelTerms: [4x1 double]
Coefficients: [-0.11492 0.0099578 0.94474 -0.0059965]
ParameterVar: [2.2444e-06 5.4231e-06 2.7675e-05 1.3305e-05]
ParameterStd: [0.0014981 0.0023288 0.0052607 0.0036476]
R2: 0.99788
RMSE: 0.029182
VarNames: {'X1'}
Polyfitn returns a standard deviation and variance for each parameter. The ratio of the coefficient to the standard deviation will give you a measure of the significance of each term.
P.Coefficients./P.ParameterStd
ans =
-76.712 4.276 179.58 -1.6439
My expectation is the linear and cubic terms will have been significant, since a sine wave has a series approximation with only odd order terms. I really want to compare these numbers to a Student t statistic, but I'm feeling lazy now. The error I've added was large enough that the quadratic term is actually potentially significant.
2 commentaires
Walter Roberson
le 8 Fév 2011
Asymptotically, the error for a polynomial of degree N fit from N+1 data points, product(x-A[i],i=1..N), will approach +/- sum(eps(A[i]),i=1..N+1) * x^(N-1) -- assuming, that is, perfect representation of x^N.
(If x itself is subject to error in representation then the question of the uncertainty in the coefficients does not arise in the same form.)
That's as x approaches infinity. Higher errors are possible within min(A[i]) to max(A[i]).
0 commentaires
Bruno Luong
le 9 Fév 2011
For Gaussian noise, here is a code to show you how to do without stat toolbox. There is a first part that uses brute force method and the second part that shows the method I propose. Both results allow to check that the estimation is correctly done:
x = linspace(0,1);
yext = 1+2*x+3.*x.^2;
sigmay = 0.1; % standard deviation of the noise data
ntest = 1000;
Ptab = zeros(ntest,3);
for n = 1:ntest
y = yext+sigmay*randn(size(yext));
P = polyfit(x,y,2);
Ptab(n,:) = P;
end
% Brute force estimate of uncertainty of P
std(Ptab,1)
% Compute the uncertainty
X = bsxfun(@power,x(:),2:-1:0);
% if not known sigmay can be estimated as std(polyval(P,x)-y)
stdP = sqrt(diag(inv(X'*X)))*sigmay
Bruno
0 commentaires
Arjen
le 9 Fév 2011
2 commentaires
Bruno Luong
le 9 Fév 2011
Yes you should have explained better. The latest attempt is no better. Bye then.
John D'Errico
le 9 Fév 2011
Explain what you want to see, else nobody will bother to try to guess your goal.
Voir également
Catégories
En savoir plus sur Logical dans Help Center et File Exchange
Produits
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!