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.  *      https://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.  */

  18. /*
  19.  * This is not the original file distributed by the Apache Software Foundation
  20.  * It has been modified by the Hipparchus project
  21.  */
  22. package org.hipparchus.distribution.continuous;

  24. import java.io.Serializable;

  26. import org.hipparchus.analysis.UnivariateFunction;
  27. import org.hipparchus.analysis.solvers.UnivariateSolverUtils;
  28. import org.hipparchus.distribution.RealDistribution;
  29. import org.hipparchus.exception.LocalizedCoreFormats;
  30. import org.hipparchus.exception.MathIllegalArgumentException;
  31. import org.hipparchus.util.FastMath;
  32. import org.hipparchus.util.MathUtils;

  34. /**
  35.  * Base class for probability distributions on the reals.
  36.  * <p>
  37.  * Default implementations are provided for some of the methods
  38.  * that do not vary from distribution to distribution.
  39.  */
  40. public abstract class AbstractRealDistribution
  41.     implements RealDistribution, Serializable {

  43.     /** Default absolute accuracy for inverse cumulative computation. */
  44.     protected static final double DEFAULT_SOLVER_ABSOLUTE_ACCURACY = 1e-9;
  45.     /** Serializable version identifier */
  46.     private static final long serialVersionUID = 20160320L;

  48.     /** Inverse cumulative probability accuracy. */
  49.     private final double solverAbsoluteAccuracy;

  51.     /** Simple constructor.
  52.      * @param solverAbsoluteAccuracy the absolute accuracy to use when
  53.      * computing the inverse cumulative probability.
  54.      */
  55.     protected AbstractRealDistribution(double solverAbsoluteAccuracy) {
  56.         this.solverAbsoluteAccuracy = solverAbsoluteAccuracy;
  57.     }

  59.     /**
  60.      * Create a real distribution with default solver absolute accuracy.
  61.      */
  62.     protected AbstractRealDistribution() {
  63.         this.solverAbsoluteAccuracy = DEFAULT_SOLVER_ABSOLUTE_ACCURACY;
  64.     }

  66.     /**
  67.      * For a random variable {@code X} whose values are distributed according
  68.      * to this distribution, this method returns {@code P(x0 < X <= x1)}.
  69.      *
  70.      * @param x0 Lower bound (excluded).
  71.      * @param x1 Upper bound (included).
  72.      * @return the probability that a random variable with this distribution
  73.      * takes a value between {@code x0} and {@code x1}, excluding the lower
  74.      * and including the upper endpoint.
  75.      * @throws MathIllegalArgumentException if {@code x0 > x1}.
  76.      * The default implementation uses the identity
  77.      * {@code P(x0 < X <= x1) = P(X <= x1) - P(X <= x0)}
  78.      */
  79.     @Override
  80.     public double probability(double x0,
  81.                               double x1) throws MathIllegalArgumentException {
  82.         if (x0 > x1) {
  83.             throw new MathIllegalArgumentException(LocalizedCoreFormats.LOWER_ENDPOINT_ABOVE_UPPER_ENDPOINT,
  84.                                                    x0, x1, true);
  85.         }
  86.         return cumulativeProbability(x1) - cumulativeProbability(x0);
  87.     }

  89.     /**
  90.      * {@inheritDoc}
  91.      *
  92.      * The default implementation returns
  93.      * <ul>
  94.      * <li>{@link #getSupportLowerBound()} for {@code p = 0},</li>
  95.      * <li>{@link #getSupportUpperBound()} for {@code p = 1}.</li>
  96.      * </ul>
  97.      */
  98.     @Override
  99.     public double inverseCumulativeProbability(final double p) throws MathIllegalArgumentException {
  100.         /*
  101.          * IMPLEMENTATION NOTES
  102.          * --------------------
  103.          * Where applicable, use is made of the one-sided Chebyshev inequality
  104.          * to bracket the root. This inequality states that
  105.          * P(X - mu >= k * sig) <= 1 / (1 + k^2),
  106.          * mu: mean, sig: standard deviation. Equivalently
  107.          * 1 - P(X < mu + k * sig) <= 1 / (1 + k^2),
  108.          * F(mu + k * sig) >= k^2 / (1 + k^2).
  109.          *
  110.          * For k = sqrt(p / (1 - p)), we find
  111.          * F(mu + k * sig) >= p,
  112.          * and (mu + k * sig) is an upper-bound for the root.
  113.          *
  114.          * Then, introducing Y = -X, mean(Y) = -mu, sd(Y) = sig, and
  115.          * P(Y >= -mu + k * sig) <= 1 / (1 + k^2),
  116.          * P(-X >= -mu + k * sig) <= 1 / (1 + k^2),
  117.          * P(X <= mu - k * sig) <= 1 / (1 + k^2),
  118.          * F(mu - k * sig) <= 1 / (1 + k^2).
  119.          *
  120.          * For k = sqrt((1 - p) / p), we find
  121.          * F(mu - k * sig) <= p,
  122.          * and (mu - k * sig) is a lower-bound for the root.
  123.          *
  124.          * In cases where the Chebyshev inequality does not apply, geometric
  125.          * progressions 1, 2, 4, ... and -1, -2, -4, ... are used to bracket
  126.          * the root.
  127.          */

  129.         MathUtils.checkRangeInclusive(p, 0, 1);

  131.         double lowerBound = getSupportLowerBound();
  132.         if (p == 0.0) {
  133.             return lowerBound;
  134.         }

  136.         double upperBound = getSupportUpperBound();
  137.         if (p == 1.0) {
  138.             return upperBound;
  139.         }

  141.         final double mu = getNumericalMean();
  142.         final double sig = FastMath.sqrt(getNumericalVariance());
  143.         final boolean chebyshevApplies;
  144.         chebyshevApplies = !(Double.isInfinite(mu) || Double.isNaN(mu) ||
  145.                              Double.isInfinite(sig) || Double.isNaN(sig));

  147.         if (lowerBound == Double.NEGATIVE_INFINITY) {
  148.             if (chebyshevApplies) {
  149.                 lowerBound = mu - sig * FastMath.sqrt((1. - p) / p);
  150.             } else {
  151.                 lowerBound = -1.0;
  152.                 while (cumulativeProbability(lowerBound) >= p) {
  153.                     lowerBound *= 2.0;
  154.                 }
  155.             }
  156.         }

  158.         if (upperBound == Double.POSITIVE_INFINITY) {
  159.             if (chebyshevApplies) {
  160.                 upperBound = mu + sig * FastMath.sqrt(p / (1. - p));
  161.             } else {
  162.                 upperBound = 1.0;
  163.                 while (cumulativeProbability(upperBound) < p) {
  164.                     upperBound *= 2.0;
  165.                 }
  166.             }
  167.         }

  169.         final UnivariateFunction toSolve = new UnivariateFunction() {
  170.             /** {@inheritDoc} */
  171.             @Override
  172.             public double value(final double x) {
  173.                 return cumulativeProbability(x) - p;
  174.             }
  175.         };

  177.         double x = UnivariateSolverUtils.solve(toSolve,
  178.                                                lowerBound,
  179.                                                upperBound,
  180.                                                getSolverAbsoluteAccuracy());

  182.         if (!isSupportConnected()) {
  183.             /* Test for plateau. */
  184.             final double dx = getSolverAbsoluteAccuracy();
  185.             if (x - dx >= getSupportLowerBound()) {
  186.                 double px = cumulativeProbability(x);
  187.                 if (cumulativeProbability(x - dx) == px) {
  188.                     upperBound = x;
  189.                     while (upperBound - lowerBound > dx) {
  190.                         final double midPoint = 0.5 * (lowerBound + upperBound);
  191.                         if (cumulativeProbability(midPoint) < px) {
  192.                             lowerBound = midPoint;
  193.                         } else {
  194.                             upperBound = midPoint;
  195.                         }
  196.                     }
  197.                     return upperBound;
  198.                 }
  199.             }
  200.         }
  201.         return x;
  202.     }

  204.     /**
  205.      * Returns the solver absolute accuracy for inverse cumulative computation.
  206.      * You can override this method in order to use a Brent solver with an
  207.      * absolute accuracy different from the default.
  208.      *
  209.      * @return the maximum absolute error in inverse cumulative probability estimates
  210.      */
  211.     protected double getSolverAbsoluteAccuracy() {
  212.         return solverAbsoluteAccuracy;
  213.     }

  215.     /**
  216.      * {@inheritDoc}
  217.      * <p>
  218.      * The default implementation simply computes the logarithm of {@code density(x)}.
  219.      */
  220.     @Override
  221.     public double logDensity(double x) {
  222.         return FastMath.log(density(x));
  223.     }
  224. }
Created with JaCoCo 0.8.12.202403310830