Inference with Multinomial Data: Why to Weaken the Prior Strength

1 Introduction

In this paper we consider the problem of inference from multinomial data with chances θ = [θ1, . . . , θk]′ . We compare the Maximum Likelihood Estimator (MLE) and the most commonly used Bayesian estimators obtained by assuming a prior Dirichlet distribution with "non-informative" prior parameters such as Laplace, Perks, Jeffreys, and Haldane. Inference in a multinomial-Dirichlet model is a recurrent problem in Artificial Intelligence and Statistics. For instance, it appears in parameter learning of probabilistic graphical models (such as Bayesian networks and some variations) [Koller and Friedman, 2009, Ch. 17], in smoothing methods in information retrieval [Zhai and Lafferty, 2001] and topic models [Mimno and McCallum, 2008], in Bayesian reliability analysis [Somerville et al., 1997], etc.

Consider c1, . . . , ck categories and θj the chance of cj to be observed, for j = 1, . . . , k. Inference about the vector of k parameters θ = [θ1, . . . , θk]′ ∈ Sθ is desired, where Sθ = {θj : 0 ≤ θj ≤ 1 for all j and θ′1 = 1} (1 denotes a row of 1s, i.e., 1 = [1, 1, . . . , 1]′ ∈ Rk ). The observed data consists in a vector of counts n = [n1, n2, . . . , nk]′ , where nj is the number of times in which the j-th category is observed and N = n′1 is the total number of observations.

The goal is thus to estimate the parameter vector θ based on the vector of observations n. Assuming that the probability of observing n, conditionally on θ, can be represented as a multinomial distribution: P (θ, n) =

be obtained by computing the Maximum Likelihood Estimator (MLE), i.e., by maximizing the likelihood P (θ, n) w.r.t. θ subject to the constraint θ′1 = 1, which gives: ˆ θMLE = n/N .

An alternative way is to follow a Bayesian approach. The multinomial is a member of the exponential family and its natural conjugate prior is the Dirichlet distribution. Hence, assuming a Dirichlet prior over θ and applying Bayes' rule to the multinomial-Dirichlet conjugate model, the following posterior density function is obtained:

where D(α, θ) ∝ k∏ j=1 θαj-1 j is the Dirichlet prior with parameters α = [α1, . . . , αk]′ and αj > 0 for j = 1, . . . , k. In the following we introduce the notation used in this paper. We assume αj = stj and s = ∑kj=1 αj, with 0 < t < 1, t′1 = 1, t = [t1, t2, . . . , tk]′ . Notice that s is the strength of the prior information (equivalent sample size or number of pseudo-counts) and tj is the prior mean. The posterior expectation of θ given n is then given by:

which gives a point estimate of θ.

The parameters s and t represent the a-priori information. In case no prior information is available, the common approach is to select these parameters to represent a noninformative prior. The most used non-informative priors select tj = 1/k for j = 1, 2, . . . , k but differ in the choice of the value of s. Bayes and Laplace suggest to use a uniform prior s = k, Perks suggests s = 1, Jeffreys suggests s = k/2, and Haldane suggests s = 0. Nevertheless, the analysis we conduct is general and applies to other choices of s too.

To compare the goodness of different point estimates, a measure of estimation performance must be defined. A popular measure of performance is the matrix mean-squared error (MSE), which is defined as

where the first term of the summation is the "squared-bias" of the estimator and the second term is its variance matrix. 1 Here the unknown parameter vector θ is assumed to be deterministic and, thus, the expectation is only over the data.

In the problem of inference from multinomial data, it is well known that the estimate ˆ θMLE is unbiased, which means that En[ ˆ θ] = θ; and achieves the Cramer-Rao Lower Bound (CRLB) for unbiased estimators, i.e. En[( ˆ θMLE - θ)( ˆ θMLE - θ)′ ] = ΣMLE, where ΣMLE is the inverse of the Fisher information matrix. These facts do not imply that MLE always provides a small MSE, especially for small data samples. In fact, since "MSE=variance + squared bias" and trading-off bias for variance, it is possible to design estimators that yield a lower MSE than the CRLB for unbiased estimators [Ghosh et al., 1983; Stein, 1956].

Since the MSE depends on the unknown θ, it is not obvious how to compare estimators in terms of MSE. However some estimators may be uniformly better than others in terms of MSE, in other words, they can be better for all possible values of θ. For this purpose, we say that an estimator ˆ θ dominates another estimator ˆ θ0 on a convex set Θ if its MSE is never greater than that of ˆ θ0 for all values of θ in Θ, and is strictly smaller for some θ in Θ. An estimator is Θ-admissible if it is not dominated by any other estimator on Θ [Berger, 1985]. Hence, it is reasonable to prefer admissible estimators. In the problem of inference frommultinomial data, it can be shown that the MLE is admissible w.r.t. the MSE criterion if Θ = Sθ [Johnson, 1971]. However, MLE might not be admissible on Θ ⊂ Sθ, because estimators that dominate MLE may exist if a proper subregion of the parameter space is considered.

ˆ

If one assumes the estimator θ = (n + st)/(N + s), obtained by a prior Dirichlet distribution with parameters s and t for θ, then the values of s and t can be designed in order to dominate MLE on Θ. The answer to this question is partially given in [Benavoli and de Campos, 2009], where the authors determine a closed-form solution for the dominance. This solution is employed there with two aims. First, for the binomial case, the authors analyze the performance of Bayesian estimators with t = 1/2 and different choices of s corresponding to the most used non-informative priors. In particular, they determine the region Θ (an interval in the binomial case) where these estimators dominates MLE. Second, assuming that Θ is given as prior information, they derive an ad-hoc criterion to choose an "optimal" value for s and t which guarantees the dominance. In this paper, we generalize the analysis to the multinomial case, where we show that the coverage of the set Θ (that is, the ratio between volume of Θ and the volume of the whole space Sθ) on which the Bayesian estimator of Eq. (2) dominatesMLE decreases at the increasing of the number of categories k. This means that, if s is kept constant, then the region where MLE is preferable to the estimator of Eq. (2) becomes larger with the increasing of k, and soon the MLE becomes the only admissable estimator for any practical scenario. However, this can be avoided if the strength s of prior is not kept constant but decreased with k. As it will be clear by the analysis, we argue that s should decrease at a rate proportional to k, or in other words, each Dirichlet parameter αj should be further corrected by dividing it to k. This corrected version of the Bayesian estimator tends quickly to MLE as k increases, having almost no practical difference already for somewhat small values of k (15 or so), but are still preferable to MLE as they avoid problems with zero counts.

1 The MSE defined in Eq. (3) is a matrix and not a scalar. In the literature, sometimes the MSE is defined as the trace of the matrix in Eq. (3), however in this paper we adopt the matrix definition. The motivation for this choice will be clarified in Section 2.

Before proceeding, we point out that the analysis performed here assumes that no additional information is available to select the prior. It is obvious that better priors can be chosen if a-priori information, for example from domain knowledge or other data source, is available. For instance, this is mostly the case in language processing [Zhai and Lafferty, 2001]. Nevertheless, the argument that prior strength should be reduced with the increase of k might still have to be taken into account, however centering the analysis on the informative prior.

2 MLE-dominating priors

In this section we summarize the results from [Benavoli and de Campos, 2009] that are used in the rest of this paper. Consider an estimator with structure as in Eq. (2). The goal is to choose the free parameters s and t so as to guarantee that: 2

for each vector θ in a convex set Θ. The right-hand side of Ineq. (4) is denoted by ΣMLE = (σij), which represents the covariance matrix of the MLE whose elements are σii = θi(1 - θi)/N and σij = -θiθj/N, for i, j = 1, 2, . . . , k and i ̸= j. The matrix domination considered in Ineq. (4) guarantees a MSE reduction for all the components of the parameter vector θ to be estimated and, thus, is stronger than a trace domination that would only guarantee an improvement for the sum of the MSEs of such components. This is the motivation behind the choice of the matrix MSE instead of the trace MSE that we have followed in this paper.

Manipulating En[(θ - ˆ θ)(θ - ˆ θ)′ ], it can be shown that:

2 Notice that Ineq. (4) is a matrix inequality. For two matrices A,B of compatible dimensions, the inequality A ≤ B means that B - A is nonnegative definite.

and Ineq. (4) becomes

The above inequality is satisfied if and only if

holds for each θ ∈ Θ. Hence, an estimator ˆ θ has MSE lower than that of MLE for all θ ∈ Θ if s and t are chosen according to Ineq. (7). If Θ is a convex polytope of vertices θv1 , θv2 , . . . , θvm , i.e. Θ = Ch{θv1 , θv2 , . . . , θvm } (Ch{·} stands for convex hull), then Ineq. (7) can be further simplified. In fact, in this case, a necessary and sufficient condition for Ineq. (7) to be satisfied for each θ ∈ Θ is to hold on the vertices of the polytope Θ. The MLE-dominance is guaranteed if s and t are chosen such that:

where θvj i denotes the i-th component of the j-th vertex. The above m-inequalities define all the values of s and t which guarantee theMLE-dominance. Using Ineq. (8), the binomial case can be analyzed by taking a Bayesian estimators with t = 1/2 and different choices of s. In particular, the set Θ becomes an interval [ϵ, 1 - ϵ] and Ineq. (8) is satisfied if

Hence, the values of the true θ1 (notice that θ2 = 1 - θ1) for which the MLE-dominance condition is satisfied when N → ∞ are: Haldane (s = 0) needs 0 ≤ θ1 ≤ 1; Jeffreys (s = 0.5) needs 0.05 ≤ θ1 ≤ 0.95; Perks (s = 1) needs 0.1 ≤ θ1 ≤ 0.9; Bayes/Laplace (s = 2) needs 0.15 ≤ θ1 ≤ 0.85; For instance, if s = 6 the Bayesian estimator has a lower MSE than MLE if the true θ is in [0, 25, 0.75], and thus it is preferable in half of the parameter space. Assuming N → ∞ leads to two properties: (i) the choice of s (if one wants to base their choice in this analysis) does not depend on the sample size; (ii) the obtained value of s is a tighter bound than that of using any finite N, which implies that such value is also a feasible choice for any finite N (just slightly smaller than it could be if the finite N was used).

3 Multinomial data

Hereafter we extend the analysis of the end of Section 2 to the multinomial case. In particular, we aim to show that the most used non-informative Bayesian estimators do have a region where they are MLE-dominant, but such region quickly reduces in size with the increase in the number of categories of the multinomial. Before performing this analysis we must define the meaning of "size" of a MLE-dominance region. There are two criteria that are defined and explored here: coverage and fitness. In the following, we assume that θ = [θ1, . . . , θk]′ , with θ′1 = 1 and θ ∈ Θk (from now on we use the superscript k on sets to indicate the dimension on which the set is embedded).

3.1 Fitness of a region

A first way to characterize MLE-dominance regions is fitness. The left-hand side of Ineq. (8) can be seen as a chi-square distributed statistics, with t being the observed frequency and θ the true distribution:

where X 2 (k) is a chi-square with k degrees of freedom. Let Θk be the subset of Sθ k where the Bayesian estimator dominates MLE. Eq. (10) suggests that the "significance" of information encoded by a set Θk ⊆ Sθ k can be evaluated by a chi-square test with k - 1 degrees of freedom. As described in Section 2, it not hard to see that the extremes of Θk will generate the most extreme values of this statistics. Hence, we define the fitness of Θk (w.r.t. prior mean t) by comparing

to the chi-square distribution X 2 (k - 1). If CDFk is the cumulative distribution function of X 2 (k), we have

defined as the fitness measurement of the set Θk w.r.t. t. In fact this is the p-value of observing frequencies t if the true is the farthest extreme of Θk , and therefore a small value of F (Θk , t) indicates that few information is encoded by Θk . For instance, F (Sθ k , t) = 0 and F (Θk , t) = 1 if Θk = {t}.

3.2 Coverage of a region

Another way to characterize the information carried out by Θk is through the ratio between its volume and the volume of the whole parameter space Sθ k . Assuming that these sets are regular (k - 1)-simplices (corresponding to points in dimension k whose coordinates sum one) with side L, their volumes are given by Lk-1 √k 2(k-1)/2 (k-1)! . Considering for instance Sθ k , which has side √2, its volume is:

We thus define the proportion of coverage λ(Θk ) for a set Θk ⊆ Sθ k to be equal to its volume divided by the volume of Sθ k , that is λ(Θk ) = V (Θk )/V (Sθ k ). An interesting property of coverage is that, under the assumption that the true θ has equal probability of being any point within Sθ, λ(Θk ) can be viewed as the chance of θ lying in Θk .

3.3 The ε-contaminated set

In order to analyze the dominance that is implied by the choice of different well-known priors, in the sequel we assume Θk to be an ε-contaminated set defined as follows:

where

ext{Sθ k } = {[θ1, . . . , θk]′ ∈ Sθ k : θi = 1, i ∈ {1, . . . , k}} are the extreme points (vertices) of the simplex. Note that this set is symmetric on Sθ k and has k vertices, namely

The set Θkε has several useful properties for the analysis of MLE-dominance: (i) it is symmetric w.r.t. t = 1k ; (ii) it is a regular simplex of side √2(1 - ε); (iii) its sides are equally distant from the border of the simplex Sθ k and touch it only if ε → 0. This latter property is very important for the MLEdominance analysis. Consider Ineq. (8): if any coordinate θi of θ is zero, then the left-hand side of the inequality will go to infinity (because the denominator is zero and the numerator is approximately 1/k), forbidding the inequality to hold for any s except zero. As the coordinates θi get farther from zero as the left-hand side of the inequality makes it easier to be satisfied.

Considering the fitness and coverage of Θkε , we have:

and

Because of the symmetry of Θkε , Ineq. (8), which has to hold for each vertex of the given set, reduces to:

As larger ε as faster the coverage of Θkε decreases, which implies that the farthest possible θ in the set becomes quickly close to 1/k (the sets are shrinking). Moreover, a quick analysis of Ineq. (13) shows that ε and the strength s of a prior whose MLE-dominance region equals to Θkε are directly correlated. For any given k and N, at the increase of s there is an increase of ε (and vice-versa).

3.4 Analysis of estimators

Figure 1 presents the estimators of Haldane, Jeffreys, Perks and Bayes/Laplace, as well as an estimator with s = 2 (all of them use t = 1k ). We see that the coverages of all estimators (apart Haldane, which gives the same estimate as MLE) quickly drop with the number of categories, meaning that the size (relative to the size of the simplex) of the region where they are preferred quickly approaches zero. At the same pace, their fitness increases, again showing their reduction in size. For these reasons, MLE becomes the preferred estimator (in the sense of better MSE on more than half of the parameter space) already with k = 3 and greater, because the coverage of the Bayesian estimators drastically reduces with k (Figure 1 shows that coverage is already small even for k = 4).

By Eq. (12) we clearly see that the coverage of Θkε considerably decreases when k increases (ε is kept fixed on k - this is equivalent to s kept fixed). Hence, a natural approach to maintain the quality of the estimators is to keep the coverage λ(Θkε) constant over k, which implies that ε (and thus the strength s) has to vary with k. If λ0 denotes the desired coverage, we obtain:

Figure 2 shows the value of εk( 12 ), which keeps a coverage of one half of the parameter space for every k. It also presents the (dashed) curve with the value of εk to keep the fitness measure F (Θkε) constant instead. We see that both curves have similar slopes, indicating that coverage and fitness of Θkε react similarly to the increase of k. At first this is slightly surprising, because chi-square distribution that is used by the fitness measure has a correction for the degrees of freedom of the multinomial. On the other hand, volumes used to compute the coverage have no such adaptive parameter corresponding to the increase in dimensionality. In spite of that, we analyze how fitness and coverage are correlated (Figure 3). We point out that this correlation is almost linear for small values of k, but becomes non-linear with its increase. More over, with the increase of the number of categories, we see that fitness goes faster and faster to zero. Still, both measures lead to similar conclusions w.r.t. the correction that has to be applied to ε (or strength s) when one varies k within practical settings - considerably large k suggests ε → 0 (i.e. s → 0) anyway.

Applying Eq. (14) on Ineq. (13) and assuming that N → ∞ (which does not considerably affect the analysis, as we discuss in the final part of this section), we have

(When λ0 is omitted, then it is assumed that λ0 = 1/2 so as to separate the parameter space in two equal parts.)

The main goal of this study is to devise rules to smartly choose the strength s of the prior distribution. If the estimator under analysis dominates MLE on (at least) half of Sθ k , then such an estimator is preferred to MLE. Hence, using Eq. (15) with λ0 = 1/2, we obtain the maximum values that are admissible for s in settings with different values of k (Table 1). Figure 5 presents the graph of sk for λ0 = 1/2 (and also the dashed curve to keep the fitness measure constant).

We point out that the numerator of Eq. (15) approaches 2 log( 1λ ) when k →∞, and it, together with the denominator (k - 1), justifies the reduction of s by a magnitude of O(k) when k increases. Figure 6 shows that an adjustment of O(k) is enough to make the coverage of Perks and the estimator with s = 2 almost constant, but still insufficient for the estimators of Jeffreys and Bayes/Laplace, which use s = k/2 and s = k, respectively. In fact any estimator using strength s = c/O(k), for any constant c and a properly chosen linear function O(k), will have its coverage kept constant with the increase of k. Therefore our suggestion of using estimators with s = c/k follows.

Finally, if we do not assume N →∞, the same values derived for s are valid, because finite values of N can only help (in the sense that finite N can only increase the upper bound of s and the maximum strength s devised by Eq. (15) would still suffice). Figure 4 shows the actual maximum value of s such that at least half of the parameter space is covered, while N varies between 1 and 30. For k = 2, s can be chosen as high as 24 if the sample size is very small, and converges to 6 as the sample size increases. For greater values of k, the convergence to the limit value as if N → ∞ is much faster, and differences in the maximum admissible value of s only occur for very small values of N (Table 1).

4 Conclusion

This paper discusses the problem of inference from multinomial data and address the problem of choosing the strength of the Dirichlet prior under a MLE-dominance criterion. This approach consists of designing free parameters of the estimator so as to guarantee, for any value of the unknown parameter vector to be estimated, an improvement of the mean-squared error with respect to MLE. Given that the true parametrization is equally probable to be any vector of the parameter space, desirable priors are those that lead to MLEdominance in at least half of the parameter space. We show that non-informative Bayesian estimators have a region of MLE-dominance that shrinks with the increase in the number of categories of the multinomial. After a careful analysis, we devise formulas that suggest how one should select the strength of their prior to avoid such problem. They are are based on the coverage of the parameter space and the fitness of the underlying MLE-dominance region. We conclude that priors must have their strength reduced by a factor proportional to the number of categories of the multinomial, otherwise the MLE becomes a preferred estimator. We emphasize that this regards even the estimators that are already "smoothed" by k, such as Laplace (which would receive s/k2 for each Dirichlet parameter αj). Finally, we point out that if one has additional information and can choose a better prior than the non-informative, then the analysis of this paper does not directly apply. Yet, the additional information could be integrated into the analysis, and the general conclusion would be similar: strengths of priors have to react to the number of categories of the multinomial. As future work, we intend to investigate other measures of quality for the estimators, such as Kullback-Leibler divergence and mean absolute error instead of mean squared error, as well as analyze the case of informative priors.

Acknowledgments

This work has been partially supported by the Computational Life Sciences - Ticino in Rete project and the Swiss NSF grant n. 200020-121785/1. The authors thank the anonymous reviewers and the program chair for their suggestions.