
% BB Sep 20 2008.
% 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

\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} \r
 \widehat{g}_N (\overline{x})\ =\ N\, g_N(N\overline{x})\ . 
\end{equation}
Proof:

\newpage
Substitute $\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})\,N d\overline{x} = 
\int_{-\infty}^{+\infty} \widehat{g}_N(\overline{x})\,d\overline{x}\ .$$
Example:
\begin{figure} \begin{center} 
\includegraphics[height=0.48\textheight]{figures/mean_u.pdf}
\caption{\small Probability densities for the arithmetic means 
 of two to six 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
{\r 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 
$\s \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} }

{\r 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. 

{\r Generalized limit theorems} exist for broad distributions, like 
those given by a probability density
$$\s f(x) = \frac{a}{2 (1+|x|)^{1+a}},\ 0<a\le 2\,. $$
This goes sometimes under the name Levy statistics, but is in the Russian
literature more closely associated with work by Khintchine.


\foilhead{\r Binning }\b \vskip -22pt

The notion of introduced here should not be confused with 
histograming! 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 $\s \sqrt{(120/\pi)}\, \exp [-120\,(x-1/2)^2]$.
Each binned data point is the average of 20 uniformly distributed 
random numbers. Assignment {\r\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 {\r\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 $\s \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\rm 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}
Our {\r workhorse} for error bar calculations from data is the
routines {\s\tt steb0.f} of {\s\tt ForLib} used via
$$\s\tt CALL\ STEB0(NDAT,DAT,DATM,DATV,DATME)\ .$$
See today's {\r classwork.}

\foilhead{\r Gaussian difference test} \vskip -34pt

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 data set are sampled from identical Gaussian distributions, 
$\s Q$ is a uniformly distributed random variable in the range 
$\s [0,1)$. Example: See today's 
{\r classwork.}

% \begin{table}[ht]
% \caption{Gaussian difference tests (assignment {\r\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}

\foilhead{\r Gosset's Student Distribution} \vskip -34pt

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} 

\begin{figure} \begin{center} 
\includegraphics[height=0.78\textheight]{figures/stud_pd.pdf}
\caption{\small Probability densities of some student distributions.} 
\end{center} \end{figure}

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]
\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
{\s \begin{eqnarray} 
\overline{x} ~~{\rm calculated\ from}~~ M ~~{\rm events},\
 {\it i.e.},~~ \sigma^2_{\overline x} &=& \sigma^2_x/M \\ 
\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. Example: Todays's
{\r classwork.}
% in the same way as the probability of the Gaussian difference test 
% is determined by the normal distribution. 
\newpage

\noindent {\bf More examples:} 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 {\r\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 $\s M=16$, 
$\s 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.

\foilhead{\r $\chi^2$ Distribution, Error of the Error Bar, 
             Variance Ratio Test} \vskip -30pt
\noindent For normally distributed data with mean zero the random 
variable
\begin{equation} \label{chi2_r} \s
 (\chi^r)^2 = \sum_{i=1}^f (y_i^r)^2\ ,
\end{equation}
defines the {\bf ${\s \chi^2}$ distribution} of $\s f$ degrees of 
freedom. The study of the variance $\s (s^r_x)^2$ of a Gaussian sample 
can be reduced to the $\s \chi^2$-distribution with $\s f=N-1$
\begin{equation} \label{chi2_s2} \s
 (\chi_f^r)^2 = {(N-1)\, (s^r_x)^2 \over \sigma^2_x } = 
 \sum_{i=1}^N { (x^r_i - \overline{x}^{\,r})^2 \over \sigma^2_x }\ .
\end{equation}
The probability density of {\bf $\s \chi^2$ per degree of freedom 
(pdf)} is
\begin{equation} \label{chi2_pdf_pd} \s
 f_N(\chi^2)\ =\ N f(N\chi^2)\ = 
{a\ e^{- a \chi^2}\ \left( a \chi^2 \right)^{a-1} \over
\Gamma (a) } ~~{\rm where}~~ a={N\over 2}\ . 
\end{equation}
% \newpage
For $\s N=1,2,\dots ,20$ this probability density is plotted in
the figure and we can see the central limit theorem at work. 
Picking the curves at $\s \chi^2/N=1$, increasing $\s f_N(\chi^2)$ 
values correspond to increasing $\s N=1,2,\dots,20$.
Assignment {\r\tt a0202\_03}.

\begin{figure} \begin{center}
\includegraphics[height=0.7\textheight]{figures/chi2_pd.pdf}
\end{center} \end{figure}

\newpage

{\bf\r The Error of the Error Bar}

{\bf For normally distributed data} the number of data alone determines
the errors of error bars, because the $\s \chi^2$ distribution is
exactly known. One does not have to rely on estimators!
Confidence intervals for variance estimates $\s s^2_x=1$ 
from {\tt NDAT} data (assignment {\r\tt a0204\_01}):

\begin{table}[tbp]  
\centering
\begin{verbatim}
 
                      q         q         q       1-q       1-q

    NDAT=2**K      .025      .150      .500      .850      .975
 
       2    1      .199      .483     2.198    27.960  1018.255
       4    2      .321      .564     1.268     3.760    13.902
       8    3      .437      .651     1.103     2.084     4.142
      16    4      .546      .728     1.046     1.579     2.395
      32    5      .643      .792     1.022     1.349     1.768
    1024   10      .919      .956     1.001     1.048     1.093
   16384   14      .979      .989     1.000     1.012     1.022
\end{verbatim}
\label{tab_var_e}
\end{table}

\foilhead{\r Variance ratio test ($\s F$-test)}

We assume that two sets of normal data are given together with
estimates of their variances $\s \sigma^2_{x_1}$ and {\s
$\sigma^2_{x_2}$: $\left( s^2_{x_1},N_1\right)$ and
$\s \left(s^2_{x_2},N_2\right).$} We would like to test 
whether the ratio $\s F={s^2_{x_1}\over s^2_{x_2}}$ 
differs from $\s F=1$.  With $\s f_i=N_i-1,\ i=1,2$
$$\s {\rm The\ probability}~~ {f_1\over f_2}\,F  < w ~~{\rm is}~~
 H(w) = 1 - B_I \left( {1\over w + 1}, {1\over 2}\,f_2,
 {1\over 2}\,f_1 \right)\ . $$
{\r This test allows us later to compare the efficiency of MC algorithms.} 

\noindent Examples: Today's {\r classwork.}

%  (assignment {\r\tt a0505\_02}):
% \begin{table}[ht] \centering
% \begin{tabular}{||c|c|c|c|c|c|c|c|c||}    \hline
% $\triangle \overline{x}_1$ 
% & 1.0& 1.0&  1.0 &1.0 & 1.0 &1.0&  1.0& 1.0\\ \hline
% $N_1$&  16&  16&   64 &1024& 2048& 32& 1024& 16 \\ \hline
% $\triangle \overline{x}_2$
%     & 1.0& 1.0&  1.0 &1.05& 1.05&2.0&  2.0& 2.0\\ \hline
% $N_2$& 16 &   8&   16 &1024& 2048&  8&  256& 16 \\ \hline
% $Q$&  1.0& 0.36& 0.005&0.12&0.027&0.90&0.98&0.01\\ \hline
% \end{tabular} \label{tab_Ftest} \end{table}

\foilhead{\r Statistical Bootstrap for Mean Values} \vskip -20pt \b

This approach is best done in three steps:

\begin{enumerate}

\item Block (bin) the data into {\tt\s NDAT} blocked data.

\item Sample {\r without replacement} from the blocked data a sample 
      of {\tt\s NBTR} mean estimates, the \r bootstrap sample. \b 
      See subroutine {\tt\s bootstrp.f}.

\item Sort the bootstrp sample and estimate confidence limits from
      the Cumulative Distribution Function (CDF). See subroutine 
      {\tt\s bootstrp.f}.

\end{enumerate}

\r Example: \b See classwork (pi estimate again).

Typically, the method is useful when one has relatively few 
independent data and/or the distribution of the estimated
mean values is not Gaussian.

\end{document}
