Source code for monte_carlo_analysis.uncertainty_metrics.EarthMoverDistanceDistributionSimilarity

"""
**Author** : Robin Camarasa

**Institution** : Erasmus Medical Center

**Position** : PhD student

**Contact** : r.camarasa@erasmusmc.nl

**Date** : 2020-10-13

**Project** : monte_carlo_analysis

**Implement class EarthMoverDistanceDistributionSimilarity**

"""
from monte_carlo_analysis.utils.numba_utils import numba_cumul_histogram
from monte_carlo_analysis.uncertainty_metrics import DistributionSimilarityUncertaintyMetric
import numpy as np
from numba import jit


[docs]class EarthMoverDistanceDistributionSimilarity(DistributionSimilarityUncertaintyMetric):
    """Implement EarthMoverDistanceDistributionSimilarity. The formula applied to the pair of distributions :math:`q_{c'}(y_j|x)`
    and :math:`q_{c'}(y_j|x)` is:

    .. math::
        S^E(q_{c'}(y_j|x), q_{c''}(y_j|x)) = \mathcal{L}_1(\int_0^t q_{c'}(y_j = t | x),\int_0^t q_{c''}(y_j = t|x))

    :param nbins: The discretization step of the integral of the :math:`\mathcal{L}_1` norm
    """
    def __init__(self, nbins=100):
        super(EarthMoverDistanceDistributionSimilarity, self)
        self.nbins = nbins
        self.transformation = self.get_transformation(nbins)

[docs]    def get_transformation(self, nbins: int) -> callable:
        """Define the transformation applied to a pair of distributions

        :param nbins: The discretization step of the integral
        :return: transformation to apply to a pair of distributions
        """
        @jit(nopython=True)
        def transformation(distribution_1: np.array, distribution_2) -> float:
            # In 1D the EarthMoverdistance is proportional to the L1 distance
            # between 2 histogram
            cum_discretized_distribution_1 = numba_cumul_histogram(
                distribution_1, bins=nbins, min_value=0, max_value=1,
                normalized=True
            )[0]

            cum_discretized_distribution_2 = numba_cumul_histogram(
                distribution_2, bins=nbins, min_value=0, max_value=1,
                normalized=True
            )[0]

            return np.sum(np.abs(
                cum_discretized_distribution_1 - cum_discretized_distribution_2
            )) / nbins
        return transformation