Use Commons Statistics

docs/change_log.rst

@@ -30,6 +30,8 @@ Minor release of GDSC SMLM.
   The results can be added to the source image as an overlay.
 * Added the :numref:`{name} <results_plugins:Trace Viewer>` plugin to display a table of traced
   results.
+* Update from Commons Math to Commons Statistics for inference testing and univariate
+  descriptive statistics.
 
 
 Version 1.1

gdsc-smlm-ij/pom.xml

@@ -94,6 +94,14 @@ Software for single molecule localisation microscopy (SMLM) in ImageJ
       <artifactId>3D_Viewer</artifactId>
       <!-- For testing <version>4.0.3-SNAPSHOT</version> -->
     </dependency>
+    <dependency>
+      <groupId>org.apache.commons</groupId>
+      <artifactId>commons-statistics-descriptive</artifactId>
+    </dependency>
+    <dependency>
+      <groupId>org.apache.commons</groupId>
+      <artifactId>commons-statistics-inference</artifactId>
+    </dependency>
 
     <dependency>
       <groupId>uk.ac.sussex.gdsc</groupId>

gdsc-smlm-ij/src/main/java/uk/ac/sussex/gdsc/smlm/ij/plugins/CmosAnalysis.java

@@ -74,13 +74,15 @@
 import org.apache.commons.lang3.concurrent.ConcurrentRuntimeException;
 import org.apache.commons.lang3.time.StopWatch;
 import org.apache.commons.math3.stat.correlation.PearsonsCorrelation;
-import org.apache.commons.math3.stat.inference.TestUtils;
-import org.apache.commons.math3.stat.inference.WilcoxonSignedRankTest;
 import org.apache.commons.rng.UniformRandomProvider;
 import org.apache.commons.rng.sampling.distribution.ContinuousSampler;
 import org.apache.commons.rng.sampling.distribution.DiscreteSampler;
 import org.apache.commons.rng.sampling.distribution.NormalizedGaussianSampler;
 import org.apache.commons.rng.sampling.distribution.SharedStateContinuousSampler;
+import org.apache.commons.statistics.inference.KolmogorovSmirnovTest;
+import org.apache.commons.statistics.inference.SignificanceResult;
+import org.apache.commons.statistics.inference.TTest;
+import org.apache.commons.statistics.inference.WilcoxonSignedRankTest;
 import uk.ac.sussex.gdsc.core.data.IntegerType;
 import uk.ac.sussex.gdsc.core.data.SiPrefix;
 import uk.ac.sussex.gdsc.core.ij.HistogramPlot;
@@ -1213,25 +1215,24 @@ private void computeError(int slice, ImageStack simulationStack, WindowOrganiser
         .setPlotLabel(result.toString()).show(wo);
 
     // Kolmogorov–Smirnov test that the distributions are the same
-    double pvalue = TestUtils.kolmogorovSmirnovTest(x, y);
-    result.append(" : Kolmogorov–Smirnov p=").append(MathUtils.rounded(pvalue)).append(' ')
-        .append(((pvalue < 0.001) ? REJECT : ACCEPT));
+    SignificanceResult r = KolmogorovSmirnovTest.withDefaults().test(x, y);
+    result.append(" : Kolmogorov–Smirnov p=").append(MathUtils.rounded(r.getPValue())).append(' ')
+        .append((r.reject(0.001) ? REJECT : ACCEPT));
 
     if (slice == 3) {
       // Paired T-Test compares two related samples to assess whether their
       // population means differ.
       // T-Test is valid when the difference between the means is normally
       // distributed, e.g. gain
-      pvalue = TestUtils.pairedTTest(x, y);
-      result.append(" : Paired T-Test p=").append(MathUtils.rounded(pvalue)).append(' ')
-          .append(((pvalue < 0.001) ? REJECT : ACCEPT));
+      r = TTest.withDefaults().test(x, y);
+      result.append(" : Paired T-Test p=").append(MathUtils.rounded(r.getPValue())).append(' ')
+          .append((r.reject(0.001) ? REJECT : ACCEPT));
     } else {
       // Wilcoxon Signed Rank test compares two related samples to assess whether their
       // population mean ranks differ
-      final WilcoxonSignedRankTest wsrTest = new WilcoxonSignedRankTest();
-      pvalue = wsrTest.wilcoxonSignedRankTest(x, y, false);
-      result.append(" : Wilcoxon Signed Rank p=").append(MathUtils.rounded(pvalue)).append(' ')
-          .append(((pvalue < 0.001) ? REJECT : ACCEPT));
+      r = WilcoxonSignedRankTest.withDefaults().test(x, y);
+      result.append(" : Wilcoxon Signed Rank p=").append(MathUtils.rounded(r.getPValue()))
+          .append(' ').append((r.reject(0.001) ? REJECT : ACCEPT));
     }
 
     ImageJUtils.log(result.toString());

gdsc-smlm-ij/src/main/java/uk/ac/sussex/gdsc/smlm/ij/plugins/DensityImage.java

@@ -46,9 +46,9 @@
 import java.util.function.DoubleUnaryOperator;
 import java.util.stream.Stream;
 import org.apache.commons.lang3.concurrent.ConcurrentRuntimeException;
-import org.apache.commons.math3.stat.descriptive.SummaryStatistics;
 import org.apache.commons.rng.JumpableUniformRandomProvider;
 import org.apache.commons.rng.simple.RandomSource;
+import org.apache.commons.statistics.descriptive.Mean;
 import uk.ac.sussex.gdsc.core.clustering.DensityManager;
 import uk.ac.sussex.gdsc.core.data.utils.TypeConverter;
 import uk.ac.sussex.gdsc.core.ij.ImageJUtils;
@@ -591,20 +591,19 @@ private int[] cropBorder(MemoryPeakResults results, int[] density) {
    * @param density the density
    * @param radiusInPixels the radius
    * @param filtered the filtered
-   * @return the summary statistics
    */
-  private SummaryStatistics logDensityResults(MemoryPeakResults results, int[] density,
-      double radiusInPixels, int filtered) {
+  private void logDensityResults(MemoryPeakResults results, int[] density, double radiusInPixels,
+      int filtered) {
     final double region =
         radiusInPixels * radiusInPixels * ((settings.useSquareApproximation) ? 4 : Math.PI);
 
     final Rectangle bounds = results.getBounds();
     final double area = (double) bounds.width * bounds.height;
     final double expected = results.size() * region / area;
-    final SummaryStatistics summary = new SummaryStatistics();
+    final Mean summary = Mean.create();
 
     for (int i = 0; i < results.size(); i++) {
-      summary.addValue(density[i]);
+      summary.accept(density[i]);
     }
 
     final DensityManager dm = createDensityManager(results);
@@ -613,19 +612,17 @@ private SummaryStatistics logDensityResults(MemoryPeakResults results, int[] den
     final double l = (settings.useSquareApproximation) ? dm.ripleysLFunction(radiusInPixels)
         : dm.ripleysLFunction(density, radiusInPixels);
 
+    final double m = summary.getAsDouble();
     String msg = String.format(
         "Density %s : N=%d, %.0fpx^2 : Radius=%spx (%s%s): L(r) - r = %s : E = %s, Obs = %s (%sx)",
-        results.getName(), summary.getN(), area, rounded(radiusInPixels), rounded(settings.radius),
+        results.getName(), results.size(), area, rounded(radiusInPixels), rounded(settings.radius),
         UnitHelper.getShortName(DistanceUnit.forNumber(settings.distanceUnit)),
-        rounded(l - radiusInPixels), rounded(expected), rounded(summary.getMean()),
-        rounded(summary.getMean() / expected));
+        rounded(l - radiusInPixels), rounded(expected), rounded(m), rounded(m / expected));
     if (settings.filterLocalisations) {
       msg += String.format(" : Filtered=%d (%s%%)", filtered,
           rounded(filtered * 100.0 / density.length));
     }
     IJ.log(msg);
-
-    return summary;
   }
 
   private static String rounded(double value) {

gdsc-smlm-ij/src/main/java/uk/ac/sussex/gdsc/smlm/ij/plugins/FilterMolecules.java

@@ -38,8 +38,8 @@
 import java.util.concurrent.ExecutorService;
 import java.util.concurrent.Executors;
 import java.util.concurrent.atomic.AtomicReference;
-import org.apache.commons.math3.exception.MathIllegalArgumentException;
-import org.apache.commons.math3.stat.descriptive.rank.Percentile;
+import org.apache.commons.statistics.descriptive.Median;
+import org.apache.commons.statistics.descriptive.NaNPolicy;
 import uk.ac.sussex.gdsc.core.data.DataException;
 import uk.ac.sussex.gdsc.core.data.utils.ConversionException;
 import uk.ac.sussex.gdsc.core.data.utils.TypeConverter;
@@ -230,12 +230,12 @@ void run(MemoryPeakResults results) {
         p.getPrecision();
 
         // Precision in nm using the median
-        precision = new Percentile().evaluate(p.precisions, 50);
+        precision = Median.withDefaults().with(NaNPolicy.ERROR).evaluate(p.precisions);
         // Convert from nm to um to raw units
         final double rawPrecision = distanceConverter.convertBack(precision / 1e3);
         // Get the localisation error (4s^2) in units^2
         error = 4 * rawPrecision * rawPrecision;
-      } catch (final MathIllegalArgumentException | DataException ex) {
+      } catch (final IllegalArgumentException | DataException ex) {
         ImageJUtils.log(TITLE + " - Unable to compute precision: " + ex.getMessage());
       }
 

gdsc-smlm-ij/src/main/java/uk/ac/sussex/gdsc/smlm/ij/plugins/Fire.java

@@ -49,6 +49,7 @@
 import java.awt.geom.Rectangle2D;
 import java.util.Arrays;
 import java.util.Collection;
+import java.util.EnumSet;
 import java.util.LinkedList;
 import java.util.List;
 import java.util.concurrent.ExecutorService;
@@ -79,8 +80,10 @@
 import org.apache.commons.math3.optim.univariate.UnivariateObjectiveFunction;
 import org.apache.commons.math3.optim.univariate.UnivariateOptimizer;
 import org.apache.commons.math3.optim.univariate.UnivariatePointValuePair;
-import org.apache.commons.math3.stat.descriptive.DescriptiveStatistics;
 import org.apache.commons.math3.util.MathArrays;
+import org.apache.commons.statistics.descriptive.DoubleStatistics;
+import org.apache.commons.statistics.descriptive.Quantile;
+import org.apache.commons.statistics.descriptive.Statistic;
 import uk.ac.sussex.gdsc.core.data.DataException;
 import uk.ac.sussex.gdsc.core.data.utils.ConversionException;
 import uk.ac.sussex.gdsc.core.ij.HistogramPlot;
@@ -2515,34 +2518,40 @@ private PrecisionHistogram calculatePrecisionHistogram() {
     double yMax = 0;
 
     // Set the min and max y-values using 1.5 x IQR
-    final DescriptiveStatistics stats = precision.getStatistics();
-    final double lower = stats.getPercentile(25);
-    final double upper = stats.getPercentile(75);
+    final double[] values = precision.values();
+    final double[] q = Quantile.withDefaults().evaluate(values, 0.25, 0.75);
+    final double lower = q[0];
+    final double upper = q[1];
     final boolean logFitParameters = IJ.debugMode;
+    final DoubleStatistics stats = DoubleStatistics
+        .of(EnumSet.of(Statistic.STANDARD_DEVIATION, Statistic.MIN, Statistic.MAX), values);
+    final double min = stats.getAsDouble(Statistic.MIN);
+    final double max = stats.getAsDouble(Statistic.MAX);
     if (Double.isNaN(lower) || Double.isNaN(upper)) {
       if (logFitParameters) {
         ImageJUtils.log("Error computing IQR: %f - %f", lower, upper);
       }
     } else {
       final double iqr = upper - lower;
 
-      yMin = Math.max(lower - iqr, stats.getMin());
-      yMax = Math.min(upper + iqr, stats.getMax());
+      yMin = Math.max(lower - iqr, min);
+      yMax = Math.min(upper + iqr, max);
 
       if (logFitParameters) {
-        ImageJUtils.log("  Data range: %f - %f. Plotting 1.5x IQR: %f - %f", stats.getMin(),
-            stats.getMax(), yMin, yMax);
+        ImageJUtils.log("  Data range: %f - %f. Plotting 1.5x IQR: %f - %f", min, max, yMin, yMax);
       }
     }
 
     if (yMin == Double.NEGATIVE_INFINITY) {
       final int n = 5;
-      yMin = Math.max(stats.getMin(), stats.getMean() - n * stats.getStandardDeviation());
-      yMax = Math.min(stats.getMax(), stats.getMean() + n * stats.getStandardDeviation());
+      final double mean = stats.getAsDouble(Statistic.MEAN);
+      final double sd = stats.getAsDouble(Statistic.STANDARD_DEVIATION);
+      yMin = Math.max(min, mean - n * sd);
+      yMax = Math.min(min, mean + n * sd);
 
       if (logFitParameters) {
-        ImageJUtils.log("  Data range: %f - %f. Plotting mean +/- %dxSD: %f - %f", stats.getMin(),
-            stats.getMax(), n, yMin, yMax);
+        ImageJUtils.log("  Data range: %f - %f. Plotting mean +/- %dxSD: %f - %f", min, max, n,
+            yMin, yMax);
       }
     }
 

gdsc-smlm-ij/src/main/java/uk/ac/sussex/gdsc/smlm/ij/plugins/Optics.java

@@ -75,7 +75,7 @@
 import java.util.concurrent.atomic.AtomicInteger;
 import java.util.function.IntFunction;
 import org.apache.commons.lang3.tuple.Pair;
-import org.apache.commons.math3.stat.descriptive.rank.Percentile;
+import org.apache.commons.statistics.descriptive.Quantile;
 import uk.ac.sussex.gdsc.core.annotation.Nullable;
 import uk.ac.sussex.gdsc.core.clustering.optics.ClusteringResult;
 import uk.ac.sussex.gdsc.core.clustering.optics.DbscanResult;
@@ -1783,17 +1783,17 @@ public Pair<OpticsSettings, SettingsList> doWork(Pair<OpticsSettings, SettingsLi
 
         if (settings.getOpticsMode() == OpticsMode.FAST_OPTICS.ordinal()) {
           // The profile may be very high. Compute the outliers and remove.
-          final Percentile p = new Percentile();
-          p.setData(profile);
+          final Quantile p = Quantile.withDefaults();
           double max;
           final boolean useIqr = true;
           if (useIqr) {
-            final double lq = p.evaluate(25);
-            final double uq = p.evaluate(75);
+            final double[] q = p.evaluate(profile, 0.25, 0.75);
+            final double lq = q[0];
+            final double uq = q[1];
             max = (uq - lq) * 2 + uq;
           } else {
             // Remove top 2%
-            max = p.evaluate(98);
+            max = p.evaluate(profile, 0.98);
           }
           if (limits[1] > max) {
             limits[1] = max;

gdsc-smlm-ij/src/main/java/uk/ac/sussex/gdsc/smlm/ij/plugins/PsfCreator.java

@@ -68,7 +68,7 @@
 import java.util.concurrent.atomic.AtomicReference;
 import org.apache.commons.lang3.exception.ExceptionUtils;
 import org.apache.commons.math3.analysis.interpolation.LoessInterpolator;
-import org.apache.commons.math3.stat.descriptive.DescriptiveStatistics;
+import org.apache.commons.statistics.descriptive.Quantile;
 import uk.ac.sussex.gdsc.core.data.DataException;
 import uk.ac.sussex.gdsc.core.data.FloatStackTrivalueProvider;
 import uk.ac.sussex.gdsc.core.data.procedures.FloatStackTrivalueProcedure;
@@ -848,9 +848,9 @@ public void execute(PeakResult peak) {
    */
   private static double[] getLimits(double[] data) {
     final double[] limits = MathUtils.limits(data);
-    final DescriptiveStatistics stats = new DescriptiveStatistics(data);
-    final double lower = stats.getPercentile(25);
-    final double upper = stats.getPercentile(75);
+    final double[] q = Quantile.withDefaults().evaluate(data, 0.25, 0.75);
+    final double lower = q[0];
+    final double upper = q[1];
     final double iqr = (upper - lower) * 2;
     limits[0] = Math.max(lower - iqr, limits[0]);
     limits[1] = Math.min(upper + iqr, limits[1]);