octave_packages/statistics-1.1.3/jackknife.m

   1 ## Copyright (C) 2011 Alexander Klein <alexander.klein@math.uni-giessen.de>
   2 ##
   3 ## This program is free software; you can redistribute it and/or modify it under
   4 ## the terms of the GNU General Public License as published by the Free Software
   5 ## Foundation; either version 3 of the License, or (at your option) any later
   6 ## version.
   7 ##
   8 ## This program is distributed in the hope that it will be useful, but WITHOUT
   9 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  10 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
  11 ## details.
  12 ##
  13 ## You should have received a copy of the GNU General Public License along with
  14 ## this program; if not, see <http://www.gnu.org/licenses/>.
  15
  16 ## -*- texinfo -*-
  17 ## @deftypefn{Function File} {@var{jackstat} =} jackknife (@var{E}, @var{x}, @dots{})
  18 ## Compute jackknife estimates of a parameter taking one or more given samples as parameters.
  19 ## In particular, @var{E} is the estimator to be jackknifed as a function name, handle,
  20 ## or inline function, and @var{x} is the sample for which the estimate is to be taken.
  21 ## The @var{i}-th entry of @var{jackstat} will contain the value of the estimator
  22 ## on the sample @var{x} with its @var{i}-th row omitted.
  23 ##
  24 ## @example
  25 ## @group
  26 ## jackstat(@var{i}) = @var{E}(@var{x}(1 : @var{i} - 1, @var{i} + 1 : length(@var{x})))
  27 ## @end group
  28 ## @end example
  29 ##
  30 ## Depending on the number of samples to be used, the estimator must have the appropriate form:
  31 ## If only one sample is used, then the estimator need not be concerned with cell arrays,
  32 ## for example jackknifing the standard deviation of a sample can be performed with
  33 ## @code{@var{jackstat} = jackknife (@@std, rand (100, 1))}.
  34 ## If, however, more than one sample is to be used, the samples must all be of equal size,
  35 ## and the estimator must address them as elements of a cell-array,
  36 ## in which they are aggregated in their order of appearance:
  37 ##
  38 ## @example
  39 ## @group
  40 ## @var{jackstat} = jackknife(@@(x) std(x@{1@})/var(x@{2@}), rand (100, 1), randn (100, 1)
  41 ## @end group
  42 ## @end example
  43 ##
  44 ## If all goes well, a theoretical value @var{P} for the parameter is already known,
  45 ## @var{n} is the sample size,
  46 ## @code{@var{t} = @var{n} * @var{E}(@var{x}) - (@var{n} - 1) * mean(@var{jackstat})}, and
  47 ## @code{@var{v} = sumsq(@var{n} * @var{E}(@var{x}) - (@var{n} - 1) * @var{jackstat} - @var{t}) / (@var{n} * (@var{n} - 1))}, then
  48 ## @code{(@var{t}-@var{P})/sqrt(@var{v})} should follow a t-distribution with @var{n}-1 degrees of freedom.
  49 ##
  50 ## Jackknifing is a well known method to reduce bias; further details can be found in:
  51 ## @itemize @bullet
  52 ## @item Rupert G. Miller: The jackknife-a review; Biometrika (1974) 61(1): 1-15; doi:10.1093/biomet/61.1.1
  53 ## @item Rupert G. Miller: Jackknifing Variances; Ann. Math. Statist. Volume 39, Number 2 (1968), 567-582; doi:10.1214/aoms/1177698418
  54 ## @item M. H. Quenouille: Notes on Bias in Estimation; Biometrika Vol. 43, No. 3/4 (Dec., 1956), pp. 353-360; doi:10.1093/biomet/43.3-4.353
  55 ## @end itemize
  56 ## @end deftypefn
  57
  58 ## Author: Alexander Klein <alexander.klein@math.uni-giessen.de>
  59 ## Created: 2011-11-25
  60
  61 function jackstat = jackknife ( anEstimator, varargin )
  62
  63
  64         ## Convert function name to handle if necessary, or throw
  65         ## an error.
  66         if ( !strcmp ( typeinfo ( anEstimator ), "function handle" ) )
  67
  68                 if ( isascii ( anEstimator ) )
  69
  70                         anEstimator = str2func ( anEstimator );
  71
  72                 else
  73
  74                         error ( "Estimators must be passed as function names or handles!" );
  75                 end
  76         end
  77
  78
  79         ## Simple jackknifing can be done with a single vector argument, and
  80         ## first and foremost with a function that does not care about
  81         ## cell-arrays.
  82         if ( length ( varargin ) == 1 && isnumeric ( varargin { 1 } ) )
  83
  84                 aSample = varargin { 1 };
  85
  86                 g = length ( aSample );
  87
  88                 jackstat = zeros ( 1, g );
  89
  90                 for k = 1 : g
  91                         jackstat ( k ) = anEstimator ( aSample ( [ 1 : k - 1, k + 1 : g ] ) );
  92                 end
  93
  94         ## More complicated input requires more work, however.
  95         else
  96
  97                 g = cellfun ( @(x) length ( x ), varargin );
  98
  99                 if ( any ( g - g ( 1 ) ) )
 100
 101                         error ( "All passed data must be of equal length!" );
 102                 end
 103
 104                 g = g ( 1 );
 105
 106                 jackstat = zeros ( 1, g );
 107
 108                 for k = 1 : g
 109
 110                         jackstat ( k ) = anEstimator ( cellfun ( @(x) x( [ 1 : k - 1, k + 1 : g ] ), varargin, "UniformOutput", false ) );
 111                 end
 112
 113         end
 114 endfunction
 115
 116
 117 %!test
 118 %! ##Example from Quenouille, Table 1
 119 %! d=[0.18 4.00 1.04 0.85 2.14 1.01 3.01 2.33 1.57 2.19];
 120 %! jackstat = jackknife ( @(x) 1/mean(x), d );
 121 %! assert ( 10 / mean(d) - 9 * mean(jackstat), 0.5240, 1e-5 );
 122
 123 %!demo
 124 %! for k = 1:1000
 125 %!  x=rand(10,1);
 126 %!  s(k)=std(x);
 127 %!  jackstat=jackknife(@std,x);
 128 %!  j(k)=10*std(x) - 9*mean(jackstat);
 129 %! end
 130 %! figure();hist([s',j'], 0:sqrt(1/12)/10:2*sqrt(1/12))
 131
 132 %!demo
 133 %! for k = 1:1000
 134 %!  x=randn(1,50);
 135 %!  y=rand(1,50);
 136 %!  jackstat=jackknife(@(x) std(x{1})/std(x{2}),y,x);
 137 %!  j(k)=50*std(y)/std(x) - 49*mean(jackstat);
 138 %!  v(k)=sumsq((50*std(y)/std(x) - 49*jackstat) - j(k)) / (50 * 49);
 139 %! end
 140 %! t=(j-sqrt(1/12))./sqrt(v);
 141 %! figure();plot(sort(tcdf(t,49)),"-;Almost linear mapping indicates good fit with t-distribution.;")