ktau.hpp

// Copyright © 2020 Thomas Nagler
//
// This file is part of the wdm library and licensed under the terms of
// the MIT license. For a copy, see the LICENSE file in the root directory
// or https://github.com/tnagler/wdm/blob/master/LICENSE.
 
#pragma once
 
#include "utils.hpp"
 
namespace wdm {
 
namespace impl {
 
inline void
normalize_weights(std::vector<double>& w)
{
  if (w.size() > 0) {
    double s = utils::sum(w);
    for (size_t i = 0; i < w.size(); i++)
      w[i] /= s;
  }
}
 
inline double
ktau(std::vector<double> x,
     std::vector<double> y,
     std::vector<double> weights = std::vector<double>())
{
  utils::check_sizes(x, y, weights);
 
  // 1.1 Sort x, y, and weights in x order; break ties in according to y.
  utils::sort_all(x, y, weights);
 
  // 1.2 Count pairs of tied x and simultaneous ties in x and y.
  double ties_x = utils::count_tied_pairs(x, weights);
  double ties_both = utils::count_joint_ties(x, y, weights);
 
  // 2.1 Sort y again and count exchanges (= number of discordant pairs).
  double num_d = 0.0;
  utils::merge_sort(y, weights, num_d);
 
  // 2.2 Count pairs of tied y.
  double ties_y = utils::count_tied_pairs(y, weights);
 
  // 3. Calculate Kendall's tau.
  if (weights.size() == 0)
    weights = std::vector<double>(x.size(), 1.0);
  double num_pairs = utils::perm_sum(weights, 2);
  double num_c = num_pairs - (num_d + ties_x + ties_y - ties_both);
  double tau = num_c - num_d;
  tau /= std::sqrt((num_pairs - ties_x) * (num_pairs - ties_y));
 
  return tau;
}
 
inline double
ktau_stat_adjust(std::vector<double> x,
                 std::vector<double> y,
                 std::vector<double> weights)
{
  utils::check_sizes(x, y, weights);
 
  if (weights.size() == 0)
    weights = std::vector<double>(x.size(), 1.0);
  // Put weights in effective-sample-size units: both their sum and squared
  // sum then equal n_eff, while their relative magnitudes remain unchanged.
  double effective_scale =
    utils::sum(weights) / utils::sum(utils::pow(weights, 2));
  for (auto& weight : weights)
    weight *= effective_scale;
 
  // 1.1 Sort x, y, and weights in x order; break ties in according to y.
  utils::sort_all(x, y, weights);
 
  // 1.2 Count pairs and triplets of tied x and simultaneous ties in x and y.
  double pair_x = utils::count_tied_pairs(x, weights);
  double trip_x = utils::count_tied_triplets(x, weights);
  double v_x = utils::count_ties_v(x, weights);
 
  // 2.1 Sort y and weights in y order; break ties according to x.
  utils::sort_all(y, x, weights);
 
  // 2.2 Count pairs and triplets of tied y.
  double pair_y = utils::count_tied_pairs(y, weights);
  double trip_y = utils::count_tied_triplets(y, weights);
  double v_y = utils::count_ties_v(y, weights);
 
  // 3. Calculate adjustment factor.
  double s = utils::sum(weights);
  double s2 = utils::perm_sum(weights, 2);
  double s3 = utils::perm_sum(weights, 3);
  double v_0 = 2 * s2 * (2 * s + 5);
  double v_1 = 2 * pair_x * 2 * pair_y / (2 * 2 * s2);
  double v_2 = 6 * trip_x * 6 * trip_y / (9 * 6 * s3);
  double v = (v_0 - v_x - v_y) / 18 + v_1 + v_2;
  return std::sqrt((s2 - pair_x) * (s2 - pair_y) / v);
}
 
}
 
}