Changeset 41340 for trunk/Ohana/src/opihi/cmd.data/test/medimage.sh
- Timestamp:
- Apr 16, 2020, 1:54:47 PM (6 years ago)
- File:
-
- 1 edited
-
trunk/Ohana/src/opihi/cmd.data/test/medimage.sh (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
trunk/Ohana/src/opihi/cmd.data/test/medimage.sh
r41164 r41340 1 1 2 macro test.mean 3 4 $Nsample = 16 5 medimage delete -q t1 6 for i 0 $Nsample 7 mgaussdev t 50 50 0.0 1.0 8 medimage add t1 t 9 imhist -q t x y -range -10 10 -delta 0.1 10 11 $C0 = 0 12 $C1 = 1.5 13 $C2 = 400 14 $C3 = 0 15 vgauss -q x y con yf 16 # echo $C1 2 macro go 3 mcreate a 30 30 4 for i 0 40 5 set a$i = zero(a) + $i 17 6 end 18 7 19 medimage calc t1 T -mean 20 21 imhist -q T x y -range -10 10 -delta 0.1 22 lim -n 1 x y; clear; box; plot -x hist x y 23 $C0 = 0 24 $C1 = 1.5 25 $C2 = 400 26 $C3 = 0 27 vgauss -q x y con yf 28 echo "expect {1/sqrt($Nsample)} : $C1" 29 plot -c red -x line x yf 8 for i 0 40 9 medimage add t1 a$i 10 end 30 11 end 31 32 macro test.median33 34 $Nsample = 1635 medimage delete -q t136 for i 0 $Nsample37 mgaussdev t 50 50 0.0 1.038 medimage add t1 t39 imhist -q t x y -range -10 10 -delta 0.140 41 $C0 = 042 $C1 = 1.543 $C2 = 40044 $C3 = 045 vgauss -q x y con yf46 end47 48 # note that median of a gaussian distributed variable is not distributed with sigma' = sigma / sqrt(N)49 # (somewhat higher scatter)50 medimage calc t1 T51 52 imhist -q T x y -range -10 10 -delta 0.153 lim -n 1 x y; clear; box; plot -x hist x y54 $C0 = 055 $C1 = 1.556 $C2 = 40057 $C3 = 058 vgauss -q x y con yf59 echo "expect {1/sqrt($Nsample)} : $C1 (actually should be a bit higher)"60 plot -c red -x line x yf61 end62 63 macro test.wtmean64 65 $Nsample = 866 $sig1 = 1.067 $sig2 = 3.068 69 medimage delete -q t170 for i 0 $Nsample71 mgaussdev t 50 50 0.0 $sig172 set v = $sig1^2 + zero(t)73 medimage add t1 t -variance v74 imhist -q t x y -range -10 10 -delta 0.175 76 $C0 = 077 $C1 = 1.578 $C2 = 40079 $C3 = 080 vgauss -q x y con yf81 # echo $C182 end83 84 for i 0 $Nsample85 mgaussdev t 50 50 0.0 $sig286 set v = $sig2^2 + zero(t)87 medimage add t1 t -variance v88 imhist -q t x y -range -10 10 -delta 0.189 90 $C0 = 091 $C1 = 1.592 $C2 = 40093 $C3 = 094 vgauss -q x y con yf95 # echo $C196 end97 98 # note that median of a gaussian distributed variable is not distributed with sigma' = sigma / sqrt(N)99 # (somewhat higher scatter)100 medimage calc t1 T -wtmean101 102 imhist -q T x y -range -10 10 -delta 0.1103 lim -n 1 x y; clear; box; plot -x hist x y104 $C0 = 0105 $C1 = 1.5106 $C2 = 400107 $C3 = 0108 vgauss -q x y con yf109 $S1 = $Nsample / $sig1^2 + $Nsample / $sig2^2110 echo "expect {1/sqrt($S1)} : $C1"111 plot -c red -x line x yf112 end113 114 macro test.irls115 medimage delete -q t1116 $Nsample = 16117 $sig = 1.0118 for i 0 $Nsample119 mgaussdev t 50 50 0.0 $sig120 set v = $sig^2 + zero(t)121 122 set bad = (rnd(t) < 0.05) ? 10*rnd(t) + 5 : zero(t)123 set ts = t + bad124 125 medimage add t1 ts -variance v126 imhist -q t x y -range -10 10 -delta 0.1127 128 $C0 = 0129 $C1 = 1.5130 $C2 = 400131 $C3 = 0132 vgauss -q x y con yf133 # echo $C1134 end135 136 # get stats for straight mean:137 medimage calc t1 Tm -mean138 139 imhist -q Tm x y -range -10 10 -delta 0.1140 lim -n 1 x y; clear; box; plot -x hist x y141 $C0 = 0142 $C1 = 1.5143 $C2 = 400144 $C3 = 0145 vgauss -q x y con yf146 echo "sigma from straight stdev: $C1"147 # stats Tm148 149 plot -c red -x line x yf150 151 # get stats for irls152 medimage calc t1 Ti -irls153 154 imhist -q Ti x y -range -10 10 -delta 0.1155 lim -n 2 x y; clear; box; plot -x hist x y156 $C0 = 0157 $C1 = 1.5158 $C2 = 400159 $C3 = 0160 vgauss -q x y con yf161 echo "sigma from irls: $C1 (ideal is {$sig/sqrt($Nsample)})"162 # stats Ti163 164 plot -c red -x line x yf165 end166 167 168 ###################33169 170 171 macro test.mean.var172 173 $Nsample = 64174 $sig = 2.0175 medimage delete -q t1176 for i 0 $Nsample177 mgaussdev t 100 100 0.0 $sig178 medimage add t1 t179 imhist -q t x y -range -10 10 -delta 0.1180 181 $C0 = 0182 $C1 = 1.5183 $C2 = 400184 $C3 = 0185 vgauss -q x y con yf186 # echo $C1187 end188 189 medimage calc t1 T -mean -variance Tv190 191 imhist -q T x y -range -10 10 -delta 0.1192 lim -n 1 x y; clear; box; plot -x hist x y193 $C0 = 0194 $C1 = 1.5195 $C2 = 400196 $C3 = 0197 vgauss -q x y con yf198 plot -c red -x line x yf199 200 imhist Tv xv yv -range -1 4 -delta 0.1201 lim -n 2 xv yv; clear; box; plot xv yv -x hist202 203 stat Tv204 echo "$C1 vs {sqrt($MEDIAN)} : expect {$sig/sqrt($Nsample)}"205 end206 207 macro test.median.var208 209 $Nsample = 64210 $sig = 2.0211 medimage delete -q t1212 for i 0 $Nsample213 mgaussdev t 50 50 0.0 $sig214 medimage add t1 t215 imhist -q t x y -range -10 10 -delta 0.1216 217 $C0 = 0218 $C1 = 1.5219 $C2 = 400220 $C3 = 0221 vgauss -q x y con yf222 # echo $C1223 end224 225 # note that median of a gaussian distributed variable is not distributed with sigma' = sigma / sqrt(N)226 # (somewhat higher scatter)227 medimage calc t1 T -variance Tv228 229 imhist -q T x y -range -10 10 -delta 0.1230 lim -n 1 x y; clear; box; plot -x hist x y231 $C0 = 0232 $C1 = 1.5233 $C2 = 400234 $C3 = 0235 vgauss -q x y con yf236 plot -c red -x line x yf237 238 imhist Tv xv yv -range -1 4 -delta 0.1239 lim -n 2 xv yv; clear; box; plot xv yv -x hist240 241 stat Tv242 echo "$C1 vs {sqrt($MEDIAN)} : expect {$sig/sqrt($Nsample)}"243 end244 245 macro test.wtmean.var246 247 $Nsample = 32248 $sig1 = 1.0249 $sig2 = 1.0250 251 medimage delete -q t1252 for i 0 $Nsample253 mgaussdev t 50 50 0.0 $sig1254 set v = $sig1^2 + zero(t)255 medimage add t1 t -variance v256 imhist -q t x y -range -10 10 -delta 0.1257 258 $C0 = 0259 $C1 = 1.5260 $C2 = 400261 $C3 = 0262 vgauss -q x y con yf263 # echo $C1264 end265 266 for i 0 $Nsample267 mgaussdev t 50 50 0.0 $sig2268 set v = $sig2^2 + zero(t)269 medimage add t1 t -variance v270 imhist -q t x y -range -10 10 -delta 0.1271 272 $C0 = 0273 $C1 = 1.5274 $C2 = 400275 $C3 = 0276 vgauss -q x y con yf277 # echo $C1278 end279 280 # note that median of a gaussian distributed variable is not distributed with sigma' = sigma / sqrt(N)281 # (somewhat higher scatter)282 medimage calc t1 T -wtmean -variance Tv283 284 imhist -q T x y -range -10 10 -delta 0.1285 lim -n 1 x y; clear; box; plot -x hist x y286 $C0 = 0287 $C1 = 1.5288 $C2 = 400289 $C3 = 0290 vgauss -q x y con yf291 plot -c red -x line x yf292 293 stat -q Tv294 $S1 = $Nsample / $sig1^2 + $Nsample / $sig2^2295 echo $C1 vs {sqrt($MEDIAN)} : expect {1/sqrt($S1)}296 end297 298 macro test.irls.var299 300 $Nsample = 16301 $sig = 1.0302 303 medimage delete -q t1304 for i 0 $Nsample305 mgaussdev t 50 50 0.0 $sig306 set v = $sig^2 + zero(t)307 308 set bad = (rnd(t) < 0.05) ? 10*rnd(t) + 5 : zero(t)309 set ts = t + bad310 311 medimage add t1 ts -variance v312 imhist -q t x y -range -10 10 -delta 0.1313 314 $C0 = 0315 $C1 = 1.5316 $C2 = 400317 $C3 = 0318 vgauss -q x y con yf319 # echo $C1320 end321 322 # get stats for straight mean:323 medimage calc t1 Tm -mean -variance Tv324 325 imhist -q Tm x y -range -10 10 -delta 0.1326 lim -n 1 x y; clear; box; plot -x hist x y327 $C0 = 0328 $C1 = 1.5329 $C2 = 400330 $C3 = 0331 vgauss -q x y con yf332 echo "sigma from straight stdev: $C1"333 # stats Tm334 335 plot -c red -x line x yf336 337 # get stats for irls338 medimage calc t1 Ti -irls -variance Tv339 340 imhist -q Ti x y -range -10 10 -delta 0.1341 lim -n 2 x y; clear; box; plot -x hist x y342 $C0 = 0343 $C1 = 1.5344 $C2 = 400345 $C3 = 0346 vgauss -q x y con yf347 echo "sigma from irls: $C1 (ideal is {$sig/sqrt($Nsample)})"348 # stats Ti349 350 plot -c red -x line x yf351 352 set dTv = sqrt(Tv)353 imhist dTv xv yv -range -1 4 -delta 0.02; lim -n 3 xv yv; clear; box; plot xv yv -x hist354 355 stat -q Tv356 echo $C1 vs {sqrt($MEDIAN)} (ideal is {$sig/sqrt($Nsample)})"357 end358 359 macro test.irls.boot.var360 361 $Nsample = 64362 $sig = 1.0363 364 medimage delete -q t1365 for i 0 $Nsample366 mgaussdev t 200 200 0.0 $sig367 set v = $sig^2 + zero(t)368 369 set bad = (rnd(t) < 0.05) ? 10*rnd(t) + 5 : zero(t)370 set ts = t + bad371 372 mgaussdev noise 200 200 0.0 0.5373 set ts = ts + noise374 375 medimage add t1 ts -variance v376 imhist -q t x y -range -10 10 -delta 0.1377 378 $C0 = 0379 $C1 = 1.5380 $C2 = 400381 $C3 = 0382 vgauss -q x y con yf383 # echo $C1384 end385 386 # get stats for straight mean:387 medimage calc t1 Tm -mean -variance Tv388 389 imhist -q Tm x y -range -10 10 -delta 0.02390 lim -n 1 x y; clear; box; plot -x hist x y391 peak -q x y392 $C0 = $peakpos393 $C1 = 1.5*$sig / sqrt($Nsample)394 $C2 = $peakval395 $C3 = 0396 vgauss -q x y con yf397 echo "sigma from straight stdev: $C1"398 # stats Tm399 400 plot -c red -x line x yf401 402 # get stats for irls403 date404 medimage calc t1 Ti -irls -variance Tv -bootstrap-iter 100405 date406 407 imhist -q Ti x y -range -10 10 -delta 0.02408 lim -n 2 x y; clear; box; plot -x hist x y409 peak -q x y410 $C0 = $peakpos411 $C1 = 1.5*$sig / sqrt($Nsample)412 $C2 = $peakval413 $C3 = 0414 vgauss -q x y con yf415 echo "sigma from irls: $C1 (ideal is {$sig/sqrt($Nsample)})"416 # stats Ti417 418 plot -c red -x line x yf419 420 set dTv = sqrt(Tv)421 imhist dTv xv yv -range 0 {5*$sig/sqrt($Nsample)} -delta 0.02; lim -n 3 xv yv; clear; box; plot xv yv -x hist422 423 stat -q Tv424 echo $C1 vs {sqrt($MEDIAN)} (ideal is {$sig/sqrt($Nsample)})"425 end426 427 ##############################428 macro test.irls.boot.test429 430 $Nsample = 100431 $sig1 = 1.0432 433 medimage delete -q t1434 for i 0 $Nsample435 mgaussdev t 100 100 0.0 $sig1436 set v = $sig1^2 + zero(t)437 438 medimage add t1 t -variance v439 end440 441 # get stats for irls442 medimage calc t1 Ti -irls -variance Tv -bootstrap443 444 imhist -q Ti x y -range {-10*$sig1/sqrt($Nsample)} {10*$sig1/sqrt($Nsample)} -delta 0.01445 lim -n 2 x y; clear; box; plot -x hist x y446 peak -q x y447 $C0 = $peakpos448 $C1 = 1.5*$sig1/sqrt($Nsample)449 $C2 = $peakval450 $C3 = 0451 vgauss x y con yf452 echo "sigma from irls: $C1 (ideal is {$sig1/sqrt($Nsample)})"453 # stats Ti454 455 plot -c red -x line x yf456 457 set dTv = sqrt(Tv)458 imhist dTv xv yv -range 0 {5*$sig1/sqrt($Nsample)} -delta 0.02; lim -n 3 xv yv; clear; box; plot xv yv -x hist459 460 stat -q irls_npt461 $Npix = $MEAN462 463 stat -q Tv464 echo "sigma of irls average: $C1, sqrt(mean) of irls variance: {sqrt($MEAN)}, (ideal is {$sig1/sqrt($Npix)})"465 end466 467 macro test.irls.range.var468 medimage delete -q t1469 for i 0 8470 mgaussdev t 50 50 0.0 1.0471 set v = 1.0 + zero(t)472 473 set bad = (rnd(t) < 0.05) ? 10*rnd(t) + 5 : zero(t)474 set ts = t + bad475 476 medimage add t1 ts -variance v477 imhist -q t x y -range -10 10 -delta 0.1478 479 $C0 = 0480 $C1 = 1.5481 $C2 = 400482 $C3 = 0483 vgauss -q x y con yf484 echo $C1485 end486 487 for i 0 8488 mgaussdev t 50 50 0.0 3.0489 set v = 3.0 + zero(t)490 491 set bad = (rnd(t) < 0.05) ? 10*rnd(t) + 5 : zero(t)492 set ts = t + bad493 494 medimage add t1 ts -variance v495 imhist -q t x y -range -10 10 -delta 0.1496 497 $C0 = 0498 $C1 = 1.5499 $C2 = 400500 $C3 = 0501 vgauss -q x y con yf502 echo $C1503 end504 505 # get stats for straight mean:506 medimage calc t1 Tm -mean -variance Tv507 508 imhist -q Tm x y -range -10 10 -delta 0.1509 lim -n 1 x y; clear; box; plot -x hist x y510 $C0 = 0511 $C1 = 1.5512 $C2 = 400513 $C3 = 0514 vgauss -q x y con yf515 echo $C1516 stats Tm517 518 plot -c red -x line x yf519 520 # get stats for irls521 medimage calc t1 Ti -irls -variance Tv522 523 imhist -q Ti x y -range -10 10 -delta 0.1524 lim -n 2 x y; clear; box; plot -x hist x y525 $C0 = 0526 $C1 = 1.5527 $C2 = 400528 $C3 = 0529 vgauss -q x y con yf530 echo $C1531 stats Ti532 533 plot -c red -x line x yf534 535 stat -q Tv536 echo $C1 vs {sqrt($MEDIAN)}537 538 set dTv = sqrt(Tv)539 imhist dTv xv yv -range -1 4 -delta 0.02; lim -n 3 xv yv; clear; box; plot xv yv -x hist540 end
Note:
See TracChangeset
for help on using the changeset viewer.
