Quit[]

(* The Jack-knife *)

For details see: B.Efron, The jackknife, the bootstrap, and other resampling plans, In Society of Industrial and Applied Mathematics
CBMS - NSF Monographs, 38, (1982).

(* Various options *)
SetDirectory[NotebookDirectory[]];
Needs["ErrorBarPlots`"];
Off[CompiledFunction::cflist];

(* Some fake data with noise *)
SeedRandom[1264645];
(* A model *)
model[x_,a_,b_]:=a x^2+b x
(* Mock data with noise *)
data=Table[{x,model[x,3,-1]+RandomVariate[NormalDistribution[0,.1]],RandomReal[{.05,.1}]},{x,0,1,0.02}];
ndata=Length[data];
(* Plot of the data and model *)
Show[ErrorListPlot[data,PlotRange->All],Plot[model[x,3,-1],{x,0,1},PlotStyle->Red],Frame->True,FrameLabel->{"x","y"},Axes->False]

(* Define the chi^2 as a function of the data and the parameters a,b *)
chi2[dat_,a_,b_]:=Sum[((dat[[i,2]]-model[dat[[i,1]],a,b])/dat[[i,3]])^2,{i,1,Length[dat]}]

(* Delete the ith point each time *)
Table[Delete[Range[5],i],{i,1,Length[Range[5]]}]//MatrixForm

(* Launch some Kernels to speed up the work *)
LaunchKernels[7];

(* Sample with replacement and make the new datasets *)
mockdata=Table[Delete[data,i],{i,1,Length[data]}];

(* The best-fit of the real data *)
fmin0=FindMinimum[chi2[data,a,b],{a,1},{b,1}]

(* The best-fit of all of the mock data // This in practice will take the most time *)
fminmock=ParallelTable[FindMinimum[chi2[mockdata[[j]],a,b],{a,1},{b,1}],{j,1,Length[mockdata]},DistributedContexts:>Automatic];

(* Gather and plot the parameters (a,b) *)
mock2D=fminmock[[All,2,{1,2},2]];
amock=fminmock[[All,2,1,2]];
bmock=fminmock[[All,2,2,2]];
ListPlot[mock2D,Frame->True,Epilog->{PointSize[0.02],Red,Point[fmin0[[2,All,2]]]},FrameLabel->{"a","b"},AspectRatio->1]

fmin0

(* The mean values *)
{Mean[fminmock[[All,1]]],am=Mean[amock],bm=Mean[bmock]}

(* ... and the errors *)
errs={(Sqrt[((ndata-1)/ndata)Sum[(amock[[i]]-am)^2,{i,1,Length[amock]}]]),(Sqrt[((ndata-1)/ndata)Sum[(bmock[[i]]-bm)^2,{i,1,Length[bmock]}]])}

(* Fisher Matrix estimation of errors *)
p={a,b};
(* Find the best-fit by minimizing the (χ^2) *)
fmin=FindMinimum[chi2[data,a,b],{a,1},{b,1}]
fmin[[1]]/(ndata-2)
(* Fisher matrix *)
Fij=Table[(1/2)D[chi2[data,a,b],p[[i]],p[[j]]]/.fmin[[2]],{i,1,Length[p]},{j,1,Length[p]}];
Print["(Subscript[F, ij])= ",Fij//MatrixForm]
(* Covariance matrix *)
Cij=Inverse[Fij];
Print["(Subscript[C, ij])= ",Cij//MatrixForm]
(* Parameter errors *)
Print["1σ errors = ",(Sqrt[Cij])//Diagonal]

Print["Fisher matrix errors = ",(Sqrt[Diagonal[Cij]])]
Print["Jack-knife errors = ",errs]

(* The agreement ain't that bad! *)

(* end *)