A Good-Turing estimator for feature allocation models

TL;DR

This paper introduces a new Good-Turing based estimator for the missing mass in feature allocation models, providing theoretical guarantees, interpretability, and practical confidence intervals, applicable to complex biological data.

Contribution

It presents a robust, nonparametric estimator for the missing mass in feature models, with proven minimax optimality and connections to existing estimators, extending applicability to real data.

Findings

Estimator is simple, robust, and theoretically sound.

Provides non-asymptotic confidence intervals for the estimator.

Demonstrates effectiveness on synthetic and genomic data.

Abstract

Feature allocation models generalize species sampling models by allowing every observation to belong to more than one species, now called features. Under the popular Bernoulli product model for feature allocation, given $n$ samples, we study the problem of estimating the missing mass $M_{n}$, namely the expected number hitherto unseen features that would be observed if one additional individual was sampled. This is motivated by numerous applied problems where the sampling procedure is expensive, in terms of time and/or financial resources allocated, and further samples can be only motivated by the possibility of recording new unobserved features. We introduce a simple, robust and theoretically sound nonparametric estimator $\hat{M}_{n}$ of $M_{n}$. $\hat{M}_{n}$ turns out to have the same analytic form of the popular Good-Turing estimator of the missing mass in species sampling models, with the difference that the two estimators have different ranges. We show that $\hat{M}_{n}$ admits a natural interpretation both as a jackknife estimator and as a nonparametric empirical Bayes estimator, we give provable guarantees for the performance of $\hat{M}_{n}$ in terms of minimax rate optimality, and we provide with an interesting connection between $\hat{M}_{n}$ and the Good-Turing estimator for species sampling. Finally, we derive non-asymptotic confidence intervals for $\hat{M}_{n}$, which are easily computable and do not rely on any asymptotic approximation. Our approach is illustrated with synthetic data and SNP data from the ENCODE sequencing genome project.

Tables

Table 1. Table 1: Bias of M ^ n subscript ^ 𝑀 𝑛 \hat{M}_{n} : p j = 1 / j s subscript 𝑝 𝑗 1 superscript 𝑗 𝑠 p_{j}=1/j^{s} for j ≤ 10 5 𝑗 superscript 10 5 j\leq 10^{5} and p j = 0 subscript 𝑝 𝑗 0 p_{j}=0 for j > 10 5 𝑗 superscript 10 5 j>10^{5} , and different values of s 𝑠 s .

$s$	$n=10$	$n=50$	$n=100$	$n=1000$
0.6	1.310	0.459	0.241	0.001
0.8	0.268	0.075	0.042	0.003
1.0	0.100	0.020	0.010	0.001
1.2	0.052	0.008	0.004	0.000
1.4	0.033	0.004	0.002	0.000
1.6	0.023	0.002	0.001	0.000

Table 2. Table 2: B ​ i ​ a ​ s ​ ( M ^ n ) 2 / R n 𝐵 𝑖 𝑎 𝑠 superscript subscript ^ 𝑀 𝑛 2 subscript 𝑅 𝑛 Bias(\hat{M}_{n})^{2}/R_{n} *100: p j = 1 / j s subscript 𝑝 𝑗 1 superscript 𝑗 𝑠 p_{j}=1/j^{s} for j ≤ 10 5 𝑗 superscript 10 5 j\leq 10^{5} and p j = 0 subscript 𝑝 𝑗 0 p_{j}=0 for j > 10 5 𝑗 superscript 10 5 j>10^{5} , and different values of s 𝑠 s .

$s$	$n=10$	$n=50$	$n=100$	$n=1000$
0.6	16.06	13.05	9.57	0.89
0.8	2.99	1.44	1.08	0.19
1.0	1.35	0.34	0.19	0.03
1.2	0.94	0.16	0.07	0.01
1.4	0.76	0.10	0.04	0.00
1.6	0.64	0.07	0.03	0.00

Table 3. Table 3: Synthetic data: missing mass estimation with p j = 1 / j s subscript 𝑝 𝑗 1 superscript 𝑗 𝑠 p_{j}=1/j^{s} for j ≤ 10 5 𝑗 superscript 10 5 j\leq 10^{5} and p j = 0 subscript 𝑝 𝑗 0 p_{j}=0 for j > 10 5 𝑗 superscript 10 5 j>10^{5} , and different values of s 𝑠 s .

$s$	$n=50$			$n=250$
	$M_n$	${\hat{M}}_n$	$95\%$ CI	$M_n$	${\hat{M}}_n$	$95\%$ CI
0.6	183.81	184.66	(174.35, 198.12)	105.38	105.61	(101.74, 110.53)
0.8	33.79	33.67	(29.45, 39.51)	26.40	26.05	(24.40, 28.30)
1.0	7.02	7.02	(4.88, 9.91)	5.41	5.43	(4.66, 6.49)
1.2	1.92	1.87	(0.54, 3.60)	1.35	1.35	(0.92, 1.93)
1.4	0.71	0.72	(0, 1.98)	0.44	0.44	(0.15, 0.81)
1.6	0.34	0.36	(0, 1.39)	0.18	0.18	(0, 0.46)

Table 4. Table 4: Real data: missing mass estimation for the ENCODE sequencing genome project.

Population	$n=20$		$n=40$
	${\hat{M}}_n$	$95\%$ CI	${\hat{M}}_n$	$95\%$ CI
CEPH	55.6	(38.7, 88.9)	22.3	(15.6, 38.9)
CHB	50.0	(37.3, 81.3)	17.6	(12.6, 32.8)
JPT	61.7	(46.9, 93.4)	26.7	(20.9, 42.2)
YRI	88.3	(65.2, 125.2)	26.6	(18.9, 44.7)

Equations

$$\begin{array}{rl}{M}_n({\mathbf{Y}}_n,(p_j{)}_{j\ge 1})& =\mathbb{E}\left(\sum _{j\ge 1}\mathbb{I}\{X_{n,j}=0,Y_{n+1,j}=1\}|{\mathbf{Y}}_n\right)\\ & =\sum _{j\ge 1}{p}_j\mathbb{I}\{X_{n,j}=0\},\end{array}$$

$$K_{n,r}=\sum _{j\ge 1}\mathbb{I}\{X_{n,j}=r\}$$

$$\mathbb{E}(M_n)$$

$$=\frac{1}{\left(\genfrac{}{}{0ex}{}{n+1}{1}\right)}\sum _{j\ge 1}\left(\genfrac{}{}{0ex}{}{n+1}{1}\right)p_j(1-p_j{)}^{(n+1)-1}=\frac{1}{n+1}\mathbb{E}(K_{n+1,1})\simeq \frac{1}{n}\mathbb{E}(K_{n,1})$$

$${\hat{M}}_n=\frac{K_{n,1}}{n}$$

$${\tilde{M}}_n(\vartheta ,\alpha ,\beta )=\frac{\vartheta \Gamma (\alpha +\beta +n)}{\Gamma (\beta +n+1)}.$$

$${\hat{\vartheta }}_n=\frac{K_{n,1}\Gamma (\beta +n+1)}{\Gamma (\beta +\alpha +n)n}$$

$${\tilde{M}}_n({\hat{\vartheta }}_n,\alpha ,\beta )=\frac{K_{n,1}}{n}={\hat{M}}_n.$$

$$R_n({\hat{T}}_n;(p_j{)}_{j\ge 1})=\mathbb{E}[(M_n-{\hat{T}}_n{)}^2]=\mathrm{Bias}({\hat{T}}_n{)}^2+\mathrm{Var}({\hat{T}}_n-M_n),$$

$$R_n^{\ast }=\underset{{\hat{T}}_n}{\inf }\underset{(p_j{)}_{j\ge 1}\in {\mathcal{P}}_W}{\sup }{R}_n({\hat{T}}_n;(p_j{)}_{j\ge 1}),$$

$$\underset{(p_j{)}_{j\ge 1}\in {\mathcal{P}}_W}{\sup }{R}_n({\hat{M}}_n;(p_j{)}_{j\ge 1})\le \frac{1}{n^2}{W}^2+\frac{2n+1}{n(n+1)}W;$$

$$R_n^{\ast }\ge \frac{2W}{9(3n+1)}-\frac{14}{n^2}.$$

$$M_n({\mathbf{Y}}_n;P)=\sum _{r=0}^{W-1}{M}_n^{(r)}({\mathbf{Y}}_n^{(r)};(p_j^{(r)}{)}_{j\ge 1}),$$

$$\mathbb{P}\left({\hat{W}}_n-W>\sqrt{\frac{4W}{n}\log (1/\delta )}+\frac{1}{n}\log (1/\delta )\right)\le \delta$$

$$\mathbb{P}\left(-({\hat{W}}_n-W)>\sqrt{\frac{2W}{n}\log (1/\delta )}\right)\le \delta .$$

$$W\le {\left(\sqrt{{\hat{W}}_n+\frac{\log (1/\delta )}{2n}}+\sqrt{\frac{\log (1/\delta )}{2n}}\right)}^2.$$

$$\sqrt{W}\ge \sqrt{{\hat{W}}_n}-\sqrt{\frac{\log (1/\delta )}{n}}$$

$${\tilde{M}}_{\tilde{n}}^{\text{GT}}=\frac{K_{n,1}}{\tilde{n}}=\frac{K_{n,1}}{{\sum }_{j\ge 1}{X}_{n,j}}.$$

$${\hat{M}}_n=\hat{W}{\tilde{M}}_{\tilde{n}}^{\text{GT}}=\frac{{\sum }_{j\ge 1}{X}_{n,j}}{n}\frac{K_{n,1}}{{\sum }_{j\ge 1}{X}_{n,j}}=\frac{K_{n,1}}{n},$$

$${\hat{M}}_n-L_{\delta }(K_{n,1},K_{n,2})\le M_n\le {\hat{M}}_n+U_{\delta }(K_n),$$

$$L_{\delta }(K_{n,1},K_{n,2})=\frac{2c_{\delta }(K_{n,2})}{n(n-1)}+\frac{\log (6/\delta )}{n}+\sqrt{2\log \left(\frac{6}{\delta }\right)\left(\frac{4c_{\delta }(K_{n,1})}{n(n-1)}+\frac{2c_{\delta }(K_{n,2})}{n^2}\right)}$$

$$U_{\delta }(K_n)=\frac{\log (6/\delta )}{n-1}+\sqrt{2\log \left(\frac{6}{\delta }\right)\frac{4c_{\delta }(K_n)}{(n-1{)}^2(1-2/n)}}.$$

$$n^{\ast }=\inf \{n\ge 0:[h(K_n+{\hat{M}}_n)-h(K_n)]\le c\}.$$

$${\tilde{M}}_n(\vartheta ,\alpha ,\beta )=E[\tilde{\mu }((0,1)\setminus \left\{F_1^{\ast },\dots ,F_{K_n}^{\ast }\right\})|{\mathbf{Y}}_n].$$

$$\tilde{\mu }|{\mathbf{Y}}_n\stackrel{d}{=}{\tilde{\mu }}_n+\sum _{\ell =1}^{K_n}{J}_{\ell }{\delta }_{F_{\ell }^{\ast }}$$

$${\tilde{M}}_n(\vartheta ,\alpha ,\beta )$$

$$\displaystyle=\mathbb{E}[\tilde{\mu}_{n}((0,1))]=-\frac{{\rm d}}{{\rm d}u}\mathbb{E}[e^{-u\tilde{\mu}_{n}((0,1))}]\Big{|}_{u=0}=\vartheta\frac{\Gamma(\alpha+\beta+n)}{\Gamma(\beta+n+1)}$$

$$\overline{\rho }(s):={\int }_x^1\rho (s)\mathrm{d}s=\frac{\vartheta }{\Gamma (1-\alpha )}{\int }_x^1{s}^{-1-\alpha }(1-s{)}^{\beta +\alpha -1}\mathrm{d}s$$

$$\underset{x\to 0}{\lim }\frac{\overline{\rho }(s)}{x^{-\alpha }}=\frac{\vartheta }{\Gamma (1-\alpha )}\underset{x\to 0}{\lim }\frac{x^{-1-\alpha }(1-x{)}^{\beta +\alpha -1}}{\alpha x^{-\alpha -1}}=\frac{\vartheta }{\Gamma (1-\alpha )\alpha }.$$

$$\mathbb{E}(M_{n,1})=\sum _{j\ge 1}{p}_j\mathbb{P}(X_{n,j}=1)=\sum _{j\ge 1}\left(\genfrac{}{}{0ex}{}{n}{1}\right)p_j^2(1-p_j{)}^{n-1}$$

$$\mathbb{E}({\hat{M}}_n-M_n)$$

Full text

Keywords: Bernoulli product model, Feature allocation model, Good-Turing estimator, Minimax rate optimality, Missing mass, Non-asymptotic uncertainty quantification, Nonparametric empirical Bayes, SNP data

1 Introduction

Feature allocation models generalize classical species sampling models by allowing every observation to belong to more than one species, now called features. In particular, every observation is endowed with a (unknown) finite set of features selected from a (possibly infinite) collection of features $(F_{j})_{j\geq 1}$. We conveniently represent each observation as a binary sequence, with entries indicating the presence (1) or absence (0) of each feature. Every feature $F_{j}$ is associated to an unknown probability $p_{j}$, and each observation displays feature $F_{j}$ with probability $p_{j}$, for $j\geq 1$. The Bernoulli product model, or binary independence model, is arguably the most popular feature allocation model. It models the $i$–th observation as a sequence $Y_{i}=(Y_{i,j})_{j\geq 1}$ of independent Bernoulli random variables with unknown success probabilities $(p_{j})_{j\geq 1}$, with the assumption that $Y_{r}$ is independent of $Y_{s}$ for any $r\neq s$. Bernoulli product models have found applications in several scientific disciplines. They have been applied extensively in ecology when animals are captured using traps and each observation is an incidence vector collecting the presence or absence of each species in the traps (e.g., Colwell et al. (2012), Chao et al. (2014) and Chao and Colwell (2017)). Besides ecology, Bernoulli product models found applications in the broad area of biosciences (e.g., Ionita-Laza et al. (2009), Gravel (2014) and Zou et al. (2016)); in the analysis of choice behaviour arising from psychology and marketing (Görür et al. (2006)); in binary matrix factorization for dyadic data (Meeds et al. (2007)); in graphical models (e.g., Wood et al. (2006) and Wood and Griffiths (2007)); in the analysis of similarity judgement matrices (Navarro and Griffiths (2007)); in network data analysis (Miller et al. (2006)).

Let $\mathbf{Y}_{n}=(Y_{1},\ldots,Y_{n})$ denote a collection of $n$ random samples collected under the Bernoulli product model with unknown feature probabilities $(p_{j})_{j\geq 1}$. Furthermore, let $X_{n,j}=\sum_{1\leq i\leq n}Y_{i,j}$ be the number of times that feature $F_{j}$ has been observed in $\mathbf{Y}_{n}$. That is, $X_{n,j}$ is a Binomial random variable with parameters $(n,p_{j})$ for any $j\geq 1$. In this paper we consider the problem of estimating the conditional expected number, given the random sample $\mathbf{Y}_{n}$, of hitherto unseen features that would be observed if one additional sample $Y_{n+1}$ was collected, i.e.

[TABLE]

where $\mathbb{I}{}$ is the indicator function. For easiness of notation, in the rest of the paper we will not highlight the dependence on $\mathbf{Y}_{n}$ and $(p_{j})_{j\geq 1}$, and simply write $M_{n}$ instead of $M_{n}(\mathbf{Y}_{n},(p_{j})_{j\geq 1})$. The statistic $M_{n}$ is referred to as the missing mass, namely the sum of the probability masses of unobserved features in the first $n$ samples.

Interest in estimating the missing mass is motivated by numerous applied problems where the sampling procedure is expensive, in terms of time and/or financial resources, and further draws can be only motivated by the possibility of recording new unobserved features. In genetics, for instance, the ambitious prospect of growing databases to encompass hundreds of thousands of human genomes, makes important to quantify the power of large sequencing projects to discover new genetic variants (Auton et al. (2015)). An accurate estimate of $M_{n}$ will enable better study design and quantitative evaluation of the potential and limitations of these datasets. Indeed one can fix a suitable threshold such that the sequencing procedure takes place until the estimate of $M_{n}$ becomes for the first time smaller than the threshold. This introduces a criterion for evaluating the effectiveness of further sampling, providing a roadmap for large-scale sequencing projects. A Bayesian parametric estimator of (1) has been introduced in Ionita-Laza et al. (2009), and it relies on a Beta prior distribution for the unknown probabilities $p_{j}$’s. A limitation of this approach is that it requires parametric forms for the distribution of variant frequencies, which requires some model of demography and selection. For instance, the Beta prior distribution is a reasonable assumption for neutrally evolving variants but may not be appropriate for deleterious mutations. To overcome this drawback, a nonparametric approach to estimate $M_{n}$ has been proposed in the recent work of Zou et al. (2016). This is a purely algorithmic type approach, whose core is a linear programming based algorithm.

In this paper, we consider an estimator $\hat{M}_{n}$ of $M_{n}$ which has the same analytic form of the popular Good-Turing estimator of the missing mass under the multinomial model for species sampling (e.g., Good (1953) and Robbins (1968)), with the difference that the two estimators have different ranges (supports). We refer to $\hat{M}_{n}$ as the Good-Turing estimator for feature allocation models. As in Zou et al. (2016), our approach is nonparametric in the sense that it does not rely on any distributional assumption on the unknown $p_{j}$’s. The use of the Good-Turing estimator for estimating $M_{n}$ was first considered in Chao and Colwell (2017), where also a bootstrap procedure is presented to approximate the variance of the estimator. Our work delves into the Good-Turing estimator for feature allocation models, thus providing theoretical guarantees for its use. We show that $\hat{M}_{n}$ admits a natural interpretation both as a jackknife (resampling) estimator (Quenouille (1956) and Tukey (1958)) and as a nonparametric empirical Bayes estimator in the sense of Efron and Morris (1973). Theoretical properties of the estimator $\hat{M}_{n}$ are investigated. Specifically, we first provide a lower bound for the minimax risk under a squared loss function and then we show that the mean squared error of the proposed estimator achieves the optimal minimax rate for the estimation of $M_{n}$. This asymptotic analysis provides with an interesting connection between the Good-Turing estimators for species sampling models and for feature allocation models in terms of the limiting Poisson regime of the Binomial distribution. Finally, we derive a non-asymptotic confidence interval for $\hat{M}_{n}$ by means of novel Bernstein type concentration inequalities which are of separate interest; the confidence interval is easy to compute and it does not rely on any asymptotic assumption or any parametric constraint on the unknown $p_{j}$’s. We illustrate the proposed methodology by the analysis of various synthetic data and SNP datasets from the ENCODE (http://www.hapmap.org/downloads/encode1.html.en) sequencing genome project.

The paper is structured as follows. In Section 2 we introduce the Good-Turing estimator $\hat{M}_{n}$ for the missing mass $M_{n}$, we present its interpretation as a jackknife estimator and a nonparametric empirical Bayes estimator, and we give provable guarantees for its performance. In Subsection 3 we derive a non-asymptotic confidence interval for $\hat{M}_{n}$, in Subsection 4 we discuss the use of $\hat{M}_{n}$ for designing cost-effective feature inventories, and in Section 5 we present a numerical illustration of $\hat{M}_{n}$ with synthetic and real data. Some concluding remarks on future works are discussed in Section 6. Proofs are deferred to Appendix A.

2 A Good-Turing estimator for $M_{n}$

Let $\mathbf{Y}_{n}=(Y_{1},\ldots,Y_{n})$ be a collection of $n$ random samples under the Bernoulli product model with unknown feature probabilities $(p_{j})_{j\geq 1}$ such that $\sum_{j\geq 1}p_{j}<+\infty$. That is, $Y_{i}=(Y_{i,j})_{j\geq 1}$ is a sequence of independent Bernoulli random variables, with $p_{j}$ being the success probability of $Y_{i,j}$ for any $i=1,\ldots,n$, and $Y_{r}$ is independent of $Y_{s}$ for any $r\neq s$. Note that the assumption $\sum_{j\geq 1}p_{j}<+\infty$ implies that each $Y_{i}$ displays finitely many features almost surely; indeed, by monotone convergence, we have $\sum_{j\geq 1}p_{j}<+\infty$ if and only if $\mathbb{E}(\sum_{j\geq 1}Y_{i,j})<+\infty$, which in turns implies that $\sum_{j\geq 1}Y_{i,j}<+\infty$ almost surely. If $X_{n,j}$ denotes the number of times that feature $F_{j}$ has been observed in the random sample $\mathbf{Y}_{n}$, then

[TABLE]

is the number of features with frequency $r$ in $\mathbf{Y}_{n}$. Let $K_{n}$ be the total number of features in $\mathbf{Y}_{n}$, i.e. $K_{n}=\sum_{r\geq 1}K_{n,r}$. For any two sequences $(a_{n})_{n\geq 1}$ and $(b_{n})_{n\geq 1}$, write $a_{n}\simeq b_{n}$ to mean that $a_{n}/b_{n}\rightarrow 1$ as $n\rightarrow+\infty$. An intuitive estimator of $M_{n}$ can be deduced from a comparison between expectations of $M_{n}$ and $K_{n,1}$. Specifically, since $X_{n,j}$ is a Binomial random variable with parameter $(n,p_{j})$,

[TABLE]

and set

[TABLE]

as an estimator of $M_{n}$, for any $n\geq 1$. The estimator $\hat{M}_{n}$ is nonparametric, in the sense that it does not rely on any distributional assumptions on the feature probabilities $(p_{j})_{j\geq 1}$. See Chao and Colwell (2017) for a somehow related derivation of $\hat{M}_{n}$.

The estimator $\hat{M}_{n}$ turns out to have the same analytic form of the popular Good-Turing estimator of the missing mass for species sampling models. The Good-Turing estimator first appeared in Good (1953) as a nonparametric empirical Bayes estimator under the classical multinomial model for species sampling, i.e., $(Y_{1},\ldots,Y_{n})$ are $n$ random samples from a population of individuals belonging to a (possibly infinite) collection of species $(S_{j})_{j\geq 1}$ with unknown proportions $(p_{j})_{j\geq 1}$ such that $\sum_{j\geq 1}p_{j}=1$. Under the multinomial model every observation is endowed with one species selected from $(S_{j})_{j\geq 1}$, and hence $K_{n,1}\in\{0,1,\ldots,n\}$. Therefore, although the estimator $\hat{M}_{n}$ has the same analytic form of the Good-Turing estimator, the two estimators have different ranges (supports): while the Good-Turing estimator for species sampling models takes values in $[0,1]$, the Good-Turing estimator for feature allocation models takes positive (finite) values.

Hereafter we show that: i) the estimator $\hat{M}_{n}$ admits natural interpretations both as a jackknife (resampling) estimator and as a nonparametric empirical Bayes estimator in the sense of Efron and Morris (1973); ii) the mean squared error of $\hat{M}_{n}$ converges to zero at the best possible rate for the estimation of the missing mass $M_{n}$; iii) the estimator $\hat{M}_{n}$ is linked to the Good-Turing estimator for species sampling models in terms of the limiting Poisson regime of the Binomial distribution.

2.1 Interpretations of $\hat{M}_{n}$

We first show that the estimator $\hat{M}_{n}$ admits a natural interpretation as a jackknife estimator in the sense of Quenouille (1956) and Tukey (1958). See also the monograph by Efron (1987) and references therein for a comprehensive account on jackknife (resampling) estimators. Let $K_{n}(j)$ denote the number of distinct features observed in the sample $\mathbf{Y}_{n}$ after the removal of the $j$-th sample $Y_{j}$. It is easy to show that the missing mass $M_{n-1}$ equals the posterior expected value of $K_{n}-K_{n-1}(j)$, given all the samples but $Y_{j}$, i.e., $M_{n-1}=\mathbb{E}[K_{n}-K_{n-1}(j)|\left\{Y_{1},\ldots,Y_{n}\right\}\setminus Y_{j}]$. Accordingly, the difference $K_{n}-K_{n-1}(j)$ provides with an unbiased estimator of $M_{n-1}$. The jackknife estimator is then $\sum_{1\leq j\leq n}(K_{n}-K_{n-1}(j))/n$ which equals $\hat{M}_{n-1}=K_{n,1}/n$; indeed, for any fixed $j$, $K_{n}-K_{n-1}(j)$ coincides with the number of features displayed only by the $j$-th sample $Y_{j}$.

The estimator $\hat{M}_{n}$ also admits a natural interpretation as a nonparametric empirical Bayes estimator in the sense of Efron and Morris (1973). A Bayesian nonparametric approach to estimate $M_{n}$ relies on the specification of a prior distribution for the unknown feature probabilities $(p_{j})_{j\geq 1}$. Here we consider the three parameters Beta process prior (e.g., Teh and Görür (2009) and James (2017)). This is a well-known generalization of the celebrated Beta process prior of Hjort (1990), and it is defined as the distribution of a completely random measure (Daley and Vere-Jones (2008)) with intensity measure $\nu({\rm d}s,{\rm d}f)=\rho(s){\rm d}s\mathbb{I}{(0,1)}(f){\rm d}f$ where $\rho(s)=\vartheta s^{-1-\alpha}(1-s)^{\beta+\alpha-1}/\Gamma(1-\alpha)$, for $\alpha\in(0,1)$, $\beta>-\alpha$ and $\vartheta>0$, with $\Gamma(\cdot)$ being the Gamma function. In particular. under the three parameters Beta process prior, the Bayesian point estimator of $M_{n}$ under squared loss function is

[TABLE]

We refer to the Appendix A.1 for further details on the three parameters Beta process prior, and on the derivation of (2). In particular, it is shown in Appendix A.1 that

[TABLE]

is a consistent estimator of $\vartheta$ for any $\alpha\in(0,1)$ and $\beta>-\alpha$, and plugging $\hat{\vartheta}_{n}$ into (2) we obtain

[TABLE]

In other terms, the proposed estimator $\hat{M}_{n}$ coincides with the nonparametric empirical Bayes estimator $\tilde{M}_{n}(\hat{\vartheta}_{n},\alpha,\beta)$ of the missing mass. $\tilde{M}_{n}(\hat{\vartheta}_{n},\alpha,\beta)$ is obtained by assigning the three parameters Beta process prior to $(p_{j})_{j\geq 1}$ and then estimating its mass parameter $\vartheta$ with an appropriately chosen consistent estimator.

2.2 Optimality of $\hat{M}_{n}$

We start by showing that the mean squared error ($L_{2}$-risk) of the estimator $\hat{M}_{n}$ goes to zero at the best possible rate. This result legitimates the use of $\hat{M}_{n}$ as an estimator of the missing mass $M_{n}$. For any $0<W<+\infty$ let $\mathscr{P}_{W}=\{(p_{j})_{j\geq 1}:p_{j}\in(0,1)\text{ and }\,\sum_{j\geq 1}p_{j}\leq W\}$ be a class of features probabilities. Also, let $\hat{T}_{n}$ denote a whichever estimator of the missing mass $M_{n}$ based on the random sample $\mathbf{Y}_{n}$, that is a map on the $n$-fold product space of observations $\{0,1\}^{\infty}\times\cdots\times\{0,1\}^{\infty}$ and taking values in $\mathds{R}^{+}$. For a specific collection of feature probabilities $(p_{j})_{j\geq 1}$ in the class $\mathscr{P}_{W}$, the $L_{2}$-risk of the estimator $\hat{T}_{n}$ is defined as follows

[TABLE]

and the minimax risk is

[TABLE]

i.e. the infimum, with respect to the set of possible estimators $\hat{T}_{n}$ of the missing mass $M_{n}$, of the worst-case risk over the class of feature probabilities $\mathscr{P}_{W}$. The next theorem provides with an upper bound for $R_{n}(\hat{M}_{n};(p_{j})_{j\geq 1})$ and a lower bound for $R^{\ast}_{n}$.

Theorem 2.1

Let $(p_{j})_{j\geq 1}$ be feature probabilities in the class $\mathscr{P}_{W}$. Then,

• i)

for any $n\geq 1$

[TABLE]

• ii)

for any $n\geq 2$

[TABLE]

See Appendix A.2 for the proof of Theorem 2.1. For any two sequences $(a_{n})_{n\geq 1}$ and $(b_{n})_{n\geq 1}$, write $a_{n}\lesssim b_{n}$ to mean that there exists a constant $C>0$ such that $a_{n}\leq Cb_{n}$ for all $n$ and $a_{n}\asymp b_{n}$ to mean that we have both $a_{n}\lesssim b_{n}$ and $b_{n}\lesssim a_{n}$. Part i) of Theorem 2.1 shows that $\sup_{(p_{j})_{j\geq 1}\in\mathscr{P}_{W}}R_{n}(\hat{M}_{n};(p_{j})_{j\geq 1})\lesssim W/n$, namely as $n$ goes to infinity the mean squared error of $\hat{M}_{n}$ goes to zero at rate $W/n$. Part ii) of Theorem 2.1 provides with a lower bound for the minimax risk, showing that $R^{\ast}\gtrsim W/n$. In other terms, this shows that the mean squared error of $\hat{M}_{n}$ goes to zero at the best possible rate. Theorem 2.1 then shows that the estimator $\hat{M}_{n}$ is essentially rate optimal due to the matching minimax lower bound in the class $\mathscr{P}_{W}$ of admissible probabilities’ masses up to a constant factor.

In Theorem 2.1 we considered the class $\mathscr{P}_{W}$ of feature probabilities having total mass bounded by $W$. We highlight the crucial role or this assumption in order to provide asymptotic results. First, notice from Theorem 2.1 that, for a fixed value of the sample size $n$, the minimax rate increases linearly in $W$. This implies that, for a fixed $n$, the estimation problem can be made as hard as desired if no bounds are imposed on the sum of the possible vectors of probabilities $(p_{j})_{j\geq 1}$. Since at first glance this result can seem slightly counter intuitive, let us consider an illustrative example which should clarify why the estimation difficulty of the problem is proportional to $W$. To simplify the argument, suppose that $W$ is a strictly positive integer and consider $W$ frequency vectors $(p^{(0)}_{j})_{j\geq 1},\ldots,(p^{(W-1)}_{j})_{j\geq 1}\in\mathscr{P}_{1}$, i.e. such that $\sum_{j}p^{(r)}_{j}\leq 1$ for all $r\in\{0,\ldots,W-1\}$. Now, let us construct another vector of frequencies $P=(p_{j})_{j\geq 1}\in\mathscr{P}_{W}$ obtained by setting $p_{qW+r}=p^{(r)}_{q}$ for $q\geq 0$ and $0\leq r\leq W-1$. The resulting vector $P$ is the concatenation of $(p^{(0)}_{j})_{j\geq 1},\ldots,(p^{(W-1)}_{j})_{j\geq 1}$, such that in the first $W$ entries of $P$ we put the first elements $(p^{(0)}_{1},\ldots,p^{(W-1)}_{1})$, followed by the second ones $(p^{(0)}_{2},\ldots,p^{(W-1)}_{2})$, the third ones $(p^{(0)}_{3},\ldots,p^{(W-1)}_{3})$ and so on. Furthermore, any observation $(Y_{i,j})_{j\geq 1}$ sampled from the Bernoulli product model with frequencies given by this $P$ can be rewritten, following a similar construction, as the concatenation of $W$ observations $Y^{(r)}_{i}=(Y^{(r)}_{i,j})_{j\geq 1}$, $r\in\{0,\ldots,W-1\}$, each of them sampled from a binomial product model with the corresponding frequencies $(p^{(r)}_{j})_{j\geq 1}$. Hence, the missing mass for a random sample from $P$ can be related to the missing masses of each of its $W$ subcomponents, by

[TABLE]

where $\mathbf{Y}^{(r)}=(Y_{1}^{(r)},\ldots,Y_{n}^{(r)})$. From this construction, we can see that by trying to estimate the missing mass on the left hand side, we are basically trying to estimate a sum of $W$ unrelated quantities, which explains why the error is linear in $W$.

In order to apply Theorem 2.1 and evaluate the performances of the Good-Turing estimator compared to the minimax rate, we would need to know an upper bound $W$, i.e. an upper bound on the total mass of the unknown vector of probabilities $(p_{j})_{j\geq 1}$. In real life applications, $W$ is unlikely to be known. However, we can easily estimate it. Specifically, since the total mass $W$ is the expected number of features displayed per observation, we can use as a consistent estimator for $W$ the estimator $\hat{W}_{n}:=\sum_{j\geq 1}X_{n,j}/n$. Besides the following concentration inequality for $\hat{W}_{n}$ around its mean $W$ may be useful to measure the probability that $W_{n}$ deviates from $W$. In particular it can be proved that for any $\delta\in(0,1)$

[TABLE]

and

[TABLE]

The proofs of these exponential tail bounds rely on suitable bounds on the corresponding log-Laplace transform of $\hat{W}_{n}-\mathbb{E}\hat{W}_{n}$. See Section A.3 for details. Equation (7) leads to an upper bound on $W$, indeed with probability $1-\delta$ one has

[TABLE]

Analogously, one may directly apply Equation (6) in order to find a lower bound for $W$. More precisely, it is straightforward to see the validity of the following inequality

[TABLE]

with probability bigger than $1-\delta$. Therefore (8)–(9) result in a confidence intervals for $W$.

Finally, we conclude this subsection by remarking that the proposed estimator displays some bias. The formula of the theoretical bias can be found in the proof of Theorem 2.1. However, the effect of this bias becomes negligible for any reasonable sample size. As a simple illustration, in Table 1, we display the theoretical bias of the estimator for different sample sizes, when the frequencies $(p_{j})_{j\geq 1}$ are chosen according to the Zipf’s Law for different values of the parameter $s$ (the same simulation setting will be reconsidered and better described in Section 5). In Table 2, we report for each combination of sample size $n$ and parameter $s$ the percentage of the risk that is due to the bias. As can be seen from Table 1 and 2, expect for very small sample size and Zipf parameter, under which most of the frequencies are high, the bias is always almost negligible and does not contribute much to the risk of estimator, which is mainly due to its variance.

2.3 Connection to the Good Turing estimator for species sampling models

We relate the Good-Turing estimator for feature allocation models with the classical Good-Turing estimator for species sampling. This link relies on the well-known limiting Poisson approximation of Binomial random variables. Theorem 2.1 states that, in the feature allocation models, the Good-Turing estimator achieves a risk of order $R_{n}\asymp\frac{W}{n}$, while it is known from Rajaraman et al. (2017) that the risk $\tilde{R}_{n}$ of the Good Turing estimator in the species sampling case asymptotically behaves as $\tilde{R}_{n}\asymp 1/n$. In order to compare the two models, we will consider the limiting scenario when $W\rightarrow 0$. Let us consider a vector of feature frequencies $(p_{j})_{j\geq 1}$, such that $\sum_{j}p_{j}=W=\frac{\lambda}{n}$ for some positive value $\lambda$ and denote $\tilde{p}_{j}=\frac{p_{j}}{\sum_{j\geq 1}p_{j}}$ the normalized probability vector. Applying the large $n$ Poisson approximation of the binomial distribution, it follows that each $X_{n,j}$ is now approximately distributed according to a Poisson distribution with mean $\lambda\tilde{p}_{j}$. Therefore, the approximated model for large $n$ boils down to sample first an effective size $\tilde{n}=\sum_{j\geq 1}X_{n,j}$, distributed according to a Poisson with parameter $\lambda$, and, conditionally on it, sample $\tilde{n}$ independent and identically distributed observations $(\tilde{Y}_{1},\ldots,\tilde{Y}_{\tilde{n}})$ from probability vector $(\tilde{p}_{j})_{j\geq 1}$. Hence it is equivalent to a species sampling model where the sample size $\tilde{n}$ is assumed to be random. Denote $\tilde{M}_{\tilde{n}}$ the missing mass in the corresponding species sampling model. Now, by noticing that $M_{n}=\frac{\lambda}{n}\tilde{M}_{\tilde{n}}$, it follows that $R_{n}\approx\frac{\lambda^{2}}{n^{2}}\mathbb{E}_{\tilde{n}}(\tilde{R}_{\tilde{n}})$. Therefore, we expect the quantity $\frac{n^{2}}{\lambda^{2}}R_{n}$ to have the same asymptotic behaviour of $\mathbb{E}_{\tilde{n}}\left(\frac{1}{\tilde{n}}\mathbb{I}{\{\tilde{n}>0\}}\right)$. Indeed, from Theorem 2.1 and $W=\frac{\lambda}{n}$, we get that $\frac{n^{2}}{\lambda^{2}}R_{n}\asymp\frac{1}{\lambda}=\frac{1}{\mathbb{E}(\tilde{n})}$.

We conclude by remarking that the above construction also provides a justification for the Good-Turing estimator in the feature allocation models. Indeed, in the context of feature allocation models, we want to estimate $M_{n}=W\ \tilde{M}_{\tilde{n}}$, where both $W$ and $\tilde{M}_{\tilde{n}}$ are unknown and need to be estimated. In order to estimate $\tilde{M}_{\tilde{n}}$, we can use the Good-Turing estimator for species models, which here turns to be

[TABLE]

However, we also need to estimate $W=\sum_{j}p_{j}$ in order to make use of $\tilde{M}^{GT}_{\tilde{n}}$ as an estimator of $M_{n}$. As pointed out at the end of the previous subsection, a consistent estimator of $W$ is $\hat{W}=\frac{\sum_{j}X_{n,j}}{n}$. Then, by combining the estimator $\hat{W}$ for $W$ with the estimator $\tilde{M}_{\tilde{n}}^{\text{GT}}$ for $\tilde{M}_{\tilde{n}}$, we obtain the following estimato of $M_{n}$,

[TABLE]

which turns out to be exactly the Good-Turing estimator for feature allocation models.

3 A confidence interval for $\hat{M}_{n}$

In this section, we consider the problem of uncertainty quantification of the proposed estimator of the missing mass. In particular, we exploit tools from concentration inequalities for sum of independent random variables (Boucheron et al. (2013)) to introduce a non-asymptotic level-$\delta$ confidence interval for $\hat{M}_{n}$, for any $\delta\in(0,1)$.

Theorem 3.1

Let $(p_{j})_{j\geq 1}$ be any sequence of feature probabilities s.t. $\sum_{j\geq 1}p_{j}<+\infty$, and set $c_{\delta}(x)=(\sqrt{\log(1/\delta)/2}+\sqrt{7\log(1/\delta)/6+x})^{2}$ for any nonnegative integer $x$ and $\delta\in(0,1)$. Then, with probability at least $1-\delta$

[TABLE]

where

[TABLE]

and

[TABLE]

This theorem can be proven in two steps. We first derive bounds on the log-Laplace transform of $D_{n}=\hat{M}_{n}-M_{n}$, from which we can obtain a confidence interval for $D_{n}$ as a function of $\mathbb{E}K_{n,1}$, $\mathbb{E}K_{n,2}$ and $\mathbb{E}K_{n}$. Then, by deriving bounds on the log-Laplace of $K_{n,1}$, $K_{n,2}$ and $K_{n}$, we are able to bound with high probability the deviation of these random variables from their mean. Theorem 3.1 then follows by plugging-in the results of step 2 in the bounds of step 1. The details are left in Appendix A.4 . For any $\delta\in(0,1)$, Theorem 3.1 provides with a level $1-\delta$ confidence interval for $\hat{M}_{n}$ that holds for every sample size $n$ and for every possible values of the feature probabilities $(p_{j})_{j\geq 1}$. Indeed, this is derived by applying finite sample concentration inequalities, without using any asymptotic approximation, and it does not rely on any distributional assumption on the feature probabilities $(p_{j})_{j\geq 1}$. Note that the proposed confidence interval can be easily evaluated without resorting to any simulation based strategy. It is enough to count the total number of features and number of features with frequency one and two in the observable sample, i.e. $K_{n}$, $K_{n,1}$ and $K_{n,2}$, and plug these quantities into $U_{\delta}(K_{n})$ and $L_{\delta}(K_{n,1},K_{n,2})$ to be added and subtracted to $\hat{M}_{n}$.

4 A Stopping Rule for the Discovery Process

As we recalled in the Introduction, interest in estimating missing mass $M_{n}$ is motivated by the design of feature inventories that are cost-effective in terms of the number of features discovered and the amount of resources allocated in the discovery process. This is a fundamental issue in numerous scientific disciplines where the sampling procedure is expensive, in terms of time and/or financial resources. Feature inventories must then be designed in such a way to redirect the search of new features to more productive sites, methods or time periods whenever the sampling effort becomes unprofitable. In such a context, the estimator $\hat{M}_{n}$ is the key ingredient for defining an adaptive sequential approach in terms a stopping rule for the discovery process. Specifically, let $h$ be an utility function, defined on the integers, such that $h(k)$ is the gain of observing $k$ features; assume that $h$ in non-decreasing and concave. If $c$ is the cost associated with each sampling step, and $K_{n}$ is the number of features in the sample $\mathbf{Y}_{n}$, then the stopping rule may be defined as

[TABLE]

This brief discussion highlights how to exploit the estimator $\hat{M}_{n}$ within the context of designing cost-effective feature inventories. In particular, it gives rises to the challenging problem of finding the time $n^{\ast}$ at which it is optimal to conclude the discovery process, i.e. the time that maximizes the expected payoff $\mathbb{E}(h(K_{n})-cn)$.

5 Numerical Illustration

In this section we illustrate the experimental performance of the estimator $\hat{M}_{n}$ by the analysis of various synthetic data, and by the analysis SNP datasets from a genome project. In particular, let $N$ denote the total number of possible features in the whole population of individuals. With regards to synthetic data, we present a numerical illustration by setting $p_{j}=1/j^{s}$ for $j\leq N$ and $p_{j}=0$ for $j>N$, with $s$ being a nonnegative parameter. Note that the feature probability masses $p_{j}$’s correspond to the unnormalized masses of the ubiquitous Zipf distribution. Recall that the parameter $s$ controls how the total mass is spread among the features: the lager $s$, the smaller is the number of features with high probability. In other terms the larger $s$, the larger is the number of features with small frequencies (rare features). Hereafter we set $N=10^{5}$, and we consider the following values for the discount parameter $s$: 0.6, 0.8, 1.0, 1.2, 1.4, 1.6. For each of these values of the parameter $s$, we take $100$ samples of size $n=50,250,1000$ from the population $(p_{j})_{1\leq j\leq N}$. Table 1 displays the true value of the missing mass $M_{n}$, the estimated value $\hat{M}_{n}$ and the corresponding $95\%$ confidence interval (CI). All the experiments are averaged over the 100 samples. Results in Table 1 show that $\hat{M}_{n}$ provides good estimates of true value of $M_{n}$ in all the scenarios that we considered. In addition, confidence intervals are quite narrow around the estimator and contain always the true value of the missing mass.

We conclude with an application to SNP data. In this context, each sample represents a genome sequence and we are interested in variations with respect to the reference genome. Single nucleotide polymorphisms (SNPs) are the most common type of genetic variation among people. Each SNP represents a difference in a single DNA building block, which is called a nucleotide. For example, at a specific position in the human genome, the C (Cytosine) nucleotide may appear in most individuals, but the position can be replaced by an A (Adenine) in a small group of people. This means that there is a SNP at that position of the DNA. SNPs are important biological markers, helping scientists locate genes that are associated with diseases. We use the SNP datasets from the ENCODE sequencing genome project (http://www.hapmap. org/downloads/encode1.html.en). The same project was analyzed in Ionita-Laza et al. (2009). Ten 500-kb regions of the genome were sequenced in 209 unrelated DNA samples: 60 Yoruba (YRI), 60 CEPH European (CEPH), 45 Han Chinese (CHB), and 44 Japanese (JPT). These regions were chosen to be representative of the genome in general, including various chromosomes, recombination rates, gene density, and values of nontranscribed conservation with mouse. To make results comparable across the 4 populations (YRI, CEPH, CHB, and JPT), we consider only $n=20$ of the sequenced individuals for each dataset. Table 2 displays the estimated values and $95\%$ confidence interval (CI). For samples of 20 individuals, the YRI population displays the highest estimate of the missing mass. This agrees with results in Ionita-Laza et al. (2009), showing that the African population is the most diverse. We also consider increasing the sample size to $n=40$ of the sequenced individuals for each dataset. Table 2 shows how the missing mass decreases with respect to the case $n=20$. This highlights the saturation effect in discovering new variants. The discovery process is very efficient in the beginning, but after many individuals are sequenced, each additional individual contributes less and less to the pool of the newly discovered variants.

6 Concluding Remarks

In this paper we delved into the Good-Turing estimator for feature allocation models, thus providing theoretical guarantees for its use. In particular, we proved that the mean squared error of $\hat{M}_{n}$ goes to zero at the best possible rate as the sample size $n$ goes to infinity. Our results are simple, intuitive and easily implementable from practitioners. It distinguishes from the approaches of Ionita-Laza et al. (2009) and Zou et al. (2016) for being the first nonparametric statistical approach for estimating the missing mass in feature allocation models. In particular, differently from Ionita-Laza et al. (2009) and Zou et al. (2016), we associated to the proposed estimator an exact (non-asymptotic) quantification of its uncertainty. Results and techniques presente in this paper paves the way to new research directions in the context of feature allocation problems. In particular, three promising directions are: i) estimating the conditional expected number, given an observable sample $\mathbf{Y}_{n}=(Y_{1},\ldots,Y_{n})$, of features with frequency $r>0$ in $\mathbf{Y}_{n}$ that would be observed in one additional sample (Chao and Colwell (2017)); ii) solving the optimal stopping problem discussed in Section 4; iii) estimating the conditional expected number, given an observable sample $\mathbf{Y}_{n}$, of new features that would be observed in $m>1$ additional (unobservable) samples (Colwell et al. (2012)).

Appendix A Appendix - Proofs

A.1 Nonparametric empirical Bayes

In the present section we derive $\hat{M}_{n}$ as a nonparametric empirical Bayes estimator in the sense of Efron and Morris (1973). Recall that $F_{j}$ denotes the label of feature $j$, in the sequel assumed to be $F_{j}\in(0,1)$. $Y_{i}=(Y_{i,j})_{j\geq 1}$ denotes the observation of $i$-th individual. Each entry $Y_{i,j}$ in this sequence is distributed according to a Bernoulli distribution with parameter $p_{j}$ and is independent of the others. $Y_{i,j}=1$ (resp. $Y_{i,j}=0$) indicates the presence (resp. absence) of feature $j$ in the $i$–th individual. The Bayesian nonparametric approach to estimate the missing mass $M_{n}$ requires a prior specification for the $p_{j}$’s: we resort to the three-parameter Beta process introduced by Teh and Görür (2009). Such a prior distribution is defined as the distribution of a Completely Random Measure (CRM) $\tilde{\mu}$ (see e.g. Daley and Vere-Jones (2008)), with Lévy intensity $\nu({\rm d}s,{\rm d}f)=\rho(s){\rm d}s\mathbb{I}{(0,1)}(f){\rm d}f$, where $\rho(s)=\vartheta/\Gamma(1-\alpha)s^{-1-\alpha}(1-s)^{\beta+\alpha-1}$, being $\alpha,\beta\in(0,1)$ and $\vartheta>0$.

Let $\mathbf{Y}_{n}$ be a random sample of size $n$ and $F_{1}^{*},\ldots,F_{K_{n}}^{*}$ be the $K_{n}$ distinct features observed in it, the posterior estimate of $M_{n}$, under a square loss function, equals

[TABLE]

We can characterize the posterior distribution of $\tilde{\mu}$ in (10) resorting to (James, 2017, Theorem 3.1 (ii)), which gives

[TABLE]

where the $J_{\ell}$’s are jumps having a density on the positive real line and $\tilde{\mu}_{n}$ is a CRM having updated Lévy intensity given by $\nu_{n}({\rm d}s\,{\rm d}f)=(1-s)^{n}\rho(s){\rm d}s\mathbb{I}{(0,1)}(f)$ ${\rm d}f$. Hence the Bayesian estimator in (10) boils down to

[TABLE]

where the second equality follows from the fact that the base measure of $\tilde{\mu}_{n}$ is diffuse and the last equality follows from the availability of the Laplace functional and from standard calculations. We are now going to show that $\hat{\vartheta}_{n}=(K_{n,1}\Gamma(\beta+n+1))/(\Gamma(\beta+\alpha+n)n)$ is a consistent estimator of $\vartheta$, then we will conclude that the empirical Bayes estimator $\tilde{M}_{n}(\hat{\vartheta}_{n},\alpha,\beta)$ coincides with $\hat{M}_{n}$. The consistency of $\hat{\vartheta}_{n}$ can be established resorting to the regular variation theory by Karlin (1967); Gnedin et al. (2007). More specifically we are able to show that the tail integral of $\rho(s)$ is a regularly varying function, having regular variation index $\alpha$, indeed

[TABLE]

fulfills

[TABLE]

Hence we just proved that $\overline{\rho}(x)\simeq x^{-\alpha}\vartheta/(\Gamma(1-\alpha)\alpha)$ as $n\to+\infty$, and therefore an application of (Gnedin et al., 2007, Corollary 21) gives the asymptotic relation $\frac{K_{n,1}}{n^{\alpha}}\rightarrow\vartheta$. Moreover, from Stirling’s approximation, $\frac{\Gamma(\alpha+\beta+n)}{\Gamma(n+1+\beta)}\simeq n^{\alpha-1}$. Therefore, $\hat{\vartheta}_{n}=\frac{K_{n,1}\Gamma(\beta+n+1)}{\Gamma(\beta+\alpha+n)n}\rightarrow\vartheta$ and $\hat{\vartheta}_{n}$ is a consistent estimator for $\vartheta$. Plugging $\hat{\vartheta}_{n}$ into the posterior estimate $\tilde{M}_{n}(\vartheta,\alpha,\beta)$, we obtain $\tilde{M}_{n}(\hat{\vartheta}_{n},\alpha,\beta)=\frac{K_{n,1}}{n}=\hat{M}_{n}$. ∎

A.2 Proof of Theorem 2.1

We prove part i) (upper bound) and ii) (lower bound) of Theorem 2.1 separately.

We first focus on the upper bound for the $L_{2}$-risk of $\hat{M}_{n}$. In the sequel we denote by $M_{n,1}$ the total mass of features with frequency 1 in a sample of size $n$, formally $M_{n,1}=\sum_{j\geq 1}p_{j}\mathbb{I}{\{X_{n,j}=1\}},$ we observe that the expectation of $M_{n,1}$ equals

[TABLE]

and obviously $\mathbb{E}[M_{n,1}]\leq W$, for any sequence $(p_{j})_{j\geq 1}$ belonging to the class $\mathscr{P}_{W}$.

In order to bound the $L_{2}$-risk of the estimator $\hat{M}_{n}$ from above for any $(p_{j})_{j\geq 1}\in\mathscr{P}_{W}$, we remind that $R_{n}(\hat{M}_{n};(p_{j})_{j\geq 1})=\mathbb{E}(\hat{M}_{n}-M_{n})^{2}$ may be seen as the sum of the squared bias and the variance. The upper bound for the bias can be easily proved as follows

[TABLE]

As for the variance bound, we note that

[TABLE]

where we have exploited the independence of the random variables $X_{n,j}$. We now observe that the events $\{X_{n,j}=1\}$ and $\{X_{n,j}=0\}$ are incompatible, hence we get

[TABLE]

Putting together the bound on the bias and the variance we immediately obtain (4).

We now focus on the proof of the lower bound of the minimax risk $R^{\ast}_{n}$, which has been defined in (3). Let us first notice that, since $(p_{j})_{j\geq 1}$ belongs to $\mathscr{P}_{W}$, then $M_{n}\leq W$ almost surely, and hence $(M_{n}-\hat{T}_{n})^{2}\geq(M_{n}-\min(\hat{T}_{n},W))^{2}$. Therefore, we can restrict the minimum in (3) over all the possible estimators less than $W$, hence the problem boils down to determine an estimate from below for $\inf_{\hat{T}_{n}\leq W}\sup_{(p_{j})_{j}\in\mathscr{P}_{W}}\mathbb{E}[(M_{n}-\hat{T}_{n})^{2}]$. The proof of this estimate relies on a Bayesian parametric approach based on the randomization of the probability masses $p_{j}$’s. More specifically we consider the randomized vector $\mathsf{P}=(\mathsf{p}_{1},\cdots,\mathsf{p}_{m},0,0,\ldots)$ where $m$ is a Poisson random variable with mean $\lambda=nW$ and $\mathsf{p}_{j}$ are independent beta random variables with parameters $a=1$ and $b=2n-1$. We denote by $\mathsf{P}_{F}:=\sum_{j=1}^{m}\mathsf{p}_{j}$ the total mass of the random sequence $(\mathsf{p}_{j})_{j\geq 1}$; we further define $X:=(X_{n,j})_{j\geq 1}$ and it is worth noticing that, conditionally on $\mathsf{P}$, the $X_{n,j}$’s are independent, having a Binomial distribution with parameters $(n,\mathsf{p}_{j})$. For generic random variables $U,V$ we will write $E_{V}$ (resp. $E_{V|U}$) to mean that the expectation is made with respect to $V$ (resp. to $V$ given $U$). It is easily shown that

[TABLE]

We bound separately the two terms on the r.h.s. of (11). We start by deriving an upper bound for the term on the right. Using Fubini’s Theorem, it comes that

[TABLE]

which leads to the more convenient form

[TABLE]

In order to provide an upper bound for the probability of the event $\{\mathsf{P}_{F}>t\}$, we use a Markov inequality with third centralized moment. We first evaluate the mean of the total mass,

[TABLE]

Now, denoting $\mu_{3}(Y)$ the centralized third moment of the r.v. $Y$, the law of total cumulant implies that

[TABLE]

$m$ is Poisson distributed, so its variance and third centralized moment are equal to its mean $\lambda=nW$. Further noticing that $\frac{2n-1}{2n+1}<1$, we upper bound the third centralized moment by $\mathbb{E}[\mu_{3}(\mathsf{P}_{F}|m)]+\frac{W}{2n^{2}}.$ Conditioning on $m$, the total mass $\mathsf{P}_{F}$ is a sum of independent and identically distributed random variables so $\mu_{3}(\mathsf{P}_{F}|m)=m\mu_{3}(p_{1})$. Since $p_{1}$ is Beta distributed with parameters $a$, and $b$, its third centralized moment is given by

[TABLE]

Using the fact that $\frac{(2n-1)(n-1)}{(2n+1)(n+1)}\leq 1\leq 2$, we successively find that $\mathbb{E}[\mu_{3}(\mathsf{P}_{F}|m)]\leq\frac{W}{2n^{2}}$ and finally $\mu_{3}(\mathsf{P}_{F})\leq\frac{W}{n^{2}},$ obtaining the inequality we need to bound $\mathbb{P}(\mathsf{P}_{F}>t).$ Indeed, Markov’s inequality implies that

[TABLE]

Now we can use the last inequality in (12) and find

[TABLE]

We now provide a lower bound for the first term on the r.h.s. of (11)

[TABLE]

Let $m\in\mathbb{N}$, with the convention that the risk is [math] when $m=0$. We will now lower bound inside the first expectation in the previous expression, using the fact that given $X$ and $m$, we know the posterior of $\mathsf{P}_{F}$.

[TABLE]

where we used the fact that $\mathsf{p}_{j}|X_{n,j}=x$ is a beta random variable ${\rm Beta}(x+a,n-x+b)$. Then,

[TABLE]

where we observed that $(1-\mathsf{p}_{j})\sim{\rm Beta}(b,a)$ to obtain (14). Therefore we can derive the bound

[TABLE]

Together with (13), the previous inequality implies that

[TABLE]

Plugging in (11) the previous results, we obtain

[TABLE]

which concludes the proof. Notice that here we are only interested in showing that the minimax rate is or order $\frac{W}{n}$. We could have obtained sharper bounds (with a better constant) by for instance taking $b=(1+\epsilon)n-1$ instead of $b=2n-1$.

∎

A.3 Details for the determiantion of (6) and (7)

The concentration inequalities (6)–(7) can be proved my means of suitable bounds on the log-Laplace transform of $\hat{W}_{n}-W$, using the techniques of Boucheron et al. (2013). To this end we recall that a random variable $Z$ is sub-Gaussian on the right tail (resp. left tail) with variance factor $v$ if for all $\lambda\geq 0$ (resp. $\lambda\leq 0$) one has

[TABLE]

We will also say that a random variable $Z$ is sub-Gamma on the right tail with variance factor $v$ and scale parameter $c$ if

[TABLE]

Furthermore $Z$ is said to be sub-Gamma on the left tail with variance factor $v$ and scale parameter $c$, if $-Z$ is sub-Gamma on the right tail with variance factor $v$ and scale parameter $c$. These types of bounds on the log-Laplace transform imply exponential tail bounds, useful to find confidence intervals, for example if $Z$ is sub-Gamma on the right tail as above then $\mathbb{P}(Z>\mathbb{E}Z+\sqrt{2vs}+cs)\leq e^{-s}$ for any $s\geq 0$. Analogously if $Z$ is sub-Gaussian on the left tail (i.e. (15) holds true for all $\lambda\leq 0$), then $\mathbb{P}(Z<\mathbb{E}Z-\sqrt{2vs})\leq e^{-s}$ for any $s\geq 0$. See Boucheron et al. (2013) for additional details on this subject. In the next proposition we prove that $\hat{W}_{n}$ is sub-Gamma on the right tail and sub-Gaussian on the left tail, providing suitable bounds on the log-Laplace transform. As just explained these bounds on the log-Laplace immediately implies the exponential tail bounds (6) and (7).

Proposition A.1

The random variable $\hat{W}_{n}:=\sum_{j\geq 1}X_{n,j}/n$ is an unbiased estimator of $W$. In addition $\hat{W}_{n}$ is sub-Gamma on the right tail with variance factor $w_{n}^{+}:=2W/n$ and scale factor $1/n$, i.e. for any $0<\lambda<n$

[TABLE]

on the left tail, $\hat{W}_{n}$ is sub-Gaussian with variance factor $w_{n}^{-}=W/n$, namely for all $\lambda\geq 0$ we have

[TABLE]

Proof.

One can immediately show that $\mathbb{E}(\hat{W}_{n})=W$. We first prove (16), that is to say $\hat{W}_{n}-W$ is sub-Gamma on the right tail, indeed, by the independence of the $X_{n,j}$’s, we get

[TABLE]

having observed that $X_{n,j}$ is a Binomial random variable with parameters $n$ and $p_{j}$. An application of the inequality $\log(z)\leq z-1$, for $z>0$, yields

[TABLE]

where $\phi(t)=e^{t}-1-t$. By a series expansion of the exponential function, we obtain

[TABLE]

and then (16) follows for any $\lambda\in(0,n)$. The proof of (17) goes along similar lines, indeed proceeding as before one can prove that for any $\lambda\geq 0$

[TABLE]

and then (17) follows by observing that $\phi(-\lambda/n)\leq\lambda^{2}/(2n^{2})$. ∎

A.4 Proof of Theorem 3.1 and related results

In this section we focus on Theorem 3.1. First of all we introduce some preliminary results which are useful to prove the theorem and are interesting in their own right. In the next proposition we derive concentration inequalities for both $K_{n,r}$ and $K_{n}$.

Proposition A.2

For any $n\geq 1$, $r\geq 1$ and $\delta\in(0,1)$ we have

[TABLE]

and

[TABLE]

Proof.

Let us focus on the derivation of (19). First of all we will determine a variance bound for the random variable $K_{n,r}$, which will be employed to derive a corresponding bound on the log-Laplace transform of $(K_{n,r}-\mathbb{E}(K_{n,r}))$. Then the result will follow by a suitable application of the Bernstein inequality.

Thanks to the independence of the random variables $X_{n,j}$’s, the variance of $K_{n,r}$ may be bounded as follows

[TABLE]

We now establish the bound on the log-Laplace. Since $K_{n,r}$ is the sum of independent random variables, for any $t\in\mathds{R}$ we can write:

[TABLE]

where we have implicitly defined the function $\phi(t)=e^{t}-1-t$ and we have applied the Bennett inequality. The validity of the previous bound on the log-Laplace implies that for any $x\geq 0$

[TABLE]

thanks to the Bernstein inequality. If one defines $s:=x^{2}/(2(\mathbb{E}(K_{n,r})+x/3))$ the previous inequality boils down to the following one

[TABLE]

which in turn implies the validity of (19), as a consequence of the change of variable $2e^{-s}=\delta$ and the elementary inequality $\sqrt{a+b}\leq\sqrt{a}+\sqrt{b}$, valid for any positive real numbers $a,b$. The proof of (20) goes along similar lines, having observed that $\operatornamewithlimits{Var}(K_{n})\leq\mathbb{E}(K_{n})$. ∎

The next Remark is a simple consequence of Proposition A.2 and will be used in the derivation of the confidence interval for $\hat{M}_{n}$ to upper bound the expected values $\mathbb{E}(K_{n,r})$ and $\mathbb{E}(K_{n})$ in terms of the two statistics $K_{n,r}$ and $K_{n}$ respectively.

Remark A.1

The concentration inequality for $K_{n,r}$ (19) implies that with probability bigger than $1-\delta$

[TABLE]

Solving the previous inequality with respect to $\sqrt{\mathbb{E}(K_{n,r})}$, we can conclude that with probability at least $1-\delta$

[TABLE]

Analogously we can employ the concentration inequality for $K_{n}$ (see (20)) to state that the event

[TABLE]

has probability at least $1-\delta$.

The confidence interval we suggest is mainly based on a concentration inequality for $\hat{M}_{n}-M_{n}$. Such concentration result is proved my means of suitable bounds on the log-Laplace transform of $\hat{M}_{n}-M_{n}$, as in Section A.3, more precisely in Proposition A.3 we will prove that $\hat{M}_{n}-M_{n}$ is a sub-Gamma random variable on the right and the left tails (see Boucheron et al. (2013)). In order to provide the bounds on the log-Laplace transform of $\hat{M}_{n}-M_{n}$, the following lemma is required.

Lemma A.1

For any $n\geq 1$, we have

[TABLE]

and

[TABLE]

Proof.

First of all, we focus on the proof of (22), straightforward calculations lead to

[TABLE]

where we have observed that $(1-p_{j})^{n+1-2}\leq(1-p_{j})^{n-2}$, which is legitimated by the fact that $n$ is a sample size and can be assumed sufficiently large.

Second, we focus on the proof of the variance bound (23). Exploiting the independence of the random variables $X_{n,j}$’s, we get

[TABLE]

We now observe that the events $\{X_{n,j}=1\}$ and $\{X_{n,j}=0\}$ are incompatible, hence

[TABLE]

and the bound on the variance now follows.

∎

To fix the notation, it is worth to define the quantities

[TABLE]

for any $k\geq 1$. It can be easily seen that $\mathbb{E}[K_{k}(n)]\leq\mathbb{E}[K(n)]$, and in addition $\mathbb{E}[K(n)]$ is an increasing function of $n$. We are now ready to prove that $\hat{M}_{n}-M_{n}$ is sub-Gamma on its tails.

Proposition A.3

The variable $D_{n}=\hat{M}_{n}-M_{n}$ is sub–gamma on the right tail with variance factor $v_{n}^{+}=2\{\mathbb{E}(K_{n,1})+2n^{2}\mathbb{E}(K_{n,2})/[n\cdot(n-1)]\}/n^{2}$ and scale factor $1/n$, i.e. for any $0<\lambda<n$

[TABLE]

On the left tail, $D_{n}$ is again sub–gamma with variance factor $v_{n}^{-}=4\mathbb{E}[K(n)]/(n-1)^{2}$ and scale parameter $1/(n-1)$, i.e. for any $0<\lambda<(n-1)$

[TABLE]

Proof.

First we prove (24), observe that when $\lambda>0$

[TABLE]

For the sake of simplifying notation, define $H_{j}:=\mathbb{I}{\{X_{n,j}=1\}}-np_{j}\mathbb{I}{\{X_{n,j}=0\}}\leq 1$, hence we have

[TABLE]

where we used the inequality $\log(z)\leq z-1$ for any positive $z$. Now, observing that the function $f(z):=\frac{e^{z}-1-z}{z^{2}}$ is increasing and that $H_{j}\leq 1$, we get

[TABLE]

now one can proceed along similar arguments as in the proof of (23) to obtain

[TABLE]

hence (24) now follows.

Finally we prove (25), analogous calculations as before show that

[TABLE]

having defined $W_{j}:=\mathbb{I}{\{X_{n,j}=0\}}-\frac{1}{np_{j}}\mathbb{I}{\{X_{n,j}=1\}}\leq 1$. Thanks to the inequality $\log(z)\leq z-1$, which is valid for any positive $z$, and applying similar considerations as those used in the first part of the proof, we have

[TABLE]

where we have observed that $(1-p)^{n}\leq e^{-np}$. Recalling the definition of $K_{k}(n)$ and that $\mathbb{E}[K_{k}(n)]\leq\mathbb{E}[K(n)]$, we get

[TABLE]

hence (25) has been now established. ∎

Finally we prove Theorem 3.1.

Proof of Theorem 3.1.

The result is a consequence of the two bounds that we are now going to prove separately:

[TABLE]

where $\hat{m}_{n},\hat{v}^{-}_{n}$ and $\hat{v}^{+}_{n}$ are suitable quantities that will be defined in the proof. First we discuss how to determine (26), we remind that $D_{n}=\hat{M}_{n}-M_{n}$ is a sub–Gamma random variable, as shown in Proposition A.3, hence the following holds (see Ben-Hamou et al. (2017))

[TABLE]

for any $s\geq 0$. Putting $e^{-s}=\delta$, we obtain

[TABLE]

On the basis of (22), we can upper bound $\mathbb{E}(D_{n})$ in (28) as

[TABLE]

and, as a consequence of Remark A.3.2, the following inequality is true with probability bigger than $1-\delta$

[TABLE]

Remark A.3.2 also applies to upper bound the expected values in $v_{n}^{+}$ with the corresponding empirical values, indeed the following bound holds true

[TABLE]

with probability bigger than $1-2\delta$. Using the inequalities (29)–(30) in (28) we can conclude

[TABLE]

which gives (26).

Now we focus on the derivation of (27), first of all recall that by Proposition A.3, $-D_{n}$ is a sub–Gamma random variable, hence

[TABLE]

for any $s\geq 0$. Define $e^{-s}=\delta$, then the previous bound implies that

[TABLE]

since $\mathbb{E}[D_{n}]\geq 0$, we get

[TABLE]

One can now apply the bound contained in (Gnedin et al., 2007, Lemma 1)

[TABLE]

which entails that $\mathbb{E}[K(n)]\leq\mathbb{E}(K_{n})/(1-2/n)$, and hence

[TABLE]

Thanks to Remark A.3.2, with probability bigger than $1-\delta$, the following holds

[TABLE]

Using (32) in (31), we get

[TABLE]

which in turn implies (27). Putting together (26)–(27), the thesis now follows. ∎

Acknowledgements

Federico Camerlenghi and Stefano Favaro received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme under grant agreement No 817257. Federico Camerlenghi and Stefano Favaro gratefully acknowledge the financial support from the Italian Ministry of Education, University and Research (MIUR), “Dipartimenti di Eccellenza” grant 2018–2022.