001/*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements.  See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License.  You may obtain a copy of the License at
008 *
009 *      http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017package org.apache.commons.math4.stat.descriptive.rank;
018
019import java.io.Serializable;
020import java.util.Arrays;
021import java.util.BitSet;
022
023import org.apache.commons.math4.exception.MathIllegalArgumentException;
024import org.apache.commons.math4.exception.NotPositiveException;
025import org.apache.commons.math4.exception.NumberIsTooLargeException;
026import org.apache.commons.math4.exception.OutOfRangeException;
027import org.apache.commons.math4.exception.util.LocalizedFormats;
028import org.apache.commons.math4.stat.descriptive.AbstractUnivariateStatistic;
029import org.apache.commons.math4.stat.ranking.NaNStrategy;
030import org.apache.commons.math4.util.FastMath;
031import org.apache.commons.math4.util.KthSelector;
032import org.apache.commons.math4.util.MathArrays;
033import org.apache.commons.math4.util.MathUtils;
034import org.apache.commons.math4.util.MedianOf3PivotingStrategy;
035import org.apache.commons.math4.util.PivotingStrategyInterface;
036import org.apache.commons.numbers.core.Precision;
037
038/**
039 * Provides percentile computation.
040 * <p>
041 * There are several commonly used methods for estimating percentiles (a.k.a.
042 * quantiles) based on sample data.  For large samples, the different methods
043 * agree closely, but when sample sizes are small, different methods will give
044 * significantly different results.  The algorithm implemented here works as follows:
045 * <ol>
046 * <li>Let <code>n</code> be the length of the (sorted) array and
047 * <code>0 &lt; p &lt;= 100</code> be the desired percentile.</li>
048 * <li>If <code> n = 1 </code> return the unique array element (regardless of
049 * the value of <code>p</code>); otherwise </li>
050 * <li>Compute the estimated percentile position
051 * <code> pos = p * (n + 1) / 100</code> and the difference, <code>d</code>
052 * between <code>pos</code> and <code>floor(pos)</code> (i.e. the fractional
053 * part of <code>pos</code>).</li>
054 * <li> If <code>pos &lt; 1</code> return the smallest element in the array.</li>
055 * <li> Else if <code>pos &gt;= n</code> return the largest element in the array.</li>
056 * <li> Else let <code>lower</code> be the element in position
057 * <code>floor(pos)</code> in the array and let <code>upper</code> be the
058 * next element in the array.  Return <code>lower + d * (upper - lower)</code>
059 * </li>
060 * </ol>
061 * <p>
062 * To compute percentiles, the data must be at least partially ordered.  Input
063 * arrays are copied and recursively partitioned using an ordering definition.
064 * The ordering used by <code>Arrays.sort(double[])</code> is the one determined
065 * by {@link java.lang.Double#compareTo(Double)}.  This ordering makes
066 * <code>Double.NaN</code> larger than any other value (including
067 * <code>Double.POSITIVE_INFINITY</code>).  Therefore, for example, the median
068 * (50th percentile) of
069 * <code>{0, 1, 2, 3, 4, Double.NaN}</code> evaluates to <code>2.5.</code></p>
070 * <p>
071 * Since percentile estimation usually involves interpolation between array
072 * elements, arrays containing  <code>NaN</code> or infinite values will often
073 * result in <code>NaN</code> or infinite values returned.</p>
074 * <p>
075 * Further, to include different estimation types such as R1, R2 as mentioned in
076 * <a href="http://en.wikipedia.org/wiki/Quantile">Quantile page(wikipedia)</a>,
077 * a type specific NaN handling strategy is used to closely match with the
078 * typically observed results from popular tools like R(R1-R9), Excel(R7).</p>
079 * <p>
080 * Since 2.2, Percentile uses only selection instead of complete sorting
081 * and caches selection algorithm state between calls to the various
082 * {@code evaluate} methods. This greatly improves efficiency, both for a single
083 * percentile and multiple percentile computations. To maximize performance when
084 * multiple percentiles are computed based on the same data, users should set the
085 * data array once using either one of the {@link #evaluate(double[], double)} or
086 * {@link #setData(double[])} methods and thereafter {@link #evaluate(double)}
087 * with just the percentile provided.
088 * </p>
089 * <p>
090 * <strong>Note that this implementation is not synchronized.</strong> If
091 * multiple threads access an instance of this class concurrently, and at least
092 * one of the threads invokes the <code>increment()</code> or
093 * <code>clear()</code> method, it must be synchronized externally.</p>
094 */
095public class Percentile extends AbstractUnivariateStatistic implements Serializable {
096    /** Serializable version identifier */
097    private static final long serialVersionUID = 20150412L;
098
099    /** Maximum number of partitioning pivots cached (each level double the number of pivots). */
100    private static final int MAX_CACHED_LEVELS = 10;
101
102    /** Maximum number of cached pivots in the pivots cached array */
103    private static final int PIVOTS_HEAP_LENGTH = 0x1 << MAX_CACHED_LEVELS - 1;
104
105    /** Default KthSelector used with default pivoting strategy */
106    private final KthSelector kthSelector;
107
108    /** Any of the {@link EstimationType}s such as {@link EstimationType#LEGACY CM} can be used. */
109    private final EstimationType estimationType;
110
111    /** NaN Handling of the input as defined by {@link NaNStrategy} */
112    private final NaNStrategy nanStrategy;
113
114    /**
115     * Determines what percentile is computed when evaluate() is activated
116     * with no quantile argument.
117     */
118    private double quantile;
119
120    /** Cached pivots. */
121    private int[] cachedPivots;
122
123    /** Weights*/
124    private double[] weights;
125    /**
126     * Constructs a Percentile with the following defaults.
127     * <ul>
128     *   <li>default quantile: 50.0, can be reset with {@link #setQuantile(double)}</li>
129     *   <li>default estimation type: {@link EstimationType#LEGACY},
130     *   can be reset with {@link #withEstimationType(EstimationType)}</li>
131     *   <li>default NaN strategy: {@link NaNStrategy#REMOVED},
132     *   can be reset with {@link #withNaNStrategy(NaNStrategy)}</li>
133     *   <li>a KthSelector that makes use of {@link MedianOf3PivotingStrategy},
134     *   can be reset with {@link #withKthSelector(KthSelector)}</li>
135     * </ul>
136     */
137    public Percentile() {
138        // No try-catch or advertised exception here - arg is valid
139        this(50.0);
140    }
141
142    /**
143     * Constructs a Percentile with the specific quantile value and the following
144     * <ul>
145     *   <li>default method type: {@link EstimationType#LEGACY}</li>
146     *   <li>default NaN strategy: {@link NaNStrategy#REMOVED}</li>
147     *   <li>a Kth Selector : {@link KthSelector}</li>
148     * </ul>
149     * @param quantile the quantile
150     * @throws MathIllegalArgumentException  if p is not greater than 0 and less
151     * than or equal to 100
152     */
153    public Percentile(final double quantile) {
154        this(quantile, EstimationType.LEGACY, NaNStrategy.REMOVED,
155             new KthSelector(new MedianOf3PivotingStrategy()));
156    }
157
158    /**
159     * Copy constructor, creates a new {@code Percentile} identical
160     * to the {@code original}
161     *
162     * @param original the {@code Percentile} instance to copy.
163     * Cannot be {@code null}.
164     */
165    public Percentile(final Percentile original) {
166        MathUtils.checkNotNull(original);
167        estimationType   = original.getEstimationType();
168        nanStrategy      = original.getNaNStrategy();
169        kthSelector      = original.getKthSelector();
170
171        setData(original.getDataRef());
172        if (original.cachedPivots != null) {
173            System.arraycopy(original.cachedPivots, 0, cachedPivots, 0, original.cachedPivots.length);
174        }
175        setQuantile(original.quantile);
176
177    }
178
179    /**
180     * Constructs a Percentile with the specific quantile value,
181     * {@link EstimationType}, {@link NaNStrategy} and {@link KthSelector}.
182     *
183     * @param quantile the quantile to be computed
184     * @param estimationType one of the percentile {@link EstimationType  estimation types}
185     * @param nanStrategy one of {@link NaNStrategy} to handle with NaNs.
186     * Cannot be {@code null}.
187     * @param kthSelector a {@link KthSelector} to use for pivoting during search
188     * @throws MathIllegalArgumentException if p is not within (0,100]
189     */
190    protected Percentile(final double quantile,
191                         final EstimationType estimationType,
192                         final NaNStrategy nanStrategy,
193                         final KthSelector kthSelector) {
194        setQuantile(quantile);
195        cachedPivots = null;
196        MathUtils.checkNotNull(estimationType);
197        MathUtils.checkNotNull(nanStrategy);
198        MathUtils.checkNotNull(kthSelector);
199        this.estimationType = estimationType;
200        this.nanStrategy = nanStrategy;
201        this.kthSelector = kthSelector;
202    }
203
204    /** {@inheritDoc} */
205    @Override
206    public void setData(final double[] values) {
207        if (values == null) {
208            cachedPivots = null;
209        } else {
210            cachedPivots = new int[PIVOTS_HEAP_LENGTH];
211            Arrays.fill(cachedPivots, -1);
212        }
213        super.setData(values);
214    }
215    /**
216     * Set the data array.
217     * @param values Data array.
218     * Cannot be {@code null}.
219     * @param sampleWeights corresponding positive and non-NaN weights.
220     * Cannot be {@code null}.
221     * @throws MathIllegalArgumentException if lengths of values and weights are not equal.
222     * @throws org.apache.commons.math4.exception.NotANumberException if any weight is NaN
223     * @throws org.apache.commons.math4.exception.NotStrictlyPositiveException if any weight is not positive
224     */
225    public void setData(final double[] values,
226                        final double[] sampleWeights) {
227        setData(values, sampleWeights, 0, values.length);
228    }
229
230    /** {@inheritDoc} */
231    @Override
232    public void setData(final double[] values, final int begin, final int length) {
233        if (values == null) {
234            cachedPivots = null;
235        } else {
236            cachedPivots = new int[PIVOTS_HEAP_LENGTH];
237            Arrays.fill(cachedPivots, -1);
238        }
239        super.setData(values, begin, length);
240    }
241    /**
242     * Set the data and weights arrays.  The input array is copied, not referenced.
243     * @param values Data array.
244     * Cannot be {@code null}.
245     * @param sampleWeights corresponding positive and non-NaN weights.
246     * Cannot be {@code null}.
247     * @param begin the index of the first element to include
248     * @param length the number of elements to include
249     * @throws MathIllegalArgumentException if lengths of values and weights are not equal or values or weights is null
250     * @throws NotPositiveException if begin or length is not positive
251     * @throws NumberIsTooLargeException if begin + length is greater than values.length
252     * @throws org.apache.commons.math4.exception.NotANumberException if any weight is NaN
253     * @throws org.apache.commons.math4.exception.NotStrictlyPositiveException if any weight is not positive
254     */
255    public void setData(final double[] values,
256                        final double[] sampleWeights,
257                        final int begin,
258                        final int length) {
259        if (begin < 0) {
260            throw new NotPositiveException(LocalizedFormats.START_POSITION, begin);
261        }
262
263        if (length < 0) {
264            throw new NotPositiveException(LocalizedFormats.LENGTH, length);
265        }
266
267        if (begin + length > values.length) {
268            throw new NumberIsTooLargeException(LocalizedFormats.SUBARRAY_ENDS_AFTER_ARRAY_END,
269                                                begin + length, values.length, true);
270        }
271
272        if (sampleWeights == null) {
273            throw new MathIllegalArgumentException(LocalizedFormats.NULL_NOT_ALLOWED);
274        }
275        cachedPivots = new int[PIVOTS_HEAP_LENGTH];
276        Arrays.fill(cachedPivots, -1);
277
278        /** Check length */
279        if (values.length != sampleWeights.length) {
280            throw new MathIllegalArgumentException(LocalizedFormats.LENGTH,
281                                                   values, sampleWeights);
282        }
283        /** Check weights > 0 */
284        MathArrays.checkPositive(sampleWeights);
285        MathArrays.checkNotNaN(sampleWeights);
286
287        super.setData(values, begin, length);
288        weights = new double[length];
289        System.arraycopy(sampleWeights, begin, weights, 0, length);
290    }
291    /**
292     * Returns the result of evaluating the statistic over the stored data.
293     * If weights have been set, it will compute weighted percentile.
294     * <p>
295     * The stored array is the one which was set by previous calls to
296     * {@link #setData(double[])} or {@link #setData(double[], double[] sampleWeights, int begin, int length)}
297     * </p>
298     * @param p the percentile value to compute
299     * @return the value of the statistic applied to the stored data
300     * @throws MathIllegalArgumentException if lengths of values and weights are not equal or values or weights is null
301     * @throws NotPositiveException if begin, length is negative
302     * @throws org.apache.commons.math4.exception.NotStrictlyPositiveException if any weight is not positive
303     * @throws org.apache.commons.math4.exception.NotANumberException if any weight is NaN
304     * @throws OutOfRangeException if p is invalid
305     * @throws NumberIsTooLargeException if begin + length is greater than values.length
306     * (p must be greater than 0 and less than or equal to 100)
307     */
308    public double evaluate(final double p) {
309        if (weights == null) {
310            return evaluate(getDataRef(), p);
311        } else {
312            return evaluate(getDataRef(), weights, p);
313        }
314    }
315
316    /**
317     * Returns an estimate of the <code>p</code>th percentile of the values
318     * in the <code>values</code> array.
319     * <p>
320     * Calls to this method do not modify the internal <code>quantile</code>
321     * state of this statistic.</p>
322     * <ul>
323     * <li>Returns <code>Double.NaN</code> if <code>values</code> has length
324     * <code>0</code></li>
325     * <li>Returns (for any value of <code>p</code>) <code>values[0]</code>
326     *  if <code>values</code> has length <code>1</code></li>
327     * <li>Throws <code>MathIllegalArgumentException</code> if <code>values</code>
328     * is null or p is not a valid quantile value (p must be greater than 0
329     * and less than or equal to 100) </li>
330     * </ul>
331     * <p>
332     * See {@link Percentile} for a description of the percentile estimation
333     * algorithm used.</p>
334     *
335     * @param values input array of values
336     * @param p the percentile value to compute
337     * @return the percentile value or Double.NaN if the array is empty
338     * @throws MathIllegalArgumentException if <code>values</code> is null or p is invalid
339     */
340    public double evaluate(final double[] values, final double p) {
341        MathArrays.verifyValues(values, 0, 0);
342        return evaluate(values, 0, values.length, p);
343    }
344    /**
345     * Returns an estimate of the <code>p</code>th percentile of the values
346     * in the <code>values</code> array with their weights.
347     * <p>
348     * See {@link Percentile} for a description of the percentile estimation
349     * algorithm used.</p>
350     * @param values input array of values
351     * @param sampleWeights weights of values
352     * @param p the percentile value to compute
353     * @return the weighted percentile value or Double.NaN if the array is empty
354     * @throws MathIllegalArgumentException if lengths of values and weights are not equal or values or weights is null
355     * @throws NotPositiveException if begin, length is negative
356     * @throws org.apache.commons.math4.exception.NotStrictlyPositiveException if any weight is not positive
357     * @throws org.apache.commons.math4.exception.NotANumberException if any weight is NaN
358     * @throws OutOfRangeException if p is invalid
359     * @throws NumberIsTooLargeException if begin + length is greater than values.length
360     */
361    public double evaluate(final double[] values, final double[] sampleWeights, final double p) {
362        MathArrays.verifyValues(values, 0, 0);
363        MathArrays.verifyValues(sampleWeights, 0, 0);
364        return evaluate(values, sampleWeights, 0, values.length, p);
365    }
366
367    /**
368     * Returns an estimate of the <code>quantile</code>th percentile of the
369     * designated values in the <code>values</code> array.  The quantile
370     * estimated is determined by the <code>quantile</code> property.
371     * <ul>
372     * <li>Returns <code>Double.NaN</code> if <code>length = 0</code></li>
373     * <li>Returns (for any value of <code>quantile</code>)
374     * <code>values[begin]</code> if <code>length = 1 </code></li>
375     * <li>Throws <code>MathIllegalArgumentException</code> if <code>values</code>
376     * is null, or <code>start</code> or <code>length</code> is invalid</li>
377     * </ul>
378     * <p>
379     * See {@link Percentile} for a description of the percentile estimation
380     * algorithm used.</p>
381     *
382     * @param values the input array
383     * @param start index of the first array element to include
384     * @param length the number of elements to include
385     * @return the percentile value
386     * @throws MathIllegalArgumentException if the parameters are not valid
387     *
388     */
389    @Override
390    public double evaluate(final double[] values, final int start, final int length) {
391        return evaluate(values, start, length, quantile);
392    }
393    /**
394     * Returns an estimate of the weighted <code>quantile</code>th percentile of the
395     * designated values in the <code>values</code> array.  The quantile
396     * estimated is determined by the <code>quantile</code> property.
397     * <p>
398     * See {@link Percentile} for a description of the percentile estimation
399     * algorithm used.</p>
400     *
401     * @param values the input array
402     * @param sampleWeights the weights of values
403     * @param start index of the first array element to include
404     * @param length the number of elements to include
405     * @return the percentile value
406     * @throws MathIllegalArgumentException if lengths of values and weights are not equal or values or weights is null
407     * @throws NotPositiveException if begin, length is negative
408     * @throws org.apache.commons.math4.exception.NotStrictlyPositiveException if any weight is not positive
409     * @throws org.apache.commons.math4.exception.NotANumberException if any weight is NaN
410     * @throws OutOfRangeException if p is invalid
411     * @throws NumberIsTooLargeException if begin + length is greater than values.length
412     */
413    public double evaluate(final double[] values, final double[] sampleWeights,
414                           final int start, final int length) {
415        return evaluate(values, sampleWeights, start, length, quantile);
416    }
417
418     /**
419     * Returns an estimate of the <code>p</code>th percentile of the values
420     * in the <code>values</code> array, starting with the element in (0-based)
421     * position <code>begin</code> in the array and including <code>length</code>
422     * values.
423     * <p>
424     * Calls to this method do not modify the internal <code>quantile</code>
425     * state of this statistic.</p>
426     * <ul>
427     * <li>Returns <code>Double.NaN</code> if <code>length = 0</code></li>
428     * <li>Returns (for any value of <code>p</code>) <code>values[begin]</code>
429     *  if <code>length = 1 </code></li>
430     * <li>Throws <code>MathIllegalArgumentException</code> if <code>values</code>
431     *  is null , <code>begin</code> or <code>length</code> is invalid, or
432     * <code>p</code> is not a valid quantile value (p must be greater than 0
433     * and less than or equal to 100)</li>
434     * </ul>
435     * <p>
436     * See {@link Percentile} for a description of the percentile estimation
437     * algorithm used.</p>
438     *
439     * @param values array of input values
440     * @param p the percentile to compute
441     * @param begin the first (0-based) element to include in the computation
442     * @param length the number of array elements to include
443     * @return the percentile value.
444     * @throws MathIllegalArgumentException if the parameters are not valid.
445     */
446    public double evaluate(final double[] values, final int begin,
447                           final int length, final double p) {
448        MathArrays.verifyValues(values, begin, length);
449        if (p > 100 || p <= 0) {
450            throw new OutOfRangeException(LocalizedFormats.OUT_OF_BOUNDS_QUANTILE_VALUE,
451                                          p, 0, 100);
452        }
453        if (length == 0) {
454            return Double.NaN;
455        }
456        if (length == 1) {
457            return values[begin]; // always return single value for n = 1
458        }
459
460        final double[] work = getWorkArray(values, begin, length);
461        final int[] pivotsHeap = getPivots(values);
462        return work.length == 0 ?
463            Double.NaN :
464            estimationType.evaluate(work, pivotsHeap, p, kthSelector);
465    }
466     /**
467     * Returns an estimate of the <code>p</code>th percentile of the values
468     * in the <code>values</code> array with <code>sampleWeights</code>, starting with the element in (0-based)
469     * position <code>begin</code> in the array and including <code>length</code>
470     * values.
471     * <p>
472     * See {@link Percentile} for a description of the percentile estimation
473     * algorithm used.</p>
474     *
475     * @param values array of input values
476     * @param sampleWeights positive and non-NaN weights of values
477     * @param begin  the first (0-based) element to include in the computation
478     * @param length  the number of array elements to include
479     * @param p  percentile to compute
480     * @return the weighted percentile value
481     * @throws MathIllegalArgumentException if lengths of values and weights are not equal or values or weights is null
482     * @throws NotPositiveException if begin, length is negative
483     * @throws org.apache.commons.math4.exception.NotStrictlyPositiveException if any weight is not positive
484     * @throws org.apache.commons.math4.exception.NotANumberException if any weight is NaN
485     * @throws OutOfRangeException if p is invalid
486     * @throws NumberIsTooLargeException if begin + length is greater than values.length
487     */
488    public double evaluate(final double[] values, final double[] sampleWeights, final int begin,
489                           final int length, final double p) {
490        /** Check length */
491        if (values.length != sampleWeights.length) {
492            throw new MathIllegalArgumentException(LocalizedFormats.LENGTH,
493                                                   values, sampleWeights);
494        }
495        if (values == null || sampleWeights == null) {
496            throw new MathIllegalArgumentException(LocalizedFormats.NULL_NOT_ALLOWED);
497        }
498        MathArrays.verifyValues(values, begin, length);
499        MathArrays.verifyValues(sampleWeights, begin, length);
500        MathArrays.checkPositive(sampleWeights);
501        MathArrays.checkNotNaN(sampleWeights);
502
503        if (p > 100 || p <= 0) {
504            throw new OutOfRangeException(LocalizedFormats.OUT_OF_BOUNDS_QUANTILE_VALUE,
505                                          p, 0, 100);
506        }
507        if (length == 0) {
508            return Double.NaN;
509        }
510        if (length == 1) {
511            // Always return single value for n = 1
512            return values[begin];
513        }
514
515        final double[] work = getWorkArray(values, begin, length);
516        final double[] workWeights = getWorkArray(values, sampleWeights, begin, length);
517        return work.length == 0 ? Double.NaN :
518                    estimationType.evaluate(work, workWeights, p);
519    }
520    /**
521     * Returns the value of the quantile field (determines what percentile is
522     * computed when evaluate() is called with no quantile argument).
523     *
524     * @return quantile set while construction or {@link #setQuantile(double)}
525     */
526    public double getQuantile() {
527        return quantile;
528    }
529
530    /**
531     * Sets the value of the quantile field (determines what percentile is
532     * computed when evaluate() is called with no quantile argument).
533     *
534     * @param p a value between 0 &lt; p &lt;= 100
535     * @throws MathIllegalArgumentException  if p is not greater than 0 and less
536     * than or equal to 100
537     */
538    public void setQuantile(final double p) {
539        if (p <= 0 || p > 100) {
540            throw new OutOfRangeException(LocalizedFormats.OUT_OF_BOUNDS_QUANTILE_VALUE,
541                                          p, 0, 100);
542        }
543        quantile = p;
544    }
545
546    /**
547     * {@inheritDoc}
548     */
549    @Override
550    public Percentile copy() {
551        return new Percentile(this);
552    }
553
554    /**
555     * Get the work array to operate. Makes use of prior {@code storedData} if
556     * it exists or else do a check on NaNs and copy a subset of the array
557     * defined by begin and length parameters. The set {@link #nanStrategy} will
558     * be used to either retain/remove/replace any NaNs present before returning
559     * the resultant array.
560     *
561     * @param values the array of numbers
562     * @param begin index to start reading the array
563     * @param length the length of array to be read from the begin index
564     * @return work array sliced from values in the range [begin,begin+length)
565     * @throws MathIllegalArgumentException if values or indices are invalid
566     */
567    private double[] getWorkArray(final double[] values, final int begin, final int length) {
568        final double[] work;
569        if (values == getDataRef()) {
570            work = getDataRef();
571        } else {
572            switch (nanStrategy) {
573                case MAXIMAL: // Replace NaNs with +INFs
574                    work = replaceAndSlice(values, begin, length, Double.NaN, Double.POSITIVE_INFINITY);
575                    break;
576                case MINIMAL: // Replace NaNs with -INFs
577                    work = replaceAndSlice(values, begin, length, Double.NaN, Double.NEGATIVE_INFINITY);
578                    break;
579                case REMOVED: // Drop NaNs from data
580                    work = removeAndSlice(values, begin, length, Double.NaN);
581                    break;
582                case FAILED: // NaN is not acceptable
583                    work = copyOf(values, begin, length);
584                    MathArrays.checkNotNaN(work);
585                    break;
586                default: // FIXED
587                    work = copyOf(values,begin,length);
588                    break;
589            }
590        }
591        return work;
592    }
593    /**
594     * Get the work arrays of weights to operate.
595     *
596     * @param values the array of numbers
597     * @param sampleWeights the array of weights
598     * @param begin index to start reading the array
599     * @param length the length of array to be read from the begin index
600     * @return work array sliced from values in the range [begin,begin+length)
601     */
602    protected double[] getWorkArray(final double[] values, final double[] sampleWeights,
603                                    final int begin, final int length) {
604        final double[] work;
605        if (values == getDataRef()) {
606            work = this.weights;
607        } else {
608            switch (nanStrategy) {
609                case REMOVED: // Drop weight if the data is NaN
610                    work = removeAndSliceByRef(values, sampleWeights, begin, length, Double.NaN);
611                    break;
612                default: // FIXED
613                    work = copyOf(sampleWeights, begin, length);
614                    break;
615            }
616        }
617        return work;
618    }
619    /**
620     * Make a copy of the array for the slice defined by array part from
621     * [begin, begin+length)
622     * @param values the input array
623     * @param begin start index of the array to include
624     * @param length number of elements to include from begin
625     * @return copy of a slice of the original array
626     */
627    private static double[] copyOf(final double[] values, final int begin, final int length) {
628        MathArrays.verifyValues(values, begin, length);
629        return Arrays.copyOfRange(values, begin, begin + length);
630    }
631
632    /**
633     * Replace every occurrence of a given value with a replacement value in a
634     * copied slice of array defined by array part from [begin, begin+length).
635     *
636     * @param values the input array
637     * @param begin start index of the array to include
638     * @param length number of elements to include from begin
639     * @param original the value to be replaced with
640     * @param replacement the value to be used for replacement
641     * @return the copy of sliced array with replaced values
642     */
643    private static double[] replaceAndSlice(final double[] values,
644                                            final int begin, final int length,
645                                            final double original,
646                                            final double replacement) {
647        final double[] temp = copyOf(values, begin, length);
648        for(int i = 0; i < length; i++) {
649            temp[i] = Precision.equalsIncludingNaN(original, temp[i]) ?
650                replacement :
651                temp[i];
652        }
653
654        return temp;
655    }
656    /**
657     * Remove the occurrence of a given value in a copied slice of array
658     * defined by the array part from [begin, begin+length).
659     * @param values the input array
660     * @param begin start index of the array to include
661     * @param length number of elements to include from begin
662     * @param removedValue the value to be removed from the sliced array
663     * @return the copy of the sliced array after removing the removedValue
664     */
665    private static double[] removeAndSlice(final double[] values,
666                                           final int begin, final int length,
667                                           final double removedValue) {
668        MathArrays.verifyValues(values, begin, length);
669        final double[] temp;
670        // Indicates where the removedValue is located
671        final BitSet bits = new BitSet(length);
672        for (int i = begin; i < begin+length; i++) {
673            if (Precision.equalsIncludingNaN(removedValue, values[i])) {
674                bits.set(i - begin);
675            }
676        }
677        // Check if empty then create a new copy
678        if (bits.isEmpty()) {
679            // Nothing removed, just copy
680            temp = copyOf(values, begin, length);
681        } else if(bits.cardinality() == length) {
682            // All removed, just empty
683            temp = new double[0];
684        } else {
685            // Some removable, so new
686            temp = new double[length - bits.cardinality()];
687            // Index from source array (i.e values)
688            int start = begin;
689            // Index in destination array(i.e temp)
690            int dest = 0;
691            // Index of bit set of next one
692            int nextOne = -1;
693            // Start index pointer of bitset
694            int bitSetPtr = 0;
695            while ((nextOne = bits.nextSetBit(bitSetPtr)) != -1) {
696                final int lengthToCopy = nextOne - bitSetPtr;
697                System.arraycopy(values, start, temp, dest, lengthToCopy);
698                dest += lengthToCopy;
699                start = begin + (bitSetPtr = bits.nextClearBit(nextOne));
700            }
701            // Copy any residue past start index till begin+length
702            if (start < begin + length) {
703                System.arraycopy(values,start,temp,dest,begin + length - start);
704            }
705        }
706        return temp;
707    }
708    /**
709     * Remove weights element if the corresponding data is equal to the given value
710     * in [begin, begin+length)
711     *
712     * @param values the input array
713     * @param sampleWeights weights of the input array
714     * @param begin start index of the array to include
715     * @param length number of elements to include from begin
716     * @param removedValue the value to be removed from the sliced array
717     * @return the copy of the sliced array after removing weights
718     */
719    private static double[] removeAndSliceByRef(final double[] values,
720                                               final double[] sampleWeights,
721                                               final int begin, final int length,
722                                               final double removedValue) {
723        MathArrays.verifyValues(values, begin, length);
724        final double[] temp;
725        //BitSet(length) to indicate where the removedValue is located
726        final BitSet bits = new BitSet(length);
727        for (int i = begin; i < begin+length; i++) {
728            if (Precision.equalsIncludingNaN(removedValue, values[i])) {
729                bits.set(i - begin);
730            }
731        }
732        //Check if empty then create a new copy
733        if (bits.isEmpty()) {
734            temp = copyOf(sampleWeights, begin, length); // Nothing removed, just copy
735        } else if(bits.cardinality() == length) {
736            temp = new double[0];                 // All removed, just empty
737        }else {                                   // Some removable, so new
738            temp = new double[length - bits.cardinality()];
739            int start = begin;  //start index from source array (i.e sampleWeights)
740            int dest = 0;       //dest index in destination array(i.e temp)
741            int nextOne = -1;   //nextOne is the index of bit set of next one
742            int bitSetPtr = 0;  //bitSetPtr is start index pointer of bitset
743            while ((nextOne = bits.nextSetBit(bitSetPtr)) != -1) {
744                final int lengthToCopy = nextOne - bitSetPtr;
745                System.arraycopy(sampleWeights, start, temp, dest, lengthToCopy);
746                dest += lengthToCopy;
747                start = begin + (bitSetPtr = bits.nextClearBit(nextOne));
748            }
749            //Copy any residue past start index till begin+length
750            if (start < begin + length) {
751                System.arraycopy(sampleWeights,start,temp,dest,begin + length - start);
752            }
753        }
754        return temp;
755    }
756    /**
757     * Get pivots which is either cached or a newly created one
758     *
759     * @param values array containing the input numbers
760     * @return cached pivots or a newly created one
761     */
762    private int[] getPivots(final double[] values) {
763        final int[] pivotsHeap;
764        if (values == getDataRef()) {
765            pivotsHeap = cachedPivots;
766        } else {
767            pivotsHeap = new int[PIVOTS_HEAP_LENGTH];
768            Arrays.fill(pivotsHeap, -1);
769        }
770        return pivotsHeap;
771    }
772
773    /**
774     * Get the estimation {@link EstimationType type} used for computation.
775     *
776     * @return the {@code estimationType} set
777     */
778    public EstimationType getEstimationType() {
779        return estimationType;
780    }
781
782    /**
783     * Build a new instance similar to the current one except for the
784     * {@link EstimationType estimation type}.
785     * <p>
786     * This method is intended to be used as part of a fluent-type builder
787     * pattern. Building finely tune instances should be done as follows:
788     * </p>
789     * <pre>
790     *   Percentile customized = new Percentile(quantile).
791     *                           withEstimationType(estimationType).
792     *                           withNaNStrategy(nanStrategy).
793     *                           withKthSelector(kthSelector);
794     * </pre>
795     * <p>
796     * If any of the {@code withXxx} method is omitted, the default value for
797     * the corresponding customization parameter will be used.
798     * </p>
799     * @param newEstimationType estimation type for the new instance.
800     * Cannot be {@code null}.
801     * @return a new instance, with changed estimation type
802     */
803    public Percentile withEstimationType(final EstimationType newEstimationType) {
804        return new Percentile(quantile, newEstimationType, nanStrategy, kthSelector);
805    }
806
807    /**
808     * Get the {@link NaNStrategy NaN Handling} strategy used for computation.
809     * @return {@code NaN Handling} strategy set during construction
810     */
811    public NaNStrategy getNaNStrategy() {
812        return nanStrategy;
813    }
814
815    /**
816     * Build a new instance similar to the current one except for the
817     * {@link NaNStrategy NaN handling} strategy.
818     * <p>
819     * This method is intended to be used as part of a fluent-type builder
820     * pattern. Building finely tune instances should be done as follows:
821     * </p>
822     * <pre>
823     *   Percentile customized = new Percentile(quantile).
824     *                           withEstimationType(estimationType).
825     *                           withNaNStrategy(nanStrategy).
826     *                           withKthSelector(kthSelector);
827     * </pre>
828     * <p>
829     * If any of the {@code withXxx} method is omitted, the default value for
830     * the corresponding customization parameter will be used.
831     * </p>
832     * @param newNaNStrategy NaN strategy for the new instance.
833     * Cannot be {@code null}.
834     * @return a new instance, with changed NaN handling strategy
835     */
836    public Percentile withNaNStrategy(final NaNStrategy newNaNStrategy) {
837        return new Percentile(quantile, estimationType, newNaNStrategy, kthSelector);
838    }
839
840    /**
841     * Get the {@link KthSelector kthSelector} used for computation.
842     * @return the {@code kthSelector} set
843     */
844    public KthSelector getKthSelector() {
845        return kthSelector;
846    }
847
848    /**
849     * Get the {@link PivotingStrategyInterface} used in KthSelector for computation.
850     * @return the pivoting strategy set
851     */
852    public PivotingStrategyInterface getPivotingStrategy() {
853        return kthSelector.getPivotingStrategy();
854    }
855
856    /**
857     * Build a new instance similar to the current one except for the
858     * {@link KthSelector kthSelector} instance specifically set.
859     * <p>
860     * This method is intended to be used as part of a fluent-type builder
861     * pattern. Building finely tune instances should be done as follows:
862     * </p>
863     * <pre>
864     *   Percentile customized = new Percentile(quantile).
865     *                           withEstimationType(estimationType).
866     *                           withNaNStrategy(nanStrategy).
867     *                           withKthSelector(newKthSelector);
868     * </pre>
869     * <p>
870     * If any of the {@code withXxx} method is omitted, the default value for
871     * the corresponding customization parameter will be used.
872     * </p>
873     * @param newKthSelector KthSelector for the new instance.
874     * Cannot be {@code null}.
875     * @return a new instance, with changed KthSelector
876     */
877    public Percentile withKthSelector(final KthSelector newKthSelector) {
878        return new Percentile(quantile, estimationType, nanStrategy, newKthSelector);
879    }
880
881    /**
882     * An enum for various estimation strategies of a percentile referred in
883     * <a href="http://en.wikipedia.org/wiki/Quantile">wikipedia on quantile</a>
884     * with the names of enum matching those of types mentioned in
885     * wikipedia.
886     * <p>
887     * Each enum corresponding to the specific type of estimation in wikipedia
888     * implements  the respective formulae that specializes in the below aspects
889     * <ul>
890     * <li>An <b>index method</b> to calculate approximate index of the
891     * estimate</li>
892     * <li>An <b>estimate method</b> to estimate a value found at the earlier
893     * computed index</li>
894     * <li>A <b> minLimit</b> on the quantile for which first element of sorted
895     * input is returned as an estimate </li>
896     * <li>A <b> maxLimit</b> on the quantile for which last element of sorted
897     * input is returned as an estimate </li>
898     * </ul>
899     * <p>
900     * Users can now create {@link Percentile} by explicitly passing this enum;
901     * such as by invoking {@link Percentile#withEstimationType(EstimationType)}
902     * <p>
903     * References:
904     * <ol>
905     * <li>
906     * <a href="http://en.wikipedia.org/wiki/Quantile">Wikipedia on quantile</a>
907     * </li>
908     * <li>
909     * <a href="https://www.amherst.edu/media/view/129116/.../Sample+Quantiles.pdf">
910     * Hyndman, R. J. and Fan, Y. (1996) Sample quantiles in statistical
911     * packages, American Statistician 50, 361–365</a> </li>
912     * <li>
913     * <a href="http://stat.ethz.ch/R-manual/R-devel/library/stats/html/quantile.html">
914     * R-Manual </a></li>
915     * </ol>
916     */
917    public enum EstimationType {
918        /**
919         * This is the default type used in the {@link Percentile}.This method
920         * has the following formulae for index and estimates<br>
921         * \( \begin{align}
922         * &amp;index    = (N+1)p\ \\
923         * &amp;estimate = x_{\lceil h\,-\,1/2 \rceil} \\
924         * &amp;minLimit = 0 \\
925         * &amp;maxLimit = 1 \\
926         * \end{align}\)
927         */
928        LEGACY("Legacy Apache Commons Math") {
929            /**
930             * {@inheritDoc}.This method in particular makes use of existing
931             * Apache Commons Math style of picking up the index.
932             */
933            @Override
934            protected double index(final double p, final int length) {
935                final double minLimit = 0d;
936                final double maxLimit = 1d;
937                return Double.compare(p, minLimit) == 0 ? 0 :
938                       Double.compare(p, maxLimit) == 0 ?
939                               length : p * (length + 1);
940            }
941            @Override
942            public double evaluate(final double[] work, final double[] sampleWeights,
943                                   final double p) {
944                throw new MathIllegalArgumentException(LocalizedFormats.UNSUPPORTED_OPERATION);
945            }
946        },
947        /**
948         * The method R_1 has the following formulae for index and estimates<br>
949         * \( \begin{align}
950         * &amp;index= Np + 1/2\,  \\
951         * &amp;estimate= x_{\lceil h\,-\,1/2 \rceil} \\
952         * &amp;minLimit = 0 \\
953         * \end{align}\)
954         */
955        R_1("R-1") {
956
957            @Override
958            protected double index(final double p, final int length) {
959                final double minLimit = 0d;
960                return Double.compare(p, minLimit) == 0 ? 0 : length * p + 0.5;
961            }
962
963            /**
964             * {@inheritDoc}This method in particular for R_1 uses ceil(pos-0.5)
965             */
966            @Override
967            protected double estimate(final double[] values,
968                                      final int[] pivotsHeap, final double pos,
969                                      final int length, final KthSelector selector) {
970                return super.estimate(values, pivotsHeap, FastMath.ceil(pos - 0.5), length, selector);
971            }
972            @Override
973            public double evaluate(final double[] work, final double[] sampleWeights,
974                                   final double p) {
975                throw new MathIllegalArgumentException(LocalizedFormats.UNSUPPORTED_OPERATION);
976            }
977        },
978        /**
979         * The method R_2 has the following formulae for index and estimates<br>
980         * \( \begin{align}
981         * &amp;index= Np + 1/2\, \\
982         * &amp;estimate=\frac{x_{\lceil h\,-\,1/2 \rceil} +
983         * x_{\lfloor h\,+\,1/2 \rfloor}}{2} \\
984         * &amp;minLimit = 0 \\
985         * &amp;maxLimit = 1 \\
986         * \end{align}\)
987         */
988        R_2("R-2") {
989
990            @Override
991            protected double index(final double p, final int length) {
992                final double minLimit = 0d;
993                final double maxLimit = 1d;
994                return Double.compare(p, maxLimit) == 0 ? length :
995                       Double.compare(p, minLimit) == 0 ? 0 : length * p + 0.5;
996            }
997
998            /**
999             * {@inheritDoc}This method in particular for R_2 averages the
1000             * values at ceil(p+0.5) and floor(p-0.5).
1001             */
1002            @Override
1003            protected double estimate(final double[] values,
1004                                      final int[] pivotsHeap, final double pos,
1005                                      final int length, final KthSelector selector) {
1006                final double low =
1007                        super.estimate(values, pivotsHeap, FastMath.ceil(pos - 0.5), length, selector);
1008                final double high =
1009                        super.estimate(values, pivotsHeap,FastMath.floor(pos + 0.5), length, selector);
1010                return (low + high) / 2;
1011            }
1012            @Override
1013            public double evaluate(final double[] work, final double[] sampleWeights,
1014                                   final double p) {
1015                throw new MathIllegalArgumentException(LocalizedFormats.UNSUPPORTED_OPERATION);
1016            }
1017        },
1018        /**
1019         * The method R_3 has the following formulae for index and estimates<br>
1020         * \( \begin{align}
1021         * &amp;index= Np \\
1022         * &amp;estimate= x_{\lfloor h \rceil}\, \\
1023         * &amp;minLimit = 0.5/N \\
1024         * \end{align}\)
1025         */
1026        R_3("R-3") {
1027            @Override
1028            protected double index(final double p, final int length) {
1029                final double minLimit = 1d/2 / length;
1030                return Double.compare(p, minLimit) <= 0 ?
1031                        0 : FastMath.rint(length * p);
1032            }
1033            @Override
1034            public double evaluate(final double[] work, final double[] sampleWeights,
1035                                   final double p) {
1036                throw new MathIllegalArgumentException(LocalizedFormats.UNSUPPORTED_OPERATION);
1037            }
1038        },
1039        /**
1040         * The method R_4 has the following formulae for index and estimates<br>
1041         * \( \begin{align}
1042         * &amp;index= Np\, \\
1043         * &amp;estimate= x_{\lfloor h \rfloor} + (h -
1044         * \lfloor h \rfloor) (x_{\lfloor h \rfloor + 1} - x_{\lfloor h
1045         * \rfloor}) \\
1046         * &amp;minLimit = 1/N \\
1047         * &amp;maxLimit = 1 \\
1048         * \end{align}\)
1049         */
1050        R_4("R-4") {
1051            @Override
1052            protected double index(final double p, final int length) {
1053                final double minLimit = 1d / length;
1054                final double maxLimit = 1d;
1055                return Double.compare(p, minLimit) < 0 ? 0 :
1056                       Double.compare(p, maxLimit) == 0 ? length : length * p;
1057            }
1058            @Override
1059            public double evaluate(final double[] work, final double[] sampleWeights,
1060                                   final double p) {
1061                throw new MathIllegalArgumentException(LocalizedFormats.UNSUPPORTED_OPERATION);
1062            }
1063        },
1064        /**
1065         * The method R_5 has the following formulae for index and estimates<br>
1066         * \( \begin{align}
1067         * &amp;index= Np + 1/2\\
1068         * &amp;estimate= x_{\lfloor h \rfloor} + (h -
1069         * \lfloor h \rfloor) (x_{\lfloor h \rfloor + 1} - x_{\lfloor h
1070         * \rfloor}) \\
1071         * &amp;minLimit = 0.5/N \\
1072         * &amp;maxLimit = (N-0.5)/N
1073         * \end{align}\)
1074         */
1075        R_5("R-5") {
1076
1077            @Override
1078            protected double index(final double p, final int length) {
1079                final double minLimit = 1d/2 / length;
1080                final double maxLimit = (length - 0.5) / length;
1081                return Double.compare(p, minLimit) < 0 ? 0 :
1082                       Double.compare(p, maxLimit) >= 0 ?
1083                               length : length * p + 0.5;
1084            }
1085            @Override
1086            public double evaluate(final double[] work, final double[] sampleWeights,
1087                                   final double p) {
1088                throw new MathIllegalArgumentException(LocalizedFormats.UNSUPPORTED_OPERATION);
1089            }
1090        },
1091        /**
1092         * The method R_6 has the following formulae for index and estimates<br>
1093         * \( \begin{align}
1094         * &amp;index= (N + 1)p \\
1095         * &amp;estimate= x_{\lfloor h \rfloor} + (h -
1096         * \lfloor h \rfloor) (x_{\lfloor h \rfloor + 1} - x_{\lfloor h
1097         * \rfloor}) \\
1098         * &amp;minLimit = 1/(N+1) \\
1099         * &amp;maxLimit = N/(N+1) \\
1100         * \end{align}\)
1101         * <p>
1102         * <b>Note:</b> This method computes the index in a manner very close to
1103         * the default Commons Math Percentile existing implementation. However
1104         * the difference to be noted is in picking up the limits with which
1105         * first element (p&lt;1(N+1)) and last elements (p&gt;N/(N+1))are done.
1106         * While in default case; these are done with p=0 and p=1 respectively.
1107         */
1108        R_6("R-6") {
1109
1110            @Override
1111            protected double index(final double p, final int length) {
1112                final double minLimit = 1d / (length + 1);
1113                final double maxLimit = 1d * length / (length + 1);
1114                return Double.compare(p, minLimit) < 0 ? 0 :
1115                       Double.compare(p, maxLimit) >= 0 ?
1116                               length : (length + 1) * p;
1117            }
1118            @Override
1119            public double evaluate(final double[] work, final double[] sampleWeights,
1120                                   final double p) {
1121                throw new MathIllegalArgumentException(LocalizedFormats.UNSUPPORTED_OPERATION);
1122            }
1123        },
1124
1125        /**
1126         * The method R_7 implements Microsoft Excel style computation has the
1127         * following formulae for index and estimates.<br>
1128         * \( \begin{align}
1129         * &amp;index = (N-1)p + 1 \\
1130         * &amp;estimate = x_{\lfloor h \rfloor} + (h -
1131         * \lfloor h \rfloor) (x_{\lfloor h \rfloor + 1} - x_{\lfloor h
1132         * \rfloor}) \\
1133         * &amp;minLimit = 0 \\
1134         * &amp;maxLimit = 1 \\
1135         * \end{align}\)
1136         * The formula to evaluate weighted percentiles is as following.<br>
1137         * \( \begin{align}
1138         * &amp;S_k = (k-1)w_k + (n-1)\sum_{i=1}^{k-1}w_i
1139         * &amp;Then find k s.t. \frac{S_k}{S_n}\leq p \leq \frac{S_{k+1}}{S_n}
1140         * \end{align}\)
1141         */
1142        R_7("R-7") {
1143            @Override
1144            protected double index(final double p, final int length) {
1145                final double minLimit = 0d;
1146                final double maxLimit = 1d;
1147                return Double.compare(p, minLimit) == 0 ? 0 :
1148                       Double.compare(p, maxLimit) == 0 ?
1149                               length : 1 + (length - 1) * p;
1150            }
1151
1152            @Override
1153            public double evaluate(final double[] work, final double[] sampleWeights,
1154                                   final double p) {
1155                MathArrays.sortInPlace(work, sampleWeights);
1156                double[] sk = new double[work.length];
1157                for(int k = 0; k < work.length; k++) {
1158                    sk[k] = 0;
1159                    for (int j = 0; j < k; j++) {
1160                        sk[k] += sampleWeights[j];
1161                    }
1162                    sk[k] = k * sampleWeights[k] + (work.length - 1) * sk[k];
1163                }
1164
1165                double qsn = (p / 100) * sk[sk.length-1];
1166                int k = searchSk(qsn, sk, 0, work.length - 1);
1167
1168                double ret;
1169                if (qsn == sk[k] && k == work.length - 1) {
1170                    ret = work[k] - (work[k] - work[k-1]) * (1 - (qsn - sk[k]) / (sk[k] - sk[k-1]));
1171                } else {
1172                    ret = work[k] + (work[k+1] - work[k]) * (qsn - sk[k]) / (sk[k+1] - sk[k]);
1173                }
1174                return ret;
1175            }
1176        },
1177
1178        /**
1179         * The method R_8 has the following formulae for index and estimates<br>
1180         * \( \begin{align}
1181         * &amp;index = (N + 1/3)p + 1/3  \\
1182         * &amp;estimate = x_{\lfloor h \rfloor} + (h -
1183           \lfloor h \rfloor) (x_{\lfloor h \rfloor + 1} - x_{\lfloor h
1184         * \rfloor}) \\
1185         * &amp;minLimit = (2/3)/(N+1/3) \\
1186         * &amp;maxLimit = (N-1/3)/(N+1/3) \\
1187         * \end{align}\)
1188         * <p>
1189         * As per Ref [2,3] this approach is most recommended as it provides
1190         * an approximate median-unbiased estimate regardless of distribution.
1191         */
1192        R_8("R-8") {
1193            @Override
1194            protected double index(final double p, final int length) {
1195                final double minLimit = 2 * (1d / 3) / (length + 1d / 3);
1196                final double maxLimit =
1197                        (length - 1d / 3) / (length + 1d / 3);
1198                return Double.compare(p, minLimit) < 0 ? 0 :
1199                       Double.compare(p, maxLimit) >= 0 ? length :
1200                           (length + 1d / 3) * p + 1d / 3;
1201            }
1202            @Override
1203            public double evaluate(final double[] work, final double[] sampleWeights,
1204                                   final double p) {
1205                throw new MathIllegalArgumentException(LocalizedFormats.UNSUPPORTED_OPERATION);
1206            }
1207        },
1208
1209        /**
1210         * The method R_9 has the following formulae for index and estimates<br>
1211         * \( \begin{align}
1212         * &amp;index = (N + 1/4)p + 3/8\\
1213         * &amp;estimate = x_{\lfloor h \rfloor} + (h -
1214           \lfloor h \rfloor) (x_{\lfloor h \rfloor + 1} - x_{\lfloor h
1215         * \rfloor}) \\
1216         * &amp;minLimit = (5/8)/(N+1/4) \\
1217         * &amp;maxLimit = (N-3/8)/(N+1/4) \\
1218         * \end{align}\)
1219         */
1220        R_9("R-9") {
1221            @Override
1222            protected double index(final double p, final int length) {
1223                final double minLimit = 5d/8 / (length + 0.25);
1224                final double maxLimit = (length - 3d/8) / (length + 0.25);
1225                return Double.compare(p, minLimit) < 0 ? 0 :
1226                       Double.compare(p, maxLimit) >= 0 ? length :
1227                               (length + 0.25) * p + 3d/8;
1228            }
1229            @Override
1230            public double evaluate(final double[] work, final double[] sampleWeights,
1231                                   final double p) {
1232                throw new MathIllegalArgumentException(LocalizedFormats.UNSUPPORTED_OPERATION);
1233            }
1234        },
1235        ;
1236
1237        /** Simple name such as R-1, R-2 corresponding to those in wikipedia. */
1238        private final String name;
1239
1240        /**
1241         * Constructor
1242         *
1243         * @param type name of estimation type as per wikipedia
1244         */
1245        EstimationType(final String type) {
1246            this.name = type;
1247        }
1248
1249        /**
1250         * Finds the index of array that can be used as starting index to
1251         * {@link #estimate(double[], int[], double, int, KthSelector) estimate}
1252         * percentile. The calculation of index calculation is specific to each
1253         * {@link EstimationType}.
1254         *
1255         * @param p the p<sup>th</sup> quantile
1256         * @param length the total number of array elements in the work array
1257         * @return a computed real valued index as explained in the wikipedia
1258         */
1259        protected abstract double index(final double p, final int length);
1260
1261        /**
1262         * Estimation based on K<sup>th</sup> selection. This may be overridden
1263         * in specific enums to compute slightly different estimations.
1264         *
1265         * @param work array of numbers to be used for finding the percentile
1266         * @param pos indicated positional index prior computed from calling
1267         *            {@link #index(double, int)}
1268         * @param pivotsHeap an earlier populated cache if exists; will be used
1269         * @param length size of array considered
1270         * @param selector a {@link KthSelector} used for pivoting during search
1271         * @return estimated percentile
1272         */
1273        protected double estimate(final double[] work, final int[] pivotsHeap,
1274                                  final double pos, final int length,
1275                                  final KthSelector selector) {
1276
1277            final double fpos = FastMath.floor(pos);
1278            final int intPos = (int) fpos;
1279            final double dif = pos - fpos;
1280
1281            if (pos < 1) {
1282                return selector.select(work, pivotsHeap, 0);
1283            }
1284            if (pos >= length) {
1285                return selector.select(work, pivotsHeap, length - 1);
1286            }
1287
1288            final double lower = selector.select(work, pivotsHeap, intPos - 1);
1289            final double upper = selector.select(work, pivotsHeap, intPos);
1290            return lower + dif * (upper - lower);
1291        }
1292
1293        /**
1294         * Evaluate method to compute the percentile for a given bounded array
1295         * using earlier computed pivots heap.<br>
1296         * This basically calls the {@link #index(double, int) index} and then
1297         * {@link #estimate(double[], int[], double, int, KthSelector) estimate}
1298         * functions to return the estimated percentile value.
1299         *
1300         * @param work array of numbers to be used for finding the percentile.
1301         * Cannot be {@code null}.
1302         * @param pivotsHeap a prior cached heap which can speed up estimation
1303         * @param p the p<sup>th</sup> quantile to be computed
1304         * @param selector a {@link KthSelector} used for pivoting during search
1305         * @return estimated percentile
1306         * @throws OutOfRangeException if p is out of range
1307         */
1308        protected double evaluate(final double[] work, final int[] pivotsHeap, final double p,
1309                                  final KthSelector selector) {
1310            MathUtils.checkNotNull(work);
1311            if (p > 100 || p <= 0) {
1312                throw new OutOfRangeException(LocalizedFormats.OUT_OF_BOUNDS_QUANTILE_VALUE,
1313                                              p, 0, 100);
1314            }
1315            return estimate(work, pivotsHeap, index(p/100d, work.length), work.length, selector);
1316        }
1317
1318        /**
1319         * Evaluate method to compute the percentile for a given bounded array.
1320         * This basically calls the {@link #index(double, int) index} and then
1321         * {@link #estimate(double[], int[], double, int, KthSelector) estimate}
1322         * functions to return the estimated percentile value. Please
1323         * note that this method does not make use of cached pivots.
1324         *
1325         * @param work array of numbers to be used for finding the percentile.
1326         * Cannot be {@code null}.
1327         * @param p the p<sup>th</sup> quantile to be computed
1328         * @return estimated percentile
1329         * @param selector a {@link KthSelector} used for pivoting during search
1330         * @throws OutOfRangeException if length or p is out of range
1331         */
1332        public double evaluate(final double[] work, final double p, final KthSelector selector) {
1333            return this.evaluate(work, null, p, selector);
1334        }
1335        /**
1336         * Evaluate weighted percentile by estimation rule specified in {@link EstimationType}.
1337         * @param work array of numbers to be used for finding the percentile
1338         * @param sampleWeights the corresponding weights of data in work
1339         * @param p the p<sup>th</sup> quantile to be computed
1340         * @return estimated weighted percentile
1341         * @throws MathIllegalArgumentException if weighted percentile is not supported by the current estimationType
1342         */
1343        public abstract double evaluate(final double[] work, final double[] sampleWeights,
1344                                        final double p);
1345        /**
1346         * Search the interval q*sn locates in.
1347         * @param qsn q*sn, where n refers to the data size
1348         * @param sk the cumulative weights array
1349         * @param lo start position to search qsn
1350         * @param hi end position to search qsn
1351         * @return the index of lower bound qsn locates in
1352         */
1353        private static int searchSk(double qsn, double[] sk, int lo, int hi) {
1354            if (sk.length == 1) {
1355                return 0;
1356            }
1357            if (hi - lo == 1) {
1358                if (qsn == sk[hi]) {
1359                    return hi;
1360                }
1361                return lo;
1362            } else {
1363                int mid = (lo + hi) / 2;
1364                if (qsn == sk[mid]) {
1365                  return mid;
1366                } else if (qsn > sk[mid]) {
1367                    return searchSk(qsn, sk, mid, hi);
1368                } else {
1369                    return searchSk(qsn, sk, lo, mid);
1370                }
1371            }
1372        }
1373        /**
1374         * Gets the name of the enum
1375         *
1376         * @return the name
1377         */
1378        String getName() {
1379            return name;
1380        }
1381    }
1382}