source: trunk/procs/stat_error.pro @ 27

FUNCTION stat_error, data, window, DATA2 = data2, NITER = niter, BAVARD = bavard, MEAN = mean

; ------------------------------------
; Program to estimate statistical error of either one time series (std)
; or significance of std dev difference between two time series

; (c) 2005, E. Guilyardi - thanks to Pascal Terray for theory

; Errors are estimated using a moving block bootstrap (of length
; 'window') to account for serial correlations

; INPUTS
; data: 1d array
; window: window size (int)
; NITER = number of iteration (999 by default to have 1% precision)

; 0. Inits
; ---------

   niter1 = 999
   IF keyword_set(niter) THEN niter1 = niter

   errstats = fltarr(niter1)

; define array: either [data] or [data, data2]


   IF keyword_set(data2) THEN BEGIN
      ; 2 arrays
      array = [data, data2]
      nblocks = floor(n_elements(data2)/window)
      IF n_elements(data2) NE nblocks*window THEN BEGIN 
         print,  ' ERROR: window must divide data2 size: ', n_elements(data2), window 
         return,  -1
      ENDIF 
      sizetest = n_elements(data2) 
   ENDIF ELSE BEGIN 
      ; 1 array
      array = data
      nblocks = floor(n_elements(data)/window)
;      IF n_elements(data) NE nblocks*window THEN BEGIN 
;         print,  ' ERROR: window must divide data size: ', n_elements(data), window 
;         return,  -1
;      ENDIF 
      sizetest = n_elements(data) 
   ENDELSE   

   size = n_elements(array) 
   nblockst = floor(size/window)

   IF keyword_set(bavard) THEN BEGIN  

      print, '  Error stats on array size/min/max: ', size, min(array), max(array)
      print, '        using moving window of ',  window, ' with N iter =',  niter1
      print, '        total number of blocks', nblockst 
      print, '        number of blocks and size of test array', nblocks, sizetest

   ENDIF 

   it = 0
   seed = 0

   idx0 = findgen(window)


   WHILE it LE niter1-1 DO BEGIN  
;
; 1. Randomly define blocks
; -------------------------
      ib = 0
      rand = fix(nblockst*randomu(seed, nblocks))

      testarr = array[(rand(0)*window)+idx0]

      WHILE ib LE nblocks - 2 DO BEGIN 

         testarr = [testarr, array((rand(ib)*window)+idx0)]
         ib = ib + 1
      ENDWHILE 
;
; 2. Compute stats
; ----------------
      IF keyword_set(mean) THEN BEGIN 
         errstats(it) = mean(testarr)
      ENDIF ELSE BEGIN  
;         errstats(it) = sqrt((moment(testarr))[1])
         errstats(it) = sqrt( total(testarr^2)/n_elements(testarr))
      ENDELSE 
; method 1 keep mean i.e = sqrt(sum(testarr(i)**2)/N)
; method 2      errstats(it) = sqrt((moment(testarr))[1])

      it = it + 1
   ENDWHILE   

;
; 3. Organise output
; ------------------
;

   return, errstats

END