Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
73,037 characters · 21 sections · 25 citation commands
\vskip.25em
\noindentKeywords: Density sharpening; $\DS(p_0,m)$ distributions; LP-Fourier analysis; Explanatory goodness-of-fit; Jaynes' dice problem; Compressive $\chi^2$; Data-efficient learning. \vskip1.7em
{.5ex} { \setcounter{tocdepth}{2} } {1.4ex} \setstretch{1.58}
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.
\vskip.34em 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.' \vskip.34em
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.
\vskip.34em
Organization. Section (ref) introduces a new family of nonparametric approximation and smoothing techniques for discrete probability distributions, which is built on the principle of `Density Sharpening.' Section (ref) 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 (ref). Additional applications and methodological details are deferred to the Supplementary Appendix to ensure the smooth flow of the main ideas.
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 (ref). Here is a bird's eye view of the core mechanism, which is a three-stage process. \vskip.35em 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 $\hp(x)$ that will eliminate the incompatibility with the data. \vskip.65em
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 $\Fmn(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) \ge u\}$ for $0<u<1$. By $\cL^2(dF_0)$ we mean the set of all square integrable functions with respect to the discrete measure $\dd F_0$, i.e, for a function $\psi \in \cL^2(dF_0)$: $\int |\psi|^2 \dd F_0 := \sum_x |\psi(x)|^2 p_0(x) < \infty$. The inner product of two functions $\psi_1$ and $\psi_2$ in $\cL^2(dF_0)$ will be denoted by $\langle \psi_1, \psi_2 \rangle_{F_0}:=\int \psi_1 \psi_2 \dd F_0$. Expectation with respect to $p_0(x)$ will be abbreviated as $\Ex_0(\psi(X)) :=\int \psi \dd F_0$.
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)$. \vskip.3em
The task is to nonparametrically approximate $d(F_0(x);F_0,F)$ to be able to apply the density sharpening equation (ref). We approximate $d \hspace{-.08em}\circ \hspace{-.08em}F_0(x) \in \cL^2({dF_0})$ by projecting it into a space of polynomials of $F_0(x)$ that are orthonormal with respect to the base measure$\dd 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.
We describe a general theory of constructing LP-polynomials---a new class of robust polynomials $\{T_j(x;F_0)\}_{j\ge 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)$.
Step 1: Define the first-order LP-basis function as standardized mid-distribution transform: \beq T_1(x;F_0) \,=\,\dfrac{\sqrt{12} \big[\Fmn(x) - 0.5\big]}{\sqrt{1-\sum_x p_0^3(x)}}. \eeq Verify that $\Ex_0[T_1(X;F_0)]=0$ and $\Ex_0[|T_1(X;F_0)|^2]=1$, since $\Ex[\Fmn(X)]=1/2$ and $\Var[\Fmn(X)]=\sqrt{(1-\sum_x p_0^3(x))}/12$. \vskip.25em
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$\dd F_0$ \[\sum\nolimits_x p_0(x) T_j(x;F_0)=0;~~\,\sum\nolimits_x p_0(x) T_j(x;F_0)T_k(x;F_0)=\delta_{jk}, ~~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 (ref). \vskip.3em Fig (ref) 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.
To estimate the unknown coefficients $\LP[j;F_0,F]$ of the $\DS(p_0,m)$ model, note the following important identity: \bea \LP[j;F_0,F]&=& \int d(F_0(x);F_0,F) T_j(x;F_0) \dd F_0(x) \nonumber \\ &=& \int T_j(x;F_0) \dd F(x) \nonumber \\ &=& \Ex_F\big[ T_j(X;F_0) \big]. \eea This immediately leads to the following “weighted mean” estimator: \beq \tLP_j\,:= \LP[j;F_0;\widetilde F]\,=\, \Ex_{\wtF}\big[T_j(X;F_0)\big]\,=\,\sum_x \tp(x) T_j(x;F_0), \eeq Using standard empirical process theory csorgHo1983quantile,parzen1998statistical one can show that the limiting distribution of sample LP-statistic is i.i.d $\cN(0,n^{-1/2})$, under the null hypothesis $H_0:p=p_0$. Thus one can quickly obtain a sparse estimated $\DS(p_0,m)$ model by retaining only the `significant' LP-coefficients, which are greater than $2/\sqrt{n}$.
{\bf Earthquake Data Example}. The first $m=10$ estimated $\tLP_j$ are shown in Fig. (ref) of the appendix, which indicates that the only interesting non-zero LP-coefficient is $\tLP_6$. The explicit form of the estimated $\DS({\rm NB},m=6)$ model for the earthquake data is given by: \beq \hp(x) = p_0(x)\big[ 1 + 0.20 T_6(x;F_0) \big],\eeq where $p_0={\rm NB}(x;\mu=19, \phi=12)$. The resulting $\hp(x)$ is plotted as a red curve in Fig. (ref).
To ensure non-negativity, we expand $\log d$ (instead of $d$ as we have done in Eq. (ref)) in LP-Fourier series, which results in the following exponential model: \beq d_{\teb}(u;F_0,F)\,=\,\exp\Big \{ \sum_{j\ge 1} \te_j S_j(u;F_0)\,-\, \Psi(\teb)\Big \}, 0<u<1 \eeq where $\Psi(\teb)=\log \int_0^1 \exp\{ \sum_j \te_j S_j(u;F_0)\}\dd u.$ This model is also called the maximum-entropy (maxent) comparison density model because it maximizes the entropy $-\int d_{\teb} \log d_{\teb}$ (flattest possible; thus promotes smoothness) under the following LP-moment constraints: \beq \Ex_{\teb}[S_j(U;F_0)]\,=\,\LP[j;F_0,\wtF], (j=1,2\ldots). \eeq LP-moments $\LP[j;F_0,\wtF]$ are `compressed measurements' (linear combinations of observed data; verify from (ref)), which are sufficient statistics for the comparison density $d_{\teb}$. \vskip.4em
\vskip.1em
{\bf Earthquake Data Example}. The estimated maxent $\DS(p_0,m)$ model for the earthquake distribution is given by \beq \hhp(x) = p_0(x)\exp\big \{ 0.195 T_6(x;F_0) - 0.02\big \},\eeq whose shape is almost indistinguishable from the LP-Fourier estimated p.m.f. (ref).
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.
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 $\DS(p_0,m=1)$ model \beq p(x)=p_0(x) \Big\{1+\LP[1;p_0,p] T_1(x;F_0)\Big\}, x=0,1.\eeq 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'_0:\LP[1;p_0,p]=0$.
Step 3. We derive an explicit formula for $\LP[1;p_0,\tp]$. It's a two step process: First, we need the analytic expression of the LP-basis $T_1(x;p_0)$ \beq T_1(x;p_0) = \left\{
\right. \eeq We then apply formula (ref) to deduce: \beq \LP(1;p_0;\tp) = (1-\tp)T_1(0;p_0) + \tp T_1(1;p_0) = \dfrac{\widetilde p - p_0}{\sqrt{p_0(1-p_0)}}. \eeq Step 4. A remarkable fact is that the test based on (ref) exactly matches with the classical Z-test, whose null distribution is: $\sqrt{n} \LP[1;p_0,\tp] \sim \cN(0,1)$ as $\nti$. 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.
We will focus now on one important special case of maxent $\DS(p_0,m)$ family of distributions (ref), where the base measure $p_0(x)$ is taken to be a negative Binomial (NB) distribution. \beq 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\ge 1} \te_j T_j(x;F_0)\,-\, \Psi(\teb)\Big \}, x \in \mathbb{N}, \eeq note that the basis functions $\{T_j(x;F_0)\}_{j\ge 1}$ are specially-designed LP-orthonormal polynomials associated with the base measure $p_0 = {\rm NB}(\phi,\mu)$. We call (ref) 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 $\texttt{XNB}$-distribution in (ref) in the context of modeling the earthquake distribution. We now turn our attention to two further real-data examples.
Given a random sample of size $n$ from the unknown population distribution $p(x)$, chi-square goodness-of-fit statistic between the sample probabilities $\tp(x)$ and the expected $p_0(x)$ can be re-written as follows: \beq \dfrac{\chi^2}{n} = \sum_x \dfrac{\big( \tp(x) - p_0(x) \big)^2}{p_0(x)} = \sum_x p_0(x) \big[ \tp(x)/p_o(x)\,-1 \big]^2= \int_0^1 \big[ d(u;p_0,\tp) - 1\big]^2\dd u,\eeq By applying Parseval's identity on the LP-Fourier expansion of $d$, we have the following important equality: \beq \dfrac{\chi^2}{n}\, =\, \sum_{j=1}^{r-1}\Big| \LP[j;p_0,\tp] \Big|^2\defeq\,\LP(p_0 \| \tp), \eeq 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.
What is an explanatory goodness-of-fit? Why do we need it? This is perhaps best answered by quoting the views of John Tukey:
To satisfactorily answer these questions we need to design a GOF-procedure that is simultaneously confirmatory and exploratory in nature: \vskip.3em $\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 $\LP(p_0 \| \tp)$, defined as \beq \LP(p_0 \| \tp)\,=\,\int_0^1 \big(d(u;p_0,\tp) - 1\big)^2 \dd u = \sum_j \Big| \LP[j;p_0,\tp] \Big|^2= \sum_j \Big|\sum_x \tp(x) T_j(x;F_0)\Big|^2, \eeq where the index $j$ runs over the significant components. See Appx. (ref) for more details.
\vskip.3em $\bullet$ On the exploratory side, the graphical visualization of comparison density $d(u;p_0,\tp)$ 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.
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_{\teb}$: \beq \KLD(p\|p_0)\,=\,\sum_x p(x) \log \Big\{\dfrac{p(x)}{p_0(x)}\Big\}\,=\,\Ex_F\Big[ \log d_{\teb}\big(F_0(X); p_0, p\big)\Big] \eeq Substituting $F_0(x)=u$, leads to the following important formula for relative entropy in terms of LP-parameters: \bea \KLD(p\|p_0)&=&\int_0^1 d(u;F_0,F) \log d(u;F_0,F) \dd u \nonumber \\ &=& \int_0^1 d(u;F_0,F) \Big\{ \sum_j \te_j S_j(u;F_0)\,-\, \Psi(\teb) \Big\} \nonumber \\ &=&\sum_j \te_j \LP_j\,-\,\Psi(\teb). \eea The second equality follows from eq. (ref) and the last one from eq. (ref). 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 (ref) with their sample estimates.
\vskip.34em
{\bf Statistical Inference}. Relative entropy-based statistical inference procedure is applied to few real datasets in the context of model validation and uncertainty quantification. \vskip.1em $\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 $\KLD(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. (ref). \vskip.2em $\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: \KLD(p\|p_0)=0$. Generate $n$ samples from $p_0(x)$; compute $\widehat{\KLD}(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.
Consider the following question aldous1986: 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.
\vskip.23em Theoretically-expected distribution: One of the classical theorems in this direction is due to Pierre de1713essay 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)$.
\vskip.23em
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)$.
\vskip.23em Modeling via goodness-of-fit: Fig. (ref) shows a CARD$(k,n)$ dataset with $k=150$ and $n=500$. There is a clear discrepancy between the observed probabilities $\tp_k$ and the theoretical $p_0$. The estimated $d$-sharpen Poisson$(1)$ is given below: \beq \widehat p_k(x)= \dfrac{e^{-1}}{x\,!} \big[ 1+ 0.130 \,T_1(x;F_0)\big], x=0,1,\ldots\eeq which shows that `first order perturbation' (location correction) is needed: $\LP[1;p_0,\tp_k]=0.130$ with pvalue $0.003$. The positive sign of $\LP_1$ indicates that the mean of the fixed points distribution with $k=150$ is larger than the postulated $\la_0=1$; more shuffling is needed to make the deck close to random. The shape of (ref) is shown in Fig. (ref).
\vskip.3em 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$ ($\la_F$) in terms of the mean of $F_0$($\la_0$). In this case, we have \bea \la_{k} = \Ex_{F_k}[X]&=& \int x \,\big[ 1+ 0.130 \,T_1(x;F_0) \big] \dd F_0(x) \nonumber\\ &=&\int_0^1 Q_0(u) \big[ 1+ 0.130 \,S_1(u;F_0) \big] \dd u \nonumber\\ &=&\int_0^1 Q_0(u) \dd u + 0.130 \,\langle Q_0, S_1 \rangle_{\cL^2[0,1]} \nonumber\\ &=& 1 + 0.130 \times 0.9596 \approx 1.125. \eea
A few additional comments:
$\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. (ref), which shows that $k=170$ shuffles is probably a safe bet to declare a deck to be fair---i.e., uniformly distributed.
$\bullet$ diaconis2018bayGOF 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.
The celebrated Jaynes' dice problem jaynes62 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$.
\vskip.25em 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 $\DS(p_0,m=1)$ model: \beq p(x) = \dfrac{1}{6}\Big\{ 1 + \LP[1;p_0,p]\,T_1(x;p_0) \Big\} \eeq The coefficient $\LP[1;p_0,\tp]$ 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 (ref) with $F_0(x)=x/6$ and $p_0(x)=1/6$ \beq 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), for\,x=1,\ldots,6. \eeq
Step 2. Compute $\LP[1;p_0,\tp]$ by applying formula (ref) \beq \LP[1;p_0,\tp] \,=\, \sum_x \tp(x) T_1(x;p_0) \,= \,\sqrt{\dfrac{12}{35}} \big( \sum_x x \tp(x) - 3.5\big) \,=\, 0.586. \eeq The non-zero $\LP[1;p_0,\tp]$ indicates it was a loaded die.
Step 3. Substitute the value of $\LP[1;p_0,\tp]$ in eq. (ref) to get the LP-Fourier $DS(p_0,m=1)$ model as \beq \widehat{p}(x) = \dfrac{1}{6}\Big\{ 1 + 0.586\,T_1(x;p_0) \Big\}, x=1,\ldots,6. \eeq This is shown as the blue curve in Fig. (ref).
Step 4. Finally, return the estimated exponential $\DS(p_0,m=1)$ probability estimates \beq \hhp(x) (x)\,=\,\dfrac{1}{6}\exp\big\{ -0.193 + 0.634\,T_1(x;F_0) \big\}, x=1,\ldots,6. \eeq This is shown as the red curve in Fig. (ref).
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 \ge 0, \sum_{i=1}^k p_i =1$. Let us look at a realistic example before discussing a general method.
\vskip.45em
\vskip.3em {\bf 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. (ref). In particular, we choose the relative entropy-based formula given in (ref) 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.
\vskip.15em
{\bf 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,\wtF)$. The left panel of Fig. (ref) displays the estimated $\whd_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)$.
\vskip.15em The exploratory graphical display of $\whd_0$ reveals a few noteworthy points. Firstly, the non-uniformity of $\whd_0$ makes us skeptical about $p_0$---this is completely in agreement with Phase-1 analysis. Secondly and more importantly, the shape of $\whd_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.
\vskip.3em {\bf 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.
Interested readers may also refer to lyons08 and the Nature news article by castelvecchi2018lhc for a clear exposition on the scientific importance of these issues.
{\bf Statistical Discovery: Inference Algorithm}. The following are the main steps of the inference algorithm whose results are summarized in Fig. (ref):
Step 1. Parametric bootstrap: To measure the natural statistical variation of $\dhat_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\{ \whd_0^{(1)}(x), \ldots, \whd_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}\, \whd_0^{(j)}(x) \ge \whd_0(x)\big\} }{B+1 }\] Fig. (ref) 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.
\vskip.3em {\bf 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 $\DS(p_0,m)$ model sharpens the parametric null (ref) to provide a nonparametrically-adjusted, parsimonious model: \beq \hp(x)\,=\,p_0(x)\Big[ 1 -\sum_{j \mathcal{J}_5} \LP[j;p_o,\tp]\,T_j(x;F_0)\Big], \eeq where the active set $\mathcal{J}_5=\{2,3,5,7,8\}$ along with the LP-coefficients are given in Table (ref).
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.
\vskip.3em 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 (ref) for a related discussion on sub-linear algorithms for big data.
\vskip.3em {\bf 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. (ref). We compared our LPgof method (ref) with two state-of-the-art algorithms (proposed by theoretical computer scientists): (i) valiant2017automatic and (ii) acharya2015optimal---interestingly, this exact test has been proposed earlier by zelterman1987, which in Statistics literature is known as Zelterman's D-statistic\footnote{`Those who ignore Statistics are condemned to reinvent it'---Brad Efron.}. Conclusion: LPgof requires {\bf 50% less} data to reach the correct conclusion with power 1.
{\bf Empirical Power Comparisons}. We compare the power of different methods under six different settings, as described in the Figs (ref) 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. (ref) and (ref) of the Appendix.
“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.”
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 $\{\tp_\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. (ref).
\vskip.35em
{\bf 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.
\vskip.35em {\bf Algorithm}. It consists of two main steps: \vskip.2em Step 1. LP-Fourier Transform Matrix: Define the LP-transform matrix $L \in \cR^{g\times m}$ by \beq L[\ell,j]=\LP[j;p_0,\tp_\ell] = \int T_j(x;F_0) \dd \wtF_\ell, \eeq for $\ell=1,\ldots,g=1000$ and $j=1,\ldots, m=10$. \vskip.25em
Step 2. DSS-plot: Perform the singular value decomposition (SVD) of $L=U\Lambda U^{T}$ $= \sum_r \la_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}(\la_1,\ldots,\la_m)$, $\la_1 \ge$ $ \cdots \la_m \ge 0$. DSS-plot is the two-dimensional graph in the right panel of Fig (ref) (b), which is formed using the points $(\la_1 u_{1\ell}, \la_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.
\vskip.4em Interpretation: DSS-plot displays and compares a large number of (empirical) distributions $\tp_\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=(\la_1 u_{1\ell} -0)^2\,+\, (\la_2 u_{2\ell} -0)^2,~~~~\ell=1,\ldots,g=1000~~\] 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.
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:
\vskip.5em
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.
\vskip.5em
2) $\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.
\vskip.5em 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.
The supplementary material includes additional theoretical and algorithmic details.