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
020import org.apache.commons.rng.UniformRandomProvider;
021
022/**
023 * Implementation of the trapezoidal distribution.
024 *
025 * <p>The probability density function of \( X \) is:
026 *
027 * <p>\[ f(x; a, b, c, d) = \begin{cases}
028 *       \frac{2}{d+c-a-b}\frac{x-a}{b-a} &amp; \text{for } a\le x \lt b \\
029 *       \frac{2}{d+c-a-b}                &amp; \text{for } b\le x \lt c \\
030 *       \frac{2}{d+c-a-b}\frac{d-x}{d-c} &amp; \text{for } c\le x \le d
031 *       \end{cases} \]
032 *
033 * <p>for \( -\infty \lt a \le b \le c \le d \lt \infty \) and
034 * \( x \in [a, d] \).
035 *
036 * <p>Note the special cases:
037 * <ul>
038 * <li>\( b = c \) is the triangular distribution
039 * <li>\( a = b \) and \( c = d \) is the uniform distribution
040 * </ul>
041 *
042 * @see <a href="https://en.wikipedia.org/wiki/Trapezoidal_distribution">Trapezoidal distribution (Wikipedia)</a>
043 */
044public abstract class TrapezoidalDistribution extends AbstractContinuousDistribution {
045    /** Lower limit of this distribution (inclusive). */
046    protected final double a;
047    /** Start of the trapezoid constant density. */
048    protected final double b;
049    /** End of the trapezoid constant density. */
050    protected final double c;
051    /** Upper limit of this distribution (inclusive). */
052    protected final double d;
053
054    /**
055     * Specialisation of the trapezoidal distribution used when the distribution simplifies
056     * to an alternative distribution.
057     */
058    private static class DelegatedTrapezoidalDistribution extends TrapezoidalDistribution {
059        /** Distribution delegate. */
060        private final ContinuousDistribution delegate;
061
062        /**
063         * @param a Lower limit of this distribution (inclusive).
064         * @param b Start of the trapezoid constant density.
065         * @param c End of the trapezoid constant density.
066         * @param d Upper limit of this distribution (inclusive).
067         * @param delegate Distribution delegate.
068         */
069        DelegatedTrapezoidalDistribution(double a, double b, double c, double d,
070                                         ContinuousDistribution delegate) {
071            super(a, b, c, d);
072            this.delegate = delegate;
073        }
074
075        @Override
076        public double density(double x) {
077            return delegate.density(x);
078        }
079
080        @Override
081        public double probability(double x0, double x1) {
082            return delegate.probability(x0, x1);
083        }
084
085        @Override
086        public double logDensity(double x) {
087            return delegate.logDensity(x);
088        }
089
090        @Override
091        public double cumulativeProbability(double x) {
092            return delegate.cumulativeProbability(x);
093        }
094
095        @Override
096        public double inverseCumulativeProbability(double p) {
097            return delegate.inverseCumulativeProbability(p);
098        }
099
100        @Override
101        public double survivalProbability(double x) {
102            return delegate.survivalProbability(x);
103        }
104
105        @Override
106        public double inverseSurvivalProbability(double p) {
107            return delegate.inverseSurvivalProbability(p);
108        }
109
110        @Override
111        public double getMean() {
112            return delegate.getMean();
113        }
114
115        @Override
116        public double getVariance() {
117            return delegate.getVariance();
118        }
119
120        @Override
121        public Sampler createSampler(UniformRandomProvider rng) {
122            return delegate.createSampler(rng);
123        }
124    }
125
126    /**
127     * Specialisation of the trapezoidal distribution used when {@code b == c}.
128     *
129     * <p>This delegates all methods to the triangular distribution.
130     */
131    private static class TriangularTrapezoidalDistribution extends DelegatedTrapezoidalDistribution {
132        /**
133         * @param a Lower limit of this distribution (inclusive).
134         * @param b Start/end of the trapezoid constant density (mode).
135         * @param d Upper limit of this distribution (inclusive).
136         */
137        TriangularTrapezoidalDistribution(double a, double b, double d) {
138            super(a, b, b, d, TriangularDistribution.of(a, b, d));
139        }
140    }
141
142    /**
143     * Specialisation of the trapezoidal distribution used when {@code a == b} and {@code c == d}.
144     *
145     * <p>This delegates all methods to the uniform distribution.
146     */
147    private static class UniformTrapezoidalDistribution extends DelegatedTrapezoidalDistribution {
148        /**
149         * @param a Lower limit of this distribution (inclusive).
150         * @param d Upper limit of this distribution (inclusive).
151         */
152        UniformTrapezoidalDistribution(double a, double d) {
153            super(a, a, d, d, UniformContinuousDistribution.of(a, d));
154        }
155    }
156
157    /**
158     * Regular implementation of the trapezoidal distribution.
159     */
160    private static class RegularTrapezoidalDistribution extends TrapezoidalDistribution {
161        /** Cached value (d + c - a - b). */
162        private final double divisor;
163        /** Cached value (b - a). */
164        private final double bma;
165        /** Cached value (d - c). */
166        private final double dmc;
167        /** Cumulative probability at b. */
168        private final double cdfB;
169        /** Cumulative probability at c. */
170        private final double cdfC;
171        /** Survival probability at b. */
172        private final double sfB;
173        /** Survival probability at c. */
174        private final double sfC;
175
176        /**
177         * @param a Lower limit of this distribution (inclusive).
178         * @param b Start of the trapezoid constant density.
179         * @param c End of the trapezoid constant density.
180         * @param d Upper limit of this distribution (inclusive).
181         */
182        RegularTrapezoidalDistribution(double a, double b, double c, double d) {
183            super(a, b, c, d);
184
185            // Sum positive terms
186            divisor = (d - a) + (c - b);
187            bma = b - a;
188            dmc = d - c;
189
190            cdfB = bma / divisor;
191            sfB = 1 - cdfB;
192            sfC = dmc / divisor;
193            cdfC = 1 - sfC;
194        }
195
196        @Override
197        public double density(double x) {
198            // Note: x < a allows correct density where a == b
199            if (x < a) {
200                return 0;
201            }
202            if (x < b) {
203                final double divident = (x - a) / bma;
204                return 2 * (divident / divisor);
205            }
206            if (x < c) {
207                return 2 / divisor;
208            }
209            if (x < d) {
210                final double divident = (d - x) / dmc;
211                return 2 * (divident / divisor);
212            }
213            return 0;
214        }
215
216        @Override
217        public double cumulativeProbability(double x)  {
218            if (x <= a) {
219                return 0;
220            }
221            if (x < b) {
222                final double divident = (x - a) * (x - a) / bma;
223                return divident / divisor;
224            }
225            if (x < c) {
226                final double divident = 2 * x - b - a;
227                return divident / divisor;
228            }
229            if (x < d) {
230                final double divident = (d - x) * (d - x) / dmc;
231                return 1 - divident / divisor;
232            }
233            return 1;
234        }
235
236        @Override
237        public double survivalProbability(double x)  {
238            // By symmetry:
239            if (x <= a) {
240                return 1;
241            }
242            if (x < b) {
243                final double divident = (x - a) * (x - a) / bma;
244                return 1 - divident / divisor;
245            }
246            if (x < c) {
247                final double divident = 2 * x - b - a;
248                return 1 - divident / divisor;
249            }
250            if (x < d) {
251                final double divident = (d - x) * (d - x) / dmc;
252                return divident / divisor;
253            }
254            return 0;
255        }
256
257        @Override
258        public double inverseCumulativeProbability(double p) {
259            ArgumentUtils.checkProbability(p);
260            if (p == 0) {
261                return a;
262            }
263            if (p == 1) {
264                return d;
265            }
266            if (p < cdfB) {
267                return a + Math.sqrt(p * divisor * bma);
268            }
269            if (p < cdfC) {
270                return 0.5 * ((p * divisor) + a + b);
271            }
272            return d - Math.sqrt((1 - p) * divisor * dmc);
273        }
274
275        @Override
276        public double inverseSurvivalProbability(double p) {
277            // By symmetry:
278            ArgumentUtils.checkProbability(p);
279            if (p == 1) {
280                return a;
281            }
282            if (p == 0) {
283                return d;
284            }
285            if (p > sfB) {
286                return a + Math.sqrt((1 - p) * divisor * bma);
287            }
288            if (p > sfC) {
289                return 0.5 * (((1 - p) * divisor) + a + b);
290            }
291            return d - Math.sqrt(p * divisor * dmc);
292        }
293
294        @Override
295        public double getMean() {
296            // Compute using a standardized distribution
297            // b' = (b-a) / (d-a)
298            // c' = (c-a) / (d-a)
299            final double scale = d - a;
300            final double bp = bma / scale;
301            final double cp = (c - a) / scale;
302            return nonCentralMoment(1, bp, cp) * scale + a;
303        }
304
305        @Override
306        public double getVariance() {
307            // Compute using a standardized distribution
308            // b' = (b-a) / (d-a)
309            // c' = (c-a) / (d-a)
310            final double scale = d - a;
311            final double bp = bma / scale;
312            final double cp = (c - a) / scale;
313            final double mu = nonCentralMoment(1, bp, cp);
314            return (nonCentralMoment(2, bp, cp) - mu * mu) * scale * scale;
315        }
316
317        /**
318         * Compute the {@code k}-th non-central moment of the standardized trapezoidal
319         * distribution.
320         *
321         * <p>Shifting the distribution by scale {@code (d - a)} and location {@code a}
322         * creates a standardized trapezoidal distribution. This has a simplified
323         * non-central moment as {@code a = 0, d = 1, 0 <= b < c <= 1}.
324         * <pre>
325         *               2             1       ( 1 - c^(k+2)           )
326         * E[X^k] = ----------- -------------- ( ----------- - b^(k+1) )
327         *          (1 + c - b) (k + 1)(k + 2) (    1 - c              )
328         * </pre>
329         *
330         * <p>Simplification eliminates issues computing the moments when {@code a == b}
331         * or {@code c == d} in the original (non-standardized) distribution.
332         *
333         * @param k Moment to compute
334         * @param b Start of the trapezoid constant density (shape parameter in [0, 1]).
335         * @param c End of the trapezoid constant density (shape parameter in [0, 1]).
336         * @return the moment
337         */
338        private static double nonCentralMoment(int k, double b, double c) {
339            // As c -> 1 then (1 - c^(k+2)) loses precision
340            // 1 - x^y == -(x^y - 1)    [high precision powm1]
341            //         == -(exp(y * log(x)) - 1)
342            // Note: avoid log(1) using the limit:
343            // (1 - c^(k+2)) / (1-c) -> (k+2) as c -> 1
344            final double term1 = c == 1 ? k + 2 : Math.expm1((k + 2) * Math.log(c)) / (c - 1);
345            final double term2 = Math.pow(b, k + 1);
346            return 2 * ((term1 - term2) / (c - b + 1) / ((k + 1) * (k + 2)));
347        }
348    }
349
350    /**
351     * @param a Lower limit of this distribution (inclusive).
352     * @param b Start of the trapezoid constant density.
353     * @param c End of the trapezoid constant density.
354     * @param d Upper limit of this distribution (inclusive).
355     */
356    TrapezoidalDistribution(double a, double b, double c, double d) {
357        this.a = a;
358        this.b = b;
359        this.c = c;
360        this.d = d;
361    }
362
363    /**
364     * Creates a trapezoidal distribution.
365     *
366     * <p>The distribution density is represented as an up sloping line from
367     * {@code a} to {@code b}, constant from {@code b} to {@code c}, and then a down
368     * sloping line from {@code c} to {@code d}.
369     *
370     * @param a Lower limit of this distribution (inclusive).
371     * @param b Start of the trapezoid constant density (first shape parameter).
372     * @param c End of the trapezoid constant density (second shape parameter).
373     * @param d Upper limit of this distribution (inclusive).
374     * @return the distribution
375     * @throws IllegalArgumentException if {@code a >= d}, if {@code b < a}, if
376     * {@code c < b} or if {@code c > d}.
377     */
378    public static TrapezoidalDistribution of(double a, double b, double c, double d) {
379        if (a >= d) {
380            throw new DistributionException(DistributionException.INVALID_RANGE_LOW_GTE_HIGH,
381                                            a, d);
382        }
383        if (b < a) {
384            throw new DistributionException(DistributionException.TOO_SMALL,
385                                            b, a);
386        }
387        if (c < b) {
388            throw new DistributionException(DistributionException.TOO_SMALL,
389                                            c, b);
390        }
391        if (c > d) {
392            throw new DistributionException(DistributionException.TOO_LARGE,
393                                            c, d);
394        }
395        // For consistency, delegate to the appropriate simplified distribution.
396        // Note: Floating-point equality comparison is intentional.
397        if (b == c) {
398            return new TriangularTrapezoidalDistribution(a, b, d);
399        }
400        if (d - a == c - b) {
401            return new UniformTrapezoidalDistribution(a, d);
402        }
403        return new RegularTrapezoidalDistribution(a, b, c, d);
404    }
405
406    /**
407     * {@inheritDoc}
408     *
409     * <p>For lower limit \( a \), start of the density constant region \( b \),
410     * end of the density constant region \( c \) and upper limit \( d \), the
411     * mean is:
412     *
413     * <p>\[ \frac{1}{3(d+c-b-a)}\left(\frac{d^3-c^3}{d-c}-\frac{b^3-a^3}{b-a}\right) \]
414     */
415    @Override
416    public abstract double getMean();
417
418    /**
419     * {@inheritDoc}
420     *
421     * <p>For lower limit \( a \), start of the density constant region \( b \),
422     * end of the density constant region \( c \) and upper limit \( d \), the
423     * variance is:
424     *
425     * <p>\[ \frac{1}{6(d+c-b-a)}\left(\frac{d^4-c^4}{d-c}-\frac{b^4-a^4}{b-a}\right) - \mu^2 \]
426     *
427     * <p>where \( \mu \) is the mean.
428     */
429    @Override
430    public abstract double getVariance();
431
432    /**
433     * Gets the start of the constant region of the density function.
434     *
435     * <p>This is the first shape parameter {@code b} of the distribution.
436     *
437     * @return the first shape parameter {@code b}
438     */
439    public double getB() {
440        return b;
441    }
442
443    /**
444     * Gets the end of the constant region of the density function.
445     *
446     * <p>This is the second shape parameter {@code c} of the distribution.
447     *
448     * @return the second shape parameter {@code c}
449     */
450    public double getC() {
451        return c;
452    }
453
454    /**
455     * {@inheritDoc}
456     *
457     * <p>The lower bound of the support is equal to the lower limit parameter
458     * {@code a} of the distribution.
459     */
460    @Override
461    public double getSupportLowerBound() {
462        return a;
463    }
464
465    /**
466     * {@inheritDoc}
467     *
468     * <p>The upper bound of the support is equal to the upper limit parameter
469     * {@code d} of the distribution.
470     */
471    @Override
472    public double getSupportUpperBound() {
473        return d;
474    }
475}