# MaxEnt-Copula-Oct2020

mesuvodeep

(September 2020)

###### Abstract

A new nonparametric model of maximum-entropy (MaxEnt) copula density function
is proposed, which offers the following advantages: (i) it is valid for mixed
random vector. By ‘mixed’ we mean the method works for any combination of
discrete or continuous variables in a fully automated manner; (ii) it yields a
bonafide density estimate with intepretable parameters. By ‘bonafide’ we mean
the estimate guarantees to be a non-negative function, integrates to $1$; and
(iii) it plays a unifying role in our understanding of a large class of
statistical methods.

A Maximum Entropy Copula Model for Mixed Data:
Representation, Estimation, and Applications
Deep Mukhopadhyay
<EMAIL_ADDRESS>

Keywords: Maximum entropy copula learning; Mixed data smoothing; Mutual
information; Categorical data analysis; Log-linear modeling; Copula-logistic
regression.

###### Contents

  1. 1 The Problem
  2. 2 Auto-Adaptive Nonparametric Models
    1. 2.1 Background Concepts and Notation
    2. 2.2 Log-bilinear Model
    3. 2.3 Log-linear Model
    4. 2.4 A Few Examples
  3. 3 Applications to Statistical Modeling
    1. 3.1 Goodman’s Association Model For Categorical Data
    2. 3.2 Logratio biplot: Graphical Exploratory Analysis
    3. 3.3 Loglinear Modeling of Large Sparse Contingency Tables
    4. 3.4 Modeling Bivariate Discrete Distributions
    5. 3.5 Mutual Information
    6. 3.6 Nonparametric Copula-Logistic Regression
  4. 4 Conclusion
  5. 5 Supplementary Appendix
    1. A.1 Nonparametric LP-Polynomials
    2. A.2 Two Cultures of Maximum Entropy Modeling
    3. A.3 Connection With Vladimir Vapnik’s Statistical Invariants
  6. A.4 Additional Figures and Tables

##  1 The Problem

A new class of nonparametric maximum-entropy (MaxEnt) approach of copula
modeling is discussed. The proposed statistical formalism has the following
advantages: First, it yields a bonafide (smooth, non-negative, and integrates
to $1$) copula density estimate with interpretable parameters that provide
insights into the nature of the association between the random variables
$(X,Y)$. Secondly, the method is data-type agnostic—which is to say that it
automatically adapts to mixed-data types (any combination of discrete,
continuous, or even categorical). Thirdly, and most notably, our copula
representation theory subsumes and unifies a wide range of statistical
learning methods using a common mathematical notation—unlocking deep,
surprising connections and insights, which were previously unknown.

##  2 Auto-Adaptive Nonparametric Models

In this section, we introduce two new maximum-entropy (MaxEnt) copula density
models. But before diving into technical details, it will be instructive to
review some basic definitions and concepts related to copula.

###  2.1 Background Concepts and Notation

Sklar’s Copula Representation Theory (Sklar, 1959). The joint cumulative
distribution function (cdf) of any pair of random variables $(X,Y)$

$F_{X,Y}(x,y)=\Pr(X\leq x,Y\leq y),\leavevmode\nobreak\ \leavevmode\nobreak\
\mathrm{for}\leavevmode\nobreak\ (x,y)\in\mathcal{R}^{2}$

can be decomposed as a function of the marginal cdfs $F_{X}$ and $F_{Y}$

$F_{X,Y}(x,y)\,=\,\operatorname{Cop}_{X,Y}\big{(}F_{X}(x),F_{Y}(y)\big{)},\leavevmode\nobreak\
\leavevmode\nobreak\ \mathrm{for}\leavevmode\nobreak\ (x,y)\in\mathcal{R}^{2}$
(2.1)

where $\operatorname{Cop}_{X,Y}$ denotes a copula distribution function with
uniform marginals. To set the stage, we start with the continuous marginals
case, which will be generalized later to allow mixed-$(X,Y)$. Taking
derivative of Eq. (2.1), we get

$f_{X,Y}(x,y)\,=\,f_{X}(x)f_{Y}(y)\operatorname{cop}_{X,Y}\big{(}F_{X}(x),F_{Y}(y)\big{)},\leavevmode\nobreak\
\leavevmode\nobreak\ \mathrm{for}\leavevmode\nobreak\ (x,y)\in\mathcal{R}^{2}$
(2.2)

which decouples the joint density into the marginals and the copula. One can
rewrite Eq. (2.2) to represent copula as a “normalized” joint density function

$\operatorname{cop}_{X,Y}\big{(}F_{X}(x),F_{Y}(y)\big{)}\,:=\,{\rm
dep}_{X,Y}(x,y)\,=\dfrac{f_{X,Y}(x,y)}{f_{X}(x)f_{Y}(y)},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ $ (2.3)

which is also known as the dependence function, pioneered by Hoeffding (1940).
To make (2.3) a proper density function (i.e., one that integrates to one) we
perform quantile transformation by substituting
$F_{X}(x)=u\leavevmode\nobreak\ \text{and}\leavevmode\nobreak\ F_{Y}(y)=v$:

$\operatorname{cop}_{X,Y}(u,v)={\rm
dep}_{X,Y}(F_{X}^{-1}(u),F_{Y}^{-1}(v))=\dfrac{f_{X,Y}\big{(}F_{X}^{-1}(u),F_{Y}^{-1}(v)\big{)}}{f_{X}(F_{X}^{-1}(u))f_{Y}(F_{Y}^{-1}(v))},\leavevmode\nobreak\
\leavevmode\nobreak\ (u,v)\in[0,1]^{2}.$ (2.4)

We are now ready to extend this copula density concept to the mixed $(X,Y)$
case. Pre-Copula: Conditional Comparison Density (Parzen and Mukhopadhyay,
2013). Before we introduce the generalized copula density, we need to
introduce a new concept—conditional comparison density (CCD). For a continuous
random variable $X$, define

$d(u;X,X|Y=y)\,=\,\dfrac{f_{X|Y}\big{(}F_{X}^{-1}(u)\mid
y\big{)}}{f_{X}\big{(}F_{X}^{-1}(u)\big{)}},\leavevmode\nobreak\
\leavevmode\nobreak\ 0<u<1$ (2.5)

For $Y$ discrete, we represent it using probability mass function (pmf):

$d(v;Y,Y|X=x)\,=\,\dfrac{p_{Y|X}\big{(}Q_{Y}(v)|x\big{)}}{p_{Y}\big{(}Q_{Y}(v)\big{)}}\,=\,\dfrac{\Pr(Y=Q_{Y}(v)|X=x)}{\Pr(Y=Q_{Y}(v))},\leavevmode\nobreak\
\leavevmode\nobreak\ 0<v<1.$ (2.6)

where $Q_{Y}(v)$ is the quantile function of $Y$. It is easy to see that the
CCDs (2.5) and (2.6) are proper densities in the sense that

$\int_{0}^{1}d(u;X,X|Y=y)\;\mathrm{d}u\,=\,\int_{0}^{1}d(v;Y,Y|X=x)\;\mathrm{d}v\,=\,1.$

Generalized Copula Representation Theory (Mukhopadhyay and Parzen, 2020). For
the mixed case, when $Y$ is discrete and $X$ is continuous the joint density
of (2.3) is defined by either side of the following identity:

$\Pr(Y\mid X=x)f_{X}(x)\,=\,f_{X|Y}(x|y)\Pr(Y=y).$ (2.7)

This can be rewritten as the ratios of conditionals and their respective
marginals:

$\text{Pre-Bayes' Rule}:\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \dfrac{\Pr(Y=y\mid
X=x)}{\Pr(Y=y)}\,=\,\dfrac{f_{X|Y}(x|y)}{f_{X}(x)}\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ $ (2.8)

This formula (2.8) can be interpreted as the slices of the mixed copula
density, since

$\operatorname{cop}_{X,Y}\big{(}F_{X}(x),F_{Y}(y)\big{)}\,=\,\dfrac{f_{X|Y}(x|y)}{f_{X}(x)}\,=\,\dfrac{\Pr(Y\mid
X=x)}{\Pr(Y=y)}.\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ $ (2.9)

Substituting $F_{X}(x)=u$ and $F_{Y}(y)=v$, we get the following definition of
the generalized copula density in terms of conditional comparison density
(CCD):

$\operatorname{cop}_{X,Y}(u,v)\,=\,d\big{(}u;X,X|Y=Q(v;Y)\big{)}\,=\,d\big{(}v;Y,Y|X=Q(u;X)\big{)},\leavevmode\nobreak\
\leavevmode\nobreak\ 0<u,v<1.$ (2.10)

Bayes’ theorem ensures the equality of two CCDs, with copula being the common
value. Based on this foundation, we are now ready to develop the MaxEnt copula
modeling theory.

###  2.2 Log-bilinear Model

An exponential Fourier series representation of copula density function is
given. The reasons for entertaining an exponential model for copula is
motivated from two different perspectives.

The Problem of Unboundedness. One peculiar aspect of copula density function
is that it can be unbounded at the corners of the unit square. In fact, many
common parametric copula families—Gaussian, Clayton, Gumbel, etc.—tend to
infinity at the boundaries. So naturally the question arises: How to develop
suitable approximation methods that can accommodate a broader class of copula
density shapes, including the unbounded ones? The first key insight: logarithm
of the copula density function is far more convenient to approximate (due to
its well-behaved nature) than the original density itself. We thus express the
logarithm of copula density $\log\operatorname{cop}_{X,Y}$ in the Fourier
series—instead of doing canonical $L_{2}$ approximation, which expands
$\operatorname{cop}_{X,Y}$ directly in an orthogonal series (Mukhopadhyay and
Parzen, 2020). Accordingly, for densities with rapidly changing tails, ‘log-
Fourier’ method leads to an improved estimate that is less wiggly and more
parsimonious than the $L_{2}$-orthogonal series model. In addition, the
resulting exponential form guarantees the non-negativity of the estimated
density function.

Choice of Orthonormal Basis. We choose nonparametric LP-polynomial basis (see
Appendix A.1) to expand log-copula density function. LP-bases’ appeal lies in
its ability to approximate the quirky shapes of mixed-copula functions in a
completely automated way; see Fig. 1. Consequently, it provides a unified way
to develop nonparametric smoothing algorithms that simultaneously hold for
mixed data types.

######  Definition 1.

The exponential copula model admits the following LP-expansion

$\operatorname{cop}_{\bm{\theta}}(u,v;X,Y)\leavevmode\nobreak\
=\leavevmode\nobreak\
\dfrac{1}{Z_{\bm{\theta}}}\exp\Big{\\{}\mathop{\sum\sum}_{j,k>0}\theta_{jk}S_{j}(u;X)S_{k}(v;Y)\Big{\\}},$
(2.11)

where $Z_{\bm{\theta}}$ is the normalization factor that ensures
$\operatorname{cop}_{\bm{\theta}}$ is a proper density

$Z_{\bm{\theta}}\,=\,\iint_{[0,1]^{2}}\exp\Big{\\{}\mathop{\sum\sum}_{j,k>0}\theta_{jk}S_{j}(u;X)S_{k}(v;Y)\Big{\\}}\;\mathrm{d}u\;\mathrm{d}v.$
(2.12)

We refer (2.11) as the log-bilinear copula model.

The Maximum Entropy Principle. Another justification for choosing the
exponential model comes from the principle of maximum entropy (MaxEnt),
pioneered by E. T. Jaynes (1957). The maxent principle defines a unique
probability distribution by maximizing the entropy
$H(\operatorname{cop})=-\int\operatorname{cop}_{X,Y}\log\operatorname{cop}_{X,Y}$
under the normalization constraint $\int\operatorname{cop}_{X,Y}=1$ and the
following LP-co-moment conditions:

$\mathbb{E}_{\operatorname{cop}_{{\bm{\theta}}}}[S_{j}(U;X)S_{k}(V;Y)]\,=\,\operatorname{LP}_{jk}.\leavevmode\nobreak\
$ (2.13)

LP-co-means are orthogonal “moments” of copula, which can be estimated by

$\widetilde{\operatorname{LP}}_{jk}\,=\,\mathbb{E}_{\widetilde{\operatorname{Cop}}}\big{[}S_{j}(U;X)S_{k}(V;Y)\big{]}\,=\,\dfrac{1}{n}\sum_{i=1}^{n}S_{j}\big{(}\widetilde{F}_{X}(x_{i});X\big{)}\,S_{k}\big{(}\widetilde{F}_{Y}(y_{i});Y\big{)}.$
(2.14)

Applying calculus of variations, one can show that the maxent constrained
optimization problem leads to the exponential (2.11) form. The usefulness of
Jaynes’ maximum entropy principle lies in providing a constructive mechanism
to uniquely identify a probability distribution that is maximally non-
committal (flattest possible) with regard to all unspecified information
beyond the given constraints.

Estimation. We fit a truncated exponential series estimator of copula density

$\operatorname{cop}_{\bm{\theta}}(u,v;X,Y)\leavevmode\nobreak\
=\leavevmode\nobreak\
\dfrac{1}{Z_{\bm{\theta}}}\exp\Big{\\{}\sum_{j=1}^{m_{1}}\sum_{k=1}^{m_{2}}\theta_{jk}S_{j}(u;X)S_{k}(v;Y)\Big{\\}}$

The task of finding the maximum likelihood estimates (MLE) of ${\bm{\theta}}$
boils down to solving the following sets of equations for $j=1,\ldots,m_{1}$
and $k=1,\ldots,m_{2}$:

$\frac{\partial\log\operatorname{cop}_{\bm{\theta}}}{\partial\theta_{jk}}\leavevmode\nobreak\
\equiv\leavevmode\nobreak\ \frac{\partial\log
Z_{\bm{\theta}}}{\partial\theta_{jk}}\,-\,\dfrac{1}{n}\sum_{i=1}^{n}S_{j}\big{(}\widetilde{F}_{X}(x_{i});X\big{)}\,S_{k}\big{(}\widetilde{F}_{Y}(y_{i});Y\big{)}\,=\,0.\leavevmode\nobreak\
\leavevmode\nobreak\ $ (2.15)

Note that the derivative of the log-partition function is equal to the
expectation of the LP-co-mean functions:

$\dfrac{\partial\log
Z_{\bm{\theta}}}{\partial\theta_{jk}}\,=\,\mathbb{E}_{\operatorname{cop}_{{\bm{\theta}}}}[S_{j}(U;X)S_{k}(V;Y)].$
(2.16)

Replacing (2.16) and (2.14) into (2.15) implies that the MLE of MaxEnt model
is same as the method of moments estimator satisfying the following moment
conditions:

$\iint_{[0,1]^{2}}S_{j}(u;X)S_{k}(v;Y)\operatorname{cop}_{\bm{\theta}}(u,v;X,Y)\;\mathrm{d}u\;\mathrm{d}v\,=\,\widetilde{\operatorname{LP}}_{jk}.$
(2.17)

At this point, one can apply any convex optimization111since the second
derivative of the log-partition function is a s positive semi-definite
covariance matrix routine (e.g., Newton’s method, gradient descent, stochastic
gradient descent, etc.) to solve for $\widehat{{\bm{\theta}}}$.

Asymptotic. Let the sequence of $m_{1}$ and $m_{2}$ increase with sample size
with an appropriate rate $\frac{(m_{1}m_{2})^{3}}{n}\to 0$ as
$n\rightarrow\infty$. Then, under certain suitable regularity conditions, the
exponential $\operatorname{cop}_{\widehat{\bm{\theta}}}$ is a consistent
estimate in the sense of Kullback-Leibler distance; see Barron and Sheu (1991)
for more details.

Determining Informative Constraints. Jayne’s maximum entropy principle assumes
that a proper set of constraints (i.e., sufficient statistic functions) are
given to the modeler, one that captures the phenomena under study. This
assumption may be legitimate for studying thermodynamic experiments in
statistical mechanics or for specifying prior distribution in Bayesian
analysis, but certainly not for building empirical models.

Which comes first: a parametric model or sufficient statistics? After all, the
identification of significant components (sufficient statistics) is a
prerequisite for constructing a legitimate probability model from the data
(Mukhopadhyay et al., 2012). Therefore the question of how to judiciously
design and select the constraints from data, seems inescapable for
nonparametrically learning maxent copula density function from data; also see
Appendix A.3, which discusses the ‘two cultures’ of maxent modeling. We
address this issue as follows: (i) compute
$\widetilde{\operatorname{LP}}_{jk}$ using the formula eq. (2.14); (ii) sort
them in descending order based on their magnitude (absolute value); (iii)
compute the penalized ordered sum of squares

${\rm PenSum}(q)\leavevmode\nobreak\ =\leavevmode\nobreak\ \text{Sum of
squares of top $q$ LP comeans}\leavevmode\nobreak\ -\leavevmode\nobreak\
\dfrac{\gamma_{n}}{n}q.\leavevmode\nobreak\ \leavevmode\nobreak\ $

For AIC penalty choose $\gamma_{n}=2$, for BIC choose $\gamma_{n}=\log n$,
etc. Further details can be found in Mukhopadhyay and Parzen (2020, Sec. 4.3).
(iv) Find the $q$ that maximizes the ${\rm PenSum}(q)$. Store the selected
indices $(j,k)$ in the set $\mathcal{I}$. (v) Carry out maxent optimization
routine based only on the selected LP-sufficient statistics-based constraints:

$\Big{\\{}S_{j}(u;X)S_{k}(v;Y)\Big{\\}},\leavevmode\nobreak\
\leavevmode\nobreak\ (j,k)\in\mathcal{I}.\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ $

This pruning strategy guards against overfitting. Finally, return the
estimated reduced-order (with effective dimension $|\mathcal{I}|$) maxent
copula model.

######  Remark 1 (Nonparametric MaxEnt).

The proposed nonparametric maxent mechanism produces a copula density
estimate, which is flexible (can adapt to the ‘shape of the data’ without
making risky a priori assumptions) and yet possesses a compact analytical
form.

###  2.3 Log-linear Model

We provide a second parameterization of copula density.

######  Definition 2.

The log-linear orthogonal expansion of LP-copula is given by:

$\operatorname{cop}_{\bm{\mu}}(u,v;X,Y)\leavevmode\nobreak\
=\leavevmode\nobreak\
\exp\big{\\{}\mu_{0}+\sum_{k>0}\mu_{k}\,\phi_{k}(u;X)\,\psi_{k}(v;Y)\big{\\}},\leavevmode\nobreak\
\leavevmode\nobreak\ 0<u,v<1,$ (2.18)

We call the parameters of this model “log-linear LP-correlations” that satisfy
for $k>0$

$\mu_{k}=\iint_{[0,1]^{2}}\phi_{k}(u;X)\psi_{k}(v;Y)\log\operatorname{cop}_{\bm{\mu}}(u,v;X,Y)\;\mathrm{d}u\;\mathrm{d}v,\leavevmode\nobreak\
\leavevmode\nobreak\ $ (2.19)

Connection. Two fundamental representations, namely the log-bilinear (2.18)
and loglinear (2.11) copula models, share some interesting connections. To see
that perform singular value decomposition (SVD) of the $\Theta$-matrix whose
$(j,k)$th entry is $\theta_{jk}$:

$\Theta\,=\,U\Omega V^{T},$

$u_{ij}$ and $v_{ij}$ are the elements of the singular vectors with singular
values $\mu_{1}\geq\mu_{2}\geq\cdots\geq 0$. Then the spectral bases can be
expressed as the linear combinations of the LP-polynomials:

$\displaystyle\phi_{k}(u;X)\,=\,\sum\nolimits_{j}u_{jk}S_{j}(u;X)\leavevmode\nobreak\
\leavevmode\nobreak\ $ (2.20)
$\displaystyle\psi_{k}(u;Y)\,=\,\sum\nolimits_{l}v_{lk}S_{l}(v;Y).\leavevmode\nobreak\
\leavevmode\nobreak\ $ (2.21)

Hence, the LP-spectral functions (2.20-2.21) satisfy the following
orthonormality conditions:

$\displaystyle\int\phi_{k}(u;X)\;\mathrm{d}u\,=\,\int\psi_{k}(v;Y)\;\mathrm{d}v\,=\,0\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ $ $\displaystyle\hskip
10.00002pt\int\phi_{j}(u;X)\psi_{k}(u;Y)\;\mathrm{d}u\,=\,\delta_{jk},\,\text{for}\leavevmode\nobreak\
j\neq k.\leavevmode\nobreak\ \leavevmode\nobreak\ \vspace{-.6em}$

###  2.4 A Few Examples

We demonstrate the flexibility of the LP-copula models using real data
examples.

######  Example 1.

Kidney fitness data (Efron and Hastie, 2016, Sec 1.1). It contains
measurements on $n=157$ healthy volunteers (potential donors). For each
volunteer, we have their age (in years) and a composite measure “tot” of
overall function of kidney function. To understand the relationship between
age and tot, we estimate the copula:

$\small{\widehat{\operatorname{cop}}_{X,Y}(u,v)\,=\,\exp\big{\\{}-.40S_{1}(u;X)S_{1}(v;Y)+.18S_{2}(u;X)S_{2}(v;Y)-0.12\big{\\}},}$
(2.22)

displayed in Fig. 1(a). At the global scale, the shape of the copula density
indicates a prominent negative ($\widehat{\theta}_{11}=-0.40$) association
between age and tot. Moreover, at the local scale, significant heterogeneity
of the strength of dependence is clearly visible, as captured by the nonlinear
asymmetric copula: the correlation between age and tot is quite high for older
(say, $>70$) donors, compared to younger ones. This allows us to gain refined
insights into how kidney function declines with age.

Figure 1: MaxEnt LP-copula density estimates for (a) Kidney fitness (age vs
tot: both continuous marginals), (b) PLOS data (title length vs number of
authors: both discrete marginals), (c) horseshoe crab data (number of
satellites vs width: mixed discrete-continuous marginals), and (d) Challenger
space shuttle data (temperature vs number of O-ring failures: mixed
continuous-discrete marginals).

######  Example 2.

PLOS data (both discrete marginals). It contains information on $n=878$
journal articles published in PLOS Medicine between 2011 and 2015. For each
article, two variables were extracted: length of the title and the number of
authors. The dataset is available in the R-package dobson. The checkerboard-
shaped estimated discrete copula

$\displaystyle\widehat{\operatorname{cop}}_{X,Y}(u,v)\,=\,\exp\big{\\{}.42S_{1}(u;X)S_{1}(v;Y)+.10S_{2}(u;X)S_{2}(v;Y)-.07S_{2}(u;X)S_{2}(v;Y)-0.12\big{\\}}.$

is shown in Fig. 1(b), which shows a strong positive nonlinear association. In
particular, the sharp lower-tail, around the $(0,0)$, indicates that the
smaller values of ($X,Y$) have a greater tendency to occur together than the
larger ones.

######  Example 3.

Horseshoe Crabs Data (mixed marginals). The study consists of $n=173$ nesting
horseshoe crabs (Agresti, 2013). For each female crab in the study, we have
its carapace width (cm) and number of male crabs residing nearby her nest. The
goal of the study is to investigate whether carapace width affects number of
male satellites for the female horseshoe crabs. If so, how–what is the shape
of the copula dependence function? The estimated copula, shown in Fig. 1(c),
is given by:

$\small{\widehat{\operatorname{cop}}_{X,Y}(u,v)\,=\,\exp\Big{\\{}0.375S_{1}(u;X)S_{1}(v;Y)-0.077\Big{\\}}.}$

This indicates a significant positive linear correlation between the width and
number of satellites of a female crab.

######  Example 4.

1986 Challenger Shuttle O-Ring data. On January 28, 1986, just after seventy-
three seconds into the flight, Challenger space shuttle broke apart, killing
all seven crew members on board. The purpose of this study is to investigate
whether the ambient temperature during the launch was related to the damage of
shuttle’s O-rings. For that we have $n=23$ previous shuttle missions data,
consisting of launch temperatures (degrees F), and number of damaged O-rings
(out of 6). The estimated LP-maxent copula density

$\small{\widehat{\operatorname{cop}}_{X,Y}(u,v)\,=\,\exp\Big{\\{}-0.37S_{1}(u;X)S_{1}(v;Y)-0.27S_{3}(u;X)S_{1}(v;Y)-0.14\Big{\\}}}$

is displayed in panel (d), and shows a strong negative association between the
temperature at the launch and the number of damaged O-rings. Moreover, the
sharp peak of the copula density around the edge $(0,1)$ further implies that
cold temperatures can excessively increase the risk of failure of the o-rings.

##  3 Applications to Statistical Modeling

The scope of the general theory of the preceding section goes far beyond
simply a tool for nonparametric copula approximation. In this section, we show
how one can derive a large class of applied statistical methods in a unified
manner by suitably reformulating them in terms of the LP-maxent copula model.
In doing so, we also provide statistical interpretations of LP-maxent
parameters under different data modeling tasks.

###  3.1 Goodman’s Association Model For Categorical Data

Categorical data analysis will be viewed through the lens of LP-copula
modeling. Let $X$ and $Y$ denote two discrete categorical variables; $X$ with
$I$ categories and $Y$ with $J$ categories. The data are summarized in an
$I\times J$ contingency table $\mathbf{F}$: $f_{kl}$ is the observed cell
count in row $k$ and column $l$ of $\mathbf{F}$ and $n=f_{++}$ is the total
frequency. The row and column totals are denoted as $f_{k+}$ and $f_{+l}$. The
observed joint $\Pr(X=k,Y=l)$ is denoted by
${\widetilde{p}}_{kl}=f_{kl}/f_{++}$; the respective row and column marginals
are given by ${\widetilde{p}}_{k+}=f_{k+}/f_{++}$ and
${\widetilde{p}}_{+l}=f_{+l}/f_{++}$.

LP Log-linear Model. We specialize our general copula model (2.18) for two-way
contingency tables. The discrete LP-copula for the $I\times J$ table is given
by

$\operatorname{cop}(F_{X}(k),F_{Y}(l))\,=\,\exp\Big{(}\mu_{0}+\sum_{j=1}^{m}\mu_{j}\phi_{jk}\psi_{jl}\Big{)},$
(3.1)

where we abbreviate the row and columns scores $\phi_{j}(F_{X}(k))=\phi_{jk}$
and $\psi_{j}(F_{Y}(l))=\psi_{jl}$ for $k=1,2,\ldots,I$ and $l=1,2,\ldots,J$.
The number of components $m\leq M=\min(I-1,J-1)$; we call the log-linear model
(3.1) ‘saturated’ (or ‘dense’) when we have $m=M$ components. The non-
increasing sequence of model parameters $\mu_{j}$’s are called “intrinsic
association parameters” that satisfy222Compare our equation (3.2) with
equation (34) of Goodman (1996).

$\mu_{j}\,=\,\sum_{k=1}^{I}\sum_{l=1}^{J}\Big{(}\log
p_{kl}\Big{)}p_{k+}p_{+l}\phi_{jk}\psi_{jl},\leavevmode\nobreak\
\leavevmode\nobreak\ \text{for}\leavevmode\nobreak\ j=1,\ldots,m.$ (3.2)

Note that the discrete LP-row and column scores, by design, satisfy (for
$j\neq j^{\prime}$):

$\displaystyle\sum_{k=1}^{I}\phi_{jk}p_{k+}\,=\leavevmode\nobreak\
\sum_{l=1}^{J}\psi_{jl}p_{+l}\leavevmode\nobreak\ =\leavevmode\nobreak\ 0$
(3.3) $\displaystyle\sum_{k=1}^{I}\phi^{2}_{jk}p_{k+}\,=\leavevmode\nobreak\
\sum_{l=1}^{J}\psi^{2}_{jl}p_{+l}\leavevmode\nobreak\ =\leavevmode\nobreak\ 1$
(3.4)
$\displaystyle\sum_{k=1}^{I}\phi_{jk}\phi_{j^{\prime}k}p_{k+}\,=\leavevmode\nobreak\
\sum_{l=1}^{J}\phi_{jl}\psi_{j^{\prime}l}p_{+l}\,=\,0.$ (3.5)

Interpretation. It is clear from (3.2) that the parameters $\mu_{j}$’s are
fundamentally different from the standard Pearsonian-type correlation

$\rho_{j}\,=\,\sum_{k=1}^{I}\sum_{l=1}^{J}p_{kl}\phi_{jk}\psi_{jl},\leavevmode\nobreak\
\leavevmode\nobreak\ \text{for}\,j=1,\ldots,m.$ (3.6)

The coefficients of the LP-MaxEnt-copula expansion for contingency tables
carry a special interpretation in terms of log-odds-ratio. To see this we
start by examining the $2\times 2$ case. The $2\times 2$ Contingency Table.
Applying (3.2) for two-by-two tables we have

$\mu_{1}\,=\,\sum_{k=0}^{1}\sum_{l=0}^{1}\Big{(}\log
p_{kl}\Big{)}p_{k+}p_{+l}\phi_{1k}\psi_{1l}.$ (3.7)

Note that for dichotomous $X$ the LP-spectral basis $\phi_{1}(F_{X}(x))$ is
equal to $T_{1}(x;F_{X})$. Consequently, we have the following explicit
formula for $\phi_{1}$ and $\psi_{1}$:

$\displaystyle\phi_{1}(F_{X}(x))$ $\displaystyle=$
$\displaystyle\dfrac{x-p_{2+}}{\sqrt{p_{1+}p_{2+}}}$ (3.8)
$\displaystyle\psi_{1}(F_{Y}(y))$ $\displaystyle=$
$\displaystyle\dfrac{y-p_{+2}}{\sqrt{p_{+1}p_{+2}}}$ (3.9)

Substituting this into (3.7) yields the following important result.

######  Theorem 1.

For 2-by-2 contingency tables, the estimate of the statistical parameter
$\mu_{1}$ of the maxent LP-copula model

$\operatorname{cop}_{\bm{\mu}}(u,v;X,Y)\,=\,e^{\mu_{0}+\mu_{1}\phi_{1}(u;X)\psi_{1}(v;Y)}$

can be expressed as follows:

$\widehat{\mu}_{1}\,=\,\log\Bigg{[}\dfrac{{\widetilde{p}}_{11}{\widetilde{p}}_{22}}{{\widetilde{p}}_{12}{\widetilde{p}}_{21}}\Bigg{]}\Big{(}{\widetilde{p}}_{1+}{\widetilde{p}}_{+1}{\widetilde{p}}_{2+}{\widetilde{p}}_{+2}\Big{)}^{1/2},$
(3.10)

where the part inside the square bracket is the sample log-odds-ratio.

######  Remark 2 (Significance of Theorem 1).

We have derived the log-odds-ratio statistic from first principles using a
copula-theoretic framework. To the best of our knowledge, no other study has
discovered this connection; see also eq.(16) of Goodman (1991). In fact, one
can view Theorem 1 as a special case of the much more general result described
next.

######  Theorem 2.

For an $I\times J$ table, consider a two-by-two subtable with rows $k$ and
$k^{\prime}$ and columns $l$ and $l^{\prime}$. Then the logarithm of odds-
ratio $\vartheta_{kl,k^{\prime}l^{\prime}}$ is connected with the intrinsic
association parameters $\mu_{j}$ in the following way:

$\log\vartheta_{kl,k^{\prime}l^{\prime}}\,=\,\sum_{j=1}^{M}\mu_{j}\big{(}\phi_{jk}-\phi_{jk^{\prime}}\big{)}\big{(}\psi_{jl}-\psi_{jl^{\prime}}\big{)}.$
(3.11)

To deduce (3.10) from (3.11), verify the following, utilizing the LP-basis
formulae (3.8)-(3.9)

$\displaystyle\phi_{11}-\phi_{10}:=\phi_{1}(F_{X}(1))-\phi_{1}(F_{X}(0))=\dfrac{p_{1+}+p_{2+}}{\sqrt{p_{1+}p_{2+}}}=\big{(}p_{1+}p_{2+}\big{)}^{-1/2}$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\
\psi_{11}-\psi_{10}:=\psi_{1}(F_{Y}(1))-\psi_{1}(F_{Y}(0))=\dfrac{p_{+1}+p_{+2}}{\sqrt{p_{+1}p_{+2}}}=\big{(}p_{+1}p_{+2}\big{)}^{-1/2}.$

Reproducing Goodman’s Association Model. Our discrete copula-based categorical
data model (3.1) expresses the logarithm of “dependence-ratios”333 Goodman
(1996, pp. 410) calls it “Pearson ratios.”

$\operatorname{cop}(F_{X}(k),F_{Y}(l))\,=\,\dfrac{p_{kl}}{p_{k+}p_{+l}}$
(3.12)

as a linear combination of LP-orthonormal row and column scores satisfying
(3.3)-(3.5). The copula-dependence ratio (3.12) measures the strength of
association between the $k$-th row category and $l$-the column category. To
make the connection even more explicit, rewrite (3.1) for two-way contingency
tables as follows:

$\log p_{kl}\,=\,\mu_{0}+\mu_{k}^{{\rm R}}+\mu_{l}^{{\rm
C}}+\sum_{j=1}^{m}\mu_{j}\phi_{jk}\psi_{jl},$ (3.13)

where $\mu_{k}^{{\rm R}}$ denotes the logarithm of row marginal $\log p_{k+}$
and $\mu_{l}^{{\rm C}}$ denotes the logarithm of column marginal $\log
p_{+l}$. Goodman (1991) called this model a “weighted association model” where
weights are marginal row and column proportions. He used the term “association
model” (to distinguish it from correlation (3.6) based model) as it studies
the relationship between rows and columns using odds-ratio.

######  Remark 3.

Log-linear models are a powerful statistical tool for categorical data
analysis (Agresti, 2013). Here we have provided a contemporary unified view of
loglinear modeling for contingency tables from discrete LP-copula viewpoint.
This newfound connection might open up new avenues of research.

###  3.2 Logratio biplot: Graphical Exploratory Analysis

We describe a graphical exploratory tool—logratio biplot, which allows a quick
visual understanding of the relationship between the categorical variables $X$
and $Y$. In the following, we describe the process of constructing logratio
biplot from the LP-copula model (3.1).

Copula-based Algorithm. Construct two scatter plots based on the top two
dominant components of the LP-copula model: the first one is associated with
the row categories, formed by the points $(\mu_{1}\phi_{1k},\mu_{2}\phi_{2k})$
for $k=1,\ldots,I$; and the second one is associated with the column
categories, formed by the points $(\mu_{1}\psi_{1l},\mu_{2}\psi_{2l})$ for
$l=1,\ldots,J$. Logratio biplot is a two-dimensional display obtained by
overlaying these two scatter plots–the prefix ‘bi’ refers to the fact that it
shares a common set of axes for both the rows and columns categories.

Interpretation. Here we offer an intuitive explanation of the logratio biplot
from the copula perspective. We start by recalling the definition of
conditional comparison density (CCD; see eq. 2.8-2.9), as the copula-slice.
For fixed $X=k$, logratio-embedding coordinates $(\mu_{j}\phi_{jk})$ can be
viewed as the LP-Fourier coefficients of the $d(F_{Y}(y);Y,Y|X=k)$, since

$d(F_{Y}(y);Y,Y|X=k)\,=\,\exp\Big{\\{}\mu_{0}+\sum_{j=1}^{m}\big{(}\mu_{j}\phi_{jk}\big{)}\psi_{jl}\Big{\\}},$
(3.14)

Similarly, the logratio coordinates $(\mu_{j}\psi_{jl})$ for fixed $Y=l$ can
be interpreted as the LP-expansion coefficients of $d(F_{X}(x);X,X|Y=l)$.
Hence, the logratio biplot can alternatively be viewed as follows: (i)
estimates the discrete LP-copula density; (ii) Extract the copula slice
$\widehat{d}(u;Y,Y|X=k)$ along with its LP-coefficients
$(\mu_{1}\phi_{1k},\mu_{2}\phi_{2k})$; (iii) similarly, get the estimated
$\widehat{d}(v;X,X|Y=l)$—the copula slice at $Y=F_{Y}(l)$ along with its LP-
coefficients $(\mu_{1}\psi_{1l},\mu_{2}\psi_{2l})$; (iv) Hence, the logratio
biplot (see Fig. 2(b)) measures the association between the row and column
categories $X=k$ and $Y=l$ by measuring the similarity between the ‘shapes’ of
$\widehat{d}(u;Y,Y|X=k)$ and $\widehat{d}(v;X,X|Y=l)$ through their LP-Fourier
coefficients.

######  Remark 4 (Historical Significance).

The following remarks are pertinent: (i) Log-ratio map traditionally taught
and practiced using matrix-algebra (Greenacre, 2018). This is in sharp
contrast with our approach, which has provided a statistical synthesis of log-
ratio biplot from a new copula-theoretic viewpoint. To the best of author’s
knowledge, this is the first work that established such a connection. (ii)
Logratio biplot has some important differences with the correspondence
analysis pioneered by the French statistician Jean-Paul Benzécri; for more
details, see Goodman (1991) and Benzecri (1991). However, in practice, these
two methods often lead to very similar conclusions (e.g., contrast Fig. 2(b)
and Fig. 7).

######  Example 5.

1970 Vietnam War-era US Draft Lottery (Fienberg, 1971). All eligible American
men aged 19 to 26 were drafted through a lottery system in 1970 to fill the
needs of the country’s armed forces. The results of the draft are given to us
in the form of a $12\times 3$ contingency table (see Table 3 in the appendix):
rows are months of the year from January to December, and columns denote three
categories of risk of being drafted—high, medium, and low. The discrete
‘staircase-shaped’ LP-copula estimate is shown in the Fig. 2 (a), whose
explicit form is given below:

$\widehat{\operatorname{cop}}_{X,Y}(u,v;X,Y)\leavevmode\nobreak\
=\leavevmode\nobreak\
\exp\big{\\{}0.26\,\phi_{1}(u;X)\,\psi_{1}(v;Y)+0.18\,\phi_{2}(u;X)\,\psi_{2}(v;Y)\,-\,0.052\big{\\}}.$

We now overlay the scatter plots $(0.26\phi_{1k},0.18\phi_{2k})$ for
$k=1,\ldots,12$ and $(0.26\psi_{1l},0.18\psi_{2l})$ for $l=1,2,3$ to construct
the logratio biplot, as displayed in Fig. 2. This easy-to-interpret two-
dimensional graph captures the essential dependence pattern between the row
(month: in blue dots) and the column (risk category: in red triangles)
variables.

Figure 2: 1970 draft lottery: (a) piecewise-constant nonlinear LP-copula
density estimate, and (b) two-dimensional logratio biplot that essentially
captures all the useful information.

###  3.3 Loglinear Modeling of Large Sparse Contingency Tables

It has been known for a long time (Fienberg and Rinaldo, 2007) that classical
maximum likelihood-based log-linear models break down when applied to large
sparse contingency tables with many zero cells. Here we discuss a new maxent
copula-based smooth method for fitting a parsimonious log-linear model to
sparse contingency tables.

######  Example 6.

Zelterman data. The dataset (Zelterman, 1987, Table 1) is summarized as a
$28\times 26$ cross-classified table that reports monthly salary and number of
years of experience since bachelor’s degree of $n=129$ women employed as
mathematicians or statisticians. The table is extremely sparse—86% cells are
empty! See Fig. 8 of Appendix A.4.

Figure 3: Estimated LP-copula for the $28\times 26$ Zelterman contingency
table data. A strong positive correlation is evident.

A Parsimonious Model. The estimated smooth log-linear LP-copula model for the
Zelterman data is given by:

$\displaystyle\widehat{\operatorname{cop}}_{X,Y}(u,v)\,=\,\exp\Big{\\{}0.52\,\phi_{1}(u;X)\,\psi_{1}(v;Y)+0.37\,\phi_{2}(u;X)\,\psi_{2}(v;Y)+0.18\,\phi_{3}(u;X)\,\psi_{3}(v;Y)\,-\,0.27\Big{\\}},$

displayed in Fig. 3. This shows a strong positive correlation between salary
and number of years of experience. However, the most notable aspect is the
effective model dimension, which can be viewed as the intrinsic degrees of
freedom (df). Our LP-maxent approach distills a compressed representation with
reduced numbers of parameters that yields a smooth estimates: only requiring
$m=3$ components to capture the pattern in the data—a radical compression with
negligible information loss! Contrast this with the dimension of the saturated
loglinear model: $(28-1)\times(26-1)=675$ —a case of a severely
overparameterized non-smooth model with inflated degrees of freedom, which
leads to an inaccurate goodness-of-fit test for checking independence between
rows and columns. More on this in Sec. 3.5.1.

Smoothing Ordered Contingency Tables. The nonparametric maximum likelihood-
based cell probability estimates ${\widetilde{p}}_{X,Y}(k,l)=f_{kl}/n$ are
very noisy and unreliable for sparse contingency tables. By sparse, we mean
tables with a large number of cells relative to the number of observations.

Using Sklar’s representation theorem, one can simply estimate the joint
probability $\widehat{p}_{X,Y}(k,l)$ by multiplying the empirical product pmf
${\widetilde{p}}_{X}(k){\widetilde{p}}_{Y}(l)$ with the smoothed LP-copula. In
other words, the copula can be viewed as a data-adaptive bivariate discrete
kernel function that smooths the independent product-density to estimate the
cell probabilities.

$\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\texttt{dKernel}(k,l)\,=\,\widehat{\operatorname{cop}}_{X,Y}\big{(}\widetilde{F}_{X}(k),\widetilde{F}_{Y}(l)\big{)},\leavevmode\nobreak\
\leavevmode\nobreak\ \text{for}\leavevmode\nobreak\
k=1,\ldots,I;\leavevmode\nobreak\ l=1,\ldots,J.$ (3.15)

where the discrete-kernel function satisfies

$\dfrac{1}{n^{2}}\sum_{k=1}^{I}\sum_{l=1}^{J}\texttt{dKernel}(k,l)f_{k+}f_{+l}\,=\,1.$
(3.16)

This approach can be generalized for any bivariate discrete distribution; see
next section.

######  Remark 5.

For a comprehensive literature on traditional kernel-based nonparametric
smoothing methods for sparse contingency tables, readers are encouraged to
consult Simonoff (1985, 1995) and references therein.

###  3.4 Modeling Bivariate Discrete Distributions

The topic of nonparametric smoothing for multivariate discrete distributions
has received far less attention than the continuous one. Two significant
contributions in this direction include: Aitchison and Aitken (1976) and
Simonoff (1983). In what follows, we discuss a new LP-copula-based procedure
for modeling correlated discrete random variables.

######  Example 7.

Shunter accident data (Arbous and Kerrich, 1951). As a motivating example,
consider the following data: we are given the number of accidents incurred by
$n=122$ shunters in two consecutive year periods, namely 1937-1942 and
1943-1947. To save space, we display the bivariate discrete data in a
contingency table format; see Table 1.

Algorithm. The main steps of our analysis are described below.

Step 1. Modeling marginal distributions. We start by looking at the marginal
distributions of $X$ and $Y$. As seen in Fig. 4, negative binomial
distributions provide excellent fit. To fix the notation, by
$G_{\mu,\phi}={\rm NB}(y;\mu,\phi)$, we mean the following probability
distribution:

${\rm
NB}(y;\mu,\phi)=\binom{y+\phi-1}{y}\,\left(\frac{\mu}{\mu+\phi}\right)^{\\!y}\,\left(\frac{\phi}{\mu+\phi}\right)^{\\!\phi}\\!,\leavevmode\nobreak\
\leavevmode\nobreak\ y\in\mathbb{N},$ (3.17)

where $\mathbb{E}(X)=\mu$ and
$\operatorname{Var}(X)=\mu+\frac{\mu^{2}}{\phi}$. Using the method of MLE, we
get: $X\sim G_{1}={\rm NB}(x;\hat{\mu}=0.97,\hat{\phi}=3.60)$, and $Y\sim
G_{2}={\rm NB}(y;\hat{\mu}=4.30,\hat{\phi}=1.27)$.

Figure 4: Accident Data. Top panel: marginal modeling—comparing the observed empirical pmf with the fitted negative binomial (NB) distribution. Bottom panel shows the estimated smooth joint pmf. | 1943-47
---|---
1937-42 | 0 | 1 | 2 | 3 | 4 | 5 | 6
0 | 21 | 18 | 8 | 2 | 1 | 0 | 0
1 | 13 | 14 | 10 | 1 | 4 | 1 | 0
2 | 4 | 5 | 4 | 2 | 1 | 0 | 1
3 | 2 | 1 | 3 | 2 | 0 | 1 | 0
4 | 0 | 0 | 1 | 1 | 0 | 0 | 0
7 | 0 | 1 | 0 | 0 | 0 | 0 | 0
Table 1: Shunter accidents data displayed compactly as XY-contingency table.

Step 2. Generalized copula density. The probability of bivariate distribution
at $(x,y)$ can be written as follows (generalizing Sklar’s Theorem):

$\Pr(X=x,Y=y):=p_{X,Y}(x,y)=g_{1}(x)g_{2}(y)\,\operatorname{cop}_{X,Y}\big{(}G_{1}(x),G_{2}(y)\big{)}$
(3.18)

where generalized log-copula density admits the following decomposition:

$\log\big{\\{}\operatorname{cop}_{X,Y}(G_{1}(x),G_{2}(y))\big{\\}}\,=\,\sum_{j=1}^{m_{1}}\sum_{k=1}^{m_{2}}\theta_{jk}T_{j}(x;G_{1})T_{k}(y;G_{2})\,-\,\log
Z_{\bm{\theta}}$ (3.19)

It is important to note that the set of LP-basis functions
$\\{T_{j}(x;G_{1})\\}$ and $\\{T_{k}(y;G_{2})\\}$ are specially designed for
the parametric marginals $G_{1}$ and $G_{2}$. We call them gLP-basis, to
distinguish them from the earlier empirical LP-polynomial systems
$\\{T_{j}(x;\widetilde{F}_{X})\\}$ and $\\{T_{k}(y;\widetilde{F}_{Y})\\}$.

Step 3. Exploratory goodness-of-fit. The estimated log-bilinear LP-copula is

$\widehat{\operatorname{cop}}_{X,Y}(u,v)\,=\,\exp\Big{\\{}0.287S_{1}(u;G_{1})S_{1}(v;G_{2})-0.043\Big{\\}}.$
(3.20)

There are three important conclusions that can be drawn from this non-uniform
copula density estimate: (i) Goodness of fit diagnostic: the independence
model (product of parametric marginals) $g_{\perp}(x,y)=g_{1}(x)g_{2}(y)$ is
not adequate for the data. (ii) Nature of discrepancy: the presence of
significant $\hat{\theta}_{11}=0.287$ in the model (3.20) implies that the
tentative independence model should be updated by incorporating the strong
(positive) ‘linear’ correlation between $X$ and $Y$. (iii) Nonparametric
repair: how to update the initial $g_{\perp}(x,y)$ to construct a “better”
model? Eq. (3.18) gives the general updating rule, which simply says: copula
provides the necessary bivariate-correction function to reduce the ‘gap’
between the starting misspecified model $g_{\perp}(x,y)$ and the true unknown
distribution $p_{X,Y}(x,y)$.

Step 4. LP-smoothed probability estimation. The bottom panel Fig. 4 shows the
final smooth probability estimate $\widehat{p}_{X,Y}(x,y)$, computed by
substituting (3.20) into (3.18). Also compare Tables 4 and 5 of Appendix A.4.

######  Remark 6.

Our procedure fits a ‘hybrid’ model: a nonparametrically corrected (through
copula) multivariate parametric density estimate.444A similar philosophy was
proposed in Mukhopadhyay (2017) for univariate continuous distribution case.
One can use any parametric distribution instead of a negative binomial. The
algorithm remains fully automatic, irrespective of the choice of parametric
marginals $G_{1}$ and $G_{2}$, which makes it a universal procedure.

###  3.5 Mutual Information

Mutual information (MI) is a fundamental quantity in Statistics and Machine
Learning, with wide-ranging applications from neuroscience to physics to
biology. For continuous random variables $(X,Y)$, mutual information is
defined as

$\operatorname{MI}(X,Y)\,=\,\iint
f_{X,Y}(x,y)\log\dfrac{f_{X,Y}(x,y)}{f_{X}(x)f_{Y}(y)}\;\mathrm{d}x\;\mathrm{d}y.$
(3.21)

Among non-parametric MI estimators, $k$-nearest-neighbor and kernel-density-
based methods (Moon et al., 1995, Kraskov et al., 2004, Zeng et al., 2018) are
undoubtedly the most popular ones. Here we are concerned with a slightly
general problem of developing a flexible MI estimation algorithm that is: (D1)
applicable for mixed555Reliably estimating MI for mixed case is notoriously
challenging task (Gao et al., 2017). $(X,Y)$; (D2) robust in the presence of
noise; and, (D3) invariant under monotone transformations666This is essential
to make the analysis less sensitive to various types of data preprocessing,
which is done routinely in applications like bioinformatics, astronomy, and
neuroscience.. To achieve this goal, we start by rewriting MI (3.21) using
copula:

$\operatorname{MI}(X,Y)\,=\,\int_{[0,1]^{2}}\operatorname{cop}_{X,Y}(u,v)\log\operatorname{cop}_{X,Y}(u,v)\;\mathrm{d}u\;\mathrm{d}v.$
(3.22)

The next theorem presents an elegant closed-form expression for MI in terms of
LP-copula parameters, which allows a fast and efficient estimation algorithm.

######  Theorem 3.

Let $(X,Y)$ be a mixed-pair of random variables. Under the LP log-bilinear
copula model (2.11), the mutual information between $X$ and $Y$ has the
following representation in terms of LP-co-mean parameters
$\operatorname{LP}_{jk}$ and maximum entropy coefficients $\theta_{jk}$

$\operatorname{MI}_{\bm{\theta}}(X,Y)\,=\,\mathop{\sum\sum}_{j,k>0}\theta_{jk}\operatorname{LP}_{jk}\,-\,\log
Z_{\bm{\theta}}$ (3.23)

Proof. Express mutual information as:

$\operatorname{MI}_{\bm{\theta}}(X,Y)\,=\,\mathbb{E}_{X,Y}\big{[}\log\operatorname{cop}_{\bm{\theta}}\big{]}\,=\,\sum_{j}\sum_{k}\theta_{jk}\mathbb{E}_{X,Y}\big{[}S_{j}(U;X)S_{k}(V;Y)\big{]}\,-\,\log
Z_{\bm{\theta}}$ (3.24)

The first equality follows from (3.22) and the second one from (2.11).
Complete the proof by replacing by $\operatorname{LP}_{jk}$ by
$\mathbb{E}[S_{j}(U;X)S_{k}(V;Y)]$ by virtue of (2.13). As a practical
consequence, we have the following efficient and direct MI-estimator,
satisfying D1-D3:

$\widehat{\rm
MI}_{\bm{\theta}}(X,Y)\,=\,\mathop{\sum\sum}_{j,k>0}\widehat{\theta}_{jk}{\widetilde{\operatorname{LP}}}_{jk}\,-\,\log
Z_{\widehat{\bm{\theta}}}.\vspace{-.5em}$ (3.25)

Continuous $(X,Y)$ example. Consider the kidney fitness data, discussed in
Example 1. The estimated LP-copula density-based (2.22) MI estimate:
$\widehat{\operatorname{MI}}=0.20$ with bootstrap pvalue is almost zero.

######  Remark 7.

MI (3.22) measures the departure of copula density from uniformity. This is
because, MI can be viewed as the Kullback-Leibler (KL) divergence between
copula and the uniform density: ${\rm MI}(X,Y)={\rm
KL}(\operatorname{cop};U_{[0,1]^{2}})$. A few immediate consequences: (i) MI
is always nonnegative, i.e., ${\rm MI}(X,Y)\geq 0$, and equality holds if and
only if variables are independent. Moreover, the stronger the dependence
between two variables, the larger the MI. (ii) MI is also invariant under
different marginalizations. Two additional applications of MI (for categorical
data and feature selection problems) are presented below.

####  3.5.1 Application 1: (X,Y) Discrete: Smooth-G2 Statistic

Given $n$ independent samples from an $I\times J$ contingency table, the
$G^{2}$-test of goodness-of-fit, also known as the log-likelihood ratio test,
is defined as

$G^{2}(X,Y)\,=\,2n\sum_{k=1}^{I}\sum_{l=1}^{J}\,{\widetilde{p}}_{kl}\log\dfrac{{\widetilde{p}}_{kl}}{{\widetilde{p}}_{k+}{\widetilde{p}}_{+l}},$
(3.26)

which under the null hypothesis of independence has asymptotic
$\chi^{2}_{(I-1)(J-1)}$ distribution. From (3.26) one can immediately conclude
the following.

######  Theorem 4.

The $G^{2}$ log-likelihood ratio statistic can be viewed as the raw
nonparametric MI-estimate

$G^{2}(X,Y)/2n\,=\,\widetilde{\operatorname{MI}}(X,Y),\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ $ (3.27)

where $\widetilde{\operatorname{MI}}(X,Y)$ is obtained by replacing the
unknown distributions in (3.21) with their empirical estimates.

######  Example 8.

Hellman’s Infant Data (Yates, 1934) We verify the identity (3.27) for the
following $2\times 2$ table, which shows cross-tabulation of $n=42$ infants
based on whether the infant was breast-fed or bottle-fed.

| Normal teeth | Malocclusion
---|---|---
Breast-fed | 4 | 16
Bottle-fed | 1 | 21
Table 2: The data table on malocclusion of the teeth in infants were obtained
by M. Hellman and reported in the classic paper by Frank Yates (1934, p.230).

The estimated LP-copula density (shown in Fig. 10 of Appendix A.4) is given by

$\widehat{\operatorname{cop}}(u,v;X,Y)=\exp\big{\\{}0.234S_{1}(u;X)S_{1}(v;Y)\,-\,0.03\big{\\}}.$
(3.28)

The empirical MI estimate is given by
$2n\times\widetilde{\operatorname{MI}}(X,Y)=2.50$, with the pvalue $0.12$
computed using the asymptotic null distribution $\chi^{2}_{1}$. This exactly
matches with the $G^{2}$-statistic value; one may use the R-function GTest.

The problem arises when we try to apply $G^{2}$-test for large sparse tables,
and it is not hard to see why: the adequacy of asymptotic
$\chi^{2}_{(I-1)(J-1)}$ distribution depends both on the sample size $n$ and
the number of cells $p=IJ$. Koehler (1986) showed that the approximation
completely breaks down when $n/p<5$, leading to erroneous statistical
inference due to significant loss of power. The following example demonstrates
this.

######  Example 9.

Zelterman Data Continued. Log-likelihood ratio $G^{2}$-test produces pvalue
$1$, firmly concluding the independence between salary and years of
experience. This directly contradicts our analysis of Sec. 3.3, where we found
a clear positive dependence between these two variables. Why $G^{2}$-test was
unable to detect that effect? Because it is based on chi-square approximation
with degrees of freedom $(28-1)\times(26-1)=675$. This inflated degrees of
freedom completely ruined the power of the test. To address this problem, we
recommend the following smoothed version:

${\rm
Smooth-}G^{2}(X,Y)/2n=\widehat{\operatorname{MI}}(X,Y),\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ $ (3.29)

where $\widehat{\operatorname{MI}}(X,Y)$ is computed based on the LP-bilinear
copula model:

$\widehat{\operatorname{cop}}(u,v;Y,X)=\exp\Big{\\{}0.50\,S_{1}(u;Y)\,S_{1}(v;X)+0.18\,S_{2}(u;Y)S_{2}(v;X)\,-\,0.167\Big{\\}},$

${\rm Smooth-}G^{2}$ analysis generates pvalue $2.48\times 10^{-11}$, thereby
successfully detecting the association. The crucial aspect of our approach
lies in its ability to provide a reduced dimensional parametrization of copula
density. For the Zelterman data, we need just two components (i.e., the
effective degrees of freedom is $2$) to capture the pattern.

######  Remark 8 (Discrete variables with many categories).

Discrete distributions over large domains routinely arise in large-scale
biomedical data such as diagnosis codes, drug compounds and genotypes (Seok
and Kang, 2015). The method proposed here can be used to jointly model such
random variables.

####  3.5.2 Application 2: (X,Y) Mixed: Feature Importance Score

We consider the two-sample feature selection problem where $Y$ is a binary
response variable, and $X$ is a predictor variable that can be either discrete
or continuous.

Figure 5: CKD data: The top panel shows the estimated copula density. The
bottom panel shows the two-sample boxplots and the conditional comparison
density (CCD) $d(v;X,X|Y=1)$. The piecewise constant shape of estimated CCD is
not an aberration of our nonparametric approximation method; it reflects the
inherent discreteness of the feature $X$ (red blood cell count).

######  Example 10.

Chronic Kidney Disease data. The goal of this study is to investigate whether
kidney function is related to red blood cell count (RBC). We have a sample of
$n=203$ participants, among whom $79$ have chronic kidney disease (ckd) and
another $124$ are non-ckd. $Y$ denotes the kidney disease status and $X$
denotes the measurements on RBC (unit in million cells per cubic millimeter of
blood). The estimated LP-copula density

$\displaystyle\widehat{\operatorname{cop}}(u,v)=\exp\Big{\\{}-0.76\,S_{1}(u;Y)\,S_{1}(v;X)+0.18\,S_{1}(u;Y)S_{2}(v;X)-0.19\,S_{1}(u;Y)S_{4}(v;X)\,-\,0.33\Big{\\}},$

(3.30)

is shown in Fig. 5. We propose mutual information-based feature importance
measure based on the formula (3.25); this yields
$\widehat{\operatorname{MI}}(Y,X)=0.36$ with pvalue almost zero, strongly
indicating that RBC is an importance risk-factor related to kidney
dysfunction.

What additionally insights can we glean from this copula? To answer that
question, let’s focus our attention on the copula-slice for $u\in[0.61,1]$.
This segment of the copula is essentially the conditional comparison density
$d(v;X,X|Y=1)$ (see Sec. 2.1), which can be easily derived by substituting
$S_{1}(\widetilde{F}_{Y}(1);Y)=\sqrt{\frac{1-\tilde{\mu}}{\tilde{\mu}}}=1.25$
into (3.30), where $\tilde{\mu}=79/203=0.389$:

$\widehat{d}(v;X,X|Y=1)=\exp\big{\\{}-0.95S_{1}(v;X)+0.23S_{2}(v;X)-0.24S_{4}(v;X)-0.33\big{\\}}.$
(3.31)

A few remarks on the interpretation of the above formula:

$\bullet$ Distributional effect-size: Bearing in mind eq. (2.5, 2.9), it is
easy to see that $d(v;X,X|Y=1)$ compares two densities: $f_{X|Y=1}(x)$ with
$f_{X}(x)$, thereby capturing the distributional difference. This has
advantages over traditional two-sample feature importance statistic (e.g.,
Student’s t or Wilcoxon statistic) that can only measure differences in
location or mean.

$\bullet$ Explainability: The estimated $\widehat{d}(v;X,X|Y=1)$ involves
three significant LP-components of $X$; the presence of the 1st order ‘linear’
$S_{1}(v;X)$ indicates location-difference; the 2nd order ‘quadratic’
$S_{2}(v;X)$ indicates scale-difference; and 4th order ‘quartic’ $S_{4}(v;X)$
indicates the presence of tail-difference in the two RBC-distributions. In
addition, the negative sign of the linear effect $\widehat{\theta}_{11}=-0.95$
implies reduced mean level of RBC in the ckd-population. Medically, this makes
complete sense, since a dysfunctional kidney cannot produce enough
Erythropoietin (EPO) hormone, which causes the RBC to drop.

###  3.6 Nonparametric Copula-Logistic Regression

We describe a new copula-based nonparametric logistic regression model. The
key result is given by the following theorem, which provides a first-principle
derivation of a robust nonlinear generalization of the classical linear
logistic regression model.

######  Theorem 5.

Let $\mu=\Pr(Y=1)$ and $\mu(x)=\Pr(Y=1|X=x)$. Then we have the following
expression for the logit (log-odds) probability model:

$\operatorname{logit}\left\\{\mu(x)\right\\}\,=\,\log\left(\dfrac{\mu(x)}{1-\mu(x)}\right)\,=\,\alpha_{0}\,+\sum_{j}\alpha_{j}T_{j}(x;F_{X}),\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ $ (3.32)

where $\alpha_{0}=\operatorname{logit}(\mu)$ and
$\alpha_{j}=\frac{\theta_{j1}}{\sqrt{\mu(1-\mu)}}$.

Proof. The proof consists of four main steps.

Step 1. To begin with, notice that for $Y$ binary and $X$ continuous, the
general log-bilinear LP-copula density function (2.11) reduces to the
following form:

$\operatorname{cop}_{\bm{\theta}}(u,v;X,Y)\leavevmode\nobreak\
=\leavevmode\nobreak\
\dfrac{1}{Z_{\theta}}\exp\Big{\\{}\sum_{j}\theta_{j1}S_{j}(u;X)S_{1}(v;Y)\Big{\\}},$
(3.33)

since we can construct at most $2-1=1$ LP-basis function for binary $Y$.

Step 2. Apply copula-based Bayes Theorem (Eq. 2.10) and express the
conditional comparison densities as follows:

$d_{1}(x)\,\equiv\,d(F_{X}(x);X,X|Y=1)\,=\,\dfrac{\mu(x)}{\mu},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ $ (3.34)

and also,

$d_{0}(x)\,\equiv\,d(F_{X}(x);X,X|Y=0)\,=\,\dfrac{1-\mu(x)}{1-\mu}.$ (3.35)

Taking logarithm of the ratio of (3.34) and (3.35), we get the following
important identity:

$\log\left(\dfrac{\mu(x)}{1-\mu(x)}\right)\,=\,\log\left(\dfrac{\mu}{1-\mu}\right)\,+\,\log
d_{1}(x)\,-\,\log d_{0}(x).$ (3.36)

Step 3. From (3.33), one can deduce the following orthonormal expansion of
maxent-conditional copula slices $\log d_{1}$ and $\log d_{0}$:

$\displaystyle\log d_{1}(x)$ $\displaystyle=$
$\displaystyle\sum_{j}\big{(}\theta_{j1}T_{1}(1;F_{Y})\big{)}\,T_{j}(x;F_{X})\,-\,\log
Z_{\theta}$ (3.37) $\displaystyle\log d_{0}(x)$ $\displaystyle=$
$\displaystyle\sum_{j}\big{(}\theta_{j1}T_{1}(0;F_{Y})\big{)}\,T_{j}(x;F_{X})\,-\,\log
Z_{\theta}$ (3.38)

Step 4. Substituting (3.37) and (3.38) into (3.36) we get:

$\log\left(\dfrac{\mu(x)}{1-\mu(x)}\right)\,=\,\log\left(\dfrac{\mu}{1-\mu}\right)\,+\,\sum_{j}\Big{\\{}\frac{\theta_{j1}}{\sqrt{\mu(1-\mu)}}\Big{\\}}T_{j}(x;F_{X}),$

since for binary $Y$ we have (see Appendix A.1):

$T_{1}(1;F_{Y})-T_{1}(0;F_{Y})\,=\,\dfrac{1-\mu}{\sqrt{\mu(1-\mu)}}\,+\,\dfrac{\mu}{\sqrt{\mu(1-\mu)}}\,=\,\dfrac{1}{\sqrt{\mu(1-\mu)}}.$
(3.39)

Substitute $\alpha_{0}=\operatorname{logit}(\mu)$ and
$\alpha_{j}=\frac{\theta_{j1}}{\sqrt{\mu(1-\mu)}}$ to complete the proof. ∎

Figure 6: UCI credit card data: We demonstrate the predictive performance and
explainability of the LP-copula-based additive logistic regression model. The
boxplots of the accuracy (AUC) are shown in the top left panel. On average,
LP-logistic regression provides a 12% boost in the accuracy. The LS-plot is
shown in the top left. For easy interpretation, we display the scaled the LS-
plot:
$(\frac{\widehat{\alpha}_{j1}}{max_{j}\widehat{\alpha}_{j1}},\frac{\widehat{\alpha}_{j2}}{max_{j}\widehat{\alpha}_{j2}})$.
We can see a tight cluster around the origin $(0,0)$, which indicates that
most of the variables are irrelevant for prediction (sparsity-assumption).
Also, it is evident that the majority of the features have either location or
scale information (differences), the only exception being the variable
pay$\\_$0—which denotes the repayment status of the last two months (-1=pay
duly, 1=payment delay for one month, 2=payment delay for two months, and so
on). This is further illustrated using three variables, as marked in the LS-
plot; see also the two-sample density estimates shown at the bottom panel.

####  3.6.1 High-dimensional Copula-based Additives Logistic Regression

Generalize the univariate copula-logistic regression model (3.32) to the high-
dimensional case as follows:

$\operatorname{logit}\big{(}\mu(x)\big{)}\,=\,\alpha_{0}+\sum_{j=1}^{p}h_{j}(x_{j}).$
(3.40)

Nonparametrically approximate the unknown smooth $h_{j}$’s by LP-polynomial
series of $X_{j}$

$h_{j}(x_{j})\,=\,\sum_{k=1}^{m}\alpha_{jk}T_{k}(x_{j};F_{X_{j}}),\leavevmode\nobreak\
\text{for}\,j=1,\ldots,p.$ (3.41)

######  Remark 9 (Estimation and Computation).

A sparse nonparametric LP-additive model of the form (3.40-3.41) can be
estimated by penalized regression techniques (lasso, elastic net, etc.), whose
R-computation is remarkably easy using glmnet function (Friedman et al.,
2010):

$\texttt{glmnet}(y\sim T_{\bf X},\text{family={binomial}})\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ $

where $T_{\bf X}$ simply is the column-wise stacked $[T_{X_{1}}\mid\cdots\mid
T_{X_{p}}]$ LP-feature matrix.

######  Example 11.

UCI Credit Card data. The dataset is available in the UCI Machine Learning
Repository. It contains records of $n=30,000$ cardholders from an important
Taiwan-based bank. For each customer, we have a response variable $Y$
denoting: default payment status (Yes = 1, No = 0), along with $p=23$
predictor variables (e.g., gender, education, age, history of past payment,
etc.). We randomly partition the data into training and test sets, with an
80-20 split, repeated $100$ times. We measure the prediction accuracy using
AUC (the area under the ROC curve). Fig. 6 compares two kinds of lasso-
logistic regressions: (i) usual version: based on feature matrix $X$; and (ii)
LP-copula version: based on feature matrix $T_{{\bf X}}$. As we can see, LP-
copula based additive logistic regression classifier significantly outperforms
the classical logistic regression model. To gain further insight into the
nature of impact of each variable, we plot the lasso-smoothed location and
scale coefficients:

$\text{LS-Feature plot}:\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\big{(}\widehat{\alpha}_{j1},\,\widehat{\alpha}_{j2}\big{)},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ {\rm
for}\leavevmode\nobreak\ j=1,\ldots,p.\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ $

L-stands for location and S-stands for scale. The purpose of the LS-feature
plot is to characterize ‘how’ each feature impacts the classification task.
For example, consider the three variables pay$\\_$0, limit$\\_$balance, and
pay$\\_$6, shown in the bottom panel Fig. 6. Each one of them contains unique
discriminatory information: pay$\\_$0 has location as well as scale
information, hence it appeared at the top-right of the LS-plot; The variable
limit$\\_$balance only shows location difference, whereas the variable
pay$\\_$6 shows contrasting scale in the two populations. In short, LS-plot
explains ‘why and how’ each variable is important using a compact diagram,
which is easy to interpret by researchers and practitioners.

##  4 Conclusion

This paper makes the following contributions: (i) we introduce modern
statistical theory and principles for maximum entropy copula density
estimation that is auto-adaptive for the mixed(X,Y)—described in Section 2.
(ii) Our general copula-based formulation provides a unifying framework of
data analysis from which one can systematically distill a number of
statistical methods by revealing some completely unexpected connections
between them. The importance of our theory in applied and theoretical
statistics is highlighted in Section 3, taking examples from different sub-
fields of statistics: Log-linear analysis of categorical data, logratio
biplot, smoothing large sparse contingency tables, mutual information,
smooth-$G^{2}$ statistic, feature selection, and copula-based logistic
regression. We hope that this new perspective will motivate new methodological
developments for cohesive statistical data analysis.

### Dedication: Two Legends from Two Different Cultures

This paper is dedicated to the birth centenary of E. T. Jaynes (1922–1998),
the originator of the maximum entropy principle.

I also like to dedicate this paper to the memory of Leo Goodman (1928–2020)—a
transformative legend of categorical data analysis, who passed away on
December 22, 2020, at the age of 92 due to COVID-19.

This paper is inspired in part by the author’s intention to demonstrate how
these two modeling philosophies can be connected and united in some ways. This
is achieved by employing a new nonparametric representation theory of
generalized copula density.

## References

  * Agresti (2013) Agresti, A. (2013).  Categorical data analysis (3rd ed.).  John Wiley & Sons.
  * Aitchison and Aitken (1976) Aitchison, J. and C. G. Aitken (1976).  Multivariate binary discrimination by the kernel method.  Biometrika 63(3), 413–420.
  * Arbous and Kerrich (1951) Arbous, A. G. and J. Kerrich (1951).  Accident statistics and the concept of accident-proneness.  Biometrics 7(4), 340–432.
  * Barron and Sheu (1991) Barron, A. R. and C. Sheu (1991).  Approximation of density functions by sequences of exponential families.  Annals of Statistics 19, 1347–1369.
  * Benzecri (1991) Benzecri, J. P. (1991).  Comment on Leo Goodman’s “Models, and graphical displays in the analysis of cross-classified data”.  Journal of the American Statistical Association. 86(416), 1112–1115.
  * Efron and Hastie (2016) Efron, B. and T. Hastie (2016).  Computer Age Statistical Inference, Volume 5.  Cambridge University Press.
  * Fienberg (1971) Fienberg, S. E. (1971).  Randomization and social affairs: the 1970 draft lottery.  Science 171(3968), 255–261.
  * Fienberg and Rinaldo (2007) Fienberg, S. E. and A. Rinaldo (2007).  Three centuries of categorical data analysis: Log-linear models and maximum likelihood estimation.  Journal of Statistical Planning and Inference 137(11), 3430–3445.
  * Friedman et al. (2010) Friedman, J., T. Hastie, and R. Tibshirani (2010).  Regularization paths for generalized linear models via coordinate descent.  Journal of statistical software 33(1), 1.
  * Gao et al. (2017) Gao, W., S. Kannan, S. Oh, and P. Viswanath (2017).  Estimating mutual information for discrete-continuous mixtures.  In Advances in neural information processing systems, pp. 5986–5997.
  * Goodman (1991) Goodman, L. A. (1991).  Measures, models, and graphical displays in the analysis of cross-classified data (with discussion).  Journal of the American Statistical association 86(416), 1085–1111.
  * Goodman (1996) Goodman, L. A. (1996).  A single general method for the analysis of cross-classified data: reconciliation and synthesis of some methods of Pearson, Yule, and Fisher, and also some methods of correspondence analysis and association analysis.  Journal of the American Statistical Association 91(433), 408–428.
  * Greenacre (2018) Greenacre, M. (2018).  Compositional data analysis in practice.  CRC Press.
  * Hoeffding (1940) Hoeffding, W. (1940).  Massstabinvariante korrelationstheorie.  Schriften des Mathematischen Seminars und des Instituts f$\ddot{u}$r Angewandte Mathematik der Universit$\ddot{a}$t Berlin 5(3), 179–233.
  * Jaynes (1957) Jaynes, E. T. (1957).  Information theory and statistical mechanics.  Physical review 106(4), 620.
  * Koehler (1986) Koehler, K. J. (1986).  Goodness-of-fit tests for log-linear models in sparse contingency tables.  Journal of the American Statistical Association 81(394), 483–493.
  * Kraskov et al. (2004) Kraskov, A., H. Stögbauer, and P. Grassberger (2004).  Estimating mutual information.  Physical review E 69(6), 066138.
  * Moon et al. (1995) Moon, Y.-I., B. Rajagopalan, and U. Lall (1995).  Estimation of mutual information using kernel density estimators.  Physical Review E 52(3), 2318.
  * Mukhopadhyay (2017) Mukhopadhyay, S. (2017).  Large-scale mode identification and data-driven sciences.  Electronic Journal of Statistics 11(1), 215–240.
  * Mukhopadhyay and Parzen (2020) Mukhopadhyay, S. and E. Parzen (2020).  Nonparametric universal copula modeling.  Applied Stochastic Models in Business and Industry, special issue on “Data Science” 36(1), 77–94.
  * Mukhopadhyay et al. (2012) Mukhopadhyay, S., E. Parzen, and S. Lahiri (2012).  From data to constraints.  Bayesian Inference And Maximum Entropy Methods In Science And Engineering: 31st International Workshop, Waterloo, Canada 1443(1), 32–39.
  * Parzen and Mukhopadhyay (2013) Parzen, E. and S. Mukhopadhyay (2013).  United statistical algorithms, LP-comoment, copula density, nonparametric modeling.  59th ISI World Statistics Congress (WSC), Hong Kong, 4719–4724.
  * Seok and Kang (2015) Seok, J. and Y. S. Kang (2015).  Mutual information between discrete variables with many categories using recursive adaptive partitioning.  Scientific reports 5(1), 1–10.
  * Simonoff (1983) Simonoff, J. S. (1983).  A penalty function approach to smoothing large sparse contingency tables.  The Annals of Statistics, 208–218.
  * Simonoff (1985) Simonoff, J. S. (1985).  An improved goodness-of-fit statistic for sparse multinomials.  Journal of the American Statistical Association 80(391), 671–677.
  * Simonoff (1995) Simonoff, J. S. (1995).  Smoothing categorical data.  Journal of Statistical Planning and Inference 47(1), 41–69.
  * Sklar (1959) Sklar, M. (1959).  Fonctions de répartition à n dimensions et leurs marges.  Publ. Inst. Statistique Univ. Paris 8, 229–231.
  * Yates (1934) Yates, F. (1934).  Contingency tables involving small numbers and the $\chi^{2}$ test.  Supplement to the Journal of the Royal Statistical Society 1(2), 217–235.
  * Zelterman (1987) Zelterman, D. (1987).  Goodness-of-fit tests for large sparse multinomial distributions.  Journal of the American Statistical Association 82(398), 624–629.
  * Zeng et al. (2018) Zeng, X., Y. Xia, and H. Tong (2018).  Jackknife approach to the estimation of mutual information.  Proceedings of the National Academy of Sciences 115(40), 9956–9961.

##  5 Supplementary Appendix

###  A.1 Nonparametric LP-Polynomials

Preliminaries. For a random variable $X$ with the associated probability
distribution $F_{X}$, define the mid-distribution function as
$F^{\rm{mid}}(x;F_{X})=F_{X}(x)-\frac{1}{2}p(x;F_{X})$ where $p(x;F_{X})$ is
probability mass function. The $F^{\rm{mid}}(X;F_{X})$ has mean
$\mathbb{E}[F^{\rm{mid}}(X;F_{X})]=.5$ and variance
$\operatorname{Var}[F^{\rm{mid}}(X;F_{X})]=\frac{1}{12}\big{(}1-\sum_{x}p^{3}(x;F_{X})\big{)}$.
Define first-order basis function $T_{1}(X;F_{X})$ by standardizing
$F^{\rm{mid}}$-transformed random variable:

$T_{1}(x;F_{X})\leavevmode\nobreak\ =\leavevmode\nobreak\
\dfrac{\sqrt{12}\big{\\{}F^{\rm{mid}}(x;F_{X})-1/2\big{\\}}}{\sqrt{1-\sum_{x}p^{3}(x;F_{X})}}.$
(6.1)

Construct the higher-order LP-polynomial bases777Here the number of LP-basis
functions $m$ is always less than $|\mathscr{U}|$, where $\mathscr{U}$ denotes
the set of all unique values of $X$. $\\{T_{j}(X;F_{X})\\}_{1\geq j\leq m}$ by
Gram-Schmidt orthonormalization of $T_{1}^{2},T_{1}^{3},\ldots,T_{1}^{m}$. We
call these specially-designed polynomials of mid-distribution transforms as
LP-basis, which by construction, satisfy the following orthonormality with
respect to the measure $F_{X}$:

$\int_{x}T_{j}(x;F_{X})\;\mathrm{d}F_{X}(x)=0,\leavevmode\nobreak\
\leavevmode\nobreak\ \text{and}\leavevmode\nobreak\ \leavevmode\nobreak\
\int_{x}T_{j}(x;F_{X})T_{k}(x;F_{X})\;\mathrm{d}F_{X}(x)=\delta_{jk},$ (6.2)

where $\delta_{jk}$ is the Kronecker delta. Two particular LP-family of
polynomials are given below: $\bullet$ gLP-Polynomial Basis. Let $X\sim G$,
where $G$ is a known distribution. Following the above recipe, construct
parametric gLP-basis $\\{T_{j}(X;G)\\}_{j\geq 1}$ for the given distribution
$G$. For $X\sim{\rm Bernoulli}(p)$, one can show that
$T_{1}(x;G)=\frac{x-p}{\sqrt{p(1-p)}},$ $x=0,1;$ see Mukhopadhyay and Parzen
(2020) for more details.

$\bullet$ eLP-Polynomial Basis. In practice, we only have access to a random
sample $X_{1},\ldots,X_{n}$ from an unknown distribution $F$. To perform
statistical data analysis, it thus becomes necessary to nonparametrically
‘learn’ an appropriate basis that is orthogonal with respect to the (discrete)
empirical measure $\widetilde{F}_{X}$. To address that need, construct LP-
basis with respect to the empirical measure
$\\{T_{j}(X;\widetilde{F}_{X})\\}_{j\geq 1}$. LP-unit Basis. We will
occasionally express the $T_{j}$’s in the quantile domain

$S_{j}(u;X)\,=\,T_{j}\big{(}Q_{X}(u);F_{X}\big{)},\leavevmode\nobreak\ 0<u<1.$
(6.3)

We call these S-functions the unit LP-bases. The ‘S-form’ and the ‘T-form’
will be used interchangeably throughout the paper, depending on the context.

LP-product-bases and Copula Approximation. The LP-product-bases
$\big{\\{}T_{j}(x;F_{X})T_{k}(y;F_{Y})\big{\\}}$ can be used to expand any
square-integrable function of the form $\Psi(F_{X}(x),F_{Y}(y))$. In
particular, the logarithm of Hoeffding’s dependence function (2.3) can be
approximated by LP-Fourier series:

$\log\operatorname{cop}_{X,Y}\big{(}F_{X}(x),F_{Y}(y)\big{)}\leavevmode\nobreak\
\sim\leavevmode\nobreak\ \sum_{j,k}\theta_{j,k}T_{j}(x;F_{X})T_{k}(y;F_{Y}).$
(6.4)

Represent (6.4) in the quantile domain by substituting $x=Q_{X}(u)$ and
$y=Q_{Y}(v)$ to get the copula density expression (2.11).

###  A.2 Two Cultures of Maximum Entropy Modeling

The form of the maximum-entropy exponential model directly depends on the form
of the set of constraints, i.e., the sufficient statistics functions. The
maxent distribution depends on the data only through the sample averages for
these functions.

$\bullet$ Parametric maxent modeling culture: The traditional practice of
maxent density modeling assumes that the appropriate sufficient statistics
functions are known or given beforehand, which, in turn, puts restrictions on
the possible ‘shape’ of the probability distribution. This parametric maxent
modeling culture was first established by Ludwig Boltzmann in 1877, and then
later popularized by E. T. Jaynes in 1960s.

$\bullet$ Nonparametric maxent modeling culture: It proceeds by identifying a
small set of most important sufficient statistics functions from data
(Mukhopadhyay et al., 2012). In the next step, we build the maxent probability
distribution that agrees with these specially-designed relevant constraints.
We have used LP-orthogonal polynomials to systematically and robustly design
the constraining functions. See, Section A.3 for more discussion.

###  A.3 Connection With Vladimir Vapnik’s Statistical Invariants

Vladimir Vapnik (2020) calls the sufficient statistics functions as
‘predicates’ and the associated moment constraints as ‘statistical
invariants.’ While describing his learning theory he acknowledged that

> “The only remaining question in the complete statistical learning theory is
> how to choose a (small) set of predicates. The choice of predicate functions
> reflects the intellectual part of the learning problem.”

The question of how to systematically design and search for informative
predicates was previously raised by Mukhopadhyay et al. (2012) in the context
of learning maxent probability models from data; also see Section 2.2 of the
main article where we discussed some concrete strategies to address this
problem. However, it is important to note that Vapnik and Izmailov (2020)
works entirely within the traditional least-square ($L_{2}$ risk) setup,
instead of maximum-entropy framework.

## References

  * Mukhopadhyay and Parzen (2020) Mukhopadhyay, S. and E. Parzen (2020).  Nonparametric universal copula modeling.  Applied Stochastic Models in Business and Industry, special issue on “Data Science” 36(1), 77–94.
  * Mukhopadhyay et al. (2012) S. Mukhopadhyay, E. Parzen. and S. N. Lahiri (2012).  From data to constraints.  Bayesian Inference And Maximum Entropy Methods In Science And Engineering: 31st International Workshop, Waterloo, Canada 1443, 32–39.
  * Vapnik and Izmailov (2020) Vapnik, V. and R. Izmailov (2020).  Complete statistical theory of learning: learning using statistical invariants.  In A. Gammerman, V. Vovk, Z. Luo, E. Smirnov, and G. Cherubin (Eds.), Proceedings of the Ninth Symposium on Conformal and Probabilistic Prediction and Applications, Volume 128 of Proceedings of Machine Learning Research, pp. 4–40.

###  A.4 Additional Figures and Tables

Figure 7: The classical correspondence analysis plot for the 1970 draft
lottery data; implemented using the R-package ca. Contrast it with Fig. 2(b).

Figure 8: Zelterman data is a $28\times 26$ table, which is extremely sparse.
Here we display the histogram of the observed cell frequencies. For more
details see Example 6.

Figure 9: CKD data. Figure 10: LP-copula density estimate for Hellman $2\times 2$ data. | High | Med | Low
---|---|---|---
Jan | 9 | 12 | 10
Feb | 7 | 12 | 10
Mar | 5 | 10 | 16
Apr | 8 | 8 | 14
May | 9 | 7 | 15
Jun | 11 | 7 | 12
Jul | 12 | 7 | 12
Aug | 13 | 7 | 11
Sep | 10 | 15 | 5
Oct | 9 | 15 | 7
Nov | 12 | 12 | 6
Dec | 17 | 10 | 4
Table 3: 1970 US draft lottery data; see Example 5. | 0 | 1 | 2 | 3 | 4 | 5 | 6
---|---|---|---|---|---|---|---
0 | 0.17 | 0.15 | 0.07 | 0.02 | 0.01 | 0.00 | 0.00
1 | 0.11 | 0.12 | 0.08 | 0.01 | 0.03 | 0.01 | 0.00
2 | 0.03 | 0.04 | 0.03 | 0.02 | 0.01 | 0.00 | 0.01
3 | 0.02 | 0.01 | 0.02 | 0.02 | 0.00 | 0.01 | 0.00
4 | 0.00 | 0.00 | 0.01 | 0.01 | 0.00 | 0.00 | 0.00
7 | 0.00 | 0.01 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00
Table 4: Accident data. The raw empirical probability estimates $\widetilde{p}(x,y)$. | 0 | 1 | 2 | 3 | 4 | 5 | 6
---|---|---|---|---|---|---|---
0 | 0.19 | 0.13 | 0.06 | 0.02 | 0.01 | 0.00 | 0.00
1 | 0.09 | 0.10 | 0.07 | 0.03 | 0.01 | 0.01 | 0.00
2 | 0.03 | 0.05 | 0.04 | 0.02 | 0.01 | 0.00 | 0.00
3 | 0.01 | 0.02 | 0.02 | 0.01 | 0.01 | 0.00 | 0.00
4 | 0.00 | 0.01 | 0.01 | 0.00 | 0.00 | 0.00 | 0.00
7 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00
Table 5: Accident data. The LP-smoothed probability estimates
$\widehat{p}(x,y)$; see Example 7.