1 // double centering utilities for the energy package
2 //
3 // Maria L. Rizzo <mrizzo@bgsu.edu>
4 // August, 2016
5
6
7
8 #include <Rcpp.h>
9 using namespace Rcpp;
10
11 NumericMatrix D_center(NumericMatrix Dx);
12 NumericMatrix U_center(NumericMatrix Dx);
13
14 // [[Rcpp::export]]
D_center(NumericMatrix Dx)15 NumericMatrix D_center(NumericMatrix Dx) {
16 /*
17 computes the double centered distance matrix for distance matrix Dx
18 for dCov, dCor, etc.
19 a_{ij} - a_{i.}/n - a_{.j}/n + a_{..}/n^2, all i, j
20 */
21 int j, k;
22 int n = Dx.nrow();
23 NumericVector akbar(n);
24 NumericMatrix A(n, n);
25 double abar = 0.0;
26
27 for (k=0; k<n; k++) {
28 akbar(k) = 0.0;
29 for (j=0; j<n; j++) {
30 akbar(k) += Dx(k, j);
31 }
32 abar += akbar(k);
33 akbar(k) /= (double) (n);
34 }
35 abar /= (double) (n * n);
36
37 for (k=0; k<n; k++)
38 for (j=k; j<n; j++) {
39 A(k, j) = Dx(k, j) - akbar(k) - akbar(j) + abar;
40 A(j, k) = A(k, j);
41 }
42
43 return A;
44 }
45
46 // [[Rcpp::export]]
U_center(NumericMatrix Dx)47 NumericMatrix U_center(NumericMatrix Dx) {
48 /*
49 computes the A_{kl}^U distances from the distance matrix (Dx_{kl}) for dCov^U
50 U-centering: if Dx = (a_{ij}) then compute U-centered A^U using
51 a_{ij} - a_{i.}/(n-2) - a_{.j}/(n-2) + a_{..}/((n-1)(n-2)), i \neq j
52 and zero diagonal
53 */
54 int j, k;
55 int n = Dx.nrow();
56 NumericVector akbar(n);
57 NumericMatrix A(n, n);
58 double abar = 0.0;
59
60 for (k=0; k<n; k++) {
61 akbar(k) = 0.0;
62 for (j=0; j<n; j++) {
63 akbar(k) += Dx(k, j);
64 }
65 abar += akbar(k);
66 akbar(k) /= (double) (n-2);
67 }
68 abar /= (double) ((n-1)*(n-2));
69
70 for (k=0; k<n; k++) {
71 for (j=k; j<n; j++) {
72 A(k, j) = Dx(k, j) - akbar(k) - akbar(j) + abar;
73 A(j, k) = A(k, j);
74 }
75 }
76 /* diagonal is zero */
77 for (k=0; k<n; k++)
78 A(k, k) = 0.0;
79
80 return A;
81 }
82
83