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

Abstract

Determining the strength of nonlinear, 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 a consistent 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 the two variables, without seeking to model the specific way in which they 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 in which the large-entropy variable registers coincidences.

1. Introduction

Inferring the statistical dependencies between two variables from a few measured samples is a ubiquitous task in many areas of study. Variables are often linked through nonlinear, 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. Replacing such estimates in the formula for the mutual information yields the so-called “plug-in” estimator of mutual information. However, unless all the states of the joint distribution are sampled extensively, this procedure typically over-estimates the mutual information [3,4,5]. In fact, when the number of samples is of the order of the effective number of joint states, even independent variables tend to appear as correlated.

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 states 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 variable has a large number of effective states, and the other only a few. Examples of asymmetric problems can be found for instance in neuroscience, when assessing the statistical dependence between the activity of a population of neurons and a few stereotyped behaviors, as “lever to the left” or “lever to the right.” If the neural activity is represented as a collection of binary strings, with zeroes and ones tagging the absence or presence of spikes in small time bins, the set of possible responses is huge, and typically remains severely undersampled in electrophysiological experiments. If the behavioral paradigm is formulated in terms of a binary forced choice, the behavioral response only has two possible states.

In our approach, 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. This prior is motivated in the framework of Bayesian estimation of information, and then tested with selected examples in which the hypotheses implied in the prior are fulfilled in varying degree, from high to low. The examples show that our estimator has very low bias, even in the severely under-sampled regime, where there are few coincidences that is, when any given state of the large-entropy variable has a low probability of being 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,\cdots ,x_{k_x}\}$ with probabilities $q_x$, and Y be a variable that varies in a small set $y\in \{y_1,\cdots ,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\left(Y\right)$ is

$${\displaystyle H\left(Y\right)=-\sum _y{q}_{y}log{q}_y,}$$

(1)

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). Since the entropy is a function of the probabilities, and not of the actual values taken by the random variable, here we follow the usual notation in which the expressions $H\left(Y\right)$ and $H\left(\mathbf{q}\right)=H(q_{y_1},\cdots ,q_{y_{k_y}})$ are taken to represent the same concept, defined in Equation (1). The conditional entropy $H\left(Y\right|X)$ is the average uncertainty of the variable Y once X is known,

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

(2)

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

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

(3)

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$. The plug-in estimators of entropy, conditional entropy and mutual information are defined as the ones in which all the probabilities appearing in Equations (1)–(3) are estimated naïvely by the frequencies obtained in a given sample. In fact, when $n_x\sim \mathcal{O}\left(1\right)$, the plug-in estimator typically underestimates $H\left(Y\right|x)$ severely [5], and often leads to an overestimation of $I(X,Y)$.

One possibility is to estimate $H\left(X\right),H\left(Y\right)$ and $H(X,Y)$ using a Bayesian estimator, and then insert the obtained values in Equation (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\left(X\right)$, but the same logic applies to $H\left(Y\right)$, or $H(X,Y)$.

We seek the Bayesian estimator of the entropy conditional to having sampled the state $x_i$ a number $n_i$ of times. We denote the collection of sampled data as a vector $\mathbf{n}=(n_1,\cdots ,n_{k_x})$, and the total number of samples $N={\sum }_i{n}_i$. The Bayesian estimator of the entropy is a function that maps each vector $\mathbf{n}$ on a non-negative real number $\widehat{H}\left(\mathbf{n}\right)$, in such a way that the posterior expected error is minimized [15,16]. In the Bayesian framework, the sampled data $\mathbf{n}$ are known, whereas the underlying distribution $\mathbf{q}$ originating the data is unknown, and must be inferred. In this framework, the estimator of the entropy is the function $\widehat{H}\left(\mathbf{n}\right)$ for which

$$〈{\left[\widehat{H}\left(\mathbf{n}\right)-H\right]}^2〉=\int {\left[\widehat{H}\left(\mathbf{n}\right)-H\left(\mathbf{q}\right)\right]}^{2}p\left(\mathbf{q}\right|\mathbf{n})\mathrm{d}\mathbf{q}$$

(4)

is minimized. The Bayesian nature of this framework is embodied in the fact that the integral weighs all candidate distributions $\mathbf{q}$ that could have originated the data with the posterior distribution $p\left(\mathbf{q}\right|\mathbf{n})$, which can be related to the multinomial distribution $p\left(\mathbf{n}\right|\mathbf{q})$ by means of Bayes rule. The estimator $\widehat{H}\left(\mathbf{n}\right)$ can be found by taking the derivative of Equation (4) with respect to $\widehat{H}$ and equating to zero. The result is [17]

$$\begin{array}{ccc}\hfill \widehat{H}\left(\mathbf{n}\right)=\langle H|\mathbf{n}\rangle & =& \int H\left(\mathbf{q}\right)p\left(\mathbf{q}\right|\mathbf{n})\mathrm{d}\mathbf{q}\hfill \\ & =& {\left[p\left(\mathbf{n}\right)\right]}^{-1}\int H\left(\mathbf{q}\right)p\left(\mathbf{n}\right|\mathbf{q})p\left(\mathbf{q}\right)\mathrm{d}\mathbf{q},\hfill \end{array}$$

(5)

where the first equation results from the minimization procedure, and the second, from using Bayes rule. As a result, the Bayesian estimator is the expected value of $H\left(\mathbf{q}\right|\mathbf{n})$.

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

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

(6)

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

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

(7)

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

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

(8)

The decomposition requires to introduce a prior $p\left(\beta \right)$. Hence, the former search for an adequate prior $p\left(\mathbf{q}\right)$ is now replaced by the search for an adequate prior $p\left(\beta \right)$. The replacement implies an assumption and also a simplification. The family of priors that can be generated by Equation (8) does not necessarily 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\left(\beta \right)$ is more restricted than the search for $p\left(\mathbf{q}\right)$ because the space of 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 Nemenman–Shafee–Bialek (NSB) [13] and Pitman–Yor Mixture (PYM) [14]. In NSB, the functions $p\left(\mathbf{q}\right|\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

$$\langle H|\mathbf{n}\rangle =\frac{1}{p\left(\mathbf{n}\right)}\int \mathrm{d}\beta p\left(\beta \right)W\left(\beta \right|\mathbf{n}),$$

(9)

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

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

(10)

When choosing the family of functions $p\left(\mathbf{q}\right|\beta )$, it is convenient to select them in such a way that the weight $W\left(\beta \right|\mathbf{n})$ can be solved analytically. However, this is not the only requirement. In order to calculate the integral in $\beta$, the prior $p\left(\beta \right)$ also plays a role. The decomposition of Equation (8) becomes most useful when the arbitrariness in the choice of $p\left(\beta \right)$ is less serious than the arbitrariness in the choice of $p\left(\mathbf{q}\right)$. This assumption is justified when $W\left(\beta \right|\mathbf{n})$ is peaked around a specific $\beta$-value, so that, in practice, the shape of $p\left(\beta \right)$ 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\left(\mathbf{q}\right|\beta )$ proposed by NSB and PYM, $W\left(\beta \right|\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\left(\mathbf{q}\right|\beta )$ favours a range of $\mathbf{q}$ values that are compatible with the measured data (as assessed by $p\left(\mathbf{n}\right|\mathbf{q})$), and also produce non-negligible entropies (Equation (10)).

When the chosen Bayesian estimates of the entropies are plugged into Equation (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\left(q_x\right)$ and $p\left(q_y\right)$ 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 negative estimates [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\left(\mathbf{q}\right|\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

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

(11)

where ${\mathbf{q}}_{y|x}$ and ${\mathbf{q}}_y$ stand for the $k_y$-dimensional vectors $(q_{y_1|x},\cdots ,q_{y_{k_y}|x})$ and $(q_{y_1},\cdots ,q_{y_{k_y}})$, and $D_{\mathrm{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. A Prior Distribution for the 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\left(X\right)}}$, since in this regime, some of the x states are expected to be sampled more than once [18,19]. 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 ${\widehat{q}}_y$ of the marginal probabilities $q_y$ can be assumed to be accurate that is,

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

(12)

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

$$\begin{array}{ccc}\hfill p\left(\left\{{\mathbf{q}}_{y|x}\right\}\right|\beta )& =& \mathsf{\Gamma }{\left(\beta \right)}^{k_x}\prod _{xy}\frac{{q_{y|x}}^{\beta {\widehat{q}}_y-1}}{\mathsf{\Gamma }\left(\beta {\widehat{q}}_y\right)}\hfill \\ & =& {\left[\frac{\mathsf{\Gamma }\left(\beta \right)}{{\prod }_y\mathsf{\Gamma }\left(\beta {\widehat{q}}_y\right)}\right]}^{k_x}{\mathrm{e}}^{-\beta k_{x}H\left({\widehat{\mathbf{q}}}_y\right)}\frac{exp\left[-\beta {\sum }_{xy}{q}_y(log{q}_y-log{q}_{y|x})\right]}{{\prod }_{xy}{q}_{y|x}}\hfill \\ & \propto & \frac{exp\left[-\beta {\sum }_x{D}_{\mathrm{KL}}(\widehat{{\mathbf{q}}_y}\left|\right|{\mathbf{q}}_{y|x})\right]}{{\prod }_{xy}{q}_{y|x}},\hfill \end{array}$$

(13)

where $\left\{{\mathbf{q}}_{y|x}\right\}$ contains the $k_x$ conditional probabilities ${\mathbf{q}}_{y|x}$ corresponding to different x values. Large $\beta$ values select conditional probabilities close to ${\widehat{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\left(\left\{{\mathbf{q}}_{y|x}\right\}\right)$ defined on the conditional probabilities $q_{y|x}$, and make no effort to model the prior probability of ${\mathbf{q}}_x$. In practice, we estimate the values of $q_x$ with the maximum likelihood estimator ${\widehat{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 Equation (11) is to weigh each of the Kullback–Leibler divergences appearing in the average. If $k_x$ is large, each $D_{KL}$-value appears in several terms of the sum, rendering the individual value of the accompanying $q_x$ irrelevant, only the sum of all those with the same $D_{\mathrm{KL}}$ matters. In Section 6, we tackle the full problem of making Bayesian inference both in ${\mathbf{q}}_x$ and $\left\{{\mathbf{q}}_{y|x}\right\}$.

The choice of prior of Equation (13) 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 Equation (13). Admittedly, this divergence is not exactly the one governing the mutual information (Equation (11)), since ${\mathbf{q}}_{y|x}$ and ${\mathbf{q}}_y$ are swapped. Yet, it is still a measure of spread. The exchange, as well as the denominator in Equation (13), 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\left(\beta \right|\mathbf{n})$ for the parameter $\beta$

$$\begin{array}{ccc}\hfill W\left(\beta \right|\mathbf{n})& =& \int \left\{\mathrm{d}{\mathbf{q}}_{y|x}\right\}I({\widehat{\mathbf{q}}}_x,\left\{{\mathbf{q}}_{y|x}\right\})p\left(\mathbf{n}\right|\widehat{{\mathbf{q}}_x},\left\{{\mathbf{q}}_{y|x}\right\})p\left(\left\{{\mathbf{q}}_{y|x}\right\}\right|\beta )\hfill \\ & =& p\left(\beta \right|\mathbf{n}\left)F\right(\beta ,\mathbf{n}),\hfill \end{array}$$

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

$$\begin{array}{ccc}\hfill p\left(\beta \right|\mathbf{n})& =& \frac{p\left(\beta \right)p\left(\mathbf{n}\right|\beta )}{p\left(\mathbf{n}\right)}=\frac{p\left(\beta \right)}{p\left(\mathbf{n}\right)}\int \left\{\mathrm{d}{\mathbf{q}}_{y|x}\right\}p\left(\mathbf{n}\right|{\widehat{\mathbf{q}}}_x,\left\{{\mathbf{q}}_{y|x}\right\})p\left(\left\{{\mathbf{q}}_{y|x}\right\}\right|\beta )\hfill \\ & =& \frac{p\left(\beta \right)}{p\left(\mathbf{n}\right)}\prod _x\frac{\mathsf{\Gamma }\left(\beta \right)}{\mathsf{\Gamma }(n_x+\beta )}\prod _{y=1}^{k_y}\frac{\mathsf{\Gamma }(n_{xy}+\beta {\widehat{q}}_y)}{\mathsf{\Gamma }\left(\beta {\widehat{q}}_y\right)}.\hfill \end{array}$$

(14)

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

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

(15)

The important point here is that the ratio of the Gamma functions of Equation (14) 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 Equation (13) renders the choice of $p\left(\beta \right)$ inconsequential.

Assuming that the marginal probability of Y is well-sampled, the entropy $H\left(Y\right)$ is well approximated by the plug-in estimator $\widehat{H}\left(Y\right)=-{\sum }_y(n_y/N)log(n_y/N)$. For each $\beta$, the expected posterior information can be calculated analytically,

$$\langle I|\mathbf{n},\beta \rangle =\widehat{H}\left(Y\right)-\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],$$

(16)

where ${\psi }_0$ is the digamma function. A code implementing the estimate of Equation (16) is publicly available in the site mentioned in the Supplementary Material below.

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, Equation (16) coincides with the plug-in estimator, which is consistent [20]. As a consequence, our estimator is also 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 $\mathit{Y}$-Variable

In this section, we analyze the case of a binary Y-variable, which, for simplicity, is assumed to be symmetric that is, $q_{y=0}=q_{y=1}=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$, such that $H\left(Y\right)=log2$ nats. In this case, the Dirichlet prior in each factor of Equation (13) becomes a Beta distribution

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

(17)

and $p\left(\left\{q_{1|x}\right\}|\beta \right)={\prod }_{x}p\left(q_{1|x}\right|\beta )$. Large values of $\beta$ mostly select conditional probabilities $q_{y|x}$ close to $\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$. 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|x}\sim 0$ or $q_{y|x}\sim 1$). In this case, $q_{y|x}$ is strongly dependent on x (see Figure 1b), so the mutual information is large. The expected prior mutual information $\langle I|\mathbf{n}=\mathbf{0},\beta \rangle$ can be calculated using the analytical approach developed by [14,17],

$$\begin{array}{cc}\hfill \langle I|\mathbf{n}=\mathbf{0},\beta \rangle & =log2-{\psi }_0(\beta +1)+{\psi }_0(\beta /2+1).\hfill \end{array}$$

(18)

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|\mathbf{n}=\mathbf{0},\beta \rangle |$,

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

(19)

When $k_y=2$, the expected posterior information (Equation (16)) becomes

$$\langle I|\mathbf{n},\beta \rangle =\widehat{H}\left(Y\right)-\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].$$

(20)

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

$$\begin{array}{ccc}\hfill p\left(\mathbf{n}\right|\beta )& =& \prod _x{\int }_0^1\mathrm{d}{q}_{1|x}p(n_{x1},n_{x0}|q_{1|x})p\left(q_{1|x}\right|\beta )\hfill \\ & \propto & \prod _x\frac{\mathsf{\Gamma }(n_{x1}+\beta /2)\mathsf{\Gamma }(n_{x0}+\beta /2)\mathsf{\Gamma }\left(\beta \right)}{\mathsf{\Gamma }(n_x+\beta ){\mathsf{\Gamma }(\beta /2)}^2}.\hfill \end{array}$$

(21)

The posterior for $\beta$ can be obtained by adding a prior $p\left(\beta \right)$, as $p\left(\beta \right|\mathbf{n})\propto p(\mathbf{n}\left|\beta \right)p\left(\beta \right)$. 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 Equation (20). Grouping together the terms of Equation (21) that are equal, the marginal likelihood can be rewritten in terms of the multiplicities $m_{n{n}^{\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\}$,

$$\begin{array}{cc}\hfill logp\left(\mathbf{n}\right|\beta )& {\displaystyle =\sum _{n\ge n^{\prime }}{m}_{n{n}^{\prime }}log\left[\frac{\mathsf{\Gamma }(n+\beta /2)\mathsf{\Gamma }(n^{\prime }+\beta /2)\mathsf{\Gamma }\left(\beta \right)}{\mathsf{\Gamma }(n+n^{\prime }+\beta )\mathsf{\Gamma }{(\beta /2)}^2}\right]}\hfill \\ & {\displaystyle =\sum _{n\ge n^{\prime }}{m}_{n{n}^{\prime }}log{p}_{n{n}^{\prime }}\left(\beta \right),}\hfill \end{array}$$

(22)

where

$$\begin{array}{cc}\hfill p_{10}\left(\beta \right)& {\displaystyle =\frac{\beta /2}{\beta }=\frac{1}{2},}\hfill \\ \hfill p_{11}\left(\beta \right)& {\displaystyle =\frac{{(\beta /2)}^2}{\beta (\beta +1)}=\frac{\beta }{4(\beta +1)},}\hfill \\ \hfill p_{20}\left(\beta \right)& {\displaystyle =\frac{(\beta /2)(\beta /2+1)}{\beta (\beta +1)}=\frac{(\beta /2+1)}{2(\beta +1)},}\hfill \\ \hfill \cdots \\ \hfill p_{n{n}^{\prime }}\left(\beta \right)& {\displaystyle =\frac{(\beta /2)(\beta /2+1)\cdots (\beta /2+n-1)(\beta /2)(\beta /2+1)\cdots (\beta /2+n^{\prime }-1)}{\beta (\beta +1)\cdots (\beta +n+n^{\prime }-1)}.}\hfill \end{array}$$

(23)

The posterior for $\beta$ is independent from states x with just a single count, as $p_{10}\left(\beta \right)=\mathrm{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 (denoted as ${\beta }^{\ast }$) that maximizes $logp\left(\mathbf{n}\right|\beta )$ in the particular case where at most two samples coincide on the same x, obtaining

$$\frac{\partial }{\partial \beta }logp\left(\mathbf{n}\right|\beta )={\left.\frac{m_{11}}{\beta }+\frac{m_{20}}{\beta +2}-\frac{m_{11}+m_{20}}{\beta +1}\right|}_{{\beta }^{\ast }}=0.$$

(24)

Denoting the fraction of 2-count states that have one count for each y-value as $f_{11}=m_{11}/(m_{11}+m_{20})$, Equation (24) implies that the likelihood $p\left(\mathbf{n}\right|\beta )$ is maximized at

$${\beta }^{\ast }=\left\{\begin{array}{ccc}+\infty \hfill & \mathrm{when}\hfill & f_{11}\ge \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.,\hfill \\ \frac{f_{11}}{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.-f_{11}}\hfill & \mathrm{when}\hfill & f_{11}<\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right..\hfill \end{array}\right.$$

(25)

If the y-values are independent of x, we expect $f_{11}\sim \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$. This case corresponds to a large ${\beta }^{\ast }$ and, consequently, to little information. On the other side, for small $f_{11}$, the parameter ${\beta }^{\ast }$ is also small and the information grows. Moreover, the width of $p\left(\mathbf{n}\right|\beta )$ is also modulated by $f_{11}$. When the information is large, the peak around ${\beta }^{\ast }$ is narrow. Low information values, instead, require more evidence, and they come about with more uncertainty around ${\beta }^{\ast }$.

In Equation (24), 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\left(\beta \right)$ is included, the $\beta$ that maximizes the posterior $p\left(\beta \right|\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, Equation (20), focusing on the conditional entropy $\langle H_{Y|X}\rangle \left(\beta \right)$. First, in terms of the multiplicities, the conditional entropy can be rewritten as

$${\displaystyle \langle H_{Y|X}\rangle \left(\beta \right)=\sum _r{f}_r\sum _{n+n^{\prime }=r}{f}_{n{n}^{\prime }}{H}_{n{n}^{\prime }}\left(\beta \right),}$$

(26)

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

$${\displaystyle H_{n{n}^{\prime }}\left(\beta \right)={\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 }.}$$

(27)

A priori, $\langle I|\mathbf{n}=0,\beta \rangle =log2-H_{00}\left(\beta \right)$, as in Figure 1d. Surprisingly, from the property ${\psi }_0(z+1)={\psi }_0\left(z\right)+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,

$$\begin{array}{c}{\displaystyle \langle H_{Y|X}\rangle =f_1{H}_{00}\left(\beta \right)+f_2{f}_{11}{H}_{11}\left(\beta \right)+f_2{f}_{20}{H}_{20}\left(\beta \right)+\cdots }\hfill \end{array}$$

(28)

In the severely under-sampled regime, these first terms are the most important ones. Moreover, when evaluating these terms in $\beta ={\beta }^{\ast }$ (Equation (25)), the conditional entropy simplifies into

$$\langle H_{Y|X}\rangle =\left\{\begin{array}{ccc}log2\hfill & \mathrm{when}\hfill & f_{11}\ge \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.,\hfill \\ (f_1+f_2)H_{00}\left({\beta }^{\ast }\right)+\mathcal{O}\left(f_3\right)\hfill & \mathrm{when}\hfill & f_{11}<\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right..\hfill \end{array}\right.$$

(29)

Typically, $(f_1+f_2)$ takes most of the weight, so the estimation is close to the prior $H_{00}$ evaluated at the $\beta$-value that maximizes the marginal likelihood (or the posterior).

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 Appendix A. The variance (Equation (A4)) 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.

If an estimator, such as ours, is guaranteed to provide a non-negative estimate for all possible sets of sampled data, it cannot be free of bias, not at least if the samples are generated by an independent distribution ($q_{y|x}=q_y$), for which the true information vanishes. In Appendix B, we show that, in this specific case, the bias decreases with the square root of the number of coincidences. This number may be large, even in the severely undersampled regime, if $exp\left[H\right(X)/2]<N\ll exp\left[H\right(X\left)\right]$. If the distribution $q_x$ is approximately uniform, the number of coincidences is proportional to the square of the total number of samples, so the bias is inversely proportional to N.

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\left(H_{XY}\right)$. 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, Equation (17). 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 in 50 different sets of samples $\mathbf{n}$ of the measured data. We present the estimate obtained with $\langle I|\mathbf{n},{\beta }^{\ast }\rangle$ from Equation (20) evaluated in the ${\beta }^{\ast }$ that maximizes the marginal likelihood $p\left(\mathbf{n}\right|\beta )$. We did not observe any improvement when integrating over the whole posterior $p\left(\beta \right|\mathbf{n})$ with the prior $p\left(\beta \right)$ of Equation (19), 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 Figure 2, the performance of our estimator is compared with that of three other methods widely employed in the literature: Plug in, NSB and PYM. In addition, two other estimators were evaluated, but not shown in the figure to avoid cramming: the one of reference [21], which is a particularly convenient case of the Schuermann family of estimators [22], and the one of reference [23], extensively used in ecology. Their estimates fell between the plug-in estimator and NSB in the first case, and between NSB and PYM in the second case.

In the first example (Figure 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$ (as described in [14]). These values correspond to a Pitman–Yor prior with a heavy tail. The conditional probabilities $q_{y|x}$ are defined by sampling a symmetric Beta distribution $q_{y|x}\sim \mathrm{Beta}(\beta /2,\beta /2)$, as in Equation (17). In Figure 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\left(H_{XY}\right)\simeq 800$. We compare our estimator to the plug-in estimator (Plug-in), 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 in the Plug-in estimator, comparable to NSB and smaller than PYM. All the other methods (Plug-in, NSB and to a lesser extend PYM) overestimate the mutual information. In Figure 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 [25]), 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\left|\mathrm{R}\right(x)}$, where R is a rotation. Using a spherical harmonic representation [24], the frequency components ${\pi }_{\ell }\left(f\left(x\right)\right)$ of the spherical spectrum are obtained, where $f\left(x\right)$ 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\left(q_{y|x}\right)$ has some level of roughness (inset in Figure 2b), due to peaks coming from the invariant classes in $\{x_1,\cdots ,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\left(H_{XY}\right)\simeq 5000$. Our estimator has little bias (Figure 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 three counts, or maybe four. 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 (Equation (18)) at the ${\beta }^{\ast }$ that maximizes the marginal likelihood $p\left(\mathbf{n}\right|\beta )$, which in turn is mainly determined by $f_{11}$. In Figure 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-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.)+0.25\delta (q-q_0)+0.25\delta (q-1+q_0)]$, with $q_0=0.1$. The delta placed in $q=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$ 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}=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$, 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 Figure 2c,f). The marginal likelihood $p\left(\mathbf{n}\right|\beta )$ contains a single peak in an intermediate value of $\beta$, coinciding with none of the deltas in $p\left(q_{1|x}\right)$, 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}$, even on average.

Finally, we check numerically the accuracy of the analytically predicted mean posterior information (Equation (20)) and variance (Equation (A4)) in the severely under-sampled regime. The test is performed in a different spirit than the numerical evaluations of Figure 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 Equations (20) and (A4), however, must be interpreted in the Bayesian sense. The square brackets in $\langle I|\mathbf{n}\rangle$ and $\langle H_{Y|X}^2\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 \mathrm{DP}\left(\alpha \right)$ (a Dirichlet Process with concentration parameter $\alpha$) and $q_{y|x}\sim \mathrm{Beta}(\beta /2,\beta /2)$. A total of 13,500 distributions $q_{xy}$ are produced, with $log\beta$ sampled from Equation (19), 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,\left\{{\mathbf{q}}_{y|x}\right\},I({\mathbf{q}}_x,\left\{{\mathbf{q}}_{y|x}\right\}),\mathbf{n},\langle I|\mathbf{n}\rangle$ and ${\sigma }^2\left(I\right|\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 $\left\{m_{n{n}^{\prime }}\right\}$ —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},\left\{{\mathbf{q}}_{y|x}\right\})$ of the corresponding class, and compare them with our predicted estimates $\langle I|\left\{m_{n{n}^{\prime }}\right\}\rangle$ and $\langle {\sigma }_I^2|\left\{m_{n{n}^{\prime }}\right\}\rangle$, using the prior $p(log\beta )$ from Equation (19). 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 $\left\{{\mathbf{q}}_{y|x}\right\}$, 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 ${\widehat{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 Equation (13) by another prior that depends on both ${\mathbf{q}}_x$ and $\left\{{\mathbf{q}}_{y|x}\right\}$.

The simplest hypothesis is to assume that the prior $p({\mathbf{q}}_x,\left\{{\mathbf{q}}_{y|x}\right\})$ factorizes as $p\left({\mathbf{q}}_x\right)p\left(\left\{{\mathbf{q}}_{y|x}\right\}\right)$, implying that the marginal probabilities $q_x$ are independent of the conditional probabilities $q_{y|x}$. We propose $q_x\sim \mathrm{DP}\left(\alpha \right)$, so that the marginal probabilities $q_x$ are drawn from a Dirichlet Process with concentration parameter $\alpha$, associated with 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

$$\begin{array}{ccc}\hfill \langle I|\mathbf{n},\alpha ,\beta \rangle & =& \frac{N}{N+\alpha }\left\{\widehat{H}\left(Y\right)-\sum _{x,n_x>0}\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]\right\}\hfill \\ & & +\frac{\alpha }{N+\alpha }\left[\widehat{H}\left(Y\right)-{\psi }_0(\beta +1)+\sum _y\widehat{q_y}{\psi }_0(\beta \widehat{q_y}+1)\right].\hfill \end{array}$$

(30)

Before including the prior $p\left({\mathbf{q}}_x\right)$, in the severely undersampled regime, the mean posterior information was approximately equal to the prior information evaluated in the best $\beta$ (Equation (16)). The new calculation (Equation (30)) 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 $\left\{{\mathbf{q}}_{y|x}\right\}$ implies that

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

(31)

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

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

(32)

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\left({\mathbf{n}}_x\right|\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 the Plug-in estimator.

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 (11) 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 (Equation (13)), which is designed to capture the same idea, and, in addition, to be analytically tractable. We choose to work with a hyper-parameter $\beta$ that regulates the scatter of $q_{y|x}$ around $q_y$, and not the scatter of $q_{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\left(\mathbf{n}\right|\beta )$.

Although our proposal is inspired in previous Bayesian studies, the procedure described here is not strictly Bayesian, since our prior (Equation (13)) requires the knowledge of ${\widehat{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 ${\widehat{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 Equation (13) by maximizing $p\left(q_y\right|\mathbf{n})$, and then replace ${\widehat{q}}_y$ in Equation (13) 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\left(\left\{{\mathbf{q}}_{y|x}\right\}\right)$ to those generated by Equation (13), 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 [26], and that is shared also by other methods [8,13,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 [26]. 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.

In the explored examples, our estimator had a small bias, even in the severely under-sampled regime, and it outperformed other estimators discussed in the literature. More importantly, the second and third examples of Section 5 showed that it even worked when the collection of true conditional probabilities $q_{y|x}$ was not contained in the family of priors generated by $p\left(q_{y|x}\right|\beta )$. In these cases, the success of the method relies on the peaked nature of the posterior distribution $p\left(\beta \right|\mathbf{n})$. Even if the selected $p\left(q_{y|x}\right|\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.

When the number of samples N is much larger than the total number of available states $k_x\times k_y$, our estimator of mutual information tends to the plug-in estimator, which is known to be consistent [20]. Consequently, our estimator is also consistent. By construction, $〈I|\mathbf{n}〉$ is equal to the mutual information averaged over all distributions ${\mathbf{q}}_x$ and $\left\{{\mathbf{q}}_{y|x}\right\}$ compatible with the measured data, each weighted by its posterior probability. As such, it can never produce a negative result, which is a desirable property. The down side is that the estimator must be positively biased, at least, when the true information vanishes. The derivation in Appendix B shows that, when the number of samples $N\in [{\mathrm{e}}^{H\left(X\right)/2},{\mathrm{e}}^{H\left(X\right)}]$, this bias is inversely proportional to the square root of the number of pairwise coincidences or, when $q_x$ is fairly uniform, to the inverse of the total number of samples. Moreover, the factor of proportionality is significantly smaller than the one obtained in the bound of the bias of other frequentist estimators [3,27]. If the number of samples N grows even further, the bias tends to zero, since the bias of all consistent estimators vanishes asymptotically [28].

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\left(\beta \right|\mathbf{n})$, as indicated by Equation (14), and only shift the posterior information towards the prior (Equation (30)). More interestingly, the x-states with just a single count are also irrelevant, both in shaping $p\left(\beta \right|\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\left(\beta \right|\mathbf{n})$, and in deviating the posterior information away from the prior.

Our method can also be extended to generalizations of mutual information designed to characterize the degree of interdependence of more than two variables [29,30,31,32,33,34,35,36,37,38,39,40,41], as long as all but one of the variables are extensively sampled. The applicability of the method to these measures will be the object of future work, once a certain degree of consensus has been reached regarding their meaning and range of applicability.

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,42,43], or in genomics, to understand associations between genes (large-entropy variable) and a few specific phenotypes [44]. The method can also shed light on the development of rate-distortion methods to be employed in situations in which only a few samples are available. In particular, it can be applied within the information bottleneck framework [45,46], aimed at extracting a maximally compressed representation of an input variable, but still preserving those features that are relevant for the prediction of an output variable. The possibility of detecting statistical dependencies 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 [47].

8. Conclusions

We have proposed a novel estimator for the mutual information $I(X;Y)$ between two variables, applicable to those cases in which the marginal distribution of one of the 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 a consistent estimator that presents very low bias, outperforming previous methods discussed in the literature. The main data statistics determining the estimated value is the inhomogeneity of the conditional distribution of the low-entropy variable in those states in which the large-entropy variable registers coincidences.