Estimating the Mutual Information between two Discrete, Asymmetric Variables with Limited Samples

TL;DR

This paper introduces a new Bayesian estimator for mutual information between two discrete variables, especially effective in limited sample scenarios where one variable is well-sampled but the other is not, reducing bias significantly.

Contribution

The authors propose a novel mutual information estimator that performs well with severely undersampled data, focusing on the inhomogeneity of the conditional distribution of the low-entropy variable.

Findings

Outperforms previous methods in low-sample regimes

Achieves very low bias even with few coincidences

Effective when one variable is well-sampled and the other is not

Abstract

Determining the strength of non-linear statistical dependencies between two variables is a crucial matter in many research fields. The established measure for quantifying such relations is the mutual information. However, estimating mutual information from limited samples is a challenging task. Since the mutual information is the difference of two entropies, the existing Bayesian estimators of entropy may be used to estimate information. This procedure, however, is still biased in the severely under-sampled regime. Here we propose an alternative estimator that is applicable to those cases in which the marginal distribution of one of the two variables---the one with minimal entropy---is well sampled. The other variable, as well as the joint and conditional distributions, can be severely undersampled. We obtain an estimator that presents very low bias, outperforming previous methods even when the sampled data contain few coincidences. As with other Bayesian estimators, our proposal focuses on the strength of the interaction between two discrete variables, without seeking to model the specific way in which the variables are related. A distinctive property of our method is that the main data statistics determining the amount of mutual information is the inhomogeneity of the conditional distribution of the low-entropy variable in those states (typically few) in which the large-entropy variable registers coincidences.

Equations

$$H(Y)=-\sum _y{q}_y\log q_y,$$

$$H(Y|X)=\sum _x{q}_x\left[-\sum _y{q}_{y|x}\log q_{y|x}\right]=\sum _x{q}_{x}H(Y|x).$$

$$I(X,Y)=H(X)+H(Y)-H(X,Y)=H(Y)-H(Y|X).$$

$$\begin{array}[]{c}\langle H|\mathbf{n}\rangle=\int{\rm d}\mathbf{q}\ H(\mathbf{q})\ p(\mathbf{q}|\mathbf{n})=[p(\mathbf{n})]^{-1}\int{\rm d}\mathbf{q}\ H(\mathbf{q})\ p(\mathbf{n}|\mathbf{q})\ p(\mathbf{q}).\end{array}$$

$$p(\mathbf{n}|\mathbf{q})=N!\prod _x\frac{q_x^{n_x}}{n_x!},$$

$$p(\mathbf{n})=\int \mathrm{d}\mathbf{q}p(\mathbf{n}|\mathbf{q})p(\mathbf{q}),$$

$$p(\mathbf{q})=\int \mathrm{d}\beta p(\beta )p(\mathbf{q}|\beta ).$$

$$⟨H|\mathbf{n}⟩=\frac{1}{p(\mathbf{n})}\int \mathrm{d}\beta p(\beta )W(\beta |\mathbf{n}),$$

$$W(\beta |\mathbf{n})=\int \mathrm{d}\mathbf{q}H(\mathbf{q})p(\mathbf{n}|\mathbf{q})p(\mathbf{q}|\beta ).$$

$$I(X;Y)=\sum _x{q}_x\sum _y{q}_{y|x}\log \frac{q_{y|x}}{q_y}=\sum _x{q}_x{D}_{\text{KL}}({\mathbf{q}}_{y|x}||{\mathbf{q}}_y),$$

$${\hat{q}}_y=\frac{n_y}{N}\approx q_y,\forall y.$$

$$p(\{{\mathbf{q}}_{y|x}\}|\beta )$$

$$W(\beta |\mathbf{n})$$

$$p(\beta |\mathbf{n})$$

$$p(\mathbf{n}|{\hat{\mathbf{q}}}_x,\{{\mathbf{q}}_{y|x}\})=N!\prod _{xy}\frac{[{\hat{q}}_x{q}_{y|x}{]}^{n_{xy}}}{n_{xy}!}.$$

$$⟨I|\mathbf{n}⟩(\beta )=\hat{H}(Y)-\sum _x\frac{n_x}{N}\left[{\psi }_0(\beta +n_x+1)-\sum _y\frac{\beta \widehat{q_y}+n_{xy}}{\beta +n_x}{\psi }_0(\beta \widehat{q_y}+n_{xy}+1)\right],$$

$$p(q_{1|x}|\beta )=\frac{\Gamma (\beta )}{\Gamma (\beta /2{)}^2}{\left[q_{1|x}(1-q_{1|x})\right]}^{\beta /2-1}.$$

$$\begin{array}[]{rl}\langle I\rangle(\beta)&=\log 2-\psi_{0}(\beta+1)+\psi_{0}(\beta/2+1).\end{array}$$

$$p(\log \beta )=\frac{\beta /2}{\log 2}|2{\psi }_1(\beta +1)-{\psi }_1(\beta /2+1)|.$$

$$⟨I|\mathbf{n}⟩(\beta )=\hat{H}(Y)-\sum _x\frac{n_x}{N}\left[{\psi }_0(\beta +n_x+1)-\sum _{y\in \{0,1\}}\left(\frac{\beta /2+n_{xy}}{\beta +n_x}\right){\psi }_0(\beta /2+n_{xy}+1)\right].$$

$$p(\mathbf{n}|\beta )$$

$$\begin{array}[]{rl}\log p(\mathbf{n}|\beta)&\displaystyle=\sum_{n\geq n^{\prime}}m_{nn^{\prime}}\log\left[\frac{\Gamma(n+\beta/2)\Gamma(n^{\prime}+\beta/2)\Gamma(\beta)}{\Gamma(n+n^{\prime}+\beta)\Gamma(\beta/2)^{2}}\right]\\ \\ &\displaystyle=\sum_{n\geq n^{\prime}}m_{nn^{\prime}}\log p_{nn^{\prime}}(\beta),\end{array}$$

$$\begin{array}[]{rl}p_{10}(\beta)&\displaystyle=\frac{\beta/2}{\beta}=\frac{1}{2}\\ \\ p_{11}(\beta)&\displaystyle=\frac{(\beta/2)^{2}}{\beta(\beta+1)}=\frac{\beta}{4(\beta+1)}\\ \\ p_{20}(\beta)&\displaystyle=\frac{(\beta/2)(\beta/2+1)}{\beta(\beta+1)}=\frac{(\beta/2+1)}{2(\beta+1)}\\ \dots\\ p_{nn^{\prime}}(\beta)&\displaystyle=\frac{(\beta/2)(\beta/2+1)\dots(\beta/2+n-1)\,(\beta/2)(\beta/2+1)\dots(\beta/2+n^{\prime}-1)}{\beta(\beta+1)\dots(\beta+n+n^{\prime}-1)}.\end{array}$$

$$\frac{\partial }{\partial \beta }\log p(\mathbf{n}|\beta )=\frac{m_{11}}{\beta }+\frac{m_{20}}{\beta +2}-\frac{m_{11}+m_{20}}{\beta +1}=0.$$

$$⟨H_{Y|X}⟩(\beta )=\sum _k{f}_k\sum _{n+n^{\prime }=k}{f}_{n{n}^{\prime }}{H}_{n{n}^{\prime }}(\beta )$$

$$H_{n{n}^{\prime }}(\beta )={\psi }_0(n+n^{\prime }+\beta +1)-\frac{(n+\beta /2){\psi }_0(n+\beta /2+1)+(n^{\prime }+\beta /2){\psi }_0(n^{\prime }+\beta /2+1)}{n+n^{\prime }+\beta }.$$

$$\begin{array}[]{l}\displaystyle\langle H_{Y|X}\rangle=f_{1}H_{00}(\beta)+f_{2}f_{11}H_{11}(\beta)+f_{2}f_{20}H_{20}(\beta)+\dots\end{array}$$

$$⟨I|\mathbf{n}⟩(\beta ,\alpha )$$

$$p(\mathbf{n}|\alpha ,\beta )=p({\mathbf{n}}_x|\alpha )p(\mathbf{n}|\beta ).$$

$$p({\mathbf{n}}_x|\alpha )\propto \frac{\Gamma (1+\alpha )}{\Gamma (N+\alpha )}{\alpha }^{k_1-1},$$

Full text

K****eywords bayesian estimation, mutual information, bias, sampling

1 Introduction

Inferring the statistical dependencies between two variables from a few measured samples is an ubiquitous task in many areas of study. Variables are often linked through non-linear relations, which contain stochastic components. The standard measure employed to quantify the amount of dependency is the mutual information, defined as the reduction in entropy of one of the variables when conditioning the other variable [1, 2]. If the states of the joint distribution are well-sampled, the joint probabilities can be estimated by the observed frequencies, yielding the maximum-likelihood estimator of mutual information. However, this procedure on average over-estimates the mutual information [3, 4, 5], so that independent variables may appear to be correlated, especially when the number of samples is small.

The search for an estimator of mutual information that remains approximately unbiased even with small data samples is an open field of research [6, 7, 8, 9, 10, 11]. Here we focus on discrete variables, and assume it is not possible to overcome the scarceness of samples by grouping elements that are close according to some metric. In addition to corrections that only work in the limit of large samples [12], the state of the art for this problem corresponds to quasi-Bayesian methods that estimate mutual information indirectly through measures of the entropies of the involved variables [8, 13, 14]. These approaches have the drawback of not being strictly Bayesian, since the linear combination of two or more Bayesian estimates of entropies does not, in general, yield a Bayesian estimator of the combination of entropies [8]. The concern is not so much to remain within theoretical Bayesian purity, but rather, to avoid frameworks that may be unnecessarily biased, or where negative estimates of information may arise.

Here we propose a new method for estimating mutual information that is valid in the specific case in which there is an asymmetry between the two variables: One of them has a large number of effective states, and the other only a few. No hypotheses are made about the probability distribution of the large-entropy variable, but the marginal distribution of the low-entropy variable is assumed to be well sampled. The prior is chosen so as to accurately represent the amount of dispersion of the conditional distribution of the low-entropy variable around its marginal distribution. The main finding is that our estimator has very low bias, even in the severely under-sampled regime where there are few coincidences, that is, when a given state of the large-entropy variable is only seldom sampled more than once. The key data statistics that determine the estimated information is the inhomogeneity of the distribution of the low-entropy variable in those states of the high-entropy variable where two or more samples are observed. In addition to providing a practical algorithm to estimate mutual information, our approach sheds light on the way in which just a few samples reveal those specific properties of the underlying joint probability distribution that determine the amount of mutual information.

2 Bayesian approaches to the estimation of entropies

We seek a low-bias estimate of the mutual information between two discrete variables. Let $X$ be a random variable with a large number $k_{x}$ of effective states $\{x_{1},\dots,x_{k_{x}}\}$ with probabilities $q_{x}$, and $Y$ be a variable that varies in a small set $y\in\{y_{1},\dots,y_{k_{y}}\}$, with $k_{y}\ll k_{x}$ . Given the conditional probabilities $q_{y|x}$, the marginal and joint probabilities are $q_{y}=\sum_{x}q_{x}\ q_{y|x}$ and $q_{xy}=q_{x}\ q_{y|x}$, respectively. The entropy $H(Y)$ is

[TABLE]

and can be interpreted as the average number of well-chosen yes/no questions required to guess the sampled value of $Y$ (when using a logarithm of base two). The conditional entropy $H(Y|X)$ is the average uncertainty of the variable $Y$ once $X$ is known,

[TABLE]

The mutual information is the reduction in uncertainty of one variable once we know the other [2]

[TABLE]

Our aim is to estimate $I(X,Y)$ when $Y$ is well sampled, but $X$ is severely undersampled, in particular, when the sampled data contain few coincidences in $X$. Hence, for most values $x$, the number of samples $n_{x}$ is too small to estimate the conditional probability $q_{y|x}$ from the frequencies $n_{xy}/n_{x}$. In fact, when $n_{x}\sim\mathcal{O}(1)$, the maximum likelihood estimator typically underestimates $H(Y|x)$ severely [5], and consequently leads to an overestimation of $I(X,Y)$.

One possibility is to estimate $H(X),H(Y)$ and $H(X,Y)$ using a Bayesian estimator, and then plug the obtained values in Eq. 3 to estimate the mutual information. We now discuss previous approaches to Bayesian estimators for entropy, to later analyze the case of information. For definiteness, we focus on $H(X)$, but the same logic applies to $H(Y)$, or $H(X,Y)$.

The Bayesian estimator is the expected value of $H(\mathbf{q}_{x}|\mathbf{n})$, where $\mathbf{q}_{x}$ are unknown probabilities $q_{x_{1}},\dots,q_{x_{k}}$, and $\mathbf{n}$ represents the number of sampled data $\mathbf{n}=(n_{1},\dots,n_{k})$ obtained in each state. That is,

[TABLE]

Since $p(\mathbf{n}|\mathbf{q})$ is the multinomial distribution

[TABLE]

and since the normalization constant $p(\mathbf{n})$ can be calculated from the integral

[TABLE]

the entire gist of the Bayesian approach is to find an adequate prior $p(\mathbf{q})$ to plug into Eqs. 4, 5 and 6. For the sake of analytical tractability, $p(\mathbf{q})$ is often decomposed into a weighted combination of distributions $p(\mathbf{q}|\beta)$ that can be easily integrated, each tagged by one or or a few parameters, here generically called $\beta$, that vary within a certain domain,

[TABLE]

The decomposition requires to introduce a prior $p(\beta)$. Hence, the former search for an adequate prior $p(\mathbf{q})$ is now replaced by the search for an adequate prior $p(\beta)$. The replacement implies an assumption and also a simplification. The family of priors that can be generated by Eq. 7 does not encompass the entire space of possible priors. The decomposition relies on the assumption that the remaining family is still rich enough to make good inference about the quantity of interest, in this case, the entropy. The simplification stems from the fact that the search for $p(\beta)$ is more restricted than the search for $p(\mathbf{q})$, because the space of possible alternatives is smaller (the dimensionality of $\mathbf{q}$ is typically high, whereas the one of $\beta$ is low). Two popular proposals of Bayesian estimators for entropies are NSB [13] and PYM [14]. In NSB, the functions $p(\mathbf{q}|\beta)$ are Dirichlet distributions, in which $\beta$ takes the role of a concentration parameter. In PYM, these functions are Pitman-Yor processes, and $\beta$ stands for two parameters: one accounting for the concentration, and the other for the so-called discount. In both cases, the Bayesian machinery implies

[TABLE]

where $W(\beta|\mathbf{n})$ is the weight of each $\beta$ in the estimation of the expected entropy

[TABLE]

When choosing the family of functions $p(\mathbf{q}|\beta)$, it is convenient to select them in such a way that the weight $W(\beta|\mathbf{n})$ can be solved analytically. However, this is not the only requirement. In order to calculate the integral in $\beta$, the prior $p(\beta)$ also plays a role. The decomposition of Eq. 7 becomes most useful when the arbitrariness in the choice of $p(\beta)$ is less serious than the arbitrariness in the choice of $p(\mathbf{q})$. This assumption is justified when $W(\beta|\mathbf{n})$ is peaked around a specific $\beta$ value, so that in practice, the shape of $p(\beta)$ hardly has an effect. In these cases, a narrow range of relevant $\beta$ values is selected by the sampled data, and all assumptions about the prior probability outside this range play a minor role. For the choices of the families $p(\mathbf{q}|\beta)$ proposed by NSB and PYM, $W(\beta|\mathbf{n})$ can be calculated analytically, and one can verify that indeed, a few coincidences in the data suffice for a peak to develop. In both cases, the selected $\beta$ is one for which $p(\mathbf{q}|\beta)$ favours a range of $\mathbf{q}$ values that are compatible with the measured data (as assessed by $p(\mathbf{n}|\mathbf{q})$), and also produce non-negligible entropies (Eq. 9).

When the chosen Bayesian estimates of the entropies are plugged into Eq. 3 to obtain an estimate of the information, each term is dominated by its own preferred $\beta$. Since the different entropies are estimated independently, the $\beta$ values selected by the data to dominate the priors $p(q_{x})$ and $p(q_{y})$ need not be compatible with the ones dominating the priors of the joint or the conditional distributions. As a consequence, the estimation of the mutual information is no longer Bayesian, and can suffer from theoretical issues, as for example, yield a negative estimate [8].

A first alternative would be to consider an integrable prior containing a single $\beta$ for the joint probability distribution $q_{xy}$, and then replace $H$ by $I$ in the equations above, to calculate $\langle I\rangle$. This procedure was tested by Archer et al. [8], and the results were only good when the collection of $q_{xy}$ values governing the data were well described by a distribution that was contained in the family of proposed priors $p(\mathbf{q}|\beta)$. The authors concluded that mixtures of Dirichlet priors do not provide a flexible enough family of priors for highly-structured joint distributions, at least for the purpose of estimating mutual information.

To make progress, we note that $I(X,Y)$ can be written as

[TABLE]

where $\mathbf{q}_{y|x}$ and $\mathbf{q}_{y}$ stand for the $k_{y}$-dimensional vectors $(q_{y_{1}|x},\dots,q_{y_{k_{y}}|x})$ and $(q_{y_{1}},\dots,q_{y_{k_{y}}})$, and $D_{\text{KL}}$ represents the Kullback-Leibler divergence. The average divergence between $\mathbf{q}_{y|x}$ and $\mathbf{q}_{y}$ captures a notion of spread. Therefore, the mutual information is sensitive not so much to the value of the probabilities $q_{y|x}$, but rather, to their degree of scatter around the marginal $q_{y}$. The parameters controlling the prior should hence be selected in order to match the width of the distribution of $q_{y|x}$ values, and not so much each probability. With this intuition in mind, in this paper we put forward a new prior for the whole ensemble of conditional probabilities $q_{y|x}$ obtained for different $x$ values. In this prior, the parameter $\beta$ controls the spread of the conditionals $q_{y|x}$ around the marginal $q_{y}$.

3 Prior distribution for conditional entropies

Our approach is valid when the total number of samples $N$ is at least of the order of magnitude of $\sqrt{e^{H(X)}}$, since in this regime, some of the $x$ states are expected to be sampled more than once [15, 16]. In addition, the marginal distribution $q_{y}$ must be well sampled. This regime is typically achieved when $X$ has a much larger set of available states than $Y$. In this case, the maximum likelihood estimators $\hat{q}_{y}$ of the marginal probabilities $q_{y}$ can be assumed to be accurate, that is,

[TABLE]

In this paper, we put forward a Dirichlet prior distribution centered at $\hat{q}_{y}$, that is,

[TABLE]

where $\{\mathbf{q}_{y|x}\}$ contains the $k_{x}$ conditional probabilities $\mathbf{q}_{y|x}$ corresponding to different $x$ values. Large $\beta$ values select conditional probabilities close to $\hat{q}_{y}$, while small values imply a large spread, that pushes the selection towards the border of the $k_{y}$-simplex.

For the moment, for simplicity we work with a prior $p(\{\mathbf{q}_{y|x}\})$ defined on the conditional probabilities $q_{y|x}$, and make no effort to model the prior probability of the vector $\mathbf{q}_{x}$. In practice, we estimate the values of $q_{x}$ with the maximum likelihood estimator $\hat{q}_{x}=n_{x}/N$. Since $X$ is assumed to be severely undersampled, this is a poor procedure to estimate $q_{x}$. Still, the effect on the mutual information turns out to be negligible, since the only role of $q_{x}$ in Eq. 10 is to weigh each of the Kullback-Leibler divergences appearing in the average. If $k_{x}$ is large, each $D_{KL}$ value will appear in several terms of the sum, rendering the individual value of the accompanying $q_{x}$ irrelevant, only the sum of them matters. In Sect. 6, we tackle the full problem of making Bayesian inference both in $\mathbf{q}_{x}$ and $\{\mathbf{q}_{y|x}\}$.

The choice of prior of Eq. 12 is inspired in three facts. First, $\beta$ captures the spread of $q_{y|x}$ around $q_{y}$, as implied by the Kullback-Leibler divergence in Eq. 12. Admittedly, this divergence is not exactly the one governing the mutual information (Eq. 10), since $\mathbf{q}_{y|x}$ and $\mathbf{q}_{y}$ are interchanged. Yet, it is still a measure of spread. The exchange, as well as the denominator in Eq. 12, were introduced for the sake of the second fact, namely, analytical tractability. The third fact regards the emergence of a single relevant $\beta$ when the sampled data begin to register coincidences. If we follow the Bayesian rationale of the previous section, now replacing the entropy by the mutual information, we can again define a weight $W(\beta|\mathbf{n})$ for the parameter $\beta$

[TABLE]

where $F(\beta,\mathbf{n})$ can be obtained analytically, and is a well behaved function of its arguments, whereas

[TABLE]

For each $x$, the vector $\mathbf{q}_{y|x}$ varies in a $k_{y}$-dimensional simplex. For $p(\mathbf{n}|\hat{\mathbf{q}}_{x},\{\mathbf{q}_{y|x}\})$ we take the multinomial

[TABLE]

The important point here, is that the ratio of the Gamma functions of Eq. 13 develops a peak in $\beta$ as soon as the collected data register a few coincidences in $x$. Hence, with few samples, the prior proposed in Eq. 12 renders the choice of $p(\beta)$ inconsequential.

Assuming that the marginal probability of $Y$ is well-sampled, the entropy $H(Y)$ is well approximated by the maximum-likelihood estimator $\hat{H}(Y)=-\sum_{y}(n_{y}/N)\ \log(n_{y}/N)$. For each $\beta$, the expected posterior information can be calculated analytically,

[TABLE]

where $\psi_{0}$ is the digamma function. When the system is well sampled, $n_{xy}\gg 1$, so the effect of $\beta$ becomes negligible, the Digamma functions tend to logarithms, and the frequencies match the probabilities. In this limit, Eq. 15 coincides with the maximal likelihood estimator, which is consistent. The rest of the paper focuses on the case in which the marginal probability of $X$ is severely undersampled.

4 A closer look on the case of a symmetric and binary $Y$-variable

In this section, for simplicity we take $q_{y=0}=q_{y=1}=\nicefrac{{1}}{{2}}$, such that $H(Y)=\log 2$ nats. In this case, the Dirichlet prior of Eq. 12 becomes a Beta distribution

[TABLE]

Large values of $\beta$ mostly select conditional probabilities $q_{y|x}$ close to $1/2$. If all conditional probabilities are similar, and similar to the marginal, the mutual information is low, since the probability of sampling a specific $y$ value hardly depends on $x$. Instead, small values of $\beta$ produce conditional probabilities $q_{y|x}$ around the borders ($q_{y|.}\sim 0$ or $q_{y|.}\sim 1$). In this case, $q_{y|x}$ is strongly dependent on $x$ (see Fig. 1 b), so the mutual information is large. The expected prior mutual information $\langle I(\beta)\rangle$ can be calculated using the analytical approach developed by [17, 14],

[TABLE]

The prior information is a slowly-varying function of the order of magnitude of $\beta$, namely of $\log\beta$. Therefore, if a uniform prior in information is desired, it suffices to choose a prior on $\log\beta$ such that $p(\log\beta)\propto|\partial_{\log\beta}\langle I\rangle(\beta)|$,

[TABLE]

When $k_{y}=2$, the expected posterior information (Eq. 15) becomes

[TABLE]

The marginal likelihood of the data given $\beta$ is also analytically tractable. The likelihood is binomial for each $x$, so

[TABLE]

The posterior for $\beta$ can be obtained by adding a prior $p(\beta)$, as $p(\beta|\mathbf{n})\propto p(\mathbf{n}|\beta)p(\beta)$. The role of the prior becomes relevant when the number of coincidences is too low for the posterior to develop a peak (see below).

In order to gain intuition about the statistical dependence between variables with few samples, we here highlight the specific aspects of the data that influence the estimator of Eq. 19. Grouping together the terms of Eq. (20) that are equal, the marginal likelihood can be rewritten in terms of the multiplicities $m_{nn^{\prime}}$, that is, the number of states $x$ with specific occurrences $\{n_{x1}=n,\,n_{x0}=n^{\prime}\}$ or $\{n_{x1}=n^{\prime},\,n_{x0}=n\}$,

[TABLE]

where

[TABLE]

The posterior for $\beta$ is independent from states $x$ with just a single count, as $p_{10}(\beta)=\text{constant}$. Only states $x$ with coincidences matter. In order to see how the sampled data favor a particular $\beta$, we search for the $\beta$ value that maximizes $\log p(\mathbf{n}|\beta)$ in the particular case where at most two samples coincide on the same $x$, obtaining

[TABLE]

Denoting the fraction of $2$-count states that have one count for each $y$ value as $f_{11}=m_{11}/(m_{11}+m_{20})$, Eq. 23 implies that $\beta\rightarrow\infty$ if $f_{11}\geq 1/2$, and $\beta=f_{11}/(1/2-f_{11})$, otherwise. If the $y$-values are independent of $x$, we expect $f_{11}\sim 1/2$. This case corresponds to a large $\beta$ and, consequently, to a low information. On the other side, for small $f_{11}$, the parameter $\beta$ is also small and the information grows.

In Eq. 23, the data only intervene through $m_{11}$ and $m_{20}$, which characterize the degree of asymmetry of the $y$ values throughout the different $x$ states. This asymmetry, hence, constitutes a sufficient statistics for $\beta$. If a prior $p(\beta)$ is included, the $\beta$ that maximizes the posterior $p(\beta|\mathbf{n})$ may shift, but the effect becomes negligible as the number of coincidences grows.

We now discuss the role of the selected $\beta$ in the estimation of information, Eq. (19), focusing on the conditional entropy $\langle H_{Y|X}\rangle(\beta)$. First, in terms of the multiplicities, the conditional entropy can be rewritten as

[TABLE]

where $f_{k}$ is the fraction of the $N$ samples that fall in states $x$ with $k$ counts, and $f_{nn^{\prime}}$ is the fraction of all states $x$ with $n+n^{\prime}$ counts that have $n$ for one $y$-value (whichever) and $n^{\prime}$ for the other. Finally, $H_{nn^{\prime}}(\beta)$ is the estimation of the entropy of a binary variable after $\{n,n^{\prime}\}$ samples,

[TABLE]

A priori, $\langle I\rangle(\beta)=\log 2-H_{00}(\beta)$, as in Fig. 1d. Surprisingly, from the property $\psi_{0}(z+1)=\psi_{0}(z)+1/z$, it turns out that $H_{00}=H_{10}$ (in fact, $H_{nn}=H_{(n+1)n}$). Hence, if only a single count breaks the symmetry between the two $y$ values, there is no effect on the conditional entropy. This is a reasonable result, since a single extra count is no evidence of an imbalance between the underlying conditional probabilities, it is just the natural consequence of comparing the counts falling on an even number of states (2) when taking an odd number of samples. Expanding the first terms for the conditional entropy,

[TABLE]

In the severely under-sampled regime, these first terms are the most important ones. Typically, $f_{1}$ takes most of the weight, and Eq. (26) implies that the estimation is close to the prior $H_{00}$ evaluated in the value of $\beta$ that maximizes the marginal likelihood (or the posterior).

Finally, we mention that when dealing with few samples, it is important to have not just a good estimate of the mutual information, but also, a confidence interval. Even a small information may be relevant, if the evidence attests that it is strictly above zero. The theory developed here also allows us to estimate the posterior variance of the mutual information, as shown in the Appendix. The variance (Eq. 33) is shown to be inversely proportional to the number of states $k_{x}$, thereby implying that our method benefits from a large number of available states $X$, even if undersampled.

5 Testing the estimator

We now analyze the performance of our estimator in three examples where the number of samples $N$ is below or in the order of the effective size of the system $\exp(H_{XY})$. In this regime, most observed $x$ states have very few samples. In each example, we define the probabilities $q_{x}$ and $q_{1|x}$ with three different criteria, giving rise to collections of probabilities that can be described with varying success by the prior proposed in this paper, Eq. 16. Once the probabilities are defined, the true value $I_{XY}$ of the mutual information can be calculated, and compared to the one estimated by our method, as well as by three other estimators employed in the literature, in 50 different sets of samples $\mathbf{n}$ of the measured data. As our estimator we use $\langle I|\mathbf{n}\rangle$ from Eq. (19) evaluated in the $\beta$ that maximizes the marginal likelihood $p(\mathbf{n}|\beta)$. We did not observe any improvement when integrating over the whole posterior $p(\beta|\mathbf{n})$ with the prior $p(\beta)$ of Eq. 18, except when $m_{20}$ or $m_{10}$ were of order 1. This fact implies the existence of a well-defined peak in the marginal likelihood.

In the first example (Fig. 2a, d), the probabilities $q_{x}$ are obtained by sampling a Pitman-Yor distribution with concentration parameter $\alpha=50$ and tail parameter $d=0.55$. These values correspond to a PYM prior with a heavy tail. The conditional probabilities $q_{y|x}$ are defined by sampling a symmetric Beta distribution $q_{y|x}\sim\text{Beta}(\beta/2,\beta/2)$, as in Eq. 16. In Fig. 2a, we use $\beta=2.3$. Once the joint probability $q_{xy}$ is defined, 50 sets of samples $\mathbf{n}$ are generated. The effective size of the system is $\exp(H_{XY})\simeq 800$. We compare our estimator to maximum likelihood (ML), NSB and PYM when applied to $H_{X}$ and $H_{XY}$ (all methods coincide in the estimation of $H_{Y}$). Our estimator has a low bias, even when the number of samples per effective state is as low as $N/e^{H_{XY}}=0.15$. The variance is larger than ML, comparable to NSB and smaller than PYM. All the other methods (ML, NSB and to a lesser extend PYM) overestimate the mutual information. In Fig. 2d, the performance of the estimators is also tested for different values of the exact mutual information $I_{XY}$, which we explore by varying $\beta\in(0.04,\,14)$. For each $\beta$, the conditional probabilities $q_{1|x}$ are sampled once. Each vector $\mathbf{n}$ contains $N=500$ samples, and $\mathbf{n}$ is sampled $50$ times. Our estimates have very low bias, even as the mutual information goes to zero —namely, for independent variables.

Secondly, we analyze an example where the statistical relation between $X$ and $Y$ is remarkably intricate (example inspired by [19]), which underscores the fact that making inference about the mutual information does not require inferences on the joint probability distribution. The variable $x$ is a binary vector of dimension $12$. Each component represents the presence or absence of one of a maximum of $12$ delta functions equally spaced on the surface of a sphere. There are $2^{12}$ possible $x$ vectors, and they are governed by a uniform prior probability: $q_{x}=2^{-12}$. The conditional probabilities are generated in such a way that they be invariant under rotations of the sphere, that is, $q_{y|x}=q_{y|\text{R}(x)}$, where $R$ is a rotation. Using a spherical harmonic representation [18], the frequency components $\pi_{\ell}(f(x))$ of the spherical spectrum are obtained, where $f(x)$ is the combination of deltas. The conditional probabilities $q_{y|x}$ are defined as a sigmoid function of $(\pi_{0}-\pi_{1}-\pi_{2})$. The offset of the sigmoid is chosen such that $q_{y=1}\simeq 0.5$, and the gain such that $I_{XY}\simeq 0.5$ nats. In this example, and unlike the Dirichlet prior implied by our estimator, $p(q_{y|x})$ has some level of roughness (inset in Fig. 2b), due to peaks coming from the invariant classes in $\{x_{1},\dots,x_{2^{12}}\}$. Hence, the example does not truly fit into the hypothesis of our method. With these settings, the effective size of the system is $\exp(H_{XY})\simeq 5000$. Our estimator has little bias (Figs. 2b, e), even with $N/e^{H_{XY}}=0.2$ samples per effective state. In this regime, around $\sim 80\%$ of the samples fall on $x$ states that occur only once ($f_{1}\simeq 0.8$), $\sim 19\%$ on states that occur twice and $\sim 1\%$ on states with $3$ counts, or maybe $4$. As mentioned above, in such cases, the value of $I_{XY}$ is very similar to the one that would be obtained by evaluating the prior information $\langle I|\mathbf{n}=\mathbf{0},\beta\rangle$ of (Eq. 17) at the $\beta$ that maximizes the marginal likelihood $p(\mathbf{n}|\beta)$, which in turn is mainly determined by $f_{11}$. In Fig. 2e, the estimator is tested with a fixed number of samples $N=2000$ for different values of the mutual information, which we explore by varying the gain of the sigmoid. The bias of the estimate is small in the entire range of mutual informations.

In the third place, we consider an example where the conditional probabilities are generated from a distribution that is poorly approximated by a Dirichlet prior. The conditional probabilities are sampled from three Dirac deltas, as $q_{y|x}\sim[0.5\,\delta(q-\nicefrac{{1}}{{2}})+0.25\,\delta(q-q_{0})+0.25\,\delta(q-1+q_{0})]$, with $q_{0}=0.1$. The delta placed in $q=\nicefrac{{1}}{{2}}$ could be approximated by a Dirichlet prior with a large $\beta$, while the other two deltas could be approximated by a small $\beta$, but there is no single value of $\beta$ that can approximate all three deltas at the same time. The $x$ states are generated as Bernoulli ($p=0.05$) binary vectors of dimension $D=40$, while the conditional probabilities $q_{1|x}$ depend on the parity of the sum of the components of the vector $x$. When the sum is even, we assign $q_{y|x}=\nicefrac{{1}}{{2}}$, and when it is odd, we assign $q_{y|x}=q_{0}$ or $q_{y|x}=1-q_{0}$, both options with equal probability. Although in this case our method has some degree of bias, it still preserves a good performance in relation to the other approaches (see Fig. 2c, f). The marginal likelihood $p(\mathbf{n}|\beta)$ contains a single peak in an intermediate value of $\beta$, coinciding with none of the deltas in $p(q_{1|x})$, but still capturing the right value of the mutual information. As in the previous examples we also test the performance of the estimator for different values of the mutual information, varying in this case the value of $q_{0}$ (with $N=2000$). Our method performs acceptably for all values of mutual information. The other methods, instead, are challenged more severely, probably because a large fraction of the $x$ states have a very low probability, and are therefore difficult to sample. Those states, however, provide a crucial contribution to the relative weight of each of the three values of $q_{1|x}$. PYM, in particular, sometimes produces a negative estimate for $I_{XY}$.

Finally, we check numerically the accuracy of the analytically predicted mean posterior information (Eq. 19) and variance (Eq. 33) in the severely under-sampled regime. The test is performed in a different spirit than the numerical evaluations of Fig. 2. There, averages were taken for multiple samples of the vector $\mathbf{n}$, from a fixed choice of the probabilities $\mathbf{q}_{x}$ and $\mathbf{q}_{y|x}$. The averages of Eqs. 19 and Eq. 33, however, must be interpreted in the Bayesian sense. The square brackets in $\langle I|\mathbf{n}\rangle$ and $\langle H^{2}_{Y|X}\rangle$ represent averages taken for a fixed data sample $\mathbf{n}$, and unknown underlying probability distributions $\mathbf{q}_{x}$ and $\mathbf{q}_{y|x}$. We generate many such distributions with $q_{x}\sim\text{DP}(\alpha)$ (a Dirichlet Process with concentration parameter $\alpha$) and $q_{y|x}\sim\text{Beta}(\beta/2,\beta/2)$. A total of $13,500$ distributions $q_{xy}$ are produced, with $\log\beta$ sampled from Eq. 18, and three equiprobable values of $\alpha=\{e^{4},e^{5},e^{6}\}$. For each of these distributions we generate five (5) sets of just $N=40$ samples, thereby constructing a list of $5\times 13,500$ cases, each case characterized by specific values of $\alpha,\beta,\mathbf{q}_{x},\{\mathbf{q}_{y|x}\},I(\mathbf{q}_{x},\{\mathbf{q}_{y|x}\}),\mathbf{n},\langle I|\mathbf{n}\rangle$ and $\sigma^{2}(I|\mathbf{n})$ . Following the Bayesian rationale, we partition this list in classes, each class containing all the cases that end up in the same set of multiplicities $\{m_{nn^{\prime}}\}$ —for example, $\{m_{10}=36,\,m_{20}=2\}$. For each of the $100$ most occurring sets of multiplicities (which together cover $70\%$ of all the cases), we calculate the mean and the standard deviation of the mutual information $I(\mathbf{q},\{\mathbf{q}_{y|x}\})$ of the corresponding class, and compare them with our predicted estimates $\langle I|\{m_{nn^{\prime}}\}\rangle$ and $\langle\sigma^{2}_{I}|\{m_{nn^{\prime}}\}\rangle$, using the prior $p(\log\beta)$ from Eq. (18). Figure 3 shows a good match between the numerical ($y$-axis) and analytical ($x$-axis) averages that define the mean information (panel a) and the standard deviation (b). The small departures from the diagonal stem from the fact that the analytical average contains all the possible $\mathbf{q}_{x}$ and $\{\mathbf{q}_{y|x}\}$, even if some of them are highly improbable for one given set of multiplicities. The numerical average, instead, includes the subset of the 13,500 explored cases that produced the tested multiplicity. All the depicted subsets contained many cases, but still, they remained unavoidably below the infinity covered by the theoretical result.

We have also tested cases where $Y$ takes more than two values, and where the marginal distribution $q_{y}$ is not uniform, observing similar performance of our estimator.

6 A prior distribution for the large entropy variable

The prior considered so far did not model the probability $\mathbf{q}_{x}$ of the large-entropy variable $X$. Throughout the calculation, the probabilities $\mathbf{q}_{x}$ were approximated by the maximum likelihood estimator $\hat{q}_{X}=n_{x}/N$. Here we justify such procedure by demonstrating that proper Bayesian inference on $\mathbf{q}_{x}$ hardly modifies the estimation of the mutual information. To that end, we replace the prior of Eq. 12 by another prior that depends on both $\mathbf{q}_{x}$ and $\{\mathbf{q}_{y|x}\}$.

The simplest hypothesis is to assume that the prior $p(\mathbf{q}_{x},\{\mathbf{q}_{y|x}\})$ factorizes as $p(\mathbf{q}_{x})\ p(\{\mathbf{q}_{y|x}\})$, implying that the marginal probabilities $q_{x}$ are independent of the conditional probabilities $q_{y|x}$. We propose $q_{x}\sim\text{DP}(\alpha)$, so that the marginal probabilities $q_{x}$ are drawn from a Dirichlet Process with concentration parameter $\alpha$, associated to the total number of pseudo-counts. After integrating in $\mathbf{q}_{x}$ and in $\mathbf{q}_{y|x}$, the mean posterior mutual information for fixed hyper-parameters $\beta$ and $\alpha$ is

[TABLE]

Before including the prior $p(\mathbf{q}_{x})$, in the severely undersampled regime the mean posterior information was approximately equal to the prior information evaluated in the best $\beta$ (Eq. 15). The new calculation (Eq. 27) contains the prior information explicitly, weighted by $\alpha/(N+\alpha)$, that is, the ratio between the number of pseudo-counts from the prior and the total number of counts. Thereby, the role of the non-observed (but still inferred) states is established.

The independence assumed between $\mathbf{q}_{x}$ and $\{\mathbf{q}_{y|x}\}$ implies that

[TABLE]

The inference over $\alpha$ coincides with the one of PYM with the tail parameter as $d=0$ [14], since

[TABLE]

where $k_{1}=\sum_{x,n_{x}>0}1$ is the number of states $x$ with at least one sample. With few coincidences in $x$, $p(\mathbf{n}_{x}|\alpha)$ develops a peak around a single $\alpha$ value that represents the number of effective states. Compared to the present Bayesian approach, maximum likelihood underestimates the number of effective states (or entropy) in $x$. Since the expected variance of the mutual information decreases with the square root of the number of effective states, the Bayesian variance is reduced with respect to the one of ML.

7 Discussion

In this work we propose a novel estimator for mutual information of discrete variables $X$ and $Y$, which is adequate when $X$ has a much larger number of effective states than $Y$. If this condition does not hold, the performance of the estimator breaks down. We inspire our proposal in the Bayesian framework, in which the core issue can be boiled down to finding an adequate prior. The more the prior is dictated by the data, the less we need to assume from outside. Equation 10 implies that the mutual information $I(X,Y)$ is the spread of the conditional probabilities of one of the variables (for example, $q_{y|x}$, but the same holds for $q_{x|y}$) around the corresponding marginal ($q_{y}$ or $q_{x}$, respectively). This observation inspires the choice of our prior (Eq. 12), which is designed to capture the same idea, and in addition, to be analytically tractable. We choose to work with an hyper-parameter $\beta$ that regulates the scatter of $q_{y|x}$ around $q_{y}$, and not the scatter of $p_{x|y}$ around $q_{x}$, because the asymmetry in the number of available states of the two variables makes the $\beta$ of the first option (and not the second) strongly modulated by the data, by the emergence of a peak in $p(\mathbf{n}|\beta)$.

Although our proposal is inspired in previous Bayesian studies, the procedure described here is not strictly Bayesian, since our prior (Eq. 12) requires the knowledge of $\hat{q}_{y}$, which depends on the sampled data. However, in the limit in which $q_{y}$ is well sampled, this is a pardonable crime, since $\hat{q}_{y}$ is defined by a negligible fraction of the measured data. Still, Bayesian purists should employ a two-step procedure to define their priors. First, they should perform Bayesian inference on the center of the Dirichlet distribution of Eq. 12 by maximizing $p(q_{y}|\mathbf{n})$, and then replace $\hat{q}_{y}$ in Eq. 12 by the inferred $q_{y}$. For all practical purposes, however, if the conditions of validity of our method hold, both procedures lead to the same result.

By confining the set or possible priors $p(\{\mathbf{q}_{y|x}\})$ to those generated by Eq. 12 we relinquish all aspiration to model the prior of, say, $q_{y|x=3}$, in terms of the observed frequencies at $x=3$. In fact, the preferred $\beta$ value is totally blind to the specific $x$ value of each sampled datum. Only the number of $x$-values containing different counts of each $y$-value matters. Hence, the estimation of mutual information is performed without attempting to infer the specific way the variables $X$ and $Y$ are related, a property named equitability [20], and that is shared also by other methods [13, 8, 14]. Although this fact may be seen as a disadvantage, deriving a functional relation between the variables can actually bias the inference on mutual information [20]. Moreover, fitting a relation is unreasonable in the severe under-sampled regime, in which not all $x$-states are observed, most sampled $x$-states contain a single count, and few $x$-states contain more than two counts. At least, without a strong assumption about the probability space. In fact, if the space of probabilities of the involved variables has some known structure or smoothness condition, other approaches that estimate information by fitting the relation first may perform well [9, 10, 11]. Part of the approach developed here could be extended to continuous variables or spaces with a determined metric. This extension is left for future work.

The main result of the paper, is that our estimator has small bias, even in the severely under-sampled regime. It outperforms other estimators discussed in the literature (at least, when the conditions of validity hold), and by construction, it never produces negative values. More importantly, it even works in cases where the collection of true conditional probabilities $q_{y|x}$ is not contained in the family of priors generated by $p(\mathbf{q}|\beta)$, as demonstrated by the second and third examples of Sect. 5. In these cases, the success of the method relies on the peaked nature of the posterior distribution for $\beta$. Even if the selected $p(\mathbf{q}|\beta)$ provides a poor description of the actual collection of probabilities, the dominant $\beta$ captures the right value of mutual information. This is the sheer instantiation of the equitability property discussed above.

Our method provides also a transparent way to identify the statistics that matter, out of all the measured data. Quite naturally, the $x$ states that have not been sampled provide no evidence in shaping $p(\beta|\mathbf{n})$, as indicated by Eq. 13, and only shift the posterior information towards the prior (Eq. 27). More interestingly, the $x$ states with just a single count are also irrelevant, both in shaping $p(\beta|\mathbf{n})$ and in modifying the posterior information away from the prior. These states are unable to provide evidence about the existence of either flat or skewed conditional probabilities $q_{y|x}$. Only the states $x$ that have been sampled at least twice contribute to the formation of a peak in $p(\beta|\mathbf{n})$, and in deviating the posterior information away from the prior.

Several fields can benefit from the application of our estimator of mutual information. Examples can be found in neuroscience, when studying whether neural activity (a variable with many possible states) correlates with a few selected stimuli or behavioral responses [12, 21, 22], or in genomics, to understand associations between genes (large-entropy variable) and a few specific phenotypes [23]. The method can also shed light into the development of rate-distortion methods to be employed in situations in which only a few samples are available [24, 25]. The possibility of detecting statistical dependences with only few samples is of key importance, not just for analyzing data sets, but also to understand how living organisms quickly infer dependencies in their environments and adapt accordingly [26].

Funding

This research was funded by CONICET, CNEA, ANPCyT Raíces 2016 grant number 1004.

Acknowledgments

We thank Ilya Nemenman for his fruitful comments and discussions.

Appendix A Expected variance for a symmetric, binary $Y$-variable

The posterior variance of the mutual information is

[TABLE]

In the first place, we demonstrate that this quantity is proportional to $k_{x}^{-1}$, implying that our estimator becomes increasingly accurate as the number of states of the $X$-variable increases. Given that

[TABLE]

with

[TABLE]

it is easy to show that

[TABLE]

In other words, if the marginal entropy is well sampled, the variance in the information is mainly due to the variance in the conditional entropy. In turn, $H_{Y|X}$ is defined as an average of $k_{x}$ terms (Eq. 32). The independence hypothesis implied in the prior of Eq. 12, and in the way different $q_{1|x}$ and $q_{x}$ factor out in $q(\mathbf{n}|\mathbf{q}_{x},\{\mathbf{q}_{y|x}\})$ (Eq. 5), imply that the different terms of (Eq. 32) are all independent of each other. The average of $k_{x}$ independent terms has a variance proportional to 1/$k_{x}$, so the estimator proposed here becomes increasingly accurate as $k_{x}$ grows.

We now derive the detailed dependence of $\langle H_{Y|X}^{2}\rangle-\langle H_{Y|X}\rangle^{2}$ on the sampled data $\mathbf{n}$. The mean conditional entropy $\langle H_{Y|X}\rangle$ can be written in terms of $\langle H_{Y|x}\rangle(n_{x1},n_{x0},\beta)$, that is, of the entropy of the variable $Y$ for a particular state $x$ with $n_{x}=n_{x0}+n_{x1}$ counts at fixed $\beta$

[TABLE]

Similarly, for the second moment,

[TABLE]

In turn, $\text{Var}[H_{Y|x}((n_{x0},n_{x1},\beta)]=\langle H_{Y|x}^{2}(n_{x0},n_{x1},\beta)\rangle-\langle H_{Y|x}(n_{x0},n_{x1},\beta)\rangle^{2}$, where the first moment $\langle H_{Y|x}\rangle(n_{x0},n_{x1},\beta)$ is given in Eq. 34. The second moment is [14],

[TABLE]

In this equation,

[TABLE]

and $\psi_{1}(z)$ is the first polygamma function.

Replacing the obtained expressions in Eq. 33, the variance of the estimated information is obtained. The two terms of Eq. (35) represent the two sources of uncertainty of the conditional entropy: the uncertainty of $\beta$ (first term), manifested in the width of $p(\beta|\mathbf{n})$, and the uncertainty of the conditional entropies for a fixed $\beta$ (second term), manifested in the width of $p(\{q_{1|x}\}|\beta)$. As the number of samples decreases, the uncertainty in $\beta$ becomes the dominant term.

Finally, we need to mention that the approximate symbol in Eq. 31 stems from the fact that we are assuming that $H(Y)$ is well approximated by its maximum-likelihood estimator. We are therefore neglecting the error in the marginal entropy $H_{Y}$, and assuming that the error in the mutual information only stems from the uncertainty in the conditional entropy $H_{Y|X}$ (Eq. 33). This assumption is well justified in the context explored in this paper, that is, when $H(X)\gg H(Y)$.