Commits

Matija Polajnar committed f50b356

Merge in Lan Umek's implementations of reliability estimation for classification.

  • Participants
  • Parent commits 7dfb468

Comments (0)

Files changed (2)

File _reliability/__init__.py

 BLENDING_ABSOLUTE = 9
 ICV_METHOD = 10
 MAHAL_TO_CENTER_ABSOLUTE = 13
+DENS_ABSOLUTE = 14
 
 # Type of estimator constant
 SIGNED = 0
                3: "BAGV absolute", 4: "CNK signed", 5: "CNK absolute",
                6: "LCV absolute", 7: "BVCK_absolute", 8: "Mahalanobis absolute",
                9: "BLENDING absolute", 10: "ICV", 11: "RF Variance", 12: "RF Std",
-               13: "Mahalanobis to center"}
+               13: "Mahalanobis to center", 14: "Density based"}
 
 select_with_repeat = Orange.core.MakeRandomIndicesMultiple()
 select_with_repeat.random_generator = Orange.misc.Random()
     t = r * (df / ((-r + 1.0 + 1e-30) * (r + 1.0 + 1e-30))) ** 0.5
     return statc.betai (df * 0.5, 0.5, df / (df + t * t))
 
+
+# Distances between two discrete probability distributions
+#TODO Document those.
+def normalize_both(p, q):
+    if not p.normalized:
+        p.normalize()
+    if not q.normalized:
+        q.normalize()
+    return p, q
+
+def minkowsky_dist(p, q, m=2):
+    p, q = normalize_both(p, q)
+    dist = 0
+    for i in range(len(p)):
+        dist += abs(p[i]-q[i])**m
+    return dist**(1./m)
+
+def manhattan_distance(p, q):
+    return minkowsky_dist(p, q, m=1)
+
+def euclidean_dist(p, q):
+    return minkowsky_dist(p, q, m=2)
+
+def variance_dist(p, q):
+    return euclidean_dist(p, q) ** 2
+
+def max_dist(p, q):
+    p, q = normalize_both(p, q)
+    return max([abs(p[i]-q[i]) for i in range(len(p))])
+
+def hellinger_dist(p, q):
+    p, q = normalize_both(p, q)
+    dist = 0
+    for i in range(len(p)):
+        dist += (math.sqrt(p[i])-math.sqrt(q[i])) ** 2
+    return dist
+
+def my_log(x):
+    return 0 if x == 0 else x * math.log(x)
+
+def kullback_leibler(p, q):
+    p, q = normalize_both(p, q)
+    dist = 0
+    for i in range(len(p)):
+        dist += my_log(p[i]-q[i])
+    return dist
+
+def cosine(p, q):
+    p, q = normalize_both(p, q)
+    p, q = [pp for pp in p], [qq for qq in q]
+    return 1 - numpy.dot(x,y) / (numpy.linalg.norm(p)*numpy.linalg.norm(q))
+
+
 class Estimate:
     """
     Reliability estimate. Contains attributes that describe the results of
     :rtype: :class:`Orange.evaluation.reliability.BaggingVarianceClassifier`
     
     :math:`m` different bagging models are constructed and used to estimate
-    the value of dependent variable for a given instance. The variance of
-    those predictions is used as a prediction reliability estimate.
+    the value of dependent variable for a given instance. In regression,
+    the variance of those predictions is used as a prediction reliability
+    estimate.
 
     :math:`BAGV = \\frac{1}{m} \sum_{i=1}^{m} (K_i - K)^2`
 
     where :math:`K = \\frac{\sum_{i=1}^{m} K_i}{m}` and :math:`K_i` are
-    predictions of individual constructed models.
+    predictions of individual constructed models. Note that a greater value
+    implies greater error.
+
+    For classification, 1 minus the average Euclidean distance between class
+    probability distributions predicted by the model, and distributions
+    predicted by the individual bagged models, is used as the BAGV reliability
+    measure. Note that in this case a greater value implies a better
+    prediction.
     
     """
     def __init__(self, m=50):
     def __call__(self, instances, learner):
         classifiers = []
 
+        if instances.domain.class_var.var_type == Orange.feature.Descriptor.Discrete:
+            classifier = learner(instances)
+        else:
+            classifier = None
+
         # Create bagged classifiers using sampling with replacement
         for _ in xrange(self.m):
             selection = select_with_repeat(len(instances))
             data = instances.select(selection)
             classifiers.append(learner(data))
-        return BaggingVarianceClassifier(classifiers)
+        return BaggingVarianceClassifier(classifiers, classifier)
 
 class BaggingVarianceClassifier:
     def __init__(self, classifiers):
         self.classifiers = classifiers
 
     def __call__(self, instance, *args):
+        def __init__(self, classifiers, classifier=None):
+            self.classifiers = classifiers
+            self.classifier = classifier
+
+    def __call__(self, instance, *args):
         BAGV = 0
 
         # Calculate the bagging variance
-        bagged_values = [c(instance, Orange.core.GetValue).value for c in self.classifiers if c is not None]
-
+        if instance.domain.class_var.var_type == Orange.feature.Descriptor.Continuous:
+            bagged_values = [c(instance, Orange.core.GetValue).value for c in self.classifiers if c is not None]
+        elif instance.domain.class_var.var_type == Orange.feature.Descriptor.Discrete:
+            estimate = self.classifier(instance, Orange.core.GetProbabilities)
+            bagged_values = [euclidean_dist(c(instance, Orange.core.GetProbabilities), estimate) for c in self.classifiers if c is not None]
         k = sum(bagged_values) / len(bagged_values)
 
         BAGV = sum((bagged_value - k) ** 2 for bagged_value in bagged_values) / len(bagged_values)
+        if instance.domain.class_var.var_type == Orange.feature.Descriptor.Discrete:
+            BAGV = 1 - BAGV
 
         return [Estimate(BAGV, ABSOLUTE, BAGV_ABSOLUTE)]
 
 class LocalCrossValidation:
     """
-    
+
     :param k: Number of nearest neighbours used in LCV estimate
     :type k: int
-    
+
+    :param distance: function that computes a distance between two discrete
+        distributions (used only in classification problems). The default
+        is Hellinger distance.
+    :type distance: function
+
+    :param distance_weighted: for classification reliability estimation,
+        use an average distance between distributions, weighted by :math:`e^{-d}`,
+        where :math:`d` is the distance between predicted instance and the
+        neighbour.
+
     :rtype: :class:`Orange.evaluation.reliability.LocalCrossValidationClassifier`
-    
+
     :math:`k` nearest neighbours to the given instance are found and put in
     a separate data set. On this data set, a leave-one-out validation is
-    performed. Reliability estimate is then the distance weighted absolute
-    prediction error.
+    performed. Reliability estimate for regression is then the distance
+    weighted absolute prediction error. In classification, 1 minus the average
+    distance between the predicted class probability distribution and the
+    (trivial) probability distributions of the nearest neighbour.
 
     If a special value 0 is passed as :math:`k` (as is by default),
     it is set as 1/20 of data set size (or 5, whichever is greater).
-    
+
+    Summary of the algorithm for regression:
+
     1. Determine the set of k nearest neighours :math:`N = { (x_1, c_1),...,
        (x_k, c_k)}`.
     2. On this set, compute leave-one-out predictions :math:`K_i` and
        prediction errors :math:`E_i = | C_i - K_i |`.
     3. :math:`LCV(x) = \\frac{ \sum_{(x_i, c_i) \in N} d(x_i, x) * E_i }{ \sum_{(x_i, c_i) \in N} d(x_i, x) }`
-    
+
     """
-    def __init__(self, k=0):
+    def __init__(self, k=0, distance=hellinger_dist, distance_weighted=True):
         self.k = k
+        self.distance = distance
+        self.distance_weighted = distance_weighted
 
     def __call__(self, instances, learner):
         nearest_neighbours_constructor = Orange.classification.knn.FindNearestConstructor()
         if self.k == 0:
             self.k = max(5, len(instances) / 20)
 
-        return LocalCrossValidationClassifier(distance_id, nearest_neighbours, self.k, learner)
+        return LocalCrossValidationClassifier(distance_id, nearest_neighbours, self.k, learner,
+            distance=self.distance, distance_weighted=self.distance_weighted)
 
 class LocalCrossValidationClassifier:
-    def __init__(self, distance_id, nearest_neighbours, k, learner):
+    def __init__(self, distance_id, nearest_neighbours, k, learner, **kwds):
         self.distance_id = distance_id
         self.nearest_neighbours = nearest_neighbours
         self.k = k
         self.learner = learner
+        for a,b in kwds.items():
+            setattr(self, a, b)
 
     def __call__(self, instance, *args):
         LCVer = 0
 
             classifier = self.learner(Orange.data.Table(train))
 
-            returned_value = classifier(knn[i], Orange.core.GetValue)
+            if instance.domain.class_var.var_type == Orange.feature.Descriptor.Continuous:
+                returned_value = classifier(knn[i], Orange.core.GetValue)
+                e = abs(knn[i].getclass().value - returned_value.value)
 
-            e = abs(knn[i].getclass().value - returned_value.value)
+            elif instance.domain.class_var.var_type == Orange.feature.Descriptor.Discrete:
+                returned_value = classifier(knn[i], Orange.core.GetProbabilities)
+                probabilities = [knn[i].get_class() == val for val in instance.domain.class_var.values]
+                e = self.distance(returned_value, Orange.statistics.distribution.Discrete(probabilities))
 
-            LCVer += e * math.exp(-knn[i][self.distance_id])
-            LCVdi += math.exp(-knn[i][self.distance_id])
+            dist = math.exp(-knn[i][self.distance_id]) if self.distance_weighted else 1.0
+            LCVer += e * dist
+            LCVdi += dist
 
         LCV = LCVer / LCVdi if LCVdi != 0 else 0
         if math.isnan(LCV):
             LCV = 0.0
+
+        if instance.domain.class_var.var_type == Orange.feature.Descriptor.Discrete:
+            LCV = 1 - LCV
+
         return [ Estimate(LCV, ABSOLUTE, LCV_ABSOLUTE) ]
 
 class CNeighbours:
     
     :param k: Number of nearest neighbours used in CNK estimate
     :type k: int
+
+    :param distance: function that computes a distance between two discrete
+        distributions (used only in classification problems). The default
+        is Hellinger distance.
+    :type distance: function
     
     :rtype: :class:`Orange.evaluation.reliability.CNeighboursClassifier`
     
-    CNK is defined for an unlabeled instance as a difference between average
-    label of its nearest neighbours and its prediction. CNK can be used as a
-    signed or absolute estimate.
+    For regression, CNK is defined for an unlabeled instance as a difference
+    between average label of its nearest neighbours and its prediction. CNK
+    can be used as a signed or absolute estimate.
     
     :math:`CNK = \\frac{\sum_{i=1}^{k}C_i}{k} - K`
     
     where :math:`k` denotes number of neighbors, C :sub:`i` denotes neighbours'
-    labels and :math:`K` denotes the instance's prediction.
+    labels and :math:`K` denotes the instance's prediction. Note that a greater
+    value implies greater prediction error.
+
+    For classification, CNK is equal to 1 minus the average distance between
+    predicted class distribution and (trivial) class distributions of the
+    $k$ nearest neighbours from the learning set. Note that in this case
+    a greater value implies better prediction.
     
     """
-    def __init__(self, k=5):
+    def __init__(self, k=5, distance=hellinger_dist):
         self.k = k
+        self.distance = distance
 
     def __call__(self, instances, learner):
         nearest_neighbours_constructor = Orange.classification.knn.FindNearestConstructor()
 
         distance_id = Orange.feature.Descriptor.new_meta_id()
         nearest_neighbours = nearest_neighbours_constructor(instances, 0, distance_id)
-        return CNeighboursClassifier(nearest_neighbours, self.k)
+        return CNeighboursClassifier(nearest_neighbours, self.k, distance=self.distance)
 
 class CNeighboursClassifier:
     def __init__(self, nearest_neighbours, k):
         knn = [ex for ex in self.nearest_neighbours(instance, self.k)]
 
         # average label of neighbors
-        for ex in knn:
-            CNK += ex.getclass().value
+        if ex.domain.class_var.var_type == Orange.feature.Descriptor.Continuous:
+            for ex in knn:
+                CNK += ex.getclass().value
+            CNK /= self.k
+            CNK -= predicted.value
 
-        CNK /= self.k
-        CNK -= predicted.value
+            return [Estimate(CNK, SIGNED, CNK_SIGNED),
+                    Estimate(abs(CNK), ABSOLUTE, CNK_ABSOLUTE)]
+        elif ex.domain.class_var.var_type == Orange.feature.Descriptor.Discrete:
+            knn_l = Orange.classification.knn.kNNLearner(k=self.k)
+            knn_c = knn_l(knn)
+            for ex in knn:
+                CNK -= self.distance(probabilities, knn_c(ex, Orange.classification.Classifier.GetProbabilities))
+            CNK /= self.k
+            CNK += 1
 
-        return [Estimate(CNK, SIGNED, CNK_SIGNED),
-                Estimate(abs(CNK), ABSOLUTE, CNK_ABSOLUTE)]
+            return [Estimate(CNK, ABSOLUTE, CNK_ABSOLUTE)]
 
 class Mahalanobis:
     """
 
         return [Estimate(value.value, SIGNED, SABIAS_SIGNED)]
 
+def gauss_kernel(x, sigma=1):
+    return 1./(sigma*math.sqrt(2*math.pi)) * math.exp(-1./2*(x/sigma)**2)
+
+class ParzenWindowDensityBased:
+    """
+    :param K: kernel function. Default: gaussian.
+    :type K: function
+
+    :param d_measure: distance measure for inter-instance distance.
+    :type d_measure: :class:`Orange.distance.DistanceConstructor`
+
+    :rtype: :class:`Orange.evaluation.reliability.ParzenWindowDensityBasedClassifier`
+
+    Returns a value that estimates a density of problem space around the
+    instance being predicted.
+    """
+    def __init__(self, K=gauss_kernel, d_measure=Orange.distance.Euclidean()):
+        self.K = K
+        self.d_measure = d_measure
+
+    def __call__(self, instances):
+
+        self.distance = self.d_measure(instances)
+
+        def density(x):
+            l, dens = len(instances), 0
+            for ex in instances:
+                dens += self.K(self.distance(x,ex))
+            return dens / l
+
+        max_density = max([density(ex) for ex in instances])
+
+        return ParzenWindowDensityBasedClassifier(density, max_density)
+
+class ParzenWindowDensityBasedClassifier:
+
+    def __init__(self, density, max_density):
+        self.density = density
+        self.max_density = max_density
+
+
+    def __call__(self, instance, *args):
+
+        DENS = self.max_density-self.density(instance)
+
+        return [Estimate(DENS, ABSOLUTE, DENS_ABSOLUTE)]
+
 class Learner:
     """
     Reliability estimation wrapper around a learner we want to test.
             return probabilities
         else:
             return predicted, probabilities
+
+# Functions for testing and plotting
+#TODO Document those.
+def get_acc_rel(method, data, learner):
+    estimators = [method]
+    reliability = Orange.evaluation.reliability.Learner(learner, estimators=estimators)
+    #results = Orange.evaluation.testing.leave_one_out([reliability], data)
+    results = Orange.evaluation.testing.cross_validation([reliability], data)
+
+    rels, acc = [], []
+
+    for res in results.results:
+        rels.append(res.probabilities[0].reliability_estimate[0].estimate)
+        acc.append(res.probabilities[0][res.actual_class])
+
+    return rels, acc
+
+def acc_rel_plot(method, data, learner, file_name="acc_rel_plot.png", colors=None):
+
+    import matplotlib.pylab as plt
+
+    plt.clf()
+
+    rels, acc = get_acc_rel(method, data, learner)
+    print "rels", rels
+    print "acc", acc
+
+    if colors is None:
+        colors = "k"
+    plt.scatter(acc, rels, c=colors)
+    plt.xlim(0.,1.)
+    plt.ylim(ymin=0.)
+    plt.savefig(file_name)
+
+def acc_rel_correlation(method, data, learner):
+    import scipy.stats
+    rels, acc = get_acc_rel(method, data, learner)
+    return scipy.stats.spearmanr(acc, rels)[0]

File docs/rst/Orange.evaluation.reliability.rst

 *************************************
 
 Reliability assessment statistically predicts reliability of single
-predictions. Most of implemented algorithms are taken from Comparison of
-approaches for estimating reliability of individual regression predictions,
-Zoran Bosnić, 2008.
+predictions. Most of implemented algorithms for regression are taken from
+Comparison of approaches for estimating reliability of individual
+regression predictions, Zoran Bosnić, 2008. Implementations for
+classification follow descriptions in Evaluating Reliability of Single
+Classifications of Neural Networks, Darko Pevec, 2011.
 
 The following example shows basic usage of reliability estimation methods:
 
 Reliability Methods
 ===================
 
+For regression, all the described measures can be used. Classification domains
+are supported by the following methods: BAGV, LCV, CNK and DENS.
+
 Sensitivity Analysis (SAvar and SAbias)
 ---------------------------------------
 .. autoclass:: SensitivityAnalysis
 
 .. autoclass:: MahalanobisToCenter
 
+Density estimation using Parzen window (DENS)
+---------------------------------------------
+
+.. autoclass:: ParzenWindowDensityBased
+
 Reliability estimation wrappers
 ===============================
 
 .org/abstract_S0269888909990154>`_ *The Knowledge Engineering Review* 25(1),
 pp. 27-47.
 
+Pevec, D., Štrumbelj, E., Kononenko, I. (2011) `Evaluating Reliability of
+Single Classifications of Neural Networks. <http://www.springerlink.com
+/content/48u881761h127r33/export-citation/>`_ *Adaptive and Natural Computing
+Algorithms*, 2011, pp. 22-30.