/usr/share/octave/packages/optim-1.3.0/test_wpolyfit.m is in octave-optim 1.3.0-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 | ## Author: Paul Kienzle <pkienzle@gmail.com>
## This program is granted to the public domain.
## Tests for wpolyfit.
##
## Test cases are taken from the NIST Statistical Reference Datasets
## http://www.itl.nist.gov/div898/strd/
1;
function do_test(n,x,y,p,dp,varargin)
[myp,s] = wpolyfit(x,y,n,varargin{:});
%if length(varargin)==0, [myp,s] = polyfit(x,y,n); else return; end
mydp = sqrt(sumsq(inv(s.R'))'/s.df)*s.normr;
if length(varargin)>0, mydp = [mydp;0]; end %origin
%[svdp,j,svddp] = svdfit(x,y,n);
disp('parameter certified value rel. error');
[myp(:), p, abs((myp(:)-p)./p)] %, svdp, abs((svdp-p)./p) ]
disp('p-error certified value rel. error');
[mydp(:), dp, abs((mydp(:) - dp)./dp)] %, svdp, abs((svddp - dp)./dp)]
input('Press <Enter> to proceed to the next test');
endfunction
## x y dy
data = [0.0013852 0.2144023 0.0020470
0.0018469 0.2516856 0.0022868
0.0023087 0.3070443 0.0026362
0.0027704 0.3603186 0.0029670
0.0032322 0.4260864 0.0033705
0.0036939 0.4799956 0.0036983 ];
x=data(:,1); y=data(:,2); dy=data(:,3);
wpolyfit(x,y,dy,1);
disp('computing parameter uncertainty from monte carlo simulation...');
fflush(stdout);
n=100; p=zeros(2,n);
for i=1:n, p(:,i)=(polyfit(x,y+randn(size(y)).*dy,1)).'; end
printf("%15s %15s\n", "Coefficient", "Error");
printf("%15g %15g\n", [mean(p'); std(p')]);
input('Press <Enter> to see some sample regression lines: ');
t = [x(1), x(length(x))];
[p,s] = wpolyfit(x,y,dy,1); dp=sqrt(sumsq(inv(s.R'))'/s.df)*s.normr;
hold off;
for i=1:15, plot(t,polyval(p(:)+randn(size(dp)).*dp,t),'-g;;'); hold on; end
errorbar(x,y,dy,"~b;;");
[yf,dyf]=polyconf(p,x,s,0.05,'ci');
plot(x,yf-dyf,"-r;;",x,yf+dyf,'-r;95% confidence interval;')
hold off;
input('Press <Enter> to continue with the tests: ');
##Procedure: Linear Least Squares Regression
##Reference: Filippelli, A., NIST.
##Model: Polynomial Class
## 11 Parameters (B0,B1,...,B10)
##
## y = B0 + B1*x + B2*(x**2) + ... + B9*(x**9) + B10*(x**10) + e
##Data:
## y x
data = [ 0.8116 -6.860120914
0.9072 -4.324130045
0.9052 -4.358625055
0.9039 -4.358426747
0.8053 -6.955852379
0.8377 -6.661145254
0.8667 -6.355462942
0.8809 -6.118102026
0.7975 -7.115148017
0.8162 -6.815308569
0.8515 -6.519993057
0.8766 -6.204119983
0.8885 -5.853871964
0.8859 -6.109523091
0.8959 -5.79832982
0.8913 -5.482672118
0.8959 -5.171791386
0.8971 -4.851705903
0.9021 -4.517126416
0.909 -4.143573228
0.9139 -3.709075441
0.9199 -3.499489089
0.8692 -6.300769497
0.8872 -5.953504836
0.89 -5.642065153
0.891 -5.031376979
0.8977 -4.680685696
0.9035 -4.329846955
0.9078 -3.928486195
0.7675 -8.56735134
0.7705 -8.363211311
0.7713 -8.107682739
0.7736 -7.823908741
0.7775 -7.522878745
0.7841 -7.218819279
0.7971 -6.920818754
0.8329 -6.628932138
0.8641 -6.323946875
0.8804 -5.991399828
0.7668 -8.781464495
0.7633 -8.663140179
0.7678 -8.473531488
0.7697 -8.247337057
0.77 -7.971428747
0.7749 -7.676129393
0.7796 -7.352812702
0.7897 -7.072065318
0.8131 -6.774174009
0.8498 -6.478861916
0.8741 -6.159517513
0.8061 -6.835647144
0.846 -6.53165267
0.8751 -6.224098421
0.8856 -5.910094889
0.8919 -5.598599459
0.8934 -5.290645224
0.894 -4.974284616
0.8957 -4.64454848
0.9047 -4.290560426
0.9129 -3.885055584
0.9209 -3.408378962
0.9219 -3.13200249
0.7739 -8.726767166
0.7681 -8.66695597
0.7665 -8.511026475
0.7703 -8.165388579
0.7702 -7.886056648
0.7761 -7.588043762
0.7809 -7.283412422
0.7961 -6.995678626
0.8253 -6.691862621
0.8602 -6.392544977
0.8809 -6.067374056
0.8301 -6.684029655
0.8664 -6.378719832
0.8834 -6.065855188
0.8898 -5.752272167
0.8964 -5.132414673
0.8963 -4.811352704
0.9074 -4.098269308
0.9119 -3.66174277
0.9228 -3.2644011];
##Certified values:
## p dP
target = [ -1467.48961422980 298.084530995537
-2772.17959193342 559.779865474950
-2316.37108160893 466.477572127796
-1127.97394098372 227.204274477751
-354.478233703349 71.6478660875927
-75.1242017393757 15.2897178747400
-10.8753180355343 2.23691159816033
-1.06221498588947 0.221624321934227
-0.670191154593408E-01 0.142363763154724E-01
-0.246781078275479E-02 0.535617408889821E-03
-0.402962525080404E-04 0.896632837373868E-05];
if 1
disp("Filippelli, A., NIST.");
do_test(10, data(:,2),data(:,1),flipud(target(:,1)),flipud(target(:,2)));
endif
##Procedure: Linear Least Squares Regression
##
##Reference: Pontius, P., NIST.
## Load Cell Calibration.
##
##Model: Quadratic Class
## 3 Parameters (B0,B1,B2)
## y = B0 + B1*x + B2*(x**2)
##Data: y x
data = [ ...
.11019 150000
.21956 300000
.32949 450000
.43899 600000
.54803 750000
.65694 900000
.76562 1050000
.87487 1200000
.98292 1350000
1.09146 1500000
1.20001 1650000
1.30822 1800000
1.41599 1950000
1.52399 2100000
1.63194 2250000
1.73947 2400000
1.84646 2550000
1.95392 2700000
2.06128 2850000
2.16844 3000000
.11052 150000
.22018 300000
.32939 450000
.43886 600000
.54798 750000
.65739 900000
.76596 1050000
.87474 1200000
.98300 1350000
1.09150 1500000
1.20004 1650000
1.30818 1800000
1.41613 1950000
1.52408 2100000
1.63159 2250000
1.73965 2400000
1.84696 2550000
1.95445 2700000
2.06177 2850000
2.16829 3000000 ];
## Certified Regression Statistics
##
## Standard Deviation
## Estimate of Estimate
target = [ ...
0.673565789473684E-03 0.107938612033077E-03
0.732059160401003E-06 0.157817399981659E-09
-0.316081871345029E-14 0.486652849992036E-16 ];
if 1
disp("Pontius, P., NIST");
do_test(2, data(:,2),data(:,1),flipud(target(:,1)),flipud(target(:,2)));
endif
#Procedure: Linear Least Squares Regression
#Reference: Eberhardt, K., NIST.
#Model: Linear Class
# 1 Parameter (B1)
#
# y = B1*x + e
#Data: y x
data =[...
130 60
131 61
132 62
133 63
134 64
135 65
136 66
137 67
138 68
139 69
140 70 ];
# Certified Regression Statistics
#
# Standard Deviation
# Estimate of Estimate
target = [ ...
0 0
2.07438016528926 0.165289256198347E-01 ];
if 1
disp("Eberhardt, K., NIST");
do_test(1, data(:,2),data(:,1),flipud(target(:,1)),flipud(target(:,2)),'origin');
endif
#Reference: Wampler, R. H. (1970).
# A Report of the Accuracy of Some Widely-Used Least
# Squares Computer Programs.
# Journal of the American Statistical Association, 65, 549-565.
#
#Model: Polynomial Class
# 6 Parameters (B0,B1,...,B5)
#
# y = B0 + B1*x + B2*(x**2) + B3*(x**3)+ B4*(x**4) + B5*(x**5)
#
# Certified Regression Statistics
#
# Standard Deviation
# Parameter Estimate of Estimate
target = [...
1.00000000000000 0.000000000000000
1.00000000000000 0.000000000000000
1.00000000000000 0.000000000000000
1.00000000000000 0.000000000000000
1.00000000000000 0.000000000000000
1.00000000000000 0.000000000000000 ];
#Data: y x
data = [...
1 0
6 1
63 2
364 3
1365 4
3906 5
9331 6
19608 7
37449 8
66430 9
111111 10
177156 11
271453 12
402234 13
579195 14
813616 15
1118481 16
1508598 17
2000719 18
2613660 19
3368421 20 ];
if 1
disp("Wampler1");
do_test(5, data(:,2),data(:,1),flipud(target(:,1)),flipud(target(:,2)));
endif
##Reference: Wampler, R. H. (1970).
## A Report of the Accuracy of Some Widely-Used Least
## Squares Computer Programs.
## Journal of the American Statistical Association, 65, 549-565.
##Model: Polynomial Class
## 6 Parameters (B0,B1,...,B5)
##
## y = B0 + B1*x + B2*(x**2) + B3*(x**3)+ B4*(x**4) + B5*(x**5)
##
## Certified Regression Statistics
## Standard Deviation
## Parameter Estimate of Estimate
target = [ ...
1.00000000000000 0.000000000000000
0.100000000000000 0.000000000000000
0.100000000000000E-01 0.000000000000000
0.100000000000000E-02 0.000000000000000
0.100000000000000E-03 0.000000000000000
0.100000000000000E-04 0.000000000000000 ];
#Data: y x
data = [ ...
1.00000 0
1.11111 1
1.24992 2
1.42753 3
1.65984 4
1.96875 5
2.38336 6
2.94117 7
3.68928 8
4.68559 9
6.00000 10
7.71561 11
9.92992 12
12.75603 13
16.32384 14
20.78125 15
26.29536 16
33.05367 17
41.26528 18
51.16209 19
63.00000 20 ];
if 1
disp("Wampler2");
do_test(5, data(:,2),data(:,1),flipud(target(:,1)),flipud(target(:,2)));
endif
##Reference: Wampler, R. H. (1970).
## A Report of the Accuracy of Some Widely-Used Least
## Squares Computer Programs.
## Journal of the American Statistical Association, 65, 549-565.
##
##Model: Polynomial Class
## 6 Parameters (B0,B1,...,B5)
##
## y = B0 + B1*x + B2*(x**2) + B3*(x**3)+ B4*(x**4) + B5*(x**5)
##
## Certified Regression Statistics
##
## Standard Deviation
## Parameter Estimate of Estimate
target = [...
1.00000000000000 2152.32624678170
1.00000000000000 2363.55173469681
1.00000000000000 779.343524331583
1.00000000000000 101.475507550350
1.00000000000000 5.64566512170752
1.00000000000000 0.112324854679312 ];
#Data: y x
data = [ ...
760. 0
-2042. 1
2111. 2
-1684. 3
3888. 4
1858. 5
11379. 6
17560. 7
39287. 8
64382. 9
113159. 10
175108. 11
273291. 12
400186. 13
581243. 14
811568. 15
1121004. 16
1506550. 17
2002767. 18
2611612. 19
3369180. 20 ];
if 1
disp("Wampler3");
do_test(5, data(:,2),data(:,1),flipud(target(:,1)),flipud(target(:,2)));
endif
##Model: Polynomial Class
## 6 Parameters (B0,B1,...,B5)
##
## y = B0 + B1*x + B2*(x**2) + B3*(x**3)+ B4*(x**4) + B5*(x**5)
##
## Certified Regression Statistics
##
## Standard Deviation
## Parameter Estimate of Estimate
target = [...
1.00000000000000 215232.624678170
1.00000000000000 236355.173469681
1.00000000000000 77934.3524331583
1.00000000000000 10147.5507550350
1.00000000000000 564.566512170752
1.00000000000000 11.2324854679312 ];
#Data: y x
data = [...
75901 0
-204794 1
204863 2
-204436 3
253665 4
-200894 5
214131 6
-185192 7
221249 8
-138370 9
315911 10
-27644 11
455253 12
197434 13
783995 14
608816 15
1370781 16
1303798 17
2205519 18
2408860 19
3444321 20 ];
if 1
disp("Wampler4");
do_test(5, data(:,2),data(:,1),flipud(target(:,1)),flipud(target(:,2)));
endif
##Model: Polynomial Class
## 6 Parameters (B0,B1,...,B5)
##
## y = B0 + B1*x + B2*(x**2) + B3*(x**3)+ B4*(x**4) + B5*(x**5)
##
## Certified Regression Statistics
##
## Standard Deviation
## Parameter Estimate of Estimate
target = [...
1.00000000000000 21523262.4678170
1.00000000000000 23635517.3469681
1.00000000000000 7793435.24331583
1.00000000000000 1014755.07550350
1.00000000000000 56456.6512170752
1.00000000000000 1123.24854679312 ];
##Data: y x
data = [ ...
7590001 0
-20479994 1
20480063 2
-20479636 3
25231365 4
-20476094 5
20489331 6
-20460392 7
18417449 8
-20413570 9
20591111 10
-20302844 11
18651453 12
-20077766 13
21059195 14
-19666384 15
26348481 16
-18971402 17
22480719 18
-17866340 19
10958421 20 ];
if 1
disp("Wampler5");
do_test(5, data(:,2),data(:,1),flipud(target(:,1)),flipud(target(:,2)));
endif
|