\noindent A possible drawback of the ordinary correlation coefficient $\rho$ for two real random variables $X$ and $Y$ is that $\rho=0$ does not imply independence. In this paper a correlation coefficient $\rho^*$ is introduced which assumes values between zero and one, equalling zero iff the two variables are independent and equalling one iff the two variables are linearly related. For a $2\times 2$ contingency table, $\rho^*$ and the squared ordinary correlation coincide. Two interpretations of $\rho^*$ are given. First, it is shown that $\rho^*$ is a weighted sum of squared correlations between marginal eigenfunctions. For uniform, logistic and exponential marginals closed form expressions for the weights and eigenfunctions are given and for other marginals it is shown how numerical approximations can be obtained. Fr{\'e}chet bounds are discussed for the component correlations. Second, it is demonstrated that $\rho^*$ can be viewed as a weighted average of phi-square coefficients. The weights are plotted for several marginal distributions, indicating how the tails of the distributions are weighted. Estimators of $\rho^*$ based on both U- and V-statistics are discussed and their asymptotic distributions under independence, which depend on the marginal distributions of the variables, are derived. It is shown that as a special case a new Cram{\'e}r-von Mises test for $K$ ordered samples is obtained, which reduces to the ordinary Cram{\'e}r-von Mises test for two samples. The paper includes an extensive consideration of grade versions of $\rho^*$, and discrete distributions are given ample attention.
abs_cor.tex
%\documentclass[11pt]{article}\input{defs}

%\documentclass{aoscm}\input{defs_aos}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Annals of Statistics - Template file for LaTeX 2e
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\documentclass{aoscm}
%\documentclass[dblspace,leqno]{aoscm}% <-- gives double-spaced document/manuscript

\usepackage{apacite}



\raggedbottom
\usepackage{apacite,latexsym,graphics,epsfig,subfigure}

\def\id{\mathop{\rm id}\nolimits}
\def\sq{\mathop{\rm sq}\nolimits}
\def\sign{\mathop{\rm sign}\nolimits}
\def\half{{\scriptstyle\frac{1}{2}}}
\def\fourth{{\scriptstyle\frac{1}{4}}}
\def\threefourth{{\scriptstyle\frac{3}{4}}}

\newcommand{\cov}{{\mbox{\rm cov}}}
\newcommand{\var}{{\mbox{\rm var}}}

\newcommand{\tr}{\mbox{\rm tr}}
\newcommand{\erf}{\mbox{\rm erf}}
\newcommand{\erfi}{\mbox{\rm erfi}}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\setcounter{page}{1} % <-- To fill at IMS or the typesetter
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%% Definitions %%%%%%%%%%%%%%%%%%
%\def\theequation{\arabic{equation}} % If you don't want the equations counter
                                     % to be reset with each section,
                                     % outcomment this definition
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\font\msbmx=MSBM10                                       % \emptyset should be changed to \varnothing
%\font\msbmvii=MSBM7                                      % in the paper and this will give the proper symbol
%\font\msbmv=MSBM5
\def\varnothing{\mathchoice{\mbox{\msbmx\char'077}}%
{\mbox{\msbmx\char'077}}{\mbox{\msbmvii\char'077}}{\mbox{\msbmv\char'077}}}%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\def\Cov{\mathop{\rm Cov}\nolimits}% % A sample of proper declaration of
\def\Argmax{\mathop{\rm Arg\,Max}\limits}%    a math operator.
\def\Var{\mathop{\rm Var}\nolimits}%
%%%%%%%%%%%% Theorems %%%%%%%%%%%%%%%%%%%%%
% The AOS style allows for theorem-like environments to be numbered throughout
% the paper or to start from one each section (carrying in that case both the
% section number and "theorem" number). If you prefer the first method, declare
% all the environments you need as:
%   \newcounter{<counter_name>}
%   \newtheorem{<counter_name>}{<printed_name>},
% for example,
%   \newcounter{prop}
%   \newtheorem{prop}{Proposition}
% and that will appear as:
%
% Proposition 2. ...
%
% If you want to reset the counter with each new section, use the following
% code:
%  \newcounter{<counter_name>}
%  \newtheorem{<counter_name>}{<printed_name>}[section],
% for example,
%   \newcounter{prob}
%   \newtheorem{prob}{Problem}[section]
% and that will appear as:
%
% Problem 3.1. ...
%
% We list just some basic counters and leave for the user to decide which
% counters will be used, and should the above mentioned hints about
% declaring the theorem-like environments be used in one or the other way.
\newtheorem{theorem}{Theorem}
\newtheorem{lemma}{Lemma}
\newtheorem{definition}{Definition}
\newtheorem{corollary}{Corollary}
\newtheorem{proposition}{Proposition}
% With proof it is slightly different, and the proof environment does
% expect an argument showing the exact text that will start the proof,
% since this considerably varies from case to case. For example
%  \proof{Proof of the Modified Cramer-Rao Thoerem.}...
%  \endproof,
% where <Proof of the Modified Cramer-Rao Thoerem.> is the actual text to
% start the proof.
% The \proof command can be used also with other similar environments set
% in roman type (not italic as theorems, lemmas, etc.) as Remark, where
% also a halmos (open box) at the end is required:
%   \proof{Remark.}...\endproof
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% For References AOS provides two possible styles. The standard one
% numbers the references setting the numbers in square brackets. The
% other style that uses only the authors names can be declared by the
% command \NONUMBIB: If you want to use that style, outcomment the
% following line:
\NONUMBIB






\newcommand{\ind}{\perp\hspace*{-2.1mm}\perp}
\def\bR{{\bf R}}
\newtheorem{example}{Example}






%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{document}
%%%%%%%%%%%%%%%%

% On the opening page AOS style requires some unmarked footnotes in the
% following order:
% 1. Received line (no footnotemark)
% 2. One or more numerated footnotes showing author(s) grants, etc.
% 3. AMS Classification line (no footnotemark)
% 4. Keywords line (no footnotemark)
% To achieve this, a macro with up to nine parameters called
% \SPECFNSYMBOL is provided, where the user has to list the used markers
% in correct order. For example, if you need two footnotes for authors
% "thanks", you should use
\SPECFNSYMBOL{1}{2}{}{}{}{}{}{}{}%
% where the fourth and fifth {} are "placeholders" for the
% unmarked footnotes. To reset the footnote style back to normal at the
% beginning of the article body, use there the command
% \BACKTONORMALFOOTNOTE{3} (see below), "3" indicating that the
% first following footnote (if any) will carry the mark "3".
% To set (almost) properly the front matter in Computer Modern fonts,
% use now the command:
\AOSMAKETITLE
% and after doing so fill in the blanks (some will be done by the
% Publisher at the typesettimg stage --- leave X's or 0's) following
% the sample below:
% \AOSyr{1999}      % <-- actual numbers
% \AOSvol{00}       % <-- will be supplied
% \AOSno{00}        % <-- at typesetting
% \AOSpp{000--000}  % <-- stage
% \AOSReceived{}
% \AOSAMS{}
% \AOSKeywords{}
% \AOStitle{}%  % with a possible "thanks" footnote
% \AOSauthor{}% % with a additional "thanks" footnotes
% \AOSaffil{}%  % Just the short version of affiliation. The full address
%               % is at the end of the article.
% \AOSlrh{}%    % Author(s) (in all caps) (see the sample below)
% \AOSrrh{}%    % Shortened title for the running head (in all caps)
% \AOSAbstract{}% The text of the abstract follows in the braces.

 \AOSyr{2006}
 \AOSvol{00}
 \AOSno{00}
 \AOSpp{000--000}
 \AOSReceived{Received }
 \AOSAMS{Primary 62H20; secondary 62E99}
 \AOSKeywords{correlation, covariance, test of independence,
spectral decomposition, eigenvalues, eigenfunctions, Hilbert-Schmidt
operator, Fr{\'e}chet bounds, contingency tables, phi-square,
canonical correlation, Cram{\'e}r-von Mises tests, rank tests,
Fredholm integral equation of the second kind.}
 \AOStitle{ON A NEW CORRELATION COEFFICIENT, ITS ORTHOGONAL DECOMPOSITION AND ASSOCIATED TESTS OF INDEPENDENCE}%
 \AOSauthor{Wicher Bergsma\footnote{Supported by The Netherlands
Organization for Scientific Research (NWO), Project Number
400-20-001. }}
 \AOSaffil{London School of Economics and Political Science}%
 \AOSlrh{Wicher Bergsma}%
 \AOSrrh{A new correlation coefficient}%
 \AOSAbstract{A possible drawback of the ordinary correlation coefficient $\rho$
for two real random variables $X$ and $Y$ is that zero correlation
does not imply independence. In this paper we introduce a new
correlation coefficient $\rho^*$ which assumes values between zero
and one, equalling zero iff the two variables are independent and
equalling one iff the two variables are linearly related. The
coefficients $\rho^*$ and $\rho^2$ are shown to be closely related
algebraically, and they coincide for distributions on a $2\times 2$
contingency table. We derive an orthogonal decomposition of $\rho^*$
as a positively weighted sum of squared ordinary correlations
between certain marginal eigenfunctions. Estimation of $\rho^*$ and
its component correlations and their asymptotic distributions are
discussed, and we develop visual tools for assessing the nature of a
possible association in a bivariate data set. The paper includes
consideration of grade (rank) versions of $\rho^*$ as well as the
use of $\rho^*$ for contingency table analysis. As a special case a
new generalization of the Cram{\'e}r-von Mises test to $K$ ordered
samples is obtained.
}%

%%%%%%%%%%
% Standard LATEX command \maketitle will call the AOS style and set the
% top of the first page of the article in the proper style (but not in
% the proper font).
\maketitle
%%%%%%%%%%

% As discussed above, this resets the footnote counter and returns to
% the normal style of automatic numbering the footnotes (they are
% particularly not very welcomed in AOS).
\BACKTONORMALFOOTNOTE{3}




%\pagestyle{headings}
% \markright{Wicher Bergsma}
% \markleft{A new correlation coefficient}


%\pagestyle{headings}


%\tableofcontents


\noindent{\bf Contents:}

%\input{abs_cor.toc}

\contentsline {section}{\numberline {1}Introduction}{2}
\contentsline {section}{\numberline {2}Properties of the kernel
function $h_F$}{4} \contentsline {subsection}{\numberline {{\rm
2}.{\rm 1}}Key properties of $h_F$}{5} \contentsline
{subsection}{\numberline {{\rm 2}.{\rm 2}}Spectral decomposition of
$h_F$}{9} \contentsline {subsection}{\numberline {{\rm 2}.{\rm
3}}Obtaining the eigensystem in the discrete case}{13} \contentsline
{subsection}{\numberline {{\rm 2}.{\rm 4}}Obtaining the eigensystem
in the continuous case}{14} \contentsline {subsection}{\numberline
{{\rm 2}.{\rm 5}}Discrete approximation of the continuous case}{17}
\contentsline {subsection}{\numberline {{\rm 2}.{\rm 6}}Relation to
Anderson-Darling kernel}{18} \contentsline {section}{\numberline
{3}Properties of $\kappa $ and $\rho ^*$}{18} \contentsline
{subsection}{\numberline {{\rm 3}.{\rm 1}}Key properties of $\kappa
$ and $\rho ^*$}{19} \contentsline {subsection}{\numberline {{\rm
3}.{\rm 2}}Orthogonal decomposition}{22} \contentsline
{subsection}{\numberline {{\rm 3}.{\rm 3}}Parameterization of the
likelihood}{24} \contentsline {subsection}{\numberline {{\rm 3}.{\rm
4}}Fr{\'e}chet bounds for component correlations}{24} \contentsline
{section}{\numberline {4}Estimation and tests of independence}{27}
\contentsline {subsection}{\numberline {{\rm 4}.{\rm 1}}U and V
statistic estimators of $\kappa $}{27} \contentsline
{subsection}{\numberline {{\rm 4}.{\rm 2}}Permutation tests}{28}
\contentsline {subsection}{\numberline {{\rm 4}.{\rm 3}}Asymptotic
distribution of estimators under independence}{29} \contentsline
{subsection}{\numberline {{\rm 4}.{\rm 4}}Bonferroni corrections for
testing significance of component correlations}{30} \contentsline
{section}{\numberline {5}Grade versions of $\kappa $ and $\rho ^*$,
copulas, and rank tests}{31} \contentsline {subsection}{\numberline
{{\rm 5}.{\rm 1}}Rank tests}{32} \contentsline
{subsection}{\numberline {{\rm 5}.{\rm 2}}A new class of $K$-sample
Cram{\'e}r-von Mises tests as a special case}{32} \contentsline
{subsection}{\numberline {{\rm 5}.{\rm 3}}$\kappa $ as a weighted
$\phi $-coefficient}{33} \contentsline {section}{\numberline {6}Data
analysis: investigating the nature of the association}{34}
\contentsline {subsection}{\numberline {{\rm 6}.{\rm 1}}Some
artificial data sets}{35} \contentsline {subsection}{\numberline
{{\rm 6}.{\rm 2}}Mental health data}{41} \contentsline
{subsection}{\numberline {{\rm 6}.{\rm 3}}Norwegian stock
exchange}{43} \contentsline {subsection}{\numberline {{\rm 6}.{\rm
4}}Discussion}{43}


\section{Introduction}\label{sec1}

We introduce a correlation coefficient $\rho^*$ which has the
potential advantage compared to the ordinary correlation $\rho$ that
it detects arbitrary forms of association between two real random
variables $X$ and $Y$. In fact $\rho^*$, to be defined below, can be
viewed as a simple modification of $\rho^2$, as we now show. The
ordinary covariance is defined as
\begin{eqnarray*}
    \cov(X,Y) = E(X-EX)(Y-EY)
\end{eqnarray*}
and the ordinary correlation as
\begin{eqnarray*}
    \rho(X,Y) = \frac{\cov(X,Y)}{\sqrt{\cov(X,X)\cov(Y,Y)}}
\end{eqnarray*}
Now suppose that $Z$, $Z_1$ and $Z_2$ are iid with distribution
function $F$, that $X$ and $Y$ have marginal distributions $F_1$ and
$F_2$, and that $(X_1,Y_1),$ and $(X_2,Y_2)$ are iid replications of
$(X,Y)$. Then with
\begin{eqnarray}
    u_F(z_1,z_2) = (z_1-EZ)(z_2-EZ) = E(z_1-Z_1)(z_2-Z_2)
    \label{udef}
\end{eqnarray}
it is easy to verify that
\begin{eqnarray}
    \cov(X,Y)^2 = E
    u_{F_1}(X_1,X_2)u_{F_2}(Y_1,Y_2)
    \label{covdef}
\end{eqnarray}
Now straightforward algebra based on the left hand side
of~(\ref{udef}) shows that we can rewrite $u_F$ as
\begin{eqnarray*}
    u_F(z_1,z_2) =
   -
   \half E
   {\mbox{\Large $($}}\,
   |z_1-z_2|^2 - |z_1-Z_2|^2 - |Z_1-z_2|^2 + |Z_1-Z_2|^2
   \,{\mbox{\Large )}}
\end{eqnarray*}
Replacing the squares in $u_F$ by absolute values then gives
\begin{eqnarray}
   h_F(z_1,z_2) =
   -
   \half
   E{\mbox{\Large $($}}\,
   |z_1-z_2| - |z_1-Z_2| - |Z_1-z_2| + |Z_1-Z_2|
   \,{\mbox{\Large )}}
   \label{hdef1}
\end{eqnarray}
and we can define a new `covariance' $\kappa$ by replacing $u$ by
$h$ in~(\ref{covdef}):
\begin{eqnarray*}
   \kappa(X,Y) = Eh_{F_1}(X_1,X_2)h_{F_2}(Y_1,Y_2)
\end{eqnarray*}
Now we can also define
\begin{eqnarray*}
   \rho^*(X,Y)
   =
   \frac
   {\kappa(X,Y)}
   {\sqrt{\kappa(X,X)\kappa(Y,Y)}}
\end{eqnarray*}
Thus, whereas the squared covariance and $\rho^2$ are based on
squared differences, $\kappa$ and $\rho^*$ are based on absolute
differences. In this paper we demonstrate the perhaps surprising
result that $0\le \rho^*(X,Y)\le 1$, such that $\rho^*(X,Y)=0$ iff
$X$ and $Y$ are independent, and $\rho^*(X,Y)= 1$ iff $X$ and $Y$
are linearly related. A further main result we give is an orthogonal
decomposition of $\rho^*$ in terms of component correlations between
eigenfunctions of $h_{F_1}$ and $h_{F_2}$.


Based on their formulas, the following statistical interpretation of
$\rho^2$ and $\rho^*$ can be given: they measure how much two $X$
observations which are `far' apart tend to occur with $Y$
observations which are `far' apart, and similarly how much two $X$
observations which are `close' together tend to occur with $Y$
observations which are `close' together.


This paper is organized as follows. In Section~\ref{sec2}, the
properties of the kernel function $h_F$ are investigated in detail.
Some general properties are given, including conditions for its
existence and a proof that it is positive, and a large part of the
section is devoted to the spectral decomposition of $h_F$. We show
that if $h_F$ is square integrable, it has a mean square convergent
spectral decomposition in terms of the eigenvalues and vectors of
$h_F$. For discrete $F$, a set of difference equations is given
which has this eigensystem as its solution, and for continuous $F$
an analogous differential equation is given. The numerical solution
of these equations is treated in some detail. Closed form solutions
are only available in some special cases, for example, if $F$
belongs to the uniform distribution on $[0,1]$, the eigenfunctions
of $h_F$ are the Fourier cosine functions.


The results of Section~\ref{sec2} are used in Section~\ref{sec3} to
derive properties of $\rho^*$. We demonstrate the aforementioned
result that $0\le \rho^*(X,Y)\le 1$, such that $\rho^*(X,Y)=0$ iff
$X$ and $Y$ are independent, and $\rho^*(X,Y)= 1$ iff $X$ and $Y$
are linearly related. Furthermore, a decomposition of $\rho^*$ is
given in terms of a sum of squared correlations between marginal
eigenfunctions of $h_{F_1}$ and $h_{F_2}$ weighted with the product
of the corresponding marginal eigenvalues. We give a
parameterization of the likelihood in terms of component
correlations of $\rho^*$ and the marginal eigenfunctions, somewhat
analogous to the well-known canonical correlation decomposition.
Fr{\'e}chet bounds for the component correlations are discussed,
which gives some insight into the possible structure of the
dependence between two random variables. Finally in this section,
component correlations for the normal distribution are discussed as
an illustration.


In Section~\ref{sec4}, we derive sample and unbiased estimators of
$\kappa(X,Y)$ and related estimates of $\rho^*(X,Y)$, which can be
calculated in time $O(n^2)$. The asymptotic distributions of the
estimators under independence, which is a mixture of chi-squares, is
derived. Finally, small sample permutation tests and Bonferroni
corrections for testing the significance of component correlations
are discussed.


Section~\ref{sec5} concerns grade versions of $\kappa$ and $\rho^*$,
copulas, and rank tests. Rank statistics, obtained from the grade
versions of $\kappa$ and $\rho^*$, are discussed. It is shown that
the two sample Cram{\'e}r-von Mises statistic is obtained as a
special case, as well as a new generalization to the case of $K$
ordered samples. Furthermore, it is shown that $\rho^*$ is a
weighted mean of phi-square coefficients obtained from collapsing
the distribution onto a $2\times 2$ table with respect to cut-points
$(x,y)$.

In Section~\ref{sec6} we propose a methodology for gaining an
understanding of the association between two variables from a data
set. The methodology is based on combining hypothesis tests with
visual tools for displaying how much individual observations
contribute to the association.


Although many of the results of the present paper are new, we have,
naturally, also borrowed much from the literature, particularly
concerning the eigensystems and orthogonal decompositions. Some
important references here are \citeA{ad52,dk72}, \citeA{dwv73} and
\citeA{dewet87}, among others. However, the focus of much of the
literature we refer to is on studying power of hypothesis tests. The
aim of this paper, on the other hand, is on providing a meaningful
coefficient for describing association, which we hope leads to a
useful methodology for gaining an understanding of the association
between two variables, and, along the way, to tests with high power
against salient alternatives, the salience of the alternatives being
determined by the size of $\rho^*$.





%Finally, in Section~\ref{sec7}, a discussion is given.

%Several of the results in this paper also apply to multivariate
%$X$ and/or $Y$. However, those cases are not investigated here.


Throughout this paper, we use the following conventions and
assumptions. We assume that $(X,Y)$, $(X_1,Y_1)$ and $(X_2,Y_2)$ are
iid with marginal distribution functions $F_1$ and $F_2$,
respectively, and joint distribution function $F_{12}$. We impose no
restrictions on the distributions, i.e., they may be continuous,
discrete, or mixed continuous-discrete. For simplification of some
of the derivations, we define distribution functions in the
following slightly non-standard way:
\begin{eqnarray*}
   F_1(x) = P(X<x) + \half P(X=x)\\
   F_2(y) = P(Y<y) + \half P(Y=y)
\end{eqnarray*}
and
\begin{eqnarray*}
   \lefteqn{F_{12}(x,y) = P(X<x,Y<y) + \half P(X<x,Y=y) + }\\
   && \half P(X=x,Y<y) + \fourth P(X=x,Y=y)
\end{eqnarray*}











\section{Properties of the kernel function $h_F$}\label{sec2}

In this section a detailed description is given of the kernel
function $h_F$ defined by~(\ref{hdef1}). Section~\ref{basichf}
concerns existence, continuity, positivity, square integrability,
existence of the trace and the shape of the graph of $h_F$. Methods
for verifying whether several of these properties hold are given. In
Section~\ref{spect}, under the assumption of square integrability of
$h_F$, its spectral decomposition is given, and some properties of
the associated eigenvalues and functions are derived. The
eigensystem is the solution to an integral equation which may be
difficult to solve. The problem is reformulated in terms of
difference equations for the discrete case in Section~\ref{obt eig
1} and in terms of a differential equation for the continuous case
in Section~\ref{obt eig 2}. Both rewrites appear much easier to
solve than the integral equation. In Section~\ref{contapprox},
efficient numerical approximation of the eigensystem of $h_F$ for
continuous $F$ is discussed. For several well-known distributions,
including the uniform and the normal, closed form solutions or
numerical approximations of (parts of) the eigensystem are given. A
new distribution $F$ is introduced which has the seemingly rare
property that $h_F$ is square integrable but has infinite trace. In
Section~\ref{sec ad} the relation between $h_F$ and a kernel
introduced by \citeA{ad52} is given. We are not aware of the kernel
$h_F$, depending on $F$, having been described previously.



\subsection{Key properties of $h_F$}\label{basichf}

%This section concerns the following properties of $h_F$: existence,
%continuity, positive definiteness, square integrability, existence
%of the trace and shape of its graph.

The kernel $h_F$ {\em exists} if $h_F(z_1,z_2)$ is finite for some
$(z_1,z_2)\in\bR^2$. The kernel $h_F$ is {\em positive} if
\begin{eqnarray}
   Eg(Z_1)g(Z_2)h_F(Z_1,Z_2) \ge 0 \label{pd1}
\end{eqnarray}
for every function $g:\bR\rightarrow\bR$ for which the expectation
exists. The kernel $h_F$ is {\em square integrable} if
\begin{eqnarray}
   Eh_F(Z_1,Z_2)^2 = \int h_F(z_1,z_2)^2dF(z_1)dF(z_2)
   \label{hfsq}
\end{eqnarray}
is finite. The kernel $h_F$ is {\em trace class} if its trace
\begin{eqnarray*}
   \tr(h_F) = Eh_F(Z,Z) = \int h_F(z,z)dF(z)
\end{eqnarray*}
is finite. In several lemmas below we give some relatively easily
verifiable conditions for checking whether these properties hold for
$h_F$. The final Lemma~\ref{lem8} concerns the shape of the graph of
$h_F$.


The next lemma may simplify verification of the existence of $h_F$,
and asserts continuity and positivity as well as giving another
integral representation of $h_F$. First we need the following
notation:
\begin{eqnarray*}
   \gamma(x,y) =
   \left\{
   \begin{array}{cc}
      0     & x>y \\
      \half & x=y \\
      1     & x<y
   \end{array}
   \right.
\end{eqnarray*}
Note that
\begin{eqnarray}
    F(z)   &=& E\gamma(Z,z) \label{gamma1}\\
    1-F(z) &=& E\gamma(z,Z) \label{gamma2}
\end{eqnarray}
\begin{lemma}\label{lem1}
If $h_F$ exists it exists on $\bR^2$. It is then continuous and
positive, with equality in~(\ref{pd1}) only for the constant
function, and has the representation
\begin{eqnarray*}
   h_F(z_1,z_2) =
   \int_{-\infty}^{\infty}[\gamma(z_1, w)-F(w)][\gamma(z_2, w)-F(w)]dw
   \hspace{6mm}\forall z_1,z_2
\end{eqnarray*}
\end{lemma}
\proof{Proof} We first show continuity of $h_F$ on its domain. Let
$\delta>0$. Then if $|z_1-z_1'|<\delta$ and $|z_2-z_2'|<\delta$,
\begin{eqnarray*}
   \lefteqn{|h_F(z_1',z_2') - h_F(z_1,z_2)|}\\
   &=& \half |E\left[ (|z_1'-z_2'|-|z_1-z_2|)-(|z_1'-Z_2|-|z_1-Z_2|)-(|z_2'-Z_1|-|z_2-Z_1|) \right]| \\
   &\le& \half E\left[ 2\delta +\delta + \delta \right] \\
   &=& 2\delta
\end{eqnarray*}
Hence, $h_F$ is continuous and bounded on any finite domain.
Therefore, if $h_F$ exists in one point it exists on $\bR^2$.


We next derive the integral representation of $h_F$. We have
\begin{eqnarray}
   |z_1-z_2| = \int_{-\infty}^{\infty}[\gamma(z_1,w)-\gamma(z_2,w)]^2dw
   \label{difint}
\end{eqnarray}
and with $z_{i:4}$ the $i$th largest number in the set
$\{z_1,z_2,z_3,z_4\}$, we have
\begin{eqnarray}
   \lefteqn{
   \int_{-\infty}^{\infty}
   \left|
   [\gamma(z_1, w)-\gamma(z_3, w)]
   [\gamma(z_2, w)-\gamma(z_4, w)]
   \right|
   dw =}\nonumber\hspace*{30mm}\\
   &&
   \left\{
   \begin{array}{ll}
      0                & \mbox{if $z_1,z_3<z_2,z_4$ or $z_1,z_3>z_2,z_4$}  \\
      z_{3:4}-z_{2:4}  & \mbox{otherwise}
   \end{array}
   \right.
   \label{fubini1}
\end{eqnarray}
Now we can derive the desired result first using~(\ref{difint}),
then applying Fubini's theorem which is justified
because~(\ref{fubini1}) is finite, and finally using~(\ref{gamma1}):
\begin{eqnarray}
   h_F(z_1,z_2)
   &=& -\half E\int_{-\infty}^{\infty}
   \left(
   [\gamma(z_1, w)-\gamma(z_2, w)]^2
   -
   [\gamma(z_1, w)-\gamma(Z_2, w)]^2-\right.
   \nonumber\\
   &&
   \left.
   [\gamma(Z_1, w)-\gamma(z_2, w)]^2
   +
   [\gamma(Z_1, w)-\gamma(Z_2, w)]^2
   \right)
   dw
   \nonumber\\
   &=&
   E\int_{-\infty}^{\infty}
   [\gamma(z_1, w)-\gamma(Z_1, w)]
   [\gamma(z_2, w)-\gamma(Z_2, w)]
   dw
   \nonumber\\
   &=&
   \int_{-\infty}^{\infty}[\gamma(z_1, w)-F(w)][\gamma(z_2,
   w)-F(w)]dw\label{intrep}
\end{eqnarray}


Finally, we show positivity of $h_F$. Let $g$ be nonconstant. Then
using~(\ref{intrep}) and Fubini's theorem,
\begin{eqnarray*}
   \lefteqn{Eg(Z_1)g(Z_2)h_F(Z_1,Z_2)}\\
   &=&
   E\int_{-\infty}^{\infty}
   g(Z_1)[\gamma(Z_1, w)-F(w)]
   g(Z_2)[\gamma(Z_2, w)-F(w)]
   dw\\
   &=&
   \int_{-\infty}^{\infty}
   \left(
   Eg(Z)[\gamma(Z, w)-F(w)]
   \right)^2
   dw > 0
\end{eqnarray*}
Hence, $h_F$ is positive. If $g$ is constant it is easily verified
that the expression is zero.
 \hspace*{\fill}$\Box$\endproof


\noindent Note that from the lemma, it follows that for checking
existence of $h_F$, it suffices to check existence of $h_F(0,0)$.
Now $h_F(0,0)$ has the following convenient representations
\begin{eqnarray}
   h_F(0,0)
   &=&
   \half E
   \left(
      |Z_1|+|Z_2|-|Z_1-Z_2|
   \right) \nonumber\\
   &=&
   \int_{-\infty}^0F(z)^2dz + \int_0^{\infty}[1-F(z)]^2dz
   \label{hf00}
\end{eqnarray}
These representations are immediately verified from~(\ref{hdef1})
and from the representation of $h_F$ given in Lemma~\ref{lem1}.


An example of a random variable for which $h_F$ does not exist is
$Z=V^2$, where $V$ has a Cauchy distribution. This can be verified
by checking that~(\ref{hf00}) does not converge.


By giving some alternative representations of~(\ref{hfsq}), the next
lemma may be helpful in the verification of square integrability of
$h_F$:
\begin{lemma}\label{lem3}
We have:
\begin{eqnarray}
   Eh_F(Z_1,Z_2)^2
   =
   \frac16E(Z_{2:4}-Z_{3:4})^2
   =
   2\int_{z_1<z_2}F(z_1)^2[1-F(z_2)^2]dz_1dz_2
   \label{sqintrepr}
\end{eqnarray}
\end{lemma}
\proof{Proof} To prove the first equality, write $a_{ij}=|Z_i-Z_j|$.
Then, since for example $Ea_{12}=Ea_{13}$ and
$Ea_{34}a_{12}=Ea_{34}a_{15}$, we obtain
\begin{eqnarray*}
   Eh_F(Z_1,Z_2)^2
   &=&   \frac14 E
   \left(a_{12}-a_{13}-a_{24}+a_{34}\right)
   \left(a_{12}-a_{15}-a_{26}+a_{56}\right)
   \\
   &=&   \frac14 E
   \left(
   a_{12}^2
   -a_{12}a_{15}
   -a_{12}a_{26}
   +a_{12}a_{56}\right.
   \\ &&
   -a_{13}a_{12}
   +a_{13}a_{15}
   +a_{13}a_{26}
   -a_{13}a_{56}
   \\
   &&
   -a_{24}a_{12}
   +a_{24}a_{15}
   +a_{24}a_{26}
   -a_{24}a_{56}
   \\ &&
   \left.
   +a_{34}a_{12}
   -a_{34}a_{15}
   -a_{34}a_{26}
   +a_{34}a_{56}
   \right)
   \\
   &=&   \frac14 E
   \left(
   a_{12}^2
   -a_{12}a_{13}
   -a_{12}a_{24}
   +a_{12}a_{34}\right)
   \\
   &=& \frac{1}{16} E
   \left(
   a_{12}-a_{13}-a_{24}+a_{34}
   \right)^2
\end{eqnarray*}
It may now be verified that $a_{12}-a_{13}-a_{24}+a_{34}$ equals 0
if $Z_1,Z_4\le Z_2,Z_3$ or $Z_1,Z_4\ge Z_2,Z_3$, and equals $\pm
2(Z_{2:4}-Z_{3:4})$ otherwise. Hence, and because
$P(a_{12}-a_{13}-a_{24}+a_{34}\ne 0)=2/3$, we have
\begin{eqnarray*}
   Eh_F(Z_1,Z_2)^2
   &=&
   \frac{1}{16}
   \left(
   a_{12}-a_{13}-a_{24}+a_{34}
   \right)^2\\
   &=&
   \frac{1}{16}\times
   \frac{2}{3}\times
   4\,E(Z_{2:4}-Z_{3:4})^2
   = \frac16E(Z_{2:4}-Z_{3:4})^2
\end{eqnarray*}
which is the first part of the lemma.


To prove the second equality, first note that
\begin{eqnarray}
   E\gamma(Z,z_1)\gamma(Z,z_2)=\min\{F(z_1),F(z_2)\}
   \label{eqmin}
\end{eqnarray}
We now have using Lemma~\ref{lem1}, Fubini's theorem (for
justification see proof of Lemma~\ref{lem1}) and~(\ref{eqmin}),
\begin{eqnarray*}
   Eh_F(Z_1,Z_2)^2
   &=& E\int[\gamma(Z_1,z_1)-F(z_1)][\gamma(Z_2,z_1)-F(z_1)]dz_1 \times\nonumber\\
   &&   \int[\gamma(Z_1,z_2)-F(z_2)][\gamma(Z_2,z_2)-F(z_2)]dz_2\nonumber\\
   &=& \int E[\gamma(Z_1,z_1)-F(z_1)][\gamma(Z_1,z_2)-F(z_2)]\times \nonumber\\
   &&E[\gamma(Z_2,z_1)-F(z_1)][\gamma(Z_2,z_2)-F(z_2)]dz_1dz_2 \nonumber\\
   &=& \int [\min\{F(z_1),F(z_2)\}-F(z_1)F(z_2)]^2dz_1dz_2\nonumber\\
   &=& 2\int_{z_1<z_2} F(z_1)^2[1-F(z_2)]^2dz_1dz_2\nonumber
\end{eqnarray*}
 \hspace*{\fill}$\Box$\endproof

\noindent An example of a distribution for which $h_F$ exists but is
not square integrable is the Cauchy distribution. In particular,
$h_F(0,0)=2\pi^{-1}\log 2$, so $h_F$ exists. In this case,
nonexistence of $Eh_F(Z_1,Z_2)^2$ can most easily be verified using
the right hand side of~(\ref{sqintrepr}) and a computer algebra
package such as provided in Mathematica.

The next lemma may be helpful in verifying whether or not $h_F$ is
trace class:
\begin{lemma}\label{lemtr}
We have
\begin{eqnarray*}
   \tr(h_F) = \half E|Z_1-Z_2| = \int F(z)[1-F(z)]dz
\end{eqnarray*}
which is finite iff $Z$ has finite mean.
%If $h_F$ has finite trace,it is square integrable.
\end{lemma}
\proof{Proof} The first equality follows directly
from~(\ref{hdef1}), the second is well-known and can be found using
similar techniques as in the proof of the second equality in
Lemma~\ref{lem3}. As is well-known, integration by parts leads to
the representation of the mean as
\begin{eqnarray*}
   EZ = \int_0^{\infty}[1-F(z)]dz - \int_{-\infty}^0F(z)dz
\end{eqnarray*}
so the mean exists iff the terms on the right hand side exists. Now
since
\begin{eqnarray*}
   F(0)\int_0^{\infty}[1-F(z)]dz \le \int_0^{\infty}F(z)[1-F(z)]dz \le \int_0^{\infty}[1-F(z)]dz
\end{eqnarray*}
and
\begin{eqnarray*}
   [1-F(0)]\int_{-\infty}^0[1-F(z)]dz \le \int_{-\infty}^0F(z)[1-F(z)]dz \le \int_{-\infty}^0F(z)dz
\end{eqnarray*}
it follows that
\begin{eqnarray*}
   \tr(h_F) = \int_0^{\infty}F(z)[1-F(z)]dz + \int_{-\infty}^0F(z)[1-F(z)]dz
\end{eqnarray*}
exists iff $EZ$ exists.
 \hspace*{\fill}$\Box$\endproof

\noindent The quantity $E|Z_1-Z_2|$ is also called {\em Gini's mean
difference}. Note that, by Lemmas~\ref{lem3} and~\ref{lemtr}, both
$Eh_F(Z_1,Z_2)^2$ and $\tr(h_F)$ can be used as measures of
dispersion.


An example of a distribution function $F$ for which $h_F$ is square
integrable but not trace class is given in Example~\ref{ex1} in the
next subsection.



We conclude this section with a lemma concerning the shape of the
graph of $h_F$. In Figure~\ref{hnorm} a representation of the graph
of $h_F$ with $F$ the CDF of the normal distribution is given. The
statements of Lemma~\ref{lem8} can be verified in the plot. Then:
\begin{lemma}\label{lem8}
Suppose $F$ is such that $h_F$ exists. Then:
\begin{enumerate}

\item For given $z_1$, $h_F(z_1,z_2)$ is strictly decreasing in $z_2$ on the domain $\{z:z\ge z_1,F(z)<1\}$ and strictly increasing in $z_2$
on the domain $\{z:z\le z_1,F(z)>0\}$.

\item $h_F(z,z)$ is strictly increasing in $z$ on the domain $\{z:F(z)>\half\}$ and strictly decreasing in $z$ on the domain $\{z:F(z)<\half\}$.

\end{enumerate}
\end{lemma}
\proof{Proof} Part 1: With $z_1,z_2,z_3$ such that $z_1<z_2\le z_3$
and $F(z_3)<1$, we obtain using Lemma~\ref{lem1} that
\begin{eqnarray*}
   h_F(z_1,z_3) - h_F(z_1,z_2)
   &=& \int[\gamma(z_1,w)-F(w)][\gamma(z_3,w)-\gamma(z_2,w)]dw \\
   &=& -\int_{z_2}^{z_3}[1-F(w)]dw\\
   &<& 0
\end{eqnarray*}
where the strict inequality holds because $F(z_3)<1$. This proves
the strict decreasingness part, the strict increasingness is proven
by appropriately reversing signs in the above.


Part 2: For $z<z'$ we have
\begin{eqnarray*}
   h_F(z,z) - h_F(z',z')
   &=&
   \int_{-\infty}^{\infty}[\gamma(z,w)-F(w)]^2-[\gamma(z',w)-F(w)]^2dw\\
   &=&
   \int_{z}^{z'}[1-2F(w)]dw
\end{eqnarray*}
which is positive if $F(z')\ge F(z)>\half$ and negative if $F(z)\le
F(z')<\half$, proving the monotonicity relations.
 \hspace*{\fill}$\Box$\endproof


\begin{figure}
\begin{center}
\includegraphics[width=.6\linewidth,clip]{hnormplot.eps}
\end{center}
%\parbox{0.95\textwidth}
{\caption{Graph of $h_F$ with $F$ the CDF for the standard normal
distribution{\label{hnorm}}}}
\end{figure}










\subsection{Spectral decomposition of $h_F$}\label{spect}

A sequence of random variables $Z^{(N)}$ is said to converge to $Z$
in mean square if
\begin{eqnarray*}
   E\left(Z - Z^{(N)}\right)^2 \rightarrow 0
   \mbox{ as } N\rightarrow\infty
\end{eqnarray*}
For a distribution function $F$, we define
\begin{eqnarray*}
   L_2(F) = \left\{g:\bR\rightarrow\bR\left|\int g(z)^2dF(x)<\infty\right.\right\}
\end{eqnarray*}
as the set of square integrable functions with respect to $F$. A set
of functions $\{g_k\}$ is said to be {\em orthonormal} if
\begin{eqnarray}
   Eg_k(Z)g_l(Z) = \delta_{kl} \label{g1}
\end{eqnarray}
(with $\delta$ the Kronecker delta) and
\begin{eqnarray*}
   Eg_k(Z)^2=1
\end{eqnarray*}
An orthonormal system $\{g_k\}$ is {\em complete} if for any
function $g\in L_2(F)$ there exist numbers $\{\alpha_k\}$ such that
\begin{eqnarray*}
   g(Z) = \sum_{k=1}^\infty\alpha_k g_k(Z)
\end{eqnarray*}
with convergence in mean square. Then the system $\{g_k\}$ is also
called a {\em basis} of $L_2(F)$. A number $\lambda$ is called an
eigenvalue of $h_F$ with corresponding eigenvector $g$ if
$Eg(Z)^2=1$ and
\begin{eqnarray}
   \lambda g(z) = Eh_F(z,Z)g(Z)  \label{g2}
\end{eqnarray}


Then we have:
\begin{theorem}\label{th2a}
Suppose $h_F$ is square integrable. Then there exists a complete
system of orthonormal functions $\{g_k\}$ of $L_2(F)$ consisting of
eigenfunctions of $h_F$, the corresponding eigenvalues
$\{\lambda_k\}$ being nonnegative. Each $g_k$ is continuous and
satisfies
\begin{eqnarray}
   Eg_k(Z) = 0 \label{degen}
\end{eqnarray}
if $\lambda_k\ne 0$. Furthermore, $h_F$ has the spectral
decomposition
\begin{eqnarray}
   h_F(Z_1,Z_2) = \sum_{k=1}^{\infty}\lambda_kg_k(Z_1)g_k(Z_2)
   \label{spectr1}
\end{eqnarray}
where convergence is in mean square.
\end{theorem}
\proof{Proof} A continuous square integrable kernel on a
sigma-finite measure space is a Hilbert-Schmidt kernel which by a
generalization of Mercer's theorem has the desired spectral
decomposition \cite{zaanen60}. It is easy to verify that $(\bR,B,F)$
is a $\sigma$-finite measure space, so the operator mapping the
function $g$ to $\int h_F(x,y)g(y)dF(y)$ is a Hilbert-Schmidt
operator. Nonnegativity of the eigenvalues follows from positivity
of $h_F$ (see Lemma~\ref{lem1}). Furthermore, from~(\ref{g2}),
$\lambda Eg_k(Z)=0$ so $Eg_k(Z)=0$ for nonzero $\lambda_k$.
 \hspace*{\fill}$\Box$\endproof


Note that if $g$ is a solution to~(\ref{g2}), then so is $-g$. To
identify the solutions, we proceed as follows. A function $g$ is
{\em initially positive} if there exists a $z$ such that $g(z)>0$
and $g(z')<0$ for all $z'<z$. Without loss of generality, we may
assume that the $g_k$ are initially positive. We also assume the
eigenvalues are ordered: $\lambda_1>\lambda_2>\ldots$. An
interesting property of eigenfunctions of homogeneous positive
Fredholm integral equations of the second kind is that they are
oscillating, in the sense that for every $k$, $g_k$ has $k$ distinct
zeroes and no more than $k-1$ local extremes.


Some further results are as follows:
\begin{lemma}\label{lem6}
For square integrable $h_F$:
\begin{enumerate}\label{lemtr2}

\item if finite, $\tr(h_F)=\sum_{k=0}^{\infty}\lambda_k$

\item $Eh_F(Z_1,Z_2)^2=\sum_{k=0}^{\infty}\lambda_k^2$.

\end{enumerate}
\end{lemma}
\proof{Proof} From the spectral decomposition~(\ref{spectr1}) and
Fubini's theorem,
\begin{eqnarray*}
   Eh_F(Z,Z) = E\sum\lambda_kg_k(Z)^2 = \sum\lambda_k
\end{eqnarray*}
For the second part, the desired result is obtained by Parseval's
theorem.
 \hspace*{\fill}$\Box$\endproof


As an example, we consider the dichotomous case which is the
simplest possible and has closed form solutions:
\begin{example}
Consider the dichotomous case that $P(Z=0)=1-P(Z=1)=p$. Then
\begin{eqnarray*}
   h_F(0,0) &=& (1-p)^2 \\
   h_F(0,1) &=& -p(1-p) \\
   h_F(1,0) &=& -p(1-p) \\
   h_F(1,1) &=& p^2
\end{eqnarray*}
The eigenvalues are $(\lambda_0,\lambda_1) = (0,p(1-p))$ and the
eigenfunctions are $g_0(z)=1$ and
\begin{eqnarray*}
   g_1(z) = \frac{1}{\sqrt{p(1-p)}}[z-(1-p)]
%   \left(g_1(0),g_1(1)\right) = \left(-\sqrt{\frac{1-p}{p}},\sqrt{\frac{p}{1-p}}\right)
\end{eqnarray*}
which has mean zero and variance one. Now $h_F(z_1,z_2) =
\lambda_1g_1(z_1)g_1(z_2)$.
\end{example}
Solutions to~(\ref{g2}) for various other $F$ are given in
Table~\ref{eigtable}, see Section~\ref{obt eig 2} for an
explanation.


Equation~(\ref{g2}) is a homogeneous Fredholm integral equation of
the second kind based on the degenerate kernel $h_F$ (see, for
example, \citeNP{tricomi85}). In general these equations are
difficult to solve. In Sections~\ref{obt eig 1} and~\ref{obt eig
2},~(\ref{g2}) is reduced to a simpler problem for the discrete and
continuous case, respectively. We conclude this section with an
alternative formula to~(\ref{g2}) for finding the eigensystem which
we use in the next two subsections. With
\begin{eqnarray*}
   G(z) = \int_{-\infty}^{z}g(y)dF(y),
\end{eqnarray*}
we have:
\begin{lemma}\label{lem5}
The nonzero eigenvalues and eigenvectors of $h_F$ are solutions to
the equation
\begin{eqnarray*}
   \lambda
   \left[
   g(z) - g(z')
   \right]
   +
   \int_{z}^{z'}G(w)dw = 0
\end{eqnarray*}
for all $z<z'$.
%subject to the side condition $G_k(z)\rightarrow 0$ as $|z|\rightarrow\infty$.
\end{lemma}
\proof{Proof} Substitution of the expression for $h_F$ of
Lemma~\ref{lem1} into~(\ref{g2}) yields
\begin{eqnarray*}
   \lambda g (z)
   &=&
   \int_{-\infty}^{\infty}
   \int_{-\infty}^{\infty}
   [\gamma(z, w)-F(w)][\gamma(z_2, w)-F(w)]g (z_2)dF(z_2)dw\\
   &=&
   \int_{-\infty}^{\infty}
   [\gamma(z, w)-F(w)]G(w)dw\\
   &=&
   -\int_{-\infty}^{z}F(w)G (w)dw +
   \int_{z}^{\infty}[1-F(w)]G (w)dw
\end{eqnarray*}
The Lemma follows from this and since by~(\ref{degen})
$G(z)\rightarrow 0$ as $|z|\rightarrow\infty$.
 \hspace*{\fill}$\Box$\endproof

\noindent We have the following interesting implication:
\begin{corollary}\label{lindens}
If $\langle a,b\rangle$ is an interval with zero probability mass,
i.e., $F(a)=F(b)$, then a solution $g(z)$ to~(\ref{g2}) is linear on
$\langle a,b\rangle$.
\end{corollary}
\proof{Proof} If $F(a)=F(b)$ then from its definition it follows
that $G$ is constant on $\langle a,b\rangle$. Therefore, for any
$z,z'\in\langle a,b\rangle$, we obtain from Lemma~\ref{lem5} that
\begin{eqnarray*}
   \lambda
   \left[
   g(z) - g(z')
   \right]
   = c(z-z')
\end{eqnarray*}
for some constant $c$. Hence, $g$ is linear on $\langle a,b\rangle$.
 \hspace*{\fill}$\Box$\endproof

\noindent Note that for discrete distributions, it follows that the
eigenfunctions $g$ are piecewise linear.





%Plots of some eigenfunctions?????





\subsection{Obtaining the eigensystem in the discrete case}\label{obt eig 1}

We consider the case that $Z$ is a.s.\ in a finite set, say
$P(Z\in\{z_1,z_2,\ldots,z_K\})=1$. We use the shorthand
$P(Z=z_i)=p_i$ and assume without loss of generality that $p_i>0$
and $z_i<z_{i+1}$ for all $i$. Then we obtain:
\begin{lemma}\label{cor2}
With $c_i=(z_i-z_{i-1})^{-1}$ the nonzero eigenvalues and
eigenvectors of $h_F$ are solutions to the equations
\begin{eqnarray*}
   &&p_1g (z_1) = \lambda c_2[g (z_1)-g (z_2)] \\
   &&p_Kg (z_K) = \lambda c_K[g (z_{K-1})-g (z_{K})]\\
   && p_ig (z_i) =
   \lambda
   \left[c_ig (z_{i-1})-(c_{i}+c_{i+1})g (z_{i})+c_{i+1}g (z_{i+1})
   \right]
   \hspace{6mm}2\le i\le K-1
\end{eqnarray*}
\end{lemma}
\proof{Proof} First note that from the definition
%\begin{eqnarray*}
%   G (z_i) = \sum_{j=1}^{i-1}p_jg (z_j)+\half p_ig (z_i)
%\end{eqnarray*}
for $z\in\langle z_i,z_{i+1}\rangle$,
\begin{eqnarray*}
   G (z) = \sum_{j=1}^{i}p_jg (z_j)
\end{eqnarray*}
It follows that
\begin{eqnarray}
   \int_{z_i}^{z_{i+1}}G (w)dw =
   (z_{i+1}-z_i)\sum_{j=1}^{i}p_jg (z_j)
   \label{int2}
\end{eqnarray}
The first two displayed equations of the lemma now follow
from~(\ref{int2}), Lemma~\ref{lem5} and from $\sum p_ig (z_i)=0$.
From~(\ref{int2}) we further obtain
\begin{eqnarray}
   c_{i+1}\int_{z_i}^{z_{i+1}}G (w)dw -c_{i} \int_{z_{i-1}}^{z_i}G (w)dw =
   p_ig (z_i)
   \label{dif1}
\end{eqnarray}
Again from Lemma~\ref{lem5}, we have
\begin{eqnarray*}
   \lambda
   \left[
   g (z_{i}) - g (z_{i+1})
   \right]
   +
   \int_{z_{i}}^{z_{i+1}}G (w)dw
   &=& 0
   \\
   \lambda
   \left[
   g (z_{i-1}) - g (z_{i})
   \right]
   +
   \int_{z_{i-1}}^{z_{i}}G (w)dw
   &=& 0
\end{eqnarray*}
Multiplying these equations by $c_{i+1}$ and $c_{i}$ respectively,
taking differences and the use of~\ref{dif1} yields
\begin{eqnarray*}
   p_ig_k(z_i) =
   \lambda
   \left[c_ig_k(z_{i-1})-(c_{i}+c_{i+1})g_k(z_{i})+c_{i+1}g_k(z_{i+1})
   \right]
\end{eqnarray*}
which completes the proof.
 \hspace*{\fill}$\Box$\endproof

More details on difference equations of the form given in
Lemma~\ref{cor2} can be found in \citeA{agarwal92}. Lemma~\ref{cor2}
allows fast and memory efficient computation of the eigenvalues and
vectors. In matrix notation, we must solve the generalized
eigenvalue problem
\begin{eqnarray*}
   D_pg = \lambda C g
\end{eqnarray*}
where $g$ is the eigenvector with corresponding eigenvalue
$\lambda$, $D_p$ is a diagonal matrix with the coordinates of $p$
on the main diagonal, and
\begin{eqnarray*}
   C =
   \left(
   \begin{array}{ccccccc}
      c_2 & c_2        & 0          & 0 \\
      c_2 & -(c_2+c_3) & c_3        & 0   &  & \ldots\\
      0   & c_3        & -(c_3+c_4) & c_4 & \\
      0   & 0          & c_4        & -(c_4+c_5) \\
          &            &            &        & \ddots \\
          & \vdots     &            &        & & -(c_{K-1}+c_K) & c_K      \\
          &            &            &        & & c_K         & c_K      \\
   \end{array}
   \right)
\end{eqnarray*}
This method can also be used for fast approximation of continuous
systems (see Section~\ref{contapprox}). The exact solution for
continuous systems is given in the next subsection.



\subsection{Obtaining the eigensystem in the continuous case}
\label{obt eig 2}

If $F$ is differentiable, the problem of finding the eigenvalues
and vectors can be reduced to a differential equation which is
sometimes easier to solve:
\begin{lemma}
Suppose $f$ is the derivative of $F$. Then the nonzero eigenvalues
and eigenvectors of $h_F$ are solutions to the equation
\begin{eqnarray}
   \lambda g''(z) + f(z)g (z) = 0
   \label{de1}
\end{eqnarray}
subject to the side condition
\begin{eqnarray*}
   g'(z) \rightarrow 0 \hspace{4mm}\mbox{ as }|z|\rightarrow\infty
\end{eqnarray*}
\end{lemma}
\proof{Proof} Letting $z'\rightarrow z$ in Lemma~\ref{lem5} we
obtain
\begin{eqnarray*}
   \lambda g'(z) + G (z) = 0
\end{eqnarray*}
Differentiating both sides with respect to $z$ yields~(\ref{de1}).
Finally, since by~(\ref{degen}) $G (z)\rightarrow 0$ as
$|z|\rightarrow\infty$, the side condition follows.
 \hspace*{\fill}$\Box$\endproof

Equation~(\ref{de1}) leads to an interesting observation on the
behavior of the eigenfunctions $g_k$: the second derivative
$g_k''(x)$ is proportional to the local density $f(x)$ times
$g_k(x)$. As mentioned in Section~\ref{spect} and as follows from
the theory of Sturm-Liouville differential equations, the
eigenfunctions oscillate. Now if the local density is low, the
oscillatory behavior will be slower. In particular, on intervals
where the local density is zero the second derivative of an
eigenfunction is zero so the eigenfunction is linear.  (NB: this
also holds for non-continuous distributions, see
Corollary~\ref{lindens}).

If $F$ is both differentiable and invertible, we can reformulate
differential equation~(\ref{de1}) in the standard Sturm-Liouville
form:
\begin{lemma}
Suppose $F$ is invertible. Let $q (u)=g [F^{-1}(u)]$ and suppose $q
$ is twice differentiable. Then the eigenvalues and eigenvectors are
solutions to the equation
\begin{eqnarray}
   \frac{d}{du}f[F^{-1}(u)]q '(u) + \lambda ^{-1}q (u) = 0
   \label{de2}
\end{eqnarray}
subject to the side condition
\begin{eqnarray*}
   f[F^{-1}(u)]q '(u) \rightarrow 0 \mbox{ as } u\downarrow 0 \mbox{ or } u\uparrow 1
\end{eqnarray*}
\end{lemma}
\proof{Proof} Note that
\begin{eqnarray*}
   \frac{d}{du}F^{-1}(u) =
   \frac{1}{f[F^{-1}(u)]}
\end{eqnarray*}
so that
\begin{eqnarray*}
   f[F^{-1}(u)]q '(u) = f[F^{-1}(u)]\frac{d}{du}g [F^{-1}(u)] =
   g '[F^{-1}(u)]
\end{eqnarray*}
From this the side condition follows. Substituting $z=F^{-1}(u)$
into the left hand side of~(\ref{de1}) yields
\begin{eqnarray*}
   \lefteqn{\lambda g ''[F^{-1}(u)] -  f(F^{-1}(u))g [F^{-1}(u)]}\\
   &=&
   \lambda \frac{du}{dF^{-1}(u)}\,\frac{d}{du}\,g '[F^{-1}(u)] -
   f(F^{-1}(u))g [F^{-1}(u)]\\
   &=&
   \lambda f[F^{-1}(u)]\,\frac{d}{du}\,f[F^{-1}(u)]q '(u) -  f[F^{-1}(u)]q (u)
\end{eqnarray*}
Hence dividing both sides of~(\ref{de1}) by $\lambda f[F^{-1}(u)]$
yields the desired result.
 \hspace*{\fill}$\Box$\endproof
Note that for the $q_k$ the orthonormality condition~(\ref{g1})
reduces to
\begin{eqnarray*}
   \int q_k(u)q_l(u)du = \delta_{kl}
\end{eqnarray*}


In general, the differential equations~(\ref{de1}) and~(\ref{de2})
do not have closed form solutions (i.e., solutions in terms of
well-known functions). We obtained the solutions for the uniform,
the logistic and the exponential distributions which are given in
Table~\ref{eigtable}, where $P_k$ is the $k$th Legendre polynomial,
$J_i$ is the $i$th Bessel function of the first kind, $\alpha_{k}$
the $k$th zero of $J_1$ and
\begin{eqnarray*}
   \beta_k=\left(J_0(\alpha_k)^2+J_1(\alpha_k)^2\right)^{-1/2}
\end{eqnarray*}
Numerical solutions were obtained for the normal, Laplace and
chi-square distribution with one degree of freedom, see the next
subsection for details on obtaining these solutions. For the normal
and Laplace distributions we obtained the exact value of
$\sum\lambda_k$ and for the normal distribution of $\sum\lambda_k^2$
using Lemmas~\ref{lem3} and~\ref{lemtr} combined with results for
order statistics by \citeA{bg59} and \citeA{govindarajulu63}
summarized in \citeA{jkb94} and using Mathematica.

We also obtained a closed form expression for the eigensystem for
the distribution introduced in the next example. It is an example of
a distribution $F$ for which $h_F$ is square integrable but not
trace class, i.e., by Lemma~\ref{lemtr2}, $\sum\lambda_k^2$ is
finite but $\sum\lambda_k=\infty$. We are not aware of any previous
studies of this distribution.
\begin{example}\label{ex1}
With $V$ standard normally distributed, let $Z$ be defined as the
following function of $V$:
\begin{eqnarray*}
   Z = \sqrt{\frac{\pi}2}\int_0^V \exp(t^2/2)dt
\end{eqnarray*}
where the convention is used that for $a>b$,
\begin{eqnarray*}
   \int_a^b f(t)dt=-\int_b^a f(t)dt
\end{eqnarray*}
Close to zero, $Z$ has approximately a normal density, but for large
values the density is much lower, that is, $Z$ has much heavier
tails than a standard normal.

With $F$ the distribution function of $Z$, we now show that $h_F$ is
square integrable but not trace class. To show this, we shall derive
an expression for $f[F^{-1}(u)]$ to be plugged into
Equation~(\ref{de2}) so that it can be solved. The derivation
involves the so-called complex error function. The {\em error
function} is defined as
\begin{eqnarray*}
  \erf(z) = \sqrt{2\pi}\int_0^z\exp(-t^2)dt
\end{eqnarray*}
Note that the CDF of the standard normal distribution is
$\Phi(v)=(1+\erf(v/\sqrt{2}))/2$. The {\em imaginary error function}
is defined as
\begin{eqnarray*}
  \erfi(z) = -i\,\erf(i\,z)
\end{eqnarray*}
(Here $i=\sqrt{-1}$. See \citeNP{weisstein:erf-erfi}, for some of
the properties of $\erf$ and $\erfi$.) We define the inverse
$\erfi^{-1}(z)$ as the unique real $y$ satisfying $z=\erfi(y)$.

We see that
\begin{eqnarray*}
  Z = \pi\,\erfi\left(\frac{V}{\sqrt{2}}\right)
\end{eqnarray*}
Now since $V$ has a standard normal distribution we obtain for the
CDF of $Z$:
\begin{eqnarray*}
  F(z)
  = P(Z<z)
  = P\left(V < \sqrt{2}\,\erfi^{-1}(z/\pi)\right)
  = \Phi\left(\sqrt{2}\,\erfi^{-1}(z/\pi)\right)
\end{eqnarray*}
From this,
\begin{eqnarray*}
  F^{-1}(u) = \pi\,\erfi\left(\Phi^{-1}(u)/\sqrt{2}\right)
\end{eqnarray*}
Some tedious but straightforward algebra then gives
\begin{eqnarray*}
  f(F^{-1}(u)) = 1/\phi[\Phi^{-1}(u)]^{2}
\end{eqnarray*}
Plugging this expression into~(\ref{de2}) leads to the solution
\begin{eqnarray*}
   \lambda_k &=& 1/k \\
   q_k(u) &\equiv& H_k\left[\Phi^{-1}(u)\right]\sqrt{\phi\left[\Phi^{-1}(u)\right]}
\end{eqnarray*}
where $H_k$ is the $k$th Hermite polynomial. The `$\equiv$' symbol
is used to indicate that the expression for $q_k$ needs to be
suitably normalized. This solution of~(\ref{de2}) was derived by
\citeA{dwv73}, who provided a method for solving differential
equations of the form $\frac{d}{du}w(u)q_k'(u) +
\lambda_k^{-1}q_k(u)=0$ for certain types of weights $w(u)$,
including $w(u)=1/\phi[\Phi^{-1}(u)]^{2}$. Note that for $g_k$ we
obtain
\begin{eqnarray*}
   g_k(x) \equiv H_k\left[\sqrt{2}\,\erfi^{-1}(\pi\,z)\right]\sqrt{\phi\left[\sqrt{2}\,\erfi^{-1}(\pi\,z)\right]}
\end{eqnarray*}
It is well-known that here $\sum\lambda_k$ is divergent and
$\sum\lambda_k^2=\pi^2/6$.
\end{example}



\begin{table}
\begin{center}
\begin{minipage}{120mm}
\begin{tabular}{l|llll}
 Distribution      & Logistic              & Uniform             & Normal                &  \\ \hline
$1/f[F^{-1}(u)]$  & $u(1-u)$                   & $1$             & $\phi[\Phi^{-1}(u)]$  &   \\
$\lambda_k$       & $\frac{1}{k(k+1)}$  & $\frac{1}{k^2\pi^2}$   & $^{\mbox{\it\scriptsize a}}$  \\
$\sum\lambda_k$   & $1$         & $\frac{1}{6}$                  & $\frac{1}{\sqrt\pi}$            \\
$\sum\lambda_k^2$ & $\frac13(\pi^2-9)$        & $\frac{1}{90}$   & $\frac13 - \frac{\sqrt{3}-1}{\pi}$           \\
$q_k(u)$          & $\sqrt{2k+1}\,P_k(2u-1)$       & $2\cos(k\pi u)$ & $^{\mbox{\it\scriptsize a}}$               \\
$\lambda_1/\sum\lambda_k$ & 0.5000        & 0.6079               & 0.5269 \\
$\lambda_2/\sum\lambda_k$ & 0.1667        & 0.1520               & 0.1635 \\
$\lambda_3/\sum\lambda_k$ & 0.0833        & 0.0675               & 0.0795 \\
$\lambda_4/\sum\lambda_k$ & 0.0500        & 0.0380               & 0.0470 \\
\\
 Distribution      & Exponential           & Laplace              & Chi-square  & Example~\ref{ex1} \\ \hline
$1/f[F^{-1}(u)]$  & $1-u$                 & $\min(u,1-u)$        & $\frac{\phi[\Phi^{-1}(\frac{u+1}{2})]}{\Phi^{-1}\left(\frac{u+1}{2}\right)}$   & $\phi[\Phi^{-1}(u)]^2$     \\
$\lambda_k$       & $\frac{4}{\alpha_k^2}$& $^{\mbox{\it\scriptsize a}}$                 & $^{\mbox{\it\scriptsize a}}$     & $\frac{1}{k}$  \\
$\sum\lambda_k$   & $\frac12$             & $\frac34$            & 0.6360  & $\infty$ \\
$\sum\lambda_k^2$ & $\frac{1}{12}$        & 0.1458               & 0.1399  & $\frac{6}{\pi^2}$  \\
$q_k(u)$          & $\beta_kJ_0(\alpha_{k}\sqrt{1-u})$ & $^{\mbox{\it\scriptsize a}}$     & $^{\mbox{\it\scriptsize a}}$     & $H_k[\Phi^{-1}(u)]\sqrt{\phi[\Phi^{-1}(u)]}^{\mbox{\it\scriptsize b}}$ \\
$\lambda_1/\sum\lambda_k$ & 0.5453        & 0.4611               & 0.5567  & 0\\
$\lambda_2/\sum\lambda_k$ & 0.1627        & 0.1816               & 0.1615  & 0\\
$\lambda_3/\sum\lambda_k$ & 0.0774        & 0.0875               & 0.0758  & 0\\
$\lambda_4/\sum\lambda_k$ & 0.0451        & 0.0542               & 0.0438  & 0\\
\end{tabular}
\end{minipage}
\caption{Eigenvalues and eigenvectors of kernel function $h_F$ for
various $F$. \protect{\newline} $\,^{\mbox{\it\scriptsize a}}$ No
closed form expression is available. \protect{\newline}
$\,^{\mbox{\it\scriptsize b}}$ Expression not normalized.
 \label{eigtable}}
\end{center}
\end{table}


The differential equation~(\ref{de2}) with $f[F^{-1}(z)]$ replaced
by a weight function $w(z)$ is given in \citeA{ad52}, see also
\citeA{dwv73} and \citeA{dewet87}. We are not aware of
equation~(\ref{de1}) having been studied previously.







%\newcommand{\erf}{\mbox{erf}}
%\newcommand{\serf}{\mbox{\scriptsize erf}}



%\subsection{Numerical approximation and estimation of eigenvalues}




\subsection{Discrete approximation of the continuous
case}\label{contapprox}

For many continuous distribution functions $F$ the differential
equations~(\ref{de1}) and~(\ref{de2}) do not have a closed form
solution and the first $t$ eigenvalues and eigenfunctions can be
approximated by using a discrete approximation of $F$ and solving
the difference equations of Lemma~\ref{cor2}. For $i=1,\ldots,t$,
let $c_i=F^{-1}\left(\frac{i-1/2}{t}\right)$ and let $Z^{(t)}$ be a
discrete random variable with
$P(Z^{(t)}=c_i)=p_i=F(i/t)-F((i-1)/t)$. The eigenvalues and
eigenvectors of $h_{F^{(t)}}$, with $F^{(t)}$ the distribution
function of $Z^{(t)}$, can then be calculated using the method of
Section~\ref{obt eig 1}. For large $t$, this method seems to give
good approximations of the eigenvalues and eigenvectors of $h_F$. An
idea of the quality of the approximations can be gained from
Table~\ref{eigtable2}. The numerical results in Table~\ref{eigtable}
were obtained using this method. For further details on discrete
approximations of eigenvalues and vectors of kernels see
\citeA{tricomi85}.

To obtain a good approximation, $t$ should of course be chosen as
large as possible. Using Mathematica~5.2 on a Pentium~IV computer at
3.0MHz, using no special routines for calculating the eigensytem of
tridiagonal matrices, calculation of a complete solution for
$t=1000$ took 29 seconds. We expect that using software with such
special routines, it is possible to obtain solutions of the first
few eigenvalues and eigenvectors much more quickly and for much
larger $t$.



\begin{table}
\begin{center}
\begin{tabular}{l|l|l|l}
                  & True value                & Estimate ($t=101$)        & Estimate ($t=1001$)         \\
$\sum\lambda_k$   & 1                         & 0.99303                   & 0.99931           \\
$\sum\lambda_k^2$ & 0.28987                   & 0.29027                   & 0.28988           \\
$\lambda_1^*$     & 0.50000                   & 0.50370                   & 0.50035   \\
$\lambda_{10}^*$  & 9.0909$\times 10^{-3}$    & 9.3093$\times 10^{-3}$    & 9.1056$\times 10^{-3}$\\
$\lambda_{100}^*$ & 9.9010$\times 10^{-5}$    & 9.9708$\times 10^{-5}$    & 9.9145$\times 10^{-5}$     \\
$\lambda_{1000}^*$& 9.9900$\times 10^{-7}$    & -                         & 9.9970$\times 10^{-7}$    \\
\end{tabular}
\caption{Eigenvalues and their estimates based on discrete
approximations for the kernel $h_F$ with $F$ the logistic
distribution function.
 \label{eigtable2}}
\end{center}
\end{table}


For calculation of the eigensystem from a sample, see
Section~\ref{uv}.


\subsection{Relation to Anderson-Darling kernel}\label{sec ad}

A related kernel was studied by \citeA{ad52} and \citeA{dwv73},
among others, in the context of Cram{\'e}r-von Mises tests. With $w$
a nonnegative weight function, they considered the kernel
\begin{eqnarray*}
   r_w(u,v) = \int_0^1
   [\gamma(u, t)-t]
   [\gamma(v, t)-t]
   w(t)
   dt
\end{eqnarray*}
The kernel $r_w$ is closely related to the kernel $h_F$: with
$w(t)=\frac{1}{f[F^{-1}(t)]}$, we obtain
\begin{eqnarray*}
   r_w(u,v) = h_F[F^{-1}(u),F^{-1}(v)]
\end{eqnarray*}
Note that this conversion does not work for discrete $F$. The
eigensystem for $r_w$ is given by the set of solutions
to~(\ref{de2}) with $f[F^{-1}(z)]$ replaced by the weight function
$w(z)$ \cite{ad52}.

%The results in Table~\ref{eigtable} for the uniform, the logistic and the distribution of Example~\ref{ex1} are based on previous results derived for $r_w$ given by \citeA{ad52} and \citeA{dwv73}.

%Anderson and Darling give the solutions for the uniform and the logistic cases. De Wet and Venter \citeyear[see also de Wet, 1987]{dwv73} give solutions for several other weight functions $w$.

Other closely related kernels have been given in the context of
two-sample tests by \citeA{za03} and \citeA{bf04}. (See
Section~\ref{cvm} for the relation between two-sample tests and
tests of independence.)







\section{Properties of $\kappa$ and $\rho^*$}\label{sec3}

We now apply the results of Section~\ref{sec2} in order to derive
properties of $\kappa$ and $\rho^*$. In Section~\ref{sec3.1}, some
key properties are derived, including that $0\le\rho^*(X,Y)\le 1$,
with $\rho^*(X,Y)=0$ iff $X\ind Y$ and $\rho^*(X,Y)=1$ iff $X$ and
$Y$ are linearly related. In Section~\ref{orthog} we give a
decomposition of $\kappa(X,Y)$ and $\rho^*(X,Y)$ as weighted sums of
squared correlations between the marginal eigenfunctions of
$h_{F_1}$ and $h_{F_2}$, weighted by functions of the eigenvalues.
In Section~\ref{parlik}, we describe a decomposition of the
likelihood in terms of marginal eigenfunctions and component
correlations of $\rho^*$. In Section~\ref{frechet}, Fr{\'e}chet
bound for the component correlations are given, which gives some
insight into their meaning.




\subsection{Key properties of $\kappa$ and $\rho^*$}\label{sec3.1}

Some key properties of $\kappa$ and of $\rho^*$, are given in the
following theorem:
\begin{theorem}\label{th1}
Suppose $X$ and $Y$ and $Z$ are real random variables for which the
marginal kernels $h_{F_1}$ and $h_{F_2}$ exist. Then:
\begin{enumerate}
\item If $a$, $b$, $c$ and $d$ are constants, then
$\kappa(aX+b,cY+d)=a\,c\,\kappa(X,Y)$ and
$\rho^*(aX+b,cY+d)=\rho^*(X,Y)$.\label{th1p1}

\item $\kappa(X,Y)\ge 0$ with equality iff $X\ind Y$.\label{th1p2}

\item If $\kappa(X,X)<\infty$ and $\kappa(Y,Y)<\infty$ then
$\kappa(X,Y)\le\sqrt{\kappa(X,X)\kappa(Y,Y)}$ with equality iff $X$
and $Y$ are a.s.\ linearly related.\label{th1p3}

\item If both $X$ and $Y$ are dichotomous then
$\kappa(X,Y)=\cov(X,Y)^2$ and $\rho^*(X,Y)=\rho(X,Y)^2$\label{th1p4}

\item With $Z_{i:n}$ the $i$th order statistic in a sample of size $n$, $\kappa(Z,Z)=\frac{1}{6}E(Z_{2:4}-Z_{3:4})^2$.\label{th1p5}

\end{enumerate}
\end{theorem}
The proof of the theorem is given at the end of this section. Note
that $\kappa(X,X)=Eh_{F_1}(X_1,X_2)^2$ and
$\kappa(Y,Y)=Eh_{F_2}(Y_1,Y_2)^2$ so the condition in Part~3 is
equivalent to square integrability of the marginal kernels. From
Theorem~\ref{th1}, Parts~2 and~3, we immediately have the following:
\begin{corollary}\label{cor1}
Suppose $\kappa(X,X)<\infty$ and $\kappa(Y,Y)<\infty$. Then
$0\le\rho^*(X,Y)\le 1$, with $\rho^*(X,Y)=0$ iff $X\ind Y$ and
$\rho^*(X,Y)=1$ iff $X$ and $Y$ are a.s.\ linearly related.
\end{corollary}
We may compare Corollary~\ref{cor1} to the related well-known result
for the ordinary correlation $\rho$: If $\var(X)<\infty$ and
$\var(Y)<\infty$ then $0\le\rho^2\le 1$ with $\rho^2=0$ if $X\ind Y$
and $\rho^2=1$ iff $X$ and $Y$ are a.s.\ linearly related. The
important difference is that $X\ind Y$ is equivalent to $\rho^*=0$
but $X\ind Y$ only implies $\rho=0$, not vice versa. By
Part~\ref{th1p5} of Theorem~\ref{th1}, $\kappa(Z,Z)$ can be used as
measures of dispersion for a real random variable $Z$. Note the
relation with the variance, which can be defined as
$E(Z_{1:2}-Z_{2:2})^2$.


Note that, even though $\rho^*(X,Y)=1$ iff $X$ and $Y$ are linearly
related, $\rho^*$ is not a measure of linear association in the
following sense: if the slope of the linear regression line of $Y$
given $X$ is zero, $\rho^*(X,Y)$ need not be equal to zero.


%\begin{example}
%If $Z$ has a Cauchy distribution then it can be verified that $\kappa(Z,Z)=\infty$. If either $X$ or $Y$ or both have a Cauchy distribution, then $\kappa(X,Y)$ may be finite or infinite. Then if $\kappa(X,Y)<\infty$, $\rho^*(X,Y)=0$ and if $\kappa(X,Y)=\infty$, $\rho^*$ is undefined.
%\end{example}


From Lemma~\ref{lem3} and Lemma~\ref{lemtr}, a sufficient condition
for $\rho^*(X,Y)$ to exist is that $EX$ and $EY$ exist. An example
showing that this condition is not necessary is Example~\ref{ex1}.
Note that the existence of the ordinary correlation $\rho$ has the
much stronger requirement of finite marginal variances. Summarizing,
we have
\begin{eqnarray*}
    \lefteqn{
    \rho(X,Y)\mbox{ exists }
    \Leftrightarrow
    \left\{\sigma^2(X) \mbox{ and }\sigma^2(Y)\mbox{ exist }\right\}
    \Rightarrow
    \left\{EX\mbox{ and }EY\mbox{ exist }\right\}
    \Rightarrow
    }\\
    &&
    \Rightarrow\left\{E(X_{2:4}-X_{3:4})^2\mbox{ and }E(Y_{2:4}-Y_{3:4})^2\mbox{ exist }\right\}
    \Leftrightarrow \rho^*(X,Y)\mbox{ exists }
\end{eqnarray*}
where the one-way implications are strict.

We now proceed to the proof of Theorem~\ref{th1}:\vspace{3mm}\\
\proof{Proof of Theorem~\ref{th1}}

Part 1 follows directly from the definition.

Part 2: From Lemma~\ref{lem1} we obtain
\begin{eqnarray*}
   h_{F_i}(a,b)
   &=& \int
   [\gamma(a,w)-F_i(w)][\gamma(b,w)-F_i(w)]
   dw
\end{eqnarray*}
Furthermore, note that
\begin{eqnarray*}
   {F_{12}(x,y)-F_1(x)F_2(y)}
   = E[\gamma(X,x)-F_1(x)][\gamma(Y,y)-F_2(y)]
\end{eqnarray*}
which is easy to verify. Using these results and the finiteness of
each side of~(\ref{fubini1}) which allows us to apply Fubini's
theorem, we obtain
\begin{eqnarray}
   \kappa(X,Y)
   &=& Eh_{F_1}(X_1,X_2)h_{F_2}(Y_1,Y_2)\nonumber\\
   &=& E\int[\gamma(X_1,x)-F_1(x)][\gamma(X_2,x)-F_1(x)]dx \times\nonumber\\
   &&   \int[\gamma(Y_1,y)-F_2(y)][\gamma(Y_2,y)-F_2(y)]dy\nonumber\\
   &=& \int E[\gamma(X_1,x)-F_1(x)][\gamma(Y_1,y)-F_2(y)]\times \nonumber\\
   &&E[\gamma(X_2,x)-F_1(x)][\gamma(Y_2,y)-F_2(y)]dxdy \nonumber\\
   &=& \int [F_{12}(x,y)-F_1(x)F_2(y)]^2dxdy\label{fxy}
   \label{intf12}
\end{eqnarray}
If $X\ind Y$, the integrand is zero so $\kappa(X,Y)=0$. It remains
to be shown that $X\not\!\!\!\ind Y$ implies $\kappa(X,Y)\ne 0$.  We
next sketch the proof.

If $X\not\!\!\!\ind Y$ then there is an $(a,b)$ such that
$D(a,b)=F_{12}(a,b)-F_1(a)F_2(b)\ne 0$. We now show that if
$D(a,b)\ne 0$, then there is an open interval, which has $(a,b)$ as
a limit point, such that $D\ne 0$ on that interval. It then follows
that~(\ref{intf12}) is nonzero. Let $\varepsilon_{12}\ge 0$ be the
probability mass in $(a,b)$, $\varepsilon_{1}\ge 0$ be the marginal
probability mass in $a$ and $\varepsilon_{2}\ge 0$ be the marginal
probability mass in $b$. Denote by $D(a^\pm,b^\pm)$ the limit
approaching from anywhere in one of four open `quadrants' defined by
$(a,b)$. Then from the definition of $F_{12}$, $F_1$ and $F_2$ given
in the introduction,
\begin{eqnarray*}
   D(a^-,b^-) &=& D(a,b) - \fourth\varepsilon_{12} + \half\varepsilon_{1}F_2(b) + \half\varepsilon_{2}F_1(a)
   - \fourth\varepsilon_{1}\varepsilon_{2} \\
   D(a^-,b^+) &=& D(a,b) + \fourth\varepsilon_{12} + \half\varepsilon_{1}F_2(b) - \half\varepsilon_{2}F_1(a)
   + \fourth\varepsilon_{1}\varepsilon_{2} \\
   D(a^+,b^-) &=& D(a,b) + \fourth\varepsilon_{12} - \half\varepsilon_{1}F_2(b) + \half\varepsilon_{2}F_1(a)
   + \fourth\varepsilon_{1}\varepsilon_{2} \\
   D(a^+,b^+) &=& D(a,b) + \threefourth\varepsilon_{12} - \half\varepsilon_{1}F_2(b) - \half\varepsilon_{2}F_1(a)
   - \fourth\varepsilon_{1}\varepsilon_{2}
\end{eqnarray*}
Now if $D(a,b)\ne 0$ these four expressions cannot all be zero, so
there must be an open set, in at least one of the four `quadrants'
and with $(a,b)$ as a limit point, where $D$ is nonzero.
Hence,~(\ref{intf12}) cannot be zero.




\newcommand{\aseq}{\stackrel{\mbox{\scriptsize a.s.}}{=}}


%\newpage
Part 3: By definition $\kappa(X,Y)^2\le\kappa(X,X)\kappa(Y,Y)$ is
equivalent to
\begin{eqnarray*}
   \left[Eh_{F_1}(X_1,X_2)h_{F_2}(Y_1,Y_2)\right]^2\le
   Eh_{F_1}(X_1,X_2)^2 Eh_{F_2}(Y_1,Y_2)^2
\end{eqnarray*}
This is a Cauchy-Schwartz inequality so it holds, and equality holds
iff
\begin{eqnarray}
   h_{F_1}(X_1,X_2) \aseq c\,h_{F_2}(Y_1,Y_2)\label{aseq1}
\end{eqnarray}
for some constant $c$. If $Y\aseq aX+b$ for certain constants $a$
and $b$ then it is immediately verified that~(\ref{aseq1}) holds
with $c=|a|$.


The reverse implication that~(\ref{aseq1}) implies linearity remains
to be shown. Suppose that~(\ref{aseq1}) holds. Then
\begin{eqnarray*}
   \lefteqn{h_{F_1}(X_1,X_2) - h_{F_1}(X_1,X_3) - h_{F_1}(X_2,X_4) - h_{F_1}(X_3,X_4)
   \aseq} \\
   && c\left[h_{F_2}(Y_1,Y_2) - h_{F_2}(Y_1,Y_3) - h_{F_2}(Y_2,Y_4) -
   h_{F_2}(X_3,X_4)\right]
\end{eqnarray*}
which reduces to
\begin{eqnarray*}
   \lefteqn{|X_1-X_2| - |X_1-X_3| - |X_2-X_4| + |X_3-X_4|
   \aseq} \\
   && c\left(|Y_1-Y_2| - |Y_1-Y_3| - |Y_2-Y_4| + |Y_3-Y_4|\right)
\end{eqnarray*}
But this is equivalent to
\begin{eqnarray}
   X_{3:4} - X_{2:4}
   \aseq
   c\left(Y_{3:4} - Y_{2:4}\right)
   \label{osdif}
\end{eqnarray}
Now without loss of generality suppose $Y_{3:4}=c X_{3:4}+b$ and
$Y_{2:4}=c X_{2:4}+b'$ for some $b$ and $b'$. Substitution
into~(\ref{osdif}) then yields $b=b'$, so the second and third order
statistics for $X$ and $Y$ are linearly related. Now the
distribution function of the second order statistic for $X$ is
\begin{eqnarray*}
   F_{1;2:4}(x) = \int_{-\infty}^x F_1(t)[1-F_1(t)]^2dF_1(t)
\end{eqnarray*}
and for $Y$
\begin{eqnarray*}
   F_{2;2:4}(y) = \int_{-\infty}^y F_2(t)[1-F_2(t)]^2dF_2(t)
\end{eqnarray*}
It is now straightforward to show that the equation
$F_{2;2:4}(y)=F_{1;2:4}(cx+b)$ leads to $F_2(y)=F_1(cx+b)$, so $X$
and $Y$ are linearly related.






Part 4: Without loss of generality assume $X\in\{0,1\}$ and
$Y\in\{0,1\}$ with probability one. Then $h_{F_1}=u_{F_1}$ and
$h_{F_2}=u_{F_21}$ (both $u$ and $h$ defined in Section~\ref{sec1}),
so
$\kappa(X,Y)=Eh_{F_1}(X_1,X_2)h_{F_2}(Y_1,Y_2)=Eu_{F_1}(X_1,X_2)u_{F_2}(Y_1,Y_2)=\cov(X,Y)$.


Part 5: Since $\kappa(Z,Z)=Eh_F(Z_1,Z_2)^2$, this follows directly
from Lemma~\ref{lem3}
 \hspace*{\fill}$\Box$\endproof



We conclude this section by giving some representations of $\kappa$
in terms of $h_F$ and the (conditional) distribution functions of
$X$ and $Y$. Let
\begin{eqnarray*}
   F_{2|1}(y|x)=P(Y<y|X=x)+\half P(Y=y|X=x)
\end{eqnarray*}
be the conditional distribution function of $Y$ given $X=x$.
%We now show how $\kappa(X,Y)$ can be written as a simple function of the distribution function $F_{12}$ and the conditional distribution function $G_y$.
\begin{lemma}\label{lem4}
The following equalities hold for $\kappa$:
\begin{enumerate}
\item  $\kappa(X,Y) = \int[F_{12}(x,y)-F_1(x)F_2(y)]^2dx dy$

\item
$\kappa(X,Y)=Eh_{F_1}(X_1,X_2)\int[F_{2|1}(y|X_1)-F_2(y)][F_{2|1}(y|X_2)-F_2(y)]dy$

%\item $\kappa_{K_1,K_2}(X,Y) = \int\frac{[F_{12}(x,y)-F_1(x)F_2(y)]^2}
%   {k_1[K_1^{-1}\circ F_1(x)]k_2[K_2^{-1}\circ F_2(y)]}dF_1(x)dF_2(y)$

%\item $\kappa(Z,Z)=\sixth E(Z_{2:4}-Z_{3:4})^2$
\end{enumerate}
\end{lemma}
\proof{Proof} Part 1: This follows from the proof of
Theorem~\ref{th1}, Part~\ref{th1p2}

Part 2: First note that $d_{xy}F_{12}(x,y)=d_xF_1(x)d_yF_{2|1}(y|x)$
which we write in shorthand $dF_{12}(x,y)=dF_1(x)dF_{2|1}(y|x)$.
Hence,
\begin{eqnarray*}
   \kappa(X,Y)
   &=& Eh_{F_1}(X_1,X_2)h_{F_2}(Y_1,Y_2)\\
   &=& \half\int h_{F_1}(x_1,x_2)\int[\gamma(y_1,y)-F_2(y)][\gamma(y_2,y)-F_2(y)]dydF_{12}(x_1,y_1)dF_{12}(x_2,y_2)\\
   &=& \int h_{F_1}(x_1,x_2)\int[\gamma(y_1,y)-F_2(y)]dF_{2|1}(y_1|x_1)[\gamma(y_2,y)-F_2(y)]dF_{2|1}(y_2|x_2)dydF_1(x_1)dF_1(x_2)\\
   &=& \int h_{F_1}(x_1,x_2)\int[F_{2|1}(y|x_1)-F_2(y)][F_{2|1}(y|x_2)-F_2(y)]dydF_1(x_1)dF_1(x_2)\\
   &=& Eh_{F_1}(X_1,X_2)\int[F_{2|1}(y|X_1)-F_2(y)][F_{2|1}(y|X_2)-F_2(y)]dy
\end{eqnarray*}
%Since $EG_X(y)=F_2(y)$, straightforward calculations yield the desired result.

%Parts 4 and 5: This follows immediately from Lemma~\ref{lem1},Part 3.
 \hspace*{\fill}$\Box$\endproof


Note the similarity of Part~1 of Lemma~\ref{lem4} with the formula
for the covariance given by \citeA{hoeffding40}:
\begin{eqnarray*}
   \cov(X,Y) =
   \int
   [F_{12}(x,y)-F_1(x)F_2(y)]
    dx dy
\end{eqnarray*}














\subsection{Orthogonal decomposition}\label{orthog}

Let us assume $h_{F_1}$ and $h_{F_2}$ are square integrable and have
the spectral decompositions
\begin{eqnarray}
   h_{F_1}(x_1,x_2) &=& \sum_{k=0}^{\infty}\lambda_kg_{1k}(x_1)g_{1k}(x_2)\label{spec1}\\
   h_{F_2}(y_1,y_2) &=&
   \sum_{k=0}^{\infty}\mu_kg_{2k}(y_1)g_{2k}(y_2)\label{spec2}
\end{eqnarray}
See Lemma~\ref{lem3} on how to check for square integrability. For
ease of notation, we write the correlations between marginal
eigenfunctions as
\begin{eqnarray*}
   \rho_{kl}(X,Y) = \rho[g_{1k}(X),g_{2l}(Y)]
\end{eqnarray*}
We now have the orthogonal decomposition given as follows:
\begin{theorem}\label{th2}
Suppose $h_{F_1}$ and $h_{F_2}$ are square integrable with
spectral decompositions as above. Then with convergence in mean
square,
\begin{eqnarray*}
   \kappa(X,Y) = {\sum_{k=1}^\infty\sum_{l=1}^\infty\lambda_k\,\mu_l}\,\rho_{kl}(X,Y)^2
\end{eqnarray*}
and
\begin{eqnarray*}
   \rho^*(X,Y) =
   \frac{1}{\sqrt{\sum\lambda_k^2}\,\sqrt{\sum\mu_l^2}}
   {\sum_{k=1}^\infty\sum_{l=1}^\infty\lambda_k\,\mu_l}\,\rho_{kl}(X,Y)^2
\end{eqnarray*}
\end{theorem}
\proof{Proof} Write
\begin{eqnarray*}
   \kappa^{(N,N)}(X,Y)
   &=&
   E\left[\sum_{k=1}^{N}\lambda_kg_{1k}(X_1)g_{1k}(X_2)\,\sum_{l=1}^{N}\mu_lg_{2l}(Y_1)g_{2l}(Y_2)
   \right]\\
   \kappa^{(.,N)}(X,Y)
   &=&
   E\left[h_{F_1}(X_1,X_2)\,\sum_{l=1}^{N}\mu_lg_{2l}(Y_1)g_{2l}(Y_2)
   \right]\\
   \kappa^{(N,.)}(X,Y)
   &=&
   E\left[\sum_{k=1}^{N}\lambda_kg_{1k}(X_1)g_{1k}(X_2)\,h_{F_2}(Y_1,Y_2)
   \right]
\end{eqnarray*}
Then straightforward algebra gives
\begin{eqnarray*}
   \kappa^{(N,N)}(X,Y)
   =
   {\sum_{k=1}^N\sum_{l=1}^N\lambda_k\,\mu_l}\,\rho[g_{1k}(X),g_{2l}(Y)]^2
\end{eqnarray*}
By the Cauchy-Schwartz inequality we obtain
\begin{eqnarray*}
   \lefteqn{
   \left(\kappa(X,Y) - \kappa^{(N,.)}(X,Y) - \kappa^{(.,N)}(X,Y) + \kappa^{(N,N)}(X,Y)\right)^2
   }
   \\
   &=&
   E
   \left(
   \left[h_{F_1}(X_1,X_2)-
   \sum_{k=1}^{N}\lambda_kg_{1k}(X_1)g_{1k}(X_2)
   \right]
   \left[h_{F_2}(Y_1,Y_2)-
   \sum_{l=1}^{N}\mu_kg_{2l}(Y_1)g_{2l}(Y_2)
   \right]
   \right)^2
   \\
   &\le&
   E
   \left[h_{F_1}(X_1,X_2)-
   \sum_{k=1}^{N}\lambda_kg_{1k}(X_1)g_{1k}(X_2)
   \right]^2
   E
   \left[h_{F_2}(Y_1,Y_2)-
   \sum_{l=1}^{N}\mu_kg_{2l}(Y_1)g_{2l}(Y_2)
   \right]^2
\end{eqnarray*}
By mean square convergence of the spectral decomposition the latter
goes to zero as $N\rightarrow\infty$ so
\begin{eqnarray*}
   \kappa(X,Y) - \kappa^{(N,.)}(X,Y) - \kappa^{(.,N)}(X,Y) + \kappa^{(N,N)}(X,Y) \rightarrow 0
\end{eqnarray*}
as $N\rightarrow\infty$. Similarly we find
\begin{eqnarray*}
   \kappa(X,Y) - \kappa^{(N,.)}(X,Y) \rightarrow 0\\
   \kappa(X,Y) - \kappa^{(.,N)}(X,Y) \rightarrow 0
\end{eqnarray*}
as $N\rightarrow\infty$. It follows that
\begin{eqnarray*}
   \kappa^{(N,N)} \rightarrow \kappa(X,Y)
\end{eqnarray*}
as $N\rightarrow\infty$, which is the desired result.
 \hspace*{\fill}$\Box$\endproof

The simplest example of a decomposition is if both variables are
dichotomous, say $P(X=0)=1-P(X=1)=p$ and $P(Y=0)=1-P(Y=1)=q$, we
obtain $\lambda_1=2p(1-p)$, $\mu_1=2q(1-q)$, and $\lambda_k=\mu_k=0$
for $k>1$ so that $\rho^*(X,Y) = \rho(X,Y)^2$, see
Theorem~\ref{th1}, Part~4. In this special case the decomposition
consists of just one component.



\subsection{Parameterization of the likelihood}\label{parlik}

Let $f_{12}$ be the joint density of $(X,Y)$ with corresponding
marginal densities $f_1$ and $f_2$. Since
\begin{eqnarray*}
   \rho_{kl} =
   \int
   \frac{f_{12}(x,y)}{f_1(x)f_2(y)}
   \,g_{1k}(x)g_{2l}(y)dF_1(x)dF_2(y)
\end{eqnarray*}
we can decompose the joint density as:
\begin{eqnarray}
   f_{12}(x,y)
   =
   f_1(x)f_2(y)
   \left(
   1+
   \sum_{i=1}^\infty
   \sum_{j=1}^\infty
   g_{1k}(x)g_{2l}(y)\rho_{kl}
   \right)
   \label{mod1}
\end{eqnarray}
A similar equation can be given for discrete distributions, and a
general treatment can be given using the Radon-Nikodym derivative.

Decomposition~(\ref{mod1}) may be compared to the well-known
canonical correlation decomposition
\begin{eqnarray*}
   f_{12}(x,y)
   =
   f_1(x)f_2(y)
   \left(
   1+
   \sum_{k=1}^\infty
   a_{1k}(x)a_{2k}(y)\rho_{k}
   \right)
\end{eqnarray*}
Here, $a_{1k}$ and $a_{2k}$ are those functions maximizing the
correlation between $X$ and $Y$, subject to the restraint (for
$k>1$) that they are orthogonal to $a_{11},\ldots,a_{1,k-1}$ and
$a_{21},\ldots,a_{2,k-1}$, respectively, and $\rho_k$ is the
correlation between $a_{1k}(X)$ and $b_{2k}(Y)$.










\subsection{Fr{\'e}chet bounds for component
correlations}\label{frechet}

Below, we discuss the  interpretation and properties of the
component correlations. In particular, we look at bounds for the
component correlations.


For two random variables $X$ and $Y$ with joint distribution
function $F_{12}$ and marginal distribution functions $F_1$ and
$F_2$, the well-known Fr{\'e}chet upper bound $F_{12}^+$ is defined
by
\begin{eqnarray*}
   F_{12}^+(x,y) =\min\{F_1(x),F_2(y)\}
\end{eqnarray*}
and the Fr{\'e}chet lower bound $F_{12}^-$ is defined by
\begin{eqnarray*}
   F_{12}^-(x,y) = \max\{0,1-F_1(x)-F_2(y)\}
\end{eqnarray*}
Then $\rho(X,Y)=1$ if and only if $F_{12}=F_{12}^+$ and
$\rho(X,Y)=-1$ if and only if $F_{12}=F_{12}^-$. A more general
question is, for functions $g$ and $h$, for which $F_{12}$ the
correlation between $g(X)$ and $h(Y)$ is maximal or minimal. Let
\begin{eqnarray*}
   S^+_{g,h}
   &=&
   \left\{
   (x,y)\in\bR^2|g(x)=h(y)
   \right\}
\end{eqnarray*}
and
\begin{eqnarray*}
   S^-_{g,h}
   &=&
   \left\{
   (x,y)\in\bR^2|g(x)=- h(y)
   \right\}
\end{eqnarray*}
Then we have:
\begin{lemma}
For functions $g$ and $h$, $\rho[g(X),h(Y)]=1$ iff the support of
the distribution of $(X,Y)$ is a subset of $S^+_{g,h}$ and and
$\rho[g(X),h(Y)]=-1$ iff the support is a subset of $S^-_{g,h}$.
\end{lemma}
\proof{Proof} Suppose for simplicity that $g$ and $h$ are
standardized. Then $\rho[g(X),h(Y)]=1$ iff $P[g(X)=h(Y)]=1$ and
$\rho[g(X),h(Y)]=-1$ iff $P[g(X)=-h(Y)]=1$, and the lemma
immediately follows.
 \hspace*{\fill}$\Box$\endproof

For the dichotomous case we obtain the following:
\begin{example}
If both variables are dichotomous, say $P(X=0)=1-P(X=1)=p$ and
$P(Y=0)=1-P(Y=1)=q$, we obtain $S_{11}^+=\{(0,0),(1,1)\}$ and
$S_{11}^-=\{(0,1),(1,0)\}$.
\end{example}

Note that the bounds need not be attainable since it may be the case
that, for example, for certain $x$, there is no $y$ such that
$(x,y)\in S^+_{g,h}$.

Here we are interested in the bounds for the component correlations
$\rho_{ij}$ of $\rho^*$. For simplicity, we write
$S^\pm_{ij}=S^\pm_{g_i,g_j}$. The next example shows that if both
$X$ and $Y$ have uniform distributions on $[0,1]$, then the
component correlations of $\rho^*$ can attain the bounds $1$ and
$-1$:
\begin{example}
Suppose $X$ and $Y$ are uniformly distributed on $[0,1]$. Then the
eigenfunctions are the Fourier cosine functions (see
Table~\ref{eigtable}). The set $S^+_{kl}$ is formed by the solutions
$(x,y)$ to the equation
\begin{eqnarray*}
   \cos (k\pi x) = \cos (l\pi y)
\end{eqnarray*}
and $S^-_{kl}$ by the solutions of
\begin{eqnarray*}
   \cos (k\pi x) = -\cos (l\pi y)
\end{eqnarray*}
The solutions are plotted in Figure~\ref{univ fr}, for
$k=1,\ldots,4$ and $l=1,\ldots,4$. The bounds for the $\rho_{kl}$
are attainable since $\rho_{kl}(X,Y)=1$ for the uniform distribution
on $S^+_{kl}$ and $\rho_{kl}(X,Y)=-1$ for the uniform distribution
on $S^-_{kl}$.



% [htbp] means: try to place the figure either (H)ere, at the (T)op of the page,
% at the (B)ottom of the page, or on a separate (P)age with figures only -
% options are in order of preference
\begin{figure}[htbp]
\begin{center}
 \subfigure[Positive bounds $S^+_{ij}$: $\rho_{ij}=1$ if support of $(X,Y)$ is subset of $S^+_{ij}$]
 {\includegraphics[width=70mm]{frechet1a.eps}}
\end{center}
\begin{center}
 \subfigure[Negative bounds $S^-_{ij}$: $\rho_{ij}=1$ if support of $(X,Y)$ is subset of $S^-_{ij}$]
 {\includegraphics[width=70mm]{frechet1b.eps}}
\end{center}
\caption{Supports of the Fr{\'e}chet bounds for the component
correlations $\rho_{ij}$ of $\rho^*(X,Y)$ when $X$ and $Y$ are
uniformly distributed on an interval.}\label{univ fr}
\end{figure}






Note that, if $\rho_{11}=1$, then $\rho_{22}=1$, since $S_{11}^+
\subset S_{22}^+$. More generally, by the same reasoning, we have
for all $i>1$,
\begin{eqnarray*}
   \rho_{kl}(X,Y)=1 \Rightarrow
      \rho_{i\times k,i\times l}(X,Y)=1
\end{eqnarray*}
and
\begin{eqnarray*}
   \rho_{kl}(X,Y)=-1 \Rightarrow
      \rho_{i\times k,i\times l}(X,Y)=(-1)^i
\end{eqnarray*}
By a symmetry argument, we also have
$\rho_{11}=1\Rightarrow\rho_{12}=0$, and there are various other
similar implications.
\end{example}

An overview of a large amount of literature on Fr{\'e}chet bounds is
given in \citeA{ruschendorf91}








\section{Estimation and tests of independence}\label{sec4}

In this section we discuss estimation of $\kappa$ and $\rho^*$ by U-
and V-statistics. Roughly speaking, the U-statistic estimator of a
parameter is an unbiased estimator based on taking averages
\cite{hoeffding48}, and the V-statistic estimator is the estimator
based on the distribution obtained by assigning a probability weight
$1/n$ to each sample point. For $\kappa$, both the U- and
V-statistic estimators are available, but for $\rho^*$ only the
latter is. However, we can estimate $\rho^*$ by a function of
U-statistic estimators.


In Section~\ref{uv}, it is shown how estimators of $\kappa$ and
$\rho^*$ by U- and V-statistics are obtained. In Section~\ref{perm}
permutation tests, useful for small samples, are described. In
Section~\ref{asymp}, the asymptotic distribution of these estimators
is derived under the null hypothesis of independence. In
Section~\ref{bonf}, Bonferroni corrections for tests of significance
of the component correlations are described.

%\subsection{Estimation of eigenvalues and eigenfunctions}

\subsection{U and V statistic estimators of $\kappa$}\label{uv}

We first give a method for calculating the U- and V-statistic
estimators of $\kappa$ based on a sample
$(X_1,Y_1),\ldots,(X_n,Y_n)$, then we give the related estimates for
$\rho^*$.


The V-statistic estimator $\hat\kappa$ is the value of $\kappa$
based on the sample distribution functions $\hat F_1$ and $\hat
F_2$, and is obtained as follows. Let
\begin{eqnarray*}
   A_{1k} &=& \frac{1}{n}\sum_{i=1}^n|X_k-X_i|\\
   A_{2k} &=& \frac{1}{n}\sum_{i=1}^n|Y_k-Y_i|
\end{eqnarray*}
and
\begin{eqnarray*}
   B_{1} &=& \frac{1}{n^2}\sum_{i=1}^n\sum_{j=1}^n|X_i-X_j|\\
   B_{2} &=& \frac{1}{n^2}\sum_{i=1}^n\sum_{j=1}^n|Y_i-Y_j|
\end{eqnarray*}
Then we have for $k,l=1,\ldots,n$,
\begin{eqnarray*}
   h_{\hat F_1}(x_1,x_2) &=& -\half
   \left(|x_1-x_2| - A_{1k} - A_{1l} + B_1\right)\\
   h_{\hat F_2}(y_1,y_2) &=& -\half
   \left(|y_1-y_2| - A_{2k} - A_{2l} + B_2\right)\\
\end{eqnarray*}
and the sample or V-statistic estimator of $\kappa$ is given as
\begin{eqnarray*}
   \hat\kappa = \frac{1}{n^2}\sum_{i,j=1}^nh_{\hat F_1}(X_i,X_j)h_{\hat F_2}(Y_i,Y_j)
\end{eqnarray*}
Now with
\begin{eqnarray*}
   \tilde h_{\hat F_1}(x_1,x_2) &=& -\half
   \left(|x_1-x_2| - \frac{n}{n-1}A_{1k} - \frac{n}{n-1}A_{1l} + \frac{n}{n-1}B_1\right)\\
   \tilde h_{\hat F_2}(y_1,y_2) &=& -\half
   \left(|y_1-y_2| - \frac{n}{n-1}A_{2k} - \frac{n}{n-1}A_{2l} + \frac{n}{n-1}B_2\right)\\
\end{eqnarray*}
for $k,l=1,\ldots,n$, the unbiased or U-statistic estimator of
$\kappa$ is given as
\begin{eqnarray*}
   \tilde\kappa = \frac{2}{n(n-1)}\sum_{i=1}^{n-1}\sum_{j=i+1}^n\tilde h_{\hat F_1}(X_i,X_j)\tilde h_{\hat F_2}(Y_i,Y_j)
\end{eqnarray*}
By Hoeffding's theory of U-statistics we have that $\tilde\kappa$ is
an unbiased estimator of $\kappa$ \cite{hoeffding48,rw79}. Note that
$\hat\kappa\ge 0$ but $\tilde \kappa$ may be negative.


The related estimators of $\rho^*$ are the following:
\begin{eqnarray*}
   \hat\rho^*(X,Y)   &=& \frac{\hat\kappa(X,Y)}{\sqrt{\hat\kappa(X,X)\hat\kappa(Y,Y)}}\\
   \tilde\rho^*(X,Y) &=& \frac{\tilde\kappa(X,Y)}{\sqrt{\tilde\kappa(X,X)\tilde\kappa(Y,Y)}}
\end{eqnarray*}



For both types of estimators, the computational complexity of the
above method is $O(n^2)$.


The marginal eigenvalues and functions can be computed numerically
from $h_{\hat F_1}$ and $h_{\hat F_2}$ or from  $\tilde h_{\hat
F_1}$ and $\tilde h_{\hat F_2}$. See also Section~\ref{contapprox}
for computational aspects.





\subsection{Permutation tests}\label{perm}

Under independence, the sample marginal distributions of $X$ and $Y$
are ancillary statistics for $\hat\rho^*$ and $\tilde\rho^*$, so by
Fisher's theory of fiducial inference we should condition on the
sample marginals when testing for independence using $\hat\rho^*$
and $\tilde\rho^*$. If $X\ind Y$, conditioning on the marginals
ensures that $\hat\rho^*$ and $\tilde\rho^*$ are distribution free,
and exact conditional $p$-values can be calculated using the
permutation method. Evaluating all permutations quickly becomes
computationally prohibitive even for moderately large sample sizes,
and we recommend using a set of random permutations. Note that
permutation tests may also be applied to the component correlations
$\hat\rho_{ij}$ and $\tilde\rho_{ij}$.

Permutation tests may be computationally intensive. Using
non-optimized software, our experience shows that (bootstrap)
permutation tests for up to a several hundred observations are
feasible: for $n=100$, the permutation test based on 1000 random
permutations took less than four minutes, and for $n=500$ it took a
bit more than one hour. Techniques for the fast {\em exact}
evaluation of permutation tests using generating functions are
described by, among others, \citeA{bps96} and \citeA{vdwdbvdl99},
but it is not clear whether these techniques extend to statistics
such as $\hat\rho^*$ which are not based on ranks.


%There is much literature on generating function techniques for speeding up the evaluation of permutation tests, see for example \citeA{bps96} and \citeA{vdwdb99}. However, this literature concentrates
%General techniques for the fast {\em exact} evaluation of permutation tests have been developed by Van de Wiel and Di Bucchianico \citeyear{vdwdib01}.

For categorical data the permutation test is better known as an {\em
exact conditional test} (where the conditioning is, again, on the
marginal distributions), the Fisher exact test being the best known
example. There is a large body of literature on fast evaluation of
exact conditional $p$-values for contingency tables, for an overview
see \citeA{agresti92} and for more recent developments see
\citeA{fms96,ds98,bb99}.


If the permutation test is too computationally intensive, an
asymptotic test may be done using the results of the next section.








\subsection{Asymptotic distribution of estimators under independence}\label{asymp}

For the asymptotic distribution of the estimators we obtain the
following:
\begin{theorem}\label{th asymp}
Suppose $h_{F_1}$ and $h_{F_2}$ are square integrable with spectral
decompositions~(\ref{spec1}) and~(\ref{spec2}). Then if $X\ind Y$
and with $Z_{ij}$ iid standard normal variables, we obtain
\begin{eqnarray*}
   n\tilde\kappa(X,Y) \rightarrow_D
   \sum_{i,j=0}^\infty \lambda_i\mu_j(Z_{ij}^2-1)
\end{eqnarray*}
If additionally $h_{F_1}$ and $h_{F_2}$ are trace class, we obtain
\begin{eqnarray*}
   n\hat\kappa(X,Y) \rightarrow_D
   \sum_{i,j=0}^\infty \lambda_i\mu_jZ_{ij}^2
\end{eqnarray*}
\end{theorem}
\proof{Proof} By the \citeA{hoeffding61} decomposition we can write
with $R_n=O(n^{-3})$
\begin{eqnarray*}
   \tilde\kappa
   &=&
   \left(\begin{array}{c}n \\ 2\end{array}\right)^{-1}
   \sum_{1\le i<j\le n} h_{F_1}(x_i,x_j)h_{F_2}(y_i,y_j) + R_n \\
   &=&
   \left(\begin{array}{c}n \\ 2\end{array}\right)^{-1}
   \sum_{1\le i<j\le n}
   \left( \sum_{k=1}^\infty \lambda_kg_{1k}(x_i)g_{1k}(x_j) \right)
   \left( \sum_{l=1}^\infty \mu_l g_{2l}(y_i)g_{2l}(y_j) \right) + R_n \\
   &=&
   \left(\begin{array}{c}n \\ 2\end{array}\right)^{-1}
   \sum_{k,l=1}^\infty \lambda_k\mu_l
   \left[
   \left( \sum_{i=1}^n g_{1k}(x_i)g_{2l}(y_i) \right)^2
   -
   \sum_{i=1}^n g_{1k}(x_i)^2g_{2l}(y_i)^2
   \right] + R_n
\end{eqnarray*}
Since $n^{-1}\sum_{i=1}^n g_{1k}(x_i)^2g_{2l}(y_i)^2\rightarrow 1$
a.s., and $n^{-1}\left( \sum_{i=1}^n g_{1k}(x_i)g_{2l}(y_i)
\right)^2\rightarrow_D Z_{kl}^2$ we obtain using the Cram{\'e}r-Wold
device that
\begin{eqnarray*}
   n\tilde\kappa \rightarrow_D
   \sum_{i,j=0}^\infty \lambda_i\mu_j(Z_{ij}^2-1)
\end{eqnarray*}


The proof for $\hat\kappa$ is similar; we have
\begin{eqnarray*}
   \hat\kappa
   &=&
   \frac1{n^2}
   \sum_{i,j=1}^n h_{F_1}(x_i,x_j)h_{F_2}(y_i,y_j) + R_n \\
   &=&
   \frac1{n^2}
   \sum_{i,j=1}^n
   \left( \sum_{k=1}^\infty \lambda_kg_{1k}(x_i)g_{1k}(x_j) \right)
   \left( \sum_{l=1}^\infty \mu_l g_{2l}(y_i)g_{2l}(y_j) \right) + R_n \\
   &=&
   \frac1{n^2}
   \sum_{k,l=1}^\infty \lambda_k\mu_l
   \left( \sum_{i=1}^n g_{1k}(x_i)g_{2l}(y_i) \right)^2
    + R_n
\end{eqnarray*}
Since $n^{-1}\left( \sum_{i=1}^n g_{1k}(x_i)g_{2l}(y_i)
\right)^2\rightarrow_D Z_{kl}^2$ we obtain using the Cram{\'e}r-Wold
device that
\begin{eqnarray*}
   n\hat\kappa \rightarrow_D
   \sum_{i,j=0}^\infty \lambda_i\,\mu_j\,Z_{ij}^2
\end{eqnarray*}
under the condition that
\begin{eqnarray*}
   \lim_{n\rightarrow\infty}E(n\hat\kappa)=\sum_{i,j=0}^\infty \lambda_i\,\mu_j = \sum_{i=1}^\infty \lambda_i\sum_{j=0}^\infty\mu_j
\end{eqnarray*}
is finite. Now by Lemma~\ref{lem6}, the two factors on the right
hand side are finite iff $h_{F_1}$ and $h_{F_2}$ are trace class,
completing the proof.
 \hspace*{\fill}$\Box$\endproof
The proof is similar to an adaptation by \citeA{dewet87} of a proof
by \citeA{eagleson79}. See also \citeA{gregory77} and
\citeA{hall79}.

Note that by Lemma~\ref{lemtr}, $h_{F_1}$ and $h_{F_2}$ are trace
class iff $EX$ and $EY$ exist. As follows from the theorem and noted
earlier by \citeA{dewet87} for related statistics, the U-statistic
estimator has an asymptotic distribution in more cases than the
V-statistic estimator. For example, if the marginal distribution of
at least one of $X$ and $Y$ is the distribution of
Example~\ref{ex1}, $n\hat\kappa$ does not have an asymptotic
distribution but $n\tilde\kappa$ does have one.




\subsection{Bonferroni corrections for testing significance of component correlations}\label{bonf}

As well as testing the significance of $\hat\rho^*$ directly, we can
test for the significance of the empirical component correlations
$\hat\rho_{ij}$. We recommend using $\hat\rho^*$ rather than
$\tilde\rho^*$ for calculating component correlations, since
$\tilde\rho^*$ may be negative in which case no component
correlations with nonnegative eigenvalues exist.


The proof of Theorem~\ref{th asymp} suggests that the component
correlations $\hat\rho_{kl}$ are asymptotically normal and
independent. Since there are many component correlations, a
simultaneous test of their significance needs a Bonferroni
correction. The ordinary Bonferroni correction, i.e., multiplying
the exceedance probabilities by the number of tests done, which in
this case is $n^2$, would be unreasonable since the multiplication
factor increases rapidly with $n$. Instead we propose dividing the
exceedance probability for $\hat\rho_{ij}$ by
\begin{eqnarray*}
   \frac{\hat\lambda_i\hat\mu_j}{\sum_{i=1}^n\hat\lambda_i\sum_{j=1}^n\hat\mu_j}
\end{eqnarray*}
Note that these numbers will converge to zero in probability if at
least one of $h_{F_1}$ and $h_{F_2}$ is not of trace class, i.e., by
Lemma~\ref{lemtr}, if at least one of $EX$ or $EY$ does not exist,
in which case the correction may not be the most appropriate one.


The idea of looking at components of a test seems to have first
appeared in \citeA{dk72}, where components of the Cram{\'e}r-von
Mises test were investigated. This test is a special case of the
tests based on $\rho^*$ described above (see Section~\ref{cvm}).


Other related work is by \citeA{kl99}, who looked at correlations
between orthogonal functions of the marginal cumulative distribution
functions, in particular, the Legendre polynomials. This work is an
extension of the so-called smooth tests of fit of \citeA{neyman37}.
Rather than looking at all correlations between the orthogonal
functions, they considered just the first few, and developed a
selection method based on Schwartz's rule for determining how many
correlations to base the overall test on.



\section{Grade versions of $\kappa$ and $\rho^*$, copulas, and rank
tests}\label{sec5}

For ordinal random variables $X$ and $Y$, any given scale is
arbitrary and it may be desirable to use scales based on the grades
$F_1(X)$ and $F_2(Y)$ of $X$ and $Y$. A general way to base $\kappa$
and $\rho^*$ on grades is as follows. For given invertible
distribution functions $K_1$ and $K_2$, we can define
\begin{eqnarray*}
   \kappa_{K_1,K_2}(X,Y) =
   \kappa[K_1^{-1}\circ F_1(X),K_2^{-1}\circ F_2(Y)]
\end{eqnarray*}
and
\begin{eqnarray*}
   \rho^*_{K_1,K_2}(X,Y) =
   \rho^*[K_1^{-1}\circ F_1(X),K_2^{-1}\circ F_2(Y)]
\end{eqnarray*}
Note that
\begin{eqnarray*}
   \kappa_{F_1,F_2}(X,Y) = \kappa(X,Y)
\end{eqnarray*}
and
\begin{eqnarray*}
   \rho^*_{F_1,F_2}(X,Y) = \rho^*(X,Y)
\end{eqnarray*}
With $K_1$ and $K_2$ uniform distribution functions,
$\rho^*_{K_1,K_2}(X,Y)$ is to $\rho^*$ what Spearman's rho is to the
ordinary correlation $\rho$.


We can use the results of Section~4 to obtain an orthogonal
decomposition of $\kappa_{K_1,K_2}$ and $\rho^*_{K_1,K_2}$ in terms
of component correlations. These component correlations then
parameterize the {\em copula}, which is defined as the joint
distribution of $(F_1(X),F_2(Y))$. From~(\ref{mod1}), and since the
marginal eigenfunctions of $h_F$ with $F$ the uniform distribution
are the cosine functions given in Table~\ref{eigtable}, we obtain
the following decomposition of $c_{12}$, the density function of the
copula:
\begin{eqnarray*}
   c_{12}(u,v)
   =
   1+
   \sum_{i=1}^\infty
   \sum_{j=1}^\infty
   \cos(k\pi u)\cos(l\pi v)\rho_{kl}
\end{eqnarray*}
where $\rho_{kl}=\int \cos(k\pi u)\cos(l\pi v)c_{12}(u,v)dudv$. This
decomposition was earlier given in \citeA{dewet80} and
\citeA{deheuvels81}. An overview of copula theory is given in
\citeA{nelsen06}. Possible drawbacks of using $\rho^*_{K_1,K_2}$ for
some given $K_1$ and $K_2$ is the arbitrariness of any choice of
$K_1$ and $K_2$ and the loss of scale information, but these issues
are hotly debated \cite{mikosch06}.


In Section~\ref{secrank}, a brief description of rank tests based on
$\kappa$ is given. In Section~\ref{cvm} a generalization of the
Cram{\'e}r-von Mises test to the case of $K$ ordered samples is
shown to be a special case, and a convenient representation is
given. In Section~\ref{secphi} we write $\kappa_{K_1,K_2}$ as a
weighted average of $\phi$-coefficients.





\subsection{Rank tests}\label{secrank}

Rank statistics which are distribution free under independence in
the continuous case are obtained as follows. For invertible
distribution functions $K_1$ and $K_2$ let
\begin{eqnarray*}
   \hat\rho^*_{K_1,K_2}(X,Y)   &=& \hat\rho^*[K_1^{-1}\circ \hat F_1(X),K_2^{-1}\circ \hat F_2(Y)] \\
   \tilde\rho^*_{K_1,K_2}(X,Y) &=& \tilde\rho^*[K_1^{-1}\circ \hat F_1(X),K_2^{-1}\circ \hat F_2(Y)]
\end{eqnarray*}
The derivation of the asymptotic distribution of these statistics is
slightly more involved than that of the asymptotic distribution of
$\hat\rho^*(X,Y)$. De Wet \citeyear{dewet80} has done this
derivation for statistics related to $\tilde\rho^*_{K_1,K_2}(X,Y)$.
He gave the weights for optimal tests in the Bahadur sense for
certain classes of alternatives, such as the bivariate normal.


With $K_1$ and $K_2$ the uniform distribution functions,
$n\hat\rho^*_{K_1,K_2}(X,Y)$ is a statistic discussed by
\citeA{bkr61}, see also \citeA{deheuvels81}. It can be viewed as a
generalization of the ordinary Cram{\'e}r-von Mises test (see next
subsection). Similarly, with $K_1$ and $K_2$ the logistic
distribution functions, $n\hat\rho_{K_1,K_2}^*(X,Y)$ can be viewed
as a generalization of the Anderson-Darling test.

\citeA{hoeffding48ind} described a related test, namely based on the
U-statistic estimator of
\begin{eqnarray*}
   \int \left[F_{12}(x,y)-F_1(x)F_2(y)\right]^2dF_{12}(x,y)
\end{eqnarray*}
which can be obtained from the representation of $\kappa$ in
Lemma~\ref{lem4}, Part~1, by replacing $dxdy$ by $dF_{12}(x,y)$.
Hoeffding's coefficient does not fall in the framework of the
present paper.






\subsection{A new class of $K$-sample Cram{\'e}r-von Mises tests as a special case}\label{cvm}

Suppose we have $K$ samples, the $k$th sample having $n_k$ iid
observations, say $\{U_{k1},\ldots,U_{kn_k}\}$. Then a test whether
the distributions of the observations in the different samples are
equal is called a {\em $K$-sample test}. A $K$-sample test can, in
fact, be viewed as a test of independence, namely, whether
`response' depends on `group membership,' the groups referring to
the different samples. Let us consider the case that the score
$c_k\in\bR$ is assigned to sample $k$ ($k=1,\ldots,K$) . With
$N_0=0$ and $N_k=\sum_{i=1}^kn_i$ let
$(X_{N_{k-1}+i_k},Y_{N_{k-1}+i_k})=(c_i,U_{k,i_k})$ for
$k=1,\ldots,K$ and $i_k=1,\ldots,n_k$. Then it can be seen that the
$K$ sample test is a test of independence of the $X$ observations
and the $Y$ observations. (Note that here the $X$ observations are
not random). A $K$-sample test can then be based on $\hat\rho^*$ or
$\tilde\rho^*$.

If samples are ordered but have no numerical scores assigned to
them, rank scores can be assigned, for example $c_k=N_k$.

In order to arrive at the Cram{\'e}r-von Mises test, we now use
Lemma~\ref{lem4}, Part~2 to give a representation of $\kappa$ in
terms of the conditional distribution functions. Let $G_k$ be the
distribution function of $U_k$, the response for sample $k$, and let
$p_k=n_k/N_K$ be the proportion of observations in sample $k$. Then
we obtain
\begin{eqnarray*}
   \kappa(X,Y) =\sum_{i,j}
   p_ip_jh_F(c_i,c_j)\int\left[G_i(y)-F_2(y)\right]\left[G_j(y)-F_2(y)\right]dy
\end{eqnarray*}
Some straightforward algebra shows that for the two-sample case this
reduces to
\begin{eqnarray*}
   \kappa(X,Y) =
   p_1^2p_2^2\int\left[G_1(y)-G_2(y)\right]^2dy
\end{eqnarray*}
A grade version of $\kappa$ is
\begin{eqnarray*}
   \kappa_{F_1,K}(X,Y) =\sum_{i,j}
   p_ip_jh_F(c_i,c_j)\int\left[G_i(y)-F_2(y)\right]\left[G_j(y)-F_2(y)\right]w[F_2(y)]dF_2(y)
\end{eqnarray*}
where
\begin{eqnarray*}
   w(u) = \frac{1}{k[K^{-1}(u)]}
\end{eqnarray*}
In the two-sample case, the sample version of $\kappa_{F_1,K}(X,Y)$
with $K$ the uniform distribution function reduces (essentially) to
the ordinary Cram{\'e}r-von Mises statistic, so we have a
generalization to the case of $K$ ordered samples. With $K$ the
logistic distribution, $\hat\kappa_{F_1,K}(X,Y)$ reduces to the
Anderson-Darling statistic.


%The above formulation was also found by \citeA{bf04}, who studied the case of possibly multivariate $U$ and $V$.


A different generalization of the two-sample Cram{\'e}r-von Mises
test was given by \citeA{kiefer59}, namely to the case of $K$ {\em
unordered} samples.










\subsection{$\kappa$ as a weighted $\phi$-coefficient}\label{secphi}

From Lemma~\ref{lem4}, Part~1, we directly obtain
\begin{eqnarray}
   \kappa_{K_1,K_2}(X,Y) =
   \int[F_{12}(x,y)-F_1(x)F_2(y)]^2\,dK_1^{-1}\circ
   F_1(x)\,dK_2^{-1}\circ F_2(y) \label{repr3}
\end{eqnarray}
This result leads to an interesting interpretation of
$\kappa_{K_1,K_2}$. The $\phi$ coefficient for measuring the
dependence in the $2\times 2$ table obtained by collapsing the
distribution with respect to the cut point $(x,y)$ is given as
\begin{eqnarray}
   \phi(x,y) = \frac{\left|F_{12}(x,y)-F_1(x)F_2(y)\right|}{\sqrt{F_1(x)[1-F_1(x)]F_2(y)[1-F_2(y)]}}
   \label{phidef1}
\end{eqnarray}
Now suppose $\psi$ is such that
\begin{eqnarray*}
   \phi(x,y) = \psi[F_1(x),F_2(y)]
\end{eqnarray*}
Then from~(\ref{repr3}) we obtain that $\kappa_{K_1,K_2}$ can be
written as a weighted average $\psi$-square:
\begin{eqnarray*}
   \kappa_{K_1,K_2}(X,Y) = \int\psi(u,v)^2\,w_{K_1}(u)\,w_{K_2}(v)dudv
\end{eqnarray*}
where the weight function $w$ is defined by
\begin{eqnarray*}
   w_K(u) = \frac{u(1-u)}{k[K^{-1}(u)]}
\end{eqnarray*}
The normalized weight function (integrating to one) is
\begin{eqnarray*}
  {\bar w}_K(u) = \frac{w_K(u)}{\int_0^1 w_K(u)du}
\end{eqnarray*}
where
\begin{eqnarray*}
   \int_0^1 w_K(u)du
   = \int_0^1{u(1-u)}dK^{-1}(u)
   = \int_{-\infty}^{\infty}K(x)[1-K(x)]dx
\end{eqnarray*}
In Figure~\ref{ffinv}, $\bar w_K$ is plotted for the distribution
functions $K$ given in Table~\ref{eigtable}.


\begin{figure}
\begin{center}
\includegraphics[width=.98\linewidth,clip]{ggplot.eps}
\end{center}
%\parbox{0.95\textwidth}
{\caption{Plots of $\bar w_K(u)$ for several distributions
$K${\label{ffinv}}}}
\end{figure}


From Figure~\ref{ffinv} we can deduce which $\phi$-coefficients are
given most weight for different marginal distributions. As the
reference marginal, we take the logistic, which is the horizontal
line in the figure, i.e., assigning uniform weights. We see that for
a uniform marginal, the weight goes to zero in the tails. The
weights for a normal marginal are for most $u$ intermediate between
the weights for the uniform and logistic marginals. In contrast to
uniform and normal marginals, the Laplace distribution gives more
weight to the tails than a logistic marginal. An exponential
marginal gives little weight to the lower tail, but much weight to
the upper tail. Finally, the chi-square distribution gives very
large weight to the lower tail and very small weight to the upper
tail. Among the distributions considered, the biggest difference is
between the chi-square and the exponential distribution.








\section{Data analysis: investigating the nature of the association}\label{sec6}

Gaining an understanding of the nature of an association between two
random variables is probably best viewed as an art rather than a
science, and in this section we present some visual tools based on
$\rho^*$ and its components which may be helpful in reaching this
objective. For an iid bivariate sample $\{(X_i,Y_i)\}$ we propose
the following two procedures.


Firstly, we calculate $\hat\rho^*$ from the sample and test its
significance. If found to be significant, then for each data point
$(X_i,Y_i)$ we calculate the weight
\begin{eqnarray*}
   W_{i} = \frac{\frac1n\sum_{j=1}^n h_{\hat F_1}(X_i,X_j)\,h_{\hat F_2}(Y_i,Y_j)}{\sqrt{\hat\kappa(X,X)\,\hat\kappa(Y,Y)}}
\end{eqnarray*}
Since
\begin{eqnarray*}
   \hat\rho^*(X,Y) = \frac{1}{n}\sum_{i=1}^nW_{i}
\end{eqnarray*}
the weights $W_{i}$ give an indication of how much the sample
element $(X_i,Y_i)$ contributes to $\hat\rho^*$, and so can be used
to discover the nature of a possible association between $X$ and
$Y$.

Secondly, we calculate the component correlations $\hat\rho_{kl}$ of
$\hat\rho^*$ and test their significance using the Bonferroni
corrections described in Section~\ref{bonf}. Then for each
significant component correlation $\hat\rho_{kl}$, we compute the
weights
\begin{eqnarray*}
   W_{i}^{(k,l)} = \hat g_{1k}(X_i)\hat g_{2l}(Y_i)
\end{eqnarray*}
where $g_{1k}$ and $g_{2l}$ are the eigenfunctions belonging to
$h_{F_1}$ and $h_{F_2}$. Since
\begin{eqnarray*}
   \hat\rho_{kl} = \frac{1}{n}\sum_{i=1}^nW_{i}^{(k,l)}
\end{eqnarray*}
the weight $W_{i}^{(k,l)}$ is the amount the sample element
$(X_i,Y_i)$ contributes to $\hat\rho_{kl}$ (conditionally on the
marginals), and so, like $W_i$, can be used to investigate the
association between $X$ and $Y$.


In this section we show how to visualize the weights $W_i$ and
$W_i^{(k,l)}$, both for continuous and categorical data, and show
how this can be used to gain an understanding of the association.
Some artificial continuous data sets are considered in
Section~\ref{sec art}, a real categorical data set is considered in
Section~\ref{sec cat} and a real time series data set is considered
in Section~\ref{sec time}


\subsection{Some artificial data sets}\label{sec art}

% [htbp] means: try to place the figure either (H)ere, at the (T)op of the page,
% at the (B)ottom of the page, or on a separate (P)age with figures only -
% options are in order of preference
\begin{figure}[tbp]
 \subfigure[Increasing location ($\hat\rho^*=.36$)]{\includegraphics[width=60mm]{newcorex1.eps}}
\hfill
 \subfigure[Decreasing then increasing location ($\hat\rho^*=.17$)]{\includegraphics[width=60mm]{newcorex2.eps}}
\\
 \subfigure[Increasing dispersion ($\hat\rho^*=.11$)]{\includegraphics[width=60mm]{newcorex3.eps}}
\hfill
 \subfigure[Increasing then decreasing dispersion ($\hat\rho^*=.02$)]{\includegraphics[width=60mm]{newcorex4.eps}}
\caption{Scatterplots of artificial data sets. The captions denote
what happens to the conditional $Y$ distribution as $X$
increases.}\label{4scats}
\end{figure}


In Figure~\ref{4scats}, four artificial data sets are plotted, each
consisting of 100 iid points. For completeness, we explain how the
data were generated. In the following, $U$ is uniformly distributed
on $[0,1]$ and $Z_1$ and $Z_2$ are iid standard normal random
variables, and $Z(u)$ has a normal distribution with mean zero and
standard deviation $u$. The data in Figure~\ref{4scats}(a) are from
a bivariate normal distribution with $\rho=2/3$. The data in
Figure~\ref{4scats}(b) are of the form
$(X,Y)=(U,(U-1/2)^2)+(Z_1/10,Z_2/10)$. The data in
Figure~\ref{4scats}(c) are of the form $(X,Y)=(U,Z(1/5+U))$.
Finally, the data in Figure~\ref{4scats}(d) are of the form
$(X,Y)=(U,Z(1/5+\min\{U,1-U\}))$.


For all four data sets, we performed permutation tests for the
significance of $\hat\rho^*$ and its component correlations
$\hat\rho_{ij}$ based on 10,000 random permutations. This took us
about 48 minutes for each data set. For the significant component
correlations we took 1 million random permutations to get a more
accurate $p$-value, and this took about 110 seconds per component
correlation. We also computed the ordinary correlation, and, not
surprisingly, only for the data in Figure~\ref{4scats}(a) it is
significantly different from zero. There we found that
$\hat\rho=.44$ ($p=.000$).


% [htbp] means: try to place the figure either (H)ere, at the (T)op of the page,
% at the (B)ottom of the page, or on a separate (P)age with figures only -
% options are in order of preference
\begin{figure}[tbp]
 \subfigure[Increasing location ($\hat\rho^*=.36$)]{\includegraphics[width=60mm]{newcorweights1.eps}}
\hfill
 \subfigure[Decreasing then increasing location ($\hat\rho^*=.17$)]{\includegraphics[width=60mm]{newcorweights2.eps}}
\\
 \subfigure[Increasing dispersion ($\hat\rho^*=.11$)]{\includegraphics[width=60mm]{newcorweights3.eps}}
\hfill
 \subfigure[Increasing then decreasing dispersion ($\hat\rho^*=.02$)]{\includegraphics[width=60mm]{newcorweights4.eps}}
\caption{Representation of weights $W_i$ for data in
Figure~\ref{4scats}. The size of each dot is proportional to
$|W_i|$; black dots represent positive $W_i$ and white represent
negative $W_i$. Total area of dots is scaled to a fixed constant for
each plot.}\label{4scats2}
\end{figure}


For the data in Figure~\ref{4scats}(a), we found that
$\hat\rho^*=.36$ ($p=.000$), i.e., there is significant association.
In Figure~\ref{4scats2}(a), the data are plotted again, this time
each sample element $(X_i,Y_i)$ is represented with size
proportional to $|W_{i}|$; black dots represent positive $W_i$ and
white represent negative $W_i$. In all these plots, the total area
of the dots is scaled to a constant to make it easier to study them.
From Figure~\ref{4scats2}(a),we see that the association seems to
consist of a linearity in the data. It may be worthwhile however to
check if there are other forms of association present by looking at
the component correlations of $\hat\rho^*$. We found two significant
components: $\hat\rho_{11}=.61$ ($p=.000$) and $\hat\rho_{22}=.38$
($p=.004$). In Figures~\ref{4scats3}(a) and~\ref{4scats3}(b) the
weights $W_{i}^{(11)}$ and $W_{i}^{(22)}$ are visualized. The
interpetation of the black and white dots is the same as above. The
gridlines correspond to the zeroes of the marginal eigenfunctions.
Therefore, within any rectangle the dots have the same color. Also
in this case, the plots point to a linearity in the data.


For the data in Figure~\ref{4scats}(b), we found $\hat\rho^*=.17$
($p=.000$) and we found two significant component correlations:
$\hat\rho_{21}=-.78$ ($p=.000$) and $\hat\rho_{42}=.54$ ($p=.000$).
The plots in Figures~\ref{4scats2}(b),~\ref{4scats3}(c)
and~\ref{4scats3}(d) all point to a curved relationship.

For the data in Figure~\ref{4scats}(c), we found that
$\hat\rho^*=.11$ ($p=.001$). There is significant association at the
5\% level, but the evidence is not as overwhelming as in the
previous two cases. We only found one significant component
correlation: $\hat\rho_{12}=-.51$ ($p=.000$).
Figures~\ref{4scats2}(c) and~\ref{4scats3}(e) both point towards and
increase in dispersion of the $Y$ variable as $X$ increases.


For the data in Figure~\ref{4scats}(d), we found that
$\hat\rho^*=.02$ ($p=.522$). This time, the test based on $\rho^*$
does not yield a significant association. However, there is one
significant component correlation: $\hat\rho_{22}=-.38$ ($p=.004$).
Figure~\ref{4scats3}(f) indicates that the association is due to an
increase and then decrease in dispersion. Since $\hat\rho^*$ is not
significant, we should refrain from giving an interpretation based
on Figure~\ref{4scats2}(d).


% [htbp] means: try to place the figure either (H)ere, at the (T)op of the page,
% at the (B)ottom of the page, or on a separate (P)age with figures only -
% options are in order of preference
\begin{figure}[tbp]
 \subfigure[$\hat\rho_{11}=.61$ (data from Figure~\ref{4scats}(a))]
 {\includegraphics[width=60mm]{compcorweights1a.eps}}
\hfill
 \subfigure[$\hat\rho_{22}=.38$  (data from Figure~\ref{4scats}(a))]
 {\includegraphics[width=60mm]{compcorweights1b.eps}}
\\
 \subfigure[$\hat\rho_{21}=-.78$  (data from Figure~\ref{4scats}(b))]
 {\includegraphics[width=60mm]{compcorweights2a.eps}}
\hfill
 \subfigure[$\hat\rho_{42}=.54$  (data from Figure~\ref{4scats}(b))]
 {\includegraphics[width=60mm]{compcorweights2b.eps}}
\\
 \subfigure[$\hat\rho_{12}=-.51$ (data from Figure~\ref{4scats}(c))]
 {\includegraphics[width=60mm]{compcorweights3a.eps}}
\hfill
 \subfigure[$\hat\rho_{22}=-.38$  (data from Figure~\ref{4scats}(d))]
 {\includegraphics[width=60mm]{compcorweights4a.eps}}
\caption{Representation of weights $W_i^{(k,l)}$ contributing to
$\hat\rho_{kl}$ for data in Figure~\ref{4scats}. The meaning of the
dots is otherwise the same as in
Figure~\ref{4scats2}.}\label{4scats3}
\end{figure}




For comparative purposes, we plotted the weights $W_i$ for data
drawn from a distribution satisfying independence in
Figure~\ref{4scats4}. Both sets consist of 100 sample elements. The
most common pattern is that of Figure~\ref{4scats4}(a), with two
diagonally opposing clusters of black dots and two diagonally
opposing clusters of white dots. In a very limited investigation,
this type of pattern occured about half of the time. Otherwise more
complex patterns were obtained, such as the one in
Figure~\ref{4scats4}(b). These figures indicate that it doesn't seem
to make sense to interpret this kind of plot if $\hat\rho^*$ is not
significant, such as Figure~\ref{4scats2}(d).


% [htbp] means: try to place the figure either (H)ere, at the (T)op of the page,
% at the (B)ottom of the page, or on a separate (P)age with figures only -
% options are in order of preference
\begin{figure}[tbp]
 \subfigure[Most common pattern]{\includegraphics[width=60mm]{newcorweightsindep1.eps}}
\hfill
 \subfigure[Less common pattern]{\includegraphics[width=60mm]{newcorweightsindep2.eps}}
\caption{Representation of weights $W_i$ for data drawn from
distributions satisfying independence.}\label{4scats4}
\end{figure}


Concluding, we see that significance tests combined with an
inspection of the two types of plots in Figures~\ref{4scats2}
and~\ref{4scats3} can give us insight into the kind of association
present between two random variables.






\subsection{Mental health data}\label{sec cat}

\begin{table}
\begin{center}
\begin{tabular}{|l|c|c|c|c|}\hline
  &\multicolumn{4}{|l|}{Mental Health Status}\\  \hline
Parents'      &      & Mild      & Moderate  &  \\
Socioeconomic &      & Symptom   & Symptom   &      \\
Status        & Well & Formation & Formation & Impaired\\ \hline %
A (high) & 64 & 94& 58& 46 \\  %
B & 57 & 94 & 64 & 40\\  %
C & 57 & 105 & 65 & 60\\ %
D & 72 & 141 & 77 & 94\\ %
F & 36 & 97 &54 & 78\\  %
G (low) & 21 & 71 & 54 & 71\\ \hline
\end{tabular}
\end{center}\caption{Cross-classification of Mental Health Status and Socioeconomic Status}
\label{table ment}
\end{table}

% [htbp] means: try to place the figure either (H)ere, at the (T)op of the page,
% at the (B)ottom of the page, or on a separate (P)age with figures only -
% options are in order of preference
\begin{figure}[htbp]
\begin{center}
 \subfigure[$\hat\rho^*=.02$]
 {\includegraphics[width=36mm]{mentalhealth1.eps}}
 \hspace{6mm}
 \subfigure[$\hat\rho_{11}=.13$]
 {\includegraphics[width=36mm]{mentalhealth2.eps}}
 \hspace{6mm}
 \subfigure[$\hat\rho_{13}=.08$]
 {\includegraphics[width=36mm]{mentalhealth3.eps}}
 \hfill
\end{center}
\caption{Representation of weights $W_{ab}$, $W_{ab}^{(11)}$ and
$W_{ab}^{(13)}$, respectively, for mental health data of
Table~\ref{table ment}. The darker the shade of gray in cell $(a,b)$
the larger $W_{ab}$; $W_{ab}=0$ is represented by an intermediate
shade of gray.} \label{table ment2}
\end{figure}


Table~\ref{table ment} describes the relationship between child's
mental impairment and parents' socioeconomic status  for a sample of
residents of Manhattan \cite[and references
therein]{goodman85,agresti02}. Goodman used this table to illustrate
various association models for categorical data, including the
so-called linear by linear association model, the row and columns
effects model, and correspondence analysis based on canonical
correlations. Here we illustrate the use of $\hat\rho^*$ and its
components as yet an alternative method for analyzing these data. We
relied on asymptotic $p$-values because with 1670 observations
approximate evaluation of the permutation tests would have been too
time consuming using our implementation.


In the categorical case, for an $I\times J$ contingency table, it
suffices to calculate weights for the cells, i.e., it is not
necessary to calculate separately a weight for each individual
observation. For an observation $(X_i,Y_i)$ in cell $(a,b)$, the
weight $W_i$ reduces to
\begin{eqnarray*}
   W_{ab} = p_{ab}\frac{\sum_{i=1}^I\sum_{j=1}^Jp_{ij}h_{\hat F_1}(i,a)h_{\hat F_2}(j,b)}{\sqrt{\hat\kappa(X,X)\hat\kappa(Y,Y)}}
\end{eqnarray*}
where $p_{ab}$ is the proportion of observations in cell $(a,b)$.
Similarly, the weights belonging to component correlation
$\rho_{kl}$ are
\begin{eqnarray*}
   W_{ab}^{(k,l)} = p_{ab}g_{1k}(a)g_{2l}(b)
\end{eqnarray*}
for $a=1,\ldots,I$, $b=1,\ldots,J$, $k=1,\ldots,I$ and
$l=1,\ldots,J$. Note that
\begin{eqnarray*}
   \rho^* = \sum_{a=1}^I\sum_{b=1}^J W_{ab}
\end{eqnarray*}
and
\begin{eqnarray*}
   \rho_{kl} = \sum_{a=1}^I\sum_{b=1}^J W_{ab}^{(k,l)}
\end{eqnarray*}


We found that $\hat\rho^*=.02$ ($p=.000$), i.e., there is
significant association in the data. The weights $W_{ab}$ for the
cells are represented in Figure~\ref{table ment2}(a). Here, the
grayscale represents the size of $W_{ab}$: the darker the cell, the
larger $W_{ab}$; $W_{ab}=0$ is represented by a fixed intermediate
shade of gray. From Figure~\ref{table ment2}(a), it can be seen that
most of the association is of a monotone nature: the higher the
parents' socioeconomic status, the better the mental health status
of their children. We also investigated the component correlations
and found two components to be significant at the 5\% level:
$\hat\rho_{11}=.13$ ($p=.000$) and $\hat\rho_{13}=.08$ ($p=.026$).
In Figures~\ref{table ment2}(b) and~\ref{table ment2}(c) we
represented the $W_{ab}^{(11)}$ and $W_{ab}^{(13)}$ using grayscales
as above. From Figure~\ref{table ment2}(b), we see that
$\hat\rho_{11}$ indicates linearity again. However, in
Figure~\ref{table ment2}(c) we see evidence of some nonlinearity in
the data, namely an apparent reversal of the association if only the
middle categories `Mild Symptom Formation' and `Moderate Symptom
Formation' are considered. Hence, it appears that the association
which is present in the data cannot be fully explained by linearity.






\subsection{Norwegian stock exchange}\label{sec time}

% [htbp] means: try to place the figure either (H)ere, at the (T)op of the page,
% at the (B)ottom of the page, or on a separate (P)age with figures only -
% options are in order of preference
\begin{figure}[htbp]
\begin{center}
 \subfigure[Time series of stock price index]
 {\includegraphics[width=7cm]{norwayindex1.eps}}
\end{center}
\hfill \\
\begin{center}
 \subfigure[Scatterplot of $(Z_t,W_t)$ showing association between successive jumps in Figure~\ref{fig norway}(a).
 $\hat\rho^*=.09$]
 {\includegraphics[width=7cm]{norwayindex2.eps}}
\end{center}
\hfill \\
\begin{center}
 \subfigure[Weights $W_i$ for data in Figure~\ref{fig norway}(b)]
 {\includegraphics[width=7cm]{norwayweights.eps}}
\end{center}
\caption{Monthly Norwegian stock price indices, 1914-2001}
\label{fig norway}
\end{figure}

The Norwegian stock exchange data represented in Figure~\ref{fig
norway}(a) yield an example with an especially interesting form of
association. In Figure~\ref{fig norway}(a) the original time series
data $Y_t$ are plotted. In Figure~\ref{fig norway}(b) we plotted
\begin{eqnarray*}
   (Z_t,W_t)=(\mbox{arcsinh}(Y_{t}-Y_{t-1}),\mbox{arcsinh}(Y_{t+1}-Y_t))
\end{eqnarray*}
The arcsinh transformation was done to make the marginal
distributions less heavy tailed. We found a highly significant
association with $\hat\rho^*(Z,W)=.09$ ($p=.000$). From
Figure~\ref{fig norway}(b) we can already interpret the association:
large jumps (up or down) in stock prices tend to be followed by
large jumps, and small jumps by small jumps, indicating periods of
volatility. In Figure~\ref{fig norway}(c), the weights $W_{i}$ are
represented by the size and color of the dots (see Section~\ref{sec
art} for further explanation). This plot points to the same
conclusion that large jumps are followed by large jumps and small
jumps by small jumps. Note that in this case, not only the positive
weights (the black dots) but also the negative weights (the white
dots) are highly indicative of association. The plot indicates that
the up-arm (with the black dots) is `heavier' than the down-arm
(with the white dots), that is, there is evidence that in the data
generating process a jump tends to be of the same sign as the
previous jump.


Seven component correlations were found to be significant at the 5\%
level after applying the Bonferroni correction: $\hat\rho_{11}=.25$,
$\hat\rho_{12}=.13$, $\hat\rho_{22}=.64$, $\hat\rho_{24}=.19$,
$\hat\rho_{33}=.19$, $\hat\rho_{44}=.38$ and $\hat\rho_{66}=.25$,
all with $p=.000$. In Figure~\ref{fig norway 2}, the weights
corresponding to the component correlations are represented.
Figures~\ref{fig norway 2}(c), (d), (f) and (g) point to the
cross-like nature of the data. Figures~\ref{fig norway 2}(a) and (e)
indicate that the up-arm is heavier than the down-arm. We did not
find a meaningful explanation for Figure~\ref{fig norway 2}(b).



\begin{figure}[htbp]
\begin{center}
 \subfigure[$\hat\rho_{11}=.25$]
 {\includegraphics[width=4cm]{norwaystocks11.eps}}
\hfill
 \subfigure[$\hat\rho_{12}=.13$]
 {\includegraphics[width=4cm]{norwaystocks12.eps}}
\hfill
 \subfigure[$\hat\rho_{22}=.64$]
 {\includegraphics[width=4cm]{norwaystocks22.eps}}
\end{center}\hfill\\
\begin{center}
 \subfigure[$\hat\rho_{24}=.19$]
 {\includegraphics[width=4cm]{norwaystocks24.eps}}
\hfill
 \subfigure[$\hat\rho_{33}=.19$]
 {\includegraphics[width=4cm]{norwaystocks33.eps}}
\hfill
 \subfigure[$\hat\rho_{44}=.38$]
 {\includegraphics[width=4cm]{norwaystocks44.eps}}
\end{center}\hfill\\
\begin{center}
 \subfigure[$\hat\rho_{66}=.25$]
 {\includegraphics[width=4cm]{norwaystocks66.eps}}
\end{center}\hfill\\
\caption{Representation of weights contributing to $\hat\rho_{ij}$
for data in Figure~\ref{fig norway}(b)} \label{fig norway 2}
\end{figure}



\subsection{Discussion}

If a researcher investigating the association between two variables
decides on the use of $\rho^*$, we recommend the following approach.
First a test of the significance of $\hat\rho^*$ should be done,
and, if found to be significant, the weights $W_i$ should be
visualized as described above in order to determine the nature of
the association. If this does not yield the desired insight, it can
be worthwhile to investigate the component correlations and
visualize the corresponding weights $W_i^{(k,l)}$. These components
form a (unique) orthogonal decomposition of the
`infinite-dimensional' object $\rho^*$ into `one-dimensional'
objects $\rho_{kl}$, and the orthogonality ensures, in a limited
sense, that the different components measure different things; by
the latter we mean that for large samples and close to independence,
the sample component correlations are approximately independent. The
question may arise: why not investigate correlations between other
sets of orthogonal functions? A sketch of an answer is as follows.
Because of the various optimality properties of the eigenfunctions
in describing the marginal kernels, these component correlations are
likely to be a better choice for investigating the deviation of
$\rho^*$ from zero than correlations between arbitrarily chosen
functions for the marginal distributions. One way to make this
intuitive is as follows: in those regions where the marginal
distributions are sparse, the eigenfunctions vary relatively slowly
(in second derivative sense, see remark after Lemma~\ref{cor2}), and
therefore power of a test based on a component correlation will be
concentrated in those regions where there are many observations.
Thus, if we believe $\rho^*$ to be a good measure of deviation from
independence, then good `one-dimensional' objects to look at are the
correlations between the marginal eigenfunctions of $h_{F_1}$ and
$h_{F_2}$.




%\subsection{Comparison with $\phi^2$ and canonical correlation analysis}

%\newpage bcjidwfjwo \newpage




\section*{Acknowledgements}

The author would like to thank the following institutes where he has
been employed for providing a stimulating environment in which to
carry out this research: The Methodology Department of the Faculty
of Social Sciences at Tilburg University (Tilburg, the Netherlands)
and EURANDOM (Eindhoven, the Netherlands).


%The author wishes to thank Alessandro Di Bucchianico and Henry Wynn
%for useful comments, discussions and encouragement.

\bibliographystyle{apacite}
\bibliography{stats}

\end{document}



\section{Conclusion}\label{sec7}

The main aim of this paper is to contribute a tool for understanding
the nature of the association between two real random variables,
rather than the {\em testing} of independence.


\subsection{Global temperature}

\begin{figure}
\begin{center}
\includegraphics[width=.6\linewidth,clip]{co2temp.eps}
\end{center}
%\parbox{0.95\textwidth}
{\caption{CO$_2$ concentration and temperature {\label{co2temp}}}}
\end{figure}

























































\newpage



\section{Misc, possibly of use}


INTERESTING RESULT:$Eh_F(Z,Z)^2=\var(Z)$?????????????????



\begin{theorem}
Let $Z_{i:k}$ be the $i$th order statistic of $k$ iid observations.
Then with $\langle g_1,g_2\rangle=Eg_1(Z_1)g_2(Z_2)h(Z_1,Z_2)$ we
have:
\begin{eqnarray*}
   ||g||^2 =
   E \left(Z_{3:4}-Z_{2:4}\right)
   \left(
   \left[g(Z_{1:4})-g(Z_{3:4})\right]
   \left[g(Z_{2:4})-g(Z_{4:4})\right]
   +
   \left[g(Z_{1:4})-g(Z_{4:4})\right]
   \left[g(Z_{2:4})-g(Z_{3:4})\right]
   \right)
\end{eqnarray*}
and
\begin{eqnarray*}
   \langle g_1,g_2\rangle =
   E \left(Z_{3:4}-Z_{2:4}\right)
   \left(
   \left[g_1(Z_{1:4})-g_1(Z_{3:4})\right]
   \left[g_2(Z_{2:4})-g_2(Z_{4:4})\right]
   +
   \left[g_1(Z_{1:4})-g_1(Z_{4:4})\right]
   \left[g_2(Z_{2:4})-g_2(Z_{3:4})\right]
   \right)
\end{eqnarray*}
\end{theorem}



\section{OLD}


\subsection{Discrete approximation of the eigensystem}

(replace by difference equation approach which is much faster!?


Suppose $Z$ is discrete with $P(Z_i=c_i)=p_i$ for constants
$c_i\in\bR$ ($i=1,\ldots,t$, $c_i\ne c_j$ if $i\ne j$), and
$\sum_{i=1}^t p_i=1$. Let $\bc=(c_1,\ldots,c_t)^T$,
$\bp=(p_1,\ldots,p_t)^T$, $\bC$ be the $t\times t$ matrix with
$(i,j)$th coordinate $|c_i-c_j|$, $\bE =\left( I -
\bone\cdot\bp^T\right)$ and
\begin{eqnarray*}
   \bH = \half \bE\,\bC\,\bE^T
\end{eqnarray*}
Then it is straightforward to verify that $h_F(c_i,c_j)$ is the
$(i,j)$th coordinate of $\bH$. The eigenvalues of $h_F$ are the
ordinary eigenvalues of $\bH\bD_\bp$, where $\bD_\bp$ is the
diagonal matrix with $\bp$ on the main diagonal. It can be verified
that f $p_i>0$ for $i=1,\ldots,t$, then $h_F$ has $t-1$ nonzero
eigenvalues. If the distribution of $Z$, $\bp$, is known, the
eigenvalues can be calculated numerically. For $t=1000$ we were able
to calculate the eigenvalues in approximately 25 seconds using {\em
Mathematica 5.1} and a Pentium IV running at 2.8 MHz with one
gigabyte of internal memory. The largest value of $t$ for which we
were able to obtain solutions was $t=6000$.


\begin{lemma}
Let $\lambda^{(t)}_k$ be the $k$th eigenvalue based on a discrete
approximation based on $t$ points. Then if $\sum\lambda_k<\infty$,
$\lambda^{(t)}_k\rightarrow\lambda_k$ as $t\rightarrow\infty$.
\end{lemma}
\proof{Proof} Possible line of proof: numerical eigenvalues converge
to finite sum...
 \hspace*{\fill}$\Box$\endproof





\subsection{Some inequalities}

From Theorem~\ref{th2} and Lemma~\ref{lem1}, we immediately obtain
\begin{eqnarray*}
   \rho^*(X,Y) \ge \lambda_1^*\mu_1^*\rho[g_{11}(X),h_1(Y)]^2
\end{eqnarray*}
and
\begin{eqnarray*}
   \rho_{K_1,K_2}^*(X,Y) \ge \lambda_1^*\mu_1^*\rho_{K_1,K_2}[g_{11}(X),h_1(Y)]^2
\end{eqnarray*}
where $\rho_{K_1,K_2}(X,Y)=\rho(K_1\circ F_1(X),K_2\circ F_2(Y))$.
With $K$ the uniform distribution function, we obtain
\begin{eqnarray*}
   \rho_{K,K}^*(X,Y)
   \,
   \ge \,\frac{90}{\pi^4}\rho_{K,K}[\cos(\pi X),\cos(\pi Y)]^2
\end{eqnarray*}
where $\frac{90}{\pi^4} \approx 0.924$. Remember that for uniform
$K$, $\rho_{K,K}$ is also called Spearman's rho. Now it is easy to
show that
\begin{eqnarray*}
   \rho_{K,K}[\cos(\pi X),\cos(\pi Y)]\,
   \ge \,\frac23\rho_{K,K}(X,Y)
\end{eqnarray*}
so that
\begin{eqnarray*}
   \rho_{K,K}^*(X,Y)
   \ge
    \frac{40}{\pi^4}\rho_{K,K}(X,Y)^2
\end{eqnarray*}
where $40/\pi^4\approx0.412$.


Alternatively, if $K$ is the logistic distribution function, we
obtain $g_{11}(x)=h_1(x)=x$ with $\lambda_1=\mu_1=\half$. From the
eigen expansion it immediately follows that
\begin{eqnarray*}
   \rho_{K,K}^*(X,Y)
   \ge
   \frac{3}{4(\pi^2-9)}\rho_{K,K}[X,Y]^2
\end{eqnarray*}
where $\frac{3}{4(\pi^2-9)}\approx 0.862$.




Especially for the field of categorical data analysis, canonical
correlation has been used for defining parsimonious models:
\begin{eqnarray*}
   \Omega
   =
   1+
   \sum_{j=1}^M
   u_{1j}(x)u_{2j}(y)\rho_{j}
\end{eqnarray*}
The same can be done with our approach, for example,
\begin{eqnarray*}
   \Omega
   =
   1+
   \sum_{k+l\le M+1}
   g^{(M)}_{1k}(x)g^{(M)}_{2l}(y)\rho^{(M)}_{kl}
\end{eqnarray*}
for some $M>0$. For $M=1$, the two models coincide. Note that the
calculation of the correlations and eigenfunctions and values is
done in a different way than above.





\begin{figure}
\begin{center}
\includegraphics[width=.25\linewidth,clip]{fr1.eps}
\end{center}
%\parbox{0.95\textwidth}
{\caption{Logistic}}
\end{figure}


\begin{figure}
\begin{center}
\includegraphics[width=.25\linewidth,clip]{fr3.eps}
\end{center}
%\parbox{0.95\textwidth}
{\caption{Exponential}}
\end{figure}




\subsection{Interpretation of components and Fr{\'e}chet bounds}

Below, we discuss the  interpretation and properties of the
component correlations. In particular, we look at bounds for the
component correlations.


We have:
\begin{theorem}
Suppose $\sum\lambda_k<\infty$ and $\sum\mu_l<\infty$ (maybe
condition not necessary????). Then $\rho^*(X,Y)=1$ iff
$\rho_{11}(X,Y)^2=1$. If $\rho_{11}(X,Y)=1$ then
$\rho_{kl}(X,Y)=\delta_{kl}$ for all $k,l\ge 1$, and if
$\rho_{11}(X,Y)=-1$ then $\rho_{kl}(X,Y)=(-1)^k\delta_{kl}$ for all
$k,l\ge 1$.
\end{theorem}
Summarizing, we have four equivalent statements (conditions!!??):
\begin{enumerate}
\item $X$ and $Y$ are a.s.\ linearly related

\item $\rho(X,Y)^2=1$

\item $\rho^*(X,Y)=1$

\item $\rho_{11}(X,Y)^2=1$
\end{enumerate}





Cite Ruschendorf on Fr{\'e}chet bounds.....




\begin{example}
Suppose $X$ and $Y$ have logistic distributions. Let $\{g_{k}\}$ be
the eigenfunctions belonging to $h_K$ with $K$ the logistic
distribution function. For $k$ and $l$ ranging from $1$ to $5$,
$S^+_{g,h}$ and $S^-_{g,h}$ are plotted in Figure~.... It can be
seen that it is possible that $\rho[g_{k}(X),g_l(Y)]=1$ {\em only}
for $k=l$, and $\rho[g_{k}(X),g_l(Y)]=-1$ is only possible for odd
$k=l$. There are fewer implications than in the uniform case. Note
that for $k>1$
\begin{eqnarray*}
   S_{g_1,g_1}^+ &\subset& S_{g_{k},g_{k}}^+ \\
   S_{g_2,g_2}^+ &\subset& S_{g_{2k},g_{2k}}^+
\end{eqnarray*}
and for $k\ge 1$,
\begin{eqnarray*}
   S_{g_1,g_1}^-&\subset& S_{g_{2k},g_{2k}}^+ \\
   S_{g_1,g_1}^-&\subset& S_{g_{2k+1},g_{2k+1}}^+
\end{eqnarray*}
By this and a symmetry argument, we have:
\begin{eqnarray*}
   \rho[g_1(X),g_1(Y)]=1 \Rightarrow
   \left\{
   \begin{array}{ll}
      \rho[g_{i}(X),g_{i}(Y)]=1 & \forall i>1 \\
      \rho[g_{i}(X),g_{j}(Y)]=0 & \mbox{otherwise}\\
   \end{array}
   \right.
\end{eqnarray*}
\begin{eqnarray*}
   \rho[g_2(X),g_2(Y)]=1 \Rightarrow
   \left\{
   \begin{array}{ll}
      \rho[g_{2i}(X),g_{2i}(Y)]=1 & \forall i>1 \\
      \rho[g_{2i}(X),g_{2j}(Y)]=0 & \mbox{otherwise}\\
   \end{array}
   \right.
\end{eqnarray*}
and
\begin{eqnarray*}
   \rho[g_1(X),g_1(Y)]=-1 \Rightarrow
   \left\{
   \begin{array}{ll}
      \rho[g_{i}(X),g_{i}(Y)]=(-1)^i & \forall i>1 \\
      \rho[g_{i}(X),g_{j}(Y)] = 0    & \mbox{otherwise}\\
   \end{array}
   \right.
\end{eqnarray*}
Furthermore, $\rho[g_1(X),g_1(Y)]=1$ iff $\rho^*(X,Y)=1$.


\end{example}




\begin{example}
Suppose $X$ and $Y$ have exponential distributions. Let $\{g_{k}\}$
be the eigenfunctions belonging to $h_K$ with $K$ the exponential
distribution function. For $k$ and $l$ ranging from $1$ to $5$,
$S^+_{g,h}$ and $S^-_{g,h}$ are plotted in Figure~.... It can be
seen that it is possible that $\rho[g_{k}(X),g_l(Y)]=1$ {\em only}
for $k=l$, and $\rho[g_{k}(X),g_l(Y)]=-1$ is not possible for any
$k$ and $l$. Note that for $k>1$
\begin{eqnarray*}
   S_{g_1,g_1}^+ &\subset& S_{g_{k},g_{k}}^+
\end{eqnarray*}
By this, we have:
\begin{eqnarray*}
   \rho[g_1(X),g_1(Y)]=1 \Rightarrow
%   \left\{
%   \begin{array}{ll}
      \rho[g_{i}(X),g_{i}(Y)]=1 \forall i>1
%      \rho[g_{i}(X),g_{j}(Y)]=0 & \mbox{otherwise}\\
%   \end{array}
%   \right.
\end{eqnarray*}
Furthermore, $\rho[g_1(X),g_1(Y)]=1$ iff $\rho^*(X,Y)=1$.
\end{example}




































































































We have introduced a new type of correlation $\rho^*$ which is
closely related to the squared ordinary correlation, having the
range $[0,1]$ and equalling one iff the variables are linearly
related, but equalling zero iff the two variables are independent.
Two interpretations of $\rho^*$ were given: one as a weighted
average of phi-square coefficients, and one as a weighted sum of
squared correlations between eigenfunctions. The decompositions
depend on the marginal distributions, and for certain marginals,
namely the uniform, logistic and exponential, the decomposition was
given explicitly. For general marginals, the closed form solutions
for the decomposition are difficult or impossible to derive, and it
was shown how accurate numerical approximations can be derived. The
decomposition applies to both discrete and continuous random
variables. It was shown that $\kappa$ can be written as a weighted
average of $\phi$-square coefficients. Fr{\'e}chet-type bounds were
given for the component correlations of $\rho^*$....


The U- and V-statistic estimators of $\kappa$ were derived, which
can be used as test statistics for the independence hypothesis.
The asymptotic distribution of the estimators under independence
were derived. The new class of tests includes a generalization of
the Cram{\'e}r-von Mises tests to $K$ ordered samples. Associated
rank tests, which coincide with a class of tests discussed by
\citeA{dewet80}, were obtained as well.































Let
\begin{eqnarray*}
   q(a,b,c,d) = |a-b|+|c-d|-|a-c|-|b-d|
\end{eqnarray*}

Suppose we have $n$ iid data points $(x_i,y_i)$. Let
\begin{eqnarray*}
   \hat F_{1}(x) &=& \frac{1}{n}\sum_{i=1}^n\gamma(x_i,x) + \half \gamma(x_i=x) \\
   \hat F_{2}(y) &=& \frac{1}{n}\sum_{j=1}^n\gamma(y_j,y) + \half \gamma(y_j=y) \\
\end{eqnarray*}
and
\begin{eqnarray*}
   \lefteqn{\hat F_{12}(x,y)
   =
   \frac{1}{n^2}\sum_{i=1}^n\sum_{j=1}^n \gamma(x_i<x,y_i<y) + }\\
   &&\half \gamma(x_i<x,y_i=y) + \half \gamma(x_i=x,y_i<y) + \fourth \gamma(x_i=x,y_i=y)
\end{eqnarray*}
Then the V-statistic estimator can be written as
\begin{eqnarray*}
   \hat\kappa(X,Y)
   &=&
   \frac{2}{n^2}\sum_{i=1}^n\sum_{j=1}^n x_iy_j[\hat F_{12}(x_i,y_j)-\hat F_1(x_i)\hat F_2(y_j)]
   -
   \\
   &&
   \frac{2}{n}\sum_{i=1}^n x_iy_i[\hat F_{12}(x_i,y_i)-\hat F_1(x_i)\hat F_2(y_i)]
\end{eqnarray*}
The computational complexity is $O(n^2)$.


\begin{lemma}\label{lem6}
Let $(X_1,Y_1),\ldots,(X_4,Y_4)$ be iid random variables
distributed as $(X,Y)$. Then
\begin{eqnarray*}
   \kappa(X,Y)
   &=& \int[F_{12}(x,y)-F_1(x)F_2(y)]^2dx dy\\
   &=& 2\int
   xy[F_{12}(x,y)-F_1(x)F_2(y)]d[F_{12}(x,y)-F_1(x)F_2(y)]\\
   &&
   uEX_1Y_2\sign(X_1-X_2)(Y_1-Y_2)]F_1(X_1)F_2(Y_2)S(X_1,Y_2)
\end{eqnarray*}
\end{lemma}


\begin{example}
Consider the trichotomous case that $P(Z=i)=p_i$ ($i=0,1,2$). The
following functions $g_1$ and $g_2$ are orthonormal
\begin{eqnarray*}
   g_1(0) &=& \frac{-p_3}{\sqrt{p_1p_3(1-p_2)}} \\
   g_1(1) &=& 0 \\
   g_1(2) &=& \frac{p_1}{\sqrt{p_1p_3(1-p_2)}}
\end{eqnarray*}
\begin{eqnarray*}
   g_2(0) &=& \sqrt{\frac{p_2}{1-p_2}} \\
   g_2(1) &=& -\sqrt{\frac{1-p_2}{p_2}} \\
   g_2(2) &=& \sqrt{\frac{p_2}{1-p_2}}
\end{eqnarray*}
\end{example}


Denote the population value of $\kappa$ by $\kappa_0$ and write
\begin{eqnarray*}
  r(x,y) = Eh_{F_1}(x,X)h_{F_2}(y,Y)
\end{eqnarray*}
By standard arguments, if $0<\kappa<1$, the asymptotic
distribution of $\sqrt{n}(\hat\kappa(X,Y)-\kappa_0)$ is normal
with mean zero and variance
\begin{eqnarray*}
  \zeta_1 = E[r(X,Y)^2] - [Er(X,Y)]^2
\end{eqnarray*}

\begin{lemma}\label{lem6a}
Suppose $P(Y=0)=1-P(Y=1)=p$. Then we have
\begin{eqnarray*}
   \rho^*(X,Y)
   &=&
   \frac{p(1-p)}{\sqrt{\sum\lambda_k^2}}
   E
   \left(
   |U_1-V_2|+ |U_2-V_1| - |U_1-U_2| - |V_1-V_2|
   \right)\\
   &=&
   \frac{p(1-p)}{\sqrt{\sum\lambda_k^2}}
   \int [G_1(z)-G_2(z)]^2dz
\end{eqnarray*}
\end{lemma}
\proof{Proof} We have
\begin{eqnarray*}
   \kappa(X,Y)
   &=&
   Eh_{F_1}(X_1,X_2)h_{F_2}(Y_1,Y_2)\\
   &=&
   p^2   Eh_{F_1}(U_1,U_2)h_{F_2}(0,0) +
   p(1-p)Eh_{F_1}(U_1,V_2)h_{F_2}(0,1) +\\
   &&
   p(1-p)Eh_{F_1}(V_1,U_2)h_{F_2}(1,0) +
   (1-p)^2Eh_{F_1}(V_1,V_2)h_{F_2}(1,1)
   \\
   &=&
   p^2(1-p)^2
   E
   \left(
   h_{F_1}(U_1,U_2) +
   h_{F_1}(U_1,V_2) +
   h_{F_1}(V_1,U_2) +
   h_{F_1}(V_1,V_2)
   \right)
   \\
   &=&
   p^2(1-p)^2
   E
   \left(
   |U_1-V_2|+ |U_2-V_1| - |U_1-U_2| - |V_1-V_2|
   \right)
\end{eqnarray*}
Since $\mu_1=p(1-p)$ we obtain the first part.

The second part is analogous to....
 \hspace*{\fill}$\Box$\endproof


\begin{lemma}\label{lem6a}
Suppose $P(Y=c_i)=p_i$ with $\sum p_i=1$ and let
\begin{eqnarray*}
    w_{ij}=\frac{p_ip_jh_{F_2}(c_i,c_j)}{\sqrt{\sum\lambda_k^2\sum\mu_k^2}}
\end{eqnarray*}
Then we have
\begin{eqnarray*}
   \rho^*(X,Y)
   &=&
   \sum_{i,j}w_{ij}
   \left(
   |U_{i1}-U_{j2}|+ |U_{i2}-U_{j1}| - |U_{i1}-U_{i2}| - |U_{j1}-U_{j2}|
   \right)\\
   &=&
   \sum_{i,j}w_{ij}
   \left(
   [G_i(z)-G_j(z)]^2dz
   \right)
\end{eqnarray*}
\end{lemma}
\proof{Proof} We have
\begin{eqnarray*}
   \kappa(X,Y)
   &=&
   Eh_{F_1}(X_1,X_2)h_{F_2}(Y_1,Y_2)\\
   &=&
   \sum_{i,j}
   p_ip_j   Eh_{F_1}(U_{i1},U_{j2})h_{F_2}(c_i,c_j)  \\
   &=&
   \sum_{i,j}w_{ij}
   \left(
   |U_{i1}-U_{j2}|+ |U_{i2}-U_{j1}| - |U_{i1}-U_{i2}| - |U_{j1}-U_{j2}|
   \right)\\
   &=&
   \sum_{i,j}w_{ij}
   \left(
   [G_i(z)-G_j(z)]^2dz
   \right)
\end{eqnarray*}
 \hspace*{\fill}$\Box$\endproof



Let $F(z)=pG_1(z)+(1-p)G_2(z)$. The class of Cram{\'e}r-von Mises
tests is based on the population quantity
\begin{eqnarray*}
   \phi(U,V) =
   \int [G_1(z)-G_2(z)]^2w(F(z))d F(z)
\end{eqnarray*}
where $w$ is a nonnegative weight function. The following Lemma
now follows from the previous:
\begin{lemma}\label{lem6}
Suppose $P(Y=0)=1-P(Y=1)=p$. Then with $w(z)=1/k(K^{-1}(z))$ we
have
\begin{eqnarray*}
   \rho_{F\circ K^{-1},I}^*(X,Y) =
   \frac{p(1-p)}{\sqrt{\sum\lambda_k^2}}
   \phi(U,V)
\end{eqnarray*}
where $I$ is the identity function $\gamma(x)=x$.
\end{lemma}
\proof{Proof}

 \hspace*{\fill}$\Box$\endproof

\begin{eqnarray*}
   \lambda_k
   &=&
   \lambda_k \int g_k(z_1)^2d F(z_1) \\
   &=&
   \int [\lambda_k g_k(z_1)]g_k(z_1)d F(z_1) \\
   &=&
   \lefteqn{
   \int_{-\infty}^{\infty} h_F(z_1,z_2)g(z_1)g(z_2)dF(z_1)dF(z_2)  }\\
   &=&
   \int_{-\infty}^{\infty}\left(\int_{-\infty}^{\infty}[\gamma(z, w)-F(z)]dF(z)\right)^2dw
   \ge 0
\end{eqnarray*}

A square integrable $h:\bR^2\rightarrow\bR$ is called positive
definite if
\begin{eqnarray}
   \int h_F(z_1,z_2)g(z_1)g(z_2)dF(z_1)dF(z_2)\ge 0\label{pd1}
\end{eqnarray}
for all measurable functions $g$. It is well known
that~(\ref{pd1}) holds if and only if  $h_F$ has nonnegative
eigenvalues.


For the Laplace distribution we obtain
\begin{eqnarray*}
   f(F^{-1}(u)) = |u-\half|
\end{eqnarray*}
Then the eigenfunctions are the normalized versions of
$\sign(u-\half)J_0\left(\alpha_{1k}\sqrt{2|u-\half|}\right)$ and
the eigenvalues are
\begin{eqnarray*}
   \lambda_k =
   \frac{\alpha_{1k}}{\sigma[J_0(\alpha_{1k}\sqrt{2|U-\half|})]}
\end{eqnarray*}
The sum of eigenvalues is conjectured to be 0.75, the first

Without loss of generality, in this section we only consider
$\rho^*$ and do not consider $\rho^*_{K_1,K_2}$ separately.


Note that ... implies the coefficients $\kappa$ and $\rho^*$
behave similarly to the covariance and the correlation.


We can rewrite $\kappa(X,Y)$ is as
\begin{eqnarray*}
   \kappa(X,Y)
   = \int \frac{[F_{12}(x,y)-F_1(x)F_2(y)]^2}{f_1(x)f_2(y)}dF_1(x)dF_2(y)
\end{eqnarray*}
With
\begin{eqnarray*}
   G_{12}(u,v) = F_{12}(F_1^{-1}(u),F_2^{-1}(v))
\end{eqnarray*}
we can rewrite
\begin{eqnarray*}
   \kappa(X,Y) = \int_{[0,1]^2} \frac{[G_{12}(u,v)-uv]^2}{f_1[F_1^{-1}(u)]f_2[F_2^{-1}(v)]}dudv
\end{eqnarray*}


In the remaining part of this subsection we give an interpretation
of the variance of $h(Z_1,Z_2)$ (Lemma~\ref{lem5}) and show that the
eigenvalues are nonnegative (Lemma~\ref{lem8}).


\begin{lemma}\label{lem5}
Let $Z_{i:n}$ be the $i$th order statistic of an iid sample of
size $n$. Then
\begin{eqnarray*}
   Eh(Z_1,Z_2)^2 = \sum_{k=0}^{\infty}\lambda_k^2 = \frac{1}{6}E(Z_{2:4}-Z_{3:4})^2
\end{eqnarray*}
\end{lemma}
\proof{Proof} The first equality follows from Parseval's theorem. To
prove the second one, note that $|z_1-z_2| + |z_3-z_4| - |z_1-z_3| -
|z_2-z_4|$ equals 0 if $z_1,z_4\le z_2,z_3$ or $z_1,z_4\ge z_2,z_3$,
and equals $\pm(z_{2:4}-z_{3:4})$ otherwise, so that
\begin{eqnarray*}
   E \left(
   |Z_1-Z_2| + |Z_3-Z_4| - |Z_1-Z_3| - |Z_2-Z_4|
   \right)^2
   = \frac23E(Z_{2:4}-Z_{3:4})^2
\end{eqnarray*}
Write $a_{ij}=|Z_i-Z_j|$. Then we obtain
\begin{eqnarray*}
   Eh(Z_1,Z_2)^2
   &=&   \frac14
   \left(a_{12}-a_{13}-a_{24}+a_{34}\right)
   \left(a_{12}-a_{15}-a_{26}+a_{56}\right)
   \\
   &=& ...\\
   &=& \frac14 E
   \left(
   a_{12}-a_{13}-a_{24}+a_{34}
   \right)^2\\
   &=& \frac16E(Z_{2:4}-Z_{3:4})^2
\end{eqnarray*}
 \hspace*{\fill}$\Box$\endproof
An important property of the eigenvalues is the following:
\begin{lemma}\label{lem1}
The eigenvalues of $h_F$ are nonnegative.
\end{lemma}
\proof{Proof} Suppose the common distribution function of $Z_1$ and
$Z_2$ is $F$. Note that since
\begin{eqnarray*}
   |z_1-z_2| = \int_{-\infty}^{\infty}[\gamma(z_1, w)-\gamma(z_2, w)]^2dw
\end{eqnarray*}
we can write
\begin{eqnarray}
   h(z_1,z_2)
   &=& -\half E\int_{-\infty}^{\infty}
   \left(
   [\gamma(z_1, w)-\gamma(z_2, w)]^2
   -
   [\gamma(z_1, w)-\gamma(Z_2, w)]^2-\right.
   \nonumber\\
   &&
   \left.
   [\gamma(Z_1, w)-\gamma(z_2, w)]^2
   -
   [\gamma(Z_1, w)-\gamma(Z_2, w)]^2
   \right)
   dw
   \nonumber\\
   &=&
   E\int_{-\infty}^{\infty}
   [\gamma(z_1, w)-\gamma(Z_1, w)]
   [\gamma(z_2, w)-\gamma(Z_2, w)]
   dw
   \nonumber\\
   &=&
   \int_{-\infty}^{\infty}[\gamma(z_1, w)-F(w)][\gamma(z_2, w)-F(w)]dw
   \label{h5}
\end{eqnarray}
Then we obtain by~(\ref{g1}),~(\ref{g2}) and~(\ref{h5})
\begin{eqnarray*}
   \lambda_k
   &=&
   \lambda_k Eg_k^2(Z_1) \\
   &=&
   Eg_k(Z_1)g_k(Z_2)h(Z_1,Z_2) \\
   &=&
   E\int_{-\infty}^{\infty}
   g(Z_1)[\gamma(Z_1, w)-F(w)]
   g(Z_2)[\gamma(Z_1, w)-F(w)]
   dw\\
   &=&
   \int_{-\infty}^{\infty}
   \left(
   Eg(Z)[\gamma(Z, w)-F(w)]
   \right)^2
   dw \ge 0\\
\end{eqnarray*}
 \hspace*{\fill}$\Box$\endproof


If $F_1=F_2$, the support of the Fr{\'e}chet upper bound is the
straight line $\{(F_1^{-1}(x),F_2^{-1}(x)):x\in\bR\}$ and the
support of the Fr{\'e}chet lower bound is the straight line
$\{(F_1^{-1}(x),-F_2^{-1}(x)):x\in\bR\}$. The concept can be
extended to bounds for $g(X)$ and $h(Y)$ such that the correlation
correlation between $g(X)$ and $h_FY)$ is maximal. We now obtain
bounds for the individual correlations. Let
\begin{eqnarray*}
   S^+_{kl}(X,Y)
   &=&
   \left\{
   (x,y)\in\bR^2|g_k(x)=g_{2l}(y)
   \right\}\\
   S^-_{kl}(X,Y)
   &=&
   \left\{
   (x,y)\in\bR^2|g_k(x)=- g_{2l}(y)
   \right\}
\end{eqnarray*}
Then it is well-known that $\rho(g_k(X),g_{2l}(Y))=1$ iff the
support of $(X,Y)$ is $S^+_{kl}(X,Y)$ and
$\rho(g_k(X),g_{2l}(Y))=-1$ iff the support of $(X,Y)$ is
$S^-_{kl}(X,Y)$. The distributions with supports $S^+_{kl}(X,Y)$ and
$S^-_{kl}(X,Y)$ are the so-called Fr{\'e}chet bounds for the
bivariate random variable $(g_k(X),g_{2l}(Y))$.

For some specific cases we next give the solution.

\begin{example}\label{ex uni}
For the uniform distribution, we obtain the functions
$g_k(x)=\cos(\pi k x)/2$ as solutions, and eigenvalues
$\lambda_k=(\pi k)^{-2}$. In this case, the Equations~(\ref{de1})
and~(\ref{de2}) are of the same form, and so $q_k=g_k$.
\end{example}

\begin{example}\label{ex log}
For the logistic distribution with distribution function $F(x) =
\frac{1}{1+e^{-x}}$ and density $f(x) =
\frac{1}{(1+e^{-x})(1+e^{x})}$ we obtain
\begin{eqnarray*}
   f(F^{-1}(u)) = u(1-u)
\end{eqnarray*}
With $P_k$ the $k$th Legendre polynomial, the solutions are given
by $q_k(u)=\sqrt{2kU}\,P_k(2u-1)$ and
$g_k(x)=\sqrt{2kU}\,P_k[2F(x)-1]$ with eigenvalues
$\lambda_k=\frac{1}{k(kU)}$.
\end{example}


\begin{example}\label{ex exp}
For the exponential distribution we obtain
\begin{eqnarray*}
   f(F^{-1}(u)) = 1-u
\end{eqnarray*}
Let $J_i$ be the $i$th Bessel function of the first kind, let
$\alpha_{k}$ be the $k$th zero of $J_1$ and let
$\beta_k=(J_0(\alpha_k)^2+J_1(\alpha_k)^2)^{-1/2}$. The
eigenfunctions are $q_k(u)=\beta_kJ_0(\alpha_{k}\sqrt{1-u})$ and
the eigenvalues are $\lambda_k = \frac{4}{\alpha_k^2}$. Numerical
evidence indicates that the sum of the eigenvalues is $1/2$. The
first is equal to 0.27245.
\end{example}



\begin{abstract} \noindent
A possible drawback of the ordinary correlation coefficient $\rho$
for two real random variables $X$ and $Y$ is that $\rho=0$ does not
imply independence. In this paper a correlation coefficient $\rho^*$
is introduced which assumes values between zero and one, equalling
zero iff the two variables are independent and equalling one iff the
two variables are linearly related. For a $2\times 2$ contingency
table, $\rho^*$ and the squared ordinary correlation coincide. Two
interpretations of $\rho^*$ are given. First, it is shown that
$\rho^*$ is a weighted sum of squared correlations between marginal
eigenfunctions. For uniform, logistic and exponential marginals
closed form expressions for the weights and eigenfunctions are given
and for other marginals it is shown how numerical approximations can
be obtained. Fr{\'e}chet bounds are discussed for the component
correlations. Second, it is demonstrated that $\rho^*$ can be viewed
as a weighted average of phi-square coefficients. The weights are
plotted for several marginal distributions, indicating how the tails
of the distributions are weighted. Estimators of $\rho^*$ based on
both U- and V-statistics are discussed and their asymptotic
distributions under independence, which depend on the marginal
distributions of the variables, are derived. It is shown that as a
special case a new Cram{\'e}r-von Mises test for $K$ ordered samples
is obtained, which reduces to the ordinary Cram{\'e}r-von Mises test
for two samples. The paper includes an extensive consideration of
grade versions of $\rho^*$, and discrete distributions are given
ample attention.
\end{abstract}


In fact, $\kappa$ is quite closely related the the covariance and
hence $\rho^*$ is quite closely related to the ordinary
correlation. This can be seen as follows. Write $a_{ij}=|X_i-X_j|$
and $b_{ij}=|Y_i-Y_j|$ so that
\begin{eqnarray*}
   \kappa(X,Y) = E(a_{12}-a_{13})(b_{12}-b_{24})
\end{eqnarray*}
Now the squared covariance is
\begin{eqnarray*}
   \cov(X,Y)^2
   &=& \fourth [E(X_1-X_2)(Y_1-Y_2)]^2 \\
   &=& \fourth E(X_1-X_4)(X_2-X_3)(Y_1-Y_4)(Y_2-Y_3) \\
   &=& \sixteenth E\left(a_{12}^2-a_{13}^2-a_{24}^2+a_{34}^2\right)
                   \left(b_{12}^2-b_{13}^2-b_{24}^2+b _{34}^2\right) \\
   &=& \sixteenth E\left(a_{12}^2-a_{13}^2\right)\left(b_{12}^2-b_{24}^2\right)+
       \sixteenth E\left(a_{12}^2-a_{13}^2\right)\left(b_{34}^2-b_{13}^2\right)+\\
   &&  \sixteenth E\left(a_{34}^2-a_{24}^2\right)\left(b_{12}^2-b_{24}^2\right)+
       \sixteenth E\left(a_{34}^2-a_{24}^2\right)\left(b_{34}^2-b_{13}^2\right) \\
   &=& E\left(\half a_{12}^2-\half a_{13}^2\right)
       \left(\half b_{12}^2-\half b_{24}^2\right)
\end{eqnarray*}
i.e., replacing the absolute differences in $\kappa$ by half the
squared differences yields the squared covariance.



Part 3: The inequality is of the Cauchy-Schwartz type so it holds.
If $Y=aX+b$ a.s.\ for certain constants $a$ and $b$ then by Part~1
$\kappa(Y,Y)=a^2\kappa(X,X)$ and $\kappa(X,Y)=a\kappa(X,X)$ a.s.
The equality holds iff
\begin{eqnarray*}
   h_{F_1}(X_1,X_2)=ah_{F_2}(Y_1,Y_2)+b \mbox{ a.s.}
\end{eqnarray*}
for certain constants $a$ and $b$. Taking expectations it follows
that $b=0$ so we have
\begin{eqnarray}
   h_{F_1}(X_1,X_2)=ah_{F_2}(Y_1,Y_2) \mbox{ a.s.}
   \label{cauchy}
\end{eqnarray}
Hence
\begin{eqnarray*}
   Eh_{F_1}(X_1,X_2)h_{F_2}(Y_1,Y_2) = aEh_{F_2}(Y_1,Y_2)^2\mbox{ a.s.}
\end{eqnarray*}
This can be rewritten as
\begin{eqnarray*}
   \int [\gamma(X_1,w)-\gamma(X_2,w)]
\end{eqnarray*}


and $Y$ are linearly related it is immediate that equality holds.
Now suppose
\begin{eqnarray}
   h_{F_1}(z_1,z_2) = h_{F_2}(z_1,z_2) \forall z_1,z_2
   \label{heq1}
\end{eqnarray}
We need to show that $X$ and $Y$ are a.s.\ linearly related. Let
$\tilde Y=aY+b$ where the constants $a$ and $b$ are chosen such that
\begin{eqnarray*}
   E|X_1-X_2| &=& E|\tilde Y_1-\tilde Y_2| \\
   E|X| &=& E|\tilde Y|
\end{eqnarray*}
Then substituting $z_1=z$ and $z_2=0$ in~(\ref{heq1}) we obtain
\begin{eqnarray*}
   E|z-X| = E|z-\tilde Y| \forall z
\end{eqnarray*}
This becomes
\begin{eqnarray*}
   E\int [\gamma(z,w)-\gamma(X,w)]^2dw =
   E\int [\gamma(z,w)-\gamma(\tilde Y,w)]^2dw\forall z
\end{eqnarray*}
so that
\begin{eqnarray*}
   \int [\gamma(z,w) + F_1(w) -2 \gamma(z,w)F_1(w)]dw =
   \int [\gamma(z,w) + \tilde F_2(w) -2 \gamma(z,w)\tilde F_2(w)]dw \forall z
\end{eqnarray*}
Some algebra now yields
\begin{eqnarray*}
   \int [F_1(w)-tilde F_2(w)](1-2\gamma(z,w)dw = 0
\end{eqnarray*}
which evaluates to
\begin{eqnarray*}
   \int[F_1(w)-tilde F_2(w)]\sign(z-w)dw = 0
\end{eqnarray*}



Now suppose equality holds. We need to show that  $Y=aX+b$ a.s.\ for
certain constants $a$ and $b$. Denote the Fr{\'e}chet upper and
lower bounds as
\begin{eqnarray*}
   F_{12}^+(x,y) &=& \min\{F_1(x),F_2(y)\}\\
   F_{12}^-(x,y) &=& \max\{0,1-F_1(x)-F_2(y)\}
\end{eqnarray*}
for which the well-known relation holds
\begin{eqnarray*}
   F_{12}^-(x,y) \le F_{12}(x,y) \le F_{12}^+(x,y)
\end{eqnarray*}
Now we have either
\begin{eqnarray*}
   \int [F_{12}(x,y)-F_1(x)F_2(y)]^2dxdy \le
   \int [F_{12}^+(x,y)-F_1(x)F_2(y)]^2dxdy
\end{eqnarray*}
or
\begin{eqnarray*}
   \int [F_{12}(x,y)-F_1(x)F_2(y)]^2dxdy \le
   \int [F_{12}^-(x,y)-F_1(x)F_2(y)]^2dxdy
\end{eqnarray*}
or both. Let's continue with the first, the second is analogous.
First consider the continuous case. Then there is a strictly
increasing function $u$ such that
\begin{eqnarray*}
   F_1(x)=F_2[u(x)]\forall x
\end{eqnarray*}
Write
\begin{eqnarray*}
   M_i(x,y) = \min\{F_i(x),F_i(y)\}-F_i(x)F_i(y)
\end{eqnarray*}
Then from~(\ref{fxy}) and a change of variables  $y\rightarrow
u(y)$ we obtain
\begin{eqnarray*}
   \kappa(X,Y)
   &\le& \int [F_{12}(x,y)-F_1(x)F_2(y)]^2dxdy\\
   &=& \int [F_{12}^-(x,u(y))-F_1(x)F_2[u(y)]]^2dxdu(y) \\
   &=& \int [\min\{F_1(x),F_1(y)-F_1(x)F_1(y)]^2\frac{du(y)}{dy}dx dy\\
   &=& \int M_1(x,y)^2\frac{du(y)}{dy}dx dy
\end{eqnarray*}
Similarly,
\begin{eqnarray*}
   \kappa(X,X)
   = \int M_1(x,y)^2dxdy
\end{eqnarray*}
and
\begin{eqnarray*}
   \kappa(Y,Y)
   &=& \int [\min\{F_2(x),F_2(y)\}-F_2(x)F_2(y)]^2dxdy \\
   &=& \int [\min\{F_2(u(x)),F_2(u(y))\}-F_2(u(x))F_2(u(y))]^2du(x)du(y) \\
   &=& \int [\min\{F_1(x),F_1(y)\}-F_1(x)F_1(y)]^2\frac{du(x)}{dx}\frac{du(y)}{dy}dxdy \\
   &=& \int M_1(x,y)^2\frac{du(x)}{dx}\frac{du(y)}{dy}dxdy
\end{eqnarray*}
Then by Cauchy-Schwartz we have
\begin{eqnarray*}
   \left(\int M_1(x,y)^2\frac{du(y)}{dy}dx dy\right)^2
   &\le& \int M_1(x,y)^2\frac{du(x)}{dx}\frac{du(y)}{dy}dxdy
\end{eqnarray*}


Now by the Cauchy-Schwartz inequality,
\begin{eqnarray*}
   \kappa(X,Y)^2 = \kappa(X,X)\kappa(Y,Y)
\end{eqnarray*}
iff for some constant $a$,
\begin{eqnarray*}
   \frac{du(x)}{dx} = a
\end{eqnarray*}
Hence, $u(x)=ax+b$. It follows that $X$ and $Y$ are a.s.\ linearly
related.


Part 4: By Part~1 it is sufficient to consider the case that $X$
and $Y$ take values in $\{0,1\}$. Write $p_{ij}=P(X=i,Y=j)$,
$p_{i+}=P(X=i)$ and $p_{+j}=P(Y=j)$. Then
\begin{eqnarray*}
   \kappa(X,X)
   &=& Eh_{F_1}(X_1,X_2)^2\\
   &=& ... = \var(X)
\end{eqnarray*}
and
\begin{eqnarray*}
   \kappa(X,Y)
   &=& Eh_{F_1}(X_1,X_2)h_{F_2}(Y_1,Y_2)\\
   &=& ... = \cov(X,Y)
\end{eqnarray*}

Let us assume $h_F(Z_1,Z_2)$ has finite variance, i.e., that
\begin{eqnarray*}
   Eh_F(Z_1,Z_2)^2<\infty
\end{eqnarray*}
It is then also said that $h_F$ is {\em square-integrable}. It
follows by linear operator theory that there exists a complete
system of orthonormal functions $\{g_k\}$ (called eigenfunctions)
satisfying
\begin{eqnarray}
   Eg_k(Z)g_l(Z) = \delta_{kl}\label{g1}
\end{eqnarray}
(with $\delta$ the Kronecker delta) and numbers $\{\lambda_k\}$
(called eigenvalues) such that
\begin{eqnarray*}
   h_F(z_1,z_2) = \sum_{k=0}^{\infty}\lambda_kg_k(z_1)g_k(z_2)
\end{eqnarray*}
where convergence is in mean square sense, i.e.,
\begin{eqnarray*}
   E\left[h_F(Z_1,Z_2) - \sum_{k=0}^{N}\lambda_kg_k(Z_1)g_k(Z_2)
   \right]^2 \rightarrow 0
\end{eqnarray*}
as $N\rightarrow\infty$. By degeneracy of $h_F$ we have
\begin{eqnarray}
   Eg_k(Z) = 0
   \label{degen}
\end{eqnarray}


We now show how we may relate the $\rho_{ij}$ to the canonical
correlations. Canonical correlation analysis is based on the
canonical decomposition
\begin{eqnarray*}
   \Omega
   =
   1+
   \sum_{k=1}^\infty
   u_{1k}(x)u_{2k}(y)\rho_{k}
\end{eqnarray*}
where $u_{11}$ and $v_{11}$ are functions maximizing the
correlation $\rho_{1}$ between $X$ and $Y$. Furthermore, $u_{1k}$
is orthogonal to $u_{11},\ldots,u_{1,k-1}$ and $v_{1k}$ are
orthogonal to  is orthogonal to $v_{11},\ldots,v_{1,k-1}$ and they
maximize the correlation $\rho^{(C)}_k$. For $I\times J$
contingency tables, $\rho_k=0$ for $k\ge\min(I,J)$. Now we have
\begin{eqnarray*}
   \phi^2 = \sum_{k=1}^\infty\rho^2_{k}
\end{eqnarray*}


Consider the set of sign functions $\cS$ with elements functions
of the form $s_w(z)=\sign(w-z)$. Then the mean and variance of
these functions are given as where
\begin{eqnarray*}
   \mu_w = E(\sign(Z-w) \\
   \sigma_w^2 = \var[\sign(Z-w)] = P(Z\ne w)-\mu_w^2
\end{eqnarray*}
Then we obtain
\begin{eqnarray*}
   \rho_K = \sup_{s_1,s_2\in\cS}\rho[s_1(X),s_2(Y)] =
   ...
   = \sup_{x,y}[F_{12}(x,y)-F_1(x)F_2(y)]^2
\end{eqnarray*}
where $c=...$. Now this reduces to Kolmogorov test... so $\rho_K$
may be called Kolmogorov maximal correlation.



A quite different class of tests is based on canonical
correlations. Canonical correlation analysis is based on $u_{11}$
and $v_{11}$ are functions maximizing the correlation $\rho_{1}$
between $X$ and $Y$. Furthermore, $u_{1k}$ is orthogonal to
$u_{11},\ldots,u_{1,k-1}$ and $v_{1k}$ are orthogonal to  is
orthogonal to $v_{11},\ldots,v_{1,k-1}$ and they maximize the
correlation $\rho^{(C)}_k$. For $I\times J$ contingency tables,
$\rho_k=0$ for $k\ge\min(I,J)$. Now we can relate this to $\phi^2$
as
\begin{eqnarray*}
   \phi^2 = \sum_{k=1}^\infty\rho^2_{k}
\end{eqnarray*}


From results of \citeA{ks78}, this approach is problematic. In
particular, they showed that for any distribution function
$F_{12}$ with continuous marginals, there is a sequence
$F_{12}^{(k)}$ ($k=1,2,\ldots$) such that
\begin{eqnarray*}
   F_{12}^{(k)} \rightarrow F_{12}
\end{eqnarray*}
for $k\rightarrow\infty$ in mean square, while the canonical
correlation for $F_{12}^{(k)}$ equals 1. For most practical
purposes, this is unacceptable. Note that the canonical
correlation is also called Renyi's maximal correlation
\cite{renyi59}.


As a solution to this problem, \citeA{ks78} proposed to maximize
over {\em monotone} functions, that is,
\begin{eqnarray*}
   \sup_{a,b\mbox{ monotone}} \rho[a(X),b(Y)]
\end{eqnarray*}
We can consider other sets of functions as well. For example, a
subset of the set of monotone functions can be formed using
$\gamma$:
\begin{eqnarray*}
   \rho_\gamma =
   \sup_{x,y}
   \left|\rho[\gamma(X,x),\gamma(Y,y)]\right|
   =
   \sup_{x,y}\phi(x,y)
\end{eqnarray*}
where $\phi$ is defined by~(\ref{phidef1}). Thus, we have that
$\rho_\gamma$ is related $\rho_{K,K}$ with $K$ the logistic
distribution function since then $\rho_{K,K}=E\phi(X,Y)$.


We now show that the sample version of $\rho_\gamma$ is a
generalized Kolmogorov-Smirnov test statistic. If $Y$ is
dichotomous with $P(Y=0)=p$ and $P(Y=1)=1-p$, we obtain the
reduction
\begin{eqnarray*}
   \phi(x,y) =
   \frac
   {\left|F_{1|2}(x|Y=0)-F_{1|2}(x|Y=1)\right|}
   {\sqrt{F_1(x)\,(1-F_1(x))\,p\,(1-p)}}
\end{eqnarray*}
Hence, the sample version of $\rho_\gamma$ then leads to a
Kolmogorov-Smirnov test. Note that the usual Kolmogorov's test is
based on the numerator of $\phi$ which is a maximal covariance:
\begin{eqnarray*}
   \sup_{x,y} \cov[\gamma(X,x),\gamma(Y,y)]
   =\sup_{x,y}\left|F_{1|2}(x|Y=0)-F_{1|2}(x|Y=1)\right|
\end{eqnarray*}


The standard Kolmogorov test for testing whether the two
conditional distribution functions $F_{2|1}(y|X=0)$ and
$F_{2|1}(y|X=1)$ is based on the sample version of
\begin{eqnarray*}
   \upsilon = \sup_{y}\left|F_{2|1}(y|X=0)-F_{2|1}(y|X=1)\right|
\end{eqnarray*}
A natural generalization of this quantity to two arbitrary real
random variables is
\begin{eqnarray*}
   \tau = \sup_{x,y}\left|F_{12}(x,y)-F_1(x)F_2(y)\right|
\end{eqnarray*}
Now if $X$ is dichotomous with $P(X=0)=1-P(X=1)=p$, it can be
verified that
\begin{eqnarray*}
   \upsilon = p(1-p)\tau
\end{eqnarray*}


In fact, a version of Kolmogorov's test statistic which is
slightly modified for the discrete case is a maximal correlation.
Note that
\begin{eqnarray*}
   \rho[\gamma(X,x),\gamma(Y,y)]
   =
   \frac
   {F_{12}(x,y) - F_1(x)F_2(y)}
   {P(X\ne x)P(Y\ne y)}
\end{eqnarray*}
Hence, with $\Gamma=\{\gamma(.,z)|z\in\bR\}$, we have
\begin{eqnarray*}
   \tilde\tau =
   \sup_{x,y}
   \frac
   {\left|F_{12}(x,y)-F_1(x)F_2(y)\right|}
   {P(X\ne x)P(Y\ne y)}
   =
   \sup_{a,b\in\Gamma}
   |\rho[a(X),b(Y)]|
\end{eqnarray*}
Note that in the continuous case $\tilde\tau=\tau$.


As a solution to this problem, \citeA{ks78} proposed to maximize
over {\em monotone} functions, that is,
\begin{eqnarray*}
   \sup_{a,b\mbox{ monotone}} \rho[a(X),b(Y)]
\end{eqnarray*}
We can consider other sets of functions as well. For example, a
subset of the set of monotone functions can be formed using
$\gamma$:
\begin{eqnarray*}
   \rho_\gamma =
   \sup_{x,y}
   \left|\rho[\gamma(X,x),\gamma(Y,y)]\right|
   =
   \sup_{x,y}\phi(x,y)
\end{eqnarray*}
where $\phi$ is defined by~(\ref{phidef1}). Thus, we have that
$\rho_\gamma$ is related $\rho_{K,K}$ with $K$ the logistic
distribution function since then $\rho_{K,K}=E\phi(X,Y)$.


We now show that the sample version of $\rho_\gamma$ is a
generalized Kolmogorov-Smirnov test statistic. If $Y$ is
dichotomous with $P(Y=0)=p$ and $P(Y=1)=1-p$, we obtain the
reduction
\begin{eqnarray*}
   \phi(x,y) =
   \frac
   {\left|F_{1|2}(x|Y=0)-F_{1|2}(x|Y=1)\right|}
   {\sqrt{F_1(x)\,(1-F_1(x))\,p\,(1-p)}}
\end{eqnarray*}
Hence, the sample version of $\rho_\gamma$ then leads to a
Kolmogorov-Smirnov test. Note that the usual Kolmogorov's test is
based on the numerator of $\phi$ which is a maximal covariance:
\begin{eqnarray*}
   \sup_{x,y} \cov[\gamma(X,x),\gamma(Y,y)]
   =\sup_{x,y}\left|F_{1|2}(x|Y=0)-F_{1|2}(x|Y=1)\right|
\end{eqnarray*}


Finally, canonical correlation is related as follows... The first
canonical correlation is Renyi's maximal correlation, etc... A
quite different class of tests is based on canonical correlations.
Canonical correlation analysis is based on $u_{11}$ and $v_{11}$
are functions maximizing the correlation $\rho_{1}$ between $X$
and $Y$. Furthermore, $u_{1k}$ is orthogonal to
$u_{11},\ldots,u_{1,k-1}$ and $v_{1k}$ are orthogonal to  is
orthogonal to $v_{11},\ldots,v_{1,k-1}$ and they maximize the
correlation $\rho^{(C)}_k$. For $I\times J$ contingency tables,
$\rho_k=0$ for $k\ge\min(I,J)$. Now we can relate this to $\phi^2$
as
\begin{eqnarray*}
   \phi^2 = \sum_{k=1}^\infty\rho^2_{k}
\end{eqnarray*}
The two approaches thus seem to have similar potential applications.
We summarize a few salient differences between them. Firstly the
canonical correlations are `multidimensional' objects, in the sense
that many forms of associations are summarized by them; in
particular, $\rho_j=0$ iff $X$ and $Y$ are independent (Lancaster
1959, Renyi 1959). In contrast, the $\rho_{kl}$ are one-dimensional.
Secondly, the $\rho_{kl}$ parameterize the association between $X$
and $Y$: knowing the marginal distributions and the $\rho_{kl}$
completely specify the joint distribution. On the other hand, the
canonical correlations do not parameterize the association, in the
sense that knowing the marginal distributions and the canonical
correlations does not specify the joint distribution completely. For
example, consider a $2\times 2$ table with uniform marginals and a
canonical correlation of 1. Then the joint distribution either has
its probability mass uniformly distributed on the main diagonal or
on the other diagonal.



\subsection{Grade versions of $\kappa$}

The results above lead to an interesting interpretation of
$\kappa_{K_1,K_2}$. Below, it is shown that Lemma~\ref{lem4}, Part
3, can be used to write $\kappa$ as an average of $\phi$
coefficients for the $2\times 2$ tables obtained by collapsing the
bivariate distribution with respect to cut-points $(x,y)$. The
$\phi$ coefficient for measuring the dependence in the $2\times 2$
table obtained by collapsing the distribution with respect to the
cut point $(x,y)$ is given as
\begin{eqnarray}
   \phi(x,y) = \frac{\left|F_{12}(x,y)-F_1(x)F_2(y)\right|}{\sqrt{F_1(x)[1-F_1(x)]F_2(y)[1-F_2(y)]}}
   \label{phidef1}
\end{eqnarray}
Suppose $K_1$ and $K_2$ have corresponding densities $k_1$ and
$k_2$, respectively. Then from Lemma~\ref{lem4}, Part 3 we obtain
that $\kappa_{K_1,K_2}$ can be written as a weighted average
$\phi$-square:
\begin{eqnarray*}
   \kappa_{K_1,K_2}(X,Y) = \int\phi(x,y)^2\,w_{K_1}[F_1(x)]\,w_{K_2}[F_2(x)]dF_1(x)dF_2(y)
\end{eqnarray*}
where the weight function $w$ is defined by
\begin{eqnarray*}
   w_K(u) = \frac{u(1-u)}{k[K^{-1}(u)]}
\end{eqnarray*}
The normalized weight function (integrating to one) is
\begin{eqnarray*}
  {\bar w}_K(u) = \frac{w_K(u)}{\int_0^1 w_K(u)du}
\end{eqnarray*}
where
\begin{eqnarray*}
   \int_0^1 w_K(u)du
   = \int_0^1{u(1-u)}dK^{-1}(u)
   = \int_{-\infty}^{\infty}K(x)[1-K(x)]dx
\end{eqnarray*}
In Figure~\ref{ffinv}, $\bar w_K$ is plotted for the distribution
functions $K$ given in Table~\ref{eigtable}.


\begin{figure}
\begin{center}
\includegraphics[width=.48\linewidth,clip]{ggplot.eps}
\end{center}
%\parbox{0.95\textwidth}
{\caption{Plots of $\bar w_K(u)$ for several distributions
$K${\label{ffinv}}}}
\end{figure}


From Figure~\ref{ffinv} we can deduce which $\phi$-coefficients
are given most weight for different marginal distributions. As the
reference marginal, we take the logistic, which is the horizontal
line in the figure, i.e., assigning uniform weights. We see that
for a uniform marginal, the weight goes to zero in the tails. The
weights for a normal marginal are for most $u$ intermediate
between the weights for the uniform and logistic marginals. In
contrast to uniform and normal marginals, the Laplace distribution
gives more weight to the tails than a logistic marginal. An
exponential marginal gives little weight to the lower tail, but
much weight to the upper tail. Finally, the chi-square
distribution gives very large weight to the lower tail and very
small weight to the upper tail. Interestingly, of the
distributions considered the biggest difference is between the
chi-square and the exponential distribution.


We have
\begin{eqnarray}
   |z_1-z_2| = \int_{-\infty}^{\infty}[\gamma(z_1,w)-\gamma(z_2,w)]^2dw
   \label{difint}
\end{eqnarray}
and with $z_{i:4}$ the $i$th largest number in the set
$\{z_1,z_2,z_3,z_4\}$, we have
\begin{eqnarray}
   \lefteqn{
   \int_{-\infty}^{\infty}
   \left|
   [\gamma(z_1, w)-\gamma(z_3, w)]
   [\gamma(z_2, w)-\gamma(z_4, w)]
   \right|
   dw =}\nonumber\hspace*{30mm}\\
   &&
   \left\{
   \begin{array}{ll}
      0                & \mbox{if $z_1,z_3<z_2,z_4$ or $z_1,z_3>z_2,z_4$}  \\
      z_{3:4}-z_{2:4}  & \mbox{otherwise}
   \end{array}
   \right.
   \label{fubini1}
\end{eqnarray}
Now we can derive the desired result first using~(\ref{difint}),
then applying Fubini's theorem which is justified because the right
hand side of~(\ref{fubini1}) can be seen to be finite, and finally
using~(\ref{gamma1}):
\begin{eqnarray}
   h_F(z_1,z_2)
   &=& -\half E\int_{-\infty}^{\infty}
   \left(
   [\gamma(z_1, w)-\gamma(z_2, w)]^2
   -
   [\gamma(z_1, w)-\gamma(Z_2, w)]^2-\right.
   \nonumber\\
   &&
   \left.
   [\gamma(Z_1, w)-\gamma(z_2, w)]^2
   -
   [\gamma(Z_1, w)-\gamma(Z_2, w)]^2
   \right)
   dw
   \nonumber\\
   &=&
   E\int_{-\infty}^{\infty}
   [\gamma(z_1, w)-\gamma(Z_1, w)]
   [\gamma(z_2, w)-\gamma(Z_2, w)]
   dw
   \nonumber\\
   &=&
   \int_{-\infty}^{\infty}[\gamma(z_1, w)-F(w)][\gamma(z_2, w)-F(w)]dw\nonumber
\end{eqnarray}


\subsection{Regression}

The most commonly used approach for detecting association is
regression. This has the advantage that the interval nature of the
data can be used. There are two disadvantages compared to the
present approach: (i) only differences in the conditional means can
be detected, not the entire conditional distributions and (ii) one
has to choose the regression functions arbitrarily. Of course the
latter could be avoided using nonparametric regression, but in that
case testing becomes more difficult.

Compared with regression, we should say that it aims at measuring
something different than what we aim at. Namely, in regression the
conditional means (or something else such as conditional medians)
are of interest, while the conditional distributions apart from
their means are of secondary interest.


For two random variables $X$ and $Y$ with joint distribution
function $F_{12}$ and marginal distribution functions $F_1$ and
$F_2$, the Fr{\'e}chet upper bound $F_{12}^+$ is defined by
\begin{eqnarray*}
   F_{12}^+(x,y) =\min\{F_1(x),F_2(y)\}
\end{eqnarray*}
and the Fr{\'e}chet lower bound $F_{12}^-$ is defined by
\begin{eqnarray*}
   F_{12}^-(x,y) = \max\{0,1-F_1(x)-F_2(y)\}
\end{eqnarray*}
Then $\rho(X,Y)=1$ if and only if $F_{12}=F_{12}^+$ and
$\rho(X,Y)=-1$ if and only if $F_{12}=F_{12}^-$. For functions $g$
and $h$, a related question is for which $F_{12}$ the correlation
between $g(X)$ and $h(Y)$ is maximal or minimal. Suppose for
simplicity that $g$ and $h$ are standardized, then
$\rho[g(X),h(Y)]=1$ iff $P[g(X)=h(Y)]=1$ and $\rho[g(X),h(Y)]=-1$
iff $P[g(X)=-h(Y)]=1$. Note that there are $g$ and $h$ such that for
no $X$ and $Y$ these bounds can be attained. Now it is easy to
verify that $\rho[g(X),h(Y)]=1$ iff the support of the distribution
of $(X,Y)$ is given by
\begin{eqnarray*}
   S^+_{g,h}
   &=&
   \left\{
   (x,y)\in\bR^2|g(x)=h(y)
   \right\}
\end{eqnarray*}
and $\rho[g(X),h(Y)]=-1$ iff the support is given by
\begin{eqnarray*}
   S^-_{g,h}
   &=&
   \left\{
   (x,y)\in\bR^2|g(x)=- h(y)
   \right\}
\end{eqnarray*}
It follows that the bound $\rho[g(X),h(Y)]=1$ can be attained iff
for almost all $x$ there is a $y$ such that $(x,y)\in S^+_{g,h}$ and
vice versa. Similarly, the bound $\rho[g(X),h(Y)]=-1$ can be
attained if and only if for almost all $x$ there is a $y$ such that
$(x,y)\in S^-_{g,h}$ and vice versa. Furthermore, if for functions
$a$ and $b$, $S^+_{a,b}\subseteq S^+_{g,h}$, then
$\rho[a(X),b(Y)]=1$ implies $\rho[g(X),h(Y)]=1$.


Finiteness of $h_F(0,0)$ follows immediately by definition from
existence of $h_F$. Suppose now $h_F(0,0)$ is finite. Then
\begin{eqnarray*}
   \lefteqn{h_F(0,0) - \half(|z_1-z_2|-|z_1|-|z_2|)}\\
   &=&
   -\half E
   \left(
      |z_1-z_2| - |z_1| - |Z_2| - |Z_1| - |z_2| + |Z_1-Z_2|
   \right)\\
   &\ge&
   -\half E
   \left(
      |z_1-z_2| - |z_1-Z_2| - |Z_1-z_2| + |Z_1-Z_2|
   \right)\\
   &=& h_F(z_1,z_2)
\end{eqnarray*}
so $h_F$ is bounded above.

\begin{enumerate}

\item $\sum_{k=0}^{\infty}\lambda_k<\infty$ iff $EZ$ exists. Then $\sum_{k=0}^{\infty}\lambda_k=Eh_F(Z,Z)$

%$\sum_{k=0}^{\infty}\lambda_k=\int_0^1\frac{u(1-u)}{f[F^{-1}(u)]}du$

\item $\sum_{k=0}^{\infty}\lambda_k^2 = Eh_F(Z_1,Z_2)^2$.

\item $Eh_F(Z,Z)^2=\var(Z)$

\item $Eh_F(Z,Z)= \frac12 E|Z_1-Z_2|=\int F(z)[1-F(z)]dz$

%\item As $k\rightarrow\infty$, $\lambda_k\rightarrow 0$.
\end{enumerate}

Now we can derive the desired result first using~(\ref{gamma1}),
then applying Fubini's theorem which is justified because the right
hand side of~(\ref{fubini1}) can be seen to be finite, and finally
using~(\ref{difint}): \newpage\begin{eqnarray*}
   \lefteqn{\int_{-\infty}^{\infty}[\gamma(z_1, w)-F(w)][\gamma(z_2,
   w)-F(w)]dw} \nonumber\\
   &=&
   E\int_{-\infty}^{\infty}
   [\gamma(z_1, w)-\gamma(Z_1, w)]
   [\gamma(z_2, w)-\gamma(Z_2, w)]
   dw \nonumber\\
   &=& -\half E\int_{-\infty}^{\infty}
   \left(
   [\gamma(z_1, w)-\gamma(z_2, w)]^2
   -
   [\gamma(z_1, w)-\gamma(Z_2, w)]^2-\right.
   \nonumber\\
   &&
   \left.
   [\gamma(Z_1, w)-\gamma(z_2, w)]^2
   -
   [\gamma(Z_1, w)-\gamma(Z_2, w)]^2
   \right)
   dw
   \nonumber\\
   &=&
   h_F(z_1,z_2)
\end{eqnarray*}

\begin{theorem}\label{th2a}
Suppose $h_F$ exists. Then there exists a complete system of
oscillating orthonormal functions $\{g_k\}$, each $g_k$ being
initially negative, continuous and satisfying
\begin{eqnarray}
   Eg_k(Z) = 0 \label{degen}
\end{eqnarray}
and nonnegative numbers $\lambda_1>\lambda_2>\ldots$ such that
\begin{eqnarray}
   h_F(Z_1,Z_2) = \sum_{k=0}^{\infty}\lambda_kg_k(Z_1)g_k(Z_2)
   \label{spectr1}
\end{eqnarray}
where convergence is in mean square.
\end{theorem}


\subsection{Estimation of the eigensystem from a sample}

The eigenvalues can also be estimated from a sample
$Z_1,\ldots,Z_n$. Some care must be taken to take tied observations
into account. To apply the above method, let $c_1<c_2<\ldots<c_t$ be
the distinct values of the observations $\{Z_i\}$ ($t\le n$) and let
$p_k$ be the proportion of the $Z_i$ equal to $c_k$.
\begin{lemma}
Let $\hat\lambda_k$ be the estimated $k$th eigenvalue based on an
iid sample of $N$ points. Then $\lambda_k$ is consistent iff
$\sum\lambda_k<\infty$.
\end{lemma}
\proof{Proof} ...
 \hspace*{\fill}$\Box$\endproof


\begin{figure}[h]
  \hfill
  \begin{minipage}[t]{.45\textwidth}
    \begin{center}
      \epsfig{file=c:/texfiles/newcorex1.eps, scale=0.5}
      \caption{Linear relationship}
      \label{fig-tc1}
    \end{center}
  \end{minipage}
  \hfill
  \begin{minipage}[t]{.45\textwidth}
    \begin{center}
      \epsfig{file=c:/texfiles/newcorex2.eps, scale=0.5}
      \caption{Curved relationship}
      \label{fig-tc2}
    \end{center}
  \end{minipage}
  \hfill \\
  \hfill
  \begin{minipage}[t]{.45\textwidth}
    \begin{center}
      \epsfig{file=c:/texfiles/newcorex3.eps, scale=0.5}
      \caption{Increasing dispersion}
      \label{fig-tc3}
    \end{center}
  \end{minipage}
  \hfill
  \begin{minipage}[t]{.45\textwidth}
    \begin{center}
      \epsfig{file=c:/texfiles/newcorex4.eps, scale=0.5}
      \caption{Increasing then decreasing dispersion}
      \label{fig-tc4}
    \end{center}
  \end{minipage}
  \hfill
\end{figure}



\begin{figure}
\begin{center}
\includegraphics[width=.6\linewidth,clip]{sim1.eps}\\~\\
\begin{tabular}{rr|rrrrr}
  & j  & i$\;\;\;\;\;\;\;\;$ 1 & 2 & 3 & 4 & 5\\ \hline
 & 1 & {\bf 0.68} & -{0.03} & -{0.08} &  -{0.03} & -{0.00} \\
& 2 & -{0.06} & {\bf 0.38} & {0.02} &    {0.08} & {0.06} \\
& 3 & {0.09} & -{0.03} & -{\bf 0.16} &    -{0.05} & -{0.04} \\
& 4 & {0.01} & {0.14} & {0.02} &    {0.07} & -{0.03} \\
& 5 & -{0.01} & -{0.03} & {0.08} &   -{0.05} & {0.07}
\end{tabular}
\end{center}
%\parbox{0.95\textwidth}
{\caption{Scatterplot of monotonely associated data and component
estimated correlations; $\rho^*=0.317$ and $\rho^2=0.???$
{\label{scat1}}}}
\end{figure}


\begin{figure}
\begin{center}
\includegraphics[width=.6\linewidth,clip]{sim2.eps}\\~\\
\begin{tabular}{rr|rrrrr}
  & j  & i$\;\;\;\;\;\;\;\;$ 1 & 2 & 3 & 4 & 5\\ \hline
 & 1 &  -{0.01} & -{0.00} & {0.01} &  {0.00} & {0.08} \\
 & 2 &  {\bf 0.81} & -{0.09} & -{0.02} &    {0.01} & {0.03} \\
 & 3 &  {0.03} & -{0.02} & {0.00} &    -{0.00} & -{0.04} \\
 & 4 &  {0.08} & {\bf 0.38} & -{0.06} &    -{0.06} & -{0.01} \\
 & 5 &  -{0.04} & {0.03} & {0.04} &    {0.02} & -{0.03}
\end{tabular}
\end{center}
%\parbox{0.95\textwidth}
{\caption{Scatterplot of quadratically associated data and component
correlations; $\rho^*=0.134$ and $\rho^2=0.0??$ {\label{scat2}}}}
\end{figure}


\begin{figure}
\begin{center}
\includegraphics[width=.6\linewidth,clip]{sim3.eps}
\end{center}
%\parbox{0.95\textwidth}
{\caption{Scatterplot of circularly associated data; $\rho^*=0.014$
and $\rho^2=0.0??$ {\label{scat3}}}}
\end{figure}



\begin{table}
\begin{center}
\begin{tabular}{rr|rrrrr}
  & j  & i$\;\;\;\;\;\;\;\;$ 1 & 2 & 3 & 4 & 5\\ \hline
 & 1 &  -{0.00} & {0.09} & -{0.01} &    {0.02} & -{0.03} \\
 & 2 &  -{0.11} & -{\bf 0.45} & -{0.01} &    -{0.07} & -{0.08} \\
 & 3 &  -{0.08} & {0.03} & {0.01} &    {0.02} & -{0.05} \\
 & 4 &  {0.04} & -{0.01} & {0.05} &    {0.03} & -{0.07} \\
 & 5 &  -{0.07} & {0.03} & {0.03} &    -{0.03} & -{0.07}
\end{tabular}
\caption{Component correlations of $\rho^*$
 \label{comp3}}
\end{center}
\end{table}


This paper introduces a correlation coefficient which has the
advantage compared to the ordinary correlation $\rho$ that it
detects arbitrary forms of association between two real random
variables $X$ and $Y$. For a distribution function $F$ let
\begin{eqnarray}
   h_F(z_1,z_2) =
   -
   \half
   E{\mbox{\Large $($}}\,
   |z_1-z_2| - |z_1-Z_2| - |Z_1-z_2| + |Z_1-Z_2|
   \,{\mbox{\Large )}}
   \label{hdef1}
\end{eqnarray}
where $Z_1$ and $Z_2$ are iid with distribution function $F$. Now
suppose $X$ and $Y$ have marginal distributions $F_1$ and $F_2$, let
$(X_1,Y_1),$ and $(X_2,Y_2)$ be iid replications of $(X,Y)$ and
consider
\begin{eqnarray*}
   \kappa(X,Y) = Eh_{F_1}(X_1,X_2)h_{F_2}(Y_1,Y_2)
\end{eqnarray*}
Now we can define
\begin{eqnarray*}
   \rho^*(X,Y)
   =
   \frac
   {\kappa(X,Y)}
   {\sqrt{\kappa(X,X)\kappa(Y,Y)}}
\end{eqnarray*}
In this paper we show that $0\le \rho^*(X,Y)\le 1$, such that
$\rho^*(X,Y)=0$ iff $X$ and $Y$ are independent, and $\rho^*(X,Y)=
1$ iff $X$ and $Y$ are linearly related.


In fact, $\kappa$ is quite closely related the the covariance and
hence $\rho^*$ is quite closely related to the ordinary correlation.
This can be seen as follows. Consider the generalized kernel
function
\begin{eqnarray}
   h^{(p)}_F(z_1,z_2) =
   -
   \half E
   {\mbox{\Large $($}}\,
   |z_1-z_2|^p - |z_1-Z_2|^p - |Z_1-z_2|^p + |Z_1-Z_2|^p
   \,{\mbox{\Large )}}
   \label{hdef2}
\end{eqnarray}
Note that $h_F=h_F^{(1)}$. Then using straightforward algebra we
obtain that
\begin{eqnarray}
   h^{(2)}_F(z_1,z_2) = (z_1-EZ)(z_2-EZ)
\end{eqnarray}
and
\begin{eqnarray}
   \cov(X,Y)^2 = Eh^{(2)}_{F_1}(X_1,X_2)h^{(2)}_{F_2}(Y_1,Y_2)
%   Eh^{(2)}_{F_1}(X_1,X_2)h^{(2)}_{F_2}(Y_1,Y_2) = E(X_1-EX)(X_2-EX)(Y_1-EY)(Y_2-EY) = [E(X-EX)(Y-EY)]^2
\end{eqnarray}
Thus, whereas $\kappa$ and $\rho^*$ are based on absolute
differences, the squared covariance and $\rho^2$ are based on
squared differences.


An intuitive interpretation of correlations based on $h^{(p)}_F$,
including $\rho^*$ and $\rho^2$ can be given. They measure how much
two $X$ observations which are ``far'' apart tend to go with $Y$
observations which are far apart, and similarly how much two $X$
observations which are ``close'' together tend to go with $Y$
observations which are close together. For all $p\ge 1$, this leads
to a maximal correlation of 1 if $X$ and $Y$ are linearly related.
For details, see Bergsma (2005).


\begin{figure}
\begin{center}
\includegraphics[width=.63\linewidth,clip]{c:/texfiles/fr2b.eps}
\end{center}
%\parbox{0.95\textwidth}
{\caption{Supports $S^+_{g_{k},g_l}$ (top) and $S^-_{g_{k},g_l}$
(bottom) of the Fr{\'e}chet bounds for the component correlations
$\rho_{ij}$ of $\rho^*(X,Y)$ when $X$ and $Y$ are uniformly
distributed on an interval.} {\label{univ fr}}}
\end{figure}


We now show that the $EZ$ does not exist. For $|V|\ge 1$, we have
\begin{eqnarray*}
   |Z| \ge \sqrt{\frac{\pi}2}\exp((|V|-1)^2/2)
\end{eqnarray*}
as is easy to verify. But then
\begin{eqnarray*}
   E|Z|
   &\ge&
   P(V\in[-1,1])\sqrt{\frac{\pi}2}\int_{[-1,1]} \exp(t^2/2)dt
   +
   P(V\not\in[-1,1])\sqrt{\frac{\pi}2}
   \int_{\bR\setminus[-1,1]}\exp(|v|-1)^2/2)\left(\frac{\exp(-v^2/2)}{\sqrt{2\pi}}\right)dv
   \\
   &\ge&
   P(V\in[-1,1])\sqrt{\frac{\pi}2}\int_{[-1,1]} \exp(t^2/2)dt
   +
   P(V\not\in[-1,1])\sqrt{\frac{\pi}2}
   \int_{\bR\setminus[-1,1]}\exp((1-2|v|)/2)\left(\frac{1}{\sqrt{2\pi}}\right)dv
   \\
   &=& \infty
\end{eqnarray*}



Since $\exp(t^2)$ is increasing in $t$, it is easy to verify that
for $0<\alpha<1$,
\begin{eqnarray*}
   |Z| \ge (1-\alpha)|V|e^{(\alpha^2V^2)}
\end{eqnarray*}
Taking $\alpha=1/\sqrt{2}$, we obtain
\begin{eqnarray*}
   E|Z|
   &\ge& (1-\half\sqrt{2})E|V|e^{V^2/2} = (1-\half\sqrt{2})\int_{-\infty}^{\infty}
   \left(|v|e^{v^2/2}\right)\left(\frac{e^{-v^2/2}}{\sqrt{2\pi}}\right)dv
   \\
   &=&\frac{1-\half\sqrt{2}}{\sqrt{2\pi}}\int_{-\infty}^{\infty}|v|dv
   =\infty
\end{eqnarray*}
In Section~\ref{spect} we show that $Eh_F(Z_1,Z_2)^2=6/\pi^2$.

On the real line, the inverse of $\erf$ is uniquely defined. For
imaginary $z$, we define the inverse $\erf^{-1}(z)$ as the unique
imaginary $y$ satisfying $z=\erf(y)$.

We see that (with $i=\sqrt{-1}$)
\begin{eqnarray*}
  Z = -\pi\,i\,\erf\left(\frac{i\,V}{\sqrt{2}}\right)
\end{eqnarray*}
Now $V$ has a standard normal distribution, so we obtain for the CDF
of $Z$:
\begin{eqnarray*}
  F(z) = P(Z<z) = P(V < -\sqrt{2}\,i\,\erf^{-1}(\pi\, i\,z))
  = \frac12
  \left(1-\erf\left(i\,\erf^{-1}\left(\frac{i\,z}{\pi}\right)\right)
  \right)
\end{eqnarray*}
From this,
\begin{eqnarray*}
  F^{-1}(u) = -\pi\,i\,\erf(i\erf^{-1}(2u-1))
\end{eqnarray*}
Some tedious but straightforward algebra then gives
\begin{eqnarray*}
  f(F^{-1}(u)) = 1/\phi[\Phi^{-1}(u)]^{2}
\end{eqnarray*}
Plugging this expression into~(\ref{de2}) leads to the solution
\begin{eqnarray*}
   \lambda_k &=& 1/k \\
   q_k(u) &\equiv& H_k[\Phi^{-1}(u)]\sqrt{\phi[\Phi^{-1}(u)]}
\end{eqnarray*}
where $H_k$ is the $k$th Hermite polynomial and the expression for
$q_k$ needs to be suitably normalized. For a derivation of
(essentially) this solution, see \citeA{dwv73}. They provided a
method for solving differential equations of the form
$\frac{d}{du}w(u)q_k'(u) + \lambda_k^{-1}q_k(u)=0$ for certain types
of weights $w(u)$, including $w(u)=1/\phi[\Phi^{-1}(u)]^{2}$.



Finally, in Section~\ref{normal}, calculations of the component
correlations for the normal distribution are given.

\subsection{Component correlations for the normal
distribution}\label{normal}

We show what the component correlations look like for the normal
distribution. Suppose $\{g_i\}$ and $\{\lambda_i\}$ are the
eigenfunctions and eigenvalues belonging to the standard normal
distribution and let
\begin{eqnarray*}
   m_{ij}^{(k)} &=& \int x^kg_i(z)g_j(z)dz
\end{eqnarray*}
Then:
\begin{theorem}
The component correlations $\rho_{ij}$ of $\rho^*(X,Y)$ where
$(X,Y)$ has a bivariate normal distribution with ordinary
correlation $\rho$ are the solutions of the infinite dimensional
linear system of equations of the form
\begin{eqnarray*}
   \rho_{ij} =
   \lambda_i\,\lambda_j\,\rho^2
   \sum_{k,l}\rho_{kl}\,
   c_{ijkl}
\end{eqnarray*}
where
\begin{eqnarray*}
   c_{ijkl}
   =
   m_{ik}^{(0)}m_{jl}^{(0)} + 2\rho m_{ik}^{(1)}m_{jl}^{(1)} +  + \rho^2m_{ik}^{(2)}m_{jl}^{(2)}
\end{eqnarray*}
\end{theorem}
\proof{Proof} qqq
 \hspace*{\fill}$\Box$\endproof

\begin{eqnarray*}
   \rho_{kl}(X,Y)=1 \Rightarrow
   \left\{
   \begin{array}{ll}
      \rho_{i\times k,i\times l}(X,Y)=1 & \forall i>1 \\
      \rho_{ij}(X,Y)=0   & \mbox{otherwise}\\
   \end{array}
   \right.
\end{eqnarray*}
and
\begin{eqnarray*}
   \rho_{kl}(X,Y)=-1 \Rightarrow
   \left\{
   \begin{array}{ll}
      \rho_{i\times k,i\times l}(X,Y)=(-1)^i & \forall i>1 \\
      \rho_{ij}(X,Y)=0   & \mbox{otherwise}\\
   \end{array}
   \right.
\end{eqnarray*}

abs_cor.bbl
\begin{thebibliography}{}

\bibitem[\protect\BCAY{Agarwal}{Agarwal}{1992}]{agarwal92}
Agarwal, R.~P.
\newblock{}\BBOP{}1992\BBCP{}.
\newblock{}\Bem{Difference equations and inequalities.}
\newblock{}New York: Marcel Dekker.

\bibitem[\protect\BCAY{Agresti}{Agresti}{1992}]{agresti92}
Agresti, A.
\newblock{}\BBOP{}1992\BBCP{}.
\newblock{}\BBOQ{}A survey of exact inference for contingency tables.\BBCQ{}
\newblock{}\Bem{Statistical Science}, \Bem{Vol. 7, No. 1}, 131-177.

\bibitem[\protect\BCAY{Agresti}{Agresti}{2002}]{agresti02}
Agresti, A.
\newblock{}\BBOP{}2002\BBCP{}.
\newblock{}\Bem{{C}ategorical {D}ata {A}nalysis, 2nd edition.}
\newblock{}New York: Wiley.

\bibitem[\protect\BCAY{Anderson \BBA{} Darling}{Anderson \BBA{}
  Darling}{1952}]{ad52}
Anderson, T.~W.\BCBT{} \BBA{} Darling, D.~A.
\newblock{}\BBOP{}1952\BBCP{}.
\newblock{}\BBOQ{}Asymptotic theory for certain `goodness of fit' criteria
  based on stochastic processes.\BBCQ{}
\newblock{}\Bem{Ann. Math. Stat.}, \Bem{47}, 193-212.

\bibitem[\protect\BCAY{Baglivo, Pagano\BCBL{} \BBA{} Spino}{Baglivo
  \BOthers{}}{1996}]{bps96}
Baglivo, J., Pagano, M.\BCBL{} \BBA{} Spino, C.
\newblock{}\BBOP{}1996\BBCP{}.
\newblock{}\BBOQ{}Permutation distributions via generating functions, with
  applications to sensitivity analysis of discrete data.\BBCQ{}
\newblock{}\Bem{J. Am. Stat. Ass.}, \Bem{91}, 1037-1046.

\bibitem[\protect\BCAY{Baringhaus \BBA{} Franz}{Baringhaus \BBA{}
  Franz}{2004}]{bf04}
Baringhaus, L.\BCBT{} \BBA{} Franz, C.
\newblock{}\BBOP{}2004\BBCP{}.
\newblock{}\BBOQ{}On a new multivariate two-sample test.\BBCQ{}
\newblock{}\Bem{Journal of multivariate analysis}, \Bem{88}, 190-206.

\bibitem[\protect\BCAY{Blum, Kiefer\BCBL{} \BBA{} Rosenblatt}{Blum
  \BOthers{}}{1961}]{bkr61}
Blum, J.~R., Kiefer, J.\BCBL{} \BBA{} Rosenblatt, M.
\newblock{}\BBOP{}1961\BBCP{}.
\newblock{}\BBOQ{}Distribution free tests of independence based on the sample
  distribution function.\BBCQ{}
\newblock{}\Bem{The annals of mathematical statistics}, \Bem{32}, 485-498.

\bibitem[\protect\BCAY{Booth \BBA{} Butler}{Booth \BBA{} Butler}{1999}]{bb99}
Booth, J.~G.\BCBT{} \BBA{} Butler, R.
\newblock{}\BBOP{}1999\BBCP{}.
\newblock{}\BBOQ{}An importance sampling algorithm for exact conditional tests
  in log-linear models.\BBCQ{}
\newblock{}\Bem{Biometrica}, \Bem{86}, 321-332.

\bibitem[\protect\BCAY{Bose \BBA{} Gupta}{Bose \BBA{} Gupta}{1959}]{bg59}
Bose, R.~C.\BCBT{} \BBA{} Gupta, S.~S.
\newblock{}\BBOP{}1959\BBCP{}.
\newblock{}\BBOQ{}Moments of order statistics from a normal population.\BBCQ{}
\newblock{}\Bem{Biometrika}, \Bem{46}, 433 - 440.

\bibitem[\protect\BCAY{{{D}e Wet}}{{{D}e Wet}}{1980}]{dewet80}
{{D}e Wet}, T.
\newblock{}\BBOP{}1980\BBCP{}.
\newblock{}\BBOQ{}Cram{\'e}r-von {M}ises tests for independence.\BBCQ{}
\newblock{}\Bem{J. Multivariate Anal.}, \Bem{10}, 38-50.

\bibitem[\protect\BCAY{{{D}e Wet}}{{{D}e Wet}}{1987}]{dewet87}
{{D}e Wet}, T.
\newblock{}\BBOP{}1987\BBCP{}.
\newblock{}\BBOQ{}Degenerate {U}- and {V}-statistics.\BBCQ{}
\newblock{}\Bem{South African Statistical Journal}, \Bem{21}, 99-129.

\bibitem[\protect\BCAY{{{D}e Wet} \BBA{} Venter}{{{D}e Wet} \BBA{}
  Venter}{1973}]{dwv73}
{{D}e Wet}, T.\BCBT{} \BBA{} Venter, J.~H.
\newblock{}\BBOP{}1973\BBCP{}.
\newblock{}\BBOQ{}Asymptotic distributions for quadratic forms with application
  to tests of fit.\BBCQ{}
\newblock{}\Bem{Annals of Statistics}, \Bem{1}, 380-387.

\bibitem[\protect\BCAY{Deheuvels}{Deheuvels}{1981}]{deheuvels81}
Deheuvels, P.
\newblock{}\BBOP{}1981\BBCP{}.
\newblock{}\BBOQ{}An asymptotic decomposition for multivariate
  distribution-free tests of independence.\BBCQ{}
\newblock{}\Bem{J. Multivariate Anal.}, \Bem{11}, 102-113.

\bibitem[\protect\BCAY{Diaconis \BBA{} Sturmfels}{Diaconis \BBA{}
  Sturmfels}{1998}]{ds98}
Diaconis, P.\BCBT{} \BBA{} Sturmfels, B.
\newblock{}\BBOP{}1998\BBCP{}.
\newblock{}\BBOQ{}Algebraic algorithms for sampling from conditional
  distributions.\BBCQ{}
\newblock{}\Bem{Ann. Stat.}, \Bem{13}, 363-397.

\bibitem[\protect\BCAY{Durbin \BBA{} Knott}{Durbin \BBA{} Knott}{1972}]{dk72}
Durbin, J.\BCBT{} \BBA{} Knott, M.
\newblock{}\BBOP{}1972\BBCP{}.
\newblock{}\BBOQ{}Components of the {C}ram{\'e}r-von {M}ises statistics,
  {I}.\BBCQ{}
\newblock{}\Bem{J. R. Statist. Soc. B}, \Bem{34}, 260-307.

\bibitem[\protect\BCAY{Eagleson}{Eagleson}{1979}]{eagleson79}
Eagleson, G.~K.
\newblock{}\BBOP{}1979\BBCP{}.
\newblock{}\BBOQ{}Orthogonal expansions and {U}-statistics.\BBCQ{}
\newblock{}\Bem{Australian Journal of Statistics}, \Bem{21}, 221-237.

\bibitem[\protect\BCAY{Forster, {McDonald}\BCBL{} \BBA{} Smith}{Forster
  \BOthers{}}{1996}]{fms96}
Forster, J.~J., {McDonald}, J.~W.\BCBL{} \BBA{} Smith, P. W.~F.
\newblock{}\BBOP{}1996\BBCP{}.
\newblock{}\BBOQ{}Monte {C}arlo exact conditional tests for log-linear and
  logistic models.\BBCQ{}
\newblock{}\Bem{J. Roy. Stat. Soc. Ser B}, \Bem{58}, 445-453.

\bibitem[\protect\BCAY{Goodman}{Goodman}{1985}]{goodman85}
Goodman, L.~A.
\newblock{}\BBOP{}1985\BBCP{}.
\newblock{}\BBOQ{}The analysis of cross-classified data having ordered and/or
  unordered categories: association models, correlation models, and asymmetry
  models for contingency tables with or without missing entries.\BBCQ{}
\newblock{}\Bem{Ann. Stat.}, \Bem{13}, 10-69.

\bibitem[\protect\BCAY{Govindarajulu}{Govindarajulu}{1963}]{govindarajulu63}
Govindarajulu, Z.
\newblock{}\BBOP{}1963\BBCP{}.
\newblock{}\BBOQ{}On moments of order statistics and quasi-ranges from normal
  populations.\BBCQ{}
\newblock{}\Bem{Ann. Math. Statist.}, \Bem{34}, 633-651.

\bibitem[\protect\BCAY{Gregory}{Gregory}{1977}]{gregory77}
Gregory, G.~G.
\newblock{}\BBOP{}1977\BBCP{}.
\newblock{}\BBOQ{}Large sample theory for {U}-statistics.\BBCQ{}
\newblock{}\Bem{Annals of Statistics}, \Bem{5}, 110-123.

\bibitem[\protect\BCAY{Hall}{Hall}{1979}]{hall79}
Hall, P.
\newblock{}\BBOP{}1979\BBCP{}.
\newblock{}\BBOQ{}On the invariance principle for {U}-statistics.\BBCQ{}
\newblock{}\Bem{Stoch. Proc. Appl.}, \Bem{9}, 163-174.

\bibitem[\protect\BCAY{Hoeffding}{Hoeffding}{1940}]{hoeffding40}
Hoeffding, W.
\newblock{}\BBOP{}1940\BBCP{}.
\newblock{}\BBOQ{}Masstabinvariante {K}orrelationtheorie.\BBCQ{}
\newblock{}\Bem{Schriften Math. Inst Univ. Berlin}, \Bem{5}, 181-233.

\bibitem[\protect\BCAY{Hoeffding}{Hoeffding}{1948\protect\BCnt{1}}]{hoeffding4%
8}
Hoeffding, W.
\newblock{}\BBOP{}1948\protect\BCnt{1}\BBCP{}.
\newblock{}\BBOQ{}A class of statistics with asymptotically normal
  distribution.\BBCQ{}
\newblock{}\Bem{Annals of Mathematical Statistics}, \Bem{19}, 293-325.

\bibitem[\protect\BCAY{Hoeffding}{Hoeffding}{1948\protect\BCnt{2}}]{hoeffding4%
8ind}
Hoeffding, W.
\newblock{}\BBOP{}1948\protect\BCnt{2}\BBCP{}.
\newblock{}\BBOQ{}A non-parametric test of independence.\BBCQ{}
\newblock{}\Bem{Annals of Mathematical Statistics}, \Bem{19}, 546-557.

\bibitem[\protect\BCAY{Hoeffding}{Hoeffding}{1961}]{hoeffding61}
Hoeffding, W.
\newblock{}\BBOP{}1961\BBCP{}.
\newblock{}\BBOQ{}The strong law of large numbers for {U}-statistics.\BBCQ{}
\newblock{}\Bem{Institute of Statistics, University of North Carolina,
  Mimeograph Series No. 302}.

\bibitem[\protect\BCAY{Johnson, Kotz\BCBL{} \BBA{} Balakrishnan}{Johnson
  \BOthers{}}{1994}]{jkb94}
Johnson, Kotz\BCBL{} \BBA{} Balakrishnan.
\newblock{}\BBOP{}1994\BBCP{}.
\newblock{}\Bem{Continuous univariate distributions: volume 1.}
\newblock{}New York: Wiley.

\bibitem[\protect\BCAY{Kallenberg \BBA{} Ledwina}{Kallenberg \BBA{}
  Ledwina}{1999}]{kl99}
Kallenberg, W. C.~M.\BCBT{} \BBA{} Ledwina, T.
\newblock{}\BBOP{}1999\BBCP{}.
\newblock{}\BBOQ{}Data driven rank tests for independence.\BBCQ{}
\newblock{}\Bem{J. Amer. Statist. Assoc.}, \Bem{94}, 285-301.

\bibitem[\protect\BCAY{Kiefer}{Kiefer}{1959}]{kiefer59}
Kiefer, J.
\newblock{}\BBOP{}1959\BBCP{}.
\newblock{}\BBOQ{}{$K$}-sample analogues of the {K}olmogorov-{S}mirnov and
  {C}ram{\'e}r-von {M}ises tests.\BBCQ{}
\newblock{}\Bem{Ann. Math. Stat.}, \Bem{30}, 420-447.

\bibitem[\protect\BCAY{Mikosch}{Mikosch}{2006}]{mikosch06}
Mikosch, T.
\newblock{}\BBOP{}2006\BBCP{}.
\newblock{}\BBOQ{}Copulas: tales and facts (with discussion).\BBCQ{}
\newblock{}\Bem{Extremes}, \Bem{To appear}.

\bibitem[\protect\BCAY{Nelsen}{Nelsen}{2006}]{nelsen06}
Nelsen, R.~B.
\newblock{}\BBOP{}2006\BBCP{}.
\newblock{}\Bem{An introduction to copulas.}
\newblock{}New York: Springer.

\bibitem[\protect\BCAY{Neyman}{Neyman}{1937}]{neyman37}
Neyman, J.
\newblock{}\BBOP{}1937\BBCP{}.
\newblock{}\BBOQ{}``{S}mooth'' tests for goodness of fit.\BBCQ{}
\newblock{}\Bem{Skand. Aktuarietidskr.}, \Bem{20}, 150-199.

\bibitem[\protect\BCAY{Randles \BBA{} Wolfe}{Randles \BBA{} Wolfe}{1979}]{rw79}
Randles, R.~H.\BCBT{} \BBA{} Wolfe, D.~A.
\newblock{}\BBOP{}1979\BBCP{}.
\newblock{}\Bem{Introduction to the theory of nonparametric statistics.}
\newblock{}New York: Wiley.

\bibitem[\protect\BCAY{R{\"u}schendorf}{R{\"u}schendorf}{1991}]{ruschendorf91}
R{\"u}schendorf, L.
\newblock{}\BBOP{}1991\BBCP{}.
\newblock{}\BBOQ{}Fr\'echet bounds and their applications.\BBCQ{}
\newblock{}In G.~Dall'Aglio, S.~Kotz\BCBL{} \BBA{} G.~Salinetti\ (\BEDS),
  \Bem{Advances in probabability distributions with given marginals}\ (\BPG\
  151-188).
\newblock{}Kluwer.

\bibitem[\protect\BCAY{Tricomi}{Tricomi}{1985}]{tricomi85}
Tricomi, F.~G.
\newblock{}\BBOP{}1985\BBCP{}.
\newblock{}\Bem{Integral equations.}
\newblock{}Dover Publications.

\bibitem[\protect\BCAY{{Van de Wiel}, {Di Bucchianico}\BCBL{} \BBA{} {Van der
  Laan}}{{Van de Wiel} \BOthers{}}{1999}]{vdwdbvdl99}
{Van de Wiel}, M.~A., {Di Bucchianico}, A.~D.\BCBL{} \BBA{} {Van der Laan}, P.
\newblock{}\BBOP{}1999\BBCP{}.
\newblock{}\BBOQ{}Symbolic computation and exact distributions of nonparametric
  test statistics.\BBCQ{}
\newblock{}\Bem{The statistician}, \Bem{48}, 507-516.

\bibitem[\protect\BCAY{Weisstein}{Weisstein}{1999}]{weisstein:erf-erfi}
Weisstein, E.~W.
\newblock{}\BBOP{}1999\BBCP{}.
\newblock{}\BBOQ{}{``Erf'' and ``Erfi''}.\BBCQ{}
\newblock{}\Bem{From MathWorld--A Wolfram Web Resource.
  http://mathworld.wolfram.com/Erfi.html}.

\bibitem[\protect\BCAY{Zaanen}{Zaanen}{1960}]{zaanen60}
Zaanen, A.
\newblock{}\BBOP{}1960\BBCP{}.
\newblock{}\Bem{Linear analysis.}
\newblock{}Amsterdam: North Holland Publishing Co.

\bibitem[\protect\BCAY{Zech \BBA{} Aslan}{Zech \BBA{} Aslan}{2003}]{za03}
Zech, G.\BCBT{} \BBA{} Aslan, B.
\newblock{}\BBOP{}2003\BBCP{}.
\newblock{}\BBOQ{}A multivariate two-sample test based on the concept of
  minimum energy.\BBCQ{}
\newblock{}\Bem{PHYSTAT2003}, 97-100.

\end{thebibliography}
