AbstractRealDistribution.java

  1. /*
  2.  * Licensed to the Apache Software Foundation (ASF) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * The ASF licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *      http://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */
  17. package org.apache.commons.math4.legacy.distribution;

  19. import org.apache.commons.statistics.distribution.ContinuousDistribution;
  20. import org.apache.commons.math4.legacy.analysis.UnivariateFunction;
  21. import org.apache.commons.math4.legacy.analysis.solvers.UnivariateSolverUtils;
  22. import org.apache.commons.math4.legacy.exception.NumberIsTooLargeException;
  23. import org.apache.commons.math4.legacy.exception.OutOfRangeException;
  24. import org.apache.commons.math4.legacy.exception.util.LocalizedFormats;
  25. import org.apache.commons.rng.UniformRandomProvider;
  26. import org.apache.commons.rng.sampling.distribution.InverseTransformContinuousSampler;
  27. import org.apache.commons.math4.core.jdkmath.JdkMath;

  29. /**
  30.  * Base class for probability distributions on the reals.
  31.  * Default implementations are provided for some of the methods
  32.  * that do not vary from distribution to distribution.
  33.  *
  34.  * <p>
  35.  * This base class provides a default factory method for creating
  36.  * a {@link org.apache.commons.statistics.distribution.ContinuousDistribution.Sampler
  37.  * sampler instance} that uses the
  38.  * <a href="http://en.wikipedia.org/wiki/Inverse_transform_sampling">
  39.  * inversion method</a> for generating random samples that follow the
  40.  * distribution.
  41.  * </p>
  42.  *
  43.  * @since 3.0
  44.  */
  45. public abstract class AbstractRealDistribution
  46.     implements ContinuousDistribution {
  47.     /** Default absolute accuracy for inverse cumulative computation. */
  48.     public static final double SOLVER_DEFAULT_ABSOLUTE_ACCURACY = 1e-6;

  50.     /**
  51.      * For a random variable {@code X} whose values are distributed according
  52.      * to this distribution, this method returns {@code P(x0 < X <= x1)}.
  53.      *
  54.      * @param x0 Lower bound (excluded).
  55.      * @param x1 Upper bound (included).
  56.      * @return the probability that a random variable with this distribution
  57.      * takes a value between {@code x0} and {@code x1}, excluding the lower
  58.      * and including the upper endpoint.
  59.      * @throws NumberIsTooLargeException if {@code x0 > x1}.
  60.      *
  61.      * The default implementation uses the identity
  62.      * {@code P(x0 < X <= x1) = P(X <= x1) - P(X <= x0)}
  63.      *
  64.      * @since 3.1
  65.      */
  66.     @Override
  67.     public double probability(double x0,
  68.                               double x1) {
  69.         if (x0 > x1) {
  70.             throw new NumberIsTooLargeException(LocalizedFormats.LOWER_ENDPOINT_ABOVE_UPPER_ENDPOINT,
  71.                                                 x0, x1, true);
  72.         }
  73.         return cumulativeProbability(x1) - cumulativeProbability(x0);
  74.     }

  76.     /**
  77.      * {@inheritDoc}
  78.      *
  79.      * The default implementation returns
  80.      * <ul>
  81.      * <li>{@link #getSupportLowerBound()} for {@code p = 0},</li>
  82.      * <li>{@link #getSupportUpperBound()} for {@code p = 1}.</li>
  83.      * </ul>
  84.      */
  85.     @Override
  86.     public double inverseCumulativeProbability(final double p) throws OutOfRangeException {
  87.         /*
  88.          * IMPLEMENTATION NOTES
  89.          * --------------------
  90.          * Where applicable, use is made of the one-sided Chebyshev inequality
  91.          * to bracket the root. This inequality states that
  92.          * P(X - mu >= k * sig) <= 1 / (1 + k^2),
  93.          * mu: mean, sig: standard deviation. Equivalently
  94.          * 1 - P(X < mu + k * sig) <= 1 / (1 + k^2),
  95.          * F(mu + k * sig) >= k^2 / (1 + k^2).
  96.          *
  97.          * For k = sqrt(p / (1 - p)), we find
  98.          * F(mu + k * sig) >= p,
  99.          * and (mu + k * sig) is an upper-bound for the root.
  100.          *
  101.          * Then, introducing Y = -X, mean(Y) = -mu, sd(Y) = sig, and
  102.          * P(Y >= -mu + k * sig) <= 1 / (1 + k^2),
  103.          * P(-X >= -mu + k * sig) <= 1 / (1 + k^2),
  104.          * P(X <= mu - k * sig) <= 1 / (1 + k^2),
  105.          * F(mu - k * sig) <= 1 / (1 + k^2).
  106.          *
  107.          * For k = sqrt((1 - p) / p), we find
  108.          * F(mu - k * sig) <= p,
  109.          * and (mu - k * sig) is a lower-bound for the root.
  110.          *
  111.          * In cases where the Chebyshev inequality does not apply, geometric
  112.          * progressions 1, 2, 4, ... and -1, -2, -4, ... are used to bracket
  113.          * the root.
  114.          */
  115.         if (p < 0.0 || p > 1.0) {
  116.             throw new OutOfRangeException(p, 0, 1);
  117.         }

  119.         double lowerBound = getSupportLowerBound();
  120.         if (p == 0.0) {
  121.             return lowerBound;
  122.         }

  124.         double upperBound = getSupportUpperBound();
  125.         if (p == 1.0) {
  126.             return upperBound;
  127.         }

  129.         final double mu = getMean();
  130.         final double sig = JdkMath.sqrt(getVariance());
  131.         final boolean chebyshevApplies;
  132.         chebyshevApplies = !(Double.isInfinite(mu) || Double.isNaN(mu) ||
  133.                              Double.isInfinite(sig) || Double.isNaN(sig));

  135.         if (lowerBound == Double.NEGATIVE_INFINITY) {
  136.             if (chebyshevApplies) {
  137.                 lowerBound = mu - sig * JdkMath.sqrt((1. - p) / p);
  138.             } else {
  139.                 lowerBound = -1.0;
  140.                 while (cumulativeProbability(lowerBound) >= p) {
  141.                     lowerBound *= 2.0;
  142.                 }
  143.             }
  144.         }

  146.         if (upperBound == Double.POSITIVE_INFINITY) {
  147.             if (chebyshevApplies) {
  148.                 upperBound = mu + sig * JdkMath.sqrt(p / (1. - p));
  149.             } else {
  150.                 upperBound = 1.0;
  151.                 while (cumulativeProbability(upperBound) < p) {
  152.                     upperBound *= 2.0;
  153.                 }
  154.             }
  155.         }

  157.         final UnivariateFunction toSolve = new UnivariateFunction() {
  158.             /** {@inheritDoc} */
  159.             @Override
  160.             public double value(final double x) {
  161.                 return cumulativeProbability(x) - p;
  162.             }
  163.         };

  165.         return UnivariateSolverUtils.solve(toSolve,
  166.                                            lowerBound,
  167.                                            upperBound,
  168.                                            getSolverAbsoluteAccuracy());
  169.     }

  171.     /**
  172.      * Returns the solver absolute accuracy for inverse cumulative computation.
  173.      * You can override this method in order to use a Brent solver with an
  174.      * absolute accuracy different from the default.
  175.      *
  176.      * @return the maximum absolute error in inverse cumulative probability estimates
  177.      */
  178.     protected double getSolverAbsoluteAccuracy() {
  179.         return SOLVER_DEFAULT_ABSOLUTE_ACCURACY;
  180.     }

  182.     /**
  183.      * {@inheritDoc}
  184.      * <p>
  185.      * The default implementation simply computes the logarithm of {@code density(x)}.
  186.      */
  187.     @Override
  188.     public double logDensity(double x) {
  189.         return JdkMath.log(density(x));
  190.     }

  192.     /**
  193.      * Utility function for allocating an array and filling it with {@code n}
  194.      * samples generated by the given {@code sampler}.
  195.      *
  196.      * @param n Number of samples.
  197.      * @param sampler Sampler.
  198.      * @return an array of size {@code n}.
  199.      */
  200.     public static double[] sample(int n,
  201.                                   ContinuousDistribution.Sampler sampler) {
  202.         final double[] samples = new double[n];
  203.         for (int i = 0; i < n; i++) {
  204.             samples[i] = sampler.sample();
  205.         }
  206.         return samples;
  207.     }

  209.     /**{@inheritDoc} */
  210.     @Override
  211.     public ContinuousDistribution.Sampler createSampler(final UniformRandomProvider rng) {
  212.         // Inversion method distribution sampler.
  213.         return InverseTransformContinuousSampler.of(rng, this::inverseCumulativeProbability)::sample;
  214.     }
  215. }
Created with JaCoCo 0.8.10.202304240956Code Coverage Report for Apache Commons Math 4.0-SNAPSHOT