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
018package org.apache.commons.math3.distribution;
019
020import org.apache.commons.math3.exception.NumberIsTooLargeException;
021import org.apache.commons.math3.exception.NumberIsTooSmallException;
022import org.apache.commons.math3.exception.OutOfRangeException;
023import org.apache.commons.math3.exception.util.LocalizedFormats;
024import org.apache.commons.math3.random.RandomGenerator;
025import org.apache.commons.math3.random.Well19937c;
026import org.apache.commons.math3.util.FastMath;
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 * @since 3.0
035 */
036public class TriangularDistribution extends AbstractRealDistribution {
037    /** Serializable version identifier. */
038    private static final long serialVersionUID = 20120112L;
039    /** Lower limit of this distribution (inclusive). */
040    private final double a;
041    /** Upper limit of this distribution (inclusive). */
042    private final double b;
043    /** Mode of this distribution. */
044    private final double c;
045    /** Inverse cumulative probability accuracy. */
046    private final double solverAbsoluteAccuracy;
047
048    /**
049     * Creates a triangular real distribution using the given lower limit,
050     * upper limit, and mode.
051     * <p>
052     * <b>Note:</b> this constructor will implicitly create an instance of
053     * {@link Well19937c} as random generator to be used for sampling only (see
054     * {@link #sample()} and {@link #sample(int)}). In case no sampling is
055     * needed for the created distribution, it is advised to pass {@code null}
056     * as random generator via the appropriate constructors to avoid the
057     * additional initialisation overhead.
058     *
059     * @param a Lower limit of this distribution (inclusive).
060     * @param b Upper limit of this distribution (inclusive).
061     * @param c Mode of this distribution.
062     * @throws NumberIsTooLargeException if {@code a >= b} or if {@code c > b}.
063     * @throws NumberIsTooSmallException if {@code c < a}.
064     */
065    public TriangularDistribution(double a, double c, double b)
066        throws NumberIsTooLargeException, NumberIsTooSmallException {
067        this(new Well19937c(), a, c, b);
068    }
069
070    /**
071     * Creates a triangular distribution.
072     *
073     * @param rng Random number generator.
074     * @param a Lower limit of this distribution (inclusive).
075     * @param b Upper limit of this distribution (inclusive).
076     * @param c Mode of this distribution.
077     * @throws NumberIsTooLargeException if {@code a >= b} or if {@code c > b}.
078     * @throws NumberIsTooSmallException if {@code c < a}.
079     * @since 3.1
080     */
081    public TriangularDistribution(RandomGenerator rng,
082                                  double a,
083                                  double c,
084                                  double b)
085        throws NumberIsTooLargeException, NumberIsTooSmallException {
086        super(rng);
087
088        if (a >= b) {
089            throw new NumberIsTooLargeException(
090                            LocalizedFormats.LOWER_BOUND_NOT_BELOW_UPPER_BOUND,
091                            a, b, false);
092        }
093        if (c < a) {
094            throw new NumberIsTooSmallException(
095                    LocalizedFormats.NUMBER_TOO_SMALL, c, a, true);
096        }
097        if (c > b) {
098            throw new NumberIsTooLargeException(
099                    LocalizedFormats.NUMBER_TOO_LARGE, c, b, true);
100        }
101
102        this.a = a;
103        this.c = c;
104        this.b = b;
105        solverAbsoluteAccuracy = FastMath.max(FastMath.ulp(a), FastMath.ulp(b));
106    }
107
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    }
116
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    }
134
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    }
166
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    }
199
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    }
209
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    }
219
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    }
231
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    }
243
244    /** {@inheritDoc} */
245    public boolean isSupportLowerBoundInclusive() {
246        return true;
247    }
248
249    /** {@inheritDoc} */
250    public boolean isSupportUpperBoundInclusive() {
251        return true;
252    }
253
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    }
264
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}