(*****************************************************************************
Conditions of Use

Of the Mathematica package Bootstrap.m

by Sergej V. Aksenov (http://aksenov.freeshell.org/)


This software (file Bootstrap.m) is copyright by Sergej V. Aksenov and 
Michael A. Savageau.

This software and its constituent parts come with no warranty, whether expressed 
or implied, that it is free of errors or suitable for any specific purpose.  It 
must not be used to solve any problem, whose incorrect solution could result in 
injury to a person, institution, or property.  It is user's own risk to use this 
software or its parts and the authors disclaim all liability for such use.

This software is distributed ``as is''. In particular, no maintenance, support, 
troubleshooting, or subsequent upgrade is implied.

The use of this software must be acknowledged in any publication that contains 
results obtained with it; please contact the author for the citation.

The free use of this software and its parts is restricted for research 
purposes; commercial use require permission and licensing from Sergej V. 
Aksenov and Michael A. Savageau.
******************************************************************************)



(*:Mathematica Version: 4.1 *)

(*:Package Version: 1.12 *)

(*:Name: Bootstrap` *)

(*:Context: Bootstrap`Bootstrap` *)

(*:Title: Confidence intervals by bootstrap *)

(*:Author: Sergej V. Aksenov (http://aksenov.freeshell.org/) *)

(*:History:
   Created: August, 2001.
   v1.01 Last modified: October, 2001.
   v1.1 Last modified 22 April, 2002; rewritten BootstrapResamples and BootstrapCI.
   v1.11 Last Modified 23 April, 2002; changed default value and corrected a bug for ConfidenceInterval; now passing options to RootSearch.
   v1.12 Last modified 29 April, 2002; If JackknifeReplications failes, it now returns a list {$Failed, <ind>}, where <ind> is the index of the point in the  data list that caused failure; fixed a bug in BootstrapReplications.
*)

(*:Copyright: Copyright 2001-2, Sergej V. Aksenov and Michael A. Savageau. *)

(*:Reference: 
Usage messages, documentation file Bootstrap.nb, article "Mathematica and C programs for minimum distance estimation of the S distribution and for calculation of goodness-of-fit by bootstrap", by Sergej V. Aksenov and Michael A. Savageau, in preparation (2002). 
Algorithm for calculation of BCa percentile intervals is in B. Efron and R. J. Tibshirani, "An introduction to the Bootstrap", Chapman and Hall, Inc., New York, London, 1993.
Algorithm for calculation of extreme percentile intervals is in S. M. S. Lee, "Nonparametric confidence intervals based on extreme bootstrap percentiles", Statistica Sinica, vol.10, pp.475-496, 2000.*)

(*:Summary:
This package provides functions for resampling a statistic of univariate random variable in parametric and nonparamteric modes, computation of confidence intervals based on bootstrap by extreme and BCa percentile methods. *)

(*:Keywords: 
Bootstrap, Jackknife, Confidence interval, bias-corrected and accelerated (BCa) percentile, extreme percentile. *)

(*:Limitations: *)

(*:Note: *)



BeginPackage["Bootstrap`Bootstrap`", "Statistics`DescriptiveStatistics`", "Statistics`Common`DistributionsCommon`", "Statistics`ContinuousDistributions`",  "Statistics`NormalDistribution`", "Statistics`ConfidenceIntervals`", "Utilities`FilterOptions`", "Enhancements`RootSearch`"]



(* Function names that will be exported. *)

BootstrapResamples::usage = 
"BootstrapResamples[obs, jack] returns a list {B1, B2} of estimated numbers of bootstrap resamples needed for construction of an extreme percentile confidence interval, given the observed value of the statistic (obs) and a list of jackknife replications for that statistic (jack)."

BootstrapCI::usage = 
"BootstrapCI[obs, boot, jack] returns a list {lower, upper} representing a bootsrap confidence interval, given the observed value of the statistic (obs), a list of its bootstrap replications (boot), and a list of its jackknife replications (jack)."

VarianceNonParametricBCaCIEndpoints::usage =
"VarianceNonParametricBCaCIEndpoints[obs, boot, jack, bootinds] returns a list {lower, upper} representing estimates of variance of BCa percentile interval endpoints using the jackknife-after-bootstrap method, given observed value of the statistic (obs), a list of nonparametric bootstrap replications of the statistic (boot), a list of jackknife replications of the statistic (jack), and a list of indices of sample data points that do not occur in each nonparametric bootstrap resample (bootinds)."

JackknifeReplications::usage = 
"JackknifeReplications[fun, dat] returns a list of jackknife replications of a statistic, given a function that computes the statistic (fun) and the sample data (dat)."

BootstrapReplications::usage = 
"BootstrapReplications[fun, B, arg, S] returns a nested list {boot, bootinds} containing a list of B bootstrap replications of a statistic (boot) and a list of B lists each corresponding to a bootstrap resample and containing indices of data points in the sample data list that are not present in that resample (bootinds).  Resamples are generated from sample data (nonparametric mode, arg is the list of data) or from a given distribution (parametric mode, arg is a completely specified distribution).  In case of parametric resamples, an optional argument S gives the size of a sample.  Statistic is computed with function fun."



(* Option names that will be exported. *)

BootstrapCIMethod::usage = 
"Option BootstrapCIMethod specifies bootstrap method for constructing confidence intervals: ExtremePercentiles uses extreme percentile method, BCaPercentiles uses BCa percentile method."

BootstrapSamplingMethod::usage = 
"Option BootstrapSamplingMethod specifies bootstrap resampling mode: NonParametric uses nonparametric resampling from observed sample with replacement, Parametric uses parametric resampling from a completely specified distribution."

ConfidenceInterval::usage =
"Option ConfidenceInterval specifies a type of an extreme percentile interval to be constructed: OneSidedUpper, OneSidedLower for intervals with the one-sided coverage and TwoSidedAsymmetrical and TwoSidedEquiTailed for intervals with the two-sided coverage."



(* Names for values of options that will be exported. *)

ExtremePercentiles::usage = 
"ExtremePercentiles is value for option BootstrapCIMethod that selects construction of extreme percentile confidence intervals."

BCaPercentiles::usage =
"BCaPercentiles is value for option BootstrapCIMethod that selects construction of BCa percentile confidence intervals."

NonParametric::usage = 
"NonParametric is value for option BootstrapSamplingMethod that selects nonparametric resampling from the observed sample with replacement."

Parametric::usage = 
"Parametric is value for option BootstrapSamplingMethod that selects parametric resampling from the completely specified distribution."

OneSidedUpper::usage =
"OneSidedUpper is value for option ConfidenceInterval that specifies one-sided upper confidence interval."

OneSidedLower::usage =
"OneSidedLower is value for option ConfidenceInterval that specifies one-sided lower confidence interval."

TwoSidedAsymmetrical::usage =
"TwoSidedAsymmetrical is value for option ConfidenceInterval that specifies two-sided asymmetrical (coverage in the tails is not equal) confidence interval."

TwoSidedEquiTailed::usage =
"TwoSidedEquiTailed is value for option ConfidenceInterval that specifies two-sided equi-tailed confidence interval."



(* Set default options. *)

Options[BootstrapResamples] = {ConfidenceLevel -> .95, ConfidenceInterval -> TwoSidedAsymmetrical}
Options[BootstrapCI] = {ConfidenceLevel -> .95, BootstrapCIMethod -> ExtremePercentiles, ConfidenceInterval -> TwoSidedAsymmetrical} 
Options[VarianceNonParametricBCaCIEndpoints] = {ConfidenceLevel -> .95}
Options[BootstrapReplications] = {BootstrapSamplingMethod -> NonParametric}



Begin["`Private`"]



(* Macros for frequently used functions (not exported). *)

Ncdf[x_]:=CDF[NormalDistribution[0,1],x]; 
Qn[x_]:=Quantile[NormalDistribution[0,1],x]; 



(* List of jackknife replications. *)

JackknifeReplications[statfun_, data_?VectorQ, opts___?OptionQ] := Module[
   {wasOn, jackresam, jackstat, try, ipos}, 
   wasOn = Head[General::stop] =!= $Off; If[wasOn, Off[General::stop]];
   jackresam=Map[Drop[#, 1] &, NestList[RotateLeft, data, Length[data] - 1]];
   jackstat=Catch[Map[(try=statfun[#,FilterOptions[statfun,opts]]; If[try===$Failed, ipos = Position[data, Complement[data, #][[1]]][[1,1]]; Throw[{$Failed, ipos}], try])&, jackresam]];
   If[wasOn, On[General::stop]]; 
   Return[jackstat];   
]



(* List of parametric or nonparametric bootstrap replications. *)

BootstrapReplications[statfun_, B_Integer?Positive, arg_, S_Integer:1, opts___?OptionQ] := Module[
   {meth, nonpar, wasOn, i, indnotthere, bootresam, bootstat, calc, NonParBootInds},
   meth = BootstrapSamplingMethod /. {opts} /. Options[BootstrapReplications];
   Which[meth===NonParametric, nonpar=True, meth===Parametric, nonpar=False, True, Message[BootstrapReplications::meth]; nonpar=True];
   If[nonpar && !VectorQ[arg], Message[BootstrapReplications::dat]; Return[$Failed]];
   wasOn = Head[General::stop] =!= $Off; If[wasOn, Off[General::stop]]; 
   bootstat = {}; i = 0; NonParBootInds = {};
   While[i < B,  
   If[nonpar, 
   bootresam = arg[[ Table[Random[Integer, {1,Length[arg]}],{Length[arg]}] ]],
   bootresam = RandomArray[arg, S]];  
   calc = statfun[bootresam,FilterOptions[statfun,opts]];  
   If[calc =!= $Failed, 
   AppendTo[bootstat,calc];
   If[nonpar, indnotthere = Complement[arg, bootresam]; 
   AppendTo[NonParBootInds, Flatten[MapThread[Position, {Table[arg, {Dimensions[indnotthere][[1]]}], indnotthere}]]]];
   i++]    ];
   If[wasOn, On[General::stop]]; If[nonpar, Protect[NonParBootInds]]; 
   Return[{bootstat,NonParBootInds}];
 ] /; Positive[S]

BootstrapReplications::meth =
"Warning: BootstrapSamplingMethod is not specified correctly, default value NonParametric is used."
BootstrapReplications::dat =
"Argument 3 is expected to be a list of sample data."



(* Number of bootstrap resamples for extreme percentiles method. *)

BootstrapResamples[obsval_, jkdat_, opts___?OptionQ] := 
  Module[{cov, sides, equ, eql, eq2a, b, B, B1, B2, term}, 
    cov = ConfidenceLevel /. {opts} /. Options[BootstrapResamples];
    If[cov ² 0 || 1 ² cov, Message[BootstrapResamples::clev]; 
      Return[$Failed]];
    sides = ConfidenceInterval /. {opts} /. Options[BootstrapResamples];
    If[Intersection[{sides}, {OneSidedUpper, OneSidedLower, 
            TwoSidedAsymmetrical, TwoSidedEquiTailed}] === {}, 
      Message[BootstrapResamples::cint]; sides = TwoSidedAsymmetrical];
    term = -Apply[Plus, (jkdat - obsval)^3]/
          6./(Apply[Plus, (jkdat - obsval)^2])^1.5;
    bsub = {B -> b/PDF[NormalDistribution[0, 1], b - 1/b]};
    equ = 1/(B + 1) + term b^3/B == 1 - cov /. bsub;
    eql = 1/(B + 1) - term b^3/B == 1 - cov /. bsub;
    eq2a = 2/(B + 1) + term^2 b^6/B == 1 - cov /. bsub;
    Off[Power::infy, \[Infinity]::indet];
    Switch[sides,
      OneSidedUpper, 
      B1 = B2 = B /. bsub /. b -> Max[b /. RootSearch[equ, {b, 0, 5}, FilterOptions[RootSearch,opts]]],
      OneSidedLower, 
      B1 = B2 = B /. bsub /. b -> Max[b /. RootSearch[eql, {b, 0, 5}, FilterOptions[RootSearch,opts]]],
      TwoSidedAsymmetrical, 
      B1 = B2 = B /. bsub /. b -> Max[b /. RootSearch[eq2a, {b, 0, 5}, FilterOptions[RootSearch,opts]]],
      TwoSidedEquiTailed, 
      B1 = B /. bsub /. b -> Max[b /. RootSearch[equ, {b, 0, 5}, FilterOptions[RootSearch,opts]]];   
      B2 = B /. bsub /. b -> Max[b /. RootSearch[eql, {b, 0, 5}, FilterOptions[RootSearch,opts]]]; If[B2<B1, {B1,B2} = {B2,B1}]
      ];
    On[Power::infy, \[Infinity]::indet];
    Return[{B1, B2}//Round];]

BootstrapResamples::clev = 
"ConfidenceLevel is expected to be between 0 and 1."
BootstrapResamples::cint =
"Warning: ConfidenceInterval is not specified correctly, default value TwoSidedAsymmetrical is used."



(* Extreme percentiles and BCa confidence intervals. *)

BootstrapCI[val_Real, boot1_List, jack1_List:{}, opts___?OptionQ] := Module[{meth, extr, cov, acap, theta, z0cap, alpha, alpha1, alpha2, bsdim, pos, indlo, indup, B1, B2, lo, up, j, bb, bbnew, boot, jack, sides},
   meth = BootstrapCIMethod /. {opts} /. Options[BootstrapCI]; Which[meth===ExtremePercentiles, extr=True, meth===BCaPercentiles, extr=False, True, Message[BootstrapCI::meth]; extr=True];
   If[boot1=={}||jack1=={}, Message[BootstrapCI::larg]; Return[$Failed]];
   cov = ConfidenceLevel /. {opts} /. Options[BootstrapCI]; If[cov <= 0 || 1 <= cov, Message[BootstrapCI::clev]; Return[$Failed]]; alpha=1-cov;
   boot = Sort[boot1];
   If[extr,
   sides = ConfidenceInterval /. {opts} /. Options[BootstrapCI];
   If[Intersection[{sides}, {OneSidedUpper, OneSidedLower, 
            TwoSidedAsymmetrical, TwoSidedEquiTailed}] === {}, 
      Message[BootstrapCI::cint]; sides = TwoSidedAsymmetrical];
   bb=BootstrapResamples[val, jack1, opts];
   If[Max[bb]!=Dimensions[boot][[1]], Message[BootstrapCI::bdim]; Return[$Failed]];
   B1=Min[bb]; B2=Max[bb];
   Switch[sides,
(* one sided and asymmetrical tw-sided *)
     OneSidedUpper, Return[{-Infinity,boot[[B1]]}],
     OneSidedLower, Return[{boot[[1]],Infinity}],
     TwoSidedAsymmetrical, Return[{boot[[1]],boot[[B1]]}],
     TwoSidedEquiTailed,
(* two-sided equi-tailed; see Lee's paper *)   
   If[B2==bb[[1]], lo=Factorial[bb[[2]]]*Factorial[bb[[1]]-bb[[2]]]/Factorial[bb[[1]]]*Sum[Factorial[j-1]/Factorial[bb[[2]]-1]/Factorial[j-bb[[2]]]*boot[[bb[[1]]-j+1]],{j,bb[[2]],bb[[1]]}]; up=boot[[bb[[1]]]]; Return[{lo,up}]]; If[B2==bb[[2]], lo=boot[[1]]; up=Factorial[bb[[1]]]*Factorial[bb[[2]]-bb[[1]]]/Factorial[bb[[2]]]*Sum[Factorial[j-1]/Factorial[bb[[1]]-1]/Factorial[j-bb[[1]]]*boot[[j]],{j,bb[[1]],bb[[2]]}]; Return[{lo,up}]]
],
(* BCa; see Efron's book *)
   jack=Sort[jack1]; theta=Apply[Plus,jack]/Dimensions[jack][[1]]; bsdim=Dimensions[boot][[1]]; acap=Apply[Plus,(theta-jack)^3]/6./(Apply[Plus,(theta-jack)^2])^1.5; pos=Position[boot,_?(#<val &)]; z0cap=Qn[N[Max[Flatten[pos],0]/bsdim]]; alpha1=Ncdf[z0cap+1./(1./(z0cap+Qn[alpha/2.])-acap)]; alpha2=Ncdf[z0cap+1./(1./(z0cap+Qn[1.-alpha/2.])-acap)]; indlo=Floor[(bsdim+1)*alpha1]; indup=bsdim+1-Floor[(bsdim+1)*(1.-alpha2)]; If[indlo<1||indlo>bsdim,indlo=Abs[indlo-1]]; If[indup<1||indup>bsdim,indup=Abs[indup-1]]; Return[{boot[[indlo]],boot[[indup]]}]
   ];
   ]

BootstrapCI::meth =
"Warning: BootstrapCIMethod is not specified correctly, default value ExtremePercentiles is used."
BootstrapCI::clev = 
"ConfidenceLevel is expected to be between 0 and 1."
BootstrapCI::larg = 
"Argument should be an nonempty list."
BootstrapCI::bdim = 
"Length of list with bootstrap replications does not match selected type of extreme percentile interval."
BootstrapCI::cint =
"Warning: ConfidenceInterval is not specified correctly, default value TwoSidedAsymmetrical is used."


(* Variance of endpoints of BCa intervals. *)

VarianceNonParametricBCaCIEndpoints[val_Real, boot1_List, jack1_List, bootin_List, opts___?OptionQ] := Module[{i, cov, test, boot2, bootci, out1, varlists, n, means, outvariance},
   If[Dimensions[boot1][[1]]!=Dimensions[bootin][[1]], Message[VarianceNonParametricBCaCIEndpoints::bdim]; Return[$Failed]];
   n = Dimensions[jack1][[1]]; out1={};
   cov = ConfidenceLevel /. {opts} /. Options[VarianceNonParametricBCaCIEndpoints]; If[cov <= 0 || 1 <= cov, Message[VarianceNonParametricBCaCIEndpoints::clev]; Return[$Failed]];
   For[i=1, i<=n, i++, 
   test = MapThread[MemberQ, {bootin, Table[i, {Dimensions[bootin][[1]]}]}]; 
   boot2 = boot1[[ Flatten[Position[test, True]] ]]; 
   bootci = BootstrapCI[val, boot2, jack1, BootstrapCIMethod -> BCaPercentiles, ConfidenceLevel -> cov]; 
   AppendTo[out1, bootci]];
   varlists=Transpose[out1]; means=Apply[Plus, varlists, 2]/n; outvariance = (n-1)/n Apply[Plus, (varlists-means)^2, 2];
   Return[outvariance];
] /; boot1!={} && bootin!={} && jack1!={}

VarianceNonParametricBCaCIEndpoints::bdim =
"Dimensions of arguments in positions 2 and 4 are expected to be the same."
VarianceNonParametricBCaCIEndpoints::clev = 
"ConfidenceLevel is expected to be between 0 and 1."



End[]
EndPackage[]