Density Sharpening: Principles and Applications to Discrete Data Analysis
Deep Mukhopadhyay
<EMAIL_ADDRESS>

###### Abstract

This article introduces a general statistical modeling principle called
“Density Sharpening” and applies it to the analysis of discrete count data.
The underlying foundation is based on a new theory of nonparametric
approximation and smoothing methods for discrete distributions which play a
useful role in explaining and uniting a large class of applied statistical
methods. The proposed modeling framework is illustrated using several real
applications, from seismology to healthcare to physics.

Keywords: Density sharpening; $\operatorname{DS}(p_{0},m)$ distributions; LP-
Fourier analysis; Explanatory goodness-of-fit; Jaynes’ dice problem;
Compressive $\chi^{2}$; Data-efficient learning.

###### Contents

  1. 1 Principles of Statistical Model Building
  2. 2 Density Sharpening: Model and Mechanism
    1. 2.1 Learning by Comparison: $d$-Sharp Density
    2. 2.2 LP-Fourier Analysis
    3. 2.3 The $\operatorname{DS}(p_{0},m)$ Model
    4. 2.4 LP-Maximum Entropy Analysis
  3. 3 Applications in Statistical Modelling
    1. 3.1 One-sample Test of Proportion
    2. 3.2 Expandable Negative Binomial Distribution
    3. 3.3 $\chi^{2}$ and Compressive-$\chi^{2}$
    4. 3.4 Explanatory Goodness-of-Fit
    5. 3.5 Relative Entropy and Model Uncertainty
    6. 3.6 Card Shuffling Problem
    7. 3.7 Jaynes Dice Problem
    8. 3.8 Compressive Learning of Big Distributions
    9. 3.9 Data-efficient Statistical Learning
    10. 3.10 Discovery-source Separation Problem
  4. 4 Discussion
  5. 5 Appendix
    1. A.1 Goodness-of-fit Methods: Automatic verses Manual
    2. A.2 LPdiscrete: A Unified Algorithm
    3. A.3 Sub-linear Algorithms for Detecting Sparse Structure
    4. A.4 Additional Figures

##  1 Principles of Statistical Model Building

> ‘Part of a meaningful quantitative analysis is to look at models and try to
> figure out their deficiencies and the ways in which they can be improved.’
>
> —-Nobel Lecture by Lars Peter Hansen (2014)

Scientific investigation never happens in a vacuum. It builds upon previously
accumulated knowledge instead of starting from scratch. Statistical modeling
is no exception to this rule.

Suppose we are given $n$ random samples $X_{1},\ldots,X_{n}$ from an unknown
discrete distribution $p(x)$. Before we jump into the statistical analysis
part, the scientist provided us a hint on what might be an initial believable
model for the data: ‘from my years of experience working in this field, I
expect the underlying distribution to be somewhat close to $p_{0}(x)$.’ This
information came with a disclaimer: ‘don’t take $p_{0}(x)$ too seriously as it
is only a simplified approximation of reality. Use it with caution and care.’

The general problem of statistical learning then aims to address the following
questions: Whether the ‘shape of the data’ is consistent with the presumed
model-0. If it is not, then what is it? How is it different from $p_{0}$?
Revealing new hidden pattern in the data is often the most essential
statistical modeling task in science and engineering. Of course, ultimately,
the aim is to search for a rich class of sensible models in an automatic and
faster manner, by appropriately changing the misspecified $p_{0}$. Knowing how
to change the anticipated $p_{0}$ is the first step towards scientific
discovery that allows scientists to re-evaluate alternative theories to
explain the data. If we succeed at this, it will provide a mechanism to build
“hybrid” knowledge-data integrated models, which are far more interpretable
than classical fully data-driven nonparametric models. Full development of
these ideas requires a new conceptual framework and mathematical tools.

Organization. Section 2 introduces a new family of nonparametric approximation
and smoothing techniques for discrete probability distributions, which is
built on the principle of ‘Density Sharpening.’ Section 3 highlights the role
of the proposed theoretical framework in the development of statistical
methods that is rich enough to include traditional as well as contemporary
statistical methods: starting from as simple as one sample Z-test for a
proportion to as sophisticated as compressive chi-square, $d$-sharp negative
Binomial distribution, universal goodness-of-fit program, relative entropy
estimation, Jaynes dice problem, sample-efficient learning of big
distributions, etc. The paper ends with a discussion and conclusion Section 4.
Additional applications and methodological details are deferred to the
Supplementary Appendix to ensure the smooth flow of the main ideas.

##  2 Density Sharpening: Model and Mechanism

We describe a method of nonparametric approximation of discrete distribution
by sharpening the initially assumed $p_{0}(x)$. The theory is remarkably
simple, yet general enough to be vastly applicable in many areas beyond
density estimation, as described in Section 3. Here is a bird’s eye view of
the core mechanism, which is a three-stage process. Stage 1. Model-0
elicitation: The modeler starts a suitable $p_{0}(x)$ by using his/her
experience or subject-matter knowledge. Often a particular parametric form of
$p_{0}(x)$ is selected keeping convenience and simplicity in mind.

Stage 2. Exploratory uncertainty analysis: Assess the uncertainty of the
presumed model $p_{0}(x)$, in a way that can explain ‘why and how’ the assumed
model-0 (i.e., $p_{0}$) is inadequate for the data.

Stage 3. Coarse-to-Refined density: Incorporate the ‘learned’ uncertainty into
$p_{0}(x)$ to produce an improved model $\widehat{p}(x)$ that will eliminate
the incompatibility with the data.

The required theory is developed in the next few sections, which heavily
relies on the following notation: let $X$ be a discrete variable with
probability mass function $p_{0}(x)$, cumulative distribution function
$F_{0}(x)$, and mid-distribution function $F_{0}^{{\rm
mid}}(x)=F_{0}(x)-\frac{1}{2}p_{0}(x)$. The associated quantile function will
be denoted by $Q_{0}(u)=\inf\\{x:F_{0}(x)\geq u\\}$ for $0<u<1$. By
$\mathcal{L}^{2}(dF_{0})$ we mean the set of all square integrable functions
with respect to the discrete measure $\;\mathrm{d}F_{0}$, i.e, for a function
$\psi\in\mathcal{L}^{2}(dF_{0})$:
$\int|\psi|^{2}\;\mathrm{d}F_{0}:=\sum_{x}|\psi(x)|^{2}p_{0}(x)<\infty$. The
inner product of two functions $\psi_{1}$ and $\psi_{2}$ in
$\mathcal{L}^{2}(dF_{0})$ will be denoted by
$\langle\psi_{1},\psi_{2}\rangle_{F_{0}}:=\int\psi_{1}\psi_{2}\;\mathrm{d}F_{0}$.
Expectation with respect to $p_{0}(x)$ will be abbreviated as
$\mathbb{E}_{0}(\psi(X)):=\int\psi\;\mathrm{d}F_{0}$.

###  2.1 Learning by Comparison: $d$-Sharp Density

We introduce a mechanism for nonparametrically estimating the density of
$X_{1},\ldots,X_{n}$ by comparing and sharpening the presumed working model
$p_{0}(x)$.

######  Definition 1 ($d$-Sharp Density).

For a discrete random variable $X$, we have the following universal density
decomposition:

$p(x)\,=\,p_{0}(x)\,d\big{(}F_{0}(x);F_{0},F\big{)},$ (2.1)

where the $d(u;F_{0},F)$ is defined as

$d(u;F_{0},F)=\dfrac{p(Q_{0}(u))}{p_{0}(Q_{0}(u))},\leavevmode\nobreak\
\,0<u<1.$ (2.2)

The function $d(u;F_{0},F)$ is called ‘comparison density’ because it compares
the assumed $p_{0}$ with the true $p(x)$ and it integrates to one:

$\int_{0}^{1}d(u;F_{0},F)\;\mathrm{d}u\,=\,\int_{x}d(F_{0}(x);F_{0},F)\;\mathrm{d}F_{0}(x)\,=\,\sum_{x}\big{(}p(x)/p_{0}(x)\big{)}p_{0}(x)\,=\,1.\leavevmode\nobreak\
\leavevmode\nobreak\ $

For brevity’s sake, we will often abbreviate $d(F_{0}(x);F_{0},F)$ as
$d_{0}(x)$ throughout the article.

######  Remark 1 (The philosophy of ‘learning by comparison’).

The density representation formula (2.1) describes a way of building a general
$p(x)$ by comparing it with the initial $p_{0}(x)$. The $d$-modulated class of
distributions is constructed by amending (instead of abandoning) the starting
imprecise model $p_{0}(x)$.

######  Remark 2 ($d$-sharp density).

Eq. (2.1) provides a formal statistical mechanism for sharpening the initial
vague $p_{0}(x)$ using the data-guided perturbation function $d_{0}(x)$. For
this reason, we call the improved $p_{0}(x)\times d_{0}(x)$ the ‘$d$-sharp’
density.

######  Example 1 (Earthquake Data).

We are given annual counts of major earthquakes (magnitude 6 and above) for
the years 1900-2006. It is available in the R package astsa. Seismic engineers
routinely use negative binomial distribution for modeling earthquake frequency
(Kagan and Jackson, 2000, Kagan, 2010). The best fitted negative binomial (NB)
with ($\mu=19$ and $\phi=12$) is shown in Fig 1, which we take as our rough
initial model $p_{0}(x)$. From the figure, it is clearly evident that the
conventional NB distribution is unable to adequately capture the shape of the
earthquake count data. Model uncertainty quantification. For earthquake
engineers it is of utmost importance to determine the uncertainty of the
assumed NB model (Bernreuter, 1981). The problem of uncertainty quantification
is the holy grail of earthquake science due to its importance in estimating
risk and in making an accurate forecast of the next big earthquake. The
comparison density $d(u;F_{0},F)$ captures the uncertainty of the assumed NB
model. The left plot in Fig 1 displays the estimated $\widehat{d}(u;F_{0},F)$
for this data, which reveals the nature of deficiency of the base NB model. A
robust nonparametric method of estimating $\widehat{d}$ from data will be
discussed in the subsequent sections. But before going into the nonparametric
approximation theory, we will spend some time on its interpretation.

Figure 1: Modeling the earthquakes distribution. Left: Estimated comparison
density $\widehat{d}(u;F_{0},F)$; Right: The fitted NB distribution and the
re-calibrated $d$-sharpened version.

Interpretable Exploratory learning. In our model (2.1), $d$ plays the role of
a data-driven correction function, measuring the discrepancy between the
initial $p_{0}$ and the unknown $p$. Thus, the non-uniformity of $\widehat{d}$
immediately tells us that there’s something more in the data than what was
expected in light of $p_{0}(x)$. In fact, the shape of $\widehat{d}$ reveals
the nature of the most prominent deviations between the data and the presumed
$p_{0}(x)$—which, in this case, are bimodality and presence of heavier-tail
than anticipated NB distribution.

######  Remark 3 (Role of $d$).

The comparison density $d$ performs dual functions: (i) its graph acts as an
exploratory diagnostic tool that exposes the unexpected, forcing decision
makers (e.g., legislators, natural security agencies, local administrations)
to think outside the box: what might have caused this bimodality? how can we
repair the old seismic hazard forecast model so that it incorporates this new
information? etc. (ii) it provides a formal process of transforming and
revising an initially misspecified model into a useful one. The red curve in
the right panel is obtained by multiplying (perturbing) the NB pmf with the
estimated comparison density, obeying the density representation formula Eq.
(2.1).

###  2.2 LP-Fourier Analysis

The task is to nonparametrically approximate $d(F_{0}(x);F_{0},F)$ to be able
to apply the density sharpening equation (2.1). We approximate $d\circ
F_{0}(x)\in\mathcal{L}^{2}({dF_{0}})$ by projecting it into a space of
polynomials of $F_{0}(x)$ that are orthonormal with respect to the base
measure$\;\mathrm{d}F_{0}$. How to construct such a system of polynomials in a
completely automatic and robust manner for any given $p_{0}(x)$? In the
section that follows, we discuss a universal construction.

####  2.2.1 Discrete LP-Basis

We describe a general theory of constructing LP-polynomials—a new class of
robust polynomials $\\{T_{j}(x;F_{0})\\}_{j\geq 1}$ that are a function of
$F_{0}(x)$ (not raw $x$) and are orthonormal with respect to user-specified
discrete distribution $p_{0}(x)$.

Figure 2: Shape of the top four LP-orthonormal basis functions
$\\{S_{j}(u;F_{0})\\}_{1\leq j\leq 4}$ for $p_{0}(x):={\rm
NB}(x;\mu=19,\phi=12)$ for the earthquake data; where recall $S_{j}\circ
F_{0}\,(x)=T_{j}(x;F_{0})$. Notice of the global nonlinearity (linear,
quadratic, cubic and so on) and the local staircase-like (piecewise-constant
unequal-length segments) shape of these specially-designed polynomials. They
are, by construction, orthonormal with respect to the chosen measure
$p_{0}(x)$, here ${\rm NB}(x;\mu=19,\phi=12)$.

Step 1: Define the first-order LP-basis function as standardized mid-
distribution transform:

$T_{1}(x;F_{0})\,=\,\dfrac{\sqrt{12}\big{[}F_{0}^{{\rm
mid}}(x)-0.5\big{]}}{\sqrt{1-\sum_{x}p_{0}^{3}(x)}}.$ (2.3)

Verify that $\mathbb{E}_{0}[T_{1}(X;F_{0})]=0$ and
$\mathbb{E}_{0}[|T_{1}(X;F_{0})|^{2}]=1$, since $\mathbb{E}[F_{0}^{{\rm
mid}}(X)]=1/2$ and $\operatorname{Var}[F_{0}^{{\rm
mid}}(X)]=\sqrt{(1-\sum_{x}p_{0}^{3}(x))}/12$.

Step 2: Apply a _weighted_ Gram-Schmidt procedure on $\\{T_{1}^{2},\ldots
T_{1}^{k-1}\\}$ to construct a higher-order LP orthogonal system
$T_{j}(x;F_{0})$ with respect to measure$\;\mathrm{d}F_{0}$

$\sum\nolimits_{x}p_{0}(x)T_{j}(x;F_{0})=0;\leavevmode\nobreak\
\leavevmode\nobreak\
\,\sum\nolimits_{x}p_{0}(x)T_{j}(x;F_{0})T_{k}(x;F_{0})=\delta_{jk},\leavevmode\nobreak\
\leavevmode\nobreak\ 1<j,k<M$

where $\delta_{jk}$ is the Kronecker delta function and the highest-degree of
the LP-polynomials $M$ is always less than the support size of the discrete
$p_{0}$. For example, if $X$ is binary, one can construct at most $2-1=1$ LP-
basis function; see Section 3.1. Fig 2 shows the top four LP-basis functions
for the earthquake data with $p_{0}$ as ${\rm NB}(x;\mu=19,\phi=12)$. Here, we
have displayed them in a unit interval as a function of $u=F_{0}(x)$, denoted
by $S_{j}(u;F_{0}):=T_{j}(Q_{0}(u);F_{0}),0<u<1$. Notice the typical shape of
these custom-constructed discrete orthonormal polynomials: globally nonlinear
(linear, quadratic, cubic, and so on) and locally piecewise-constant with
unequal step size.

######  Remark 4 (Role of LP-coordinate system in unification of statistical
methods).

LP-bases play a unique role in statistical modeling—they provide an efficient
coordinate (data-representation) system that is fundamental to developing
unified statistical algorithms.

###  2.3 The $\operatorname{DS}(p_{0},m)$ Model

######  Definition 2 (LP-canonical Expansion).

Expand comparison density in the LP-orthogonal series

$d(F_{0}(x);F_{0},F)\,=\,1+\sum_{j}\operatorname{LP}[j;F_{0},F]\,T_{j}(x;F_{0}),$
(2.4)

where the $j$th LP-Fourier coefficient satisfies the following identity:

$\operatorname{LP}[j;F_{0},F]=\big{\langle}d\circ
F_{0},T_{j}\big{\rangle}_{F_{0}}.$ (2.5)

A change-of-basis perspective: The conventional way to represent a discrete
distribution is through indicator basis (histogram representation):

$\eta_{j}(i)\,=\,\mathbbm{1}\big{\\{}X_{i}\in[x_{j},x_{j+1})\big{\\}},\leavevmode\nobreak\
\leavevmode\nobreak\ {\rm for}\,j=1,2,\ldots,r\leavevmode\nobreak\ $ (2.6)

where $r$ is the domain size (number of unique values) of the empirical
distribution ${\widetilde{p}}(x)$. In (2.4), we have performed a “change of
basis” from the amorphous indicator-basis to a more structured LP-basis
$\\{T_{j}(x;F_{0})\\}$, where the expansion coefficients
$\operatorname{LP}[j;F_{0},F]$ act as the coordinates of $p(x)$ relative to
assumed $p_{0}(x)$:

$\big{[}F\big{]}_{F_{0}}:=\Big{(}\operatorname{LP}[1;F_{0},F],\ldots,\operatorname{LP}[m;F_{0},F]\Big{)},\leavevmode\nobreak\
\leavevmode\nobreak\ 1\leq m<r.$

For that reason, one may call these coefficients the discrete LP-Fourier
Transform (LPT) of $p(x)$ relative to $p_{0}(x)$.

######  Definition 3.

$\operatorname{DS}(p_{0},m)$ denotes a class of distributions with the
following representation:

$p(x)\,=\,p_{0}(x)\Big{[}1\,+\,\sum_{j=1}^{m}\operatorname{LP}[j;F_{0},F]\,T_{j}(x;F_{0})\Big{]},$
(2.7)

obtained by replacing (2.4) into (2.1). Here $\operatorname{DS}(p_{0},m)$
stands for Density-Sharpening of $p_{0}(x)$ using $m$-term LP-series
approximated $d_{0}(x)$. $\operatorname{DS}(p_{0},m)$ is a class of
nonparametrically-designed parametric models that are flexible enough to
capture various _shapes of discrete_ $p(x)$, like multi-modality, excess-
variation, long-tailed, and sharp peaks.

To estimate the unknown coefficients $\operatorname{LP}[j;F_{0},F]$ of the
$\operatorname{DS}(p_{0},m)$ model, note the following important identity:

$\displaystyle\operatorname{LP}[j;F_{0},F]$ $\displaystyle=$
$\displaystyle\int d(F_{0}(x);F_{0},F)T_{j}(x;F_{0})\;\mathrm{d}F_{0}(x)$
(2.8) $\displaystyle=$ $\displaystyle\int T_{j}(x;F_{0})\;\mathrm{d}F(x)$
$\displaystyle=$ $\displaystyle\mathbb{E}_{F}\big{[}T_{j}(X;F_{0})\big{]}.$

This immediately leads to the following “weighted mean” estimator:

$\widetilde{\operatorname{LP}}_{j}\,:=\operatorname{LP}[j;F_{0};\widetilde{F}]\,=\,\mathbb{E}_{\widetilde{F}}\big{[}T_{j}(X;F_{0})\big{]}\,=\,\sum_{x}{\widetilde{p}}(x)T_{j}(x;F_{0}),\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ $ (2.9)

Using standard empirical process theory (Csörgő, 1983, Parzen, 1998) one can
show that the limiting distribution of sample LP-statistic is i.i.d
$\mathcal{N}(0,n^{-1/2})$, under the null hypothesis $H_{0}:p=p_{0}$. Thus one
can quickly obtain a sparse estimated $\operatorname{DS}(p_{0},m)$ model by
retaining only the ‘significant’ LP-coefficients, which are greater than
$2/\sqrt{n}$.

Earthquake Data Example. The first $m=10$ estimated
$\widetilde{\operatorname{LP}}_{j}$ are shown in Fig. 13 of the appendix,
which indicates that the only interesting non-zero LP-coefficient is
$\widetilde{\operatorname{LP}}_{6}$. The explicit form of the estimated
$\operatorname{DS}({\rm NB},m=6)$ model for the earthquake data is given by:

$\widehat{p}(x)=p_{0}(x)\big{[}1+0.20T_{6}(x;F_{0})\big{]},$ (2.10)

where $p_{0}={\rm NB}(x;\mu=19,\phi=12)$. The resulting $\widehat{p}(x)$ is
plotted as a red curve in Fig. 1.

###  2.4 LP-Maximum Entropy Analysis

To ensure non-negativity, we expand $\log d$ (instead of $d$ as we have done
in Eq. (2.4)) in LP-Fourier series, which results in the following exponential
model:

$d_{\bm{\theta}}(u;F_{0},F)\,=\,\exp\Big{\\{}\sum_{j\geq
1}\theta_{j}S_{j}(u;F_{0})\,-\,\Psi(\bm{\theta})\Big{\\}},\leavevmode\nobreak\
\leavevmode\nobreak\ 0<u<1$ (2.11)

where
$\Psi(\bm{\theta})=\log\int_{0}^{1}\exp\\{\sum_{j}\theta_{j}S_{j}(u;F_{0})\\}\;\mathrm{d}u.$
This model is also called the maximum-entropy (maxent) comparison density
model because it maximizes the entropy $-\int d_{\bm{\theta}}\log
d_{\bm{\theta}}$ (flattest possible; thus promotes smoothness) under the
following LP-moment constraints:

$\mathbb{E}_{\bm{\theta}}[S_{j}(U;F_{0})]\,=\,\operatorname{LP}[j;F_{0},\widetilde{F}],\leavevmode\nobreak\
\leavevmode\nobreak\ (j=1,2\ldots).$ (2.12)

LP-moments $\operatorname{LP}[j;F_{0},\widetilde{F}]$ are ‘compressed
measurements’ (linear combinations of observed data; verify from (2.9)), which
are sufficient statistics for the comparison density $d_{\bm{\theta}}$.

######  Definition 4 (Maxent $\operatorname{DS}(p_{0},m)$ model).

Replacing (2.11) into (2.1), we have the following maxent
$\operatorname{DS}(p_{0},m)$ model

$p(x)\,=\,p_{0}(x)\exp\Big{\\{}\sum_{j\geq
1}\theta_{j}T_{j}(x;F_{0})\,-\,\Psi(\bm{\theta})\Big{\\}},\vspace{-.5em}$
(2.13)

To estimate a sparse maxent comparison density model, we carry out the
optimization routine by choosing only the ‘significant’ LP-moments in (2.12).

Earthquake Data Example. The estimated maxent $\operatorname{DS}(p_{0},m)$
model for the earthquake distribution is given by

$\widehat{\widehat{p\,}}(x)=p_{0}(x)\exp\big{\\{}0.195T_{6}(x;F_{0})-0.02\big{\\}},$
(2.14)

whose shape is almost indistinguishable from the LP-Fourier estimated p.m.f.
(2.10).

##  3 Applications in Statistical Modelling

We describe how the general principle of ‘density sharpening’ acts as a
unified framework for the analysis of discrete data with a wide variety of
applications, ranging from basic introductory methods to more advanced
statistical modeling techniques.

###  3.1 One-sample Test of Proportion

Given $n$ samples from a binary $X$, the one-sample proportion test is
concerned with testing whether the population proportion $p$ is equal to the
hypothesized proportion $p_{0}$. We approach this problem by reformulating it
in our mathematical notation:

Step 1. We start with the $\operatorname{DS}(p_{0},m=1)$ model

$p(x)=p_{0}(x)\Big{\\{}1+\operatorname{LP}[1;p_{0},p]T_{1}(x;F_{0})\Big{\\}},\leavevmode\nobreak\
\leavevmode\nobreak\ x=0,1.$ (3.1)

where the null model $p_{0}(x)=xp_{0}+(1-x)(1-p_{0}),$ for $x=0,1$.

Step 2. We rewrite the original hypothesis $H_{0}:p=p_{0}$ in terms of LP-
parameter as $H^{\prime}_{0}:\operatorname{LP}[1;p_{0},p]=0$.

Step 3. We derive an explicit formula for
$\operatorname{LP}[1;p_{0},{\widetilde{p}}]$. It’s a two step process: First,
we need the analytic expression of the LP-basis $T_{1}(x;p_{0})$

$T_{1}(x;p_{0})=\left\\{\begin{array}[]{rl}-\dfrac{p_{0}}{\sqrt{p_{0}(1-p_{0})}}&\mbox{for
$x=0$}\\\ \dfrac{1-p_{0}}{\sqrt{p_{0}(1-p_{0})}}&\mbox{for
$x=1$.}\end{array}\right.$ (3.2)

We then apply formula (2.9) to deduce:

$\operatorname{LP}(1;p_{0};{\widetilde{p}})=(1-{\widetilde{p}})T_{1}(0;p_{0})+{\widetilde{p}}T_{1}(1;p_{0})=\dfrac{\widetilde{p}-p_{0}}{\sqrt{p_{0}(1-p_{0})}}.\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ $ (3.3)

Step 4. A remarkable fact is that the test based on (3.3) exactly matches with
the classical Z-test, whose null distribution is:
$\sqrt{n}\operatorname{LP}[1;p_{0},{\widetilde{p}}]\sim\mathcal{N}(0,1)$ as
$n\rightarrow\infty$. This shows how the LP-theoretical device provides a
transparent first-principle derivation of the one-sample proportion test, by
shedding light on its genesis.

###  3.2 Expandable Negative Binomial Distribution

We will focus now on one important special case of maxent
$\operatorname{DS}(p_{0},m)$ family of distributions (2.13), where the base
measure $p_{0}(x)$ is taken to be a negative Binomial (NB) distribution.

$p(x)\,=\,\binom{x+\phi-1}{x}\,\left(\frac{\mu}{\mu+\phi}\right)^{\\!x}\,\left(\frac{\phi}{\mu+\phi}\right)^{\\!\phi}\\!\exp\Big{\\{}\sum_{j\geq
1}\theta_{j}T_{j}(x;F_{0})\,-\,\Psi(\bm{\theta})\Big{\\}},\leavevmode\nobreak\
\leavevmode\nobreak\ x\in\mathbb{N},\leavevmode\nobreak\ \leavevmode\nobreak\
$ (3.4)

note that the basis functions $\\{T_{j}(x;F_{0})\\}_{j\geq 1}$ are specially-
designed LP-orthonormal polynomials associated with the base measure
$p_{0}={\rm NB}(\phi,\mu)$. We call (3.4) the $m$th-order expandable NB
distributions, denoted by XNB(m). A few practical advantages of XNB(m)
distributions are: the computational ease they afford for estimating
parameters; their compactly parameterizable yet shape-flexible nature; and,
finally, their ability to provide explanatory insights into how $p(x)$ is
different from the standard NB distribution. Due to their simplicity and
flexibility, they have the potential to be a ‘default choice’ for modeling
count data. We have already seen an example of XNB-distribution in (2.14) in
the context of modeling the earthquake distribution. We now turn our attention
to two further real-data examples.

Figure 3: Nonparametric modeling of NMES 1988 data. Top left: For clarity, we
only show the data over the domain $[0,20]$. Otherwise, the ultra long-
tailedness of the data will make it harder to see what’s going on in the
interesting zones; the arrows point out some specific aspects of the data that
were missed by NB distribution. Top right: compares two boxplots—one is the
observed data and the other is the simulated data from ${\rm NB}(1,5.7)$,
which captures the difference in the tail behavior. The tail of NB has to be
stretched to match the shape of the data. Bottom left: $\widehat{d}$ captures
the ‘unexplained shape of the data,’ which tells us how to change the initial
model-0 to fit the data. Bottom right: The red curve is the rectified XNB
probability distribution.

######  Example 2 (NMES 1988 Data).

This is a part of the US National Medical Expenditure Survey (NMES) conducted
in 1987 and 1988. It is available in the R-package AER. We have $n=4,406$
observations of a discrete random variable $X$ which denotes how many times an
individual, aged 66 and covered by Medicare, visited physician’s office. As
displayed in the Fig. 3 boxplot, the distribution has a large support size
(varies between $0$ and $89$), with some regions being extremely data-sparse.

The blue curve in the top left plot shows the ${\rm
NB}(\hat{\phi}=1,\hat{\mu}=5.7)$, where the parameters are maximum-likelihood
estimates. Next, we estimate the LP-maxent $\widehat{d}_{\bm{\theta}}$, using
the theory of Sec. 2.4. At this point, it is strongly advisable to pay
attention to the shape $\widehat{d}_{\bm{\theta}}$. Why? Because, it
efficiently extracts and exposes ‘unanticipated’ aspects in the data that
cannot be explained by the initial NB distribution. The bottom-left Fig. 3
immediately reveals a few things: (i) NB underestimates the probability at
$x=0$; (ii) it overestimates the probability mass around $x=2$ and $3$; (iii)
there seems to be an excess probability mass (‘bump’ structure) around $x=4$;
(iv) NB clearly has a shorter tail than what we see in the data—this can be
seen from the sharply increasing right tail of the comparison density. To
better understand the tail-behavior, we have simulated $n$ samples from ${\rm
NB}(1,5.7)$ and contrasted the two boxplots in the top-right panel, which
strongly indicates the long-tailedness of $p(x)$ relative to NB distribution.
Any reader will agree that without the help of $\widehat{d}_{\bm{\theta}}$,
even experienced eyes could have easily missed these subtle patterns. Finally,
multiply $\widehat{d}_{\bm{\theta}}$ by the ${\rm NB}(1,5.7)$, following eq.
(3.4), to get the estimated XNB distribution—the red p.m.f curve, shown in the
bottom-right panel of Fig. 3.

######  Example 3 (Computer Breaks Data).

We are given the number of times a DEC-20 computer broke down at Open
University in each of $n=128$ consecutive weeks of operation, starting in late
1983. The data shows positive skewness with a slightly longer tail; see Fig.
16 in the appendix. The mechanics of XNB modeling proceed as follows: (i) We
start by estimating the parameters of $p_{0}(x)$, which in this case are MLE-
fitted (one can use any other method of estimation) ${\rm
NB}(\hat{\phi}=1.7,\hat{\mu}=4)$. (ii) The next step is estimation of
$\widehat{d}_{\bm{\theta}}$, which in this case is just the uniform
distribution—none of the LP-maxent parameters were large enough to be
selected. This is depicted in the left panel of supplementary Fig. 16. This
graphical diagnostic indicates that the initial $p_{0}$ fits the data
satisfactorily; no repairing is necessary. (iii) Accordingly, our ‘density
sharpening’ principle returns the $\texttt{XNB}(m=0)$ as the final model,
which is simply the starting parametric model ${\rm
NB}(\hat{\phi}=1.7,\hat{\mu}=4)$.

It is interesting to contrast our finding with Saulo et al. (2020), where the
authors fit a highly specialized nonparametric discrete distribution to this
data. The beauty of our approach is that it performs nonparametric correction
(through $d$) only when it is warranted. When the reality is already simple,
we don’t complicate it unnecessarily.

###  3.3 $\chi^{2}$ and Compressive-$\chi^{2}$

Given a random sample of size $n$ from the unknown population distribution
$p(x)$, chi-square goodness-of-fit statistic between the sample probabilities
${\widetilde{p}}(x)$ and the expected $p_{0}(x)$ can be re-written as follows:

$\dfrac{\chi^{2}}{n}=\sum_{x}\dfrac{\big{(}{\widetilde{p}}(x)-p_{0}(x)\big{)}^{2}}{p_{0}(x)}=\sum_{x}p_{0}(x)\big{[}{\widetilde{p}}(x)/p_{o}(x)\,-1\big{]}^{2}=\int_{0}^{1}\big{[}d(u;p_{0},{\widetilde{p}})-1\big{]}^{2}\;\mathrm{d}u,$
(3.5)

By applying Parseval’s identity on the LP-Fourier expansion of $d$, we have
the following important equality:

$\dfrac{\chi^{2}}{n}\,=\,\sum_{j=1}^{r-1}\Big{|}\operatorname{LP}[j;p_{0},{\widetilde{p}}]\Big{|}^{2}\mathrel{\vbox{\hbox{\scriptsize.}\hbox{\scriptsize.}}}=\,\operatorname{LP}(p_{0}\|{\widetilde{p}}),$
(3.6)

where $r$ is the number of unique values in our sample $X_{1},\ldots,X_{n}$.
This shows that chi-square information statistic is a “saturated” raw-
nonparametric measure with $r-1$ components.

######  Example 4 (The Gambler’s Die).

A gambler rolls a die $n=60$ times and gets the following observed counts:

Table 1: The observed frequencies Number on die | $1$ | $2$ | $3$ | $4$ | $5$ | $6$
---|---|---|---|---|---|---
Observed ${\widetilde{p}}$ | 4/60 | 6/60 | 17/60 | 16/60 | 8/60 | 9/60
Hypothesized $p_{0}$ | 1/6 | 1/6 | 1/6 | 1/6 | 1/6 | 1/6

The gambler wishes to determine whether the die is fair. If it is fair, we
would expect the outcomes of $1$ to $6$ are equally likely, with probability
$1/6$. Pearsonian chi-square and a full-rank LP-analysis, both lead to the
same answer:

$\chi_{{\rm
obs}}^{2}\,=\,n\times\sum_{j=1}^{6-1}\Big{|}\operatorname{LP}[j;p_{0},{\widetilde{p}}]\Big{|}^{2}=14.2,\leavevmode\nobreak\
\,\text{with degrees of freedom $5$}$

with the resulting $p$-value $0.0143$. Note that the sum of squares of the
$6-1=5$ LP-Fourier coefficients “exactly” reproduces (numerically) the
Pearson’s chisquare statistic! This further verifies the mathematical fact
elucidated in eq. (3.6). Conclusion: The die is loaded at 5% significance
level.

Figure 4: The estimated $\widehat{d}(u;F_{0},F)$ for examples 4 and 5. It
helps to identify ‘where’ the interesting differences between data
(${\widetilde{p}}$) and hypothesized model ($p_{0}$) lie.

Exploratory insight. Here we want to go beyond classical confirmatory
analysis, with the goal to understand how the die is loaded. The answer is
hidden in the shape of the comparison density $\widehat{d}$. Fig. 4(a) firmly
suggests that the die was loaded heavily in the middle —especially on the
sides $3$ and $4$, where it landed most frequently.

######  Example 5 (Sparse Dice problem).

It is an example of sparse count data with many groups. Imagine a $k=20$
dimensional dice is rolled $n=20$ times:

  * •

The hypothesized model: $p_{0}=(1/4,1/4,1/36,\ldots,1/36)$

  * •

The observed probabilities: ${\widetilde{p}}=(3/4,1/4,0,\ldots,0)$.

We would like to know whether the postulated model $p_{0}(x)$ actually
reflects the data or not. If it does not, then we want to know how the
hypothesized model differs from the observed probabilities.

Pearsonian chi-square111R-function chisq.test() generates the message that
“Chi-squared approximation may be incorrect!.” yields value $\chi_{{\rm
obs}}^{2}=30$, with degrees of freedom $19$ and $p$-value $0.052$. Conclusion:
there is no evidence of discrepancy at 5% level, even though there is a
glaring difference between ${\widetilde{p}}(1)=3/4$ and $p_{0}(1)=1/4$. The
legacy $\chi^{2}$ loses its power because of ‘inflated degrees of freedom’ for
large sparse problems: Large value of $k$ increases the critical value
$\chi^{2}_{\alpha;k-1}$, making it harder to detect ‘small’ but important
changes.

LP-Analysis of the sparse dice problem: (i) Construct the discrete LP-
polynomials $\\{T_{j}(x;p_{0})\\}$ that are specially-designed for the given
$p_{0}(x)$. Appendix Fig. 14 displays the shape of those basis functions. (ii)
Compute the LP-Fourier coefficients
$\operatorname{LP}[j;p_{0},\,{\widetilde{p}}]$ by
$\sum_{x}{\widetilde{p}}(x)T_{j}(x;p_{0})$. The Appx. figure 15 identifies the
first two LP-parameters as significant components. (iii) We now compute the
compressive-$\chi^{2}$ based on these interesting components:

$\operatorname{LP}(p_{0}\|{\widetilde{p}})\,=\,\big{|}\operatorname{LP}[1;p_{0},\,{\widetilde{p}}]\big{|}^{2}\,+\,\big{|}\operatorname{LP}[2;p_{0},\,{\widetilde{p}}]\big{|}^{2}=1.49,\leavevmode\nobreak\
\,\leavevmode\nobreak\ \text{with degrees of freedom $2$}.$ (3.7)

and $p$-value $3.4\times 10^{-7}$. Also noteworthy is the fact that
compressive LP-chisquare is numerically almost same as the raw $\chi^{2}$:

$\vspace{-1em}\chi^{2}_{{\rm obs}}\,=\,30\leavevmode\nobreak\
\approx\leavevmode\nobreak\
n\times\operatorname{LP}(p_{0}\|{\widetilde{p}})=29.8.$

######  Remark 5 (Auto-adaptability).

LP-goodness-of-fit shows an impressive adaptability property: under the usual
scenario (like in example 4) it reduces to the classical $\chi^{2}$ analysis,
and for large-sparse problems, it automatically produces a chi-square
statistic with the fewest possible degrees of freedom, which boosts its power.
Our reformulation (in terms of modern LP-nonparametric language) allowed a
better way of doing chi-square goodness-of-fit analysis that applies to a much
broader class of applied statistics problems. In John Tukey’s (1954) words:
“Do we need to find new techniques, or to use old ones better?”

######  Remark 6 (Ungrouped case).

What if we have ungrouped data: given a random sample of counts
$X_{1},\ldots,X_{n}$, check (confirm) whether the data is compatible with the
hypothesised $p_{0}(x)$, i.e., to test the hypothesis $H_{0}:p=p_{0}$. One way
to approach this problem is to forcefully group the data points into different
categories and then apply Pearson’s $\chi^{2}$ test on it. This is (almost)
always a bad strategy, since grouping leaks information. Moreover, the
arbitrariness involved in choosing the groups makes it an even less attractive
option. However, our LP-divergence measure
$\operatorname{LP}(p_{0}\|{\widetilde{p}})$ (3.6) can be applied to ungrouped
data, with no trouble. The next section expands on this.

###  3.4 Explanatory Goodness-of-Fit

What is an explanatory goodness-of-fit? Why do we need it? This is perhaps
best answered by quoting the views of John Tukey:

> “What are we trying to do with goodness of fit tests? (Surely not to test
> whether the models fits exactly, since we know that no model fits exactly!)
> What then? How should we express the answer of such test?”
>
> —John Tukey (1954)

To satisfactorily answer these questions we need to design a GOF-procedure
that is simultaneously confirmatory and exploratory in nature: $\bullet$ On
the confirmatory side, it aims to develop a universal GOF statistic that is
easy to use and fully automated for any user-specified discrete $p_{0}(x)$.
One such universal GOF measure is $\operatorname{LP}(p_{0}\|{\widetilde{p}})$,
defined as

$\operatorname{LP}(p_{0}\|{\widetilde{p}})\,=\,\int_{0}^{1}\big{(}d(u;p_{0},{\widetilde{p}})-1\big{)}^{2}\;\mathrm{d}u=\sum_{j}\Big{|}\operatorname{LP}[j;p_{0},{\widetilde{p}}]\Big{|}^{2}=\sum_{j}\Big{|}\sum_{x}{\widetilde{p}}(x)T_{j}(x;F_{0})\Big{|}^{2},\leavevmode\nobreak\
\leavevmode\nobreak\ $ (3.8)

where the index $j$ runs over the significant components. See Appx. A.1 for
more details.

$\bullet$ On the exploratory side, the graphical visualization of comparison
density $d(u;p_{0},{\widetilde{p}})$ provides explanations as to why
$p_{0}(x)$ is inadequate for the data (if so) and how to rectify it to reduce
its incompatibility with the data. This has special significance for data-
driven discovery and hypothesis generation. In particular, the non-zero LP-
coefficients indicate the “main sources of discrepancies.” In the following,
we illustrate this method using three real data examples, each of which
contains a different degree of lack-of-fit.

######  Example 6 (Spiegel Family Data).

Spiegel (1972) reported a survey data of $n=320$ families with five children.
The numbers of families with $0,1,2,3,4$ and $5$ girls were $18,56,110,88,40$
and $8$. As an obvious model for $p_{0}$ we choose
Binomial$(5,\hat{\pi}=.463)$. Estimated LP-Fourier coefficients are

$n\times\operatorname{LP}(p_{0}\|{\widetilde{p}})\,=\,n\times\sum_{j=1}^{5}\left|\operatorname{LP}[j;p_{0},{\widetilde{p}}]\right|^{2}=1.489\leavevmode\nobreak\
\leavevmode\nobreak\ $ (3.9)

with $p$-value $0.92$ under the chi-square null with df $5$. In other words,
we have just shown that the comparison density is flat uniform $d_{0}(u)=1$,
hence the binomial distribution is completely acceptable for this data; no
further density sharpening is needed.

######  Remark 7.

Note that in our analysis, the prime object of interest is the shape of the
estimated $\widehat{d}_{0}(x)$ (because it addresses Tukey’s concern about the
practical utility of goodness-of-fit), not how big or small the $p$-value is.
But if a data analyst is habituated to using a threshold $p$-value as a basis
for decision making (not a good practice), then we recommend ‘double
parametric bootstrap’ (Beran, 1988) to compute the $p$-value— admittedly a
computationally demanding task. This adjusted $p$-value takes into account the
fact that the null-parameters are not given (e.g., here the binomial
proportion $\pi$); they are estimated from the data.

######  Example 7 (Rutherford-Geiger polonium data).

Rutherford, Geiger, and Bateman (1910) presented experimental data on groups
of alpha particles emitted by Polonium, a radioactive element, in 1/8 minute
intervals. On the whole, $n=2608$ time intervals were considered, in which $k$
($k=0,1,\ldots,14$) decays were observed. The following table summarizes the
data.

Table 2: Observed number of collisions of alpha particles emitted from polonium. 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 13 | 14
---|---|---|---|---|---|---|---|---|---|---|---|---|---
57 | 203 | 383 | 525 | 532 | 408 | 273 | 139 | 45 | 27 | 10 | 4 | 1 | 1

In the 1910 article, Bateman showed that (appended as a note at the end of the
original 1910 paper) the theoretical distribution of alpha particles observed
in a small interval follows Poisson law, which we select as our model-0. The
estimated $\operatorname{DS}(p_{0},m)$ model is given by

$\widehat{p}(x)\,=\,e^{-\lambda_{0}}\dfrac{\lambda_{0}^{x}}{x!}\,\Big{[}1-0.03T_{2}(x;F_{0})-0.04T_{3}(x;F_{0})\Big{]},\leavevmode\nobreak\
\leavevmode\nobreak\ \text{with}\leavevmode\nobreak\ \lambda_{0}=3.88$ (3.10)

which is displayed in Fig. 18 in the Appendix. The model (3.10) indicates that
there is a ‘gap’ between the theoretically predicted Poisson model and the
experimental result— second-order (under‐dispersed) and third-order (less
skewed) corrections are needed. To quantify the lack-of-fit, compute:

$n\times\operatorname{LP}(p_{0}\|{\widetilde{p}})\,=\,n\times\sum_{j\in\\{2,3\\}}\left|\operatorname{LP}[j;p_{0},{\widetilde{p}}]\right|^{2}=6.82,\leavevmode\nobreak\
\leavevmode\nobreak\ \text{with pvalue $0.033$}.$

This is a borderline case, where it is important to consult subject matter
specialists before choosing sides—scientific significance is as important as
statistical significance. Hoaglin (1980) came to a similar conclusion using an
exploratory diagnostic tool called “Poissonness plot,” shown in the Appendix
Fig. 19.

######  Example 8 (Sparrow data).

This data composed of numbers of sparrow nests found in plots of area one
hectare, the sample average being $\bar{x}=1.10$. Zar (1974) previously
analyzed this dataset. We choose Poisson($1.10$) as our $p_{0}(x)$ for
$\operatorname{DS}(p_{0},m)$ model. The second-order
$\operatorname{LP}[2;p_{0},{\widetilde{p}}]=-0.328$ (pvalue=.03) turns out to
be the only significant component, which indicates that the data exhibit less-
dispersion (due to the negative sign) than postulated ${\rm Poisson}(1.10)$.
Our finding is in agreement with Gürtler and Henze (2000). Finally, return the
$d$-modified under-dispersed Poisson model for the data:

$\widehat{p}(x)\,=\,e^{-\lambda_{0}}\dfrac{\lambda_{0}^{x}}{x!}\,\big{[}1-0.33T_{2}(x;F_{0})\big{]},\leavevmode\nobreak\
\leavevmode\nobreak\ \text{with}\leavevmode\nobreak\ \lambda_{0}=1.10$

which is displayed in Fig. 20 of the Appendix.

###  3.5 Relative Entropy and Model Uncertainty

How to quantify the uncertainty of the chosen model $p_{0}(x)$? A general
information-theoretic formula of model uncertainty is derived based on
relative entropy between the true (unknown) $p(x)$ and the hypothesized
$p_{0}(x)$. Express relative entropy (also called Kullback–Leibler divergence)
as a functional of maxent comparison density $d_{\bm{\theta}}$:

$\operatorname{KL}(p\|p_{0})\,=\,\sum_{x}p(x)\log\Big{\\{}\dfrac{p(x)}{p_{0}(x)}\Big{\\}}\,=\,\mathbb{E}_{F}\Big{[}\log
d_{\bm{\theta}}\big{(}F_{0}(X);p_{0},p\big{)}\Big{]}$ (3.11)

Substituting $F_{0}(x)=u$, leads to the following important formula for
relative entropy in terms of LP-parameters:

$\displaystyle\operatorname{KL}(p\|p_{0})$ $\displaystyle=$
$\displaystyle\int_{0}^{1}d(u;F_{0},F)\log d(u;F_{0},F)\;\mathrm{d}u$ (3.12)
$\displaystyle=$
$\displaystyle\int_{0}^{1}d(u;F_{0},F)\Big{\\{}\sum_{j}\theta_{j}S_{j}(u;F_{0})\,-\,\Psi(\bm{\theta})\Big{\\}}$
$\displaystyle=$
$\displaystyle\sum_{j}\theta_{j}\operatorname{LP}_{j}\,-\,\Psi(\bm{\theta}).$

The second equality follows from eq. (2.8) and the last one from eq. (2.8).
Based on a random sample $X_{1},\ldots,X_{n}$, a nonparametric estimate of
relative entropy is obtained by replacing the unknown LP-parameters in (3.12)
with their sample estimates.

Statistical Inference. Relative entropy-based statistical inference procedure
is applied to few real datasets in the context of model validation and
uncertainty quantification. $\bullet$ Estimation and Standard error: For the
earthquake data, we like to quantify the uncertainty of $p_{0}={\rm
NB}(12,19)$. The estimated value of $\operatorname{KL}(p\|p_{0})$ is $0.070\pm
0.020$ (bootstrap standard error, based on $B=1000$), indicating serious lack-
of-fit of the starting NB model—which matches with our previous conclusion;
see Fig. 1. $\bullet$ Testing: For the Spiegel family data, the estimated
relative entropy we get is $0.0087$, quite small. Naturally, we perform
(parametric bootstrap-based) testing to check if
$H_{0}:\operatorname{KL}(p\|p_{0})=0$. Generate $n$ samples from $p_{0}(x)$;
compute $\widehat{\operatorname{KL}}(p\|p_{0})$; repeat, say, $1000$ times;
return the $p$-value based on the bootstrap null-distribution. For this
example, the $p$-value we get is $0.093$, which reaffirms that binomial
distribution explains the data fairly well.

###  3.6 Card Shuffling Problem

Consider the following question (Aldous and Diaconis, 1986): How many times
must a deck of cards be shuffled until it is close to random? To check whether
a deck of $52$ cards is uniformly shuffled we use fixed point statistic, which
is defined as the number of cards in the same position after a random
permutation. Large values of fixed points (‘too many cards left untouched’) is
an indication that the deck is not well-mixed.

Theoretically-expected distribution: One of the classical theorems in this
direction is due to Pierre de Montmort (1713) who showed that the distribution
of the number of fixed-points under $H_{0}$ (random permutation of
$\\{1,2,\ldots,52\\})$ is approximately $p_{0}={\rm Poisson}(1)$.

Data and notation: Let $n$ denotes sample size and $k$ number of shuffles.
Then CARD$(k,n)$ stands for a dataset $X_{1},X_{2},\ldots,X_{n}$, where
$X_{i}$ is the number of fixed points of a $k$-shuffled deck. By
$\widetilde{p}_{k}$, we mean the sample distribution of fixed-point statistic
after $k$ random permutations. The goal is to find the minimum value of $k$,
such that it is safe to accept $H_{0}:p_{k}=p_{0}$, where, we should recall,
the null-distribution $p_{0}$ is ${\rm Poisson}(1)$.

Modeling via goodness-of-fit: Fig. 5 shows a CARD$(k,n)$ dataset with $k=150$
and $n=500$. There is a clear discrepancy between the observed probabilities
${\widetilde{p}}_{k}$ and the theoretical $p_{0}$. The estimated $d$-sharpen
Poisson$(1)$ is given below:

$\widehat{p}_{k}(x)=\dfrac{e^{-1}}{x\,!}\big{[}1+0.130\,T_{1}(x;F_{0})\big{]},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ x=0,1,\ldots$
(3.13)

which shows that ‘first order perturbation’ (location correction) is needed:
$\operatorname{LP}[1;p_{0},{\widetilde{p}}_{k}]=0.130$ with pvalue $0.003$.
The positive sign of $\operatorname{LP}_{1}$ indicates that the mean of the
fixed points distribution with $k=150$ is larger than the postulated
$\lambda_{0}=1$; more shuffling is needed to make the deck close to random.
The shape of (3.13) is shown in Fig. 5.

New updated mean. A curious reader might want to know precisely how large the
mean of $p_{k}$ is compared to $1$. For a distribution $F\sim{\rm
DS}(p_{0},m)$, we can write an expression for the mean of $F$ ($\lambda_{F}$)
in terms of the mean of $F_{0}$($\lambda_{0}$). In this case, we have

$\displaystyle\vspace{-1em}\lambda_{k}\leavevmode\nobreak\
=\leavevmode\nobreak\ \mathbb{E}_{F_{k}}[X]$ $\displaystyle=$
$\displaystyle\int
x\,\big{[}1+0.130\,T_{1}(x;F_{0})\big{]}\;\mathrm{d}F_{0}(x)$ (3.14)
$\displaystyle=$
$\displaystyle\int_{0}^{1}Q_{0}(u)\big{[}1+0.130\,S_{1}(u;F_{0})\big{]}\;\mathrm{d}u$
$\displaystyle=$
$\displaystyle\int_{0}^{1}Q_{0}(u)\;\mathrm{d}u+0.130\,\langle
Q_{0},S_{1}\rangle_{\mathcal{L}^{2}[0,1]}$ $\displaystyle=$ $\displaystyle
1+0.130\times 0.9596\approx 1.125.$

A few additional comments:

Figure 5: Distribution of fixed points: black is empirical and blue is
theoretical—a clear mismatch between data and the theory. The rectified
$\operatorname{DS}(p_{0},m)$ model is displayed in red.

$\bullet$ Thus far, we have verified that $k=150$ is not enough to produce a
uniformly shuffled deck. However, we came to this conclusion based on a single
sample of size $n$. So, we generate (through computer simulation) several
replications of CARD$(k,n)$ data with different $(k,n)$.

$\bullet$ To reach a confident decision, we perform the experiment with
$n=500$ and $k=150,160,170,180,190,200$. The analysis was done based on
$B=250$ datasets from each $(n,k)$-pair. The results are summarized in
appendix Fig. 21, which shows that $k=170$ shuffles is probably a safe bet to
declare a deck to be fair—i.e., uniformly distributed.

$\bullet$ Diaconis and Wang (2018) describes a Bayesian approach to this
problem. In contrast, we offered a completely nonparametric solution that
leverages the additional knowledge of the expected $p_{0}(x)$ and provides
more insights into the nature of discrepancies.

###  3.7 Jaynes Dice Problem

The celebrated Jaynes’ dice problem (Jaynes, 1962) is as follows: Suppose a
die has been tossed $N$ times (unknown) and we are told only that the average
number of the faces was $4.5$—not $3.5$, as we might expect from a fair die.
Given this information (and nothing else), the goal is to determine
probability assignment, i.e., what is the probability that the next throw will
result in face $k$, for $k=1,\ldots,6$.

Solution of Jaynes’ dice problem using density sharpening principle. The
initial $p_{0}(x)$ is selected as the discrete uniform distribution
$p_{0}(x)=1/6$ for $x=1,\ldots,6$, which reflects the null hypothesis of
‘fair’ dice. As we are given only the first-order location information (mean
is $4.5$) we consider the following $\operatorname{DS}(p_{0},m=1)$ model:

$p(x)=\dfrac{1}{6}\Big{\\{}1+\operatorname{LP}[1;p_{0},p]\,T_{1}(x;p_{0})\Big{\\}}$
(3.15)

The coefficient $\operatorname{LP}[1;p_{0},{\widetilde{p}}]$ has to be
estimated, and for that we also need to know the basis function
$T_{1}(x;F_{0})$.

Step 1. To find an explicit formula for the discrete basis $T_{1}(x;F_{0})$,
apply (2.3) with $F_{0}(x)=x/6$ and $p_{0}(x)=1/6$

$T_{1}(x;F_{0})\,=\,\sqrt{12}\dfrac{(x/6-.5)}{\sqrt{1-\sum_{x=1}^{6}(1/6^{3})}}\,=\,\sqrt{\dfrac{12}{35}}\big{(}x-3.5\big{)},\leavevmode\nobreak\
\leavevmode\nobreak\ \text{for}\,x=1,\ldots,6.$ (3.16)

Step 2. Compute $\operatorname{LP}[1;p_{0},{\widetilde{p}}]$ by applying
formula (2.9)

$\operatorname{LP}[1;p_{0},{\widetilde{p}}]\,=\,\sum_{x}{\widetilde{p}}(x)T_{1}(x;p_{0})\,=\,\sqrt{\dfrac{12}{35}}\big{(}\sum_{x}x{\widetilde{p}}(x)-3.5\big{)}\,=\,0.586.$
(3.17)

The non-zero $\operatorname{LP}[1;p_{0},{\widetilde{p}}]$ indicates it was a
loaded die.

Step 3. Substitute the value of $\operatorname{LP}[1;p_{0},{\widetilde{p}}]$
in eq. (3.15) to get the LP-Fourier $DS(p_{0},m=1)$ model as

$\widehat{p}(x)=\dfrac{1}{6}\Big{\\{}1+0.586\,T_{1}(x;p_{0})\Big{\\}},\leavevmode\nobreak\
\leavevmode\nobreak\ x=1,\ldots,6.$ (3.18)

This is shown as the blue curve in Fig. 6.

Step 4. Finally, return the estimated exponential
$\operatorname{DS}(p_{0},m=1)$ probability estimates

$\widehat{\widehat{p\,}}(x)(x)\,=\,\dfrac{1}{6}\exp\big{\\{}-0.193+0.634\,T_{1}(x;F_{0})\big{\\}},\leavevmode\nobreak\
\leavevmode\nobreak\ x=1,\ldots,6.$ (3.19)

This is shown as the red curve in Fig. 6.

Figure 6: The black dots denote the starting null-model—die with $6$ equally likely outcomes. The blue curve is our LP-Fourier density estimate given in eq. (3.18). The exponential $\operatorname{DS}(p_{0},m)$ density estimate, given in eq. (3.19), is shown in red line. Table 3: The estimated LP skew distribution for the Jaynes’ dice problem. $p_{0}(x)$ | 1/6 | 1/6 | 1/6 | 1/6 | 1/6 | 1/6
---|---|---|---|---|---|---
LP-Fourier $\widehat{p}(x)$ | .025 | .080 | .140 | .195 | .250 | .310
LP-MaxEnt $\widehat{\widehat{p\,}}(x)$ | 0.054 | 0.079 | 0.114 | 0.165 | 0.240 | 0.347
Jaynes’Answer | 0.054 | 0.079 | 0.114 | 0.165 | 0.240 | 0.347

######  Remark 8.

It is a quite remarkable fact that our density-sharpening principle-based
probability assignment exactly matches with Jaynes’ maxent answer; see Table
3.

###  3.8 Compressive Learning of Big Distributions

An important class of learning problem that has recently attracted researchers
from various disciplines—including high-energy physics, neuroscience,
theoretical computer science, and machine learning—can be viewed as a modeling
problem based on samples from a distribution over large ordered domains. Let
${\bf p}=(p_{1},\ldots,p_{k})$ be a probability distribution over a very large
domain $k$, where $p_{i}\geq 0,\sum_{i=1}^{k}p_{i}=1$. Let us look at a
realistic example before discussing a general method.

######  Example 9 (HEP data).

This is an example from high-energy physics222High-energy physics is not the
only discipline where this kind of very large and sparse histogram-like data
appears. It frequently arises in many modern scientific domains: inter-spike
interval data (neuronal firing patterns); relative abundance/intensity data
(mass spectrometry data); DNA methylation and ChIP-seq data (genomics); pixel
histogram data (astronomical images of stars, galaxies etc); histograms of
activity intensities (biosignals from wearable sensor devices, mental
illnesses studies by NIMH); photometric redshift data (photo-z spectra in
Cosmology), just to name a few. There is an outstanding interest in developing
new computational tool that allows rapid and approximate statistical learning
for big-histogram-like datasets. (HEP), motivated by the PHYSTAT 2011 Banff
bump-hunting challenge task (Junk, 2011). In HEP counting experiments (e.g.,
in Large Hadron Collider) one observes data in the following form: $n$ samples
from unknown $p(x)$ as event counts (number of collisions) at $k$ finely
binned energy-cells, which we denote by HEP$(k,n)$. Fig. 7 displays one such
data with $n=10,000$ and $k=250$, with the postulated background model
(dictated by the known Standard Model) as discretized exponential distribution
$f_{0}(x)=\lambda e^{-\lambda x}$ with $\lambda=1/20$:

$p_{0}(i)\,\doteq\,\int_{x_{i}}^{x_{i+1}}f_{0}(x)\;\mathrm{d}x,\leavevmode\nobreak\
\leavevmode\nobreak\ i=1,\ldots,k.$ (3.20)

Particle physicists are interested in discovering new particles that go beyond
the known Standard model described by the background model $p_{0}$. We present
a four-step algorithmic program to address the general problem of data-driven
‘Learning and Discovery.’

Phase 1. Testing. The first question a scientist would like answered is
whether the data is consistent with the background-only hypothesis, i.e.,
$H_{0}:p=p_{0}$. We perform the information-theoretic test described in Sec.
3.5. In particular, we choose the relative entropy-based formula given in
(3.12) as our test statistic. The $p$-value obtained using parametric
bootstrap (with $B=50,000$) is almost zero—strongly suggesting that the data
contain some surprising new information in light of the known physics model
$p_{0}$. But to figure out whether that information is actually useful for
physicists, we have to dig deeper.

Figure 7: HEP$(k,n)$ data analysis. Left: The estimated comparison density
$\widehat{d}(u;p_{0},p)$, which indicates there may be a bump (around
$u=0.64$) that went unnoticed by the theoretical model. Right: The data
generated from a mixture of $p_{0}$ and $\mathcal{N}(125,2)$ with mixing
proportion $0.1$. The theory-driven model $p_{0}$ is the blue curve and the
red curve is the $d$-sharpen $p_{0}$. The green triangle denotes the true peak
at the mass position 125 Gev. The shaded area under ${\widehat{d}}_{0}(u)$
denotes the amount of excess mass on the top of the smooth background.

Phase 2. Exploration and Discovery. By definition, new discoveries can be made
only by “contrasting” data with the existing model. This is what is achieved
through $d(u;F_{0},\widetilde{F})$. The left panel of Fig. 7 displays the
estimated $\widehat{d}_{0}(u)$ for the HEP-data, which compactly encodes all
the structure of the data that cannot be described by the assumed $p_{0}(x)$.

The exploratory graphical display of $\widehat{d}_{0}$ reveals a few
noteworthy points. Firstly, the non-uniformity of $\widehat{d}_{0}$ makes us
skeptical about $p_{0}$—this is completely in agreement with Phase-1 analysis.
Secondly and more importantly, the shape of $\widehat{d}_{0}$ provides a
refined understanding of the nature of new physics that is hidden in the data,
which, in this case, revealed itself as a bump above the smooth background
$p_{0}$. The word ‘above’ is important because we are not interested in bumps
on $p$ itself, but on $d_{0}$, which is the unanticipated ‘excess mass.’ Hunt
for new physics is the problem of bump hunting on $d(u;p_{0},p)$, not on
$p(x)$. For the HEP-data, we see a prominent bump in $d_{0}(u)$ around
$u=0.64$, which (in the original data domain) corresponds to near
$Q_{0}(.64)\approx$ 125 GeV; the green triangle in the above figure.

######  Remark 9 (The discovery function).

Since, $d_{0}(x)$ encapsulates what’s new in the data by separating the
unknown from the known, we also call it the “discovery function.” It is the
“missing piece” that glues together the known $p_{0}(x)$ and the unknown
$p(x)$. It provides a graphical diagnostic tool that exposes the unexpected.
These clues can guide domain-scientists to carry out more targeted follow-up
studies.

Phase 3. Inference and Excess Mass Problem. Where is the interesting excess
mass hiding? Is it a statistical fluke or something real? How substantial is
the evidence? The real issue is: can we let the data confidently tell us where
to look next for new particles? This will result in a complete paradigm shift
because traditionally the HEP searches (for locating excess events) were
guided by theoretical considerations only.

> ‘One may feel uneasy that we may therefore only find new processes if a
> theorist has been clever enough to propose the corresponding theory ahead of
> time.’
>
> —-Glen Cowan (2007)

Interested readers may also refer to Lyons (2008) and the Nature news article
by Castelvecchi (2018) for a clear exposition on the scientific importance of
these issues.

Statistical Discovery: Inference Algorithm. The following are the main steps
of the inference algorithm whose results are summarized in Fig. 8:

Step 1. Parametric bootstrap: To measure the natural statistical variation of
${\widehat{d}}_{0}(x)$ under the null hypothesis: simulate $n$ samples from
$p_{0}(x)$ and estimate the comparison density. Repeat the whole process for a
large number of iterations (say $B=10,000$ times) to get a bundle of
comparison density curves, all of which fluctuate around the flat uniform
line.

Step 2. Pointwise $p$-value computation: At a fixed grid point
$x\in[100,250]$, we have the following values of the test statistic

$\Big{\\{}\widehat{d}_{0}^{(1)}(x),\ldots,\widehat{d}_{0}^{(B)}(x)\Big{\\}}$

calculated from the $B$ bootstrap samples. Compute the bootstrap $p$-value at
the point $x$ by

${\rm Pval}(x)\,=\,\dfrac{1+\big{\\{}\\#\,{\rm
of}\,\widehat{d}_{0}^{(j)}(x)\geq\widehat{d}_{0}(x)\big{\\}}}{B+1}$

Fig. 8 draws the curve $-\log_{10}({\rm Pval}(x))$ as a function of $x$. The
$5\sigma$ discovery region $(121.5,129.5)$ is highlighted in yellow, which
includes the true excess mass point $125$ GeV. This is how modern
nonparametric modeling based on ‘density sharpening principle’ can
convincingly guide researchers on where to look for evidence of a deeper
theory of physics.

Figure 8:  The $5\sigma$ discovery interval $(121.5,129.5)$ correctly captures
the true excess-mass point indicated by the green triangle.

Phase 4. Sharpen Scientific-Model. Finally, the goal is to sharpen the initial
scientific model $p_{0}(x)$ to achieve a more precise description of what is
loosely known or suspected. The estimated $\operatorname{DS}(p_{0},m)$ model
sharpens the parametric null (3.20) to provide a nonparametrically-adjusted,
parsimonious model:

$\widehat{p}(x)\,=\,p_{0}(x)\Big{[}1-\sum_{j\mathcal{J}_{5}}\operatorname{LP}[j;p_{o},{\widetilde{p}}]\,T_{j}(x;F_{0})\Big{]},$
(3.21)

where the active set $\mathcal{J}_{5}=\\{2,3,5,7,8\\}$ along with the LP-
coefficients are given in Table 4.

Table 4: The selected no-zero LP-coefficients $\mathcal{J}_{5}$ | 2 | 3 | 5 | 7 | 8
---|---|---|---|---|---
$\widehat{\operatorname{LP}}_{j}$ | -.097 | -.090 | 0.117 | -.095 | -.060

######  Remark 10.

A few remarks:

1\. The part in the square brackets of (3.21) shows how to change the prior
scientific model $p_{0}(x)$ to make it consistent with the data. Knowing the
nature of the deficiency of the assumed model, is an important step towards
data-driven discovery. As George Box (2001) said: “discovery usually means
learning how to change the model.”

2\. LP-parameterization requires only $5$-dimensional sufficient statistics to
approximately capture the shape of the distribution! The “compressiveness” of
the LP-transformation—the ability to extract a low-dimensional
representation—makes it less data-hungry, as demonstrated in the next section.

3\. Our model (3.21) is a ‘hybrid’ between theory-driven and data-driven
model, which decouples the overall density into two components: expected
$p_{0}(x)$ and unexpected $d_{0}(x)$.

######  Remark 11 (An Appeal to Physicists: Hypothesis Testing $\neq$
Discovery Science).

Classical statistical inference puts too much emphasis on testing, $p$-value,
standard error, and confidence intervals, etc. This ideology is reflected in
the practice of high-energy physicists—which entirely revolves around antique
tools of hypothesis testing, likelihood ratio, and $p$-value. It’s time to
break the shackles of outdated data analysis technology that starts with
hypothesis testing and ends with a $p$-value. George Box (2001) expressed a
similar sentiment, arguing that the reason why engineering and the physical
sciences rarely use statistics is: “Much of what we have been doing is
adequate for testing but not adequate for discovery.”

In this section my purpose has been to introduce some modern statistical tools
and concepts that can help scientists with their everyday tasks of discovery
and deeper exploration of data. After all, one of the main goals of data
analysis is to sharpen the scientists’ mental model by revealing the
unexpected—a continuous cycle of knowledge refinement:

Theory Measurement Discovery $d_{0}(x)$ Better theory
Figure 9: Continuous learning by iterative model refinement: It develops
increasingly ‘better’ theory that explains new phenomena by broadening the
scope of the previous theory.

###  3.9 Data-efficient Statistical Learning

We are interested in the following statistical learning problem: Given a small
handful of samples $n\ll k$ from a big probability distribution of size $k$,
how to perform time-and-storage-efficient learning? When designing such
algorithms we have to keep in mind that they must be (i) data-efficient: can
learn from limited sample data $n\ll k$; and (ii) statistically powerful: can
detect “small” deviations.

Classical nonparametric methods characterize big probability distributions
using high-dimensional sufficient statistics based on histogram counts
$\\{N_{j}\\}_{j=1}^{k}$, which, obviously, requires a very large sample for
efficient learning. And as the required sample sizes increases, this slows
down the algorithm running-time which scales with $n$. Hence, most of the
‘legacy’ statistical algorithms (e.g., Pearson’s chi-square, Freeman-Tukey
statistic, etc.) become unusable on large-scale discrete data as ‘the minimum
number of data points required to obtain an acceptable answer is too large to
be practical.’ Indeed, detecting sparse structure in a data-efficient manner
from big distributions, such as from HEP($k,n$), is a challenging problem and
requires a “new breed” of computational algorithms. There has been some
impressive progress on this front by the Theoretical Computer Science
community; see Appendix A.3 for a related discussion on sub-linear algorithms
for big data.

Figure 10: HEP$(n,k=500)$ data for varying sample size $n$. Event counts
(number of collisions) at $k=500$ finely binned energy-cells. LPgof (shown in
red) requires 50% less data to reach the correct conclusion with power 1.

HEP Data Example. We generate HEP($k=500,n$) data for varying sample size $n$.
We used $350$ null simulated data sets to estimate the 95% rejection cutoffs
for all methods at the significance level $0.05$, and used $350$ simulated
data sets from alternative to approximate the power, as displayed in Fig. 10.
We compared our LPgof method (3.8) with two state-of-the-art algorithms
(proposed by theoretical computer scientists): (i) Valiant and Valiant (2017)
and (ii) Acharya et al. (2015)—interestingly, this exact test has been
proposed earlier by Zelterman (1987), which in Statistics literature is known
as Zelterman’s D-statistic333‘Those who ignore Statistics are condemned to
reinvent it’—Brad Efron.. Conclusion: LPgof requires 50% less data to reach
the correct conclusion with power 1.

Figure 11: The red power curve is LPgof method with $m=8$. First row: Null
distribution is $p_{0}=(1/4,1/4,1/2(k-2)\ldots,1/2(k-2))$. Alternatives are
$p_{\delta}=(1/4+\delta,1/4-\delta,1/2(k-2)\ldots,1/2(k-2))$, with
$\delta=1/8,1/10$. Second row: Null distribution is $p_{0}=U_{[k]}$.
Alternatives are $.95U_{[k]}+.05{\rm Zipf}(\alpha)$, with $\alpha=1.25,1.15$.
Third row: Null distribution is $p_{0}=U_{[k]}$. Alternative probabilities are
computed using the increments $\Delta D(j/k;\mu,\pi)$, $j=1,\ldots,k$ where
$D(u)=F(\Phi^{-1}(u);X)$, and $F(x;X)=(1-\pi)\Phi(x)+\pi\Phi(x-\mu)$ for
$\mu=-1.5$ and $\pi=\\{.2,.1\\}$. In all cases the dimension $k=5000$; and
sample sizes are sublinear in dimension, i.e., $n\asymp\sqrt{k}$

Empirical Power Comparisons. We compare the power of different methods under
six different settings, as described in the Figs 10 description, and will not
repeat this here. The overall conclusion is pretty clear: LPgof emerged as the
most powerful data-efficient test—it can detect new discoveries quickly and
more reliably. The prime reason for achieving this level of performance is
fully attributable to the good sparsity (energy compaction) property of “LP-
domain data analysis.” The specially-designed discrete LP-transformation basis
provides an efficient coordinate system that requires far fewer parameters
than the size of the distribution to capture the essential information. For
additional examples see Figs. 22 and 23 of the Appendix.

###  3.10 Discovery-source Separation Problem

“At this scale it is not possible to keep all the data (the LHC produces up to
a Petabyte of data per second) and it is essential to have efficient data-
filtering mechanisms so that we can separate the wheat from the chaff.”

— Bob Jones, Project Leader at CERN.

Imagine a typical scenario, where massive experimental data are distributed
across thousands of different computation units (say, servers). Each data
source has its own ‘personal’ distribution over a large domain of size
$k=1000$, from which, we only have access to $n=1000$ samples. Thus, the
observed data can be summarized as a collection of highly noisy and sparse
empirical distributions $\\{{\widetilde{p}}_{\ell}(x)\\}$. Suppose we are
given $900$ such sparsely-sampled empirical distributions from
$p_{0}=U_{[k]}$, $50$ from $.9U_{[k]}+.1{\rm Zipf}(1.25)$, and $50$ from a
mixture of $.9U_{[k]}+.1\Delta{\rm Beta}(50,50)$, where $\Delta{\rm
Beta}(50,50)$ denotes the increments of the cdf of ${\rm Beta}(50,50)$. Shapes
of these three source-distributions are shown in the left of Fig. 12.

Discovery-source Separation (DSS). The goal here is to design algorithms that
can quickly filter out the ‘interesting’ data sources, having potential for a
new physics discovery. A more ambitious goal would be to classify different
data sources based on their ‘nature’ of discoverability.

Algorithm. It consists of two main steps: Step 1. LP-Fourier Transform Matrix:
Define the LP-transform matrix $L\in\mathcal{R}^{g\times m}$ by

$L[\ell,j]=\operatorname{LP}[j;p_{0},{\widetilde{p}}_{\ell}]=\int
T_{j}(x;F_{0})\;\mathrm{d}\widetilde{F}_{\ell},$ (3.22)

for $\ell=1,\ldots,g=1000$ and $j=1,\ldots,m=10$.

Step 2. DSS-plot: Perform the singular value decomposition (SVD) of
$L=U\Lambda U^{T}$ $=\sum_{r}\lambda_{r}u_{r}u_{r}^{T}$, where $u_{ij}$ are
the elements of the singular vector matrix $U=(u_{1},\ldots,u_{m})$, and
$\Lambda={\rm diag}(\lambda_{1},\ldots,\lambda_{m})$, $\lambda_{1}\geq$
$\cdots\lambda_{m}\geq 0$. DSS-plot is the two-dimensional graph in the right
panel of Fig 12 (b), which is formed using the points
$(\lambda_{1}u_{1\ell},\lambda_{2}u_{2\ell})$, for $\ell=1,\ldots,g$ by taking
the top two dominant singular vectors. Different data sources are shown as
points, which captures the heterogeneity in terms of discoverability.

Figure 12: Discovery-source Separation (DSS) plot. Each data-batch is
HEP$(k=1000,n=1000)$ with different shapes, among them $900$ were generated
from $U_{[k]}$ (grey dotted line), $50$ were from $.9U_{[k]}+.1{\rm
Zipf}(1.15)$ (red dotted line), and $50$ from a mixture of $U_{[k]}$ and
increments of ${\rm Beta}(50,50)$ with mixing proportion $.1$ (blue dotted
line). In the right panel, we show the DSS-plot. The chunks of points around
zero indicate the null data sources. The sources are classified into three
groups based on their nature of discoverability.

Interpretation: DSS-plot displays and compares a large number of (empirical)
distributions ${\widetilde{p}}_{\ell}$ by embedding them in 2D Euclidean
space. The cluster of points (data sources) that are near to the origin are
the ones with background distribution. The distance from the origin

$\text{discovery-
index}_{\ell}=(\lambda_{1}u_{1\ell}-0)^{2}\,+\,(\lambda_{2}u_{2\ell}-0)^{2},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\ell=1,\ldots,g=1000\leavevmode\nobreak\ \leavevmode\nobreak\ $

can be interpreted as the “degree of newness” of that dataset. The DSS-plot
successfully separates different sources based on the statistical nature of
their signal components. Researchers (like Bob Jones) can use this tool to
_quickly identify_ interesting data sets for careful investigation.

##  4 Discussion

Compared to the rich and well-developed tools for continuous data modeling
problems, many companion discrete data modeling problems are still quite
challenging. This was the main motivation for undertaking this research, which
aimed at developing a widely-applicable general theory of discrete data
modeling. This paper makes three broad contributions to the field of
nonparametric statistical modeling:

1) Model-Sharpening Principle: We have introduced a new principle of
statistical model building, called ‘Density-Sharpening,’ which performs three
tasks in one step: model verification, model exploration, and model
rectification. This was the guiding principle behind our systematic theory of
discrete data modeling. As future work, we plan to explore how model-
sharpening principle can help developing ‘auto-adaptable’ machine-learning
models.

2) $\operatorname{DS}(p_{0},m)$ Model: We have introduced a new class of
nonparametric discrete probability model and robust estimation techniques that
has the capability to leverage researchers’ vague (misspecified) prior
knowledge. It comes with novel exploratory graphical methods for ‘discovering’
new knowledge from the data that investigators neither knew nor expected.

3) Unified Statistical Framework: Our modern nonparametric treatment of the
analysis of discrete data is shown to be rich enough to subsume a large class
of statistical learning methods as a special case. This inclusivity of the
general theory has some serious implications: Firstly, from a theoretical
angle, it deepens our understanding of the ties between different statistical
methods. Secondly, it simplifies practice by developing unified algorithms
with expanded capabilities. And finally, it is also expected to be beneficial
for modernizing Statistics curriculum in a way that is applicable to small- as
well as large-scale problems.

## Supplementary Appendix

The supplementary material includes additional theoretical and algorithmic
details.

## References

  * Acharya et al. (2015) Acharya, J., C. Daskalakis, and G. C. Kamath (2015).  Optimal testing for properties of distributions.  In Advances in Neural Information Processing Systems, pp. 3591–3599.
  * Aldous and Diaconis (1986) Aldous, D. and P. Diaconis (1986).  Shuffling cards and stopping times.  The American Mathematical Monthly 93(5), 333–348.
  * Beran (1988) Beran, R. (1988).  Prepivoting test statistics: a bootstrap view of asymptotic refinements.  Journal of the American Statistical Association 83(403), 687–697.
  * Bernreuter (1981) Bernreuter, D. (1981).  Seismic hazard analysis application of methodology, results, and sensitivity studies. volume 4.  Technical report, Lawrence Livermore National Lab., CA.
  * Box (2001) Box, G. (2001).  Statistics for discovery.  Journal of Applied Statistics 28(3), 285–299.
  * Castelvecchi (2018) Castelvecchi, D. (2018).  LHC physicists embrace brute-force approach to particle hunt.  Nature 560(7718), 293–295.
  * Cowan (2007) Cowan, G. (2007).  The small-n problem in high energy physics.  In Statistical Challenges in Modern Astronomy IV ASP Conference Series; Edited by G. Jogesh Babu and Eric D. Feigelson., Volume 371, pp. 75.
  * Csörgő (1983) Csörgő, M. (1983).  Quantile processes with statistical applications.  SIAM.
  * de Montmort (1713) de Montmort, P. R. (1713).  Essay d’analyse sur les jeux de hazard.  Jacques Quillau, Paris. Reprinted 1980 by Chelsea, New York.
  * Diaconis and Wang (2018) Diaconis, P. and G. Wang (2018).  Bayesian goodness of fit tests: a conversation for david mumford.  Annals of Mathematical Sciences and Applications (arXiv:1803.11251).
  * Gürtler and Henze (2000) Gürtler, N. and N. Henze (2000).  Recent and classical goodness-of-fit tests for the poisson distribution.  Journal of Statistical Planning and Inference 90(2), 207–225.
  * Hansen (2014) Hansen, L. P. (2014).  Nobel lecture: Uncertainty outside and inside economic models.  Journal of Political Economy 122(5), 945–987.
  * Hoaglin (1980) Hoaglin, D. C. (1980).  A poissonness plot.  The American Statistician 34(3), 146–149.
  * Jaynes (1962) Jaynes, E. T. (1962).  Information theory and statistical mechanics (notes by the lecturer).  Statistical Physics 3, Lectures from Brandeis Summer Institute 1962\. New York: WA Benjamin, Inc.
  * Junk (2011) Junk, T. R. (2011).  Banff Challenge 2.  Proceedings of the PHYSTAT 2011 Workshop and Fermi National Accelerator Laboratory Report.
  * Kagan (2010) Kagan, Y. Y. (2010).  Statistical distributions of earthquake numbers: consequence of branching process.  Geophysical Journal International 180(3), 1313–1328.
  * Kagan and Jackson (2000) Kagan, Y. Y. and D. D. Jackson (2000).  Probabilistic forecasting of earthquakes.  Geophysical Journal International 143(2), 438–453.
  * Lyons (2008) Lyons, L. (2008).  Open statistical issues in particle physics.  The Annals of Applied Statistics, 887–915.
  * Parzen (1998) Parzen, E. (1998).  Statistical methods mining, two sample data analysis, comparison distributions, and quantile limit theorems.  In B. Szyszkowicz (Ed.), Asymptotic Methods in Probability and Statistics, pp. 611–617. Amsterdam: North-Holland.
  * Rutherford et al. (1910) Rutherford, E., H. Geiger, and H. Bateman (1910).  The probability variations in the distribution of $\alpha$ particles.  The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 20(118), 698–707.
  * Saulo et al. (2020) Saulo, H., R. Vila, L. Paiva, and N. Balakrishnan (2020).  On a family of discrete log-symmetric distributions.  arXiv preprint arXiv:2005.09744.
  * Spiegel (1972) Spiegel, M. R. (1972).  Statistics.  McGraw-Hill.
  * Tukey (1954) Tukey, J. W. (1954).  Unsolved problems of experimental statistics.  Journal of the American Statistical Association 49(268), 706–731.
  * Valiant and Valiant (2017) Valiant, G. and P. Valiant (2017).  An automatic inequality prover and instance optimal identity testing.  SIAM Journal on Computing 46(1), 429–455.
  * Zar (1974) Zar, J. H. (1974).  Biostatistical analysis.  Prentice Hall, Englewood Cliffs, N.J.
  * Zelterman (1987) Zelterman, D. (1987).  Goodness-of-fit tests for large sparse multinomial distributions.  Journal of the American Statistical Association 82(398), 624–629.

##  5 Appendix

It consists of four sections.

###  A.1 Goodness-of-fit Methods: Automatic verses Manual

Traditional GOF methods, mainly inspired by Lancaster (1953), construct
orthonormal polynomials on a _case-by-case basis_ for each parametric discrete
distribution $p_{0}(x)$. This is generally done by solving the heavy-duty
Emerson recurrence (Emerson, 1968, Best and Rayner, 1999a, 2003a) relation,
which could be quite complicated for non-standard distributions. Few concrete
examples: The proposal of Best and Rayner (1997) for testing Binomial
distribution is based on Krawtchouk polynomial; Best and Rayner (2003b)
develop test for geometric distribution using Meixner orthonormal polynomials;
Best and Rayner (1999b) construct Poisson GOF method based on Poisson-Charlier
orthonormal polynomials; For a theory of classical (distribution-specific)
orthonormal polynomials refer to Chihara (2011).

The problem with this line of thinking is that it makes the whole procedure
hard to automate and generalize, which may turn off practitioners.

Instead of bookkeeping a long list of esoteric polynomials for each parametric
distribution, we want to devise a universal mechanism that works for any
arbitrary $p_{0}$—which will make it easy to use and compute. This is achieved
in Section 3.4 of our paper by employing LP-orthonormal polynomial bases. Its
universality (3.8) lies in the fact that, unlike traditional methods, we have
the ability to construct discrete orthonormal polynomials (Sec. 2.2) of the
user-specified $p_{0}(x)$ in a completely automatic and robust manner.

###  A.2 LPdiscrete: A Unified Algorithm

We summarize some of the key steps of our generic discrete data modeling
approach, which starts with data $\\{X_{1},\ldots,X_{n}\\}$ and an approximate
(misspecified) working model $p_{0}(x)$.

LPdiscrete: Algorithm Based on ‘Density-Sharpening’ Principle

1\. Input: Discrete sample distribution ${\widetilde{p}}(x)$ and the assumed
reference model $p_{0}(x)$. The corresponding cdfs are denoted by
$\widetilde{F}(x)$ and $F_{0}(x)$. 2\. Compute the first-order LP-basis
$T_{1}(x;F_{0})$ by standardizing mid-distribution transform

$T_{1}(x;F_{0})\,=\,\dfrac{\sqrt{12}\big{[}F_{0}^{{\rm
mid}}(x)-0.5\big{]}}{\sqrt{1-\sum_{x}p_{0}^{3}(x)}}.$

3\. Apply _weighted_ Gram-Schmidt procedure on
$\\{T_{1}^{2},T_{1}^{3},\ldots\\}$ to construct empirical orthonormal
polynomials $T_{j}(x;X)$ with respect to measure $F_{0}$

$\sum\nolimits_{x}p_{0}(x)T_{j}(x;F_{0})=0;\leavevmode\nobreak\
\leavevmode\nobreak\
\,\sum\nolimits_{x}p_{0}(x)T_{j}(x;F_{0})T_{k}(x;F_{0})=\delta_{jk},\leavevmode\nobreak\
\leavevmode\nobreak\ 1<j,k<M$

Define the Unit LP-bases $S_{j}(u;F_{0})=T_{j}(Q_{0}(u);F_{0})$, $0\leq u\leq
1$.

4\. Estimate the discrete LP-Fourier transform coefficients of $p(x)$ with
respect to $p_{0}(x)$

$\widetilde{\operatorname{LP}}_{j}\,:=\operatorname{LP}[j;p_{0};{\widetilde{p}}]\,=\,\mathbb{E}\big{[}T_{j}(X;F_{0});\widetilde{F}\big{]}\,=\,\sum_{x}{\widetilde{p}}(x)T_{j}(x;F_{0}),$

They are the coordinates of the true unknown distribution $F$ relative to
$F_{0}$

$\big{[}F\big{]}_{F_{0}}:=\Big{(}\operatorname{LP}[1;F_{0},F],\ldots,\operatorname{LP}[m;F_{0},F]\Big{)}.$

5\. Sparse LP-transform. Identify the ‘significant’ non-zero LP-coefficients
with $|\widetilde{\operatorname{LP}}_{j}|>2/\sqrt{n}$. A more sophisticated
method: identify indices $j$ using AIC (or BIC) model selection criterion
applied to $\widetilde{\operatorname{LP}}_{j}$ arranged in decreasing
magnitude

${\rm AIC}(m)\leavevmode\nobreak\ =\leavevmode\nobreak\ \textrm{Sum of squares
of first m coefficients}\,–\,2m/n.$

Choose $m$ to maximize $\operatorname{AIC}(m)$. Store the selected indicates
in the set $\mathcal{J}$. 6\. Compute the lack-of-fit of the assumed model
$p_{0}$

$\operatorname{LP}(p_{0}\|{\widetilde{p}})\,=\sum_{j\in\mathcal{J}}\Big{|}\operatorname{LP}[j;p_{0},{\widetilde{p}}]\Big{|}^{2}=\sum_{j\in\mathcal{J}}\Big{|}\sum_{x}{\widetilde{p}}(x)T_{j}(x;F_{0})\Big{|}^{2}.$

7\. Since $\operatorname{LP}(j;p_{0},{\widetilde{p}})$ is the inner product of
comparison density $\;\mathrm{d}\widetilde{F}/\;\mathrm{d}F_{0}$ with LP-basis
$T_{j}(x;F_{0})$, they are the coefficient of the orthogonal expansion of
$d_{0}(x)$ in LP-polynomial bases. Thus, the estimated LP-Fourier comparison
density is given by:

$\widehat{d}_{0}(u):=d\big{(}u;F_{0},\widetilde{F}\big{)}\,-\,1\,=\,\sum_{j\in\mathcal{J}}\widetilde{\operatorname{LP}}_{j}\,S_{j}(u;F_{0}).$

8\. Construct MaxEnt comparison density model based on selected LP-sufficient
statistics

$\widehat{d}_{\bm{\theta}}(u;F_{0},F)\,=\,\exp\Big{\\{}\sum_{j\in\mathcal{J}}\widehat{\theta}_{j}S_{j}(u;F_{0})\,-\,\Psi(\bm{\theta})\Big{\\}},\leavevmode\nobreak\
\leavevmode\nobreak\ 0<u<1$

9\. Estimate the relative entropy ${\rm KL}(p_{0}\|p)$ using the LP-parameter-
based formula:

$\sum_{j\in\mathcal{J}}\widehat{\theta}_{j}\widetilde{\operatorname{LP}}_{j}\,-\,\Psi(\bm{\theta})$

This quantifies the uncertainty that surrounds the initial tentative model
$p_{0}(x)$. 10\. Exploratory uncertainty analysis. The graphical visualization
of $\widehat{d}(u;F_{0},F)$ plays a central role in data-driven knowledge
discovery. It conveys the missing knowledge that was not anticipated by the
a-priori assumed $p_{0}(x)$—tool for detecting the unexpected. It also
provides insights as to how to change the misspecified $p_{0}$ to get a better
model that fits the data well. The purpose of this exploratory interface is to
encourages interactive data analysis and visual reasoning. 11\. Model
amendment via density-sharpening: Construct a class of reasonable models in
the neighborhood of the reference model $p_{0}(x)$

$\widehat{p}(x)\leavevmode\nobreak\ =\leavevmode\nobreak\
p_{0}(x)\times\widehat{d}_{0}(x),$

for $\widehat{d}_{0}(x)$ either use LP-Fourier canonical estimate or the LP-
maxent estimate. The bottom line: Density-sharpening provides a systematic
statistical model developmental process, describing how a relatively simple
model $p_{0}(x)$ can evolved into a more mature model $\widehat{p}(x)$ when it
comes in contact with new data. The $\widehat{p}(x)$ is not a conventional
nonparametric model, it is a ‘hybrid’ model that supplements domain-knowledge
with data-knowledge.

###  A.3 Sub-linear Algorithms for Detecting Sparse Structure

The problem of developing ‘efficient’ algorithms for detecting sparse
structure in big discrete distributions has recently attracted lots of
attention from the Theoretical Computer Science (TCS) community. By
‘efficient’ we mean algorithms that require less data to correctly detect the
presence of signals. We call an algorithm sub-linear if the sample complexity
($n$) is sub-linear in size of the domain ($k$).

Some interesting new ideas and results recently appeared in the Theoretical
Computer Science (TCS) literature under the banner of ‘Sub-linear Algorithms
for Big data,’ ‘Big data Property testing,’ etc; see, for example, Canonne
(2015), Indyk et al. (2012), Valiant and Valiant (2011), Rubinfeld (2012),
Diakonikolas (2016), Orlitsky et al. (2003), Orlitsky and Suresh (2015), Jiao
et al. (2015). They promise to give 20th-century pioneering statistical
modeling ideas a 21st-century makeover with new ideas from TCS. They deserve
the credit for asking the right question. However, all of these methods are
_too specialized_ and develops “new tricks” for each new related learning
problem, which makes them practically inconvenient. This paper presents a more
systematic and adaptable strategy, which starts with the following question
(the holy grail of data analysis): _How to find an efficient coordinate system
adapted to the given data sets sitting on a high-dimensional discrete domain
to deliver time and sample efficient computation?_ Section 3.9 of the main
paper, performed an extensive numerical comparison of the power of different
methods. Here we perform six further experiments. The detailed settings and
results are summarized in Fig. 23, which further re-establish the power of LP-
domain statistical learning for detecting sparse structure in big discrete
distributions.

###  A.4 Additional Figures

Some figures of the main article have been moved to the Appendix to save
space.

Figure 13: Earthquake Data: Plot of
$\sqrt{n}\big{|}\widetilde{\operatorname{LP}}_{j}\big{|}$, $j=1,2,\ldots,10$.
The only significant component is
$\sqrt{n}\widetilde{\operatorname{LP}}_{6}=2.10$, marked in red circle. see
Sec. 2.3 for more details. .

Figure 14: Sparse dice problem: Shape of the top four LP-orthonormal basis
functions $\\{S_{j}(u;F_{0})\\}_{1\leq j\leq 4}$ for $p_{0}(x)$ given in
Example 5. Figure 15: Sparse dice problem: Plot of
$\sqrt{n}\big{|}\widetilde{\operatorname{LP}}_{j}\big{|}$, $j=1,2,\ldots,10$.
The two significant LP-Fourier coefficients are marked with red circles.

Figure 16: Shows the analysis of $X$, which denotes a discrete random variable
= number of times computer broke down in each of 128 consecutive weeks of
operation. See Example 3 of Section 3.2 for more details. The flat comparison
density indicates that the negative binomial (with $\hat{\mu}=4$ and
$\hat{\phi}=1.70$) adequately captures the pattern in the data. No further
adjustments are needed. Figure 17: Spiegel Family Data: Plot of
$\sqrt{n}\big{|}\widetilde{\operatorname{LP}}_{j}\big{|}$, $j=1,2,\ldots,5$
where the $p_{0}(x)$ is taken to be ${\rm Binomial}(5,p=0.4625)$. None of the
LP-coefficients are significant, thus the Binomial-model is accepted. see
Example 6 for more details. . Figure 18: Rutherford-Geiger Polonium Data: The
estimated ${\rm DS}(p_{0},m)$ model is shown in red. Here $p_{0}(x)$ is ${\rm
Poisson}(\lambda=3.87)$ shown in blue. see Example 7 for more details.
Figure 19: Rutherford’s polonium data. Hoaglin’s (1980) diagnostic plot to
measure the “Poissonness” of the data. Note that when the data truly generated
from ${\rm Poisson}(\lambda)$ we would have
$\widetilde{\phi}(x):=\log(x!{\widetilde{p}}(x))$ will be close to a linear
function of $x$, given by $\lambda+(\log\lambda)x$. The empirical function
$\widetilde{\phi}(x)$ is called ‘metameter,’ denoted by open circles in the
above plot. The dashed lines denote the confidence intervals, computed
following Hoaglin and Tukey (1985, Ch. 9); centers of those intervals are
marked with dark circles. The length of the CIs are large for data-sparse
points $x\geq 11$; see Table 2. The Poisson model is slightly dubious for the
polonium data because the fitted line is outside the confidence interval at
$x=8$. However, this is not a make-or-break factor; According to Tukey and
Hoaglin “we need to recall that one miss out of 15 is very common at 5%.” This
perfectly matches with our assessments; see Example 7 of Sec. 3.4. Figure 20:
Sparrow data: $\operatorname{DS}(p_{0},m)$ finds a under-dispersed Poisson
model, which has a distinct peak at $x=1$. Figure 21: Card shuffling data
(Sec. 3.6): $P$-value for CARD($k,n$) data with $n=500$ as a function of
number of shuffles $k$. At each $(n,k)$ combination, we generated $B=250$ sets
of random sample of size $n=500$. The average of those $B$ $p$-values are
displayed, which shows $k=170$ random transpositions may be enough to produce
a well-mixed deck. Figure 22: To examine whether the type-I error is
controlled, we simulated data under the null $p_{0}=U_{[k]}$ for $k=5000$ with
sample size $n$ ranging from $15$ to $700$. All the methods maintain type-I
error rate quite well.

Figure 23: First row: The null distribution is $p_{0}=U_{[k]}$ and step-
function alternative $k^{-1}\alpha\,\mathbbm{1}_{1\leq j\leq
k/2}+k^{-1}(2-\alpha)\,\mathbbm{1}_{k/2<j\leq k}$ with $\alpha=.5,.7$. Second
row: Construct discrete distribution as increments $\Delta F(j/n)$,
$j=1,.\ldots,k$. For null model we choose $p_{0}$ $={\rm
Beta}(\alpha=2,\beta=2)$, and for alternatives $p_{\theta}$ is ${\rm
Beta}(\alpha=2,\beta=2+\theta)$ with $\theta=1,.5$. Third row: Null
distribution is $p_{0}=U_{[k]}$. The ‘bumpy’ alternative distribution is
constructed by taking the mixture of $U_{[k]}$ and increments of ${\rm
Beta}(10,10)$ with mixing proportion $.3$ and $.2$. The red power curve
denotes the LPgof method (see Sec 3.9) with $m=8$.

## References

  * Best and Rayner (1997) Best, D. and J. Rayner (1997).  Goodness of fit for the binomial distribution.  Australian Journal of Statistics 39(3), 355–364.
  * Best and Rayner (1999a) Best, D. and J. Rayner (1999a).  Goodness of fit for the poisson distribution.  Statistics & probability letters 44(3), 259–265.
  * Best and Rayner (1999b) Best, D. and J. Rayner (1999b).  Goodness of fit for the poisson distribution.  Statistics & probability letters 44(3), 259–265.
  * Best and Rayner (2003a) Best, D. and J. Rayner (2003a).  Tests of fit for the geometric distribution.  Communications in Statistics-Simulation and Computation 32(4), 1065–1078.
  * Best and Rayner (2003b) Best, D. and J. Rayner (2003b).  Tests of fit for the geometric distribution.  Communications in Statistics-Simulation and Computation 32(4), 1065–1078.
  * Canonne (2015) Canonne, C. L. (2015).  A survey on distribution testing: Your data is big. but is it blue?  Electronic Colloquium on Computational Complexity (ECCC) 22(63).
  * Chihara (2011) Chihara, T. S. (2011).  An introduction to orthogonal polynomials.  Courier Corporation.
  * Diakonikolas (2016) Diakonikolas, I. (2016).  Learning structured distributions.  Handbook of Big Data, 267.
  * Emerson (1968) Emerson, P. L. (1968).  Numerical construction of orthogonal polynomials from a general recurrence formula.  Biometrics, 695–701.
  * Hoaglin (1980) Hoaglin, D. C. (1980).  A poissonness plot.  The American Statistician 34(3), 146–149.
  * Hoaglin and Tukey (1985) Hoaglin, D. C. and J. W. Tukey (1985).  Checking the shape of discrete distributions.  Exploring Data Tables, Trends and Shapes, 345–416.
  * Indyk et al. (2012) Indyk, P., R. Levi, and R. Rubinfeld (2012).  Approximating and testing k-histogram distributions in sub-linear time.  Proceedings of the 31st symposium on Principles of Database Systems, 15–22.
  * Jiao et al. (2015) Jiao, J., K. Venkat, Y. Han, and T. Weissman (2015).  Minimax estimation of functionals of discrete distributions.  Information Theory, IEEE Transactions on 61(5), 2835–2885.
  * Lancaster (1953) Lancaster, H. (1953).  A reconciliation of $\chi^{2}$, considered from metrical and enumerative aspects.  Sankhyā: The Indian Journal of Statistics (1933-1960) 13(1/2), 1–10.
  * Orlitsky et al. (2003) Orlitsky, A., N. P. Santhanam, and J. Zhang (2003).  Always good turing: Asymptotically optimal probability estimation.  Science 302(5644), 427–431.
  * Orlitsky and Suresh (2015) Orlitsky, A. and A. T. Suresh (2015).  Competitive distribution estimation: Why is good-turing good.  Advances in Neural Information Processing Systems, 2134–2142.
  * Rubinfeld (2012) Rubinfeld, R. (2012).  Taming big probability distributions.  XRDS: Crossroads, The ACM Magazine for Students 19(1), 24–28.
  * Valiant and Valiant (2011) Valiant, G. and P. Valiant (2011).  Estimating the unseen: an n/log (n)-sample estimator for entropy and support size, shown optimal via new clts.  Proceedings of the forty-third annual ACM symposium on Theory of computing, 685–694.