---
abstract: |
    The Dirichlet prior is widely used in estimating discrete distributions and functionals of discrete distributions.  In terms of Shannon entropy estimation, one approach is to plug-in the Dirichlet prior smoothed distribution into the entropy functional, while the other one is to calculate the Bayes estimator for entropy under the Dirichlet prior for squared error, which is the conditional expectation.  We show that in general they do *not* improve over the maximum likelihood estimator, which plugs-in the empirical distribution into the entropy functional.  No matter how we tune the parameters in the Dirichlet prior, this approach cannot achieve the minimax rates in entropy estimation, as recently characterized by Jiao, Venkat, Han, and Weissman [@Jiao--Venkat--Han--Weissman2014minimax], and Wu and Yang [@Wu--Yang2014minimax].  The performance of the minimax rate-optimal estimator with $n$ samples is essentially *at least* as good as that of the Dirichlet smoothed entropy estimators with $n\ln n$ samples. 

    We harness the theory of approximation using positive linear operators for analyzing the bias of plug-in estimators for general functionals under arbitrary statistical models, thereby further consolidating the interplay between these two fields, which was thoroughly developed and exploited by Jiao, Venkat, Han, and Weissman [@Jiao--Venkat--Han--Weissman2014MLEnew].  We establish new results in approximation theory, and apply them to analyze the bias of the Dirichlet prior smoothed plug-in entropy estimator.  This interplay between bias analysis and approximation theory is of relevance and consequence far beyond the specific problem setting in this paper. 
author:
- 'Yanjun Han, , Jiantao Jiao, , and Tsachy Weissman, '
bibliography:
- 'IEEEfull.bib'
- 'strings.bib'
- 'Reference.bib'
- 'refs.bib'
- 'wu\_journal.bib'
- 'spca.bib'
- 'wu\_conf.bib'
- 'refs\_SPCA.bib'
- 'entropy.bib'
- 'wu\_big.bib'
- 'di.bib'
- 'wu\_new.bib'
- 'refs-entropy.bib'
- 'refs-opnorm.bib'
- 'bstcontrol.bib'
title: 'Does Dirichlet Prior Smoothing Solve the Shannon Entropy Estimation Problem?' 
---

entropy estimation, approximation theory using positive linear operators, Dirichlet prior smoothing, minimax optimality, maximum risk, high dimensional statistics, large alphabet

Introduction
============

The Shannon entropy of a discrete distribution, defined by $$\label{eqn.entropy}
H(P) \triangleq \sum_{i = 1}^S p_i \ln \frac{1}{p_i},$$ emerged in Shannon’s 1948 masterpiece [@Shannon1948] in the answers to the most fundamental questions of compression and communication.  In addition to its prominent roles in information theory, Shannon entropy has been widely applied in such disciplines as genetics [@Olsen--Meyer--Bontempi2009impact], image processing [@Pluim--Maintz--Viergever2003mutual], computer vision [@Viola--Wells1997alignment], secrecy [@Batina--Gierlichs--Prouff--Rivain--Standaert--Veyrat2011mutual], ecology [@Hill1973diversity], and physics [@Franchini--Its--Korepin2008Renyi].  In most real-world applications, the true distribution underlying the generation of the data is unknown, making it impossible for us to compute the Shannon entropy exactly.  Hence, in nearly every practical problem involving Shannon entropy, we need to estimate it from data. 

Classical theory is mainly concerned with the case where the number of i.i.d.  samples $n\to \infty$, while the alphabet size $S$ is fixed.  In that scenario, the maximum likelihood estimator (MLE), which plugs in the empirical distribution into the definition of entropy, is *asymptotically efficient* [@Vandervaart2000 Thm.  8.11, Lemma 8.14] in the sense of the Hájek convolution theorem [@Hajek1970characterization] and the Hájek–Le Cam local asymptotic minimax theorem [@Hajek1972local].  In contrast, various modern data-analytic applications deal with datasets which do not fall into the regime of fixed $S$ and $n\to \infty$. In fact, in many applications the alphabet size $S$ is comparable to or even larger than the number of samples $n$. Indeed, the Wikipedia page on the Chinese characters shows that the alphabet of the Chinese language is at least $80,000$ characters large.  As another case in point, half of the words in the Shakespearean canon appeared only once [@Efron--Thisted1976]. 

The problem of entropy estimation in the large alphabet regime (or non-asymptotic analysis) has been investigated extensively in various disciplines.  We refer to [@Jiao--Venkat--Han--Weissman2014minimax] for a thorough review.  One recent breakthrough in this direction came from Valiant and Valiant [@Valiant--Valiant2011], who constructed the first explicit entropy estimator of sample complexity $n \asymp \frac{S}{\ln S}$, which they also proved to be necessary.  On the other hand, it was shown in [@Paninski2003][@Jiao--Venkat--Han--Weissman2014MLEnew] that the MLE requires $n \asymp S$ samples, implying that MLE is strictly sub-optimal in terms of sample complexity. 

However, the aforementioned estimators have not been shown to achieve the minimax $L_2$ rates.  In light of this, Jiao et al.[@Jiao--Venkat--Han--Weissman2014minimax], and Wu and Yang in [@Wu--Yang2014minimax] independently developed schemes based on approximation theory, and obtained the minimax $L_2$ convergence rates for the entropy.  Further, Jiao et al.[@Jiao--Venkat--Han--Weissman2014minimax] proposed a general methodology for estimating functionals, and showed that for a wide class of functionals (including entropy and mutual information), their methodology yields minimax rate-optimal estimators whose performance with $n$ samples is commensurate with that of the MLE with $n\ln n$ samples.  They also obtained minimax $L_2$ rates for estimating a large class of functionals.  On the practical side, Jiao et al.[@Jiao--Venkat--Han--Weissman2014beyond] showed that the minimax rate-optimal estimators introduced in [@Jiao--Venkat--Han--Weissman2014minimax] can lead to consistent and substantial performance boosts in various machine learning applications.  It was argued in [@Jiao--Venkat--Han--Weissman2014minimax] that the “only” approach that can achieve the minimax rates for entropy must either implicitly or explicitly conduct best polynomial approximation as [@Jiao--Venkat--Han--Weissman2014minimax] did.  A question that arises naturally then is whether modifications of the plug-in approach, such as the Dirichlet prior smoothing ideas, can improve and result in something which does achieve the minimax rates.  This paper answers this question negatively. 

Dirichlet smoothing may have two connotations in the context of entropy estimation:

-   [@Schurmann--Grassberger1996entropy; @Schober2013some] One first obtains 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)$.

-   [@Wolpert--Wolf1995; @Holste--Grosse--Herzel1998bayes] One calculates the Bayes estimate for entropy $H(P)$ under Dirichlet prior for squared error.  The estimator is the conditional expectation $\bE[H(P)|\mathbf{X}]$, where $\mathbf{X}$ represents the samples. 

We show in the present paper that neither approach results in improvements over the MLE in the large alphabet regime.  Specifically, these approaches require at least $n\gg S$ to be consistent, while the minimax rate-optimal estimators such as the ones in [@Jiao--Venkat--Han--Weissman2014minimax][@Wu--Yang2014minimax] only need $n \gg \frac{S}{\ln S}$ to achieve consistency. 

A main motivation for the present paper, beyond that discussed above, is to demonstrate the power of *approximation theory using positive linear operators* for bounding the bias of plug-in estimators for functionals of parameters under arbitrary statistical models.  It was first shown in Jiao et al.  [@Jiao--Venkat--Han--Weissman2014MLEnew] that under mild conditions, the problem of bias analysis of plug-in estimators for functionals from arbitrary finite dimensional statistical models is equivalent to approximation theory using positive linear operators, a subfield of approximation theory which has been developing for more than a century.  Applying advanced tools from positive linear operator theory [@Paltanea2004], Jiao et al.  [@Jiao--Venkat--Han--Weissman2014MLEnew] obtained tight non-asymptotic characterizations of maximum $L_2$ risks for MLE in estimating a variety of functionals of probability distributions.  In this paper, we contribute to the general positive linear operator theory [@Paltanea2004], and use the Dirichlet smoothing prior plug-in estimator as an example to demonstrate the efficacy of this general theory in dealing with analysis of the bias in estimation problems.  We believe this connection has far reaching implications beyond analyzing bias in statistical estimation, which itself is an important problem. 

Dirichlet Smoothing
-------------------

The Dirichlet smoothing is widely used in practice to overcome the undersampling problem, i.e., one observes too few samples from a distribution $P$. The Dirichlet distribution with order $S\geq 2$ with parameters $\alpha_1,\ldots,\alpha_S>0$ has a probability density function with respect to Lebesgue measure on the Euclidean space $\mathbb{R}^{S-1}$ given by $$f \left(x_1,\cdots, x_{S}; \alpha_1,\cdots, \alpha_S \right) = \frac{1}{\mathrm{B}(\boldsymbol\alpha)} \prod_{i=1}^S x_i^{\alpha_i - 1}$$ on the open $S-1$-dimensional simplex defined by: $$\begin{aligned}
&x_1, \cdots, x_{S-1} > 0 \\
&x_1 + \cdots + x_{S-1} < 1 \\ 
&x_S = 1 - x_1 - \cdots - x_{S-1} \end{aligned}$$ and zero elsewhere.  The normalizing constant is the multinomial Beta function, which can be expressed in terms of the Gamma function: $$\mathrm{B}(\boldsymbol\alpha) = \frac{\prod_{i=1}^S \Gamma(\alpha_i)}{\Gamma\left(\sum_{i=1}^S \alpha_i\right)},\qquad\boldsymbol{\alpha}=(\alpha_1,\cdots,\alpha_S).$$

Denote by $\mathcal{M}_S$ the set of discrete distributions supported on $S$ elements.  Assuming the unknown distribution $P$ follows prior distribution $P \sim \text{Dir}(\boldsymbol\alpha)$, and we observe a vector $\mathbf{X} = (X_1,X_2,\ldots,X_S)$ with Multinomial distribution $\mathsf{Multi}(n;p_1,p_2,\ldots,p_S)$, then one can show that the posterior distribution $P_{P|\mathbf{X}}$ is also a Dirichlet distribution with parameters $${\boldsymbol\alpha} + \mathbf{X} = \left( \alpha_1 + X_1,\alpha_2 + X_2,\ldots,\alpha_K + X_K \right).$$

Furthermore, the posterior mean (conditional expectation) of $p_i$ given $\mathbf{X}$ is given by [@Lehmann--Casella1998theory Example 5.4.4] $$\delta_i(\mathbf{X}) \triangleq \bE[p_i|\mathbf{X}] = \frac{\alpha_i + X_i}{n + \sum_{i = 1}^S \alpha_i}.$$

The estimator $\delta_i(\mathbf{X})$ is widely used in practice for various choices of $\mathbf{\alpha}$. For example, if $\alpha_i = \frac{\sqrt{n}}{S}$, then the corresponding $\left( \delta_1(\mathbf{X}),\delta_2(\mathbf{X}),\ldots,\delta_S(\mathbf{X})\right)$ is the minimax estimator for $P$ under squared loss [@Lehmann--Casella1998theory Example 5.4.5].  However, it is no longer minimax under other loss functions such as $\ell_1$ loss, which was investigated in [@Han--Jiao--Weissman2014minimaxdistribution]. 

Note that the estimator $\delta_i(\mathbf{X})$ subsumes the MLE $\hat{p}_i = \frac{X_i}{n}$ as a special case, since we can take the limit $\boldsymbol\alpha\to 0$ for $\delta_i(\mathbf{X})$ to obtain MLE. 

The Dirichlet prior smoothed distribution estimate is denoted as $\hat{P}_B$, where $$\hat{P}_B = \frac{n}{n + \sum_{i = 1}^S \alpha_i} P_n + \frac{\sum_{i = 1}^S \alpha_i}{n+\sum_{i = 1}^S \alpha_i}  \frac{\boldsymbol\alpha}{\sum_{i = 1}^S \alpha_i}.$$ Note that the *smoothed* distribution $\hat{P}_B$ can be viewed as a convex combination of the empirical distribution $P_n$ and the *prior* distribution $\frac{\boldsymbol\alpha}{\sum_{i = 1}^S \alpha_i}$. We call the estimator $H(\hat{P}_B)$ the *Dirichlet prior smoothed plug-in estimator*. 

Another way to apply Dirichlet prior in entropy estimation is to compute the Bayes estimator for $H(P)$ under squared error, given that $P$ follows Dirichlet prior.  It is well known that the Bayes estimator under squared error is the conditional expectation.  It was shown in Wolpert and Wolf [@Wolpert--Wolf1995] that $$\label{eqn.bayesestimator}
\hat{H}^{\mathsf{Bayes}} \triangleq \bE[H(P)|\mathbf{X}] = \psi\left( \sum_{i = 1}^S (\alpha_i + X_i)+1\right)  - \sum_{i = 1}^S \left( \frac{\alpha_i + X_i}{\sum_{i = 1}^S (\alpha_i + X_i)} \right) \psi(\alpha_i + X_i + 1),$$ where $\psi(z) \triangleq \frac{\Gamma'(z)}{\Gamma(z)}$ is the digamma function.  We call the estimator $\hat{H}^{\mathsf{Bayes}}$ the *Bayes estimator under Dirichlet prior*. 

Non-asymptotic analysis of $L_2$ risk
-------------------------------------

We adopt the conventional statistical decision theoretic framework [@Wald1950statistical] in analyzing the performance of any entropy estimator $\hat{H} = \hat{H}(\mathbf{X})$. Denote by $\mathcal{M}_S$ all discrete distributions with support size $S$. The $L_2$ risk of estimator $\hat{H}$ is defined as $$\bE_P \left( H(P) - \hat{H} \right)^2,$$ where the expectation is taken with respect to the distribution $P$ that generates the observations used by $\hat{H}$. Apparently, the $L_2$ risk is a function of both the *unknown* distribution $P$ and the estimator $\hat{H}$, and in order to obtain a single score that evaluates how well the estimator $\hat{H}$ in the worst possible case, we may want to adopt the minimax criterion [@Wald1950statistical][@Lehmann--Casella1998theory], and evaluate the *maximum risk* $$\sup_{P\in \mathcal{M}_S}  \bE_P \left( H(P) - \hat{H} \right)^2.$$

The $L_2$ risk can be decomposed into two parts: one is the bias, the other is the variance.  They are defined by $$\begin{aligned}
\mathsf{Bias}(\hat{H}) & \triangleq \bE \hat{H} - H(P) \\
\mathsf{Var}(\hat{H}) & \triangleq \bE \left( \hat{H} - \bE \hat{H} \right)^2.  \end{aligned}$$ We have $$\bE_P \left( H(P) - \hat{H} \right)^2 = \left( \mathsf{Bias}(\hat{H}) \right)^2 + \mathsf{Var}(\hat{H}).$$

The literature on concentration inequalities [@Boucheron--Lugosi--Massart2013] provides us with effective techniques for controlling the variance, and we will show it suffices to apply the Efron–Stein inequality to obtain tight bounds on $\mathsf{Var}(H(\hat{P}_B))$. The focus of this paper is on bias analysis, rather than variance.  How hard can it be to bound the bias? 

The reader may be tempted to use the Taylor expansion to analyze the bias of $H(\hat{P}_B)$. Unfortunately, Taylor expansion cannot give satisfactory results for this problem.  To demonstrate this point simply, let us try to apply Taylor expansion to analyze a special case of $H(\hat{P}_B)$, which is the MLE $H(P_n)$.

Considerable effort has been devoted to understanding the performance of the bias of MLE $H(P_n)$. One of the earliest investigations in this direction is due to Miller [@Miller1955], who applied the Taylor expansion and showed that, for any fixed distribution $P$, $$\label{eqn.miller}
\bE H(P_n) = H(P) - \frac{S-1}{2n} + O\left(\frac{1}{n^2}\right).$$ Equation (\[eqn.miller\]) was later refined by Harris [@Harris1975] using higher order Taylor series expansions to yield $$\label{eqn.harris}
\bE H(P_n) = H(P) - \frac{S-1}{2n} + \frac{1}{12 n^2} \left( 1- \sum_{i = 1}^S \frac{1}{p_i} \right) + O\left( \frac{1}{n^3} \right).$$

Harris’s result reveals an undesirable consequence of the Taylor expansion method: one cannot obtain uniform bounds on the bias of the MLE.  Indeed, the term $\sum_{i = 1}^S \frac{1}{p_i}$ can be arbitrarily large for distributions sufficiently close to the boundary of the simplex of probabilities.  However, it is evident that both $H(P_n)$ and $H(P)$ are bounded above by $\ln S$, since the entropy of any distribution supported on $S$ elements is upper bounded by $\ln S$. Conceivably, for such a distribution $P$ that would make $\sum_{i = 1}^S \frac{1}{p_i}$ very large, we need to compute even higher order Taylor expansions to obtain more accuracy, but even with such efforts we can never obtain a uniform bias bound for all $P$.

Harris’s analysis demonstrates a severe limitation of the Taylor expansion method.  The situation is even worse if we consider general functionals.  Indeed, if we want to analyze the bias of the plug-in estimator $F(P_n)$ for $F(P)$, where $F(P)$ is an arbitrary functional of $P$, which may not be differentiable, Taylor expansion may not even be valid. 

Approximation theory using positive linear operators
----------------------------------------------------

An operator $L$ defined on a linear space of functions, $V$, is called *linear* if $$L(\alpha f + \beta g) = \alpha L(f) + \beta L(g), \quad \text{for all }f,g\in V, \alpha,\beta \in \mathbb{R},$$ and is called *positive*, if $$L(f)\geq 0\quad \text{for all }f\in V, f\geq 0.$$

It was first shown in Jiao et al.  [@Jiao--Venkat--Han--Weissman2014MLEnew] that analyzing the bias of *any* plug-in estimator for functionals of parameters from *any* parametric families can be recast as a problem of approximation theory using positive linear operators [@Paltanea2004], which is a subarea of approximation theory.  Paltanea [@Paltanea2004] provides a comprehensive account of the state-of-the-art theory in this subject. 

For a concrete example, the classical Bernstein operator $B_n(f)$ maps a continuous function $f\in C[0,1]$ to another continuous function $B_n(f)\in C[0,1]$ such that $$B_n(f)(x) = \sum_{j = 0}^n f\left(\frac{j}{n}\right)\binom{n}{j} x^j(1-x)^{n-j},$$ where $C[0,1]$ denotes the space of continuous functions on $[0,1]$. One can easily verify that this operator is positive and linear.  Bernstein in 1912 showed that the function $B_n$ uniformly approximates any continuous function $f\in C[0,1]$. In other words, $$\label{eqn.bernsteincanapproximate}
\sup_{x\in [0,1]}|B_n(f)(x) - f(x)| \to 0,\quad n\to \infty.$$ The function $B_n(f)(x)$ is called Bernstein polynomial.  On the other hand, the operator $B_n(f)$ has a probabilistic interpretation: it can be interpreted as the expectation of random variable $f\left(\frac{X}{n}\right)$ where $X$ follows a Binomial distribution $X\sim \mathsf{B}(n,x)$. In other words, the uniform approximation property (\[eqn.bernsteincanapproximate\]) is equivalent to the fact that $$\sup_{x\in [0,1]} \left| \bE_x f\left(\frac{X}{n}\right) - f(x)\right| = \sup_{x\in [0,1]} \left | \mathsf{Bias}\left(f\left(\frac{X}{n}\right)\right) \right| \to 0,\quad n\to \infty,$$ where the expectation is taken with respect to distribution $X\sim \mathsf{B}(n,x)$. Hence, the bias of the plug-in estimator $f(X/n)$ is precisely the approximation error of the Bernstein polynomial in approximating $f$. More generally, as Paltanea [@Paltanea2004 Remark 1.1.2.]  argued using the Riesz–Markov–Kakutani representation theorem, positive linear operators of continuous functions can be “essentially” interpreted as the expectation of some plug-in estimator.  Specifically, let $X$ be a locally compact Hausdorff space, such as $\mathbb{R}^n$. For any positive linear functional $L$ on $C_c(X)$, which denotes the space of continuous compactly supported real/complex valued functions on $X$, the Riesz–Markov–Kakutani representation theorem implies that there is a unique regular Borel measure $\mu$ on $X$ such that $$L(f) = \int_X f(x)d\mu(x),$$ for all $f\in C_c(X)$. Hence, if we assume the functional $L$ maps the indicator function over $X$ to the constant one, then the measure $\mu$ is a probability measure, implying that $L(f)$ can be interpreted as the expectation of a plug-in estimator that plugs-in an estimate of $x$ with distribution $\mu$ into the continuous function $f$.

The rest of the paper is organized as follows.  Section \[sec.mainresults\] discusses the main results of this paper.  Section \[sec.applyapproximationtheory\] develops new results in approximation theory, and applies them to analyze the bias of the Dirichlet prior smoothed plug-in estimator.  Section \[sec.varianceanalysis\] presents a non-asymptotic upper bound on the variance of the Dirichlet prior smoothed plug-in estimator, and section \[sec.prooflowerbound\] proves the lower bound on the maximum risk of the Dirichlet prior smoothed plug-in estimator as well as the Bayes estimator under Dirichlet prior.  The remaining proofs are deferred to the appendices. 

Main Results {#sec.mainresults}
============

For simplicity, we restrict attention to the case where the parameter $\boldsymbol\alpha$ in the Dirichlet distribution takes the form $(a,a,\ldots,a)$. We remark that for general $\boldsymbol\alpha$ the analysis goes through seamlessly to arrive at similar conclusions. 

In comparison to MLE $H(P_n)$, where $P_n$ is the empirical distribution, the Dirichlet smoothing scheme $H(\hat{P}_B)$ has a disadvantage: it requires the knowledge of the alphabet size $S$ in general.  We define $$\hat{p}_{B,i} = \frac{n \hat{p}_i + a}{n + Sa},$$ and $$p_{B,i} = \bE[\hat{p}_{B,i}] = \frac{np_i + a}{n+Sa}.$$ Apparently, $$\label{eqn. Pbconstruction}
\hat{P}_B = \frac{n}{n+Sa} P_n + \frac{Sa}{n+Sa} \text{Uniform distribution},$$ where $P_n$ stands for the empirical distribution.  Analogously, we have $$P_B = \frac{n}{n+Sa} P + \frac{Sa}{n+Sa} \text{Uniform distribution},$$ where $P$ is the true distribution.  Throughout, we adopt the following standard notations: $a_n\lesssim b_n$ means $\sup_n a_n/b_n <\infty$, $a_n \gtrsim b_n$ means $b_n\lesssim a_n$, $a_n \asymp b_n$ means $a_n\lesssim b_n$ and $a_n \gtrsim b_n$, or equivalently, there exist two universal positive constants $c,C$ such that $$\begin{aligned}
  0<c<\liminf_{n\to\infty} \frac{a_n}{b_n} \le \limsup_{n\to\infty} \frac{a_n}{b_n} < C<\infty.\end{aligned}$$

The main results of this paper can now be stated as follows. 

\[thm.upperbound\] If $n\geq Sa$, then the maximum $L_2$ risk of $H(\hat{P}_B)$ in estimating $H(P)$ is upper bounded as $$\sup_{P \in \mathcal{M}_S} \bE_P \left( H(\hat{P}_B) - H(P)\right)^2 \leq \left( \ln\left(1+\frac{S-1}{n+Sa}\right) + \frac{2Sa}{n+Sa} \ln \left( \frac{n+Sa}{2a}\right) \right)^2 + \frac{2n}{(n+Sa)^2}\left[3+\ln\left(\frac{n+Sa}{a+1}\wedge S\right)\right]^2,$$ where $a\wedge b = \min\{a,b\}$. Here the first term bounds the squared bias, and the second term bounds the variance. 

The following corollary is immediate. 

If $n\gg S$ and $a$ is upper bounded by a constant, then the maximum $L_2$ risk of $H(\hat{P}_B)$ vanishes. 

\[thm.lowerbound\] If $n\ge\max\{15S,Sa\}$, then the maximum $L_2$ risk of $H(\hat{P}_B)$ in estimating $H(P)$ is lower bounded as $$\sup_{P \in \mathcal{M}_S} \bE_P \left( H(\hat{P}_B) - H(P)\right)^2 \geq \frac{1}{2}\left[\frac{(S-3)a}{4(n+Sa)}\ln \left( \frac{n+Sa}{a}\right) + \frac{S-1}{8n} + \frac{S^2}{80n^2} - \frac{1}{48n^2}\right]^2 + c \frac{\ln^2 S}{n},$$ where $c>0$ is a universal constant that does not depend on $\hat{P}_B,P$.

If $n<Sa$, then we have $$\sup_{P \in \mathcal{M}_S} \bE_P \left( H(\hat{P}_B) - H(P)\right)^2 \geq \frac{\ln^2 S}{16}.$$ If $n<15S$, then we have $$\sup_{P \in \mathcal{M}_S} \bE_P \left( H(\hat{P}_B) - H(P)\right)^2 \geq \left[\left( \frac{(S-3)a}{4(n+Sa)}\ln  \left( \frac{n+Sa}{a}\right) + \frac{\lfloor n/15 \rfloor}{8n} - \frac{1}{16n} \right)_+ \right]^2,$$ where $\lfloor x \rfloor$ is the largest integer that does not exceed $x$, and $(x)_+ = \max\{x,0\}$ represents the positive part of $x$.

We have the following corollary. 

If $n \lesssim S$, then the maximum $L_2$ risk of $H(\hat{P}_B)$ is bounded away from zero. 

The case $a = 0$ corresponds to the analysis of MLE, which was conducted in [@Jiao--Venkat--Han--Weissman2014MLEnew].  For the MLE, it was shown in [@Jiao--Venkat--Han--Weissman2014MLEnew] that we have, for any $n$, $$\sup_{P \in \cM_S} \bE_P \left( H(P_n) - H(P) \right)^2 \leq  \left( \ln \left( 1+ \frac{S-1}{n} \right) \right)^2 + \left(\frac{(\ln n)^2}{n}\wedge \frac{2(\ln S + 2)^2}{n}\right).$$ and if $n \geq 15S$, $$\sup_{P \in \mathcal{M}_S} \bE_P \left( H(P_n) - H(P) \right)^2 \geq \frac{1}{2}\left(\frac{S-1}{2n} + \frac{S^2}{20n^2} - \frac{1}{12n^2}\right)^2 + c \frac{\ln^2 S}{n},$$ where $c>0$ is the same universal constant as in Theorem \[thm.lowerbound\]. 

The next theorem presents a lower bound on the maximum risk of the Bayes estimator under Dirichlet prior.  Since we have assumed that all $\alpha_i = a,1\leq i\leq S$, the Bayes estimator under Dirichlet prior is $$\hat{H}^{\mathsf{Bayes}} = \psi(Sa +n + 1) - \sum_{i = 1}^S \frac{a + X_i}{Sa + n} \psi(a + X_i + 1).$$

\[thm.bayeslowerbound\] If $S \geq e(2n+1)$ and $n\geq Sa$, then $$\sup_{P\in \mathcal{M}_S} \bE_P \left( \hat{H}^{\mathsf{Bayes}} - H(P) \right)^2 \geq \left( \ln \left( \frac{S}{e(2n+1)} \right) \right)^2.$$

If $n<Sa$, then $$\sup_{P\in \mathcal{M}_S} \bE_P \left( \hat{H}^{\mathsf{Bayes}} - H(P) \right)^2 \geq \left[ \left (\ln \left( \frac{Sa + n}{e(a+n+1)} \right)\right)_+ \right]^2.$$

Evident from Theorem \[thm.upperbound\], \[thm.lowerbound\], and \[thm.bayeslowerbound\] is the fact that in the best situation (i.e. $a$ not too large), both the Dirichlet prior smoothed plug-in estimator and the Bayes estimator under Dirichlet prior still require at least $n\gg S$ samples to be consistent, which is the same as MLE.  In contrast, the minimax rate-optimal estimator in Jiao et al.  [@Jiao--Venkat--Han--Weissman2014minimax] is consistent if $n\gg \frac{S}{\ln S}$, which is the best possible rate for consistency.  Thus, we can conclude that the Dirichlet smoothing technique does *not* solve the entropy estimation problem.  From an intuitive point of view it is also clear: both the Dirichlet prior smoothed plug-in estimator and the Bayes estimator under Dirichlet prior do not exploit the special properties of the entropy functional $p\ln(1/p)$, i.e. the functional has a nondifferentiable point at $p = 0$. The analysis in [@Jiao--Venkat--Han--Weissman2014minimax] demonstrates that the minimax rate-optimal estimator has to exploit the special structure of the entropy function. 

Approximation theory for bias analysis {#sec.applyapproximationtheory}
======================================

Denote by $e_j,j\in \mathbb{N}_+ \cup \{0\}$, the monomial functions $e_j(y) = y^j,y\in I$. For a linear positive functional $F$, we adopt the following notation $$\begin{aligned}
    B_F(x) = \left|F(e_1)-xF(e_0)\right|,\qquad V_F = F\left((e_1-F(e_1)e_0)^2\right),
  \end{aligned}$$ which represent the “bias” and “variance” of a positive linear functional $F$. Define the first order and second order Ditzian–Totik modulus of smoothness [@Ditzian--Totik1987] by $$\begin{aligned}
    \omega_1^\varphi(f,2h) &\triangleq \sup\left\{|f(u)-f(v)|:u,v\in[0,1],|u-v|\le 2h\varphi\left(\frac{u+v}{2}\right)\right\} \\
    \omega_2^\varphi(f,h) &\triangleq \sup\left\{\left|f(u)-2f\left(\frac{u+v}{2}\right)+f(v)\right|:u,v\in[0,1],|u-v|\le 2h\varphi\left(\frac{u+v}{2}\right)\right\}. 
  \end{aligned}$$

A contribution to approximation theory
--------------------------------------

First we recall the following result, which is a direct corollary of [@Paltanea2004 Thm.  2.5.1]. 

\[lemma\_ori\] If $F: C[0,1]\to \mathbb{R}$ is a linear positive functional and $F(e_0)=1$, then we have $$\begin{aligned}
\label{eq:inequality1}
      |F(f)-f(x)| \le \frac{B_F(x)}{2h_1\varphi(x)}\cdot \omega_1^\varphi(f,2h_1) + \frac{5}{2}\omega_2^\varphi(f,h_1),
    \end{aligned}$$ for all $f\in C[0,1]$ and $0<h_1\le \frac{1}{2}$, where $\varphi(x)=\sqrt{x(1-x)}$ and $h_1=\sqrt{F\left((e_1-xe_0)^2\right)}/\varphi(x)=\sqrt{V_F+(B_F(x))^2}/\varphi(x)$.

We remark that Lemma \[lemma\_ori\] cannot yield the desired result for $f(p)=-p\ln p$ and $$\begin{aligned}
    F(f) = \sum_{k=0}^n f\left(\frac{k+a}{n+Sa}\right)\cdot \binom{n}{k}p^k(1-p)^{n-k}. 
  \end{aligned}$$ Specifically, it is easy to show that $$\begin{aligned}
    B_F(p) = \left|\frac{np+a}{n+Sa} - p\right| = \frac{|1-pS|a}{n+Sa},\qquad V_F = \frac{np(1-p)}{(n+Sa)^2}. 
  \end{aligned}$$ One can also easily show that $\omega_1^\varphi(f,2h) \asymp h$ for $f(x) = -x\ln x$. Regarding the second order Ditzian–Totik modulus of smoothness, the following lemma from [@Jiao--Venkat--Han--Weissman2014MLEnew] gives a precise characterization. 

\[lemma. DT2ndentropy\] For $f(x) = -x\ln x, x\in [0,1]$, the second order Ditzian–Totik modulus of smoothness $\omega_2^\varphi(f,h)$ satisfies $$\omega_2^\varphi(f,h) = \frac{h^2 \ln 4}{1+h^2}, \quad h\leq 1.$$

Hence, when $x\to0$, we conclude that $$\begin{aligned}
    h_1 = \frac{\sqrt{V_F+(B_F(x))^2}}{\varphi(x)} \ge \frac{|B_F(x)|}{\varphi(x)} = \frac{a}{n+Sa}\cdot \frac{|1-xS|}{\sqrt{x(1-x)}}
  \end{aligned}$$ is unbounded as $x\to 0$, thus does not satisfy the condition $h_1\le 1/2$. Thus, we cannot directly use the result of Lemma \[lemma\_ori\]. 

It turns out that the general result in Lemma \[lemma\_ori\] can be strictly improved in a general fashion.  The result is given by the following lemma. 

\[lemma\_mod\] If $F: C[0,1]\to \mathbb{R}$ is a linear positive functional and $F(e_0)=1$, then $$\begin{aligned}
\label{eq:inequality2}
      |F(f)-f(x)| \le \omega_1(f,B_F(x);x) + \frac{5}{2}\omega_2^\varphi(f,h_2)
    \end{aligned}$$ for all $f\in C[0,1]$ and $0<h_2\le \frac{1}{2}$, where $\varphi(x)=\sqrt{x(1-x)}$ and $h_2=\sqrt{V_F}/\varphi(x)$, and $$\begin{aligned}
      \omega_1(f,h;x) \triangleq \sup\left\{|f(u)-f(x)|:u\in[0,1],|u-x|\le h\right\}. 
    \end{aligned}$$

Applying Lemma \[lemma\_ori\] to $x=F(e_1)$ we have $$\begin{aligned}
      |F(f)-f(F(e_1))| \le \frac{5}{2}\omega_2^\varphi(f,h_2)
    \end{aligned}$$ and then (\[eq:inequality2\]) is the direct result of the triangle inequality $|F(f)-f(x)|\le |F(f)-f(F(e_1))| + |f(F(e_1))-f(x)|$.

We show that Lemma \[lemma\_mod\] is indeed stronger than Lemma \[lemma\_ori\].  Firstly, due to $h_1\ge h_2$, we have $\omega_2^\varphi(f,h_2)\le \omega_2^\varphi(f,h_1)$. Second, for $x\le 1/2$, we have $$\begin{aligned}
    \frac{B_F(x)}{2h_1\varphi(x)}\cdot \omega_1^\varphi(f,2h_1)
    &\approx \frac{B_F(x)}{2h_1\varphi(x)}\cdot \sup_{0\le s\le1}2h_1\varphi(s)f'(s)\\
    &\ge B_F(x)\cdot  \sup_{x\le s\le1-x}f'(s)\\
    &\approx \sup_{x\le s\le 1-x}\omega_1(f,B_F(x);s)
  \end{aligned}$$ which is almost the supremum of $\omega_1(f,|F(e_1-xe_0)|;s)$ over $s\in[x,1-x]$ and is no less than the pointwise result $\omega_1(f,|F(e_1-xe_0)|;x)$, and here we have used the inequality $\varphi(s)\ge\varphi(x)$ for $x\le s\le 1-x$. A similar argument also holds for $x>1/2$. Hence, Lemma \[lemma\_mod\] transforms the first order term from the norm result in Lemma \[lemma\_ori\] to a pointwise result, which may exhibit great advantages when we applied it to specific problems. 

Application of the improved general bound to our problem
--------------------------------------------------------

\[thm\_2\] If $n\ge \max\{ Sa,4\}$, then $$\begin{aligned}
  \sup_{P\in\mathcal{M}_S}\bE_P|H(\hat{P}_B) - H(P)| \le  \frac{5nS\ln2}{(n+Sa)^2} + \frac{2Sa}{n+Sa}\ln\left(\frac{n+Sa}{2a}\right).\end{aligned}$$

Note that Theorem \[thm\_2\] implies a slightly weaker bias bound than Theorem \[thm.upperbound\], but it is only sub-optimal up to a multiplicative constant.  The bias bound in Theorem \[thm.upperbound\] is obtained using another technique which is tailored for the entropy function, and is presented in the appendix. 

Now we give the proof of Theorem \[thm\_2\] using the approximation theoretic machinery we just established.  Note that $h_2 = \frac{\sqrt{n}}{n+Sa}$. In order to ensure that $h_2 \leq 1/2$, it suffices to take $n\geq 4$.

In light of Lemma \[lemma\_mod\], we have $$\begin{aligned}
    \bE_P|H(\hat{P}_B) - H(P)| &\le \sum_{i=1}^S \left(\omega_1\left(f,\frac{|1-p_iS|a}{n+Sa};p_i\right) + \frac{5n\ln 2}{(n+Sa)^2}\right)\\
    &\le -\left(\sum_{i=1}^S\frac{|1-p_iS|a}{n+Sa}\right)\ln\left(\frac{1}{S}\sum_{i=1}^S\frac{|1-p_iS|a}{n+Sa}\right) + \frac{5nS\ln2}{(n+Sa)^2}\\
    &\le \frac{2Sa}{n+Sa}\ln\left(\frac{n+Sa}{2a}\right) + \frac{5nS\ln2}{(n+Sa)^2}
  \end{aligned}$$ where we have used the fact that if $|x-y|\leq 1/2, x,y\in [0,1]$, then $|x\ln x - y\ln y|\leq -|x-y|\ln |x-y|$. The readers are referred to the proof of Cover and Thomas [@Cover--Thomas2006 Thm.  17.3.3] for details.  We also utilized the fact that if $n\geq Sa$, then for any $i, 1\leq i\leq S$ $$\begin{aligned}
    \frac{|1-p_iS|a}{n+Sa} \le \frac{Sa}{n+Sa} \le \frac{1}{2}. 
  \end{aligned}$$

Variance Analysis {#sec.varianceanalysis}
=================

The following theorem gives the variance bound in Theorem \[thm.upperbound\]. 

The variance of $H(\hat{P}_B)$ is upper bounded as follows: $$\begin{aligned}
  \mathsf{Var}\left(H(\hat{P}_B)\right) \le \frac{2n}{(n+Sa)^2}\left[3+\ln\left(\frac{n+Sa}{a+1}\wedge S\right)\right]^2.\end{aligned}$$

We recall the following bounded differences inequality first. 

[@Boucheron--Lugosi--Massart2013 Cor.  3.2] If function $f: \cZ^n\to\bR$ has the bounded differences property, i.e., for some non-negative constants $c_1,c_2,\cdots,c_n$, $$\begin{aligned}
      \sup_{z_1,z_2,\cdots,z_n,z_i'\in\cZ} |f(z_1,z_2,\cdots,z_n) - f(z_1,\cdots,z_{i-1},z_i',z_{i+1},\cdots,z_n)| \le c_i, \qquad 1\le i\le n
    \end{aligned}$$ then $$\begin{aligned}
      \mathsf{Var}\left(f(Z_1,Z_2,\cdots,Z_n)\right) \le \frac{1}{4}\sum_{i=1}^n c_i^2,
    \end{aligned}$$ when $Z_1,Z_2,\cdots,Z_n$ are independent random variables. 

In our case, apparently $F(\hat{P}_B)$ is a function of $n$ independent random variables $\{Z_i\}_{1\le i\le n}$ taking values in $\cZ=\{1,2,\cdots,S\}$. Changing one location of the sample would make some symbol with count $j$ to have count $j+1$, and another symbol with count $i$ to have count $i-1$. Then the absolute value of the total change in the functional estimator is $$\begin{aligned}
    \left|f\left(\frac{j+1+a}{n+Sa}\right) - f\left(\frac{j+a}{n+Sa}\right) - f\left(\frac{i+a}{n+Sa}\right) + f\left(\frac{i-1+a}{n+Sa}\right)\right|
    \le 2\max_{1\le k\le n} \left|f\left(\frac{k+a}{n+Sa}\right) - f\left(\frac{k-1+a}{n+Sa}\right)\right|. 
  \end{aligned}$$

In light of the Taylor expansion with integral form residue, we have that for $1\ge x\ge t>0$, $$\begin{aligned}
   (x-t)\ln(x-t) = x\ln x -t (\ln x+1) + \int_x^{x-t} \frac{x-t-u}{u}du
 \end{aligned}$$ so $$\begin{aligned}
\label{eq:entropy_inequality}
   |(x-t)\ln(x-t) - x\ln x| \le t|\ln x+1| + \left|\int_x^{x-t} \frac{x-t}{u}du\right| + t\le t|\ln x+1| + 2t \le t(3-\ln x). 
 \end{aligned}$$ As a result, $$\begin{aligned}
    \max_{1\le k\le n}\left|f\left(\frac{k+a}{n+Sa}\right) - f\left(\frac{k-1+a}{n+Sa}\right)\right|
    &\le \max_{1\le k\le n}\frac{1}{n+Sa}\left(3-\ln\left(\frac{k+a}{n+Sa}\right) \right) \\
    &\le \frac{1}{n+Sa}\left(3+\ln \left(\frac{n+Sa}{a+1}\right)\right). 
  \end{aligned}$$ Hence, the bounded differences inequality shows that $$\begin{aligned}
    \mathsf{Var}\left(H(\hat{P}_B)\right) \le n\max_{2\le k\le n}\left(f\left(\frac{k+a}{n+Sa}\right) - f\left(\frac{k-1+a}{n+Sa}\right)\right)^2\le \frac{n}{(n+Sa)^2}\left(3+\ln\left(\frac{n+Sa}{a+1}\right)\right)^2
  \end{aligned}$$ which completes the proof of the first part. 

For the second inequality, the Efron-Stein inequality gives a general upper bound on the variance. 

[@Boucheron--Lugosi--Massart2013 Thm.  3.1] Let $Z_1,\cdots,Z_n$ be independent random variables and let $f(Z_1,Z_2,\cdots,Z_n)$ be a squared integrable function.  Moreover, if $(Z_1',Z_2',\cdots,Z_n')$ are independent copies of $(Z_1, Z_2,\cdots,Z_n)$ and if we define, for every $i=1,2,\cdots,n$, $$\begin{aligned}
      f_i' = f(Z_1,Z_2,\cdots,Z_{i-1},Z_i',Z_{i+1},\cdots,Z_n)
    \end{aligned}$$ then $$\begin{aligned}
      \mathsf{Var}(f) \le \frac{1}{2}\sum_{i=1}^n \bE(f_i-f_i')^2. 
    \end{aligned}$$

Since $H(\hat{P}_B)=\hat{H}_B(Z_1,\cdots,Z_n)$ is invariant to any permutation of $(Z_1,Z_2,\cdots,Z_n)$, we know that the Efron-Stein inequality implies $$\begin{aligned}
   \mathsf{Var}\left(H(\hat{P}_B)\right) \le \frac{n}{2}\bE\left(\hat{H}_B(Z_1',Z_2,\cdots,Z_n)-\hat{H}_B(Z_1,Z_2,\cdots,Z_n)\right)^2
 \end{aligned}$$ where $Z_1'$ is an i.i.d.  copy of $Z_1$.

Now define $$\begin{aligned}
   X_i = \sum_{j=1}^n \mathbbm{1}(Z_j=i),\qquad 1\le i\le S.
 \end{aligned}$$

For brevity, we denote the $S$-tuple $(X_1,\cdots,X_S)$ as $X_1^S$, and the $n$-tuple $(Z_1,\cdots,Z_n)$ as $Z_1^n$. A specific realization of $(X_1,\cdots,X_S)$ is denoted by $x_1^S=(x_1,\cdots,x_S)$, and a specific realization of $(Z_1,\cdots,Z_n)$ is denoted by $z_1^n=(z_1,\cdots,z_n)$. Then we have $$\begin{aligned}
   &\bE\left(\hat{H}_B(Z_1',Z_2,\cdots,Z_n)-\hat{H}_B(Z_1,Z_2,\cdots,Z_n)\right)^2\\ &=
   \sum_{x_1^S} \bP(X_1^S=x_1^S)\bE\left[\left(\hat{H}_B(Z_1',Z_2,\cdots,Z_n)-\hat{H}_B(Z_1,Z_2,\cdots,Z_n)\right)^2\left|X_1^S=x_1^S\right.\right]. 
 \end{aligned}$$

In light of [@Jiao--Venkat--Han--Weissman2014MLEnew Lemma B.1], we know that the conditional distribution of $Z_1$ conditioned on $(X_1,\cdots,X_S)$ is the discrete distribution $(X_1/n,X_2/n,\cdots,X_S/n)$. Denoting $r(p)=f(\frac{np+a}{n+Sa})$, we can rewrite $$\begin{aligned}
   \hat{H}_B(Z_1',Z_2,\cdots,Z_n)-\hat{H}_B(Z_1,Z_2,\cdots,Z_n) = D_- + D_+
 \end{aligned}$$ where $$\begin{aligned}
   D_- &= r\left(\frac{X_{Z_1}-1}{n}\right) - r\left(\frac{X_{Z_1}}{n}\right)\\
   D_+ &= \begin{cases}
     r\left(\frac{X_{Z_1'}+1}{n}\right) - r\left(\frac{X_{Z_1'}}{n}\right)&Z_1\neq Z_1'\\
     r\left(\frac{X_{Z_1'}}{n}\right) - r\left(\frac{X_{Z_1'}-1}{n}\right)&Z_1= Z_1'
   \end{cases}. 
 \end{aligned}$$

Here, $D_-$ is the change in $\hat{H}_B$ that occurs when $Z_1$ is removed according to the distribution $(X_1/n,X_2/n,\cdots,X_S/n)$, and $D_+$ is the change in $\hat{H}$ that occurs when $Z_1'$ is added back according to the true distribution $P$. Now we have $$\begin{aligned}
   \bE[D_-^2|X_1^S] &= \sum_{i=1}^S \frac{X_i}{n}\left(r\left(\frac{X_{i}-1}{n}\right)-r\left(\frac{X_{i}}{n}\right)\right)^2\\
   \bE[D_+^2|X_1^S] &= \sum_{i=1}^S p_i\frac{X_i}{n}\left(r\left(\frac{X_{i}}{n}\right)-r\left(\frac{X_{i}-1}{n}\right)\right)^2
   + \sum_{i=1}^S p_i\left(1-\frac{X_i}{n}\right)\left(r\left(\frac{X_{i}+1}{n}\right)-r\left(\frac{X_{i}}{n}\right)\right)^2
 \end{aligned}$$ where we define $r(x)=0$ when $x\notin[0,1]$. Then, by the law of iterated expectation, we know that $$\begin{aligned}
   \bE[D_-^2] &= \sum_{i=1}^S \sum_{j=0}^n\frac{j}{n}\left(r\left(\frac{j-1}{n}\right)-r\left(\frac{j}{n}\right)\right)^2\times\bP(\mathsf{B}(n,p_i)=j)\\
   \bE[D_+^2] &= \sum_{i=1}^S \sum_{j=0}^n\left(\frac{j}{n}\left(r\left(\frac{j-1}{n}\right)-r\left(\frac{j}{n}\right)\right)^2 + \left(1-\frac{j}{n}\right) \left(r\left(\frac{j+1}{n}\right)-r\left(\frac{j}{n}\right)\right)^2 \right)\times p_i\bP(\mathsf{B}(n,p_i)=j)
 \end{aligned}$$

After some algebra we can show that $\bE[D_-^2]=\bE[D_+^2]$. Hence, we have $$\begin{aligned}
   \mathsf{Var}\left(H(\hat{P}_B)\right) &\le \frac{n}{2}\bE\left(D_-+D_+\right)^2\le n\bE\left(D_-^2+D_+^2\right) = 2n\bE D_-^2 \\
   &= 2n\sum_{i=1}^S\bE P_n(i)\left(r(P_n(i)-\frac{1}{n})-r(P_n(i))\right)^2\\
   &\le 2n\sum_{i=1}^S \bE P_n(i)\left(\frac{1}{n+Sa}\left[3-\ln \left(\frac{nP_n(i)+a}{n+Sa}\right)\right]\right)^2\\
   &= \frac{2n}{(n+Sa)^2}\sum_{i=1}^S \bE P_n(i)\left(3-\ln \left(\frac{nP_n(i)+a}{n+Sa}\right)\right)^2\\
   &\le \frac{2n}{(n+Sa)^2}\sum_{i=1}^S p_i\left(3-\ln \left(\frac{np_i+a}{n+Sa}\right)\right)^2\\
   &\le \frac{2n}{(n+Sa)^2}S\cdot \frac{1}{S}\left(3-\ln \left(\frac{n/S+a}{n+Sa}\right)\right)^2\\
   &= \frac{2n(3+\ln S)^2}{(n+Sa)^2}
 \end{aligned}$$ where we have used the inequality (\[eq:entropy\_inequality\]) and Jensen’s inequality due to $$\begin{aligned}
   \frac{d^2}{dx^2}\left[x\left(\ln \left(\frac{nx+a}{n+Sa}\right)-3\right)^2\right] = \frac{n}{nx+a}\left[3\left(\ln \left(\frac{nx+a}{n+Sa}\right)-3\right)+\frac{nx}{nx+a}\left(4-\ln \left(\frac{nx+a}{n+Sa}\right)\right)\right]<0. 
 \end{aligned}$$

Proof of the lower bound {#sec.prooflowerbound}
========================

Theorem \[thm.lowerbound\] follows upon applying the inequality $$\sup_{P\in\mathcal{M}_S} \bE_P\left(H(\hat{P}_B)-H(P)\right)^2 \geq \max\{C,D\}\geq \frac{C+D}{2},$$ where expressions $C$ and $D$ are given, respectively, in the following two theorems. 

[@Wu--Yang2014minimax] There exists a universal constant $c>0$ such that $$\begin{aligned}
    \sup_{P\in\mathcal{M}_S} \bE_P\left(H(\hat{P}_B)-H(P)\right)^2 \ge c\frac{\ln^2 S}{n}. 
  \end{aligned}$$

If $n\ge\max\{15S,Sa\}$, $$\sup_{P\in\mathcal{M}_S}\left |\bE_P H(\hat{P}_B) - H(P) \right| \ge \frac{(S-3)a}{4(n+Sa)}\ln \left( \frac{n+Sa}{a}\right) + \frac{S-1}{8n} + \frac{S^2}{80n^2} - \frac{1}{48n^2}$$ If $n<Sa$, then $$\sup_{P\in\mathcal{M}_S}\left |\bE_P H(\hat{P}_B) - H(P) \right| \geq  \frac{\ln S}{4}.$$ If $n<15S$, then $$\sup_{P\in\mathcal{M}_S}\left |\bE_P H(\hat{P}_B) - H(P) \right| \geq \frac{(S-3)a}{4(n+Sa)}\ln \left( \frac{n+Sa}{a}\right) + \frac{\lfloor n/15 \rfloor}{8n} - \frac{1}{16n}.$$

By setting $P=(1,0,0,\cdots,0)$, we have $H(P)=0$ and $$\begin{aligned}
    H(\hat{P}_B) = -\frac{(S-1)a}{n+Sa}\ln\left(\frac{a}{n+Sa}\right) -\frac{n+a}{n+Sa}\ln\left(\frac{n+a}{n+Sa}\right) \ge \frac{(S-1)a}{n+Sa}\ln\left(\frac{n+Sa}{a}\right)\quad a.s. 
  \end{aligned}$$ hence we have obtained the first lower bound $$\begin{aligned}
  \sup_{P\in\mathcal{M}_S}\left|\bE_P H(\hat{P}_B) - H(P)\right| \ge \frac{(S-1)a}{n+Sa}\ln\left(\frac{n+Sa}{a}\right). 
  \end{aligned}$$

Note that if $n<Sa$, then $$\sup_{P\in\mathcal{M}_S}\left|\bE_P H(\hat{P}_B) - H(P)\right| \ge \frac{(S-1)a}{2Sa}\ln S \geq \frac{\ln S}{4}.$$ From now on we assume $n\geq Sa$. For $n\ge 15S$, it follows from [@Jiao--Venkat--Han--Weissman2014MLEnew Lemma 2.11] that $$\begin{aligned}
 \label{eqn.mlepaperlowerbound}
    \sup_{P\in\mathcal{M}_S}\left|\bE_P H(\hat{P}) - H(P)\right| \ge \frac{S-1}{2n} + \frac{S^2}{20n^2} - \frac{1}{12n^2}. 
  \end{aligned}$$ If $n<15S$, then it follows from the proof of [@Jiao--Venkat--Han--Weissman2014MLEnew Lemma 2.11] that one can essentially take $S = \lfloor n/15 \rfloor$ in (\[eqn.mlepaperlowerbound\]), and obtain $$\sup_{P\in\mathcal{M}_S}\left|\bE_P H(\hat{P}) - H(P)\right| \ge \frac{\lfloor n/15 \rfloor}{2n} - \frac{1}{4n}.$$

It follows from a refinement result of Cover and Thomas [@Cover--Thomas2006 Thm.  17.3.3] that when $| \hat{p}_{B,i} - \hat{p}_i | \leq 1/2$ for all $i$ (which is ensured by condition $n\geq Sa$), we have $$\begin{aligned}
  |H(\hat{P}_B) - H(\hat{P})| \leq -\frac{2Sa}{n+Sa} \ln \left( \frac{2a}{n+Sa}\right) \le \frac{2Sa}{n+Sa}\ln\left(\frac{n+Sa}{a}\right).\end{aligned}$$ A combination of these two inequalities yield the second lower bound $$\begin{aligned}
  \sup_{P\in\mathcal{M}_S}\left|\bE_P H(\hat{P}_B) - H(P)\right | \ge  \frac{S-1}{2n} + \frac{S^2}{20n^2} - \frac{1}{12n^2} - \frac{2Sa}{n+Sa}\ln\left(\frac{n+Sa}{a}\right)\end{aligned}$$ when $n\geq 15S$, and the second lower bound $$\sup_{P\in\mathcal{M}_S}\left|\bE_P H(\hat{P}_B) - H(P)\right | \ge  \frac{\lfloor n/15 \rfloor}{2n} - \frac{1}{4n} - \frac{2Sa}{n+Sa}\ln\left(\frac{n+Sa}{a}\right)$$ when $n<15S$.

Hence we are done by using these two lower bounds and the inequality $\max\{a,b\}\ge \frac{3a+b}{4}$.

We prove Theorem \[thm.bayeslowerbound\] below. 

We first assume that $n\geq Sa$. Applying the recursive formula $\psi(x+1) = \psi(x) + 1/x$, we know that $\psi(a + X_i + 1)\geq \psi(a+1)$, since $X_i$ is an non-negative integer.  Hence, we have $$\hat{H}^{\mathsf{Bayes}} \leq \psi(Sa + n + 1) - \sum_{i = 1}^S \frac{a+X_i}{Sa + n}\psi(a+1) = \psi(Sa + n+1) - \psi(a+1).$$

We know that the digamma function $\psi(x)$ is strictly increasing from $-\infty$ to $\infty$ on $(0,\infty)$, and that for any $x\geq 1, x\in \mathbb{R}$, $$\ln \lfloor x \rfloor -\gamma \leq \psi(x) \leq 1+ \ln(\lfloor x \rfloor) -\gamma,$$ where $\gamma \approx 0.57721$ is the Euler–-Mascheroni constant. 

Hence, we know $$\hat{H}^{\mathsf{Bayes}} \leq 1-\gamma + \ln(Sa + n + 1) - (-\gamma) = 1 + \ln(Sa + n +1) \leq 1+\ln(2n+1),$$ where we used the assumption that $n\geq Sa$. Since the Bayes estimator under Dirichlet prior is upper bounded by $1+\ln(2n+1)$ for all possible realizations, the squared error it incurs for uniform distribution is at least $\left( \ln S - 1 - \ln(2n+1)\right)^2$ when $S\geq e(2n+1)$. The first part is proved. 

Regarding the second part, we assume $n<Sa$. Take $P = (1,0,\ldots,0)\in \mathcal{M}_S$. Thus, $$\hat{H}^{\mathsf{Bayes}} = \psi(Sa + n + 1) - \frac{(S-1)a}{Sa + n}\psi(a + 1) - \frac{a+n}{Sa + n}\psi(a+n+1),\quad a.s.$$ and the corresponding true entropy $H(P) = 0$.

Using the monotonicity of $\psi(x)$ on $(0,\infty)$, we have $$\begin{aligned}
\hat{H}^{\mathsf{Bayes}} & \geq \psi(Sa + n + 1) - \psi(a+n+1) \\
& \geq \ln(Sa + n) - \gamma - (1+\ln(a+n+1)-\gamma)\\
& = \ln(Sa + n) -1 -\ln(a+n+1) \\
& = \ln \left( \frac{Sa+n}{e(a+n+1)} \right).  \end{aligned}$$

Hence, the squared error incurred by the Bayes estimator under Dirichlet prior for $P = (1,0,\ldots,0)\in \mathcal{M}_S$ is at least $\left( \ln \left( \frac{Sa + n}{e(a+n+1)} \right)\right)_+^2$. The second part is proved. 

Acknowledgments
===============

The authors would like to thank Maya Gupta for raising the question answered in this paper. 

Proof of the bias analysis in Theorem \[thm.upperbound\]
========================================================

The following theorem gives the upper bound on the squared bias used in Theorem \[thm.upperbound\]. 

If $n\ge Sa$, $$\begin{aligned}
  \sup_{P\in\mathcal{M}_S}\bE_P|H(\hat{P}_B) - H(P)| \le \ln\left(1+\frac{S-1}{n+Sa}\right) + \frac{2Sa}{n+Sa} \ln \left( \frac{n+Sa}{2a}\right).\end{aligned}$$

We have $$H(\hat{P}_B) = \sum_{i = 1}^S-\hat{p}_{B,i} \ln \hat{p}_{B,i} = H(P_B) + \sum_{i =1}^S (p_{B,i} - \hat{p}_{B,i}) \ln p_{B,i} - \sum_{i = 1}^S \hat{p}_{B,i} \ln \frac{\hat{p}_{B,i}}{p_{B,i}}.$$

Taking expectations on both sides, we have $$\bE H(\hat{P}_B) - H(P) = H(P_B) - H(P) - \bE D(\hat{P}_B \| P_B),$$ where $D(P\|Q) = \sum_{i = 1}^S p_i \ln \frac{p_i}{q_i}$ is the KL divergence between distributions $P$ and $Q$. In order to analyze the bias, it suffices to analyze the two terms separately.  We first analyze $\bE D(\hat{P}_B \| P_B)$.

It follows from Jensen’s inequality that $$D(P \| Q) = \sum_{i = 1}^S p_i \ln \frac{p_i}{q_i} \leq \ln \left( \sum_{i = 1}^S \frac{p_i^2}{q_i} \right),$$ whose derivation here follows from Tsybakov [@Tsybakov2008 Lemma 2.7]. 

By Jensen’s inequality, we have $$\bE D(\hat{P}_B \| P_B) \leq \bE \ln \left( \sum_{i = 1}^S \frac{\hat{p}_{B,i}^2}{p_{B,i}} \right) \leq \ln \left( \sum_{i = 1}^S \frac{\bE \hat{p}_{B,i}^2 }{p_{B,i}} \right).$$

We also have $$\sum_{i = 1}^S \frac{p_i^2}{q_i} = 1 + \sum_{i =1}^S \frac{(p_i - q_i)^2}{q_i},$$ and that $$\bE (\hat{p}_{B,i} - p_{B,i})^2 = \frac{n^2}{(n+Sa)^2} \bE (\hat{p}_i - p_i)^2 = \frac{n p_i(1-p_i)}{(n+Sa)^2}.$$

Hence, $$\bE D(\hat{P}_B \| P_B) \leq \ln \left( 1+ \sum_{i = 1}^S \frac{n p_i(1-p_i)}{(n+Sa)(np_i+a)} \right) = \ln \left( 1+ \sum_{i = 1}^S \frac{n p_i(1-p_i)}{(n+Sa)np_i} \frac{np_i}{np_i+a} \right) \leq \ln \left( 1+ \sum_{i = 1}^S \frac{1-p_i}{(n+Sa)}  \right),$$ which implies that $$\bE D(\hat{P}_B \| P_B)  \leq \ln \left( 1+\frac{S-1}{n+Sa} \right) \leq \frac{S-1}{n+Sa}.$$

Now we consider the deterministic gap $H(P_B) - H(P)$. It follows from a refinement result of Cover and Thomas [@Cover--Thomas2006 Thm.  17.3.3] that when $| p_{B,i} - p_i | \leq 1/2$ for all $i$, we have $$\label{eq:near_lips}
|H(P_B) - H(P)| \leq \| P_B - P \|_1 \ln \frac{\| P_B - P \|_1}{S}.$$

Note that the condition $n\ge Sa$ ensures that $| p_{B,i} - p_i | \leq 1/2$. We compute $$\| P_B - P\|_1 = \frac{Sa}{n+Sa} \sum_{i = 1}^S \left|p_i -\frac{1}{S}\right|.$$ Taking a naive bound $\sum_{i = 1}^S |p_i -\frac{1}{S}| \leq 2$, we have $$\| P_B - P\|_1  \leq \frac{2Sa}{n+Sa}.$$ Hence, $$|H(P_B) - H(P)| \leq -\frac{2Sa}{n+Sa} \ln \left( \frac{2a}{n+Sa}\right).$$
