"""
The Kullback-Leibler divergence.
"""
import numpy as np
from ..helpers import as_linear, normalize_pmfs, normalize_rvs
from ..utils import flatten, unitful
from .kullback_leibler_divergence import kullback_leibler_divergence
__all__ = (
"double_power_sum",
"hellinger_sum",
"alpha_divergence",
"hellinger_divergence",
"renyi_divergence",
"tsallis_divergence",
)
### References for Divergence Formulas ###
## http://arxiv.org/pdf/1105.3259v1.pdf
## http://arxiv.org/pdf/1206.2459.pdf
## http://mitran-lab.amath.unc.edu:8082/subversion/grants/Proposals/2013/DOE-DataCentric/biblio/LieseVajdaDivergencesInforTheory.pdf
## Crooks: http://threeplusone.com/on_information.pdf
def double_power_sum(dist1, dist2, exp1=1, exp2=1, rvs=None):
"""
A common generalization of the sums needed to compute the Hellinger
and alpha divergences below.
Parameters
----------
dist1 : Distribution
The first distribution in the sum.
The second distribution in the sum.
exp1 : float, 1
First exponent used in the power sum.
exp2 : float, 1
Second exponent used in the power sum.
rvs : list, None
The indexes of the random variable used to calculate the sum.
If None, then the sum is calculated over all random variables.
Returns
-------
dkl : float
The specified sum between `dist1` and `dist2`.
Raises
------
ditException
Raised if either `dist1` or `dist2` doesn't have `rvs` or, if `rvs` is
None, if `dist2` has an outcome length different than `dist1`.
"""
ps, qs = _union_pmfs(dist1, dist2, rvs)
div = np.sum(_power(ps, exp1) * _power(qs, exp2))
return div
def _kl_nats(dist1, dist2, rvs):
"""KL divergence in nats, independent of the distributions' bases."""
kl_bits = kullback_leibler_divergence(as_linear(dist1), as_linear(dist2), rvs=rvs)
return getattr(kl_bits, "magnitude", kl_bits) * np.log(2)
def _union_pmfs(dist1, dist2, rvs):
"""
Linear pmfs of `dist1` and `dist2` over `rvs`, aligned on the union of
their supports.
"""
rvs = _validate_rvs(dist1, dist2, rvs)
return normalize_pmfs(dist1.marginal(rvs), dist2.marginal(rvs))
def _validate_rvs(dist1, dist2, rvs):
"""
Normalize `rvs` against `dist1` and check that `dist2` has them too.
"""
rvs, _ = normalize_rvs(dist1, rvs, None)
rvs = list(flatten(rvs))
normalize_rvs(dist2, rvs, None)
return rvs
def _power(x, e):
"""
``x**e`` with ``0**e = inf`` for ``e < 0`` and ``0**0 = 0``, so that an
outcome outside a distribution's support contributes nothing at ``e = 0``.
"""
with np.errstate(divide="ignore"):
return np.where(x > 0, np.power(x, e), np.inf if e < 0 else 0.0)
def hellinger_sum(dist1, dist2, alpha=1.0, rvs=None):
"""
The Hellinger sum/integral of `dist1` and `dist2`, used to define other
divergences.
Parameters
----------
dist1 : Distribution
The first distribution in the sum.
dist2 : Distribution
The second distribution in the sum.
alpha : float, 1
Exponent used in the sum.
rvs : list, None
The indexes of the random variable used to calculate the sum.
If None, then the sum is calculated over all random variables.
Returns
-------
dkl : float
The Hellinger sum between `dist1` and `dist2`.
Raises
------
ditException
Raised if either `dist1` or `dist2` doesn't have `rvs` or, if `rvs` is
None, if `dist2` has an outcome length different than `dist1`.
"""
return double_power_sum(dist1, dist2, alpha, 1 - alpha, rvs=rvs)
[docs]
def hellinger_divergence(dist1, dist2, alpha=1.0, rvs=None):
# http://mitran-lab.amath.unc.edu:8082/subversion/grants/Proposals/2013/DOE-DataCentric/biblio/LieseVajdaDivergencesInforTheory.pdf
"""
The Hellinger divergence of `dist1` and `dist2`.
Parameters
----------
dist1 : Distribution
The first distribution in the Hellinger divergence.
dist2 : Distribution
The second distribution in the Hellinger divergence.
alpha : float, 1
The divergence is a one parameter family in alpha.
rvs : list, None
The indexes of the random variable used to calculate the
Hellinger divergence between. If None, then the Hellinger
divergence is calculated over all random variables.
Returns
-------
dkl : float
The Hellinger divergence between `dist1` and `dist2`.
Raises
------
ditException
Raised if either `dist1` or `dist2` doesn't have `rvs` or, if `rvs` is
None, if `dist2` has an outcome length different than `dist1`.
"""
if alpha == 1:
return _kl_nats(dist1, dist2, rvs)
s = hellinger_sum(dist1, dist2, rvs=rvs, alpha=alpha)
return (s - 1) / (alpha - 1)
[docs]
def tsallis_divergence(dist1, dist2, alpha=1.0, rvs=None):
"""
The Tsallis divergence of `dist1` and `dist2`.
Parameters
----------
dist1 : Distribution
The first distribution in the Tsallis divergence.
dist2 : Distribution
The second distribution in the Tsallis divergence.
alpha : float, 1
The divergence is a one parameter family in alpha.
rvs : list, None
The indexes of the random variable used to calculate the
Tsallis divergence between. If None, then the Tsallis
divergence is calculated over all random variables.
Returns
-------
dkl : float
The Tsallis divergence between `dist1` and `dist2`.
Raises
------
ditException
Raised if either `dist1` or `dist2` doesn't have `rvs` or, if `rvs` is
None, if `dist2` has an outcome length different than `dist1`.
"""
# D_T = (D_alpha -1) / (alpha-1)
if alpha == 1:
div = _kl_nats(dist1, dist2, rvs)
else:
s = hellinger_sum(dist1, dist2, rvs=rvs, alpha=alpha)
div = (s - 1) / (alpha - 1)
return div
[docs]
@unitful
def renyi_divergence(dist1, dist2, alpha=1, rvs=None):
"""
The Renyi divergence of `dist1` and `dist2`.
Parameters
----------
dist1 : Distribution
The first distribution in the Renyi divergence.
dist2 : Distribution
The second distribution in the Renyi divergence.
alpha : float, 1
The divergence is a one parameter family in alpha.
rvs : list, None
The indexes of the random variable used to calculate the
Renyi divergence between. If None, then the Renyi
divergence is calculated over all random variables.
Returns
-------
dkl : float
The Renyi divergence between `dist1` and `dist2`.
Raises
------
ditException
Raised if either `dist1` or `dist2` doesn't have `rvs` or, if `rvs` is
None, if `dist2` has an outcome length different than `dist1`.
"""
# D_R = log D_alpha / (alpha-1)
if alpha == 1:
div = kullback_leibler_divergence(dist1, dist2, rvs=rvs)
else:
s = hellinger_sum(dist1, dist2, rvs=rvs, alpha=alpha)
div = np.log2(s) / (alpha - 1.0)
return div
[docs]
def alpha_divergence(dist1, dist2, alpha=1.0, rvs=None):
"""
The alpha divergence of `dist1` and `dist2`, as used in Information
Geometry. Note there is more than one inequivalent definition of
"alpha divergence" in the literature, this one comes
from http://en.wikipedia.org/wiki/Information_geometry .
Parameters
----------
dist1 : Distribution
The first distribution in the alpha divergence.
dist2 : Distribution
The second distribution in the alpha divergence.
alpha : float, 1
The divergence is a one parameter family in alpha.
rvs : list, None
The indexes of the random variable used to calculate the
alpha divergence between. If None, then the alpha
divergence is calculated over all random variables.
Returns
-------
dkl : float
The alpha divergence between `dist1` and `dist2`.
Raises
------
ditException
Raised if either `dist1` or `dist2` doesn't have `rvs` or, if `rvs` is
None, if `dist2` has an outcome length different than `dist1`.
"""
if alpha == 1:
rvs = _validate_rvs(dist1, dist2, rvs)
return _kl_nats(dist2, dist1, rvs)
if alpha == -1:
return _kl_nats(dist1, dist2, rvs)
s = double_power_sum(dist1, dist2, (1 - alpha) / 2, (1 + alpha) / 2, rvs=rvs)
return 4 * (1 - s) / (1 - alpha * alpha)
[docs]
def f_divergence(dist1, dist2, f, rvs=None):
"""
The Csiszar f-divergence of `dist1` and `dist2`. Note that it is typically
more accurate to use a specialized divergence function when available
due to roundoff errors and small probability effects.
Parameters
----------
dist1 : Distribution
The first distribution in the f-divergence.
dist2 : Distribution
The second distribution in the f-divergence.
f : function
The auxillary function defining the f-divergence
rvs : list, None
The indexes of the random variable used to calculate the
f-divergence between. If None, then the
f-divergence is calculated over all random variables.
Returns
-------
dkl : float
The f-divergence between `dist1` and `dist2`.
Raises
------
ditException
Raised if either `dist1` or `dist2` doesn't have `rvs` or, if `rvs` is
None, if `dist2` has an outcome length different than `dist1`.
Notes
-----
Outcomes where `dist2` is zero but `dist1` is not are omitted: their
contribution, ``p * lim_{t -> inf} f(t) / t``, depends on `f` and is not
computed.
"""
ps, qs = _union_pmfs(dist1, dist2, rvs)
vfunc = np.vectorize(f)
with np.errstate(divide="ignore", invalid="ignore"):
div = np.nansum(qs * vfunc(np.divide(ps, qs)))
return div