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.statistics.distribution;
019
020/**
021 * Implementation of the triangular real distribution.
022 *
023 * @see <a href="http://en.wikipedia.org/wiki/Triangular_distribution">
024 * Triangular distribution (Wikipedia)</a>
025 */
026public class TriangularDistribution extends AbstractContinuousDistribution {
027    /** Lower limit of this distribution (inclusive). */
028    private final double a;
029    /** Upper limit of this distribution (inclusive). */
030    private final double b;
031    /** Mode of this distribution. */
032    private final double c;
033
034    /**
035     * Creates a distribution.
036     *
037     * @param a Lower limit of this distribution (inclusive).
038     * @param b Upper limit of this distribution (inclusive).
039     * @param c Mode of this distribution.
040     * @throws IllegalArgumentException if {@code a >= b}, if {@code c > b}
041     * or if {@code c < a}.
042     */
043    public TriangularDistribution(double a,
044                                  double c,
045                                  double b) {
046        if (a >= b) {
047            throw new DistributionException(DistributionException.TOO_LARGE,
048                                            a, b);
049        }
050        if (c < a) {
051            throw new DistributionException(DistributionException.TOO_SMALL,
052                                            c, a);
053        }
054        if (c > b) {
055            throw new DistributionException(DistributionException.TOO_LARGE,
056                                            c, b);
057        }
058
059        this.a = a;
060        this.c = c;
061        this.b = b;
062    }
063
064    /**
065     * Gets the mode.
066     *
067     * @return the mode of the distribution.
068     */
069    public double getMode() {
070        return c;
071    }
072
073    /**
074     * {@inheritDoc}
075     *
076     * For lower limit {@code a}, upper limit {@code b} and mode {@code c}, the
077     * PDF is given by
078     * <ul>
079     * <li>{@code 2 * (x - a) / [(b - a) * (c - a)]} if {@code a <= x < c},</li>
080     * <li>{@code 2 / (b - a)} if {@code x = c},</li>
081     * <li>{@code 2 * (b - x) / [(b - a) * (b - c)]} if {@code c < x <= b},</li>
082     * <li>{@code 0} otherwise.
083     * </ul>
084     */
085    @Override
086    public double density(double x) {
087        if (x < a) {
088            return 0;
089        }
090        if (x < c) {
091            final double divident = 2 * (x - a);
092            final double divisor = (b - a) * (c - a);
093            return divident / divisor;
094        }
095        if (x == c) {
096            return 2 / (b - a);
097        }
098        if (x <= b) {
099            final double divident = 2 * (b - x);
100            final double divisor = (b - a) * (b - c);
101            return divident / divisor;
102        }
103        return 0;
104    }
105
106    /**
107     * {@inheritDoc}
108     *
109     * For lower limit {@code a}, upper limit {@code b} and mode {@code c}, the
110     * CDF is given by
111     * <ul>
112     * <li>{@code 0} if {@code x < a},</li>
113     * <li>{@code (x - a)^2 / [(b - a) * (c - a)]} if {@code a <= x < c},</li>
114     * <li>{@code (c - a) / (b - a)} if {@code x = c},</li>
115     * <li>{@code 1 - (b - x)^2 / [(b - a) * (b - c)]} if {@code c < x <= b},</li>
116     * <li>{@code 1} if {@code x > b}.</li>
117     * </ul>
118     */
119    @Override
120    public double cumulativeProbability(double x)  {
121        if (x < a) {
122            return 0;
123        }
124        if (x < c) {
125            final double divident = (x - a) * (x - a);
126            final double divisor = (b - a) * (c - a);
127            return divident / divisor;
128        }
129        if (x == c) {
130            return (c - a) / (b - a);
131        }
132        if (x <= b) {
133            final double divident = (b - x) * (b - x);
134            final double divisor = (b - a) * (b - c);
135            return 1 - (divident / divisor);
136        }
137        return 1;
138    }
139
140    /**
141     * {@inheritDoc}
142     *
143     * For lower limit {@code a}, upper limit {@code b}, and mode {@code c},
144     * the mean is {@code (a + b + c) / 3}.
145     */
146    @Override
147    public double getMean() {
148        return (a + b + c) / 3;
149    }
150
151    /**
152     * {@inheritDoc}
153     *
154     * For lower limit {@code a}, upper limit {@code b}, and mode {@code c},
155     * the variance is {@code (a^2 + b^2 + c^2 - a * b - a * c - b * c) / 18}.
156     */
157    @Override
158    public double getVariance() {
159        return (a * a + b * b + c * c - a * b - a * c - b * c) / 18;
160    }
161
162    /**
163     * {@inheritDoc}
164     *
165     * The lower bound of the support is equal to the lower limit parameter
166     * {@code a} of the distribution.
167     *
168     * @return lower bound of the support
169     */
170    @Override
171    public double getSupportLowerBound() {
172        return a;
173    }
174
175    /**
176     * {@inheritDoc}
177     *
178     * The upper bound of the support is equal to the upper limit parameter
179     * {@code b} of the distribution.
180     *
181     * @return upper bound of the support
182     */
183    @Override
184    public double getSupportUpperBound() {
185        return b;
186    }
187
188    /**
189     * {@inheritDoc}
190     *
191     * The support of this distribution is connected.
192     *
193     * @return {@code true}
194     */
195    @Override
196    public boolean isSupportConnected() {
197        return true;
198    }
199
200    /** {@inheritDoc} */
201    @Override
202    public double inverseCumulativeProbability(double p) {
203        if (p < 0 ||
204            p > 1) {
205            throw new DistributionException(DistributionException.INVALID_PROBABILITY, p);
206        }
207        if (p == 0) {
208            return a;
209        }
210        if (p == 1) {
211            return b;
212        }
213        if (p < (c - a) / (b - a)) {
214            return a + Math.sqrt(p * (b - a) * (c - a));
215        }
216        return b - Math.sqrt((1 - p) * (b - a) * (b - c));
217    }
218}