    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.distribution;
    019
    020    import org.apache.commons.math3.exception.NumberIsTooLargeException;
    021    import org.apache.commons.math3.exception.NumberIsTooSmallException;
    022    import org.apache.commons.math3.exception.OutOfRangeException;
    023    import org.apache.commons.math3.exception.util.LocalizedFormats;
    024    import org.apache.commons.math3.util.FastMath;
    025    import org.apache.commons.math3.random.RandomGenerator;
    026    import org.apache.commons.math3.random.Well19937c;
    027
    028    /**
    029     * Implementation of the triangular real distribution.
    030     *
    031     * @see <a href="http://en.wikipedia.org/wiki/Triangular_distribution">
    032     * Triangular distribution (Wikipedia)</a>
    033     *
    034     * @version $Id: TriangularDistribution.java 1416643 2012-12-03 19:37:14Z tn $
    035     * @since 3.0
    036     */
    037    public class TriangularDistribution extends AbstractRealDistribution {
    038        /** Serializable version identifier. */
    039        private static final long serialVersionUID = 20120112L;
    040        /** Lower limit of this distribution (inclusive). */
    041        private final double a;
    042        /** Upper limit of this distribution (inclusive). */
    043        private final double b;
    044        /** Mode of this distribution. */
    045        private final double c;
    046        /** Inverse cumulative probability accuracy. */
    047        private final double solverAbsoluteAccuracy;
    048
    049        /**
    050         * Creates a triangular real distribution using the given lower limit,
    051         * upper limit, and mode.
    052         *
    053         * @param a Lower limit of this distribution (inclusive).
    054         * @param b Upper limit of this distribution (inclusive).
    055         * @param c Mode of this distribution.
    056         * @throws NumberIsTooLargeException if {@code a >= b} or if {@code c > b}.
    057         * @throws NumberIsTooSmallException if {@code c < a}.
    058         */
    059        public TriangularDistribution(double a, double c, double b)
    060            throws NumberIsTooLargeException, NumberIsTooSmallException {
    061            this(new Well19937c(), a, c, b);
    062        }
    063
    064        /**
    065         * Creates a triangular distribution.
    066         *
    067         * @param rng Random number generator.
    068         * @param a Lower limit of this distribution (inclusive).
    069         * @param b Upper limit of this distribution (inclusive).
    070         * @param c Mode of this distribution.
    071         * @throws NumberIsTooLargeException if {@code a >= b} or if {@code c > b}.
    072         * @throws NumberIsTooSmallException if {@code c < a}.
    073         * @since 3.1
    074         */
    075        public TriangularDistribution(RandomGenerator rng,
    076                                      double a,
    077                                      double c,
    078                                      double b)
    079            throws NumberIsTooLargeException, NumberIsTooSmallException {
    080            super(rng);
    081
    082            if (a >= b) {
    083                throw new NumberIsTooLargeException(
    084                                LocalizedFormats.LOWER_BOUND_NOT_BELOW_UPPER_BOUND,
    085                                a, b, false);
    086            }
    087            if (c < a) {
    088                throw new NumberIsTooSmallException(
    089                        LocalizedFormats.NUMBER_TOO_SMALL, c, a, true);
    090            }
    091            if (c > b) {
    092                throw new NumberIsTooLargeException(
    093                        LocalizedFormats.NUMBER_TOO_LARGE, c, b, true);
    094            }
    095
    096            this.a = a;
    097            this.c = c;
    098            this.b = b;
    099            solverAbsoluteAccuracy = FastMath.max(FastMath.ulp(a), FastMath.ulp(b));
    100        }
    101
    102        /**
    103         * Returns the mode {@code c} of this distribution.
    104         *
    105         * @return the mode {@code c} of this distribution
    106         */
    107        public double getMode() {
    108            return c;
    109        }
    110
    111        /**
    112         * {@inheritDoc}
    113         *
    114         * <p>
    115         * For this distribution, the returned value is not really meaningful,
    116         * since exact formulas are implemented for the computation of the
    117         * {@link #inverseCumulativeProbability(double)} (no solver is invoked).
    118         * </p>
    119         * <p>
    120         * For lower limit {@code a} and upper limit {@code b}, the current
    121         * implementation returns {@code max(ulp(a), ulp(b)}.
    122         * </p>
    123         */
    124        @Override
    125        protected double getSolverAbsoluteAccuracy() {
    126            return solverAbsoluteAccuracy;
    127        }
    128
    129        /**
    130         * {@inheritDoc}
    131         *
    132         * For lower limit {@code a}, upper limit {@code b} and mode {@code c}, the
    133         * PDF is given by
    134         * <ul>
    135         * <li>{@code 2 * (x - a) / [(b - a) * (c - a)]} if {@code a <= x < c},</li>
    136         * <li>{@code 2 / (b - a)} if {@code x = c},</li>
    137         * <li>{@code 2 * (b - x) / [(b - a) * (b - c)]} if {@code c < x <= b},</li>
    138         * <li>{@code 0} otherwise.
    139         * </ul>
    140         */
    141        public double density(double x) {
    142            if (x < a) {
    143                return 0;
    144            }
    145            if (a <= x && x < c) {
    146                double divident = 2 * (x - a);
    147                double divisor = (b - a) * (c - a);
    148                return divident / divisor;
    149            }
    150            if (x == c) {
    151                return 2 / (b - a);
    152            }
    153            if (c < x && x <= b) {
    154                double divident = 2 * (b - x);
    155                double divisor = (b - a) * (b - c);
    156                return divident / divisor;
    157            }
    158            return 0;
    159        }
    160
    161        /**
    162         * {@inheritDoc}
    163         *
    164         * For lower limit {@code a}, upper limit {@code b} and mode {@code c}, the
    165         * CDF is given by
    166         * <ul>
    167         * <li>{@code 0} if {@code x < a},</li>
    168         * <li>{@code (x - a)^2 / [(b - a) * (c - a)]} if {@code a <= x < c},</li>
    169         * <li>{@code (c - a) / (b - a)} if {@code x = c},</li>
    170         * <li>{@code 1 - (b - x)^2 / [(b - a) * (b - c)]} if {@code c < x <= b},</li>
    171         * <li>{@code 1} if {@code x > b}.</li>
    172         * </ul>
    173         */
    174        public double cumulativeProbability(double x)  {
    175            if (x < a) {
    176                return 0;
    177            }
    178            if (a <= x && x < c) {
    179                double divident = (x - a) * (x - a);
    180                double divisor = (b - a) * (c - a);
    181                return divident / divisor;
    182            }
    183            if (x == c) {
    184                return (c - a) / (b - a);
    185            }
    186            if (c < x && x <= b) {
    187                double divident = (b - x) * (b - x);
    188                double divisor = (b - a) * (b - c);
    189                return 1 - (divident / divisor);
    190            }
    191            return 1;
    192        }
    193
    194        /**
    195         * {@inheritDoc}
    196         *
    197         * For lower limit {@code a}, upper limit {@code b}, and mode {@code c},
    198         * the mean is {@code (a + b + c) / 3}.
    199         */
    200        public double getNumericalMean() {
    201            return (a + b + c) / 3;
    202        }
    203
    204        /**
    205         * {@inheritDoc}
    206         *
    207         * For lower limit {@code a}, upper limit {@code b}, and mode {@code c},
    208         * the variance is {@code (a^2 + b^2 + c^2 - a * b - a * c - b * c) / 18}.
    209         */
    210        public double getNumericalVariance() {
    211            return (a * a + b * b + c * c - a * b - a * c - b * c) / 18;
    212        }
    213
    214        /**
    215         * {@inheritDoc}
    216         *
    217         * The lower bound of the support is equal to the lower limit parameter
    218         * {@code a} of the distribution.
    219         *
    220         * @return lower bound of the support
    221         */
    222        public double getSupportLowerBound() {
    223            return a;
    224        }
    225
    226        /**
    227         * {@inheritDoc}
    228         *
    229         * The upper bound of the support is equal to the upper limit parameter
    230         * {@code b} of the distribution.
    231         *
    232         * @return upper bound of the support
    233         */
    234        public double getSupportUpperBound() {
    235            return b;
    236        }
    237
    238        /** {@inheritDoc} */
    239        public boolean isSupportLowerBoundInclusive() {
    240            return true;
    241        }
    242
    243        /** {@inheritDoc} */
    244        public boolean isSupportUpperBoundInclusive() {
    245            return true;
    246        }
    247
    248        /**
    249         * {@inheritDoc}
    250         *
    251         * The support of this distribution is connected.
    252         *
    253         * @return {@code true}
    254         */
    255        public boolean isSupportConnected() {
    256            return true;
    257        }
    258
    259        @Override
    260        public double inverseCumulativeProbability(double p)
    261            throws OutOfRangeException {
    262            if (p < 0 || p > 1) {
    263                throw new OutOfRangeException(p, 0, 1);
    264            }
    265            if (p == 0) {
    266                return a;
    267            }
    268            if (p == 1) {
    269                return b;
    270            }
    271            if (p < (c - a) / (b - a)) {
    272                return a + FastMath.sqrt(p * (b - a) * (c - a));
    273            }
    274            return b - FastMath.sqrt((1 - p) * (b - a) * (b - c));
    275        }
    276    }