
% BB Nov 21 2006.
% BB Mar 08 2004.

\documentclass[20pt,landscape,footrule]{foils}
\usepackage[dvips]{graphicx}
\usepackage[pdftex]{geometry}
\geometry{headsep=3ex,hscale=0.9}
\usepackage[pdftex]{color} % for the colored block
\date{\color{blue} Singapore MCMC Meeting, March 2004}
\usepackage[ignore]{pause}
%\usepackage{pause}
\usepackage{background}
%\usepackage{pp4slide}
%\MyLogo{AV HEP,CSIT}
\LogoOff
\begin{document}
\def\r{\color{red}}
\def\b{\color{blue}}
\def\s{\color{black}}

\definecolor{bgblue}{rgb}{0.04,0.39,0.53}
\definecolor{pictcol}{rgb}{0.04,0.29,0.44}
\vpagecolor{bgblue}

\foilhead{\r Lecturenotes Statistics II -- Contents}\b \vskip 0.0pt
\begin{enumerate}
\item The Central Limit Theorem and Binning
\item Gaussian Error Analysis for Large and Small Samples
\item The Jackknife Approach

\end{enumerate}

\foilhead{\r The Central Limit Theorem and Binning }\b
How is the sum of two independent random variables
\begin{equation} \label{two_r} \s
  y^r\ =\ x_1^r + x_2^r\ . 
\end{equation}
distributed? We denote the probability density of $\s y^r$ by $\s g(y)$. 
The corresponding cumulative distribution function is given by 
$$\s G(y) = \int_{x_1+x_2 \le y} f_1(x_1)\ f_2(x_2)\ dx_1\ dx_2
        = \int_{-\infty}^{+\infty} f_1(x)\ F_2(y-x)\ dx $$
where $\s F_2(x)$ is the distribution function of the random variable 
$\s x_2^r$. We take the derivative and obtain the probability 
density of $\s y^r$
\begin{equation} \label{two_pd} \s
 g(y)\ =\ {dG(y)\over dy}\  
       =\ \int_{-\infty}^{+\infty} f_1(x)\ f_2(y-x)\ dx\ . 
\end{equation}
The probability density of a sum of two independent random variables 
is the {\r convolution of the probability densities} of these random 
variables.  

Example: Sums of uniform random numbers, corresponding to the sums 
of an uniformly distributed random variable $\s x^r\in (0,1]$:
(a) Let $\s y^r=x^r+x^r$, then
\begin{equation} \label{two_upd} \s
 g_2(y)\ =\ \cases{ y ~~~~~{\rm for}~~~ 0\le y\le 1, \cr
                    2-y  ~~{\rm for}~~~ 1\le y\le 2, \cr
                    0   ~~~{\rm elsewhere} . \cr}
\end{equation}
(b) Let $\s y^{\,r}=x^r+x^r+x^r$, then
\begin{equation} \label{three_upd} \s
 g_3(y)\ =\ \cases{ y^2/2 ~~~~~~{\rm for}~~ 0\le y\le 1, \cr
              (-2y^2+6y-3)/2 ~~{\rm for}~~ 1\le y\le 2, \cr
                 (y-3)^2/2 ~~~~{\rm for}~~ 2\le y\le 3, \cr
                  0   ~~~{\rm elsewhere} . \cr}
\end{equation}

\newpage
The convolution~(\ref{two_pd}) takes on a simple form in 
{\r Fourier space}. In statistics the {\r Fourier transformation} 
of the probability density is known as {\r characteristic function},
defined as the expectation value of $\s e^{itx^r}$:
\begin{equation} \label{chfct} \s
 \phi(t)\ =\ \langle e^{itx^r} \rangle\ =\ 
 \int_{-\infty}^{+\infty} e^{itx}\, f(x)\, dx\ .
\end{equation}
The characteristic function is particularly useful for investigating
sums of random variables, $\s y^r = x_1^r + x_2^r$:
$$\s \phi_y(t) = \left\langle e^{itx_1^r+itx_2^r} \right\rangle 
  =  \int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} 
  e^{itx_1}\,e^{itx_2}\,f_1(x_1)\,f_2(x_2)\ dx_1\,dx_2
  = \phi_{x_1}(t)\,\phi_{x_2}(t)\,.$$
{\r The characteristic function of a sum of random variables is the 
product of their characteristic functions.} The result generalizes
immediately to $\s N$ random variables
\begin{equation} \label{sum_r} \s
 y^r\ =\ x_1^r + \dots + x_N^r\ . 
\end{equation}
The characteristic function of $\s y^r$ is
\begin{equation} \label{N_cf} \s
 \phi_y(t)\ =\ \prod_{i=1}^N \phi_{x_i}(t)
\end{equation}
and the probability density of $\s y^r$ is the Fourier back-transformation
of this characteristic function
\begin{equation} \label{sum_pd} \s
  g(y)\ =\ {1\over 2\pi} \int_{-\infty}^{+\infty} dt\, e^{-ity}\,
  \phi_y(t)\ .
\end{equation}
The {\r probability density of the sample mean} is obtained as follows:
The arithmetic mean of $\s y^r$ is $\s \overline{x}^{\, r}=y^r/N$.
We denote the probability density of $\s y^r$ by $\s g_N(y)$ and the 
probability density of the arithmetic mean by $\s \widehat{g}_N
(\overline{x})$. They are related by 
\begin{equation} \label{mean_two_pd} \s
 \widehat{g}_N (\overline{x})\ =\ N\, g_N(N\overline{x})\ . 
\end{equation}

\newpage
This follows by substituting $\s y=N\overline{x}$ into $\s g_N(y)\,dy$:
$$\s 1 = \int_{-\infty}^{+\infty} g_N(y)\,dy = 
\int_{-\infty}^{+\infty} g_N(N\overline{x})\,2d\overline{x} = 
\int_{-\infty}^{+\infty} \widehat{g}_N(\overline{x})\,d\overline{x}\ .$$
Example:
\begin{figure} \begin{center} 
\includegraphics[height=0.45\textheight]{figures/mean_u.pdf}
\caption{\small Probability densities for the arithmetic means 
 of two and three uniformly distributed random variables, 
 $\s \widehat{g}_2(\overline{x})$ and $\s \widehat{g}_3(\overline{x})$, 
 respectively.} 
\end{center} \end{figure}

\newpage
This suggests that sampling leads to convergence of the mean
by reducing its variance. We use the characteristic function to 
understand the general behavior. The characteristic function of a 
sum of independent random variables is the product of their individual 
characteristic functions
\begin{equation} \label{sum_chfct} \s
 \phi_y(t)\ =\ \left[ \phi_x(t) \right]^N\ .
\end{equation}
The characteristic function for the corresponding arithmetic
average is
$$\s \phi_{{\overline x}}(t)\ =\ \int_{-\infty}^{+\infty} 
   d\overline{x}\, e^{it\overline{x}}\, \widehat{g}_N(\overline{x})\ 
=\ \int_{-\infty}^{+\infty} N d\overline{x}\, 
    e^{it\overline{x}}\, g_N(N\overline{x}) $$ 
$$\s =\ \int_{-\infty}^{+\infty} dy\, \exp \left( i\,{t\over N}\, y\right)\,
       g_N (y)\ . $$
Hence,
\begin{equation} \label{mean_chf} \s
 \phi_{{\overline x}}(t)\ =\ \phi_y \left( {t\over N} \right)\
 =\ \left[ \phi_x \left( {t\over N} \right) \right]^N\ .
\end{equation}

\newpage
{Example: The normal distribution.}

The characteristic function is obtained by Gaussian integration
\begin{equation} \label{gchf} \s
 \phi (t)\ =\ \exp \left( - {1\over 2} \sigma^2 t^2 \right)\ .
\end{equation}
Defining $\s y^r=x^r+x^r$ we have
\begin{equation} \label{2gchf} \s
 \phi_y(t)\ =\ \left[ \phi(t) \right]^2\ =\ 
 \exp \left( - {1\over 2} 2 \sigma^2 t^2 \right)\ .
\end{equation}
This is the characteristic function of a Gaussian with
variance $\s 2 \sigma^2$. We obtain the characteristic function of the
arithmetic average $\s {\overline{x}^{\,r}} = y^r/2$ by the
substitution $\s t\to t/2$:
\begin{equation} \label{mean_2gchf} \s
 \phi_{{\overline x}}(t)\ =\ 
 \exp \left( -{1\over 2} {\sigma^2\over 2} t^2 \right)\ .
\end{equation}
The variance is reduced by a factor of two.

\newpage
{\r\bf The Central Limit Theorem}

To simplify the equations we restrict ourselves to $\s {\widehat x}=0$. 
Let us consider a probability density $\s f(x)$ and assume that its 
moment exists, implying that the characteristic function is a least 
two times differentiable, so that
\begin{equation} \label{approx_chf} \s
 \phi_x (t)\ =\ 1\ -\ {\sigma_x^2 \over 2}\ t^2\ +\ {\cal O} (t^3)\ .
\end{equation}
The leading term reflects the the normalization of the probability density
% which makes the zeroth moment equal to one.
and the first moment is $\s \phi' (0)=\widehat{x}=0$. The characteristic 
function of the mean becomes
$$\s \phi_{\overline x} (t)\ =\ \left[ 1\ -\ {\sigma_x^2 \over 2 N^2} t^2\
 +\ {\cal O} \left( {t^3 \over N^3} \right) \right]^N\
 =\ \exp \left[ - {1\over 2}\,{\sigma_x^2 \over N}\, t^2 \right]
 +\ {\cal O} \left( {t^3 \over N^2} \right) \ .$$
\newpage

{\r {\bf\large The probability density of the arithmetic mean 
$\overline{x}^{\,r}$ converges towards the 
Gaussian probability density with variance}
\begin{equation} \label{mean_var}
  \sigma^2 (\overline{x}^{\,r})\ =\ {\sigma^2 (x^r) \over N}\ .
\end{equation} }

{A Counter example:} The Cauchy distribution provides an 
instructive, case for which the central limit theorem does not work. 
This is expected as its second moment does not exist. 

Nevertheless, the characteristic function of the Cauchy distribution 
exists.  For simplicity we take $\s \alpha=1$ and get 
\begin{equation} \label{cchf} \s
  \phi (t) =\ \int_{-\infty}^{+\infty} dx\, {e^{itx} \over
  \pi\, (1+x^2) }\ =\ \exp (-|t|)\ .
\end{equation}
The integration involves the residue theorem. Using 
equation~(\ref{mean_chf}) for the characteristic function 
of the mean of $\s N$ random variables, we find
\begin{equation} \label{mean_cchf} \s
  \phi_{\overline{x}} (t) =\ \left[ \exp \left( -{|t|\over N}\right)
  \right]^n\ =\ \exp (-|t|)\ .
\end{equation}
The surprisingly simple result is that the probability distribution for 
the mean values of $\s N$ independent Cauchy random variables agrees 
with the probability distribution of a single Cauchy random variable. 
Estimates of the Cauchy mean cannot be obtained by sampling. Indeed, 
the mean does not exist. 

\newpage
{\r \bf Binning}

The notion of introduced here should not be confused with 
histogramming! Binning means here that we group {\tt NDAT} data 
into {\tt NBINS} bins, where each binned data point is the arithmetic 
average of
$$\s {\tt NBIN = [NDAT/NBINS]} ~~({\rm Fortran\ integer\ division.})$$
original data points. Preferably {\tt NDAT} is a multiple of {\tt NBINS}. 
The purpose of the binning procedure is twofold:

\begin{enumerate}
\item When the the central limit theorem applies, the binned data
will become practically Gaussian, as soon as {\tt NBIN} becomes large
enough. This allows to apply Gaussian error analysis methods even
when the original are not Gaussian.

\item When data are generated by a Markov process subsequent events 
are correlated. For binned data these correlations are reduced
and can be neglected, once {\tt NBIN} is sufficiently large compared 
to the autocorrelation time.
\end{enumerate}

\newpage
{\bf Example:}

\begin{figure} \begin{center}
\includegraphics[height=0.7\textheight]{figures/bin_his.pdf}
\caption{\small Comparison of a histogram of 500 binned data with 
the normal distribution $\sqrt{(120/\pi)}\, \exp [-120\,(x-1/2)^2]$.
Each binned data point is the average of 20 uniformly distributed 
random numbers. Assignment {\tt a0108\_02}.}
\end{center} \end{figure}

\foilhead{\r Gaussian Error Analysis for Large and Small Samples}
\b \vskip -30.0pt

The central limit theorem underlines the importance of the normal 
distribution. Assuming we have a large enough sample, the 
arithmetic mean of a suitable expectation value becomes normally 
distributed and the calculation of the confidence intervals is reduced 
to studying the normal distribution. It has become the convention to 
use the {\r standard deviation of the sample mean}
\begin{equation} \label{mean_sdev} \s
  \sigma = \sigma (\overline{x}^{\,r})~~{\rm with}~~~
  \overline{x}^{\,r} = \frac{1}{N} \sum_{i=1}^N x_i^r
\end{equation}
to indicate its confidence intervals
$\s [{\widehat x}-n\sigma ,{\widehat x}+n\sigma ]$ (the dependence 
of $\s \sigma$ on $\s N$ is suppressed). For a Gaussian distribution 
the probability content $\s p$ of the confidence intervals
(\ref{mean_sdev}) to be
\begin{equation} \label{p_gauss} \s
 p = p(n) = G(n\sigma) - G(-n\sigma) = 
  {1\over \sqrt{2\pi}} \int_{-n}^{+n} dx\, e^{-{1\over 2} x^2}
  = {\rm erf} \left( {n\over \sqrt{2}} \right)\ .
\end{equation}

\vskip -20pt
\begin{table}[tbp]
\caption{Probability content $\s p$ of Gaussian confidence intervals
$\s [\widehat{x}-n\sigma ,\widehat{x}+n\sigma]$, $\s n=1,\dots ,6$, and
$\s q=(1-p)/2$. Assignment {\tt a0201\_01}.}
\centering
\begin{verbatim}

  n      1          2          3          4          5    
  p    .68        .95        1.0        1.0        1.0    
  q    .16        .23E-01    .13E-02    .32E-04    .29E-06 
\end{verbatim}
\label{tab_gcints}
\end{table}
In practice the roles of $\s \overline{x}$ and $\s \widehat{x}$ are 
interchanged: One would like to know the likelihood that the 
{\r unknown} exact expectation value $\s {\widehat x}$ will be 
in a certain confidence interval around the measured sample mean. 
The relationship
\begin{equation} \label{interchange} \s
 {\overline x}\in [{\widehat x}-n\sigma ,{\widehat x}+n\sigma]\
 \Longleftrightarrow\ {\widehat x} \in
 [ {\overline x}-n\sigma ,{\overline x}+n\sigma ]
\end{equation}
solves the problem. Conventionally, these estimates are quoted as
\begin{equation} \label{error_bar} \s
 {\widehat x}\ =\ {\overline x} \pm \triangle {\overline x}
\end{equation}
where the {\r error bar} $\s \triangle \overline{x}$ is often an 
{\r estimator} of the exact standard deviation.

\newpage
An obvious estimator for the variance 
$\sigma^2_x$ is
\begin{equation} \label{sample_variance1} \s
  (s'^{\,r}_x)^2\ =\ {1\over N} \sum_{i=1}^N (x_i^r
  - \overline{x}^{\,r})^2
\end{equation}
where the prime indicates that we shall not be happy with it, because
we encounter a {\bf bias}. An estimator is said to be biased when its 
expectation value does not agree with the exact result. In our case
\begin{equation} \label{variance_bias} \s
 \langle (s'^{\,r}_x)^2\rangle\ \ne\ \sigma^2_x\ . 
\end{equation}
An estimator whose expectation value agrees with the true expectation 
value is called {\r unbiased}. 
For the variance it is rather straightforward to construct an unbiased 
estimator $\s (s^r_x)^x$. The bias of the definition (\ref{sample_variance1}) 
comes from replacing the exact mean $\s {\widehat x}$ by its estimator
$\s \overline{x}^{\,r}$. The latter is a random variable, whereas the
former is just a number. 

\newpage
Some algebra shows that the desired {\r unbiased estimator of the 
variance} is given by
\begin{equation} \label{sample_variance} \s
 (s^r_x)^2\ =\ {N\over N-1}\, (s'^{\,r}_x)^2\ =\
{1\over N-1} \sum_{i=1}^N ( x_i^r - \overline{x}^{\,r} )^2\ .
\end{equation}
Correspondingly, the unbiased estimator of the variance of the 
sample mean is 
\begin{equation} \label{sample_mean_variance} \s
 (s^r_{\overline x})^2\ =\ {1\over N (N-1)} 
 \sum_{i=1}^N ( x_i^r - \overline{x}^{\,r} )^2\ .
\end{equation}

{\r\bf Gaussian difference test:}

In practice one is often faced with the problem to compare two 
different empirical estimates of some mean. How large must 
$\s D = {\overline x} - {\overline y}$ be in order to indicate a 
real difference? The quotient 
\begin{equation} \label{gdr} \s
  d^r\ =\ {D^r\over \sigma_D}
\end{equation}
is normally distributed with expectation zero and variance one,
so that
\begin{equation} \label{Pgd} \s
  P\ =\ P(|d^r| \le d)\ =\ G_0(d) - G_0(-d)\ =\ 1 - 2\, G_0(-d)\
     =\ {\rm erf} \left( {d\over\sqrt{2}} \right)\ .
\end{equation}
The {\r likelihood that the observed difference $\s 
|{\overline x}-{\overline y}|$ is due to chance} is defined to be
\begin{equation} \label{Qgood} \s
  Q\ =\ 1-P\ =\ 2\,G_0(-d)\ =\ 
  1-{\rm erf} \left( {d\over\sqrt{2}}\right)\ .
\end{equation}
If the assumption is correct, then $\s Q$ is a uniformly distributed 
random variable in the range $\s [0,1)$. Examples: 

\begin{table}[ht]
\caption{Gaussian difference tests (assignment {\tt a0201\_06}).} 
\medskip
\centering
\begin{tabular}{||c|c|c|c|c|c||}                            \hline
 ${\overline x}_1 \pm \sigma_{{\overline x}_1}$ & $1.0\pm 0.1$
 & $1.0\pm 0.1$ &$1.0\pm 0.1$&$1.0\pm 0.05$ & $1.000\pm 0.025$
\\ \hline
 ${\overline x_2} \pm \sigma_{{\overline x}_2}$ & $1.2\pm 0.2$
 &$1.2\pm 0.1$ &$1.2\pm 0.0$&$1.2\pm 0.00$ & $1.200\pm 0.025$
\\ \hline
 $Q$ & 0.37 & 0.16 &~~~0.046 & 0.000063 & $0.15\times 10^{-7}$
\\ \hline
\end{tabular}
\label{tab_gdt}
\end{table}

\newpage
{\r\bf Gosset's Student Distribution}

We ask the question: What happens with the Gaussian confidence 
limits when we replace the variance $\s \sigma^2_{\overline x}$ by 
its estimator $\s s^2_{\overline x}$ in statements like
$$\s { |\overline{x} - \widehat{x} | \over \sigma_{\overline x} }
  < 1.96 ~~{\rm with}~~ 95\% ~~{\rm probability.} $$
For sampling from a Gaussian distribution the answer was given by
Gosset, who published his article 1908 under the pseudonym 
{\it Student} in Biometrika. He showed that the distribution of 
the random variable
\begin{equation} \label{student_rv} \s
 t^r = { \overline{x}^{\,r} - \widehat{x} \over s^r_{\overline x} }
\end{equation}
is given by the probability density % \vskip -40pt
\begin{equation} \label{stud_pd} \s
   f(t) = {1\over (N-1)\, B(1/2,(N-1)/2)}\,
  \left( {1+{t^2\over N-1}} \right)^{-{N\over 2}}\ .
\end{equation} 
Here $\s B(x,y)$ is the beta function. The fall-off is a power law 
$\s |t|^{-f}$ for $\s |t|\to \infty$. Confidence probabilities of the 
Student distribution are:
\begin{table}[tbp]
% \caption{ Confidence probabilities for given confidence intervals, 
% see assignment~\ref{assign_student}.}
\centering
\begin{verbatim}  
     N \ S  1.0000     2.0000     3.0000     4.0000     5.0000    
     2      .50000     .70483     .79517     .84404     .87433    
     3      .57735     .81650     .90453     .94281     .96225    
     4      .60900     .86067     .94233     .97199     .98461    
     8      .64938     .91438     .98006     .99481     .99843    
    16      .66683     .93605     .99103     .99884     .99984    
    64      .67886     .95018     .99614     .99983     1.0000    
 INFINITY:  .68269     .95450     .99730     .99994     1.0000    
\end{verbatim}
\label{tab_stud1}
\end{table}
For $\s N\le 4$ we find substantial deviations from the Gaussian
confidence levels. Up to two standard deviations reasonable 
approximations of Gaussian confidence limits are obtained for 
$\s N\ge 16$ data. If desired, the Student distribution function
can always be used to calculate the exact confidence limits. 
% When the central limit theorem applies, we can {\bf bin a large 
% set of non-Gaussian data into 16 almost Gaussian data} and reduce 
% the error analysis to Gaussian methods.

\newpage
{\r \bf Student difference test}

This test takes into account that only a finite number of events are 
sampled.  Let the following (normal distributed) data be given
\begin{eqnarray} \s
\overline{x} ~~{\rm calculated\ from}~~ M ~~{\rm events},\
 {\it i.e.},~~ \sigma^2_{\overline x} &=& \sigma^2_x/M \\ \s
\overline{y} ~~{\rm calculated\ from}~~ N ~~{\rm events},\
 {\it i.e.},~~ \sigma^2_{\overline y} &=& \sigma^2_y/N
\end{eqnarray}
and an unbiased estimators of the variances 
\begin{equation} \label{partial_data} \s
s^2_{\overline x} = s^2_x/M = {\sum_{i=1}^M (x_i-\overline{x})^2
\over M\,(M-1)} ~~{\rm and}~~
s^2_{\overline y} = s^2_y/N = {\sum_{j=1}^N (y_j-\overline{y})^2
\over N\,(N-1)}\ .
\end{equation}
Under the {\r additional assumption} $\s \sigma^2_x=\sigma^2_y$ the 
probability 
\begin{equation} \s
 P(|\overline{x}-\overline{y}|>d)
\end{equation}
is is determined by the Student distribution function.
% in the same way as the probability of the Gaussian difference test 
% is determined by the normal distribution. 
\newpage

Examples for the Student difference test for  
$\s \overline{x}_1 = 1.00 \pm 0.05$ from $\s M$ data and 
$\s \overline{x}_2 = 1.20 \pm 0.05$ from $\s N$ data
(assignment {\tt a0203\_03}):

%
\begin{table}[ht]
\medskip
\centering
\begin{tabular}{||c|c|c|c|c|c|c|c||}    \hline
$M$& 512  & 32   & 16   & 16  & 4   & 3   & 2  \\
\hline
$N$& 512  & 32   & 16   &  4  & 4   & 3   & 2  \\
\hline
$Q$&0.0048&0.0063&0.0083&0.072&0.030&0.047&0.11\\
\hline
\end{tabular}
\label{tab_sdt}
\end{table}
%

The Gaussian difference test gives $\s Q=0.0047$. For $\s M=N=512$ the 
Student $\s Q$ value is practically identical with the Gaussian result, 
for $\s M=N=16$ it has almost doubled. Likelihoods above a {\s 5\% 
cut-off,} are only obtained for {\s $M=N=2$ (11\%) and $M=16$, $N=4$ 
(7\%).} The latter result looks a bit surprising, because its $\s Q$ 
value is smaller than for $\s M=N=4$.  The explanation is that for 
$\s M=16$, $\s N=4$ data one would expect the $\s N=4$ error bar to 
be two times larger than the $\s M=16$ error bar, whereas the estimated 
error bars are identical.

This leads to the question: Assume data are sampled from the same 
normal distribution, {\r when are two measured error bars consistent 
and when not?} Answered later by: $\s \chi^2$ Distribution, Error of 
the Error Bar, and Variance Ratio Test.

\newpage
{\r \bf The Jackknife Approach} 

Jackknife estimators allow to correct for the bias and the error 
of the bias. The method was introduced in the 1950s in papers by 
Quenouille and Tukey. 
The jackknife method is {\r recommended as the standard} for 
error bar calculations.  In unbiased situations the jackknife and 
the usual error bars agree. Otherwise the jackknife estimates are 
improvements, so that one cannot loose. \vskip -9pt

The unbiased estimator of the expectation value $\s \widehat{x}$ is
$$\s {\overline x} = {1\over N} \sum_{i=1}^N x_i $$
Normally bias problems occur when one estimates a non-linear 
function of $\s {\widehat x}$: 
\begin{equation} \s
 {\widehat f} = f({\widehat x})\ . 
\end{equation}
Typically, the bias is of order $\s 1/N$:
\begin{equation} \label{f_bias} \s
 {\rm bias}\ ({\overline f})\ =\ {\widehat f} -  
 \langle {\overline f}\rangle\ =\ {a_1 \over N} + {a_2 \over N^2} 
 + O({1\over N^3})\ . 
\end{equation}
Unfortunately, we lost the ability to estimate the variance 
$\s \sigma^2 ({\overline f}) = \sigma^2 (f) / N$ via the standard equation
\begin{equation} \label{f_bad_variance} \s
 s^2 ({\overline f})\ =\ {1\over N} s^2 (f)\ =\
 {1\over N\, (N-1)} \sum_{i=1}^N (f_i - {\overline f})^2\ ,  
\end{equation}
because $\s f_i=f(x_i)$ is not a valid estimator of $\s {\widehat f}$. 
Also it is in non-trivial applications almost always a bad idea to to 
use standard error propagation formulas with the aim to deduce  
$\s \triangle {\overline f}$ from $\s \triangle {\overline x}$. 
Jackknife methods are not only easier to implement, but also more 
precise and far more {\r robust}.

The error bar problem for the estimator $\s \overline{f}$ is 
conveniently overcome by using {\r jackknife estimators} 
$\s {\overline f^J}$, $\s f_i^J$, defined by
\begin{equation} \label{jackknife_estimators} \s
 {\overline f^J}\ =\ {1\over N} \sum_{i=1}^N f^J_i
 ~~{\rm with}~~ f_i^J\ =\ f(x_i^J) ~~{\rm and}~~ 
   x_i^J\ =\ {1\over N-1} \sum_{k\ne i} x_k\ . 
\end{equation}
The estimator for the variance $\s \sigma^2 ({\overline f^J})$ is
$$\s s^2_J ({\overline f^J})\ =\ {N-1 \over N} 
  \sum_{i=1}^N (f^J_i - {\overline f^J})^2\ .   $$
Straightforw algebra shows that in the unbiased case the estimator 
of the jackknife variance reduces to the normal variance. 

Notably only of order $\s N$ (not $\s N^2$) operations are needed to
construct the jackknife averages $\s x_i^J,\ i=1,\dots ,N$ from the
original data.

\noindent The bias of each jackknife estimator $\s {\overline f^J}$ is 
also of order $\s 1/N$:
$$\s {\rm bias}\ ({\overline f^J}) =
  {\widehat f} - \langle {\overline f^J} \rangle\ = {a_1\over N-1} 
  + {a_2 \over (N-1)^2} + O \left({1\over N^3}\right)\ .  $$
As the bias converges faster than the statistical error and is then
normally ignored. In exceptional situations the constant $a_1$ may 
be large. Bias corrections should be done when the bias reaches the 
order of magnitude of the statistical error bar.
An estimator for the bias is easily constructed:
$$\s a_1 = N\,(N-1)\, \left(\langle {\overline f} \rangle 
  - \langle{\overline f^J}\rangle\right)+O\left({1\over N}\right) $$
and, therefore,
$$\s {\overline b}\ =\ (N-1)\ ({\overline f}-{\overline f^J}) $$
defines an estimator for the bias up to $\s O\left({1/N^2}\right)$
corrections. 


\noindent In our Fortran code the calculation of jackknife error bars 
is implemented by the subroutines {\s\tt datjack.f} and {\s\tt stebj0.f} 
of {\s\tt ForLib}. The jackknife data $\s x^J_i$, $\s (i=1,\dots ,N)$ 
of equation (\ref{jackknife_estimators}) are generated by a call to
$$\s  {\tt DATJACK(N,X,XJ)} $$
where {\s \tt X} is the array of the original data and {\s \tt XJ} the 
array of the corresponding jackknife data (estimators). This is done in 
$\s N$ ({\r not} $\s N^2$) steps. An array {\s \tt FJ}, corresponding 
to $\s f^J_i$ has then to be defined by the user. Relying on the 
jackknife variance definition  a subsequent call to
$$\s  {\tt STEBJ0(N,FJ,FM,FV,FE)} $$
uses the array {\s \tt FJ} as input to calculate the jackknife mean 
{\s\tt FJM}, variance {\tt FJV} and error bar {\tt FJME} of the mean. 
The subroutine {\s\tt bias.f} of {\s\tt ForLib} returns the bias.
 
\end{document}
