1 ## Author: Paul Kienzle <pkienzle@gmail.com>
2 ## This program is granted to the public domain.
6 ## Test cases are taken from the NIST Statistical Reference Datasets
7 ## http://www.itl.nist.gov/div898/strd/
11 function do_test(n,x,y,p,dp,varargin)
12 [myp,s] = wpolyfit(x,y,n,varargin{:});
13 %if length(varargin)==0, [myp,s] = polyfit(x,y,n); else return; end
14 mydp = sqrt(sumsq(inv(s.R'))'/s.df)*s.normr;
15 if length(varargin)>0, mydp = [mydp;0]; end %origin
16 %[svdp,j,svddp] = svdfit(x,y,n);
17 disp('parameter certified value rel. error');
18 [myp(:), p, abs((myp(:)-p)./p)] %, svdp, abs((svdp-p)./p) ]
19 disp('p-error certified value rel. error');
20 [mydp(:), dp, abs((mydp(:) - dp)./dp)] %, svdp, abs((svddp - dp)./dp)]
21 input('Press <Enter> to proceed to the next test');
25 data = [0.0013852 0.2144023 0.0020470
26 0.0018469 0.2516856 0.0022868
27 0.0023087 0.3070443 0.0026362
28 0.0027704 0.3603186 0.0029670
29 0.0032322 0.4260864 0.0033705
30 0.0036939 0.4799956 0.0036983 ];
31 x=data(:,1); y=data(:,2); dy=data(:,3);
33 disp('computing parameter uncertainty from monte carlo simulation...');
36 for i=1:n, p(:,i)=(polyfit(x,y+randn(size(y)).*dy,1)).'; end
37 printf("%15s %15s\n", "Coefficient", "Error");
38 printf("%15g %15g\n", [mean(p'); std(p')]);
39 input('Press <Enter> to see some sample regression lines: ');
40 t = [x(1), x(length(x))];
41 [p,s] = wpolyfit(x,y,dy,1); dp=sqrt(sumsq(inv(s.R'))'/s.df)*s.normr;
43 for i=1:15, plot(t,polyval(p(:)+randn(size(dp)).*dp,t),'-g;;'); hold on; end
44 errorbar(x,y,dy,"~b;;");
45 [yf,dyf]=polyconf(p,x,s,0.05,'ci');
46 plot(x,yf-dyf,"-r;;",x,yf+dyf,'-r;95% confidence interval;')
48 input('Press <Enter> to continue with the tests: ');
51 ##Procedure: Linear Least Squares Regression
52 ##Reference: Filippelli, A., NIST.
53 ##Model: Polynomial Class
54 ## 11 Parameters (B0,B1,...,B10)
56 ## y = B0 + B1*x + B2*(x**2) + ... + B9*(x**9) + B10*(x**10) + e
60 data = [ 0.8116 -6.860120914
145 target = [ -1467.48961422980 298.084530995537
146 -2772.17959193342 559.779865474950
147 -2316.37108160893 466.477572127796
148 -1127.97394098372 227.204274477751
149 -354.478233703349 71.6478660875927
150 -75.1242017393757 15.2897178747400
151 -10.8753180355343 2.23691159816033
152 -1.06221498588947 0.221624321934227
153 -0.670191154593408E-01 0.142363763154724E-01
154 -0.246781078275479E-02 0.535617408889821E-03
155 -0.402962525080404E-04 0.896632837373868E-05];
157 disp("Filippelli, A., NIST.");
158 do_test(10, data(:,2),data(:,1),flipud(target(:,1)),flipud(target(:,2)));
161 ##Procedure: Linear Least Squares Regression
163 ##Reference: Pontius, P., NIST.
164 ## Load Cell Calibration.
166 ##Model: Quadratic Class
167 ## 3 Parameters (B0,B1,B2)
168 ## y = B0 + B1*x + B2*(x**2)
214 ## Certified Regression Statistics
216 ## Standard Deviation
217 ## Estimate of Estimate
219 0.673565789473684E-03 0.107938612033077E-03
220 0.732059160401003E-06 0.157817399981659E-09
221 -0.316081871345029E-14 0.486652849992036E-16 ];
224 disp("Pontius, P., NIST");
225 do_test(2, data(:,2),data(:,1),flipud(target(:,1)),flipud(target(:,2)));
230 #Procedure: Linear Least Squares Regression
231 #Reference: Eberhardt, K., NIST.
252 # Certified Regression Statistics
255 # Estimate of Estimate
258 2.07438016528926 0.165289256198347E-01 ];
262 disp("Eberhardt, K., NIST");
263 do_test(1, data(:,2),data(:,1),flipud(target(:,1)),flipud(target(:,2)),'origin');
267 #Reference: Wampler, R. H. (1970).
268 # A Report of the Accuracy of Some Widely-Used Least
269 # Squares Computer Programs.
270 # Journal of the American Statistical Association, 65, 549-565.
272 #Model: Polynomial Class
273 # 6 Parameters (B0,B1,...,B5)
275 # y = B0 + B1*x + B2*(x**2) + B3*(x**3)+ B4*(x**4) + B5*(x**5)
277 # Certified Regression Statistics
280 # Parameter Estimate of Estimate
282 1.00000000000000 0.000000000000000
283 1.00000000000000 0.000000000000000
284 1.00000000000000 0.000000000000000
285 1.00000000000000 0.000000000000000
286 1.00000000000000 0.000000000000000
287 1.00000000000000 0.000000000000000 ];
315 do_test(5, data(:,2),data(:,1),flipud(target(:,1)),flipud(target(:,2)));
318 ##Reference: Wampler, R. H. (1970).
319 ## A Report of the Accuracy of Some Widely-Used Least
320 ## Squares Computer Programs.
321 ## Journal of the American Statistical Association, 65, 549-565.
322 ##Model: Polynomial Class
323 ## 6 Parameters (B0,B1,...,B5)
325 ## y = B0 + B1*x + B2*(x**2) + B3*(x**3)+ B4*(x**4) + B5*(x**5)
327 ## Certified Regression Statistics
328 ## Standard Deviation
329 ## Parameter Estimate of Estimate
331 1.00000000000000 0.000000000000000
332 0.100000000000000 0.000000000000000
333 0.100000000000000E-01 0.000000000000000
334 0.100000000000000E-02 0.000000000000000
335 0.100000000000000E-03 0.000000000000000
336 0.100000000000000E-04 0.000000000000000 ];
365 do_test(5, data(:,2),data(:,1),flipud(target(:,1)),flipud(target(:,2)));
371 ##Reference: Wampler, R. H. (1970).
372 ## A Report of the Accuracy of Some Widely-Used Least
373 ## Squares Computer Programs.
374 ## Journal of the American Statistical Association, 65, 549-565.
376 ##Model: Polynomial Class
377 ## 6 Parameters (B0,B1,...,B5)
379 ## y = B0 + B1*x + B2*(x**2) + B3*(x**3)+ B4*(x**4) + B5*(x**5)
381 ## Certified Regression Statistics
383 ## Standard Deviation
384 ## Parameter Estimate of Estimate
386 1.00000000000000 2152.32624678170
387 1.00000000000000 2363.55173469681
388 1.00000000000000 779.343524331583
389 1.00000000000000 101.475507550350
390 1.00000000000000 5.64566512170752
391 1.00000000000000 0.112324854679312 ];
418 do_test(5, data(:,2),data(:,1),flipud(target(:,1)),flipud(target(:,2)));
421 ##Model: Polynomial Class
422 ## 6 Parameters (B0,B1,...,B5)
424 ## y = B0 + B1*x + B2*(x**2) + B3*(x**3)+ B4*(x**4) + B5*(x**5)
426 ## Certified Regression Statistics
428 ## Standard Deviation
429 ## Parameter Estimate of Estimate
431 1.00000000000000 215232.624678170
432 1.00000000000000 236355.173469681
433 1.00000000000000 77934.3524331583
434 1.00000000000000 10147.5507550350
435 1.00000000000000 564.566512170752
436 1.00000000000000 11.2324854679312 ];
464 do_test(5, data(:,2),data(:,1),flipud(target(:,1)),flipud(target(:,2)));
469 ##Model: Polynomial Class
470 ## 6 Parameters (B0,B1,...,B5)
472 ## y = B0 + B1*x + B2*(x**2) + B3*(x**3)+ B4*(x**4) + B5*(x**5)
474 ## Certified Regression Statistics
476 ## Standard Deviation
477 ## Parameter Estimate of Estimate
479 1.00000000000000 21523262.4678170
480 1.00000000000000 23635517.3469681
481 1.00000000000000 7793435.24331583
482 1.00000000000000 1014755.07550350
483 1.00000000000000 56456.6512170752
484 1.00000000000000 1123.24854679312 ];
511 do_test(5, data(:,2),data(:,1),flipud(target(:,1)),flipud(target(:,2)));