ajs_jackknife.pro (Anthony Smith's IDL routines)

; docformat = 'rst'
;+
; This function returns the jackknife method-estimated standard
; deviation (and optionally the jackknife mean) of n resamplings
; of a quantity, where each of the n estimates of the quantity
; is generated by removing (1/n) of the original sample data points.
;-

;+
; This function returns the jackknife method-estimated standard
; deviation (and optionally the jackknife mean) of n resamplings
; of a quantity, where each of the n estimates of the quantity
; is generated by removing (1/n) of the original sample data points.
; :Returns:
;    Return value is the jackknife standard deviation or, if
;    jackknife_estimates is [m, n], returns jackknife covariance matrix [m, m]
; :Params:
;    jackknife_estimates : in, required
;       [n] or [m, n] array of n estimates of a quantity, or n
;       estimates of m quantities
; :Keywords:
;    original_estimate : in, optional
;       Supply this value to get bias correction (see below)
;    bias_correction : out, optional
;       Bias-corrected estimate is x_c = x + b where x is the original
;       estimate and b is the bias_correction
;    jackknife_mean : out, optional
;       Mean of the jackknife samples (NB different from
;       bias-corrected estimate)
; :History:
;    24 Jan 2008 Created, Anthony Smith
;
;    7 Apr 2008 Added covariance output for [m, n] input
;-
FUNCTION ajs_jackknife, jackknife_estimates, $
                        original_estimate=original_estimate, $
                        bias_correction=bias_correction, $
                        jackknife_mean=jackknife_mean
  compile_opt idl2
  debug = ajs_debug()
  IF debug GE 2 THEN message, 'Jackknife', /inf

  IF size(reform(jackknife_estimates), /n_dimensions) EQ 1 THEN BEGIN
      IF debug GE 2 THEN message, 'Single array of values', /inf
      n_jack = n_elements(jackknife_estimates)
      ;; Variance returns sum(deviations^2) / (n-1)
      jackknife_std = sqrt((n_jack - 1.) * variance(jackknife_estimates) $
                           * (n_jack - 1.) / n_jack)
      
      IF arg_present(jackknife_mean) THEN $
         jackknife_mean = mean(jackknife_estimates)
      
      IF arg_present(bias_correction) $
         AND n_elements(original_estimate) GT 0 THEN $
            bias_correction = (n_jack - 1.) $
                              * (original_estimate $
                                 - mean(jackknife_estimates))
      
      return, jackknife_std
  ENDIF ELSE IF size(reform(jackknife_estimates), /n_dimensions) EQ 2 $
  THEN BEGIN
      IF debug GE 2 THEN message, '2d array of values: covariance', /inf
      dims = size(jackknife_estimates, /dimensions)
      n_vars = dims[0]
      n_jack = dims[1]
      
      cov_mat = correlate(jackknife_estimates, /covariance)
      jackknife_cov_mat = (n_jack - 1.) * cov_mat $
                          * (n_jack - 1.) / n_jack

      IF arg_present(jackknife_mean) THEN $
         jackknife_mean = total(jackknife_estimates, 2) / float(n_jack)

      IF arg_present(bias_correction) $
         AND n_elements(original_estimate) GT 0 THEN $
            bias_correction = (n_jack - 1.) $
                              * (original_estimate $
                                 - (total(jackknife_estimates, 2) $
                                    / float(n_jack)))
      
      return, jackknife_cov_mat
  ENDIF

END

;+
; Test ajs_jackknife
;-
PRO ajs_jackknife_test
  compile_opt idl2

  ;; Random points from Gaussian
  n_points = 10000               ; Multiple of 100
  x = randomn(seed, n_points)
  indices = indgen(n_points)
  mean_x = mean(x)
  std_mean_x = sqrt(1. / n_points)
  print, 'Mean of', n_points, ' points drawn from N(0,1)'
  print, 'True mean:', mean_x
  print, 'Expected standard deviation on mean:', std_mean_x

  ;; Jackknife removing one point at a time
  jackknife_estimates = fltarr(n_points)
  FOR i = 0L, n_points - 1 DO BEGIN
      jackknife_estimates[i] = mean(x[where(indices NE i)])
  ENDFOR
  jackknife_std = ajs_jackknife(jackknife_estimates, $
                                original_estimate=mean_x, $
                                bias_correction=bias_correction)
  print, 'Jackknife one point at a time:'
  print, '  Bias-corrected mean:', mean_x + bias_correction
  print, '  Jackknife standard deviation:', jackknife_std

  ;; Jackknife removing 1/100 of the points at a time
  jackknife_estimates = fltarr(n_points / 100)
  FOR i = 0L, (n_points / 100) - 1 DO BEGIN
      jackknife_estimates[i] = mean(x[where(indices LT i * 100 $
                                            OR indices GE (i + 1) * 100)])
  ENDFOR
  jackknife_std = ajs_jackknife(jackknife_estimates, $
                                original_estimate=mean_x, $
                                bias_correction=bias_correction)
  print, '100 jackknife resamplings:'
  print, '  Bias-corrected mean:', mean_x + bias_correction
  print, '  Jackknife standard deviation:', jackknife_std

  ;; Jackknife removing 1/20 of the points at a time
  jackknife_estimates = fltarr(n_points / 20)
  FOR i = 0L, (n_points / 20) - 1 DO BEGIN
      jackknife_estimates[i] = mean(x[where(indices LT i * 20 $
                                            OR indices GE (i + 1) * 20)])
  ENDFOR
  jackknife_std = ajs_jackknife(jackknife_estimates, $
                                original_estimate=mean_x, $
                                bias_correction=bias_correction)
  print, '20 jackknife resamplings:'
  print, '  Bias-corrected mean:', mean_x + bias_correction
  print, '  Jackknife standard deviation:', jackknife_std

END