\section{Introduction}


Developing effective methods for quantifying and comparing empirical diversity of various biological populations is one of the fundamental problems  of modern life sciences as it has  direct impact on our understanding of the basic  operating principles of our planet's  ecosystem and its evolution \citep[cf., eg.,][]{Berkov:2014fk}.
In the  course  of its 3.5 billion years of evolutionary history, nature has developed an outstanding bio- and molecular diversity among  the Earth's  species of plants and animals.   Indeed,  it is  estimated that there are currently about 8.7 million eukaryotic species on earth,  both  marine and 
terrestrial,    88\% of which are still waiting to be described \citep{Mora:2011xy}.  The diversity  at the molecular level is 
perhaps even more spectacular,   as it occurs at different levels of biological organization: 
within one individual (e.g., through RNA, DNA, proteins, and metabolites), between individuals of the same and related species,  within and between species  and ecosystems, as well as  throughout evolution \cite[see, e.g.,][]{campbell2003save}.
For instance,   the  number of different molecular types of human T-cells is estimated at $10^{18}$ \citep{Jane05} which  only slightly less that the currently estimated number of stellar objects in the known universe (the latter believed to be of  the order $10^{21}$).  
 
Whereas the power of modern  computing has  allowed us  to make  steady progress towards  building  ever more  robust empirical measures of biodiversity  based on  a variety of  considerations \cite[see, e.g.,][]{Presley:2014fr},  the most relevant to our discussion here are the  measures  borrowed from  the field of information theory.   They  include among others  the   {\em Hill number}  (or the effective number of species) and the related concept of  the {\em Renyi  entropy}   (see, e.g., the recent review \cite{chiu2014phylogenetic} and references therein). Although  originally proposed  for  quantifying ecological diversity in the macro-scale ecosystems \citep{Chao:2010rt},   the use  of the empirical Renyi entropy as a descriptor of diversity  was also adopted  for  molecular populations in \citet{DeAndrade11}.     Since then   the Renyi-type  measures were   applied to  problems of molecular populations ranging from  analyzing   regulatory variants   and   testing   genome-wide associations \citep{Sun13,Sadee:2014aa} to comparing  different   T-cell populations  \citep{Cebula2013,RS12}.  Despite their growing usage in biodiversity studies of  both macro- and molecular- level populations,  it appears that  some important   statistical properties  of  the Renyi-type  measures have not been yet sufficiently understood, especially in the context of frequency-based   analysis and large sample behavior.

Currently,   standard  methods of obtaining  molecular level data on  the {\em transcriptome} (RNA)  abundance  rely on the so-called  next-generation sequencing (NGS) technology and especially    the high-throughput  RNA sequencing or RNA-seq \citep{Wang:2009aa}. However, the   molecular count  data  from NGS often elude  standard  statistical  analysis  due to the fact that exhaustive sampling of the DNA and RNA fragments for the purpose of sequence reconstruction  is  not feasible and that the  sequencing errors increase  with  sampling intensity or  {\em sequencing depth} \citep{ORawe:2015kq}.    It has been therefore generally conceded  \citep{Oh:2014aa}  that  the  standard,   fixed-dimension,  non-parametric frequency/contingency table  analysis   (see, e.g., \citealt{AA02}) does not readily apply  to the NGS data and that a different, {\em infinite-size}  contingency table framework,   more  reflective of the current sequencing technology,  appears  necessary.  Due to the  nature of the  NGS methods,   such  framework  should be based on the large sample (high-throughput)  considerations  but,  at the same time,  should also  account   for  the  increase in   the number of sequencing errors   with  increasing sample size  as well as   for the  under-sampling bias.   

 Motivated by  the questions on comparing biodiversity in molecular data (especially arriving from  the NGS experiments)  in  the current paper we establish   some  large sample results for the empirical Renyi entropy and divergence  in order to  bridge the gap between  current heuristic approaches   and  a more formal statistical  theory of large samples.   To this end, we  derive herein several central  limit theorems (CLTs)  which  yield approximate confidence bounds  for the  (Renyi) entropy-based  measures of diversity  and similarity in  the setting of  an infinite contingency table. Our CLT results   complement both the  law of large number theorems in \cite{RS12}  as well as the  CLT   for the plugin estimates of the Shannon entropy \cite{zhang2012normal} and the Kullback-Leibler divergence estimates  \citep{paninski2003estimation,zhang2014nonparametric}.  Since  in the  NGS experiments   one typically expects to under-sample the transcriptome, we focus here on   the  Renyi entropy exponent  (which below is denoted by $\alpha$) less than one, so as to   up-weight the  contributions of the lower  counts and  our  CLT results are restricted to this case. The extensions to  arbitrary exponents are straightforward but  not considered here.  In order to provide examples of the types of  applications  motivating   the mathematical results, we  analyze   two real  biological datasets from two different types of  NGS experiments.  In the first experiment, described  in   the  study \cite{Cebula2013},   one compares   multiple  T-cell receptors populations taken from mice before and after treatment with antibiotics.  The goal of  the  second  experiment  is the  elucidation  of   differences in  gene expression profiles between cancer and  control tissues  in individuals with hepatocellular carcinoma,  as described  in \cite{HCC14}. In both presented examples  the NGS datasets are analyzed  and de-noised by applying   a multi-stage process  developed on the basis of our theoretical results. 
 
As already indicated above, the problem of empirically estimating entropy and divergence  has been extensively studied in the statistical and machine learning literature  over past several decades, both  in the context of  discrete  and continuous   distributions.  See, for instance,  the monograph by \cite{pardo2005statistical} or the review  in \cite{Krish14} for more details.  In the  general case of   Renyi's  entropy and closely related Tsallis' entropy  of  a fixed continuos  distribution $f$ in ${\mathbb R}^m$,   a class of consistent estimators   was  proposed in \cite{leonenko2008class}  based on the $k$-th nearest-neighbor distances computed from the appropriate random samples of size $n$ from $f$.  The idea was later also extended to the Renyi entropy functionals in \cite{kallberg2012statistical} and it appears that similar results could be expected  to hold in the discrete case as well.  The main difference between these types of results and what is considered here is that in our setting  the discrete density function $f$ is allowed to change as the  sample size $n$ increases. Additionally, although   in the current and that we only analyze the basic  empirical frequency (the so-called plug-in) estimates.  

 
The paper is organized as follows. In the next section  (Section~2) we outline    the relevant mathematical   concepts  along with the necessary notation. In  Section~3 we state the main theoretical results of the paper, namely the CLTs for the Hill number (or the Tsallis entropy) and  the Renyi entropy and divergence  in the  asymptotic regime when the diversity of the population (i.e., the number of different types) grows with the sample size.  The results for the simpler case (Theorems 1 and 2) when  Renyi entropy statistics admit linear approximations are  established via the intermediate CLT results for the corresponding power sums which are closely related to the CLTs for  Hill's numbers and Tsallis' entropies. These results are also  included  as parts of formulations of Theorems 1 and 2. In case of  the uniform distribution for  the  Renyi entropy  as well as   the equal-marginals bivariate distribution for  the Renyi divergence, the power sum CLTs are no longer valid (there is no linear approximation available) and other methods are required to establish weak convergence to Gaussian variates  under slightly more stringent conditions. These results  are  presented  as Theorems 3 and 4 in Section~3. As it turns out, the key ingredient needed to establish Theorems 3 and 4 is the CLT result for two  Pearson-type  chi-square statistics in an  infinite contingency table. This latter result is of interest in itself and is presented as  Lemma~2 in Section~3.   In the following Section~4, we provide some simulation-based  examples  of the asymptotic behavior  of estimates from Section~3 in the case (relevant for our applications)  of  power law distributions under   various  sampling scenarios.  These examples illustrate  in particular how the  CLTs of Section~3  may  hold or not, depending on the relations between the dimensions of the relevant contingency tables and the empirical sample sizes.  In the second part  of Section~4  we  also discuss in detail the two biological examples of NGS data analysis and  show how the results of Section~3 may be used to analyze biodiversity of T-cell receptors  and to  profile the multiple sets of transcriptomes.  The final Section~5 offers a summary and brief conclusions. The proofs of all more complicated results are provided in the appendix along with some auxiliary technical lemmas. 


\section{Power Sums,  Entropy and Divergence}\label{sec:2}
Consider a triangular array of bivariate row-wise independent  random variables  $Z_{n,k}$ for $k=1,\ldots,n$ which in each row are equidistributed with the random variable  $Z_n=(X_n,Y_n)$  such that $P(X_n=i,Y_n=j)=p_{ij}^{(n)}$ for $i,j=1,\ldots,m_n$. Below we suppress the index $n$  when possible,  writing e.g., $m, Z_k, Z, p_{ij}$, etc.  for simplicity.  

Let  $\alpha>0$ and for any  probability distribution  ${\bm p}=(p_i)_{i=1}^m$ 
 define
 \begin{equation}\label{eq:1} 
 {\cal S}_\a({\bm p}) =\sum_{i=1}^m {p}_i^\alpha.
 \end{equation}
 Similarly, for any pair  of distributions ${\bm p}=(p_i)_{i=1}^m$ and ${\bm q}=(q_i)_{i=1}^m$
 define 
  \begin{equation}\label{eq:2} {\cal S}_\a({\bm p},{\bm q}) =\sum_{i=1}^m {p}_i^\alpha q_i^{1-\alpha}.
  \end{equation} 
   (Note that  ${\cal S}_1\equiv 1$). 
   The well-known special case of the above is  $\alpha=1/2$,  which results in  a symmetric index  ${\cal S}_{1/2}({\bm p},{\bm q})={\cal S}_{1/2}({\bm q},{\bm p})$ often referred to as  the  Bhattacharyya coefficient  \citep[see, e.g.,][]{nielsen2011burbea}.
   
    Recall   \citep{Renyi63} that for a given distribution ${\bm p}$ its  Renyi entropy $\mathcal {H}_\a$  is defined as
   $$\mathcal {H}_\a({\bm p}) =\frac{1}{1-\alpha} \log\left( \sum p_i^\alpha\right)=\frac{1}{1-\alpha} \log {\cal S}_\a({\bm p}) $$ and that for  a pair 
   of  distributions $({\bm p},{\bm q})$  their   Renyi divergence $\mathcal {D}_\a$ is  defined  as 
   $$\mathcal {D}_\a({\bm p},{\bm q}) = \frac{1}{\alpha-1}\log {\cal S}_\a({\bm p},{\bm q}).$$  Note that the sign change in the normalizing constant is needed in order to ensure non-negativity of $\mathcal {H}_\a$ and $\mathcal {D}_\a$. The special case of $\mathcal {D}_\a$ with $\alpha=1/2$ is referred to as the Bhattacharyya distance,  and may be expressed in terms of the   Mahalanobis distance \citep[see, e.g.,][]{nielsen2011burbea}, whereas      the linear approximation of $\mathcal {H}_\a({\bm p})$  given by 
\begin{equation}\label{eq:tse}\mathcal{T}({\bm p})=\frac{1}{1-\alpha}({\cal S}_\a({\bm p})-1).\end{equation} 
 is sometimes referred to as the  Tsallis entropy  and has important applications in the field of statistical mechanics  \citep{tsallis1988possible}.   
 Note that for our current purposes,  we  will only consider the  quantities $\mathcal {D}_\a,\mathcal {H}_\a$, and $\mathcal{T}$  for  $\alpha$ satisfying $0<\alpha<1$. 
 
 In what follows the summation symbol  without subscripts ($\sum $) will indicate  summation with respect to the index $i$   ($i=1,\ldots,m$)  whereas ${\bm p}=(p_i)_{i=1}^m$ and ${\bm q}=(q_i)_{i=1}^m$ will  (typically) denote the marginal distributions  of  the bivariate variable $Z=(X,Y)$ whose distribution is denoted by $(p_{ij})_{i,j=1}^m$.  Additionally, the  uniform distribution on $m$ points will be denoted by ${\bm u}$.  An important   relation  between  the Renyi entropy and the Renyi divergence   is 
 \begin{equation}\label{eq:hud}
 \mathcal {H}_\a({\bm p})=\log m -\mathcal {D}_\a({\bm p},{\bm u}).
 \end{equation} 
 We  note   also the following monotonicity property of $\mathcal {D}_\a$ and $\mathcal {H}_\a$ with respect to the index $\alpha$. 
 \begin{lemma} For $0<\alpha< \beta<1$ we have 
 $\mathcal {D}_\a({\bm p},{\bm q})\le {\cal D}_\beta({\bm p},{\bm q}) $ and thus, in view of \eqref{eq:hud}, also $ \mathcal {H}_\a({\bm p})\ge {\cal H}_{\beta}({\bm p})$.
 \end{lemma} 
 \begin{proof} Note that  for $x \ge 0$ the function $x \rightarrow x^\frac{\alpha-1}{\beta-1}$ is strictly convex for $0<\alpha<\beta<1$. Therefore, by Jensen's inequality  \begin{align*}\mathcal {D}_\a({\bm p},{\bm q}) &=\frac{1}{\alpha-1} \log \sum p_i^\alpha q_i^{1-\alpha} =\frac{1}{\alpha-1} \log \sum p_i\left(\frac{q_i}{p_i}\right)^{(1-\beta)\frac{\alpha-1}{\beta-1}}\\ &\le  \frac{1}{\beta-1} \log \sum p_i\left(\frac{q_i}{p_i}\right)^{(1-\beta)}={\cal D}_\beta({\bm p},{\bm q}).\end{align*} 
 \end{proof}
\begin{example}[Hill's Number] For given $0<\alpha<1$  the measure of diversity of a   distribution  ${\bm p}$  also known as the  {\em effective number of classes } may be defined as  \citep[see, e.g.,][]{jost2007partitioning,Chao:2012mz,RS12}  $ENC_\alpha({\bm p})=\exp (\mathcal {H}_\a({\bm p}))={\cal S}_\a({\bm p})^{1/(1-\alpha)}$. It follows then from Lemma~1 that for any $0<\alpha< \beta<1$ we have 
 $ENC_\alpha({\bm p}) \ge ENC_\beta({\bm p})$. (As it turns out, this inequality may be in fact extended to arbitrary positive $\alpha<\beta$). \end{example}

\subsection{Low Diversity  Condition and Projection Variables}\label{sec:ld}
The notion of an infinite-dimension contingency table  brought up in the introduction may be now  formally introduced simply  as a requirement that for $n$-size sample from $(p_{ij})_{i,j=1}^m$ we have $m\to \infty$ as $n\to \infty$. Throughout the paper, let $a\wedge b$ denote $\min(a,b)$ for any real $a,b$ and let $a_n\sim b_n$ (resp. $a_n\sim O(b_n)$) denote $a_n/b_n\to 1$ (resp. $A<\limsup_n a_n/b_n<B$ for some finite  $A,B$) as $n\to \infty$   for any  real sequences $a_n, b_n$. 
Throughout the paper  we  consider  only  the  {\em  low diversity} (LD) schemes in which the  marginals ${\bm p},{\bm q}$,  of $Z$   satisfy  the following {\em LD condition}. \begin{equation}\label{eq:ld}   (n p_\ast)^{-1}=o(n^{-\tau})  \quad \text{for some}\quad \tau>0,  \end{equation} where $p_\ast=\min_i(p_i)\wedge\min_i(q_i)$.
Note that  since $p_\ast\le 1/m$  then  \eqref{eq:ld}  implies  in particular $m/n= o(n^{-\tau})$. As it turns out,  for many distributions  ${\bm p}$ the two conditions are in fact  equivalent, as seen  in the  following. 

\begin{example}[Power Law Model] Let ${\bm p}={\bm q}$ and 
assume that $p_i=H^{-1}(\beta,m)/(i^{\beta}l(i))$, $(i=1,\ldots,m)$ where $\beta>0 $, $l(x)$ is a non-decreasing slowly varying function (see, e.g., \citealt{soulier2009some}, chapter 1), and $H^{-1}(\beta,m)=1/\sum_{i=1}^m ({i^{\beta}l(i)})^{-1}$ is the normalizing constant. Note that if  $0<\beta<1$ then   $H^{-1}(\beta,m) \sim (1-\beta) l(m)/m^{1-\beta}$ and \eqref{eq:ld} 
is implied by  $m/n=o(n^{-\tau})$ since 
$$(n \min_i p_i)^{-1} \sim (1-\beta)^{-1}\frac{ m^\beta l(m)}{n m^{\beta-1} l(m)} = (1-\beta)^{-1}\frac{m}{n}.$$
\end{example}

For any $0<\alpha<1$ and a given pair $(m, n)$, let  us define  two random variables which will play an important role in the following section.  Let  $W_n^{(\alpha)}$ be   defined  as  
\begin{equation}\label{eq:w} P(W_n^{(\alpha)}=\alpha p_{i}^{\alpha-1})=p_{i}\end{equation} for $i=1,\ldots,m$. Similarly,  define also  $V_n^{(\alpha)}$  as
\begin{equation}\label{eq:v}
P\left(V_n^{(\alpha)}=\alpha \left(\frac{q_i}{p_i}\right)^{1-\alpha} + (1-\alpha) \left(\frac{p_j}{q_j}\right)^{\alpha}\right)=p_{ij} \end{equation} for $i,j=1,\ldots,m$. 
In the following,  for the reasons discussed below, we refer to \eqref{eq:w} and \eqref{eq:v} as the {\em projection variables} or simply {\em projections}. 
\begin{remark}\label{rem1}   Note that $$E W_n^{(\alpha)}=\alpha{\cal S}_\a({\bm p}) $$ and $Var W_n^{(\alpha)}=0$ iff $p_i=1/m$ for all $i$, that is,  ${\bm p}=(p_i)={\bm u}$ is a uniform distribution   on  $m$ support points (this case  is often referred to as  a maximal diversity model or a pure noise model).  Similarly, $$E V_n^{(\alpha)} = {\cal S}_\a({\bm p},{\bm q})$$ and it is also easy to see that 
$Var V_n^{(\alpha)}=0$ iff $p_{i}=q_i$ for all $i$, that is, ${\bm p}={\bm q}$.
\end{remark}
As it turns out, both   cases ${\bm p}={\bm u}$ and ${\bm p}={\bm q}$  require special consideration in  the asymptotic analysis of $\mathcal {H}_\a$ and $\mathcal {D}_\a$. In view of the remark above they may be referred to as the cases of ``degenerate" (zero variance)  projections.  
\begin{example}[Noise--and--Signal and Pure Noise Models]\label{ns1}
 A distribution concentrated on $m+1$ support points, such that $p_0>0$ and $p_i=(1-p_0)/m$ for $1\le i\le m$, may be considered as a simple model of  signal contamination. Note that in this case we have 
$P(W_n^{(\alpha)}=\alpha p_0^{\alpha-1})=p_0$, $P(W_n^{(\alpha)}=\alpha m^{1-\alpha}(1-p_0)^{\alpha-1})=1-p_0$ and 
$$ Var W_n^{(\alpha)} = \alpha^2\left( m^{1-\alpha}(1-p_0)^\alpha\left(\frac{p_0}{1-p_0} \right)^{1/2}-p_0^\alpha\left(\frac{1-p_0}{p_0} \right)^{1/2}\right)^2.$$ 
For the pure noise model $p_0=0$,  in which case the support reduces to $m$ points, and the above formula is not valid. However,  as already pointed out  before, in this case  we may show directly that $Var W_n^{(\alpha)}=0$. \end{example}
 \section{Limit Theorems}\label{sec:lims} Let $N(0,1)$ denote the standard Gaussian random variable and $\Rightarrow$ denote the usual weak convergence in the  space of  probability  distributions.   Define also the plug-in $n$-sample estimates of ${\bm p}$ and ${\bm q}$ as, respectively,  
$\bhp=(\hat{p}_i)_{i=1}^m,$ where $ \hat{p}_i=\sum_{k=1}^n I(X_k=i)/n$ and  $\bhq=(\hat{q}_i)_{i=1}^m,$ where $ \hat{q}_i=\sum_{k=1}^n I(Y_k=i)/n$. Here and elsewhere in the paper $I(\cdot)$ denotes the indicator function.
As it turns out,  
two distinct sets of CLTs may be derived depending on whether the variables $W_n^{(\alpha)}$ and $V_n^{(\alpha)}$ are degenerate  (that is, their respective variances vanish) or not.  For the non-degenerate case  the appropriate CLTs may be established  by expanding on  the   usual projection and Taylor's expansion arguments \citep[see, e.g.,][chapter 1]{Shao04}. This is the simpler case to consider and we discuss it first. 
\subsection{CLTs for Non-Degenerate Projections}
The first two CLT results  for the empirical (plug-in) Renyi entropy and divergence  and their  corresponding power sums  are provided in Theorems~1 and 2 below. Their respective hypotheses  $(iii)$  may be viewed  as complementing  the analogous results established for the Shannon entropy and  the Kullback-Leibler divergence \citep{paninski2003estimation,zhang2012normal,zhang2014nonparametric}.
 Note also that   
${\cal S}_\a=(ENC_\alpha)^{1-\alpha}$ where the Hill number $ENC_\alpha$ is defined in Example~1.  The proofs are deferred  to the appendix. 

Recall that for any  square integrable random variable $X$, such that $EX\ne 0$,  we define its coefficient of variation as ${\cal CV}(X)=(Var X)^{1/2}|E X|^{-1}$.  \begin{theorem}[Renyi Entropy CLT]\label{thm:1}  Let $W_n^{(\alpha)}$ be a sequence of random variables defined by \eqref{eq:w} with  $\inf_n \mathcal{CV}(W_n^{(\alpha)})>0$ and   let \begin{equation}\label{eq:pc} \sum p_i^{\alpha-1}(n Var W_n^{(\alpha)})^{-1/2}\to 0\quad\text{for}\   m,n\to \infty.\end{equation}
Then, under the LD condition \eqref{eq:ld}, as $m,n\to \infty$
\begin{itemize}
\item[(i)] ${\cal S}_\a(\hat{\bm p})/{\cal S}_\a({\bm p}) \to 1$ in probability, 
\item[(ii)] 
$\sqrt{n}({\cal S}_\a(\hat{\bm p})-{\cal S}_\a({\bm p}))/(Var W_n^{(\alpha)})^{1/2}\Rightarrow N(0,1),$
\item [(iii)] $\sqrt{n}\,(1/\alpha-1)(\mathcal {H}_\a(\hat{\bm p})-\mathcal {H}_\a({\bm p}))/\mathcal{CV} (W_n^{(\alpha)})\Rightarrow N(0,1).$
\end{itemize}
\end{theorem}
\begin{remark}Note that the first two assertions of the theorem may be equivalently stated in terms of the convergence of the Tsallis plug-in entropy defined by \eqref{eq:tse}.
\end{remark}
\begin{remark}\label{rem2}
Note that the  condition  \eqref{eq:pc} is typically stronger than \eqref{eq:ld}. Indeed, taking $\alpha>1/2$ and  the power law model from Example~2 with $0<\beta<1$ we obtain
 $\sum p_i^\alpha\sim (1-\beta)^\alpha\, m^{1-\alpha}/(1-\alpha\beta)$ and $\sum p_i^{2\alpha-1}\sim (1-\beta)^{2\alpha-1}\, m^{2-2\alpha}/(1-2\alpha\beta+\beta)$. Consequently, for some constant  $C>1$ 
$$ \frac{C\sum p_i^{\alpha-1} }{\sqrt{n (\sum p_i^{2\alpha-1} -(\sum p_i^\alpha)^2) }}\ge \frac{m}{\sqrt{n}} \frac{(\max_i p_i)^{\alpha-1}}{m^{1-\alpha}}\ge  \frac{m}{\sqrt{n}}$$   for large $m,n$ and   \eqref{eq:pc} implies  \eqref{eq:ld} with $\tau=1/2$. Similarly, (possibly for different $C>1$) 
$$ \frac{\sum p_i^{\alpha-1} }{\sqrt{n (\sum p_i^{2\alpha-1} -(\sum p_i^\alpha)^2) }}\le  \frac{C m}{\sqrt{n}} \frac{(\min_i p_i)^{\alpha-1}}{m^{1-\alpha}}\le C \frac{m}{\sqrt{n}} $$
and therefore  in this case \eqref{eq:pc} is seen to be actually equivalent to \eqref{eq:ld} with $\tau=1/2$.
\end{remark}
\begin{remark}[Plug-in Bias]\label{rem2a} Note that, in view of Jensen's inequality applied to the  strictly concave function $x\rightarrow x^\alpha$ for $x>0$ and  $0<\alpha<1$, we have 
$ E{\cal S}_\a(\hat{\bm p})/{\cal S}_\a({\bm p})\le 1$. This and  the assertion $(i)$ above imply together that under the assumptions of Theorem~1 the {\em relative bias}  of  ${\cal S}_\a(\hat{\bm p})$ satisfies   $E{\cal S}_\a(\hat{\bm p})/{\cal S}_\a({\bm p})-1\to~0$
as $n,m\to \infty$. The standard  inequality  $\log x\le x-1$ valid for $x>0$ implies then that 
the bias of the plug-in entropy estimate satisfies \begin{equation}\label{eq:bias} E \mathcal {H}_\a(\hat{\bm p})-\mathcal {H}_\a({\bm p})\to 0\qquad \text{as}\quad n,m\to\infty. \end{equation} \end{remark}
Unfortunately,  as may be  seen from the proof of Theorem~1 in  the appendix,  a more careful analysis of the tail events for  the plug-in estimate  than the one currently performed   is needed in order to actually establish  a convergence rate in \eqref{eq:bias}.

Turning now to our second result, note that the relation \eqref{eq:hud} suggests that  CLT of Theorem~1 could be also extended to the Renyi divergence.  The proof is again based on the Taylor expansion method where now the projection variable \eqref{eq:w} is  replaced by \eqref{eq:v}. 
\begin{theorem}[Renyi Divergence CLT]\label{thm:2} Let $V_n^{(\alpha)}$ be a sequence of random variables defined by \eqref{eq:v} with $\inf_n\mathcal{CV} (V_n^{(\alpha)})>0$ and let 
\begin{equation}\label{eq:pc2}\left( \sum \left({q_i}/{p_i}\right)^{1-\alpha} +\sum \left({p_i}/{q_i}\right)^{\alpha}\right) (nVar V_n^{(\alpha)})^{-1/2}\to 0\quad\text{for}\ m,n\to \infty.\end{equation}
Then, under the LD condition \eqref{eq:ld}, as $m,n\to \infty$
\begin{itemize}
\item[(i)] $\displaystyle{{\cal S}_\a(\hat{\bm p},\hat{\bq})/{\cal S}_\a({\bm p},{\bm q}) \to 1}$ in probability, 
\item[(ii)] $\displaystyle{\sqrt{n}({\cal S}_\a(\hat{\bm p},\hat{\bq})-{\cal S}_\a({\bm p},{\bm q}))/(Var V_n^{(\alpha)})^{1/2}\Rightarrow N(0,1)},$
\item [(iii)] $\displaystyle{\sqrt{n}\,(\alpha-1)(\mathcal {D}_\a(\hat{\bm p},\hat{\bq})-\mathcal {D}_\a({\bm p},{\bm q}))/\mathcal{CV} (V_n^{(\alpha)})\Rightarrow N(0,1).}$
\end{itemize}
\end{theorem}

\begin{remark}[Plug-in Bias]\label{rem2b} Note that, similarly as in Remark~\ref{rem2a}, we have  
$ E{\cal S}_\a(\hat{\bm p},\hat{\bq})/{\cal S}_\a({\bm p},{\bm q})\le 1$  and, by a similar argument as before,   Theorem~2$(i)$ implies $$E\mathcal {D}_\a(\hat{\bm p},\hat{\bq})-\mathcal {D}_\a({\bm p},{\bm q})\to~0\qquad\text{as}\ n,m\to \infty.$$ \end{remark}
\begin{example}[Symmetric Divergence for Power Laws] Consider the symmetric divergence ${\cal D}_{1/2}({\bm p},{\bm q})$  with independent marginals,  which often is the  case  of interest in  NGS applications. Note that in this situation $Var V_n^{(1/2)}=1/2-(\sum \sqrt{p_i q_i})^2/2$. Suppose additionally that   $p_i=H^{-1}(\beta_1,m)/(i^{\beta_1}l_1(i))$ and  $q_i=H^{-1}(\beta_2,m)/(i^{\beta_2}l_2(i))$, $(i=1,\ldots,m)$ where the notation is as in Example~2 with $0<\beta_1\ne\beta_2<1$. Then $$Var V_n^{(1/2)}\sim \frac{1}{2}-\frac{\sqrt{(1-\beta_1)(1-\beta_2)}}{2-\beta_1-\beta_2}$$ and, consequently,  \eqref{eq:pc2} is seen as  equivalent to $m/\sqrt{n}\to 0$  (cf. also Remark~\ref{rem2} above). 
\end{example}

 With some additional effort, the two CLT results of this section may be extended to degenerate projections. This is discussed in the next section. 
 
\subsection{CLTs for Degenerate Projections }
In case of  a  degenerate projection,    the linear term of the power sum Taylor's expansion disappears (cf. formula 
(B.6) in the appendix) and  the condition \eqref{eq:pc} is no longer needed. However,  the  LD  assumption \eqref{eq:ld} has  to be slightly strengthened in order to establish the asymptotic results  for the leading (quadratic) term of the appropriate expansion. 
\subsubsection{Chi-Square Statistic CLT}
The following lemma describing  the chi-square statistic CLT  may be of independent interest for models of  sparse contingency tables. For a recent discussion of a normal approximation to the chi-square statistic in such settings, see, e.g.,  \cite{chi2approx2013}.  Here  we  apply the chi-square  CLT  formulated below     to obtain   weak limits  for the quadratic  terms in the  entropy and divergence Taylor's expansions  leading to Theorems~3 and 4 described in the next subsection. To begin, consider    
a pair of distributions $({\bm p}, {\bm q})$  and a set of positive weights ${\bm r}=(r_i)_{i=1}^m$ and  define the corresponding  chi-square ($\chi^2$) distance function as   $${\cal X}^2_{\bm r}({\bm p},{\bm q})=n \sum \frac{(p_i-q_i)^2}{r_i}.$$
Note that,  for instance,   the $\chi^2$-distance  statistic between the empirical marginals $(\hat{\bm p},\hat{\bq})$ is  obtained by setting  $r_i=p_i+q_i$ 
$${\cal X}^2_{\bm r}(\hat{\bm p},\hat{\bq})=n\sum \frac{(\hat{p}_i-\hat{q}_i)^2}{p_i+q_i}$$ 
and  the  Pearson $\chi^2$-statistic  is  obtained by setting $r_i=p_i$ 
\begin{equation}\label{eq:ch1}
\ch_\bp(\hat{\bm p},{\bm p})=n\sum\frac{(\hat{p}_i-p_i)^2}{p_i}.
\end{equation} 
  Below we denote   ${\cal X}^2_{\bm u}(\hat{\bm u},{\bm u})=:{\cal X}^2_{\bm u}$.  


\begin{lemma}\label{lem:ch1}   Let $(p_{ij})_{i,j=1}^m$ be the bivariate distribution of $Z=(X,Y)$  with  $X$ and $Y$  having   marginals  $(p_i)_{i=1}^m$ and $(q_i)_{i=1}^m$  where $p_i=q_i>0$. Assume $m\to \infty $ as $n\to \infty $ and 
\begin{align}
& (mn)^{-1}\sum \max (p_i^{-1},p_i^{-2}m^{-1})\to 0,  \label{eq:con2}
\end{align}
Then  as  $n\to \infty$
\begin{itemize}
\item[(i)] $\displaystyle {\frac{{\cal X}^2_{\bm p}(\hat{\bm p},{\bm p})-m}{\sqrt{2m}}\Rightarrow N(0,1)},$ \\  and if additionally 
\begin{equation}\label{eq:bd}
\sup_n \max_{ij} \frac{p_{ij}}{p_ip_j}=B<\infty 
\end{equation} then also \vspace{.071in}
\item[(ii)] $\displaystyle {\frac{{\cal X}^2_{2{\bm p}}(\hat{\bm p},\hat{\bq})-\mu_n }{\sqrt{2}\gamma_n}\Rightarrow N(0,1)},$ \\ where \begin{align} \mu_n &=\sum_i (1-p_{ii}/p_i)\nonumber \\  \gamma_n^2 &=\sum_i \frac{(p_i-p_{ii})^2}{p_i^2}+\sum_{1\le i\ne j\le m} \frac{(p_{ij}+p_{ji})^2}{4p_ip_j}.\label{eq:sn}\end{align}
\end{itemize}
\end{lemma}

 \begin{remark}\label{rem3}  Note that for ${\cal X}_{\bu}^2$ the condition \eqref{eq:con2} simplifies  to $m/n\to 0$.
 \end{remark}
 
  \begin{remark}\label{rem4}  Note that under the assumption \eqref{eq:bd} we have $m -2B\le \gamma_n^2\le m+B^2$ and therefore $\gamma_n^2\sim m$. In particular,  if $p_{ij}=p_ip_j$ then  $\mu_n=\gamma_n^2=m-1$.
\end{remark} The proof of the result may be found in the appendix. 
Its application  is  discussed next.  
\subsubsection{Pure Noise and Equal Marginals CLTs}
  The first result  covers  the case of Renyi entropy  when ${\bm p}={\bm u}$. The proof is  outlined in the appendix. Recall that for real $a$ and integer $k$ we define $\binom{a}{k}=a(a-1)\cdots (a-k+1)/k!$ 
\begin{theorem}[Uniform Entropy CLT]\label{u1} Assume $m\to \infty $ as $n\to \infty $ and  $m^2/n =o (n^{-\tau})$ for $\tau>0$. Then \begin{itemize}
\item[(i)] $ \frac{n\binom{\alpha}{2}^{-1}[m^{\alpha-1}{\cal S}_\a(\hat{\bm u})-1]-m}{\sqrt{2m}}\Rightarrow N(0,1)$
\item[(ii)] $   \frac{n[\mathcal {H}_\a(\hat{\bm u}) -\log m- (1-\alpha)^{-1}{\log(1+\binom{\alpha}{2}\frac{m}{n}})]}{\alpha\sqrt{m/2}}\Rightarrow N(0,1).$
\end{itemize}
\end{theorem}


Our second CLT result is  the following theorem for Renyi divergence  when  ${\bm p}={\bm q}$. The proof is  again deferred to  the appendix.

 \begin{theorem}[Degenerate Divergence  CLT]\label{u2} Let $(p_{ij})_{i,j=1}^m$ be the bivariate distribution of $Z=(X,Y)$  with  $X$ and $Y$  having    marginals  ${\bm p}=(p_i)_{i=1}^m$ and ${\bm q}=(q_i)_{i=1}^m$  where $p_i=q_i>0$. Let $\mu_n$ and $\gamma_n^2$ be given by \eqref{eq:sn}. Assume $m\to \infty $ as $n\to \infty $ and that \eqref{eq:bd} holds,  as well as that 
\begin{equation}\label{eq:cond3}  \max \left\{\frac{1}{nm\min p_i^2},\frac{m}{n\min p_i}\right\}=o(n^{-\tau}). 
\end{equation} Then \begin{itemize}
\item[(i)] $ \frac{n(\alpha(\alpha-1))^{-1}[{\cal S}_\a(\hat{\bm p},\hat{\bq})-1]-\mu_n}{\sqrt{2}\gamma_n}\Rightarrow N(0,1)$
\item[(ii)] $   \frac{n[\mathcal {D}_\a(\hat{\bm p},\hat{\bq})- (\alpha-1)^{-1}\log(1+\alpha(\alpha-1)\frac{\mu_n}{n})]}{\alpha\sqrt{2}\gamma_n}\Rightarrow N(0,1).$
\end{itemize}
\end{theorem}
\begin{remark} Note that for ${\bm p}={\bm q}={\bm u}$ the condition \eqref{eq:cond3} reduces to $m^2/n=o(n^{-\tau})$ required in  Theorem~3.

\end{remark}



\subsubsection{Random Sample Size}
When analyzing  NGS  data   some  part of the sequences reads   is frequently removed for technical reasons, for instance, due to poor amplification or reading  errors (see next section). In such cases  one effectively  deals with  a molecular  sample of random size.    Our  CLT   results derived  earlier may be extended to this case as well, with  the help of following simple result described in Theorem~5 below. Its  various versions have been discussed, for instance,  in the context of random allocations \citep[see, e.g.,][]{Kolchin78}. 

\begin{theorem}[Randomized  Sample CLT]\label{rcpt}  Let $(Z_n)_{n=1}^\infty$ be  a  sequence of bivariate variables    supported on an $m_n\times m_n$ integer lattice with distribution  $(p_{ij})_{i,j=1}^{m_n}$.   Let 
$(\hat{Z}_n)=(\hat{p}_{ij})_{i,j=1}^{m_n}$ $(n=1,2,3, \ldots,)$ be the sequence of the  empirical estimates, each  based on an iid sample of (deterministic) size $n$. Suppose that   the statistic  ${\cal G}_n={\cal G}_n(\hat{p}_{ij})$  satisfies   $b_n({\cal G}_n-a_n)\Rightarrow N(0,1)$ as $n\to \infty$ with some non-random  $(a_n, b_n)$.   Let    $(\nu_n)_{n=1}^\infty$ be a sequence of random variables independent of $(\hat{Z}_n)_{n=1}^\infty$  and  following  the binomial distributions $bin(n,\tau_n)$ with $0<\inf_n \tau_n\le \sup_n \tau_n<1$.  Then  also  $$b_{\nu_n}({\cal G}_{\nu_n}-a_{\nu_n})\Rightarrow N(0,1).$$
  \end{theorem}
  \begin{proof}
  Denote by  ${\cal G}_{n_k}$ the random variable ${\cal G}_{\nu_k}$   conditional on the event  $\nu_k=n_k$ and by  $\Phi$  the distribution function of the standard normal random variable. 
 By assumption,  for  any  real $x$   we have $ P({\cal G}_{n_k}\le x) \to \Phi(x)$ provided that $n_k\to \infty$ as $k\to\infty$.
 Let $\varepsilon>0$ be sufficiently small and define $C_\varepsilon(k_0)=\{n_k:  k(\tau_k-\varepsilon)\le n_k\le k(\tau_k+\varepsilon), k>k_0\}$. Note that  by the weak law of large numbers  $P(\nu_k \in C_\varepsilon(k_0))\to 1$ as $k_0\to \infty$. Therefore 
\begin{align*} 
P({\cal G}_{\nu_k}\le x, \nu_k\in C_\varepsilon(k_0))  &= \sum_{n_k\in C_\varepsilon(k_0)} P({\cal G}_{n_k}\le x) P(\nu_k=n_k) \\
&= (\Phi(x)+\delta(k_0)) P(\nu_k \in C_\varepsilon(k_0))
\end{align*} where $\delta(k_0)\to 0$ as $k_0\to \infty$. Accordingly,  as $k_0\to \infty$ the left-hand side converges to $\lim_k P({\cal G}_{\nu_k}\le x)$ and  the right-hand side to $\Phi(x)$ and the result follows.
   \end{proof}

\section{Examples and NGS Applications}
We start by   providing some  numerical examples illustrating   that,  in general, the CLT results discussed above do not  hold without  assumptions on the relative rate of $m$ and  $n$. Next, we  show  two examples of  applicability of our results  to analyzing biodiversity of NGS data.  The first one is concerned with  comparing the diversity of T-cell receptor populations in transgenics  mice, whereas the second one aims at identifying the hepatocellular carcinoma  transcription profiles  in humans.  For the purpose of the   T-cell receptors   example, we  propose  a  sequential statistical procedure of NGS signal filtering based on our CLT results from the previous sections. We begin by  pointing out  to some subtleties   in   the CLT results discussed in Section~3.
\subsection{Power Law and Pure Noise Models}
Consider the power law model from Example~2 in Section~\ref{sec:ld} with $\beta=1$ and $l(x)\equiv 1$. Note that  in this case 
$ (n \min_i p_i)^{-1} \sim m \log m/\,n$ as well as  $\sum p_i^{\alpha-1}(n VarW_n^{(\alpha)})^{-1/2} \sim O(m(\log^{2\alpha} m/n)^{1/2})$ and therefore the assumptions of Theorem~1 are satisfied as soon as 
\begin{equation}\label{eq:10}
n^{\tau-1} m\to 0
\end{equation} for some $\tau>1/2$. Similarly, the assumption \eqref{eq:con2} of Lemma~2 is satisfied as soon as  
\begin{equation}\label{eq:11}
\log^2\!m\,\frac{m}{n}\to 0.
\end{equation}
In Figure~\ref{fig1} we illustrate the convergence results of Theorem~1$(iii)$ and Lemma~2$(i)$ for   this power law model and $\alpha=0.5$.  The panels of Figure~\ref{fig1} presents the  sample vs standard normal quantile  (QQ) plots  for the normalized Renyi entropy statistic and  the normalized Pearson statistic \eqref{eq:ch1} based on  $B=5000$ samples from the power law distribution, each  with $m=1000$ and three different values of $n=m^{1+\varepsilon}$ ($\varepsilon=-0.5, 0.5, 1.5$). As seen from the plots, in the absence of  \eqref{eq:10} the CLT result for the Renyi entropy (cf. Theorem~1$(iii)$) does not hold. Moreover, the middle panel QQ plot  indicates that  for  large $m,n$  satisfying $n=m^{3/2}$ the discrepancy between distribution of  the entropy function  and its plug-in estimate appears in a form of deterministic shift, indicating the presence of substantial asymptotic bias and hence the lack of  convergence  \eqref{eq:bias}. Similarly,  when  \eqref{eq:11} is not satisfied than   the Pearson statistic CLT given in Lemma~2$(i)$ fails with the middle panel again indicating that the bias of the estimate does not vanish when $m$ is too large relative to $n$. 
\begin{figure}[htp]
 \centering{
\includegraphics[scale=0.6,trim=0 10 0 10]{Fig1}}
\caption{{\bf Projection CLTs.} Normal QQ  plots  for  the normalized Renyi entropy (Theorem~\ref{thm:1}$(iii)$, lower (green) curve) and normalized Pearson $\chi^2$ statistic (Lemma~2$(i)$,  upper (blue) curve)  for the power law distribution $p_i=1/i$. The panels shows quantile plots  with different values of $n=m^{1+\varepsilon}$ ($\varepsilon=-0.5, 0.5, 1.5$) and $m=1000$. The   solid  (red) line gives  quantiles of the standard normal distribution for reference. }\label{fig1}
\end{figure}

For comparison,  we also considered the uniform distribution (pure noise) model $p_i=1/m$. Note that it  may be  viewed as a degenerate power law where $\beta=0$ and $l(x)\equiv1$.  Recall that according to Theorem~3 $(ii)$ and Lemma~2 $(i)$, the sufficient conditions for the respective CLTs are  $m^2/n^{1-\tau}\to 0$ and $m/n\to 0$ (see Remark~\ref{rem3} for the latter one).  The necessity  of  these conditions is  illustrated in the panels of  Figure~\ref{fig2}  where we again present the (normal) QQ plots for the Renyi ($\alpha=0.5$) and the Pearson statistics   for the same values of $B, n$ and $m$ as in Figure~\ref{fig1}. As seen from these plots,  only in the last panel, when $m^2/n\approx  0$, we get good  CLT approximation for both statistics. These  results appear  consistent with our theoretical results from Theorem~3 and Lemma~2.
\begin{figure}[htp]
 \centering{
\includegraphics[ scale=0.6,trim=0 10 0 10]{Fig2}
\caption{{\bf  Degenerate Projection CLTs.} Normal QQ plots for  the normalized uniform Renyi entropy (Theorem~\ref{u1}$(ii)$, represented by the lower  (green) curve)  and the normalized  Pearson $\chi^2$-statistic (Lemma~2 $(i)$, represented by the upper (blue) curve) with $p_i=m^{-1}$. The panels shows quantile plots with  different values of $n=m^{1+\varepsilon}$ ($\varepsilon=-0.5, 0.5, 1.5$) and $m=1000$. The   solid  (red) line gives  the quantiles of the standard normal distribution for reference.  Note that the normalized Renyi entropy is undefined for   the first panel. }\label{fig2}
\end{figure}

Although not presented here due to space considerations, similar examples based on the  bivariate power laws  may be used to  illustrate  the necessity of the assumptions of type  \eqref{eq:pc2} and \eqref{eq:cond3} in the CLT results for divergence in Theorems~2$(iii)$ and 4$(ii)$. 

\subsection{Applications to NGS Data}
 Our CLT results described in Section~\ref{sec:lims} were originally motivated by questions rising in NGS data analysis. Below  we describe   two  examples which adhere to the following basic framework.   Denote by  $\bm \varepsilon_1,\bm \varepsilon_2$  two independent noise distributions  each on $m$ support points,    and assume   that a pair $({\bm p}, {\bm q})$ of marginal distributions may be represented as \begin{equation}\label{eq:nm}({\bm p},{\bm q}) =\lambda (\tilde{\bm p}, \tilde{\bm q})+(1-\lambda)(\bm \varepsilon_1,\bm \varepsilon_2)\end{equation} where $(\tilde{\bm p},\tilde{\bm q})$ is a pair of  marginal   distributions having no common support points with  $(\bm \varepsilon_1,\bm \varepsilon_2)$ and $\lambda$ is the mixing proportion (or prior probability of signal).   We  assume that  each $\bm \varepsilon$ is a simple finite mixture of  $K$  uniform distributions on separate support. Note that  the noise-and-signal model from Example~\ref{ns1} in Section~\ref{sec:ld} may be viewed as a (univariate) special case of \eqref{eq:nm} with $K=1$.  In the  first  example below we took $K=2$.   
 \begin{table}
\begin{center}
\begin{tabular}{ccc}
& Antibiotic ($\tilde{\bm p}$)& Control  ($\tilde{\bm q}$)\\
\hline
\tabT$n$ & 39,084&39,084\\
$m$ & 165&165\\
$k_m$ & 17 & 17 \\
$\hat{\lambda}$ &0.46 &0.46\\
$\hat{\beta}$ & 0.869(0.05) & 0.971(0.05) \\
${\cal H}_{1/2}$ &4.81 (4.79, 4.82) & 4.64 (4.63, 4.67) \\
$ENC_{1/2}$ & 122.73 (120.30, 123.97) & 103.54 (102.51, 106.70)\\
${\cal D}_{1/2}$ & \multicolumn{2}{c}{0.155 (0.147, 0.163)}  \\
\hline

\end{tabular}
\caption{{\bf Results of  TCR data analysis}. The mixture model \eqref{eq:nm} with   heavy-tailed power laws  fitted to two sets of TCR  counts derived from mouse MLN before and after an antibiotic treatment as described in \cite{Cebula2013}. The empirical Renyi entropy, the Hill number and the Renyi diversity CIs (in parenthesis)  are obtained from the CLT results of Theorems~1 and 2.}\label{tab1}
\end{center}
\end{table}

\vspace{.1in}\noindent {\bf Algorithm 1(NGS  Diversity Analysis with $\mathcal {D}_\a$ or ${\cal S}_\a$)}\vspace{-.1in}
\begin{cenum}
\item {\em Exponent ($\alpha$) selection}. Use  problem-specific criteria (e.g. sample coverage, see \cite{RS12}) to identify the appropriate $\alpha$ value. If no  prior knowledge  exist,  the value  $\alpha=1/2$  (the  Bhattacharyya distance) may be  often used. 
\item {\em Noise filtering}. Identify the number of mixture components $K$ and the cut-off count(s) $k_m$ for  the support of  $\bm \varepsilon_i$ in \eqref{eq:nm} with   a sequential (starting from the lowest empirical frequency) procedure  based on  Lemma~2$(i)$ with  ${\bm p}=\bm \varepsilon_i$ ($i=1,2$).    The values of $\lambda$ is then estimated as the proportion of a sample falling into the $m$ 'noise' categories.  
\item{\em Equality testing}. For a pre-determined value of $\alpha$,  test the hypothesis $H_0: \tilde{\bm p}=\tilde{\bm q}$ by  comparing the  observed value of $\mathcal {D}_\a$  (alternatively, ${\cal S}_\a$) with  the asymptotic normal distribution in  Theorem~4. 
\item{\em Difference quantification}. If $H_0$ is  not rejected,  conclude that $\mathcal {D}_\a\equiv 0$ (${\cal S}_\a\equiv 1$). Otherwise,   apply Theorem~2 to obtain confidence bounds 
for $\mathcal {D}_\a$ (${\cal S}_\a$).\end{cenum}

\subsubsection{T-Cell Receptor Populations} In this  example we apply Algorithm~1 to  measure similarity between   a pair of  T-cell receptor (TCR) populations based on the observed  NGS counts of  receptor-specific nucleotide sequences.   With the  current NGS  technology,  the two main difficulties in  comparing  TCR populations are  to adjust the under-sampling  bias due to unobserved rare types  and  the `ghost`  types  created   due to  the sequencing errors  \citep{Wang:2014aa}.    The first  problem may be often alleviated  by   applying  diversity criteria, like the Renyi entropy and divergence, which allow for the sample-based up-weighting of  rare  counts (see  \citealt{RS12}).   The second one  requires typically additional assumptions, in order to perform  analysis as outlined  in Algorithm~1$(ii)$.  A recent  detailed overview  of the   TCR diversity analysis  methods was presented  by   \cite{RS12} and  earlier on, in a more general context of biodiversity,  by  \cite{Hsieh:2006sy} and  \cite{Magurran:2005by}.   For illustration, we analyze here  two populations derived from the mesenteric lymph nodes (MLN) of a  TCR mini-mouse before and after an antibiotic treatment.  The details of the experiments and   a dataset description are  given in   \cite{Cebula2013}.  For the current analysis it is important to note that,  since the experimental groups  consisted of  different animals,  we may consider  two experimental groups  as independent.  The total combined  sample size (or sequencing depths)  was $n=72,030$,  with initial $m_0=6,336$ receptor types. After performing  step $(ii)$ of   Algorithm~1 $m=165$  types were  identified  as ``signal"  based on the  cut-off  $k_m=17$ in both populations. The signal population corresponded to the remaining sample size of $38,896$   or about  54\% of the original NGS counts. We used  $\mathcal {D}_\a$ with  $\alpha=1/2$ as the diversity measure  in step $(iii)$-$(iv)$  of Algorithm~1.   Based on Theorem~2,  the asymptotic  $P$-value for testing $H_0: \tilde{\bm p}=\tilde{\bm q}$ was found to  be less than $10^{-4}$ and hence the hypothesis of equal diversity of the two populations was rejected (see Algorithm~1$(iii)$).  

\begin{figure}[htp!] \centering{ \
\includegraphics[scale=0.4,trim=15 15 15 15]{Fig3}}
\caption{{\bf Power law fit for TCR data}. QQ plot of the TCR data against quantiles of a power law distribution with $\beta_1=0.87$ ($SE=.05$)  and $\beta_2=0.97$ ($SE=.05$) values fitted via  the least squares method. }\label{fig:qq} 
\end{figure}

To compare this finding with a more standard  parametric analysis, we additionally  fitted, with the least squares method, the counts of $165$ receptor types in two populations to the   power law distributions.  Since the respective  exponent  values  for the two  fitted populations were found  to be different, with    $\beta_1=.87$ (for antibiotic treated mice) and $\beta_2=.97$ (for untreated), the parametric analysis confirmed the findings of Algorithm~1.  For illustration,  the plots of the fitted power law quantiles versus  the empirical ones  are presented in Figure~\ref{fig:qq}.  Additionally, the diversity of each of the TCR populations in terms of its  respective Renyi entropy ${\cal H}_{1/2}$ and the Hill number $ENC_{1/2}$ as well as the diversity difference  measured by the Renyi divergence  ${\cal D}_{1/2}$  are listed  in Table~\ref{tab1}, along with the  corresponding asymptotic  confidence intervals  obtained via Theorems~1 and 2.  
 As seen from the values in Table~\ref{tab1},  although the diversity of each of the NGS populations was relatively similar  in terms of the two populations count patterns, it differed in terms of the specific TCR types expressed.    
 

\subsubsection{Gene  Expression Profiling}

Beyond Algorithm~1, the results of Section~3   may be applied to facilitate various other biodiversity analysis,   for instance, in    simultaneous comparison of  several pairs of  molecular  samples.  We  illustrate this  with an NGS   data example from  the  recent hepatocellular carcinoma (HCC)  study in  \cite{HCC14} which we obtained  through the  gene expression omnibus (GEO)  database. The GEO dataset consists of HCC tumor-infected ($T$) and healthy liver  ($N$)  tissue samples from three individuals denoted below  as follows  in relation to their original database designations   $T1=HCC448T, T2=HCC473T, T3=HCC510T$  and  $N1=HCC448N, N2=HCC473N, N3=HCC510N$.
For this dataset one of the questions  of research interest was whether the expression profiles of  genes associated with regulation of cell proliferation and programmed cell death  differ  across $T$ and $N$ samples as well as across individuals  (cf., e.g., \citealt{Kong:2013aa}). To address this specific  question, in contrast with the  previous TCR example,  we were  thus only interested in  a pre-selected subset of  the NGS counts. The final values of  $m=1332 $ and $n$ between 1.2 and 1.9 million reads \footnote{Based on these  values, the empirical versions of the conditions for the relevant  theorems in Section~3 were  considered satisfied.}  were obtained after aligning the   pre-selected NGS fragments   to the  HG19  reference genome with the Tophat2/Bowtie2 software \citep{Kim:2013aa} and performing the  transcript annotation with the  Ensembl genome browser (\url{www.ensembl.org}).  After  the final fragments-to-counts  conversion, our  data analysis was performed in three steps. 
First, the null hypothesis of the tissue homogeneity $H_0^{all}= \{T_1=N_1=T_2=N_2=T_3=N_3\}$ was tested (and rejected) based on the result of  Theorem~4 and the  corresponding asymptotic $p$-value obtained from the $\chi^2(3)$ distribution. 
Next,   the  hypothesis of the across-individuals homogeneity was tested by evaluating three pairwise null hypothesis $H_0^{ij}=\{{\cal D}_{1/2}(T_i,N_i)={\cal D}_{1/2}(T_j,N_j)]\},$ $1\le i<j\le 3$  (each  rejected) based on Theorem~4.  Finally, having rejected the homogeneity hypothesis  we have used the result of Theorem~2 to quantify the differences between the three sets of  $T$ and $N$ tissue samples. The details of the analysis are presented in Table~\ref{tab3}. As seen from the numerical results,   it seems that despite the large  individual differences between patients, the set of $m=1332$ genes   associated with cell proliferation and  death may be used to distinguish between  T-type  and N-type samples in HCC patients. 


\begin{table}
\begin{center}
\begin{tabular}{cccccc}
 Hypothesis  &Statistic   & $P$-value & ${\cal D}_{1/2}$ Value (CI)   \\
\hline
\tabT $H_0^{all}$&$\sum w_i[{\cal D}_{1/2}(i)-\mu_i]^2$  & $< 0.001$  & NA \\
\tabT  $H_0^{1,2}$  &${\cal D}_{1/2}(1)-{\cal D}_{1/2}(2)$& $< 0.01$ & ${\cal D}_{1/2}(1)$=0.553 (0.551, 0.555)  \\
\tabT  $H_0^{2,3}$   &${\cal D}_{1/2}(2)-{\cal D}_{1/2}(3)$ & &${\cal D}_{1/2}(2)$=0.292 (0.291, 0.294) \\
\tabT  $H_0^{3,1}$  & ${\cal D}_{1/2}(3)-{\cal D}_{1/2}(1)$ &  &${\cal D}_{1/2}(3)$= 0.346 (0.345 0.348)\\
\tabB\\
\hline
\end{tabular}
\caption{The $95\%$ confidence intervals for the pairwise symmetric Renyi Divergence ${\cal D}_{1/2}$ between the tumor and control (healthy) tissues from three individuals based on the profile of expression of pre-selected  $m=1332$ transcripts  related to cell proliferation. Here ${\cal D}_{1/2}(i)$ denotes  ${\cal D}_{1/2}(T_i,N_i).$}\label{tab3}
\end{center}
\end{table}

\section{Summary and Conclusions}

  We derived  two sets of  limit theorems for the Renyi entropy and divergence statistics. The first set of results  holds  for lineralizeable statistics (their first order Taylor approximations exist) whereas the second one holds in the degenerate case (when  the first order approximations  vanish)  and requires   analyzing the quadratic terms in the Taylor expansions.   Our  Renyi entropy limit theorems   complement  those obtained elsewhere for the Shannon entropy and divergence.  
  
 Based on the CLT results we have proposed here  a new framework for  analyzing molecular diversity of molecular (especially NGS) data  based on  the idea of  analyzing the   frequency/contingency tables  where cell counts are highly unbalanced (for instance, as arriving from  mixtures   of  heavy tailed, power-law type and uniform distributions) and the number of cells or, equivalently,  the counts distribution support size $m$, increases with the sample size $n$. For analyzing such tables,   we suggested using the  empirical Renyi entropy and divergence as  the statistical measures of, respectively,  diversity and pairwise similarity of different molecular sub-populations.    


 In  the two  examples of NGS analysis we have  shown how  the Renyi entropy  methods may be used for filtering out low frequency noise and for establishing valid confidence bounds in  pairwise divergence analysis for pre-selected transcripts. However,  it was also seen that  in order to apply our CLT results  the number of  transcripts had to be small relative to the sequencing depth.  For the special class of  heavy-tailed  power law distributions, our  results in particular  indicate that the appropriate entropy CLTs  are valid  (and thus so is our  proposed analysis framework)  when, roughly speaking,  $m/\sqrt{n}\to 0$ and not otherwise.  As such  restriction may be often limiting in very high diversity NGS data,   other  statistics  beyond  those  discussed here  and  not requiring such  condition  could be also  of interest.  We hope to  pursuing  this matter  further in our future 
work. 




 


 