TriangularDistribution.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.  */

  18. package org.apache.commons.math3.distribution;

  20. import org.apache.commons.math3.exception.NumberIsTooLargeException;
  21. import org.apache.commons.math3.exception.NumberIsTooSmallException;
  22. import org.apache.commons.math3.exception.OutOfRangeException;
  23. import org.apache.commons.math3.exception.util.LocalizedFormats;
  24. import org.apache.commons.math3.random.RandomGenerator;
  25. import org.apache.commons.math3.random.Well19937c;
  26. import org.apache.commons.math3.util.FastMath;

  28. /**
  29.  * Implementation of the triangular real distribution.
  30.  *
  31.  * @see <a href="http://en.wikipedia.org/wiki/Triangular_distribution">
  32.  * Triangular distribution (Wikipedia)</a>
  33.  *
  34.  * @since 3.0
  35.  */
  36. public class TriangularDistribution extends AbstractRealDistribution {
  37.     /** Serializable version identifier. */
  38.     private static final long serialVersionUID = 20120112L;
  39.     /** Lower limit of this distribution (inclusive). */
  40.     private final double a;
  41.     /** Upper limit of this distribution (inclusive). */
  42.     private final double b;
  43.     /** Mode of this distribution. */
  44.     private final double c;
  45.     /** Inverse cumulative probability accuracy. */
  46.     private final double solverAbsoluteAccuracy;

  48.     /**
  49.      * Creates a triangular real distribution using the given lower limit,
  50.      * upper limit, and mode.
  51.      * <p>
  52.      * <b>Note:</b> this constructor will implicitly create an instance of
  53.      * {@link Well19937c} as random generator to be used for sampling only (see
  54.      * {@link #sample()} and {@link #sample(int)}). In case no sampling is
  55.      * needed for the created distribution, it is advised to pass {@code null}
  56.      * as random generator via the appropriate constructors to avoid the
  57.      * additional initialisation overhead.
  58.      *
  59.      * @param a Lower limit of this distribution (inclusive).
  60.      * @param b Upper limit of this distribution (inclusive).
  61.      * @param c Mode of this distribution.
  62.      * @throws NumberIsTooLargeException if {@code a >= b} or if {@code c > b}.
  63.      * @throws NumberIsTooSmallException if {@code c < a}.
  64.      */
  65.     public TriangularDistribution(double a, double c, double b)
  66.         throws NumberIsTooLargeException, NumberIsTooSmallException {
  67.         this(new Well19937c(), a, c, b);
  68.     }

  70.     /**
  71.      * Creates a triangular distribution.
  72.      *
  73.      * @param rng Random number generator.
  74.      * @param a Lower limit of this distribution (inclusive).
  75.      * @param b Upper limit of this distribution (inclusive).
  76.      * @param c Mode of this distribution.
  77.      * @throws NumberIsTooLargeException if {@code a >= b} or if {@code c > b}.
  78.      * @throws NumberIsTooSmallException if {@code c < a}.
  79.      * @since 3.1
  80.      */
  81.     public TriangularDistribution(RandomGenerator rng,
  82.                                   double a,
  83.                                   double c,
  84.                                   double b)
  85.         throws NumberIsTooLargeException, NumberIsTooSmallException {
  86.         super(rng);

  88.         if (a >= b) {
  89.             throw new NumberIsTooLargeException(
  90.                             LocalizedFormats.LOWER_BOUND_NOT_BELOW_UPPER_BOUND,
  91.                             a, b, false);
  92.         }
  93.         if (c < a) {
  94.             throw new NumberIsTooSmallException(
  95.                     LocalizedFormats.NUMBER_TOO_SMALL, c, a, true);
  96.         }
  97.         if (c > b) {
  98.             throw new NumberIsTooLargeException(
  99.                     LocalizedFormats.NUMBER_TOO_LARGE, c, b, true);
  100.         }

  102.         this.a = a;
  103.         this.c = c;
  104.         this.b = b;
  105.         solverAbsoluteAccuracy = FastMath.max(FastMath.ulp(a), FastMath.ulp(b));
  106.     }

  108.     /**
  109.      * Returns the mode {@code c} of this distribution.
  110.      *
  111.      * @return the mode {@code c} of this distribution
  112.      */
  113.     public double getMode() {
  114.         return c;
  115.     }

  117.     /**
  118.      * {@inheritDoc}
  119.      *
  120.      * <p>
  121.      * For this distribution, the returned value is not really meaningful,
  122.      * since exact formulas are implemented for the computation of the
  123.      * {@link #inverseCumulativeProbability(double)} (no solver is invoked).
  124.      * </p>
  125.      * <p>
  126.      * For lower limit {@code a} and upper limit {@code b}, the current
  127.      * implementation returns {@code max(ulp(a), ulp(b)}.
  128.      * </p>
  129.      */
  130.     @Override
  131.     protected double getSolverAbsoluteAccuracy() {
  132.         return solverAbsoluteAccuracy;
  133.     }

  135.     /**
  136.      * {@inheritDoc}
  137.      *
  138.      * For lower limit {@code a}, upper limit {@code b} and mode {@code c}, the
  139.      * PDF is given by
  140.      * <ul>
  141.      * <li>{@code 2 * (x - a) / [(b - a) * (c - a)]} if {@code a <= x < c},</li>
  142.      * <li>{@code 2 / (b - a)} if {@code x = c},</li>
  143.      * <li>{@code 2 * (b - x) / [(b - a) * (b - c)]} if {@code c < x <= b},</li>
  144.      * <li>{@code 0} otherwise.
  145.      * </ul>
  146.      */
  147.     public double density(double x) {
  148.         if (x < a) {
  149.             return 0;
  150.         }
  151.         if (a <= x && x < c) {
  152.             double divident = 2 * (x - a);
  153.             double divisor = (b - a) * (c - a);
  154.             return divident / divisor;
  155.         }
  156.         if (x == c) {
  157.             return 2 / (b - a);
  158.         }
  159.         if (c < x && x <= b) {
  160.             double divident = 2 * (b - x);
  161.             double divisor = (b - a) * (b - c);
  162.             return divident / divisor;
  163.         }
  164.         return 0;
  165.     }

  167.     /**
  168.      * {@inheritDoc}
  169.      *
  170.      * For lower limit {@code a}, upper limit {@code b} and mode {@code c}, the
  171.      * CDF is given by
  172.      * <ul>
  173.      * <li>{@code 0} if {@code x < a},</li>
  174.      * <li>{@code (x - a)^2 / [(b - a) * (c - a)]} if {@code a <= x < c},</li>
  175.      * <li>{@code (c - a) / (b - a)} if {@code x = c},</li>
  176.      * <li>{@code 1 - (b - x)^2 / [(b - a) * (b - c)]} if {@code c < x <= b},</li>
  177.      * <li>{@code 1} if {@code x > b}.</li>
  178.      * </ul>
  179.      */
  180.     public double cumulativeProbability(double x)  {
  181.         if (x < a) {
  182.             return 0;
  183.         }
  184.         if (a <= x && x < c) {
  185.             double divident = (x - a) * (x - a);
  186.             double divisor = (b - a) * (c - a);
  187.             return divident / divisor;
  188.         }
  189.         if (x == c) {
  190.             return (c - a) / (b - a);
  191.         }
  192.         if (c < x && x <= b) {
  193.             double divident = (b - x) * (b - x);
  194.             double divisor = (b - a) * (b - c);
  195.             return 1 - (divident / divisor);
  196.         }
  197.         return 1;
  198.     }

  200.     /**
  201.      * {@inheritDoc}
  202.      *
  203.      * For lower limit {@code a}, upper limit {@code b}, and mode {@code c},
  204.      * the mean is {@code (a + b + c) / 3}.
  205.      */
  206.     public double getNumericalMean() {
  207.         return (a + b + c) / 3;
  208.     }

  210.     /**
  211.      * {@inheritDoc}
  212.      *
  213.      * For lower limit {@code a}, upper limit {@code b}, and mode {@code c},
  214.      * the variance is {@code (a^2 + b^2 + c^2 - a * b - a * c - b * c) / 18}.
  215.      */
  216.     public double getNumericalVariance() {
  217.         return (a * a + b * b + c * c - a * b - a * c - b * c) / 18;
  218.     }

  220.     /**
  221.      * {@inheritDoc}
  222.      *
  223.      * The lower bound of the support is equal to the lower limit parameter
  224.      * {@code a} of the distribution.
  225.      *
  226.      * @return lower bound of the support
  227.      */
  228.     public double getSupportLowerBound() {
  229.         return a;
  230.     }

  232.     /**
  233.      * {@inheritDoc}
  234.      *
  235.      * The upper bound of the support is equal to the upper limit parameter
  236.      * {@code b} of the distribution.
  237.      *
  238.      * @return upper bound of the support
  239.      */
  240.     public double getSupportUpperBound() {
  241.         return b;
  242.     }

  244.     /** {@inheritDoc} */
  245.     public boolean isSupportLowerBoundInclusive() {
  246.         return true;
  247.     }

  249.     /** {@inheritDoc} */
  250.     public boolean isSupportUpperBoundInclusive() {
  251.         return true;
  252.     }

  254.     /**
  255.      * {@inheritDoc}
  256.      *
  257.      * The support of this distribution is connected.
  258.      *
  259.      * @return {@code true}
  260.      */
  261.     public boolean isSupportConnected() {
  262.         return true;
  263.     }

  265.     @Override
  266.     public double inverseCumulativeProbability(double p)
  267.         throws OutOfRangeException {
  268.         if (p < 0 || p > 1) {
  269.             throw new OutOfRangeException(p, 0, 1);
  270.         }
  271.         if (p == 0) {
  272.             return a;
  273.         }
  274.         if (p == 1) {
  275.             return b;
  276.         }
  277.         if (p < (c - a) / (b - a)) {
  278.             return a + FastMath.sqrt(p * (b - a) * (c - a));
  279.         }
  280.         return b - FastMath.sqrt((1 - p) * (b - a) * (b - c));
  281.     }
  282. }
Created with JaCoCo 0.7.2.201409121644