Minimax Estimation of Functionals of Discrete Distributions

I Introduction and main results

Given $n$ independent samples from an unknown discrete probability distribution $P=(p_{1},p_{2},\ldots,p_{S})$, with unknown support size $S$, consider the problem of estimating a functional of the distribution of the form:

which plays significant roles in information theory . Another information theoretic quantity which is closely related to the entropy is the mutual information, which for discrete random variables can be defined as

where $P_{X},P_{Y},P_{XY}$ denote, respectively, the distributions of random variables $X$, $Y$, and the pair $(X,Y)$.

We are also interested in the family of information measures $F_{\alpha}(P)$:

The significance of functional $F_{\alpha}(P)$ can be seen via the connection $H_{\alpha}(P)=\frac{\ln F_{\alpha}(P)}{1-\alpha}$, where $H_{\alpha}(P)$ is the Rényi entropy , emerging in answering fundamental questions in information theory ,. The functional $1-F_{2}(P)$ is also called Gini impurity, which is widely used in machine learning .

Over the years, the use of information theoretic measures, especially entropy and mutual information, has extended far beyond the information theory community, and is deeply imbued in fundamental concepts from various disciplines. In statistics, one of the popular criteria for objective Bayesian modeling is to design a prior on the parameter to maximize the mutual information between the parameter and the observations. In machine learning, the so-called infomax criterion states that the function that maps a set of input values to a set of output values should be chosen or learned so as to maximize the mutual information between the input and output, subject to a set of specified constraints. This principle has been widely adopted in practice, for example, in decision tree based algorithms in machine learning such as C4.5 , one tries to select the feature at each step of tree splitting to maximize the mutual information (called information gain principle ) between the output and the feature conditioned on previous chosen features. Other measures in feature selection have been proposed, such as the Gini impurity (used in CART ), variance reduction , and many of them can be incorporated as special cases of what we study in this paper. We emphasize that in some applications, mutual information arises naturally as the only answer, for example, the well known Chow–Liu algorithm for learning tree graphical models relies on estimation of the mutual information, which is a natural consequence of maximum likelihood estimation. Recently, it was shown that mutual information is the unique measure of relevance for inference in the presence of side information to satisfy a natural data processing property.

We also mention genetics , image processing , computer vision , secrecy , ecology , and physics as fields in which information theoretic measures are widely used. There are some other functionals that can be loosely categorized as information theoretic measures, such as the association measures quantifying certain dependency relations of random variables , and divergence measures .

In most applications, the underlying distribution is unknown, so we cannot compute these information theoretic measures exactly. Hence, in nearly every problem that uses information theoretic measures, we need to estimate these quantities from the data, which is what we study in this paper. Our contributions are threefold. (i) We show that when the number of observations $n$ is comparable to the parameter dimension (a relevant regime in the “big data” era), the prevailing approaches (such as plug-in of the maximum likelihood estimator) can be highly sub-optimal. (ii) We propose new and computationally efficient algorithms that are essentially optimal in terms of the worst case squared error risk. That is, we both characterize the fundamental limits on estimation performance, and propose practical algorithms that essentially achieve them. Our results establish that for such functional estimation scenarios, replacing the plug-in (maximum likelihood) estimator by our practical and essentially minimax optimal estimators yields an effective enlargement of the sample size from $n$ to $n\ln n$, which can make a significant difference in practice. (iii) We demonstrate the efficacy of our schemes by a comparison with existing procedures in the literature, as well as illustrate performance boosts over traditional schemes on both real and simulated data.

Notation: We use the notation $a_{\gamma}\lesssim b_{\gamma}$ to denote that there exists a universal constant $C$ such that $\sup_{\gamma}\frac{a_{\gamma}}{b_{\gamma}}\leq C$. Notation $a_{\gamma}\asymp b_{\gamma}$ is equivalent to $a_{\gamma}\lesssim b_{\gamma}$ and $b_{\gamma}\lesssim a_{\gamma}$. Notation $a_{\gamma}\gg b_{\gamma}$ means that $\liminf_{\gamma}\frac{a_{\gamma}}{b_{\gamma}}=\infty$, and $a_{\gamma}\ll b_{\gamma}$ is equivalent to $b_{\gamma}\gg a_{\gamma}$. The sequences $a_{\gamma},b_{\gamma}$ are non-negative.

Our main goal in this work is to present a general approach to the construction of minimax rate-optimal estimators for functionals of the form (1) under $L_{2}$ loss. To illustrate our approach, we describe and analyze explicit constructions for the specific cases of entropy $H(P)$ and $F_{\alpha}(P)$, from which the construction for any other functional of the form (1) will be clear. Our estimators for each of these two functionals are agnostic with respect to the support size $S$, and achieve the minimax $L_{2}$ rates (i.e. the performance of our approaches when we do not know the support size $S$ does not degrade compared with the case where the support size $S$ is known).

Our approach is to tackle the estimation problem separately for the cases of “small $p$” and “large $p$” in $H(P)$ and $F_{\alpha}(P)$ estimation, corresponding to treating regions where the functional is nonsmooth and smooth in different ways. As we describe in detail in the sections to follow, where we give a full account of our estimators, in the nonsmooth region, we rely on the best polynomial approximation of the function $f$, by employing an unbiased estimator for this approximation. The best polynomial approximation for a function $f(x)$ on domain $A$ with order no more than $K$ is defined as

where $\mathsf{poly}_{K}$ is the collection of polynomials with order at most $K$ on $A$. The part pertaining to the smooth region is estimated by a bias-corrected maximum likelihood estimator. We apply this procedure coordinate-wise based on the empirical distribution of each observed symbol, and finally sum the respective estimates.

We now look at the specific cases of entropy and $F_{\alpha}(P)$ separately. For the entropy, after we obtain the empirical distribution $P_{n}$, for each coordinate $P_{n}(i)$, if $P_{n}(i)\ll\ln n/n$, we (i) compute the best polynomial approximation for $-p_{i}\ln p_{i}$ in the regime $0\leq p_{i}\ll\ln n/n$, (ii) use the unbiased estimators for integer powers $p_{i}^{k}$ to estimate the corresponding terms in the polynomial approximation for $-p_{i}\ln p_{i}$ up to order $K_{n}\sim\ln n$, and (iii) use that polynomial as an estimate for $-p_{i}\ln p_{i}$. If $P_{n}(i)\gg\ln n/n$, we use the estimator $-P_{n}(i)\ln P_{n}(i)+\frac{1}{2n}$ to estimate $-p_{i}\ln p_{i}$. Then, we add the estimators corresponding to each coordinate. Our estimator for $F_{\alpha}(P)$ is very similar to that of entropy, with the only difference that we conduct polynomial approximation for $x^{\alpha}$ with order $K_{n}\sim\ln n$, and use the estimator $\left(1+\frac{\alpha(1-\alpha)}{2nP_{n}(i)}\right)P_{n}^{\alpha}(i)$ when $P_{n}(i)\gg\ln n/n$.

We remark that our estimator is both conceptually and algorithmically simple, with complexity linear in the number of samples $n$. Indeed, the only non-trivial computation required is the best polynomial approximation for functions, which is data independent and can be done offline before obtaining any samples from the experiment. Moreover, the coefficients of the best polynomial approximation of different orders can be preprocessed and stored in advance in the implementation of our approach. We demonstrate in Section V that the best polynomial approximation step can be performed efficiently using modern machinery from approximation theory and numerical analysis.

I-B Main results

Simple as our estimators are to describe and implement, they can be shown to be near “optimal” in the strong sense we now describe. We adopt the conventional statistical decision theoretic framework . Regarding the task of estimating functional $F(P)$, the $L_{2}$ risk of an arbitrary estimator $\hat{F}$ is defined as

where the expectation is taken with respect to the distribution $P$ that generates the observations used by $\hat{F}$. Apparently, the $L_{2}$ risk is a function of both the unknown distribution $P$ and the estimator $\hat{F}$, and our goal is to minimize this risk. Since $P$ is unknown, we cannot directly minimize it, but if we want to do well no matter what the true distribution $P$ is, we may want to adopt the minimax criterion , and try to minimize the maximum risk

where $\mathcal{M}_{S}$ denotes the set of all discrete distributions with support size $S$. The estimator that minimizes the maximum risk above is called the minimax estimator, and the corresponding risk is called the minimax risk. The exact computation of the minimax risk and the minimax estimator for general $F(P)$ seems intractable. Although the maximum risk in (7) is a convex function of $\hat{F}$ (supremum of convex functions is convex), minimizing this function involves computation of the objective function via $\sup_{P\in\mathcal{M}_{S}}$, which is a non-convex optimization problem. Moreover, even if we can compute it exactly, the minimax estimator will surely depend on the support size $S$, which is unknown to the statistician in many applications.

Hence, we slightly relax the requirement, and seek minimax rate-optimal estimators $\hat{F}^{*}$ with maximum (worst-case) risk equal to the minimax risk up to a multiplicative constant. In other words, we want to design estimator $\hat{F}^{*}$ such that there exist two universal positive constants $0<C_{1}\leq C_{2}<\infty$ that do not depend on the problem configuration (such as the support size $S$ and sample size $n$), for which

As it turns out, it is possible to construct estimators $\hat{F}^{*}$ for a wide class of functionals, which do not rely on the knowledge of support size $S$. A brief description of the constructions is given in Section I-A. We find it intriguing that our estimators, which are minimax rate-optimal, are intimately connected to the problem of best (minimax) polynomial approximation, which is a convex optimization problem. In some sense, we have transformed the difficult-to-solve minimax and convex problem of minimizing the maximum risk in (7) into another efficently solvable minimax and convex problem of minimizing the maximum deviation of a polynomial from a given function, which turns out to have been studied extensively in approximation theory for more than a century.

To ease the presentation, we consider the “Poissonized” observation model [22, Pg. 508], since we can show that the minimax risks under the Multinomial model and Poisson model are essentially the same (cf. Lemma 16). Moreover, adopting the Poisson model significantly reduces the length of the proofs, and we emphasize that similar analysis can also go through for Multinomial settings, with more nuanced analysis. In the Poisson setting, we first draw a Poisson random number $N\sim\mathsf{Poi}(n)$, and then conduct the sampling $N$ times. Consequently, the observed number of occurrences of each symbol are independent [23, Thm. 5.6].

We have the following characterization of the minimax risk for entropy estimation.

Suppose $n\gtrsim\frac{S}{\ln S}$. Then the minimax risk of estimating entropy $H(P)$ satisfies

Our estimator achieves this bound without knowledge of the support size $S$ under the Poisson model.

The following is an immediate consequence of Theorem 1.

For our entropy estimator, the maximum $L_{2}$ risk vanishes provided $n\gg\frac{S}{\ln S}$. Moreover, if $n\lesssim\frac{S}{\ln S}$, then the maximum risk of any estimator for entropy is bounded from zero.

It was first shown in that one must have $n\gg\frac{S}{\ln S}$ for consistently estimating the entropy. However, the entropy estimators based on linear programming proposed in Valiant and Valiant have not been shown to achieve the minimax risk. Another estimator proposed by Valiant and Valiant has only been shown to achieve the minimax risk in the restrictive regime of $\frac{S}{\ln S}\lesssim n\lesssim\frac{S^{1.03}}{\ln S}$. Wu and Yang independently applied the idea of best polynomial approximation to entropy estimation, and obtained its minimax $L_{2}$ rates. The minimax lower bound part of Theorem 1 follows from Wu and Yang . We also remark that, unlike the estimator we propose, the estimator in Wu and Yang relies on knowledge of the support size $S$, which generally may not be known.

For the functional $F_{\alpha}(P),0<\alpha<1$, we have the following.

Suppose $n\gtrsim\frac{S^{1/\alpha}}{\ln S}$ when we estimate $F_{\alpha}(P),0<\alpha<1$. Then we have the following characterizations of the minimax risk.

$0<\alpha\leq 1/2$. If we also have $\ln n\lesssim\ln S$, then

Our estimators $\hat{F}_{\alpha}$ achieve this bound without knowledge of the support size $S$ under the Poisson model.

One immediate corollary of Theorem 2 is the following.

For our estimators of $F_{\alpha}$, the maximum $L_{2}$ risk vanishes provided $n\gg\frac{S^{1/\alpha}}{\ln S},0<\alpha<1$. Moreover, if $n\lesssim\frac{S^{1/\alpha}}{\ln S}$, then the maximum risk of any estimator for $F_{\alpha}$ is bounded from zero.

The minimax lower bound we present in Theorem 2In a previous version of the manuscript, there is a $\sqrt{\ln S}$ gap between our minimax lower bound and the achievability in Theorem 2. Partially inspired by Wu and Yang , we modified the proof by using an argument similar to that in the lower bound proof of , thereby closing the gap. significantly improves on Paninski’s lower bound in , which states that if $n\lesssim S^{1/\alpha-1}$, then the maximum $L_{2}$ risk of any estimator for $F_{\alpha}(P),0<\alpha<1$, is bounded from zero.

The next two theorems correspond to estimation of $F_{\alpha}(P)$, $\alpha>1$.

Suppose $1<\alpha<\frac{3}{2}$. Under the Poissonized model, our estimator $\hat{F}_{\alpha}$ satisfies

In other words, our estimator $\hat{F}_{\alpha},1<\alpha<3/2$ achieves an $L_{2}$ convergence rate of $(n\ln n)^{-2(\alpha-1)}$ regardless of the support size. This also turns out to be the minimax rate, as shown by the following result.

Suppose $1<\alpha<\frac{3}{2}$. There exists a universal constant $c_{0}>0$ such that if $S=c_{0}n\ln n$ then

where the infimum is taken over all possible estimators $\hat{F}$.

Table I summarizes the minimax $L_{2}$ rates and the $L_{2}$ convergence rates of the MLE in estimating $F_{\alpha}(P),\alpha>0$ and $H(P)$. When the $L_{2}$ rates have two terms, the first and second terms represent respectively the contributions of the bias and the variance. When there is a single term, only the dominant term is retained. Conditions for these results are presented in parentheses.

From a sample complexity perspective (i.e. how should the number of samples $n$ scale with the support size $S$ to achieve consistent estimation), Table I implies the results in Table II.

Our work (including the companion paper ) is the first to obtain the minimax rates, minimax rate-optimal estimators, and the maximum risk of MLE for estimating $F_{\alpha}(P),0<\alpha<3/2$, and entropy $H(P)$ in the most comprehensive regime of $(S,n)$ pairs. Evident from Table I is the fact that the MLE cannot achieve the minimax rates for estimation of $H(P)$, and $F_{\alpha}(P)$ when $0<\alpha<3/2$. In these cases, our estimators have performance with $n$ samples essentially the same as the MLE with $n\ln n$ samples, and it is the best possible. In other words, the minimax rate-optimal schemes enlarge the “effective sample size” from $n$ to $n\ln n$. Furthermore, all the improvements we have are in the bias, which is the dominating factor in the risk. This observation suggests a simple way to obtain the minimax $L_{2}$ rates from the $L_{2}$ rates of the MLE. One need merely find the bias term in the expression of MLE $L_{2}$ rates, and replace the term $n$ by $n\ln n$. This simple rule is intimately connected to the rationale behind the construction and analysis of our estimators, on which we elaborate in Section II.

We also note that Table II is a “lossy compression” of Table I. Indeed, it did not reflect the important improvement of our estimator over the MLE in estimating $F_{\alpha}(P),1<\alpha<3/2$. However, it is more transparent about the increase of difficulty in estimation when we decrease $\alpha$. Indeed, when $\alpha\to 0^{+}$, the sample complexity in estimating $F_{\alpha}(P)$, $S^{1/\alpha}/\ln S$ becomes super-polynomial in $S$, which implies that the problem has become extremely challenging. Indeed, the limiting case of $\alpha=0$ corresponds to estimating the support size of a discrete distribution, which has long been known impossible to do consistently without additional assumptions .

I-C Discussion of main results

Within the scope of estimating entropy of discrete distributions from i.i.d. samples, the reader should be aware of different problem formulations, so as not to be confused by seemingly contradictory results. The problem we consider in this paper is to estimate the entropy $H(P)$ for all possible distributions $P$ supported on $S$ elements, a setting for which we obtain Theorem 1. However, one may impose some additional structure on the distribution $P$, and thus restrict attention to smaller uncertainty sets of distributions. One may expect different answers depending on the size and nature of the uncertainty sets. One of the popular alternative settings is to assume that the distribution $P$ comes from a uniform distribution with unknown support size $S$. Note that it is a great simplification of the problem, indeed, there is only one parameter $S$ to estimate. Correspondingly, we only need $n\gg\sqrt{S}$ samples to consistently estimate the support size, and the entropy under this setting , which is much smaller than the required $n\gg\frac{S}{\ln S}$ samples in our setting, cf. Theorem 1.

Some readers may be concerned that the minimax decision theoretic framework we adopt is too pessimistic. In some sense, it characterizes the worst-case performance over all possible distributions $P\in\mathcal{M}_{S}$, and it would be disappointing if our estimator fails to behave reasonably for distributions lying in a strict subset of $\mathcal{M}_{S}$ not including the worst case distribution. Regarding this question, Brown argued that the minimax idea has been an essential foundation for advances in many areas of statistical research, including general asymptotic theory and methodology, hierarchical models, robust estimation, optimal design, and nonparametric function analysis. Second, the statistics community in general uses the adaptive estimation framework to alleviate the pessimism of minimaxity . Specifically, one specifies a nested sequence of subsets of $\mathcal{M}_{S}$, and tries to construct an estimator that achieves simultaneously the minimax rates over each of the subsets without knowing the subset to which the active parameter $P$ actually belongs. It was shown recently in another related paper that along the nested subsets $\mathcal{M}_{S}(H)=\{P:H(P)\leq H,P\in\mathcal{M}_{S}\}$, our estimator (without knowing $H$ nor $S$) simultaneously achieves the minimax rates over $P\in\mathcal{M}_{S}(H)$ for all $H\leq\ln S$. Most surprisingly, the maximum risk of our estimator over $\mathcal{M}_{S}(H)$ for every $S$ and $H$ with $n$ samples is still essentially that of the MLE with $n\ln n$ samples, further reinforcing the effectiveness of our estimator.

It is instructive to consider our results in the context of the intriguing connections and differences between three important problems in information theory: entropy estimation, estimating a discrete distribution under relative entropy loss, and minimax redundancy in compressing i.i.d. sources. Table III summarizes the known results.

Table III conveys several important messages. First, in the asymptotic regime, there is a logarithmic factor between the redundancy of the compression and distribution estimation problems. Indeed, since compression requires use of a coding distribution $Q$ that does not depend on the data, the redundancy of compression will definitely be larger than the risk under relative entropy in estimating the distribution. However, in the large alphabet setting, the problems are equally difficult - the phase transition of vanishing risk for both compression and distribution estimation happen when $n$ is linear in the support size $S$.

Second, the large alphabet setting shows that estimation of entropy is considerably easier than both estimating the corresponding distribution, and compression. It is somewhat surprising and enlightening, since there has been a well-received tradition to apply data compression techniques to estimate entropy, even beyond the information theory community, e.g. , whereas one of the implications of Table III is that the approach of entropy estimation via compression can be highly sub-optimal.

If we plot the phase transitions of $\ln n/\ln S$ for estimating $F_{\alpha}(P)$ using $F_{\alpha}(P_{n})$ with respect to $\alpha$, we obtain Figure 1.

consistent estimationnot achievable(Theorem 2)1012$\alpha<span class="katex-display"><span class="katex"><span class="katex-mathml"><math xmlns="http://www.w3.org/1998/Math/MathML" display="block"><semantics><mrow><mfrac><mrow><mi>ln</mi><mo>⁡</mo><mi>n</mi></mrow><mrow><mi>ln</mi><mo>⁡</mo><mi>S</mi></mrow></mfrac></mrow><annotation encoding="application/x-tex">\frac{\ln n}{\ln S}</annotation></semantics></math></span><span class="katex-html" aria-hidden="true"><span class="base"><span class="strut" style="height:2.0574em;vertical-align:-0.686em;"></span><span class="mord"><span class="mopen nulldelimiter"></span><span class="mfrac"><span class="vlist-t vlist-t2"><span class="vlist-r"><span class="vlist" style="height:1.3714em;"><span style="top:-2.314em;"><span class="pstrut" style="height:3em;"></span><span class="mord"><span class="mop">ln</span><span class="mspace" style="margin-right:0.1667em;"></span><span class="mord mathnormal" style="margin-right:0.0576em;">S</span></span></span><span style="top:-3.23em;"><span class="pstrut" style="height:3em;"></span><span class="frac-line" style="border-bottom-width:0.04em;"></span></span><span style="top:-3.677em;"><span class="pstrut" style="height:3em;"></span><span class="mord"><span class="mop">ln</span><span class="mspace" style="margin-right:0.1667em;"></span><span class="mord mathnormal">n</span></span></span></span><span class="vlist-s">​</span></span><span class="vlist-r"><span class="vlist" style="height:0.686em;"><span></span></span></span></span></span><span class="mclose nulldelimiter"></span></span></span></span></span></span>1/\alpha$consistent estimationachievable via both MLEand our scheme(resp. , Theorem 2,3) We observe a sharp phase transition at $\alpha=1$, as the sample size requirement shifts from $n\gg S^{\frac{1}{\alpha}}/\ln S$ to $n\gg 1$, depending on whether $\alpha$ is in the left or right neighborhood of 1, respectively. Hence, $\alpha=1$ is a critical point in that consistent estimation requires a number of measurements super-linear or constant in the size of the alphabet according to whether $\alpha<1$ or $\alpha>1$.

Combining Table III and Figure 1 leads to the interesting observation that, in high dimensional asymptotics, estimating a functional of a distribution could be easier (e.g. $H(P),F_{\alpha}(P),\alpha>1$) or harder (e.g. $F_{\alpha}(P),0<\alpha<1$) than estimating the distribution itself. This observation taps into another interesting interpretation of the functional $F_{\alpha}(P)$. In information theory, the random variable $\imath(X)=\ln\frac{1}{P(X)}$ is known as the information density, and plays important roles in characterizing higher order fundamental limits of coding problems . The functional $F_{\alpha}(P)$ can be interpreted as the moment generating function for random variable $\imath(X)$ as

It is shown in Valiant and Valiant that the distribution of $\imath(X)$ can be estimated using $n\gg S/\ln S$ samples. Since moment generating functions can determine the distribution under some conditions, it is indeed plausible to see that the problem of estimating $F_{\alpha}(P)$, or the moment generating function of $\imath(X)$, is either easier or harder than estimating the distribution of $\imath(X)$ itself for various values of $\alpha$.

For any $\alpha\in(1,3/2)$ and any $\delta>0,\epsilon\in(0,1)$, there exists a constant $c=c_{\alpha}(\delta,\epsilon)>0$ such that,

where $\hat{F}_{\alpha}$ is the estimator from Theorem 3.

In words, with high probability $(\ln\hat{F}_{\alpha})/(1-\alpha)$ is close to the Rényi entropy provided $n\gtrsim S/\ln S$. In contrast, the MLE requires $n\gtrsim S$ samples for estimating $H_{\alpha}(P),1<\alpha<\frac{3}{2}$, as is implied by the following theorem.

For any $\alpha\in(1,3/2)$ and any constant $c>0$, there exist some $\delta=\delta_{\alpha}(c)>0$ such that the MLE $H_{\alpha}(P_{n})$ satisfies

To conclude this discussion, we conjecture that plugging in our minimax rate-optimal estimators for $F_{\alpha}(P)$ into the definition of $H_{\alpha}(P)$ results in minimax rate-optimal estimators for $H_{\alpha}(P)$ for all $\alpha>0$.

II Motivation, methodology, and related work

Existing theory proves inadequate for addressing the problem of estimating functionals of probability distributions. A natural estimator for functionals of the form (1) is the maximum likelihood estimator (MLE), or plug-in estimator, which simply evaluates $F({P}_{n})$, where ${P}_{n}$ is the empirical distribution of the data. How well does the MLE perform? Interestingly, if $f\in C^{1}(0,1]$ and we focus on $n$ i.i.d. observations from a distribution with support size $S$, then the problem of estimating $F(P)$ becomes a classical problem when $S$ is fixed, and the number of observations $n\to\infty$. This maximum likelihood estimator is asymptotically efficient [48, Thm. 8.11, Lemma 8.14] in the sense of the Hájek convolution theorem and the Hájek–Le Cam local asymptotic minimax theorem . It is therefore not surprising to encounter the following quote from the introduction of Wyner and Foster who considered entropy estimation:

“The plug-in estimate is universal and optimal not only for finite alphabet i.i.d. sources but also for finite alphabet, finite memory sources. On the other hand, practically as well as theoretically, these problems are of little interest. ”

In light of this, is it fair to say that the entropy estimation problem is solved in the finite alphabet setting? It was observed in Paninski that the maximum of $\mathsf{Var}(-\ln P(X))$ over distributions with support size $S$ is of order $(\ln S)^{2}$ (a tight bound is also given by Lemma 15 in the appendix). Since classical asymptotics (with the Delta method [48, Chap. 3]) show that

a naive interpretation of (17) might be that it suffices to take $n\gg(\ln S)^{2}$ samples to guarantee the consistency of $H(P_{n})$. Such an interpretation could however be very misleading. It was already observed in Paninski that if $n\lesssim S^{1-\delta},\delta>0$, then the maximum $L_{2}$ risk of any entropy estimator would be unbounded as $S$ grows.

This apparent discrepancy shows that (17) is not valid when $S$ might be growing with $n$, and it is of utmost importance to obtain risk bounds for estimators of entropy and other functionals of distributions in the latter regime. Indeed, in the modern era of high dimensional statistics, we often encounter situations where the support size is comparable to, or much larger than the number of observations. For example, half of the words in the collected works by Shakespeare appeared only once .

It was shown in the companion paper that for $n\gtrsim S$, the maximum risk of the MLE $H(P_{n})$ can be written as

The fact that the bias dominates the risk in entropy estimation in the large alphabet setting has been known, see . However, a general recipe to overcome the bias has defied many attempts. We briefly review some of the approaches in the literature.

One of the earliest investigations on reducing the bias of MLE in entropy estimation is due to Miller , who showed that, for any fixed distribution $P$ supported on $S$ elements, for each symbol $i,1\leq i\leq S$, we have

Another popular approach in estimating entropy is based on using the Dirichlet prior smoothing. Dirichlet smoothing may carry two different meanings in terms of entropy estimation:

One first obtain a Bayes estimate for the discrete distribution $P$, which we denote by $\hat{P}_{B}$, and then plugs it in the entropy functional to obtain the entropy estimate $H(\hat{P}_{B})$.

It was shown in that both approaches require at least $n\gg S$ to be consistent.

The jackknife is another popular technique in reducing the bias. However, it was shown by Paninski that the jackknifed MLE also requires $n\gg S$ samples to consistently estimate entropy.

Given the fact that $n\gg\frac{S}{\ln S}$ is both necessary and sufficient, the approaches mentioned above are far from optimal. We note that the problem of estimating entropy of discrete distributions from i.i.d. observations on large alphabets has been investigated in various disciplines by many authors, with many approaches difficult to analyze theoretically. Among them we mention the Miller–Madow bias-corrected estimator and its variants , the jackknified estimator , the shrinkage estimator , the Bayes estimator under various priors , the coverage adjusted estimator , the Best Upper Bound (BUB) estimator , the B-Splines estimator , and etc.

In what follows, we explain in detail our step by step approach to this problem, and arrive at minimax rate-optimal estimators for the entire range of functional estimation problems considered.

II-B How did we come up with our scheme?

Existing literature implied that it is possible to come up with consistent entropy estimators that require sublinear $n\ll S$ samples. The earliest indication to this effect appeared in Paninski , but only an existential proof based on the Stone–Weierstrass theorem was provided. It was therefore a breakthrough when Valiant and Valiant introduced the first explicit entropy estimator requiring a sublinear number of samples. They showed that $n\gg S/\ln S$ samples are both necessary and sufficient to consistently estimate the entropy of a discrete distribution. However, the entropy estimators based on linear programming proposed in Valiant and Valiant have not been shown to achieve the minimax rate. Another estimator proposed by Valiant and Valiant has only been shown to achieve the minimax rate in the restrictive regime of $\frac{S}{\ln S}\lesssim n\lesssim\frac{S^{1.03}}{\ln S}$. Moreover, the scheme of can only be applied to functionals that are Lipschitz continuous with respect to a Wasserstein metric, which can be roughly understood as those functionals that are equally “smooth” or “smoother” than entropy. Notably, this does not include the functional $F_{\alpha},\alpha<1$ and other interesting nonsmooth functionals of distributions. Also, it is not clear whether these techniques generally lead to minimax rate-optimal estimators. Readers are referred to Valiant’s thesis for more details.

Conceivably, there is a fundamental connection between the smoothness of a functional, and the hardness of estimating it. The ideal solution to this problem would be systematic and capture this trade-off for nearly every functional. Such a comprehensive view of functional estimation has yet to be realized. George Pólya commented that “the more general problem may be easier to solve than the special problem”. This motivated our present work, in which we provide a general framework and procedure for minimax estimation of functionals with non-asymptotic performance guarantees. To make things transparent, let us now start from scratch and demonstrate how our solution has a natural construction.

Suppose we would like to propose a general method to construct minimax rate-optimal estimators for functionals of the form (1). What are the prerequisites that any method must satisfy? Based on our analysis above, the following criteria appear natural:

Asymptotic efficiency. As modern asymptotic statistics tells us, if the function $f(p)$ in (1) is differentiable on $(0,1]$, then the MLE $F(P_{n})$ is asymptotically efficient. In other words, no matter how we adjust the MLE in finite sample settings, we have to ensure that when the number of samples $n$ go to infinity while $S$ remains fixed, our estimator is very similar to the MLE.

Bias reduction. As our analysis of MLE indicates, the MLE usually has large bias and small variance in high dimensions. Hence, the general method has to reduce bias in finite samples.

Let us attempt to understand an estimator’s bias more carefully. In the simplest setting, consider a Binomial random variable $X\sim\mathsf{B}(n,p)$, and suppose we wish to estimate the scalar $f(p)$ based on $X$. Denote by $g(X)$ an arbitrary estimator for $f(p)$. The bias of $g(X)$ can be written as

One may initially be tempted to use the Taylor series to approximate $f(p)$. However, a more careful inspection indicates that the Taylor polynomial is inappropriate for approximating general continuous functions. Even setting aside questions of convergence, Taylor polynomials are not defined for functions that are not infinitely differentiable. Even for functions that are analytic (such as $e^{x},x\in$), one can show that truncating the Taylor series up to order $n$ results in maximum error $\sim\frac{1}{(n+1)!}$ on $$, but there exists a polynomial with order$n$whose maximum approximation error is asymptotically$\frac{1}{2^{n}(n+1)!}$, which is the so called best approximation polynomial. The best polynomial approximation is targeted at computing the polynomial that minimizes the maximum deviation of the polynomial from the function$f(p)$. It is known that for any continuous function on a compact interval, there exists a unique best approximation polynomial for any order. Adopting this rationale, we may try to solve the following problem:

where we seek $g^{*}$ minimizing the maximum value of $|\mathsf{Bias}_{p}(g(X))|$. It gives us the best uniform control of the bias since we do not know $p$ a priori.

Applying advanced tools from approximation theory, Paninski tried the idea mentioned above, which unfortunately did not result in improved estimators. It turns out that this idea, while improving significantly in the bias, results in a blowing up of the variance term. Indeed, the squared bias of the estimator designed above can be shown to be $S^{2}/n^{4}$. Taking $n\gg\sqrt{S}$, the bias term will vanish, but the variance term will diverge, because Paninski already showed that if $n\lesssim S^{1-\delta}$, for any $\delta>0$, the maximum $L_{2}$ risk of any estimators for entropy will be bounded from zero.

In fact, there is a simpler way to understand why the global scale polynomial approximation idea of the form (22) does not work. It is destined to fail because it violates the first prerequisite of any general method to improve MLE in functional estimation. Indeed, this scheme does not behave like the MLE even if $n\gg S$.

This observation leads us to combine some core ideas that finally constitute our scheme. First, one needs to use approximation theory to reduce bias. Second, one cannot do approximation on a global scale (such as $p\in$), but can only approximate the function $f(p)$ locally. Fortunately, the measure concentration phenomenon allows us to do approximation locally. For example, upon observing $\hat{p}=X/n,X\sim\mathsf{B}(n,p)$, we have $\mathsf{Var}(\hat{p})=\frac{p(1-p)}{n}$, which vanishes as $n\to\infty$. Finally, where should we approximate? Intuitively, the bias is mainly due to the set of points where the function $f(p)$ changes abruptly. For $f(p)=-p\ln p$ or $p^{\alpha},\alpha>0$, the most “nonsmooth” point is $p=0$.

To sum up, we need to approximate locally around the “nonsmooth” points to reduce bias. Natural as this statement may seem, there are some parameters to be carefully specified. For this subsection we only consider $f(p)=-p\ln p$ or $p^{\alpha},\alpha>0$. We detail the construction of our scheme by posing the following natural questions:

If we approximate function $f(p)$ in interval $[0,\Delta_{n}]$, how should we choose $\Delta_{n}$?

If we use a polynomial with order $K_{n}$ to approximate $f(p)$ in $[0,\Delta_{n}]$, how should we choose $K_{n}$? What should we do after obtaining the polynomial?

What should we do in interval $p\in[\Delta_{n},1]$?

Let us now answer these questions in the order in which they were asked. The value $\Delta_{n}$ should always be chosen to be the smallest number such that we can localize the parameter $p$. In other words, suppose we observe $X\sim\mathsf{B}(n,p)$. Then, $\Delta_{n}$ should be chosen to ensure that if $X/n\leq c\Delta_{n}$, $c>0$ is a constant, then $p\in[0,\Delta_{n}]$ with high probability. Similarly, if $X/n\geq C\Delta_{n}$, $C>0$ is a constant, then $p\in[\Delta_{n},1]$ with high probability. It turns out that $\Delta_{n}\asymp\frac{\ln n}{n}$ fulfills this goal (cf. Lemma 21).

Regarding the second question, the value $K_{n}$ should always be chosen to be the largest number such that the increased variance does not exceed the bias. Indeed, if we use order $n$ approximation, then we essentially go back to the idea (22) that increases the variance too much such that the resulting estimator does not have vanishing risk. It turns out for $H(P)$ and $F_{\alpha}(P)$, $K_{n}\asymp\ln n$ is the correct order for which we need to conduct the best polynomial approximation. Suppose we have obtained the best polynomial approximation of order $K_{n}$ for function $f(p)$ over regime $[0,\Delta_{n}]$. Noting that any polynomial of order no more than $n$ can be estimated without bias using the estimator in (21), we use the corresponding unbiased estimator to estimate this $K_{n}$-order polynomial, thereby ensuring that the bias of this estimator when $p\in[0,\Delta_{n}]$ is exactly the polynomial approximation error in approximating $f(p)$ over $[0,\Delta_{n}]$.

The third question refers to the scheme in the “smooth” regime. Interestingly, it was already observed in 1969 by Carlton that Miller’s bias correction formula (19) should only be applied when $p_{i}\gg\frac{1}{n}$. In other words, Miller’s formula (19) is relatively accurate when $p_{i}\gg\frac{1}{n}$. In our case, since we have already chosen $\Delta_{n}\asymp\frac{\ln n}{n}$, in the smooth regime we have $p_{i}\gtrsim\frac{\ln n}{n}$. We use the first order bias-correction in this regime inspired by Miller, whose rationale is the following.

For Binomial random variable $X\sim\mathsf{B}(n,p)$, denote the empirical frequency by $\hat{p}=\frac{X}{n}$. Then it follows from Taylor’s theorem that

where $f^{\prime\prime}(p)$ is the second derivative of $f(p)$. We define the first order bias-corrected estimator of $f(\hat{p})$ by

Taking $f(p)=-p\ln p$, $f^{c}(\hat{p})=-\hat{p}\ln\hat{p}+\frac{1-\hat{p}}{2n}$, which is exactly the Miller–Madow bias corrected entropy estimator. Taking $f(p)=p^{\alpha}$, we have the corresponding bias corrected estimator

Figure 2 demonstrates the estimators for $H(P)$ and $F_{\alpha}(P)$ pictorially, where $\hat{p}_{i}=P_{n}(i)$ is the empirical frequency of $i$-th symbol. An important observation is that our estimator naturally satisfies the first prerequisite of any improved method for functional estimation as discussed above. Indeed, as $n\to\infty$, all the observations will fall in the “smooth” regime, and in the smooth regime our estimators are very similar to the MLE, which naturally implies that they are also asymptotically efficient in the sense of Hájek and Le Cam .

We conclude this subsection by comparing any minimax rate-optimal estimator with our estimator. If we consider entropy, Theorem 1 demonstrates that when $n$ is not too large, the risk is dominated by the first term $\frac{S^{2}}{(n\ln n)^{2}}$, which corresponds to the squared bias of our estimator in the “nonsmooth” regime. Further, it is shown in Wu and Yang that in the worst case, the risk contributed by the “nonsmooth” regime (i.e. $[0,\frac{\ln n}{n}]$) is at least of order $\frac{S^{2}}{(n\ln n)^{2}}$. These observations together imply that the gist of any successful scheme should contribute squared bias nearly $\frac{S^{2}}{(n\ln n)^{2}}$. However, the bias always corresponds to a polynomial approximation error, and in the interval $[0,\frac{\ln n}{n}]$, it roughly corresponds to a polynomial with order $\ln n$. The squared bias $\frac{S^{2}}{(n\ln n)^{2}}$ corresponds to a polynomial whose error in approximating $f(p)=-p\ln p$ in $[0,\frac{\ln n}{n}]$ is nearly the same as the best approximation polynomial with order $\ln n$. The theory of strong uniqueness in approximation theory states that any polynomial whose approximation property is close to the best approximation must be close to the best approximation polynomial. Thus, we conclude that any successful scheme must inherently conduct near-best polynomial approximation in the “nonsmooth” regime, which is what we do in our scheme. Similar arguments also explain $F_{\alpha}(P)$ and Theorem 2, 3, and 4.

II-C Related work

The problem of estimating functionals of parameters is one that has been studied extensively in such fields as statistics, information theory, computer science, physics, neuroscience, psychology, and ecology, to name a few. Different communities have focused on different aspects of this general problem, and some seemingly different problems can be recast as functional estimation ones. Below we review some of the core ideas in various communities.

In light of these shortcomings, there seems to be a perception that estimation under finite sample optimality criteria is not amenable to a general mathematical theory . Classical asymptotic theory is usually the refuge. The beautiful theory of Hájek and Le Cam showed that, under mild conditions, there exist systematic methods to construct an asymptotically efficient estimator $\hat{\theta}_{n}$ for the finite dimensional parameter $\theta$, where if $\varphi(\theta)$ is differentiable at $\theta$, $\varphi(\hat{\theta}_{n})$ is also asymptotically efficient for estimating $\varphi(\theta)$ [48, Lemma 8.14]. Furthermore, it is also known that if the functional is non-differentiable, then it is nearly impossible to get an elegant mathematical theory .

The question of estimating functionals of finite dimensional parameters being satisfactorily answered under classical asymptotics, functional estimation in various nonparametric settings has been a strong area of focus since. There are several profound contributions in this area, of which we only mention a few. The most developed theory deals with linear functionals, for example, see . Another well studied situation deals with the case of “smooth” functionals, see , among others. Estimation of non-smooth functionals is an extremely difficult problem, and is still largely open . In particular, the problem of estimating differential entropy $\int-f\ln f$, where $f$ is a density, remains fertile ground for research, cf. . Similar situations are also true for estimating the entropy of a discrete distribution supported on a countably infinite alphabet, cf. .

II-C2 Information theory

In the information theory community, following the seminal work of Shannon , the focus has been on estimating entropy rates of general stationary ergodic processes with fixed (usually small) support (alphabet) sizes. Outside of the favored binary alphabet, printed English contributed the other interesting example of support size $27$ (including the “space”). Cover and King gave an overview of the entropy rate estimation literature until 1978. Soon after the appearance of universal data compression algorithms proposed by Ziv and Lempel , the information theory community started applying these ideas in entropy rate estimation, e.g. Wyner and Ziv , and Kontoyiannis et al. . Verdú provides an overview of universal estimation of information measures until 2005. Jiao et al. constructed a general framework for applying data compression algorithms to establish near-optimal estimators for information rates, with a focus on directed information.

II-C3 Computer science, physics, neuroscience, psychology, ecology, etc

Much of the efforts in computer science, physics, neuroscience, psychology, ecology, and related fields have focused on some special functionals of particular interest. For example, the problem of estimating Shannon entropy $H(P)$ from a finite alphabet source with i.i.d. observations has been investigated extensively. Section II-A and II-B summarize some of the efforts.

II-C4 Modern era: high dimensions and non-asymptotics

The current era of “big data” abounds with applications in which we no longer operate in the asymptotic regime of large sample sizes. This sample scarcity regime necessitates going beyond classical asymptotic analysis and considering finitely many samples in high dimensions. Indeed, the recent successes of finite-blocklength analysis in information theory , and compressed sensing in statistics have demonstrated the benefit of carefully analyzing practical sample sizes. There are also ample recent examples in statistics approaching classical questions from a high dimensional perspective, cf. . The machine learning community has the tradition of favoring non-asymptotic analysis, and usually pose the question of the sample complexity for achieving $\epsilon$ accuracy with $1-\delta$ probability, cf. . The information theoretic counterpart of high dimensional statistics might be the large alphabet setting, with exciting recent advances (cf. ).

With the above as context, our work revisits the framework of functional estimation for finite dimensional models, with a focus on high dimensional and non-asymptotic analysis.

II-D General methodology for functional estimation

We begin by reviewing the existing general approaches to estimation. Maximum likelihood is the most widely used statistical estimation technique, which emerged in modern form 90 years ago in a series of remarkable papers by Fisher . As evidence of its ubiquity, the Google Scholar search query “Maximum Likelihood Estimation” yields approximately $2,570,000$ articles, patents and books. Indeed, in his response to Berkson in 1980, Efron explains the popularity of maximum likelihood:

“The appeal of maximum likelihood stems from its universal applicability, good mathematical properties, by which I refer to the standard asymptotic and exponential family results, and generally good track record as a tool in applied statistics, a record accumulated over fifty years of heavy usage. ”

Over the years, the following folk theorem seems to have been tacitly accepted by applied scientists:

For a finite dimensional parametric estimation problem, it is “good” to employ the MLE.

From the perspective of mathematical statistics, however, maximum likelihood is by no means sacrosanct. As early as in 1930, in his letters to Fisher, Hotelling raised the possibility of the MLE performing poorly . Subsequently, various examples showing that the performance of the MLE can be significantly improved upon have been proposed in the literature, cf. Le Cam for an excellent overview. However, as Stigler [150, Sec. 12] discussed in his 2007 survey, while these early examples created a flurry of excitement, for the most part they were not seen as debilitating to the fundamental theory. Perhaps because these examples did not provide a systematic methodology for improving the MLE.

In 1956, Stein observed that in the Gaussian location model $X\sim\mathcal{N}(\theta,I_{p})$ (where $I_{p}$ is the $p\times p$ identity matrix), the MLE for $\theta$, $\hat{\theta}^{\textrm{MLE}}=X$ is inadmissible [7, Chap. 1] when $p\geq 3$. Later, James and Stein showed that an estimator that appropriately shrinks the MLE towards zero achieves uniformly lower $L_{2}$ risk compared to the risk of the MLE. The shrinkage idea underlying the James–Stein estimator has proven extremely fruitful for statistical methodology, and has motivated further milestone developments in statistics, such as wavelet shrinkage , and compressed sensing .

One interpretation of the shrinkage idea is that, when one desires to estimate a high dimensional parameter, the MLE may have a relatively small bias compared to the variance. Shrinking the MLE introduces an additional bias, but reduces the overall risk by reducing the variance substantially. A natural question now arises: what about situations wherein the bias is the dominating term? Does there exist an analogous methodology for improving over the performance of the MLE in such scenarios? A precedent to this line of questioning can be found in the 1981 Wald Memorial Lecture by Efron entitled “Maximum Likelihood and Decision Theory”:

“…the MLE can be non-optimal if the statistician has one specific estimation problem in mind. Arbitrarily bad counterexamples, along the line of estimating $e^{\theta}$ from $X\sim\mathcal{N}(\theta,1)$, are easy to construct. Nevertheless the MLE has a good reputation, acquired over 60 years of heavy use, for producing reasonable point estimates. Useful general improvements on the MLE, such as robust estimation, and Stein estimation, are all the more impressive for their rarity. ”

For the aforementioned example, Efron argued that the reason the MLE $e^{X}$ may not be a good estimate for $e^{\theta}$, is that it has a large bias. In particular, the statistician may prefer the uniform minimum variance unbiased estimator (UMVUE), $e^{X-\frac{1}{2}}$ to estimate $e^{\theta}$. As we discussed in the presentation of our main results, the bias is usually the dominating term in estimation of functionals of high-dimensional parameters. Notably, the two general improvements of the MLE, namely robust estimation and shrinkage estimation, are not designed to handle functional estimation problems such as the one presented by Efron. Also, as Efron himself observed, the statistician cannot always rely on the UMVUE to save the day, since these are generally very hard to compute, and may not always exist [155, Remark C, Sec. 7]. Thus, there is a need to address, both in scope and methodology, the improvement over the MLE for problems where the bias is the leading term. Such a solution could be considered the dual of the idea of shrinkage, since the trade-off between bias and variance is now reversed, i.e., one might want to sacrifice the variance to reduce the bias.

II-D2 Approximation: dual of shrinkage

Our main results in this paper imply that Theorem 7 is far from true in high-dimensional non-asymptotic settings. Now, we aim to abstract our scheme in estimating functionals of type (1), and distill a general methodology for estimating functionals of parameters of any finite dimensional parametric families.

We propose to conduct the following two-step procedure in estimating $G(\theta)$.

Classify Regime: Compute $\hat{\theta}_{n}$, and declare that we are operating in the “nonsmooth” regime if $\hat{\theta}_{n}$ is “close” enough to $\Theta_{0}$. Note that $G(\theta)$ is not analytic at any $\theta\in\Theta_{0}$. Otherwise declare we are in the “smooth” regime;

If $\hat{\theta}_{n}$ falls in the “smooth” regime, use an estimator “similar” to $G(\hat{\theta}_{n})$ to estimate $G(\theta)$;

If $\hat{\theta}_{n}$ falls in the “nonsmooth” regime, replace the functional $G(\theta)$ in the “nonsmooth” regime by an approximation $G_{\text{appr}}(\theta)$ (another functional) which can be estimated without bias, then apply an unbiased estimator for the functional $G_{\text{appr}}(\theta)$.

II-D3 Details of “Approximation”

While this general recipe appears clean in its description, there are several problem-dependent features that one needs to design carefully – namely

How to determine “nonsmooth” regime? What is the size of it?

What approximation should we choose to approximate $G(\theta)$ in the “nonsmooth” regime?

What does “ ‘similar’ to $G(\hat{\theta}_{n})$” mean precisely? What exactly do we do in the “smooth” regime?

The careful reader may have realized that Questions 1,2, and 3 resemble the questions we asked in Section II-B. Answers to these questions draw on additional problems we investigated beyond these in the present paper.

We should always choose the “nonsmooth” regime to be the smallest regime such that we can still localize the parameter $\theta$. In other words, when we observe $\hat{\theta}_{n}$ in the “nonsmooth” regime, we should be able to infer with high probability that $\theta$ is also in the “nonsmooth” regime. Similarly, we should also be able to localize the parameter in the “smooth” regime. A concrete case would be the following. Say we observe $X\sim\mathsf{B}(n,p)$, and we would like to estimate a functional which is not analytic at $p_{0}=0.2$. How should we define the “nonsmooth” regime? Noting that $\mathsf{Var}(X/n)=\frac{p(1-p)}{n}$, it turns out we can set the “nonsmooth” regime to be $\left[p_{0}-\sqrt{\frac{p_{0}(1-p_{0})\ln n}{n}},p_{0}+\sqrt{\frac{p_{0}(1-p_{0})\ln n}{n}}\right]$ (cf. Lemma 21).

We should always choose an approximation $G_{\text{appr}}(\theta)$ that can be estimated without bias. This requirement leads us to the general theory of unbiased estimation, which was pioneered by Halmos and Kolmogorov . For a comprehensive survey the readers are referred to the monograph by Voinov and Nikulin .

There is a delicate trade-off: the approximation $G_{\text{appr}}(\theta)$ should be estimated without bias, but also should approximate the functional $G(\theta)$ well, and at the same time not incur too much additional variance. These three requirements yield a highly non-trivial interplay between approximation theory and statistics, of which our understanding is as yet incomplete.

where $\mathsf{poly}_{n}$ is the collection of polynomials with order at most $n$ on $A$, is a crucial object in approximation theory as well as our general methodology. Quantifying $E_{n}[f]_{A}$ and obtaining the polynomial that achieves it turned out to be extremely challenging. Remez in 1934 proposed an efficient algorithm for computing the best polynomial approximation, and it was recently implemented and highly optimized in Matlab by the Chebfun team . Regarding the theoretical understanding of $E_{n}[f]_{A}$, de la Vallée-Poussin, Bernstein, Ibragimov, Markov, Kolmogorov and others have made significant contributions, and it is still an active research area. Among others, Bernstein and Ibragimov showed various exact limiting results for some important classes of functions like $|x|^{p}$ and $|x|^{m}\ln|x|^{n}$. For example, we have

The following limit exists for all $p>0$:

where $\Gamma(\cdot)$ denotes the Gamma function.

Regarding bounds on $E_{n}[f]$ for any finite $n$, Korneichuk [166, Chap. 6] provides a comprehensive study. For a comprehensive treatment of modern approximation theory, DeVore and Lorentz , Ditzian and Totik provide excellent references. For the most up-to-date review of polynomial approximation, we refer the readers to Bustamante .

We emphasize that the discussions above refer to approximation in dimension one. The general multivariate case is extremely complicated. Rice wrote:

“The theory of Chebyshev approximation (a.k.a. best approximation) for functions of one real variable has been understood for some time and is quite elegant. For about fifty years attempts have been made to generalize this theory to functions of several variables. These attempts have failed because of the lack of uniqueness of best approximations to functions of more than one variable. ”

We do not know whether in general polynomial approximation can achieve the minimax rates in general settings. Probably other approximation bases need be chosen for certain problems.

Note that we have assumed $G(\theta)$ is analytic in the “smooth” regime. For various statistical models (like Gaussian and Poisson), any analytic functional admits unbiased estimators. We propose to use Taylor series bias correction in the “smooth” regime, where the order of the Taylor series may vary between problems.

II-E Remaining content

The rest of the paper is organized as follows. Section III details the construction of our estimators $\hat{H}$ and $\hat{F}_{\alpha}$ and their analysis. In Section IV we present our general approach for proving minimax lower bounds and apply it to establish Theorems 2 and 4. Section V presents a few experiments demonstrating the practical advantages of our estimators in entropy estimation, mutual information estimation, entropy rate estimation, and learning graphical models. Complete proofs of the remaining theorems and lemmas are provided in the appendices.

III Estimator construction and analysis

Throughout our analysis, we utilize the Poisson sampling model, which is equivalent to having a $S$-dimensional random vector $\mathbf{Z}$ such that each component $Z_{i}$ in $\mathbf{Z}$ has distribution $\mathsf{Poi}(np_{i})$, and all coordinates of $\mathbf{Z}$ are independent. For simplicity of analysis, we conduct the classical “splitting” operation on the Poisson random vector $\mathbf{Z}$, and obtain two independent identically distributed random vectors $\mathbf{X}=[X_{1},X_{2},\ldots,X_{S}]^{T},\mathbf{Y}=[Y_{1},Y_{2},\ldots,Y_{S}]^{T}$, such that each component $X_{i}$ in $\mathbf{X}$ has distribution $\mathsf{Poi}(np_{i}/2)$, and all coordinates in $\mathbf{X}$ are independent. For each coordinate $i$, the splitting process generates a random variable $T_{i}$ such that $T_{i}|\mathbf{Z}\sim\mathsf{B}(Z_{i},1/2)$, and assign $X_{i}=T_{i},Y_{i}=Z_{i}-T_{i}$. All the random variables $\{T_{i}:1\leq i\leq S\}$ are conditionally independent given our observation $\mathbf{Z}$.

For simplicity, we re-define $n/2$ as $n$, and denote

Our estimator $\hat{F}_{\alpha},\alpha>0$, is constructed as follows.

We explain each equation in detail as follows.

Note that $\hat{p}_{i,1}$ and $\hat{p}_{i,2}$ are i.i.d. random variables such that $n\hat{p}_{i,1}\sim\mathsf{Poi}(np_{i})$. We use $\hat{p}_{i,2}$ to determine whether we are operating in the “nonsmooth” regime or not. If $\hat{p}_{i,2}\leq 2\Delta$, we declare we are in the “nonsmooth” regime, and plug in $\hat{p}_{i,1}$ into function $L_{\alpha}(\cdot)$. If $\hat{p}_{i,2}>2\Delta$, we declare we are in the “smooth” regime, and plug in $\hat{p}_{i,1}$ into $U_{\alpha}(\cdot)$.

The coefficients $g_{k,\alpha},0\leq k\leq K$ are coefficients of the best polynomial approximation of $x^{\alpha}$ over $$up to degree$K$, i.e.,

where $\mathsf{poly}_{K}$ denotes the set of algebraic polynomials up to order $K$. Note that in general $g_{k,\alpha}$ depends on $K$, which we do not make explicit for brevity. Lemma 4 shows that for $nX\sim\mathsf{Poi}(np)$,

Thus, we can understand $S_{K,\alpha}(X),nX\sim\mathsf{Poi}(np)$ as a random variable whose expectation is nearly Note that we have removed the constant term from the best polynomial approximation. It is to ensure that we assign zero to symbols we do not see. the best approximation of function $x^{\alpha}$ over $[0,4\Delta]$.

Any reasonable estimator for $p_{i}^{\alpha}$ should be upper bounded by the value one. We cut off $S_{K,\alpha}(x)$ by upper bound $1$, and define the function $L_{\alpha}(x)$, which means “lower part”.

The function $U_{\alpha}(x)$ (standing for “upper part”) is nothing but a product of an interpolation function $I_{n}(x)$The usage of the interpolation function was partially inspired by Valiant and Valiant . and the bias-corrected MLE. The careful reader may note that the bias-corrected MLE is not exactly the same as what we used before in (25). It is because here we are using the Poisson model instead of the Multinomial model. In the Poisson model, the bias correction formula should be modified to

The interpolation function $I_{n}(x)$ is designed to make $U_{\alpha}(x)$ a smooth function on $$. Indeed, when$0<\alpha<1$, were it not for the interpolation function,$U_{\alpha}(x)$would be unbounded for$x$close to zero. Note that$L_{\alpha}(x)$and$U_{\alpha}(x)$are dependent on$n$. We omit this dependence in notation for brevity. The interpolation function$I_{n}(x)$ is defined as follows:

The following lemma characterizes the properties of the function $g(x;a)$ appearing in the definition of $I_{n}(x)$. In particular, it shows that $I_{n}(x)$ is four times continuously differentiable.

For the function $g(x;a)$ on $[0,a]$ defined as follows,

The function $g(x;1)$ is depicted in Figure 3.

Similarly, we define our estimator for entropy $H(P)$ as

The coefficients $\{g_{k,H}\}_{1\leq k\leq K}$ are defined as follows. We first define

Lemma 20 shows that for $nX\sim\mathsf{Poi}(np)$,

is a near-best polynomial approximation for $-p\ln p$ on $[0,4\Delta]$.

III-B Estimator analysis

We demonstrate our analysis techniques via the proof of Theorem 2 and 3, and note that similar techniques allow us to establish Theorem 1.

The next two lemmas show that the estimators $U_{\alpha}(x),U_{H}(x)$ have desirable bias and variance properties when the true probability $p$ is not too small.

Suppose $nX\sim\mathsf{Poi}(np),p\geq\Delta,c_{1}\ln n\geq 1$. For $0<\alpha<3/2$, we have

The following two lemmas characterize the performance of $S_{K,\alpha}(X)$ and $S_{K,H}(X),nX\sim\mathsf{Poi}(np)$ when $p$ is not too large.

If $nX\sim\mathsf{Poi}(np),p\leq 4\Delta,\alpha>0$, we have

and for $n$ large enough, we can take $c_{3}=\frac{2\mu(2\alpha)c_{1}^{\alpha}}{c_{2}^{2\alpha}}$, where $c_{3}$ is the constant appearing in Lemma 19. If we also have $c_{2}\leq 4c_{1}$, then

For the entropy, if $p\leq 4\Delta$, we have

When $n$ is large enough, $C$ can be taken to be $\frac{4c_{1}\nu_{1}(2)}{c_{2}^{2}}$, which is given in Lemma 20. If we also have $c_{2}\leq 4c_{1}$, then

If $nX\sim\mathsf{Poi}(np),p\leq\frac{1}{n\ln n},1<\alpha<3/2$, then for $c_{2}\leq 4c_{1}$,