1  #region License Information


2  /* HeuristicLab


3  * Copyright (C) 20022015 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 


22  using System;


23  using System.Collections.Generic;


24  using System.Linq;


25 


26  namespace 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.rproject.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.rproject.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, sidestep 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  }

