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 }