[
https://issues.apache.org/jira/browse/STATISTICS-85?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=17854364#comment-17854364
]
Alex Herbert commented on STATISTICS-85:
----------------------------------------
h2. Median
Note for the median that an odd length input array requires finding one index
(k = n/2); an even length input requires finding indices (k-1, k). Data were
created using a given length and extended with an additional range of 1 (e.g.
n=[1000, 1001]). For the BM data all distributions at n and n+1 are created.
For random data the length is randomly selected from the range. This creates a
50:50 mixture of odd and even length data. The data samples are randomly
shuffled per iteration to reduce branch prediction effects.
Tested on the BM distribution.
||n||CM3||CM4||JDK||Statistics||
|1000|3568053.381|4778367.635|7302762.526|2295992.742|
|Relative|0.489|0.654|1.0|0.314|
|100000|6242554116|6350143191|1375191908|314671031.5|
|Relative|4.539|4.618|1.0|0.229|
The CM implementation handles the small data. For larger data it is
significantly worse than a full sort. The Statistics implementation is robust
to the BM data distributions and 3-4x faster than a sort.
Tested using random data.
||n||CM3||CM4||JDK||Statistics||Relative to JDK||Relative to CM3||
|100|4.289E+04|4.277E+04|1.277E+05|3.424E+04|0.27|0.80|
|1000|1.007E+06|9.976E+05|2.900E+06|6.872E+05|0.24|0.68|
|10000|1.020E+07|1.001E+07|3.944E+07|7.215E+06|0.18|0.71|
|100000|9.613E+07|9.483E+07|5.040E+08|6.346E+07|0.13|0.66|
!MedianRandomData.png|width=1000,height=800!
On random data the CM implementation is faster than sorting. The statistics
implementation is consistently faster than CM and the relative difference
increases with length.
> Quantile implementation
> -----------------------
>
> Key: STATISTICS-85
> URL: https://issues.apache.org/jira/browse/STATISTICS-85
> Project: Commons Statistics
> Issue Type: New Feature
> Components: descriptive
> Reporter: Alex Herbert
> Priority: Major
> Attachments: 100QuantilesRandomDataLength10000.png,
> MedianRandomData.png
>
>
> Add a quantile implementation. This will interpolate the value of a sorted
> array of data for probability p in [0, 1].
> Replace the legacy API from Commons Math Percentile with an updated API. The
> new API should:
> * Decouple estimation of quantile positions inside data of length n; and the
> selection of correctly-ordered indices in array data.
> * Support multiple data types.
> * Support pre-sorted data.
> * Avoid performance issues observed in the CM Percentile implementation.
> h2. Proposed API
> {code:java}
> org.apache.commons.statistics.descriptive
> public final class Quantile {
> // overwrite=true; EstimationMethod.HF8; NaNPolicy.ERROR
> public static Quantile withDefaults();
> public Quantile withOverwrite(boolean);
> public Quantile with(EstimationMethod);
> // Could support NaN handling ... see below for NaNPolicy
> public Quantile with(NaNPolicy);
> // Create n uniform probabilities in range [p1, p2]
> public static double[] probabilities(int n);
> public static double[] probabilities(int n, double p1, double p2);
> // Quantiles on sorted data a of size n
> public double evaluate(int n, java.util.function.IntToDoubleFunction a,
> double p);
> public double[] evaluate(int n, java.util.function.IntToDoubleFunction a,
> double... p);
> // Quantiles on the primitive types that cannot be easily sorted
> public double evaluate(double[] a, double p);
> public double[] evaluate(double[] a, double... p);
> public double evaluate(int[] a, double p);
> public double[] evaluate(int[] a, double... p);
> public double evaluate(long[] a, double p);
> public double[] evaluate(long[] a, double... p);
> public double evaluate(float[] a, double p);
> public double[] evaluate(float[] a, double... p);
> // Provide the 9 methods in Hyndman and Fan (1996)
> // Sample Quantiles in Statistical Packages.
> // The American Statistician, 50, 361-365.
> public abstract class Quantile$EstimationMethod extends
> java.lang.Enum<Quantile$EstimationMethod> {
> public static final Quantile$EstimationMethod HF1;
> public static final Quantile$EstimationMethod HF2;
> public static final Quantile$EstimationMethod HF3;
> public static final Quantile$EstimationMethod HF4;
> public static final Quantile$EstimationMethod HF5;
> public static final Quantile$EstimationMethod HF6;
> public static final Quantile$EstimationMethod HF7;
> public static final Quantile$EstimationMethod HF8;
> public static final Quantile$EstimationMethod HF9;
> }
> }
> Note: The CM API used the 9 methods from Hyndman and Fann but labelled them
> as R1-9; this may be derived from the same names used in the R language. I
> propose to rename as HF1-9 to reflect the origin.
> {code}
> h2. NaNPolicy
> There are multiple options here. For reference R and Python's numpy only
> provide the option to exclude NaN:
> * R: quantile errors if NaN is present. median returns NaN. They is an
> option to exclude NaN.
> * numpy: two methods are provided: median/nanmedian + quantile/nanquantile
> (the non-nan versions will return NaN if any NaNs are present)
> Commons Math provides a remapping. Note the Statistics ranking module has the
> same NaNStrategy as that in CM:
> * MINIMAL: map to -infinity
> * MAXIMAL: map to +infinity
> * REMOVED: ignore from the data
> * FIXED: leave in place. This makes no sense for quantiles. It is done by
> moving to the end following the order imposed by Double.compare.
> * FAILED: raise an exception
> I favour the simpler option of: treating NaN so they are above/below all
> other values; removing them from the data; or raising an exception. I do not
> see the requirement to remap NaN to infinity. This can be done by the user.
> The API can be simplified by using:
> {code:java}
> public final class NaNPolicy extends java.lang.Enum<NaNPolicy> {
> public static final NaNPolicy LAST; // Move to end of data
> public static final NaNPolicy FIRST; // Move to start of data
> public static final NaNPolicy REMOVE; // Remove from data
> public static final NaNPolicy ERROR; // Raise exception
> }
> {code}
> Note that the FIRST option is not strictly required. But if there is an
> option to order the NaNs (i.e. LAST) then both orders should be supported.
> Which option to use as the default is not clear. As a drop in substitute for
> Arrays.sort the default would be handle NaN with a move to the end (LAST). As
> an API to signal to the user that the quantiles are possibly corrupted then
> ERROR would be the default. A user can then decide what to do if they receive
> errors during their analysis. Note that a separation of the Quantile API from
> the partitioning API (see below) may result in doing NaN processing twice
> introducing a performance overhead. If this will occur by default it should
> be documented so a user can override the NaN behaviour if they do not expect
> NaNs.
> Commons Math places the Percentile class in the rank package. I propose to
> move this implementation to place it in the descriptive package where it will
> sit beside other descriptive statistics for data such as mean and standard
> deviation.
> h2. Examples
> {code:java}
> double[] data = ...
> double q = Quantile.withDefaults().evaluate(data, 0.25);
> double[] quantiles = Quantile.withDefaults().evaluate(data, 0.25, 0.5, 0.75);
> // Use cases:
> // 1. Median / other quantile estimation
> // 2. Box plot of data
> // 3. Interquartile range analysis
> double[] p = Quantile.probabilities(100, 0.01, 0.99);
> double[] quantiles = Quantile.withDefaults().evaluate(data, p);
> // Use cases:
> // 1. plot p vs quantiles
> // 2. use p with an expected (inverse) probability distribution to create a
> QQ plot
> // Sorted input / unsupported datatype
> short[] data = ...
> Arrays.sort(data);
> double[] quantiles = Quantile.withDefaults().evaluate(data.length, i ->
> data[i], p);
> {code}
> h2. Implementation
> The Quantile API and the underlying implementation can be separated.
> Performing quantile estimation requires knowing the value of elements (i,
> i+1) in a sorted array a. Note that i+1 is used for interpolation:
> {noformat}
> value = a[i] + alpha * (a[i+1] - a[i]) ; alpha in [0, 1)
> {noformat}
> Note: The alpha factor and the index i are chosen from the percentile p using
> the EstimationMethod.
> The array does not have to be fully sorted. It can be partitioned so that
> each required element i is equal to the value in a fully sorted array. A
> partitioning API can be placed in Commons Numbers:
> {code:java}
> org.apache.commons.numbers.arrays
> public final class Selection {
> // Operate in-place and support a range as per Arrays.sort
> public void select(double[] a, int k);
> public void select(double[] a, int... k);
> public void select(int from, int to, double[] a, int k);
> public void select(int from, int to, double[] a, int... k);
> // Extend to other primitive types where sorting is slow
> // Sorting of 16-bit data is fast using a histogram so selection
> // has no major speed gain above a small length:
> // byte[] does counting sort at length 64 (temurin JDK 11)
> // short[]/char[] at length 1750 (temurin JDK 11)
> }
> {code}
> h3. Examples
> {code:java}
> // Find the bottom 100 (with 0-based indexing)
> Selection.select(data, 99);
> double[] bottom100 = Arrays.copyOf(data, 100);
> // Find the bottom and top 10
> Selection.select(data, 9, data.length - 10);
> double[] bottom10 = Arrays.copyOf(data, 10);
> double[] top10 = Arrays.copyOfRange(data, data.length - 10, data.length);
> {code}
--
This message was sent by Atlassian Jira
(v8.20.10#820010)
