[2] | 1 | FUNCTION stat_error, data, window, DATA2 = data2, NITER = niter, BAVARD = bavard, MEAN = mean |
---|
| 2 | |
---|
| 3 | ; ------------------------------------ |
---|
| 4 | ; Program to estimate statistical error of either one time series (std) |
---|
| 5 | ; or significance of std dev difference between two time series |
---|
| 6 | |
---|
| 7 | ; (c) 2005, E. Guilyardi - thanks to Pascal Terray for theory |
---|
| 8 | |
---|
| 9 | ; Errors are estimated using a moving block bootstrap (of length |
---|
| 10 | ; 'window') to account for serial correlations |
---|
| 11 | |
---|
| 12 | ; INPUTS |
---|
| 13 | ; data: 1d array |
---|
| 14 | ; window: window size (int) |
---|
| 15 | ; NITER = number of iteration (999 by default to have 1% precision) |
---|
| 16 | |
---|
| 17 | ; 0. Inits |
---|
| 18 | ; --------- |
---|
| 19 | |
---|
| 20 | niter1 = 999 |
---|
| 21 | IF keyword_set(niter) THEN niter1 = niter |
---|
| 22 | |
---|
| 23 | errstats = fltarr(niter1) |
---|
| 24 | |
---|
| 25 | ; define array: either [data] or [data, data2] |
---|
| 26 | |
---|
| 27 | |
---|
| 28 | IF keyword_set(data2) THEN BEGIN |
---|
| 29 | ; 2 arrays |
---|
| 30 | array = [data, data2] |
---|
| 31 | nblocks = floor(n_elements(data2)/window) |
---|
| 32 | IF n_elements(data2) NE nblocks*window THEN BEGIN |
---|
| 33 | print, ' ERROR: window must divide data2 size: ', n_elements(data2), window |
---|
| 34 | return, -1 |
---|
| 35 | ENDIF |
---|
| 36 | sizetest = n_elements(data2) |
---|
| 37 | ENDIF ELSE BEGIN |
---|
| 38 | ; 1 array |
---|
| 39 | array = data |
---|
| 40 | nblocks = floor(n_elements(data)/window) |
---|
| 41 | ; IF n_elements(data) NE nblocks*window THEN BEGIN |
---|
| 42 | ; print, ' ERROR: window must divide data size: ', n_elements(data), window |
---|
| 43 | ; return, -1 |
---|
| 44 | ; ENDIF |
---|
| 45 | sizetest = n_elements(data) |
---|
| 46 | ENDELSE |
---|
| 47 | |
---|
| 48 | size = n_elements(array) |
---|
| 49 | nblockst = floor(size/window) |
---|
| 50 | |
---|
| 51 | IF keyword_set(bavard) THEN BEGIN |
---|
| 52 | |
---|
| 53 | print, ' Error stats on array size/min/max: ', size, min(array), max(array) |
---|
| 54 | print, ' using moving window of ', window, ' with N iter =', niter1 |
---|
| 55 | print, ' total number of blocks', nblockst |
---|
| 56 | print, ' number of blocks and size of test array', nblocks, sizetest |
---|
| 57 | |
---|
| 58 | ENDIF |
---|
| 59 | |
---|
| 60 | it = 0 |
---|
| 61 | seed = 0 |
---|
| 62 | |
---|
| 63 | idx0 = findgen(window) |
---|
| 64 | |
---|
| 65 | |
---|
| 66 | WHILE it LE niter1-1 DO BEGIN |
---|
| 67 | ; |
---|
| 68 | ; 1. Randomly define blocks |
---|
| 69 | ; ------------------------- |
---|
| 70 | ib = 0 |
---|
| 71 | rand = fix(nblockst*randomu(seed, nblocks)) |
---|
| 72 | |
---|
| 73 | testarr = array[(rand(0)*window)+idx0] |
---|
| 74 | |
---|
| 75 | WHILE ib LE nblocks - 2 DO BEGIN |
---|
| 76 | |
---|
| 77 | testarr = [testarr, array((rand(ib)*window)+idx0)] |
---|
| 78 | ib = ib + 1 |
---|
| 79 | ENDWHILE |
---|
| 80 | ; |
---|
| 81 | ; 2. Compute stats |
---|
| 82 | ; ---------------- |
---|
| 83 | IF keyword_set(mean) THEN BEGIN |
---|
| 84 | errstats(it) = mean(testarr) |
---|
| 85 | ENDIF ELSE BEGIN |
---|
| 86 | ; errstats(it) = sqrt((moment(testarr))[1]) |
---|
| 87 | errstats(it) = sqrt( total(testarr^2)/n_elements(testarr)) |
---|
| 88 | ENDELSE |
---|
| 89 | ; method 1 keep mean i.e = sqrt(sum(testarr(i)**2)/N) |
---|
| 90 | ; method 2 errstats(it) = sqrt((moment(testarr))[1]) |
---|
| 91 | |
---|
| 92 | it = it + 1 |
---|
| 93 | ENDWHILE |
---|
| 94 | |
---|
| 95 | ; |
---|
| 96 | ; 3. Organise output |
---|
| 97 | ; ------------------ |
---|
| 98 | ; |
---|
| 99 | |
---|
| 100 | return, errstats |
---|
| 101 | |
---|
| 102 | END |
---|