| /* |
| * R : A Computer Language for Statistical Data Analysis |
| * Copyright (C) 1999-2016 The R Core Team. |
| * |
| * This program is free software; you can redistribute it and/or modify |
| * it under the terms of the GNU General Public License as published by |
| * the Free Software Foundation; either version 2 of the License, or |
| * (at your option) any later version. |
| * |
| * This program is distributed in the hope that it will be useful, |
| * but WITHOUT ANY WARRANTY; without even the implied warranty of |
| * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
| * GNU General Public License for more details. |
| * |
| * You should have received a copy of the GNU General Public License |
| * along with this program; if not, a copy is available at |
| * https://www.R-project.org/Licenses/ |
| */ |
| |
| /* Kendall's rank correlation tau and its exact distribution in case of no ties |
| */ |
| |
| #include <string.h> |
| #include <R.h> |
| #include <math.h> // for floor |
| #include <Rmath.h> |
| |
| /* |
| and the exact distribution of T = (n * (n - 1) * tau + 1) / 4, |
| which is -- if there are no ties -- the number of concordant ordered pairs |
| */ |
| |
| static double |
| ckendall(int k, int n, double **w) |
| { |
| int i, u; |
| double s; |
| |
| u = (n * (n - 1) / 2); |
| if ((k < 0) || (k > u)) return(0); |
| if (w[n] == 0) { |
| w[n] = (double *) R_alloc(u + 1, sizeof(double)); |
| memset(w[n], '\0', sizeof(double) * (u+1)); |
| for (i = 0; i <= u; i++) w[n][i] = -1; |
| } |
| if (w[n][k] < 0) { |
| if (n == 1) |
| w[n][k] = (k == 0); |
| else { |
| s = 0; |
| for (i = 0; i < n; i++) |
| s += ckendall(k - i, n - 1, w); |
| w[n][k] = s; |
| } |
| } |
| return(w[n][k]); |
| } |
| |
| #if 0 |
| void |
| dkendall(int *len, double *x, int *n) |
| { |
| int i; |
| double **w; |
| |
| w = (double **) R_alloc(*n + 1, sizeof(double *)); |
| |
| for (i = 0; i < *len; i++) |
| if (fabs(x[i] - floor(x[i] + 0.5)) > 1e-7) { |
| x[i] = 0; |
| } else { |
| x[i] = ckendall((int)x[i], *n) / gammafn(*n + 1, w); |
| } |
| } |
| #endif |
| |
| static void |
| pkendall(int len, double *Q, double *P, int n) |
| { |
| int i, j; |
| double p, q; |
| double **w; |
| |
| w = (double **) R_alloc(n + 1, sizeof(double *)); |
| memset(w, '\0', sizeof(double*) * (n+1)); |
| |
| for (i = 0; i < len; i++) { |
| q = floor(Q[i] + 1e-7); |
| if (q < 0) |
| P[i] = 0; |
| else if (q > (n * (n - 1) / 2)) |
| P[i] = 1; |
| else { |
| p = 0; |
| for (j = 0; j <= q; j++) p += ckendall(j, n, w); |
| P[i] = p / gammafn(n + 1); |
| } |
| } |
| } |
| |
| #include <Rinternals.h> |
| SEXP pKendall(SEXP q, SEXP sn) |
| { |
| q = PROTECT(coerceVector(q, REALSXP)); |
| int len = LENGTH(q), n = asInteger(sn); |
| SEXP p = PROTECT(allocVector(REALSXP, len)); |
| pkendall(len, REAL(q), REAL(p), n); |
| UNPROTECT(2); |
| return p; |
| } |