source: trunk/sources/HeuristicLab.Problems.DataAnalysis/3.4/OnlineCalculators/DependencyCalculator/HoeffdingsDependenceCalculator.cs @ 12012

Last change on this file since 12012 was 12012, checked in by ascheibe, 6 years ago

#2212 merged r12008, r12009, r12010 back into trunk

File size: 5.3 KB
Line 
1#region License Information
2/* HeuristicLab
3 * Copyright (C) 2002-2015 Heuristic and Evolutionary Algorithms Laboratory (HEAL)
4 *
5 * This file is part of HeuristicLab.
6 *
7 * HeuristicLab is free software: you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation, either version 3 of the License, or
10 * (at your option) any later version.
11 *
12 * HeuristicLab is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with HeuristicLab. If not, see <http://www.gnu.org/licenses/>.
19 */
20#endregion
21
22using System;
23using System.Collections.Generic;
24using System.Linq;
25
26namespace HeuristicLab.Problems.DataAnalysis {
27  public class HoeffdingsDependenceCalculator : IDependencyCalculator {
28
29    public double Maximum { get { return 1.0; } }
30
31    public double Minimum { get { return -0.5; } }
32
33    public string Name { get { return "Hoeffdings Dependence"; } }
34
35    public double Calculate(IEnumerable<double> originalValues, IEnumerable<double> estimatedValues, out OnlineCalculatorError errorState) {
36      return HoeffdingsDependenceCalculator.CalculateHoeffdings(originalValues, estimatedValues, out errorState);
37    }
38
39    public static double CalculateHoeffdings(IEnumerable<double> originalValues, IEnumerable<double> estimatedValues, out OnlineCalculatorError errorState) {
40      double d = HoeffD(originalValues, estimatedValues, out errorState);
41      if (errorState != OnlineCalculatorError.None) return double.NaN;
42      return d;
43    }
44
45    /// <summary>
46    /// computes Hoeffding's dependence coefficient.
47    /// Source: hoeffd.r from R package hmisc http://cran.r-project.org/web/packages/Hmisc/index.html
48    /// </summary>
49    private static double HoeffD(IEnumerable<double> xs, IEnumerable<double> ys, out OnlineCalculatorError errorState) {
50      double[] rx = TiedRank(xs);
51      double[] ry = TiedRank(ys);
52      if (rx.Length != ry.Length) throw new ArgumentException("The number of elements in xs and ys does not match");
53      double[] rxy = TiedRank(xs, ys);
54
55      int n = rx.Length;
56      double q = 0, r = 0, s = 0;
57      double scaling = 1.0 / (n * (n - 1));
58      for (int i = 0; i < n; i++) {
59        q += (rx[i] - 1) * (rx[i] - 2) * (ry[i] - 1) * (ry[i] - 2) * scaling;
60        r += (rx[i] - 2) * (ry[i] - 2) * rxy[i] * scaling;
61        s += rxy[i] * (rxy[i] - 1) * scaling;
62      }
63      errorState = OnlineCalculatorError.None;
64      // return 30.0 * (q - 2 * (n - 2) * r + (n - 2) * (n - 3) * s) / n / (n - 1) / (n - 2) / (n - 3) / (n - 4);
65      double t0 = q / (n - 2) / (n - 3) / (n - 4);
66      double t1 = 2 * r / (n - 3) / (n - 4);
67      double t2 = s / (n - 4);
68      return 30.0 * (t0 - t1 + t2);
69    }
70
71    private static double[] TiedRank(IEnumerable<double> xs) {
72      var xsArr = xs.ToArray();
73      var idx = Enumerable.Range(1, xsArr.Length).ToArray();
74      Array.Sort(xsArr, idx);
75      CRank(xsArr);
76      Array.Sort(idx, xsArr);
77      return xsArr;
78    }
79
80    /// <summary>
81    /// Calculates the joint rank with midranks for ties. Source: hoeffd.r from R package hmisc http://cran.r-project.org/web/packages/Hmisc/index.html
82    /// </summary>
83    /// <param name="xs"></param>
84    /// <param name="ys"></param>
85    /// <returns></returns>
86    private static double[] TiedRank(IEnumerable<double> xs, IEnumerable<double> ys) {
87      var xsArr = xs.ToArray();
88      var ysArr = ys.ToArray();
89      var r = new double[xsArr.Length];
90      int n = r.Length;
91      for (int i = 0; i < n; i++) {
92        var xi = xsArr[i];
93        var yi = ysArr[i];
94        double ri = 0.0;
95        for (int j = 0; j < n; j++) {
96          if (i != j) {
97            double cx;
98            if (xsArr[j] < xi) cx = 1.0;
99            else if (xsArr[j] > xi) cx = 0.0;
100            else cx = 0.5;  // eq
101            double cy;
102            if (ysArr[j] < yi) cy = 1.0;
103            else if (ysArr[j] > yi) cy = 0.0;
104            else cy = 0.5; // eq
105            ri = ri + cx * cy;
106          }
107        }
108        r[i] = ri;
109      }
110      return r;
111    }
112
113    /// <summary>
114    /// Calculates midranks. Source: Numerical Recipes in C. p 642
115    /// </summary>
116    /// <param name="w">Sorted array of elements, replaces the elements by their rank, including midranking of ties</param>
117    /// <returns></returns>
118    private static void CRank(double[] w) {
119      int i = 0;
120      int n = w.Length;
121      while (i < n - 1) {
122        if (w[i + 1] > w[i]) {    // w[i+1] must be larger or equal w[i] as w must be sorted
123          // not a tie
124          w[i] = i + 1;
125          i++;
126        } else {
127          int j;
128          for (j = i + 1; j < n && w[j] <= w[i]; j++) ; // how far does it go (<= effectively means == as w must be sorted, side-step equality for double values)
129          double rank = 1 + 0.5 * (i + j - 1);
130          int k;
131          for (k = i; k < j; k++) w[k] = rank; // set the rank for all tied entries
132          i = j;
133        }
134      }
135
136      if (i == n - 1) w[n - 1] = n;   // if the last element was not tied, this is its rank
137    }
138  }
139}
Note: See TracBrowser for help on using the repository browser.