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    }