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     */
017    
018    package org.apache.commons.math3.random;
019    
020    import java.io.BufferedReader;
021    import java.io.File;
022    import java.io.FileInputStream;
023    import java.io.IOException;
024    import java.io.InputStream;
025    import java.io.InputStreamReader;
026    import java.net.URL;
027    import java.nio.charset.Charset;
028    import java.util.ArrayList;
029    import java.util.List;
030    
031    import org.apache.commons.math3.distribution.AbstractRealDistribution;
032    import org.apache.commons.math3.distribution.NormalDistribution;
033    import org.apache.commons.math3.distribution.RealDistribution;
034    import org.apache.commons.math3.exception.MathIllegalStateException;
035    import org.apache.commons.math3.exception.MathInternalError;
036    import org.apache.commons.math3.exception.NullArgumentException;
037    import org.apache.commons.math3.exception.OutOfRangeException;
038    import org.apache.commons.math3.exception.ZeroException;
039    import org.apache.commons.math3.exception.util.LocalizedFormats;
040    import org.apache.commons.math3.stat.descriptive.StatisticalSummary;
041    import org.apache.commons.math3.stat.descriptive.SummaryStatistics;
042    import org.apache.commons.math3.util.FastMath;
043    import org.apache.commons.math3.util.MathUtils;
044    
045    /**
046     * <p>Represents an <a href="http://http://en.wikipedia.org/wiki/Empirical_distribution_function">
047     * empirical probability distribution</a> -- a probability distribution derived
048     * from observed data without making any assumptions about the functional form
049     * of the population distribution that the data come from.</p>
050     *
051     * <p>An <code>EmpiricalDistribution</code> maintains data structures, called
052     * <i>distribution digests</i>, that describe empirical distributions and
053     * support the following operations: <ul>
054     * <li>loading the distribution from a file of observed data values</li>
055     * <li>dividing the input data into "bin ranges" and reporting bin frequency
056     *     counts (data for histogram)</li>
057     * <li>reporting univariate statistics describing the full set of data values
058     *     as well as the observations within each bin</li>
059     * <li>generating random values from the distribution</li>
060     * </ul>
061     * Applications can use <code>EmpiricalDistribution</code> to build grouped
062     * frequency histograms representing the input data or to generate random values
063     * "like" those in the input file -- i.e., the values generated will follow the
064     * distribution of the values in the file.</p>
065     *
066     * <p>The implementation uses what amounts to the
067     * <a href="http://nedwww.ipac.caltech.edu/level5/March02/Silverman/Silver2_6.html">
068     * Variable Kernel Method</a> with Gaussian smoothing:<p>
069     * <strong>Digesting the input file</strong>
070     * <ol><li>Pass the file once to compute min and max.</li>
071     * <li>Divide the range from min-max into <code>binCount</code> "bins."</li>
072     * <li>Pass the data file again, computing bin counts and univariate
073     *     statistics (mean, std dev.) for each of the bins </li>
074     * <li>Divide the interval (0,1) into subintervals associated with the bins,
075     *     with the length of a bin's subinterval proportional to its count.</li></ol>
076     * <strong>Generating random values from the distribution</strong><ol>
077     * <li>Generate a uniformly distributed value in (0,1) </li>
078     * <li>Select the subinterval to which the value belongs.
079     * <li>Generate a random Gaussian value with mean = mean of the associated
080     *     bin and std dev = std dev of associated bin.</li></ol></p>
081     *
082     * <p>EmpiricalDistribution implements the {@link RealDistribution} interface
083     * as follows.  Given x within the range of values in the dataset, let B
084     * be the bin containing x and let K be the within-bin kernel for B.  Let P(B-)
085     * be the sum of the probabilities of the bins below B and let K(B) be the
086     * mass of B under K (i.e., the integral of the kernel density over B).  Then
087     * set P(X < x) = P(B-) + P(B) * K(x) / K(B) where K(x) is the kernel distribution
088     * evaluated at x. This results in a cdf that matches the grouped frequency
089     * distribution at the bin endpoints and interpolates within bins using
090     * within-bin kernels.</p>
091     *
092     *<strong>USAGE NOTES:</strong><ul>
093     *<li>The <code>binCount</code> is set by default to 1000.  A good rule of thumb
094     *    is to set the bin count to approximately the length of the input file divided
095     *    by 10. </li>
096     *<li>The input file <i>must</i> be a plain text file containing one valid numeric
097     *    entry per line.</li>
098     * </ul></p>
099     *
100     * @version $Id: EmpiricalDistribution.java 1422350 2012-12-15 20:47:47Z psteitz $
101     */
102    public class EmpiricalDistribution extends AbstractRealDistribution {
103    
104        /** Default bin count */
105        public static final int DEFAULT_BIN_COUNT = 1000;
106    
107        /** Character set for file input */
108        private static final String FILE_CHARSET = "US-ASCII";
109    
110        /** Serializable version identifier */
111        private static final long serialVersionUID = 5729073523949762654L;
112    
113        /** List of SummaryStatistics objects characterizing the bins */
114        private final List<SummaryStatistics> binStats;
115    
116        /** Sample statistics */
117        private SummaryStatistics sampleStats = null;
118    
119        /** Max loaded value */
120        private double max = Double.NEGATIVE_INFINITY;
121    
122        /** Min loaded value */
123        private double min = Double.POSITIVE_INFINITY;
124    
125        /** Grid size */
126        private double delta = 0d;
127    
128        /** number of bins */
129        private final int binCount;
130    
131        /** is the distribution loaded? */
132        private boolean loaded = false;
133    
134        /** upper bounds of subintervals in (0,1) "belonging" to the bins */
135        private double[] upperBounds = null;
136    
137        /** RandomDataGenerator instance to use in repeated calls to getNext() */
138        private final RandomDataGenerator randomData;
139    
140        /**
141         * Creates a new EmpiricalDistribution with the default bin count.
142         */
143        public EmpiricalDistribution() {
144            this(DEFAULT_BIN_COUNT);
145        }
146    
147        /**
148         * Creates a new EmpiricalDistribution with the specified bin count.
149         *
150         * @param binCount number of bins
151         */
152        public EmpiricalDistribution(int binCount) {
153            this(binCount, new RandomDataGenerator());
154        }
155    
156        /**
157         * Creates a new EmpiricalDistribution with the specified bin count using the
158         * provided {@link RandomGenerator} as the source of random data.
159         *
160         * @param binCount number of bins
161         * @param generator random data generator (may be null, resulting in default JDK generator)
162         * @since 3.0
163         */
164        public EmpiricalDistribution(int binCount, RandomGenerator generator) {
165            this(binCount, new RandomDataGenerator(generator));
166        }
167    
168        /**
169         * Creates a new EmpiricalDistribution with default bin count using the
170         * provided {@link RandomGenerator} as the source of random data.
171         *
172         * @param generator random data generator (may be null, resulting in default JDK generator)
173         * @since 3.0
174         */
175        public EmpiricalDistribution(RandomGenerator generator) {
176            this(DEFAULT_BIN_COUNT, generator);
177        }
178    
179        /**
180         * Creates a new EmpiricalDistribution with the specified bin count using the
181         * provided {@link RandomDataImpl} instance as the source of random data.
182         *
183         * @param binCount number of bins
184         * @param randomData random data generator (may be null, resulting in default JDK generator)
185         * @since 3.0
186         * @deprecated As of 3.1. Please use {@link #EmpiricalDistribution(int,RandomGenerator)} instead.
187         */
188        @Deprecated
189        public EmpiricalDistribution(int binCount, RandomDataImpl randomData) {
190            this(binCount, randomData.getDelegate());
191        }
192    
193        /**
194         * Creates a new EmpiricalDistribution with default bin count using the
195         * provided {@link RandomDataImpl} as the source of random data.
196         *
197         * @param randomData random data generator (may be null, resulting in default JDK generator)
198         * @since 3.0
199         * @deprecated As of 3.1. Please use {@link #EmpiricalDistribution(RandomGenerator)} instead.
200         */
201        @Deprecated
202        public EmpiricalDistribution(RandomDataImpl randomData) {
203            this(DEFAULT_BIN_COUNT, randomData);
204        }
205    
206        /**
207         * Private constructor to allow lazy initialisation of the RNG contained
208         * in the {@link #randomData} instance variable.
209         *
210         * @param binCount number of bins
211         * @param randomData Random data generator.
212         */
213        private EmpiricalDistribution(int binCount,
214                                      RandomDataGenerator randomData) {
215            super(null);
216            this.binCount = binCount;
217            this.randomData = randomData;
218            binStats = new ArrayList<SummaryStatistics>();
219        }
220    
221        /**
222         * Computes the empirical distribution from the provided
223         * array of numbers.
224         *
225         * @param in the input data array
226         * @exception NullArgumentException if in is null
227         */
228        public void load(double[] in) throws NullArgumentException {
229            DataAdapter da = new ArrayDataAdapter(in);
230            try {
231                da.computeStats();
232                // new adapter for the second pass
233                fillBinStats(new ArrayDataAdapter(in));
234            } catch (IOException ex) {
235                // Can't happen
236                throw new MathInternalError();
237            }
238            loaded = true;
239    
240        }
241    
242        /**
243         * Computes the empirical distribution using data read from a URL.
244         *
245         * <p>The input file <i>must</i> be an ASCII text file containing one
246         * valid numeric entry per line.</p>
247         *
248         * @param url url of the input file
249         *
250         * @throws IOException if an IO error occurs
251         * @throws NullArgumentException if url is null
252         * @throws ZeroException if URL contains no data
253         */
254        public void load(URL url) throws IOException, NullArgumentException, ZeroException {
255            MathUtils.checkNotNull(url);
256            Charset charset = Charset.forName(FILE_CHARSET);
257            BufferedReader in =
258                new BufferedReader(new InputStreamReader(url.openStream(), charset));
259            try {
260                DataAdapter da = new StreamDataAdapter(in);
261                da.computeStats();
262                if (sampleStats.getN() == 0) {
263                    throw new ZeroException(LocalizedFormats.URL_CONTAINS_NO_DATA, url);
264                }
265                // new adapter for the second pass
266                in = new BufferedReader(new InputStreamReader(url.openStream(), charset));
267                fillBinStats(new StreamDataAdapter(in));
268                loaded = true;
269            } finally {
270               try {
271                   in.close();
272               } catch (IOException ex) { //NOPMD
273                   // ignore
274               }
275            }
276        }
277    
278        /**
279         * Computes the empirical distribution from the input file.
280         *
281         * <p>The input file <i>must</i> be an ASCII text file containing one
282         * valid numeric entry per line.</p>
283         *
284         * @param file the input file
285         * @throws IOException if an IO error occurs
286         * @throws NullArgumentException if file is null
287         */
288        public void load(File file) throws IOException, NullArgumentException {
289            MathUtils.checkNotNull(file);
290            Charset charset = Charset.forName(FILE_CHARSET);
291            InputStream is = new FileInputStream(file);
292            BufferedReader in = new BufferedReader(new InputStreamReader(is, charset));
293            try {
294                DataAdapter da = new StreamDataAdapter(in);
295                da.computeStats();
296                // new adapter for second pass
297                is = new FileInputStream(file);
298                in = new BufferedReader(new InputStreamReader(is, charset));
299                fillBinStats(new StreamDataAdapter(in));
300                loaded = true;
301            } finally {
302                try {
303                    in.close();
304                } catch (IOException ex) { //NOPMD
305                    // ignore
306                }
307            }
308        }
309    
310        /**
311         * Provides methods for computing <code>sampleStats</code> and
312         * <code>beanStats</code> abstracting the source of data.
313         */
314        private abstract class DataAdapter{
315    
316            /**
317             * Compute bin stats.
318             *
319             * @throws IOException  if an error occurs computing bin stats
320             */
321            public abstract void computeBinStats() throws IOException;
322    
323            /**
324             * Compute sample statistics.
325             *
326             * @throws IOException if an error occurs computing sample stats
327             */
328            public abstract void computeStats() throws IOException;
329    
330        }
331    
332        /**
333         * <code>DataAdapter</code> for data provided through some input stream
334         */
335        private class StreamDataAdapter extends DataAdapter{
336    
337            /** Input stream providing access to the data */
338            private BufferedReader inputStream;
339    
340            /**
341             * Create a StreamDataAdapter from a BufferedReader
342             *
343             * @param in BufferedReader input stream
344             */
345            public StreamDataAdapter(BufferedReader in){
346                super();
347                inputStream = in;
348            }
349    
350            /** {@inheritDoc} */
351            @Override
352            public void computeBinStats() throws IOException {
353                String str = null;
354                double val = 0.0d;
355                while ((str = inputStream.readLine()) != null) {
356                    val = Double.parseDouble(str);
357                    SummaryStatistics stats = binStats.get(findBin(val));
358                    stats.addValue(val);
359                }
360    
361                inputStream.close();
362                inputStream = null;
363            }
364    
365            /** {@inheritDoc} */
366            @Override
367            public void computeStats() throws IOException {
368                String str = null;
369                double val = 0.0;
370                sampleStats = new SummaryStatistics();
371                while ((str = inputStream.readLine()) != null) {
372                    val = Double.valueOf(str).doubleValue();
373                    sampleStats.addValue(val);
374                }
375                inputStream.close();
376                inputStream = null;
377            }
378        }
379    
380        /**
381         * <code>DataAdapter</code> for data provided as array of doubles.
382         */
383        private class ArrayDataAdapter extends DataAdapter {
384    
385            /** Array of input  data values */
386            private double[] inputArray;
387    
388            /**
389             * Construct an ArrayDataAdapter from a double[] array
390             *
391             * @param in double[] array holding the data
392             * @throws NullArgumentException if in is null
393             */
394            public ArrayDataAdapter(double[] in) throws NullArgumentException {
395                super();
396                MathUtils.checkNotNull(in);
397                inputArray = in;
398            }
399    
400            /** {@inheritDoc} */
401            @Override
402            public void computeStats() throws IOException {
403                sampleStats = new SummaryStatistics();
404                for (int i = 0; i < inputArray.length; i++) {
405                    sampleStats.addValue(inputArray[i]);
406                }
407            }
408    
409            /** {@inheritDoc} */
410            @Override
411            public void computeBinStats() throws IOException {
412                for (int i = 0; i < inputArray.length; i++) {
413                    SummaryStatistics stats =
414                        binStats.get(findBin(inputArray[i]));
415                    stats.addValue(inputArray[i]);
416                }
417            }
418        }
419    
420        /**
421         * Fills binStats array (second pass through data file).
422         *
423         * @param da object providing access to the data
424         * @throws IOException  if an IO error occurs
425         */
426        private void fillBinStats(final DataAdapter da)
427            throws IOException {
428            // Set up grid
429            min = sampleStats.getMin();
430            max = sampleStats.getMax();
431            delta = (max - min)/(Double.valueOf(binCount)).doubleValue();
432    
433            // Initialize binStats ArrayList
434            if (!binStats.isEmpty()) {
435                binStats.clear();
436            }
437            for (int i = 0; i < binCount; i++) {
438                SummaryStatistics stats = new SummaryStatistics();
439                binStats.add(i,stats);
440            }
441    
442            // Filling data in binStats Array
443            da.computeBinStats();
444    
445            // Assign upperBounds based on bin counts
446            upperBounds = new double[binCount];
447            upperBounds[0] =
448            ((double) binStats.get(0).getN()) / (double) sampleStats.getN();
449            for (int i = 1; i < binCount-1; i++) {
450                upperBounds[i] = upperBounds[i-1] +
451                ((double) binStats.get(i).getN()) / (double) sampleStats.getN();
452            }
453            upperBounds[binCount-1] = 1.0d;
454        }
455    
456        /**
457         * Returns the index of the bin to which the given value belongs
458         *
459         * @param value  the value whose bin we are trying to find
460         * @return the index of the bin containing the value
461         */
462        private int findBin(double value) {
463            return FastMath.min(
464                    FastMath.max((int) FastMath.ceil((value - min) / delta) - 1, 0),
465                    binCount - 1);
466        }
467    
468        /**
469         * Generates a random value from this distribution.
470         * <strong>Preconditions:</strong><ul>
471         * <li>the distribution must be loaded before invoking this method</li></ul>
472         * @return the random value.
473         * @throws MathIllegalStateException if the distribution has not been loaded
474         */
475        public double getNextValue() throws MathIllegalStateException {
476    
477            if (!loaded) {
478                throw new MathIllegalStateException(LocalizedFormats.DISTRIBUTION_NOT_LOADED);
479            }
480    
481            // Start with a uniformly distributed random number in (0,1)
482            final double x = randomData.nextUniform(0,1);
483    
484            // Use this to select the bin and generate a Gaussian within the bin
485            for (int i = 0; i < binCount; i++) {
486               if (x <= upperBounds[i]) {
487                   SummaryStatistics stats = binStats.get(i);
488                   if (stats.getN() > 0) {
489                       if (stats.getStandardDeviation() > 0) {  // more than one obs
490                           return randomData.nextGaussian(stats.getMean(),
491                                                          stats.getStandardDeviation());
492                       } else {
493                           return stats.getMean(); // only one obs in bin
494                       }
495                   }
496               }
497            }
498            throw new MathIllegalStateException(LocalizedFormats.NO_BIN_SELECTED);
499        }
500    
501        /**
502         * Returns a {@link StatisticalSummary} describing this distribution.
503         * <strong>Preconditions:</strong><ul>
504         * <li>the distribution must be loaded before invoking this method</li></ul>
505         *
506         * @return the sample statistics
507         * @throws IllegalStateException if the distribution has not been loaded
508         */
509        public StatisticalSummary getSampleStats() {
510            return sampleStats;
511        }
512    
513        /**
514         * Returns the number of bins.
515         *
516         * @return the number of bins.
517         */
518        public int getBinCount() {
519            return binCount;
520        }
521    
522        /**
523         * Returns a List of {@link SummaryStatistics} instances containing
524         * statistics describing the values in each of the bins.  The list is
525         * indexed on the bin number.
526         *
527         * @return List of bin statistics.
528         */
529        public List<SummaryStatistics> getBinStats() {
530            return binStats;
531        }
532    
533        /**
534         * <p>Returns a fresh copy of the array of upper bounds for the bins.
535         * Bins are: <br/>
536         * [min,upperBounds[0]],(upperBounds[0],upperBounds[1]],...,
537         *  (upperBounds[binCount-2], upperBounds[binCount-1] = max].</p>
538         *
539         * <p>Note: In versions 1.0-2.0 of commons-math, this method
540         * incorrectly returned the array of probability generator upper
541         * bounds now returned by {@link #getGeneratorUpperBounds()}.</p>
542         *
543         * @return array of bin upper bounds
544         * @since 2.1
545         */
546        public double[] getUpperBounds() {
547            double[] binUpperBounds = new double[binCount];
548            for (int i = 0; i < binCount - 1; i++) {
549                binUpperBounds[i] = min + delta * (i + 1);
550            }
551            binUpperBounds[binCount - 1] = max;
552            return binUpperBounds;
553        }
554    
555        /**
556         * <p>Returns a fresh copy of the array of upper bounds of the subintervals
557         * of [0,1] used in generating data from the empirical distribution.
558         * Subintervals correspond to bins with lengths proportional to bin counts.</p>
559         *
560         * <p>In versions 1.0-2.0 of commons-math, this array was (incorrectly) returned
561         * by {@link #getUpperBounds()}.</p>
562         *
563         * @since 2.1
564         * @return array of upper bounds of subintervals used in data generation
565         */
566        public double[] getGeneratorUpperBounds() {
567            int len = upperBounds.length;
568            double[] out = new double[len];
569            System.arraycopy(upperBounds, 0, out, 0, len);
570            return out;
571        }
572    
573        /**
574         * Property indicating whether or not the distribution has been loaded.
575         *
576         * @return true if the distribution has been loaded
577         */
578        public boolean isLoaded() {
579            return loaded;
580        }
581    
582        /**
583         * Reseeds the random number generator used by {@link #getNextValue()}.
584         *
585         * @param seed random generator seed
586         * @since 3.0
587         */
588        public void reSeed(long seed) {
589            randomData.reSeed(seed);
590        }
591    
592        // Distribution methods ---------------------------
593    
594        /**
595         * {@inheritDoc}
596         * @since 3.1
597         */
598        public double probability(double x) {
599            return 0;
600        }
601    
602        /**
603         * {@inheritDoc}
604         *
605         * <p>Returns the kernel density normalized so that its integral over each bin
606         * equals the bin mass.</p>
607         *
608         * <p>Algorithm description: <ol>
609         * <li>Find the bin B that x belongs to.</li>
610         * <li>Compute K(B) = the mass of B with respect to the within-bin kernel (i.e., the
611         * integral of the kernel density over B).</li>
612         * <li>Return k(x) * P(B) / K(B), where k is the within-bin kernel density
613         * and P(B) is the mass of B.</li></ol></p>
614         * @since 3.1
615         */
616        public double density(double x) {
617            if (x < min || x > max) {
618                return 0d;
619            }
620            final int binIndex = findBin(x);
621            final RealDistribution kernel = getKernel(binStats.get(binIndex));
622            return kernel.density(x) * pB(binIndex) / kB(binIndex);
623        }
624    
625        /**
626         * {@inheritDoc}
627         *
628         * <p>Algorithm description:<ol>
629         * <li>Find the bin B that x belongs to.</li>
630         * <li>Compute P(B) = the mass of B and P(B-) = the combined mass of the bins below B.</li>
631         * <li>Compute K(B) = the probability mass of B with respect to the within-bin kernel
632         * and K(B-) = the kernel distribution evaluated at the lower endpoint of B</li>
633         * <li>Return P(B-) + P(B) * [K(x) - K(B-)] / K(B) where
634         * K(x) is the within-bin kernel distribution function evaluated at x.</li></ol></p>
635         *
636         * @since 3.1
637         */
638        public double cumulativeProbability(double x) {
639            if (x < min) {
640                return 0d;
641            } else if (x >= max) {
642                return 1d;
643            }
644            final int binIndex = findBin(x);
645            final double pBminus = pBminus(binIndex);
646            final double pB = pB(binIndex);
647            final double[] binBounds = getUpperBounds();
648            final double kB = kB(binIndex);
649            final double lower = binIndex == 0 ? min : binBounds[binIndex - 1];
650            final RealDistribution kernel = k(x);
651            final double withinBinCum =
652                (kernel.cumulativeProbability(x) -  kernel.cumulativeProbability(lower)) / kB;
653            return pBminus + pB * withinBinCum;
654        }
655    
656        /**
657         * {@inheritDoc}
658         *
659         * <p>Algorithm description:<ol>
660         * <li>Find the smallest i such that the sum of the masses of the bins
661         *  through i is at least p.</li>
662         * <li>
663         *   Let K be the within-bin kernel distribution for bin i.</br>
664         *   Let K(B) be the mass of B under K. <br/>
665         *   Let K(B-) be K evaluated at the lower endpoint of B (the combined
666         *   mass of the bins below B under K).<br/>
667         *   Let P(B) be the probability of bin i.<br/>
668         *   Let P(B-) be the sum of the bin masses below bin i. <br/>
669         *   Let pCrit = p - P(B-)<br/>
670         * <li>Return the inverse of K evaluated at <br/>
671         *    K(B-) + pCrit * K(B) / P(B) </li>
672         *  </ol></p>
673         *
674         * @since 3.1
675         */
676        public double inverseCumulativeProbability(final double p) throws OutOfRangeException {
677            if (p < 0.0 || p > 1.0) {
678                throw new OutOfRangeException(p, 0, 1);
679            }
680    
681            if (p == 0.0) {
682                return getSupportLowerBound();
683            }
684    
685            if (p == 1.0) {
686                return getSupportUpperBound();
687            }
688    
689            int i = 0;
690            while (cumBinP(i) < p) {
691                i++;
692            }
693    
694            final RealDistribution kernel = getKernel(binStats.get(i));
695            final double kB = kB(i);
696            final double[] binBounds = getUpperBounds();
697            final double lower = i == 0 ? min : binBounds[i - 1];
698            final double kBminus = kernel.cumulativeProbability(lower);
699            final double pB = pB(i);
700            final double pBminus = pBminus(i);
701            final double pCrit = p - pBminus;
702            if (pCrit <= 0) {
703                return lower;
704            }
705            return kernel.inverseCumulativeProbability(kBminus + pCrit * kB / pB);
706        }
707    
708        /**
709         * {@inheritDoc}
710         * @since 3.1
711         */
712        public double getNumericalMean() {
713           return sampleStats.getMean();
714        }
715    
716        /**
717         * {@inheritDoc}
718         * @since 3.1
719         */
720        public double getNumericalVariance() {
721            return sampleStats.getVariance();
722        }
723    
724        /**
725         * {@inheritDoc}
726         * @since 3.1
727         */
728        public double getSupportLowerBound() {
729           return min;
730        }
731    
732        /**
733         * {@inheritDoc}
734         * @since 3.1
735         */
736        public double getSupportUpperBound() {
737            return max;
738        }
739    
740        /**
741         * {@inheritDoc}
742         * @since 3.1
743         */
744        public boolean isSupportLowerBoundInclusive() {
745            return true;
746        }
747    
748        /**
749         * {@inheritDoc}
750         * @since 3.1
751         */
752        public boolean isSupportUpperBoundInclusive() {
753            return true;
754        }
755    
756        /**
757         * {@inheritDoc}
758         * @since 3.1
759         */
760        public boolean isSupportConnected() {
761            return true;
762        }
763    
764        /**
765         * {@inheritDoc}
766         * @since 3.1
767         */
768        @Override
769        public double sample() {
770            return getNextValue();
771        }
772    
773        /**
774         * {@inheritDoc}
775         * @since 3.1
776         */
777        @Override
778        public void reseedRandomGenerator(long seed) {
779            randomData.reSeed(seed);
780        }
781    
782        /**
783         * The probability of bin i.
784         *
785         * @param i the index of the bin
786         * @return the probability that selection begins in bin i
787         */
788        private double pB(int i) {
789            return i == 0 ? upperBounds[0] :
790                upperBounds[i] - upperBounds[i - 1];
791        }
792    
793        /**
794         * The combined probability of the bins up to but not including bin i.
795         *
796         * @param i the index of the bin
797         * @return the probability that selection begins in a bin below bin i.
798         */
799        private double pBminus(int i) {
800            return i == 0 ? 0 : upperBounds[i - 1];
801        }
802    
803        /**
804         * Mass of bin i under the within-bin kernel of the bin.
805         *
806         * @param i index of the bin
807         * @return the difference in the within-bin kernel cdf between the
808         * upper and lower endpoints of bin i
809         */
810        @SuppressWarnings("deprecation")
811        private double kB(int i) {
812            final double[] binBounds = getUpperBounds();
813            final RealDistribution kernel = getKernel(binStats.get(i));
814            return i == 0 ? kernel.cumulativeProbability(min, binBounds[0]) :
815                kernel.cumulativeProbability(binBounds[i - 1], binBounds[i]);
816        }
817    
818        /**
819         * The within-bin kernel of the bin that x belongs to.
820         *
821         * @param x the value to locate within a bin
822         * @return the within-bin kernel of the bin containing x
823         */
824        private RealDistribution k(double x) {
825            final int binIndex = findBin(x);
826            return getKernel(binStats.get(binIndex));
827        }
828    
829        /**
830         * The combined probability of the bins up to and including binIndex.
831         *
832         * @param binIndex maximum bin index
833         * @return sum of the probabilities of bins through binIndex
834         */
835        private double cumBinP(int binIndex) {
836            return upperBounds[binIndex];
837        }
838    
839        /**
840         * The within-bin smoothing kernel.
841         *
842         * @param bStats summary statistics for the bin
843         * @return within-bin kernel parameterized by bStats
844         */
845        private RealDistribution getKernel(SummaryStatistics bStats) {
846            // For now, hard-code Gaussian (only kernel supported)
847            return new NormalDistribution(
848                    bStats.getMean(), bStats.getStandardDeviation());
849        }
850    }