gpuMi, gpuRowTtest, gpuCorPerm, Spearman ranks: the method fixes from widestats

NEWS.md

@@ -1,5 +1,27 @@
 # gputools 2.1.1
 
+- `gpuMi()` gave `NaN` (or worse) for a constant variable, built its
+  B-spline weights in single precision when `precision = "single"` (a
+  value on a knot, as discrete data put them, could land in the
+  neighbouring bin), accepted `splineOrder > bins`, and computed every
+  pair's joint distribution in a per-pair loop. It now builds the weights
+  in double from the double data, takes every pair's joint distribution
+  from one matrix product per tile, computes the information in log form
+  (no large entropies subtracted; never below 0), returns exactly 0 for a
+  constant variable, checks `1 <= splineOrder <= bins`, and defaults to
+  Daub et al.'s `bins = 10, splineOrder = 3` instead of plain binning into
+  two bins. Its help page no longer claims the data must lie in [0, 1].
+- `gpuRowTtest()` with equal groups could count a permutation equivalent
+  to the observed labelling as more extreme, or fail to count its swap,
+  giving p-values too small for strongly separated rows: labellings are
+  now canonical and the observed t is recomputed with every block of
+  permutations. A row constant within each group with the means apart is
+  `Inf`/`-Inf` with p 0 rather than `NaN`; rows with no test no longer
+  enter the maxT maxima.
+- `gpuCorPerm()`'s rounding allowance is relative, as `gpuRowTtest()`'s.
+- Spearman ranks (`gpuCor`, `gpuCorPairs`, `gpuCorPerm`) are taken in
+  double before any narrowing to single, so distinct values no longer tie.
+
 - `gpuGranger()` reported F = 0, p = 1, "no causality", for a pair it
   could not resolve (collinear lags, a constant series, a cause identical
   to the effect), computed in single precision by default on uncentred

R/gpuCorPerm.R

@@ -41,7 +41,7 @@ gpuCorPerm <- function(x, y = NULL, B = 1000L, method = "pearson",
     perms <- vapply(seq_len(B), function(b) sample.int(n), integer(n)) - 1L
     dim(perms) <- c(n, B)
     # a permuted correlation equal to the observed one counts as at least
-    # as extreme; this allows for rounding in the comparison
+    # as extreme; a relative allowance for rounding in the comparison
     eps <- if(precision == "double") 1e-10 else 1e-5
     ans <- .Call(C_gpuCorPerm, x, y, perms, eps, precision)
 
@@ -52,7 +52,7 @@ gpuCorPerm <- function(x, y = NULL, B = 1000L, method = "pearson",
     # Westfall and Young's single-step maxT: how often the permutation's
     # largest |r| over all the tests reaches this one
     maxima <- sort(ans$maxima)
-    reach <- B - findInterval(abs(r) - eps, maxima, left.open = TRUE)
+    reach <- B - findInterval(abs(r) * (1 - eps), maxima, left.open = TRUE)
     adjusted <- (reach + 1) / (B + 1)
     p[is.na(r)] <- NA
     adjusted[is.na(r)] <- NA

R/gpuMi.R

@@ -1,37 +1,36 @@
-gpuMi <- function(x, y = NULL, bins = 2, splineOrder = 1,
-	precision = getOption("gputools.precision", "single"))
+gpuMi <- function(x, y = NULL, bins = 10, splineOrder = 3,
+    precision = getOption("gputools.precision", "single"))
 {
-	precision <- .precision(precision)
-	x <- as.matrix(x)
-	if(!is.null(y)) {
-		y <- as.matrix(y)
-	}
+    precision <- .precision(precision)
+    x <- as.matrix(x)
+    storage.mode(x) <- "double"
+    if (!is.null(y)) {
+        y <- as.matrix(y)
+        storage.mode(y) <- "double"
+        if (nrow(y) != nrow(x)) stop("x and y have different numbers of rows")
+    }
+    whole <- function(v) length(v) == 1 && is.numeric(v) && !is.na(v) &&
+        v == round(v)
+    if (!whole(bins) || !whole(splineOrder) || splineOrder < 1 ||
+        splineOrder > bins) {
+        stop("bins and splineOrder must be whole numbers with ",
+             "1 <= splineOrder <= bins")
+    }
+    if (bins > 64) stop("bins must be at most 64")
+    if (nrow(x) < 2) stop("need at least 2 rows")
+    if (anyNA(x) || (!is.null(y) && anyNA(y))) {
+        stop("x and y must have no missing values")
+    }
+    bins <- as.integer(bins)
+    splineOrder <- as.integer(splineOrder)
+    nsamples <- as.integer(nrow(x))
+    na <- as.integer(ncol(x))
+    nb <- if (is.null(y)) na else as.integer(ncol(y))
 
-	bins <- as.integer(bins)
-	splineOrder <- as.integer(splineOrder)
-
-	nsamples <- as.integer(nrow(x))
-	na <- as.integer(ncol(x))
-	a <- as.double(x)
-
-	b <- a
-	nb <- na
-	row_labels <- colnames(x)
-
-	if(is.null(y)) {
-		nb <- as.integer(ncol(x))
-		b <- as.double(x)
-		row_labels <- colnames(x)
-	} else {
-		nb <- as.integer(ncol(y))
-		b <- as.double(y)
-		row_labels <- colnames(y)
-	}
-
-	mi <- .Call(C_gpuMutualInfo, bins, splineOrder, nsamples, a, na, b, nb,
-		precision)
-	mi <- matrix(mi, nb, na)
-	rownames(mi) <- row_labels
-	colnames(mi) <- colnames(x)
-	return(mi)
+    mi <- .Call(C_gpuMutualInfo, bins, splineOrder, nsamples, x, na, y, nb,
+                precision)
+    mi <- matrix(mi, nb, na)
+    rownames(mi) <- if (is.null(y)) colnames(x) else colnames(y)
+    colnames(mi) <- colnames(x)
+    mi
 }

R/gpuRowTtest.R

@@ -40,17 +40,34 @@ gpuRowTtest <- function(x, group, var.equal = FALSE, B = 0L,
         se <- sqrt(a + b)
         df <- (a + b)^2 / (a^2 / (n1 - 1) + b^2 / (n2 - 1))
     }
+    # a zero standard error with the means apart is a complete separation
+    # of the groups: t infinite, p 0 (as SciPy reports it); with the means
+    # equal as well, or a constant row, there is no test
     t <- (m1 - m2) / se
-    t[se == 0] <- NaN
+    t[se == 0 & m1 == m2] <- NaN
+    p.value <- 2 * pt(-abs(t), df)
+    # (Welch's df is 0/0 there, so the p-value is set directly)
+    p.value[is.infinite(t)] <- 0
+    if(!var.equal) df[is.infinite(t)] <- NA
     res <- data.frame(mean1 = m1, mean2 = m2, statistic = t, df = df,
-                      p.value = 2 * pt(-abs(t), df), row.names = rownames(x))
+                      p.value = p.value, row.names = rownames(x))
 
     if(B > 0) {
         # the permutations, drawn here so set.seed() reproduces them, on
         # the GPU: every row's t under every relabelling at once
         labels <- matrix(0L, n, B + 1L)
         labels[, 1] <- in1
         for(b in seq_len(B)) labels[, b + 1L] <- in1[sample.int(n)]
+        if(n1 == n2) {
+            # With equal groups a labelling and its swap give the same |t|
+            # in exact arithmetic but are computed by different roundings,
+            # which for strongly separated rows can differ by more than the
+            # allowance: a swap of the observed labelling could then fail
+            # to count, making p-values too small. Putting sample 1 in
+            # group 1 always makes the two the same column, bit for bit.
+            flip <- labels[1, ] == 0L
+            labels[, flip] <- 1L - labels[, flip]
+        }
         eps <- if(precision == "double") 1e-10 else 1e-5
         ans <- .Call(C_gpuRowTtest, x, labels, n1, isTRUE(var.equal), eps,
                      precision)