# Nearest-Neighbours Estimators for Conditional Mutual Information

Jake Witter  Conor Houghton  School of Engineering Mathematics and Technology
University of Bristol, Bristol, United Kingdom
(jake.witter /<EMAIL_ADDRESS>
###### Abstract

The conditional mutual information quantifies the conditional dependence of
two random variables. It has numerous applications; it forms, for example,
part of the definition of transfer entropy, a common measure of the causal
relationship between time series. It does, however, require a lot of data to
estimate accurately and suffers the curse of dimensionality, limiting its
application in machine learning and data science. However, the Kozachenko-
Leonenko approach can address this problem: it is possible, in this approach
to define a nearest-neighbour estimator which depends only on the distance
between data points and not on the dimension of the data. Furthermore, the
bias can be calculated analytically for this estimator. Here this estimator is
described and is tested on simulated data.

††preprint: APS/123-QED

##  I Introduction

The mutual information is a fundamental and common approach to assessing the
relationship between random variables. It calculates the amount of information
the value of one random variable contains, on average, about the value of
another. The conditional mutual information is a natural extension of this, it
is the average information in one variable about another, given knowledge of a
third.

The conditional mutual information has many applications. In particular it is
used to calculate the transfer entropy [1, 2]; the transfer entropy is a
powerful approach to measuring causality which extends in a model-free manner
the insight provided by the autoregressive approach to commonly-used Granger
causality [3]. Conditional mutual information is also used to calculate the
interaction information, the natural generalization of mutual information to
more than two random variables [4, 5].

Mutual information is often difficult to calculate in a model-free way: there
is typically a huge data requirement when estimating it and it suffers from
the curse of dimensionality, with the data required increasing rapidly as the
dimension increases. This often renders conditional mutual information and
transfer entropy impractical for problems involving the sort of
multidimensional variables typically encountered in scientific applications
and in data science and machine learning. In Tobin and Houghton [6], Houghton
[7] a Kozachenko-Leonenko (KL) [8, 9] estimator is presented, this gives a
method for estimating mutual information which is defined on metric spaces.
Since it uses only the distance between data points it wards off the curse of
dimensionality. Here, this metric-space KL approach is extended to estimating
conditional mutual information.

This gives a formula for mutual information similar in spirit to the one
derived in Kraskov et al. [9]. They are both nearest neighbour formulae
derived using the KL estimator. The estimator described in Kraskov et al. [9]
has been adapted to conditional mutual information [10, 11] and transfer
entropy [12]; here these are referred to as the ‘KSG’ estimators. The novel
estimator that is introduced here is referred to as the ‘new’ estimator.

There are two significant differences between the KSG and new estimators; the
first is that the derivation of the new estimator relies only on metric space
properties of the data. Secondly, the new estimator allows the bias to be
calculated and the bias-correction gives a useful approach to finding a value
for $h$, the number of nearest neighbours. Although we believe these
differences are both advantages that the new estimator has over the KSG
estimator, one of the aims in this paper is to use simulated data to compare
different approaches to estimating conditional mutual information ‘in action’.

##  II Methods

The conditional mutual information between variables $X$ and $Y$, conditioned
on $Z$, is given by

$I(X,Y|Z)=\int_{\mathcal{X}}dx\int_{\mathcal{Y}}dy\int_{\mathcal{Z}}dz\log\left[\frac{p_{XY|Z}(x,y|z)}{p_{X|Z}(x|z)p_{Y|Z}(y|z)}\right]p_{XYZ}(x,y,z)$
(1)

where $X$, $Y$ and $Z$ take values in $\mathcal{X}$, $\mathcal{Y}$ and
$\mathcal{Z}$. Writing the integral as an expectation value and expanding the
probability mass function to give an equivalent form that proves more
convenient:

$I(X,Y|Z)=\left\langle\log\frac{p_{XYZ}(x,y,z)p_{Z}(z)}{p_{XZ}(x,z)p_{YZ}(y,z)}\right\rangle_{p_{XYZ}(x,y,z)}$
(2)

The challenge with calculating this quantity lies in estimating the
probabilities. The histogramme method estimates these probability mass
functions by constructing histograms and counting the data points in each bin
to estimate the likelihood of the corresponding outcomes. This will give the
correct answer in the infinite limit of available data and infinitely small
bin sizes, but with three random variables the number of bins needed to
produce an accurate answer is often far larger than what the quantity of data
will support. Of course, this problem becomes more acute if $X$, $Y$ or $Z$ is
a vector; the required number of points increases as a power of the number of
dimensions. Similarly, in transfer entropy, the random variables $X$ and $Z$
are vectors over recent time steps, making high-dimensionality unavoidable.

###  II.1 Estimator

Here the KL approach is described. It is intended here to compare both the KSG
and new estimators to a ground truth. By using simulated data, it is possible
to generate enough data to accurately calculate conditional information using
histogram-based methods. The new estimator is described in detail but there is
only a short outline of the KSG estimator since the details can be found in
Kraskov et al. [9], Frenzel and Pompe [10], Lizier [11], Gómez-Herrero et al.
[12].

A KL approach exploits the proximity structure of the data. Typically this
proximity structure is based on coordinate distances but for the new estimator
it is only assumed that the random variables, $X$, $Y$ and $Z$ take values in
metric spaces $\mathcal{X}$, $\mathcal{Y}$ and $\mathcal{Z}$. Let

$\mathcal{D}=\\{(x_{1},y_{1},z_{1}),(x_{2},y_{2},z_{2}),\ ...\
,(x_{n},y_{n},z_{n})\\}$ (3)

denote the data. Depending on the problem being modelled the outcomes $x_{i}$,
$y_{i}$ or $z_{i}$ could be high-dimensional vectors, or even values in a
space without a manifold structure. However it is assumed that it is possible
to calculate distances on the three spaces. Furthermore each set of triples
$(x_{i},y_{i},z_{i})$ is assumed to be drawn from $p_{XYZ}(x,y,z)$. Following
from the definition of the conditional mutual information, Eq. 2, $I(X,Y|Z)$
can be approximated, for suitably large $n$, by the Monte-Carlo estimate

$I(X,Y|Z)\approx\frac{1}{n}\sum_{i=1}^{n}\log\left[\frac{p_{XYZ}(x_{i},y_{i},z_{i})p_{Z}(z_{i})}{p_{XZ}(x_{i},z_{i})p_{YZ}(y_{i},z_{i})}\right].$
(4)

This reduces the problem to estimating $p_{XYZ}(x_{i},y_{i},z_{i})$,
$p_{XZ}(x_{i},z_{i})$, $p_{YZ}(y_{i},z_{i})$ and $p_{Z}(z_{i})$ at the data
points in $\mathcal{D}$.

In the KL approach the probability at a point is approximated by counting the
number of points in a nearby region. For simplicity consider temporarily a
single random variable $V$ with data $\\{v_{1},v_{2},\ldots,v_{n}\\}$. The KL
estimate of $p_{V}(v_{i})$ is

$p_{V}(v_{i})\approx\frac{h_{V}(v_{i})}{n\operatorname{vol}{[B_{V}(v_{i})]}}$
(5)

where $B_{V}(v_{i})$ is a region in $\mathcal{V}$ around $v_{i}$,
$h_{V}(v_{i})=|B_{V}(v_{i})|$ is the number of data points in the region and
$\operatorname{vol}{[B_{V}(v_{i})]}$ is the volume of the region
$B_{V}(v_{i})$. The estimate follows from the definition of the probability
mass $p_{V}(v)$:

$\langle h_{V}(v_{i})\rangle=n\int_{B_{V}(v_{i})}dv\,p_{V}(v)$ (6)

by using two approximations. Firstly approximating $p_{V}(v)$ by assuming it
is constant throughout $B_{V}(v_{i})$ and secondly, approximately evaluating
the integral by counting the data points: $h_{V}(v_{i})\approx\langle
h_{V}(v_{i})\rangle$. Clearly these two approximations favour different sizes
of $B_{V}(v_{i})$: approximating $p_{V}(v)$ as constant favours a small
region, approximating the expected number of data points by the actual size of
the sample found in $B_{V}(v_{i})$ favours a large region.

The other important ingredient is the calculation of the volume
$\operatorname{vol}{B_{V}(v_{i})}$. Since this is a metric space approach, to
avoid any assumptions about a coordinate structure, this calculation follows
[7] and calculates volumes using probability densities. In essence, the volume
of a region is given by the probability mass it contains. This is easier to
discuss in the context where there are different probability densities, as in
the problem being considered here where the goal is to estimate conditional
information and this requires the estimate of a number of different
probability mass functions: $p_{XYZ}(x_{i},y_{i},z_{i})$,
$p_{XZ}(x_{i},z_{i})$, $p_{YZ}(y_{i},z_{i})$ and $p_{Z}(z_{i})$. In fact, the
key to estimating the volume is to use different distributions on the same
space. This is required because using the same distribution to estimate volume
and probability mass trivializes the estimate: to see this consider, for
example, $p_{Z}(z_{i})$. To define the region $B_{Z}(z_{i})$ a volume is
fixed, say $h/n$; where $h$ counts the number of points and so can be thought
of as a smoothing parameter. Now, for each point $z_{i}$ a ball is expanded
around this seed point until it has volume $h/n$, since the probability mass
is being used to measure the volume, this is estimated by counting data
points:

$\operatorname{vol}{B_{z}(i)}=\frac{1}{n}\int_{B_{Z(i)}}dz\,p_{Z}(z)=\frac{\langle
h_{Z}(i)\rangle}{n}\approx\frac{h_{Z}(i)}{n}$ (7)

where, for brevity, the shorthand $B_{Z}(i)=B_{Z}(z_{i})$ and
$h_{Z}(i)=h_{Z}(z_{i})$ has been used. Thus, using the metric, a ball in
$\mathcal{Z}$ is expanded around the seed point until it contains $h$ points,
including $z_{i}$ itself. Now, using the KL estimate of the probability:

$p_{Z}(z_{i})\approx\frac{h_{Z}(i)}{n\operatorname{vol}{B_{Z}(i)}}=1,$ (8)

and thus gives a trivial estimate.

To avoid this problem a different distribution is used to estimate volumes.
For joint distributions this is provided by the marginalized distribution.
First, consider $p_{XZ}(x_{i},z_{i})$; in this case $p_{X}(x_{i})p_{Z}(z_{i})$
can be used to measure the volume. At first this is confusing, the volume of
the ball around $(x_{i},z_{i})$ depends on the number of $x_{j}$ points and
the number of $z_{j}$ points, not on the number of $(x_{j},z_{j})$ pairs, as
it would if $p_{XZ}(x,z)$ was being used as the volume measure. As such, the
ball in $\mathcal{X}\times\mathcal{Z}$ is formed by expanding a ball around
$x_{i}$ until it contains $h$ points and a ball around $z_{i}$, again, until
it also contains $h$ points. The ball $B_{XZ}(x_{i},z_{i})$ is then the
intersection of these two balls; it has volume estimated as $h^{2}/n^{2}$.
Writing $h_{XZ}(i)$ as the number of points it contains, the KL estimate is

$p_{XZ}(x_{i},z_{i})\approx\frac{nh_{XZ}(x_{i},z_{i})}{h^{2}}.$ (9)

In the case of $p_{XYZ}(x,y,z)$ the alternative distribution on the same space
is again provided by the marginal distribution, $p_{X}(x)p_{Y}(y)p_{Z}(z)$ and
the ball $B_{XYZ}(x_{i},y_{i},z_{i})$ is formed the same way as for the
$p_{XZ}(x,z)$ example. It is the intersection of $B_{X}(i)$, $B_{Y}(i)$ and
$B_{Z}(i)$ and has volume $h^{3}/n^{3}$. This gives

$p_{XYZ}(x_{i},y_{i},z_{i})\approx\frac{n^{2}h_{XYZ}(i)}{h^{3}}$ (10)

where $h_{XYZ}(i)$ is the number of points in $B_{XYZ}(x_{i},y_{i},z_{i})$.
The regions have been referred to as balls, though, in some ways, they are
more like squares or cubes; Fig. 1 provides an illustration. From now on, the
shorthand introduced earlier is extended to these regions, so
$B_{XYZ}(i)=B_{XYZ}(x_{i},y_{i},z_{i})$ and so on.

Figure 1: An illustration of the different intersections. In this cartoon the
three sets of outcomes, $\mathcal{X}$, $\mathcal{Y}$ and $\mathcal{Z}$, are
marked $X$, $Y$ and $Z$; here these are one-dimensional whereas, of course,
they can be high-dimensional or even spaces without a dimension. In A the
three balls around the seed point, $B_{X}$, $B_{Y}$ and $B_{Z}$, are drawn in
blue, yellow and red. Each of these balls contain $h$ points, this determines
their radius. In B attention is restricted $B_{Z}$; the intersection this has
with $B_{X}$ gives $B_{XZ}$ and with $B_{Y}$ gives $B_{YZ}$. All three
intersect to give the cyan region $B_{XYZ}$. C shows $B_{Z}$ again, from
above; the points are also marked, with a star for the seed point. Here,
$h_{XYZ}=4$, $h_{XZ}=6$, $h_{YZ}=9$ and all the points lie in $B_{Z}$ so $h$
in this illustration is 22.

Using these point probability estimates, the conditional mutual information,
Eq. 4, can be rewritten as

$\begin{split}I(X,Y|Z)\approx\frac{1}{n}\sum_{i=1}^{n}\log\frac{h_{XYZ}(i)h}{h_{XZ}(i)h_{YZ}(i)}.\end{split}$
(11)

Thus, the conditional mutual information is estimated by finding the nearest-
neighbours for each data point in each of $\mathcal{X}$, $\mathcal{Y}$ and
$\mathcal{Z}$ and then counting the points in the intersections $B_{XZ}(i)$,
$B_{YZ}(i)$ and $B_{XYZ}(i)$. One complication arises when there are draws;
this is discussed in Appendix 1.

###  II.2 Bias

As is typical with estimators of mutual information, this approach has an
upward bias. That is, when applying the approach to conditionally independent
$X$ and $Y$:

$X\rightarrow Z\rightarrow Y$ (12)

there will be a non-zero estimated conditional mutual information, when, of
course, the true value is zero. To remove this bias the expected value of the
estimator is calculated, assuming $X$ and $Y$ are conditionally independent,
given $Z$. This follows similar calculations in [13, 14] but is more
complicated in this case because of the additional random variable.

Calculating the bias requires a formula for the probability that the
intersection $B_{XYZ}(i)$ contains $r$ points given the values $h_{XZ}(i)$ and
$h_{YZ}(i)$, that is the sizes of $B_{XZ}(i)$ and $B_{YZ}(i)$, something that
will vary from point to point. This bias therefore varies for each data point
triple, each $i$, unlike similar formulae for the bias of the estimator in
[13, 14]. The reason for this difference is that the independence is
conditional and so the bias has to include the dependence that $X$ and $Y$
might have while calculating the expected size of $B_{XYZ}(i)$ under the
assumption of conditional independence. Thus, the bias includes a sum over all
data points, for $(x_{i},y_{i},z_{i})$ the contribution to the bias is

$I_{b}(i,h)=\sum^{h}_{r=1}\mathbb{P}(h_{XYZ}(i)=r)\log\left[\frac{rh}{h_{XZ}(i)h_{YZ}(i)}\right].$
(13)

In order to calculate $\mathbb{P}(h_{XYZ}(i)=r)$ consider the $h_{XZ}(i)$
points in $B_{XZ}(i)$; these points are chosen from the $h$ points in
$B_{Z}(i)$. The number of these points that are also among the $h_{YZ}(i)$
points in $B_{YZ}(i)$ gives $r$. This can be thought of as an urn problem.
First remove the seed point $(x_{i},y_{i},z_{i})$ itself, now there are $h-1$
balls in an urn, $h_{XZ}(i)-1$ of which are ‘red’, meaning they are in
$B_{XZ}(i)$. Next $h_{YZ}(i)-1$ balls are drawn out of that urn, representing
the number in $B_{YZ}(i)$, $r-1$ is the number of these balls that are ‘red’.
This is described by the hypergeometric distribution, so the probability
$p(r)$ that $h_{XYZ}(i)=r$ is

$p(r)=\frac{\left(\begin{array}[]{c}h_{XZ}(i)-1\\\
r-1\end{array}\right)\left(\begin{array}[]{c}h-h_{XZ}(i)\\\
h_{YZ}(i)-r\end{array}\right)}{\left(\begin{array}[]{c}h-1\\\
h_{YZ}(i)-1\end{array}\right)}\\\ $ (14)

or

$r\sim\mbox{Hypergeometric}(h-1,h_{XZ}(i)-1,h_{YZ}(i)-1).$ (15)

Now, it is possible to calculate $I_{b}(i,h)$ and hence the total bias

$I_{b}(h)=\frac{1}{n}\sum_{i=1}^{n}I_{b}(i,h)$ (16)

It is clear that the result of the estimator will be a function of $h$. As
mentioned above, a larger $h$ will mean a larger region, increasing the error
introduced by assuming the probability density is constant throughout the
region, but, conversely, makes it more likely that the regions contain a
representative number of points. While it is unclear how to choose $h$,
practically, maximising the bias-corrected information over $h$ works well and
gives a good estimate of the information, as in [13]; it balances the two
sources of error, both of which might be expected to reduce the de-biased
estimate. Hence

$I(X,Y|Z)\approx\max_{h}\frac{1}{n}\sum_{i=1}^{n}\left[\log\frac{h_{XYZ}(i)h}{h_{XZ}(i)h_{YZ}(i)}-I_{b}(i,h)\right]$
(17)

Thus, the procedure is to maximize the bias-corrected mutual information over
$h$; this is performed here using golden section search [15]. This procedure
does cause an upward bias when the mutual information is low: when there is
little mutual information the noise in the estimate is larger compared to the
true value and so maximizing over $h$ effectively selects the value of $h$
whose corresponding estimate is largest because of the noise, rather than
because this is the optimal value of $h$. However, for the more typical case,
where the estimator is being used with variables that are not conditionally
independent, the procedure appears to work well.

###  II.3 The KSG estimator

The KSG estimator, [9] also adopts a KL approach in that it also estimates the
probabilities by counting data points in small regions. The most obvious
difference between the KSG estimator and the new estimator lies in the
calculation of volumes. In the new estimator, as described above, the volume
is also estimated by counting data points whereas the KSG estimator assumes a
coordinate system and calculates volumes by integration. There is a second
difference, in the KSG estimator the probability density is estimated by
calculating the probability of finding the $k$th point at a particular
distance, a more sophisticated approach to the one used in the new estimator,
but one that depends on the space having a coordinate system relevant to the
calculation of mutual information. Ultimately this leads to a formula that
contains a digamma function, [16], rather than a logarithm, the digamma
function is typically denoted $\psi(x)$.

There are two, similar, versions of the KSG estimator:

$\displaystyle I(X,Y)\approx$ $\displaystyle\psi(k)+\psi(n)$ (18)
$\displaystyle-\frac{1}{n}\sum_{i}[\psi(k_{X}(i)+1)+\psi(k_{Y}(i)+1)]$

and

$\displaystyle I(X,Y)\approx$ $\displaystyle\psi(k)-\frac{1}{k}+\psi(n)$ (19)
$\displaystyle-\frac{1}{n}\sum_{i}[\psi(k_{X}(i))+\psi(k_{Y}(i))]$

where $k$ is a smoothing parameter, similar to $h$ in the new estimator and
$n$ is the total number of points. Like the new estimator, the formulae
involve counting points, with $k_{X}(i)$ and $k_{Y}(i)$ giving the number of
points in particular regions. They are, however, potentially confusing in that
there is actually a subtly different definition used for each of the two
formula. The definition is also different from the one used for quantities
such as $h_{XY}(i)$ in the new estimator. $k_{X}(i)$ and $k_{Y}(i)$ are what
could be thought of as “outer” counts whereas $h_{XY}(i)$ is an “inner” count.
Loosely, in an inner count two regions, $B_{X}(i)$ and $B_{Y}(i)$, are given
with specified size $h$ and the number of points in the intersection,
$h_{XY}(i)$, is calculated. In an outer count this is reversed: the size of an
intersection, $k$, is specified and counting is performed on the two
corresponding regions. More specifically, to form the outer count used in the
KSG estimator, the distance in the product space $\cal{X}\times\cal{Y}$ from a
seed point $(x_{i},y_{i})$ to its $k$th nearest neighbour is found. This
defines the intersection. The subtly arises because this requires a metric on
$\cal{X}\times\cal{Y}$ when the metrics are only defined on $\cal{X}$ and
$\cal{Y}$ separately. To lift these metrics to a metric on the product, the
maximum is used, with the obvious notation

$d_{XY}[(x_{i},y_{i}),(x_{j},y_{j})]=\text{max}[d_{X}(x_{i},x_{j}),d_{Y}(y_{i},y_{j})].$
(20)

Let $d^{*}$ denote the distance to the $k$th nearest point. There are now two
options, the first uses this distance to define $B_{X}(i)$ and $B_{Y}(i)$, so,
for example, $B_{X}(i)$ is the set of all points $x_{j}$ within $d^{*}$ of the
seed point $x_{i}$, that is $B_{X}(i)=\\{x_{j}:d_{X}(x_{i},x_{j})<d^{*}\\}$.
The second does something more complicated, it uses $d^{*}_{X}$ to define
$B_{X}(i)$ and $d^{*}_{Y}$ to define $B_{Y}(i)$ where
$d^{*}=\text{max}[d_{X}^{*},d_{Y}^{*}]$. In either case, the outer counts are
then $k_{X}(i)=|B_{X}(i)|$ and $k_{Y}(i)=|B_{Y}(i)|$. The two KSG estimators
above, Eq. 18 and Eq. 19, correspond to these two options.

It is interesting that these formula do not include terms that depend on the
dimension and in this way, although the derivation of the KSG estimator relies
on the existence of coordinates, the formula itself does not. Beyond the
theoretical objection that these estimators are not derived for spaces without
coordinates, even if they can be applied to spaces of that type, the use of
the maximum metric on $\mathcal{X}\times\mathcal{Y}$ creates practical
difficulties in using these estimators relative to using the new estimator.
Most obviously it relies on some compatibility of the metrics on $\mathcal{X}$
or $\mathcal{Y}$ since they are compared when calculating the distance on
$\mathcal{X}\times\mathcal{Y}$. There is a milder version of this issue in the
case of the new estimator where the choice of the same volume for both
$B_{X}(i)$ and $B_{Y}(i)$ also makes an implicit assumption of similarity. The
other challenge is that it is more difficult to calculate the bias. A
consequence of using the maximum metric is that counts $k_{X}(i)$ and
$k_{Y}(i)$ depend on the distributions of points in $\cal{X}$ and $\cal{Y}$
even if these two distributions are independent, preventing the calculation of
the bias without making additional assumptions.

These estimators are in turn used to derive estimators of the conditional
mutual information in [10, 11] and the transfer entropy in [12]. These
estimators, like the new estimator proposed here, are similar to the mutual
information estimators but contain more complicated counts of data points. The
two estimators are

$I(X,Y|Z)\approx\psi(k)+\frac{1}{n}\sum_{i}[\psi(k_{Z}+1)-\psi(k_{XZ}(i)+1)-\psi(k_{YZ}(i)+1)]$
(21) $\displaystyle I(X,Y|Z)\approx$
$\displaystyle\psi(k)-\frac{2}{k}+\frac{1}{n}\sum_{i}\left[\psi(k_{Z}(i))-\psi(k_{XZ}(i))-\frac{1}{k_{XZ}(i)}-\psi(k_{YZ}(i))-\frac{1}{k_{YZ}(i)}\right]$
(22)

where, again, variables like $k_{XZ}$ represents an outer count. In short the
$B_{X}(i)$, $B_{Y}(i)$ and $B_{Z}(i)$ balls are expanded to include $k$
points, $k_{Z}(i)=|B_{Z}(i)|$ whereas $k_{XZ}(i)$ is the number of $(x,z)$
data points in $B_{X}(i)\times B_{Y}(i)$.

These estimators differ from the new estimator in relying on coordinates in
their derivation, if not their implementation; the derivation is, in general,
more subtle and the estimation of the bias for zero information is not
provided, nor is there an approach to selecting a value of $K$.

###  II.4 Transfer Entropy

One important application of this new estimator for conditional mutual
information is the calculation of transfer entropy: a directed measure of
information transfer between two processes [2]. This quantity can be useful as
it gives a direction to the relationship between processes, as well as an
insight in to whether a relationship is causal. Consider $X_{t}$, $Y_{t}$ to
be two random processes with $t\in\mathbb{N}$ a discrete time index. The
transfer entropy from process $X$ to process $Y$ is the mutual information
between past values of $X$ and the present value of $Y$, conditioned on the
past values of $Y$. That is, as described in [17],

$T(X\rightarrow Y)=I(Y_{t},X_{t-1:t-\ell}|Y_{t-1:t-\ell})$ (23)

where $X_{t-1:t-\ell}=(X_{t-1},X_{t-2},\ldots,X_{t-\ell})$ is the past $\ell$
values of $X_{t}$, with a similar notation for $Y_{t-1:t-\ell}$. Thus, roughly
speaking, the transfer entropy measures the extra information about the
present value of $Y$ contained in the past values of $X$ beyond what is
already accounted for in the past of $Y$. Transfer entropy is widely used in
neuroscience, see for example [18, 19] and finance, see [20], but the data
requirement for estimating conditional mutual information is an impediment to
its practical use. The transfer entropy can be thought of as a non-parametric
version of the commonly-used Granger causality; or, rather, transfer entropy
can be considered the ideal approach to estimating causal relationship when
there is no data-constraint, with the Granger causality interpreted as an
approximation to the transfer entropy which is exact for Gaussian random
variables [3].

##  III Example

###  III.1 A simple Markov tree

To demonstrate the estimation of conditional mutual information and to compare
approaches, a simple example is considered here. This example is designed to
produce data with a conditional mutual information that can be varied by
changing a parameter. It is also designed to allow for the easy generation of
a large quantity of data, which means that histogram methods can be used to
give something approaching a ground truth. The example consists of four
Gaussian random variables; $W$, $X$, $Y$ and $Z$, distributed as

$\displaystyle w_{a}$ $\displaystyle\sim\text{Normal}(0,\sigma_{w})$ (24)
$\displaystyle x_{a}$ $\displaystyle\sim\text{Normal}(w_{a},\sigma_{x})$ (25)
$\displaystyle y_{a}$ $\displaystyle\sim\text{Normal}(w_{a},\sigma_{y})$ (26)
$\displaystyle z_{a}$ $\displaystyle\sim\text{Normal}(w_{a},\sigma_{z})$ (27)

where the $a$-index corresponds to the different dimensions of the data, so,
for example, for two dimensions $X=(X_{1},X_{2})$. Here one- and two-
dimensional examples will be considered. This relationship is illustrated in
Fig. 2.

$z_{a}$$w_{a}$$x_{a}$$y_{a}$ Figure 2: The Markov network for our simulated
data. This illustrates the relationship between the three observed variables
$X$, $Y$ and $Z$ and the unobserved variable $W$. $X$ and $Y$ are
conditionally independent given $W$; since $Z$ is a noisy version of $W$ the
mutual independence of $X$ and $Y$ given $Z$ depends on the amount of that
noise, $\sigma_{z}$.

Clearly $X$ and $Y$ are related, $I(X,Y)>0$, but they are conditionally
independent when conditioned on $W$, so $I(X,Y|W)=0$. However, here $I(X,Y|Z)$
is estimated and $Z$ is not $W$, in fact $Z$ can be thought of as a noisy
version of $W$. Thus, $I(X,Y|Z)$ is non-zero and will depend on how noisy $Z$
is, as determined by $\sigma_{z}$. As $\sigma_{z}$ approaches zero, $Z$
becomes closer to an exact copy of $W$ and $I(X,Y|Z)$ approaches zero.
Conversely, as $\sigma_{z}$ becomes large, conditioning on $Z$ becomes
increasingly meaningless and $I(X,Y|Z)$ approaches $I(X,Y)$.

\begin{overpic}[width=303.53267pt,tics=10]{figs/example_results_figure.pdf}
\put(0.0,98.0){\huge{a}} \put(0.0,58.0){\huge{b}} \put(0.0,28.0){\huge{c}}
\end{overpic} Figure 3: Results from one-dimensional Markov tree. In a the
relationship between the estimated conditional mutual information and
$\sigma_{z}$ is shown. The shaded area shows the middle 50% of estimates.
Here, the KSG and new methods use 3500 points, while the histogram uses
$5\times 10^{6}$ points to establish the ground truth. b shows how this
estimate scales over the number of data points used in the estimate. Again,
the histogram method uses $5\times 10^{6}$ points. In c, the relationship
between the $k$ used in the KSG estimates, and the estimated information is
shown. This demonstrates that choice of $k$ does strongly influence the
estimated information. Note the change in scale in the vertical axis between a
and the other two rows.
\begin{overpic}[width=303.53267pt,tics=10]{figs/example_results_figure2d.pdf}
\put(0.0,98.0){\huge{a}} \put(0.0,58.0){\huge{b}} \put(0.0,28.0){\huge{c}}
\end{overpic} Figure 4: Results from two dimensional Markov tree In a the
relationship between the estimated conditional mutual information and
$\sigma_{z}$ is shown. The shaded area shows the middle 50% of estimates.
Here, the KSG and new methods use 3500 points, while the histogram uses
$5\times 10^{6}$ points to establish the ground truth. b shows how this
estimate scales over the number of data points used in the estimate. Again,
the histogram method uses $5\times 10^{6}$ points. In c, the relationship
between the $k$ used in the KSG estimates, and the estimated information is
shown. This demonstrates that choice of $k$ does strongly influence the
estimated information. Note the change in scale in the vertical axis between c
and the other two rows.

The estimator is applied to the example in one and two dimensions, as seen in
Fig. 3 and Fig. 4. A histogram approach was also applied; this gives something
akin to a ground truth.

Considering first the one-dimensional case, Fig. 3a graphs the mutual
information as $\sigma_{z}$ varies, the histogram estimate and a numerical
show show the mutual information increasing, as expected, as $\sigma_{z}$
increased. For the estimators a fixed $n$ of 3500; while all three give a
reasonable estimate, the new estimator is closer to the correct value. In Fig.
3b the behaviour of the estimators for differing value of $n$ is considered
for three different values of $\sigma_{z}$, this shows that all three
estimators appear to require about 3500 samples to reach their best estimate.
The two KSG estimators require a choice of smoothing parameter $k$, in Fig. 3c
it is seen that this choice does affect the estimate for KSG2 and the
estimated value varies from too large to too small as $k$ increases; KSG1
consistently gives too small a value.

The same set of measures are exhibited for two dimension in Fig. 4. Here, the
results are somewhat similar, but more extreme. That is, the advantage of the
new estimator becomes more pronounced, the advantage of using the new
estimator becomes more extreme and convergence to the best estimate as the
number of sample $n$ increases is slower in all cases.

The histogram approach was applied with a huge number of data points; this
indicates that the histogram approach and the estimator described in this
paper approach similar values, but vastly more data is required for the
histogram approach. This discrepancy is more pronounced in the two-dimensional
example, and, in fact, becomes even more so as dimension increases further.

###  III.2 XY Model

In order to assess the ability of the new estimator to calculate transfer
entropy, a model is needed which combines noisy dynamics and causal
relationships. Here, the XY model is used.

Figure 5: XY model example. An example state of the XY model as described
here. The fixed pixel is in the centre, at $(25,25)$, with colour representing
the value of $\theta$ at each site.

The XY model is a well studied model in statistical mechanics which gives a
mathematical description of interacting spins. It is a lattice model in which
individual particles at lattice points have a spin described by a point on a
circle. The spins of adjacent particles interact in a way that models both the
tendency of spins to align and thermodynamics effects which allow for
fluctuations away from alignment. Here, a two-dimensional square-lattice XY
model is considered, with periodic boundary conditions in both directions,
giving a torus. Using

$\textbf{s}_{x}=(\cos\theta_{x},\ \sin\theta_{x})$ (28)

to denote the spin at the lattice point $x$, the energy is given

$\displaystyle H(\sigma)$
$\displaystyle=-J\sum_{x\in\mathcal{G}}\sum_{y\in\mathcal{N}(x)}\mathbf{s}_{x}\cdot\mathbf{s}_{y}$
$\displaystyle=-J\sum_{x\in\mathcal{G}}\sum_{y\in\mathcal{N}(x)}\cos(\theta_{x}-\theta_{y})$
(29)

where $\mathcal{N}(x)$ is the set made up of the four neighbours to $x$ and
$\mathcal{G}$ is the set of all points in the lattice. This corresponds the
simplest version of the model where coupling strengths are constant and there
is no external field. Here the ‘Metropolis’ version of the dynamics is used,
this is one of a number of different, ultimately equivalent, approaches to the
statistical dynamics. In this implementation, at each time step, all sites are
potentially updated. A change is proposed for a site by drawing a new value
from the von Mises distribution centered on the current value:

$\theta^{\prime}_{x}\sim\operatorname{VonMises}(\theta_{x},r)$ (30)

where $r$ is a concentration parameter; for small $r$ the distribution
approaches a uniform distribution, for large $r$ the value of
$\theta^{\prime}_{x}$ will tend to be very close to the value of $\theta_{x}$.
Now, the change in energy $\delta H$ that would result from switching
$\theta_{x}$ to the proposed new value $\theta^{\prime}_{x}$ is calculated. If
$\delta H<0$ the proposal is accepted corresponding to the dynamical tendency
of spins to align. If, however, $\delta H>0$ the proposal may still be
accepted with probability $\exp{(-\delta H/T)}$; this corresponds to thermal
fluctuations. $T$ is a temperature, an adjustable parameter which determines
the significance of the thermal effects. In practice, $r$, the value of the
concentration in the proposal distribution, is varied to keep the overall
acceptance rate equal to some target value $d$. In simulations, this is taken
to be $0.5$.

Here, the focus is on accessing causal structure. To this end, the XY-model
has been modified so there is one site on the grid that affects the dynamics
of the other sites, but is itself unaffected by its neighbours. In short, this
special site, which will be referred to as the causal site, changes according
to a Bernoulli random walk, it changes its angle by $\pm 0.2$ at each time
step. This is indepdent of the states of the surrounding sites, however the
causal point still ‘informs’ its neighbour in that it is incorporated in the
usual way in their update.

A similar story is shown in the case of transfer entropy. The histogramming
method becomes unworkable for larger values of $\ell$, the length of the past
considered. It can again be seen that the method here converges to the
eventual value of the histogramming method for far less data.

\begin{overpic}[width=99.73074pt,tics=10]{figs/different_pasts_1.pdf}
\put(-8.0,90.0){\huge{a}}
\end{overpic}\begin{overpic}[width=99.73074pt,tics=10]{figs/different_distances_1.pdf}
\put(-8.0,90.0){\huge{b}} \end{overpic}

Figure 6: Results for calculating transfer entropy from the XY model. In a the
relationship between number of samples used, and transfer entropy estimate is
shown. Figures show for past lengths of one and two respectively. Here,
distance is fixed at one. The dotted histogram line is after the histogram
estimate has converged, here using $5\times 10^{6}$ samples. b makes a similar
comparison, but varying distances of one and two respectively, with fixed past
length. Again, the dotted line represents the converged histogram method,
using a far greater, $5\times 10^{6}$, number of samples.

##  IV Discussion

This paper describes and demonstrates a model-free estimator for conditional
mutual information. As a metric-based approach it performs well for high-
dimensional data in a situation where other model-free approaches require huge
data sets. The estimator is very straight-forward, it relies on finding sets
of nearest points, in practice by sorting rows of a distance matrix, and
counting the sizes of the intersections of these sets.

This approach relies on a metric. It needs a ‘meaningful’ distance measure on
each of the outcome spaces for the random variable. If the metric is
unsuitable information may be lost: the estimator only estimates the
information in the distance matrix. In some applications appropriate distance
measures are well understood, this approach, for example, came from
neuroscience. In neuroscience there is an interest in the information content
of spike trains, the sequences of voltage pulses that propagate information
from neuron to neuron. The space of spike trains has good metrics but no
coordinates [21, 22]. For other data an appropriate metric may be less
obvious. Nonetheless, for data with coordinates there are always induced
metrics, such as the obvious $L^{2}$-metric and this proximity structure is
often implicitly exploited in estimation approaches, such as using a
histogram. In the estimator here, there is a broader choice of metric, since
the metric is not linked to the coordinates.

It has been assumed throughout that the smoothing parameter $h$ is the same
for each random variable so that each of the balls $B_{X}$, $B_{Y}$ and
$B_{Z}$ has the same volume. In fact, this is not required, separate values
$h_{X}$, $h_{Y}$ and $h_{Z}$ could be used. Maximising over all three could
only increase the estimated conditional mutual information; in cases where the
the information is very low, this could lead to an over estimate but this
freedom to have different smoothing parameters may be needed to get a good
estimate in examples where the three distributions are very different.

Another important factor in the application of this method is the
computational expense. The estimator appears to work effectively on a
comparably low amount of data but it does require the calculation of the
distance matrices for the observations in $\mathcal{X}$, $\mathcal{Y}$ and
$\mathcal{Z}$ and, significantly, it involves sorting each row of these
matrices. This quickly becomes the most time consuming step as the amount of
data increases so the computation time increases as $n\log{n}$.

### Applications

As well as transfer entropy, another measure based on conditional mutual
information is interaction information [4]. The interaction information
generalizes the mutual information to more than two variables, it has been
applied to problems in biology [23, 24, 25] and has interesting mathematical
properties with potential application to the study of intelligence [5]. For
three random variables it is

$I(X,Y,Z)=I(X,Y)-I(X,Y|Z)$ (31)

with the definition for higher numbers of variables given recursively:

$\displaystyle I(X_{1},X_{2},\ldots,X_{n})=$ $\displaystyle
I(X_{1},X_{2},\ldots,X_{n-1})$ (32)
$\displaystyle-I(X_{1},X_{2},\ldots,X_{n-1}|X_{n})$

Since the formula for $I(X,Y,Z)$ includes the conditional mutual information,
the estimator describe here gives an estimator for interaction information:

$I(X,Y,Z)\approx\frac{1}{n}\sum_{i}\log{\frac{nh_{XY}(i)h_{YZ}(i)h_{ZX}(i)}{h_{XYZ}(i)}}$
(33)

with similar, but increasingly complicated, formulae for interaction
information with more and more variables.

It is hoped that by reducing the amount of data required to calculate
conditional mutual information, and by implication, the transfer entropy, this
new estimator will unlock other applications of these quantities, including
applications in data science and machine learning.

##  Appendix A Dealing with draws when counting points

One issue that does arise is how to deal with draws: what happens if there are
multiple points equally far from the seed point and they cannot all be
included without the volume exceeding $h/n$? It might seem unlikely that there
will be draws given that the distance is generally a real number, however,
that does not account for the idiosyncrasy of data. In the example of spike
trains in neuroscience, some of the data points might correspond to empty
spike trains and these will all be the same distance from each other. As
another example, again related to neuroscience, one crude metric gives the
difference in the number of spikes as the distance between two spike trains,
this discrete metric will obviously lend itself to draws.

Draws can be dealt with by counting points fractionally. Consider drawing a
ball around $x_{i}$, $y_{i}$ or $z_{i}$ aiming to include $h$ points. Let $b$
be the number of points on the boundary: the draws or points that are
equidistant. Let $c$ be the points in the interior of the ball. The issue
arises when $b+c>h$. This is solved by including all of the interior $c$
points, as well as all of the boundary $b$ points. The boundary points however
are weighted at $(h-c)/b$, while interior points, as in the regular case, are
weighted one. These weights are taken in to account when counting the sizes of
intersections.

##  Appendix B Variance of estimates

Variance of the 50 conditional mutual information estimates in Figures 3, 4.

## References

  * Schreiber [2000] Thomas Schreiber.  Measuring information transfer.  _Physical Review Letters_ , 85(2):461, 2000\.
  * Palus̆ et al. [2001] Milan Palus̆, Vladimír Komárek, Zbynĕk Hrnc̆ír̆, and Katalin S̆tĕrbová.  Synchronization as adjustment of information rates: Detection from bivariate time series.  _Physical Review E_ , 63(4):046211, 2001.
  * Barnett et al. [2009] Lionel Barnett, Adam B Barrett, and Anil K Seth.  Granger causality and transfer entropy are equivalent for Gaussian variables.  _Physical Review Letters_ , 103(23):238701, 2009\.
  * McGill [1954] William McGill.  Multivariate information transmission.  _Transactions of the IRE Professional Group on Information Theory_ , 4(4):93–111, 1954.
  * Bell [2003] Anthony J Bell.  The co-information lattice.  In _Proceedings of the Fifth International Workshop on Independent Component Analysis and Blind Signal Separation: ICA_ , volume 2003\. Citeseer, 2003.
  * Tobin and Houghton [2013] R Joshua Tobin and Conor J Houghton.  A kernel-based calculation of information on a metric space.  _Entropy_ , 15(10):4540–4552, 2013.
  * Houghton [2015] Conor Houghton.  Calculating mutual information for spike trains and other data with distances but no coordinates.  _Royal Society Open Science_ , 2(5):140391, 2015\.
  * Kozachenko and Leonenko [1987] LF Kozachenko and Nikolai N Leonenko.  Sample estimate of the entropy of a random vector.  _Problemy Peredachi Informatsii_ , 23(2):9–16, 1987.
  * Kraskov et al. [2004] Alexander Kraskov, Harald Stögbauer, and Peter Grassberger.  Estimating mutual information.  _Physical Review E_ , 69(6):066138, 2004.
  * Frenzel and Pompe [2007] Stefan Frenzel and Bernd Pompe.  Partial mutual information for coupling analysis of multivariate time series.  _Physical review letters_ , 99:204101, Nov 2007.  doi: 10.1103/PhysRevLett.99.204101.
  * Lizier [2014] Joseph T Lizier.  Jidt: An information-theoretic toolkit for studying the dynamics of complex systems.  _Frontiers in Robotics and AI_ , 1:11, 2014.
  * Gómez-Herrero et al. [2015] Germán Gómez-Herrero, Wei Wu, Kalle Rutanen, Miguel C Soriano, Gordon Pipa, and Raul Vicente.  Assessing coupling dynamics from an ensemble of time series.  _Entropy_ , 17(4):1958–1970, 2015.
  * Houghton [2019] Conor Houghton.  Calculating the mutual information between two spike trains.  _Neural Computation_ , 31(2):330–343, 2019\.
  * Witter and Houghton [2021] Jake Witter and Conor Houghton.  A note on the unbiased estimation of mutual information.  _arXiv 2105.08682_ , 2021.
  * Press et al. [2007] William H Press, Saul A Teukolsky, William T Vetterling, and Brian P Flannery.  _Numerical recipes 3rd edition: The art of scientific computing_.  Cambridge university press, 2007.
  * Abramowitz and Stegun [1968] Milton Abramowitz and Irene A Stegun.  _Handbook of mathematical functions with formulas, graphs, and mathematical tables_ , volume 55.  US Government printing office, 1968.
  * Hlavác̆ková-Schindler et al. [2007] Katerina Hlavác̆ková-Schindler, Milan Palus̆, Martin Vejmelka, and Joydeep Bhattacharya.  Causality detection based on information-theoretic approaches in time series analysis.  _Physics Reports_ , 441(1):1–46, 2007.
  * Vicente et al. [2011] Raul Vicente, Michael Wibral, Michael Lindner, and Gordon Pipa.  Transfer entropy—a model-free measure of effective connectivity for the neurosciences.  _Journal of Computational Neuroscience_ , 30(1):45–67, 2011.
  * Staniek and Lehnertz [2008] Matthäus Staniek and Klaus Lehnertz.  Symbolic transfer entropy.  _Physical Review Letters_ , 100(15):158101, 2008\.
  * Bossomaier et al. [2016] Terry Bossomaier, Lionel Barnett, Michael Harré, and Joseph T Lizier.  _An introduction to transfer entropy_.  Springer, 2016.
  * Victor and Purpura [1997] Jonathan D Victor and Keith P Purpura.  Metric-space analysis of spike trains: theory, algorithms and application.  _Network: Computation in Neural Systems_ , 8(2):127–164, 1997.
  * Rossum [2001] MCW van Rossum.  A novel spike distance.  _Neural Computation_ , 13(4):751–763, 2001\.
  * Moore et al. [2006] Jason H Moore, Joshua C Gilbert, Chia-Ti Tsai, Fu-Tien Chiang, Todd Holden, Nate Barney, and Bill C White.  A flexible computational framework for detecting, characterizing, and interpreting statistical patterns of epistasis in genetic studies of human disease susceptibility.  _Journal of Theoretical Biology_ , 241(2):252–261, 2006.
  * Margolin et al. [2010] Adam A Margolin, Kai Wang, Andrea Califano, and Ilya Nemenman.  Multivariate dependence and genetic networks inference.  _IET Systems Biology_ , 4(6):428–440, 2010\.
  * LeVine and Weinstein [2014] Michael V LeVine and Harel Weinstein.  Nbit-a new information theory-based analysis of allosteric mechanisms reveals residues that underlie function in the leucine transporter leut.  _PLoS Computational Biology_ , 10(5):e1003603, 2014.