1 #include <Rcpp.h>
2 using namespace Rcpp;
3
4 NumericMatrix U_center(NumericMatrix);
5
6 //[[Rcpp::export]]
dcovU_stats(NumericMatrix Dx,NumericMatrix Dy)7 NumericVector dcovU_stats(NumericMatrix Dx, NumericMatrix Dy) {
8 // x and y must be square distance matrices
9 NumericMatrix A = U_center(Dx);
10 NumericMatrix B = U_center(Dy);
11 double ab = 0.0, aa = 0.0, bb = 0.0;
12 double V, dcorU = 0.0;
13 double eps = std::numeric_limits<double>::epsilon(); //machine epsilon
14 int n = Dx.nrow();
15 int n2 = n * (n - 3);
16
17 for (int i=0; i<n; i++)
18 for (int j=0; j<i; j++) {
19 // U-centered is symmetric, with zero diagonal
20 ab += A(i, j) * B(i, j);
21 aa += A(i, j) * A(i, j);
22 bb += B(i, j) * B(i, j);
23 }
24 ab = 2.0 * ab / (double) n2;
25 aa = 2.0 * aa / (double) n2;
26 bb = 2.0 * bb / (double) n2;
27 V = aa * bb;
28 if (V > eps)
29 dcorU = ab / sqrt(V);
30
31 return NumericVector::create(
32 _["dCovU"] = ab,
33 _["bcdcor"] = dcorU,
34 _["dVarXU"] = aa,
35 _["dVarYU"] = bb
36 );
37 }
38
39