All Downloads are FREE. Search and download functionalities are using the official Maven repository.

org.apache.commons.statistics.distribution.TrapezoidalDistribution Maven / Gradle / Ivy

/*
 * Licensed to the Apache Software Foundation (ASF) under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * The ASF licenses this file to You under the Apache License, Version 2.0
 * (the "License"); you may not use this file except in compliance with
 * the License.  You may obtain a copy of the License at
 *
 *      http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

package org.apache.commons.statistics.distribution;

import org.apache.commons.rng.UniformRandomProvider;

/**
 * Implementation of the trapezoidal distribution.
 *
 * 

The probability density function of \( X \) is: * *

\[ f(x; a, b, c, d) = \begin{cases} * \frac{2}{d+c-a-b}\frac{x-a}{b-a} & \text{for } a\le x \lt b \\ * \frac{2}{d+c-a-b} & \text{for } b\le x \lt c \\ * \frac{2}{d+c-a-b}\frac{d-x}{d-c} & \text{for } c\le x \le d * \end{cases} \] * *

for \( -\infty \lt a \le b \le c \le d \lt \infty \) and * \( x \in [a, d] \). * *

Note the special cases: *

    *
  • \( b = c \) is the triangular distribution *
  • \( a = b \) and \( c = d \) is the uniform distribution *
* * @see Trapezoidal distribution (Wikipedia) */ public abstract class TrapezoidalDistribution extends AbstractContinuousDistribution { /** Lower limit of this distribution (inclusive). */ protected final double a; /** Start of the trapezoid constant density. */ protected final double b; /** End of the trapezoid constant density. */ protected final double c; /** Upper limit of this distribution (inclusive). */ protected final double d; /** * Specialisation of the trapezoidal distribution used when the distribution simplifies * to an alternative distribution. */ private static class DelegatedTrapezoidalDistribution extends TrapezoidalDistribution { /** Distribution delegate. */ private final ContinuousDistribution delegate; /** * @param a Lower limit of this distribution (inclusive). * @param b Start of the trapezoid constant density. * @param c End of the trapezoid constant density. * @param d Upper limit of this distribution (inclusive). * @param delegate Distribution delegate. */ DelegatedTrapezoidalDistribution(double a, double b, double c, double d, ContinuousDistribution delegate) { super(a, b, c, d); this.delegate = delegate; } @Override public double density(double x) { return delegate.density(x); } @Override public double probability(double x0, double x1) { return delegate.probability(x0, x1); } @Override public double logDensity(double x) { return delegate.logDensity(x); } @Override public double cumulativeProbability(double x) { return delegate.cumulativeProbability(x); } @Override public double inverseCumulativeProbability(double p) { return delegate.inverseCumulativeProbability(p); } @Override public double survivalProbability(double x) { return delegate.survivalProbability(x); } @Override public double inverseSurvivalProbability(double p) { return delegate.inverseSurvivalProbability(p); } @Override public double getMean() { return delegate.getMean(); } @Override public double getVariance() { return delegate.getVariance(); } @Override public Sampler createSampler(UniformRandomProvider rng) { return delegate.createSampler(rng); } } /** * Specialisation of the trapezoidal distribution used when {@code b == c}. * *

This delegates all methods to the triangular distribution. */ private static class TriangularTrapezoidalDistribution extends DelegatedTrapezoidalDistribution { /** * @param a Lower limit of this distribution (inclusive). * @param b Start/end of the trapezoid constant density (mode). * @param d Upper limit of this distribution (inclusive). */ TriangularTrapezoidalDistribution(double a, double b, double d) { super(a, b, b, d, TriangularDistribution.of(a, b, d)); } } /** * Specialisation of the trapezoidal distribution used when {@code a == b} and {@code c == d}. * *

This delegates all methods to the uniform distribution. */ private static class UniformTrapezoidalDistribution extends DelegatedTrapezoidalDistribution { /** * @param a Lower limit of this distribution (inclusive). * @param d Upper limit of this distribution (inclusive). */ UniformTrapezoidalDistribution(double a, double d) { super(a, a, d, d, UniformContinuousDistribution.of(a, d)); } } /** * Regular implementation of the trapezoidal distribution. */ private static class RegularTrapezoidalDistribution extends TrapezoidalDistribution { /** Cached value (d + c - a - b). */ private final double divisor; /** Cached value (b - a). */ private final double bma; /** Cached value (d - c). */ private final double dmc; /** Cumulative probability at b. */ private final double cdfB; /** Cumulative probability at c. */ private final double cdfC; /** Survival probability at b. */ private final double sfB; /** Survival probability at c. */ private final double sfC; /** * @param a Lower limit of this distribution (inclusive). * @param b Start of the trapezoid constant density. * @param c End of the trapezoid constant density. * @param d Upper limit of this distribution (inclusive). */ RegularTrapezoidalDistribution(double a, double b, double c, double d) { super(a, b, c, d); // Sum positive terms divisor = (d - a) + (c - b); bma = b - a; dmc = d - c; cdfB = bma / divisor; sfB = 1 - cdfB; sfC = dmc / divisor; cdfC = 1 - sfC; } @Override public double density(double x) { // Note: x < a allows correct density where a == b if (x < a) { return 0; } if (x < b) { final double divident = (x - a) / bma; return 2 * (divident / divisor); } if (x < c) { return 2 / divisor; } if (x < d) { final double divident = (d - x) / dmc; return 2 * (divident / divisor); } return 0; } @Override public double cumulativeProbability(double x) { if (x <= a) { return 0; } if (x < b) { final double divident = (x - a) * (x - a) / bma; return divident / divisor; } if (x < c) { final double divident = 2 * x - b - a; return divident / divisor; } if (x < d) { final double divident = (d - x) * (d - x) / dmc; return 1 - divident / divisor; } return 1; } @Override public double survivalProbability(double x) { // By symmetry: if (x <= a) { return 1; } if (x < b) { final double divident = (x - a) * (x - a) / bma; return 1 - divident / divisor; } if (x < c) { final double divident = 2 * x - b - a; return 1 - divident / divisor; } if (x < d) { final double divident = (d - x) * (d - x) / dmc; return divident / divisor; } return 0; } @Override public double inverseCumulativeProbability(double p) { ArgumentUtils.checkProbability(p); if (p == 0) { return a; } if (p == 1) { return d; } if (p < cdfB) { return a + Math.sqrt(p * divisor * bma); } if (p < cdfC) { return 0.5 * ((p * divisor) + a + b); } return d - Math.sqrt((1 - p) * divisor * dmc); } @Override public double inverseSurvivalProbability(double p) { // By symmetry: ArgumentUtils.checkProbability(p); if (p == 1) { return a; } if (p == 0) { return d; } if (p > sfB) { return a + Math.sqrt((1 - p) * divisor * bma); } if (p > sfC) { return 0.5 * (((1 - p) * divisor) + a + b); } return d - Math.sqrt(p * divisor * dmc); } @Override public double getMean() { // Compute using a standardized distribution // b' = (b-a) / (d-a) // c' = (c-a) / (d-a) final double scale = d - a; final double bp = bma / scale; final double cp = (c - a) / scale; return nonCentralMoment(1, bp, cp) * scale + a; } @Override public double getVariance() { // Compute using a standardized distribution // b' = (b-a) / (d-a) // c' = (c-a) / (d-a) final double scale = d - a; final double bp = bma / scale; final double cp = (c - a) / scale; final double mu = nonCentralMoment(1, bp, cp); return (nonCentralMoment(2, bp, cp) - mu * mu) * scale * scale; } /** * Compute the {@code k}-th non-central moment of the standardized trapezoidal * distribution. * *

Shifting the distribution by scale {@code (d - a)} and location {@code a} * creates a standardized trapezoidal distribution. This has a simplified * non-central moment as {@code a = 0, d = 1, 0 <= b < c <= 1}. *

         *               2             1       ( 1 - c^(k+2)           )
         * E[X^k] = ----------- -------------- ( ----------- - b^(k+1) )
         *          (1 + c - b) (k + 1)(k + 2) (    1 - c              )
         * 
* *

Simplification eliminates issues computing the moments when {@code a == b} * or {@code c == d} in the original (non-standardized) distribution. * * @param k Moment to compute * @param b Start of the trapezoid constant density (shape parameter in [0, 1]). * @param c End of the trapezoid constant density (shape parameter in [0, 1]). * @return the moment */ private static double nonCentralMoment(int k, double b, double c) { // As c -> 1 then (1 - c^(k+2)) loses precision // 1 - x^y == -(x^y - 1) [high precision powm1] // == -(exp(y * log(x)) - 1) // Note: avoid log(1) using the limit: // (1 - c^(k+2)) / (1-c) -> (k+2) as c -> 1 final double term1 = c == 1 ? k + 2 : Math.expm1((k + 2) * Math.log(c)) / (c - 1); final double term2 = Math.pow(b, k + 1); return 2 * ((term1 - term2) / (c - b + 1) / ((k + 1) * (k + 2))); } } /** * @param a Lower limit of this distribution (inclusive). * @param b Start of the trapezoid constant density. * @param c End of the trapezoid constant density. * @param d Upper limit of this distribution (inclusive). */ TrapezoidalDistribution(double a, double b, double c, double d) { this.a = a; this.b = b; this.c = c; this.d = d; } /** * Creates a trapezoidal distribution. * *

The distribution density is represented as an up sloping line from * {@code a} to {@code b}, constant from {@code b} to {@code c}, and then a down * sloping line from {@code c} to {@code d}. * * @param a Lower limit of this distribution (inclusive). * @param b Start of the trapezoid constant density (first shape parameter). * @param c End of the trapezoid constant density (second shape parameter). * @param d Upper limit of this distribution (inclusive). * @return the distribution * @throws IllegalArgumentException if {@code a >= d}, if {@code b < a}, if * {@code c < b} or if {@code c > d}. */ public static TrapezoidalDistribution of(double a, double b, double c, double d) { if (a >= d) { throw new DistributionException(DistributionException.INVALID_RANGE_LOW_GTE_HIGH, a, d); } if (b < a) { throw new DistributionException(DistributionException.TOO_SMALL, b, a); } if (c < b) { throw new DistributionException(DistributionException.TOO_SMALL, c, b); } if (c > d) { throw new DistributionException(DistributionException.TOO_LARGE, c, d); } // For consistency, delegate to the appropriate simplified distribution. // Note: Floating-point equality comparison is intentional. if (b == c) { return new TriangularTrapezoidalDistribution(a, b, d); } if (d - a == c - b) { return new UniformTrapezoidalDistribution(a, d); } return new RegularTrapezoidalDistribution(a, b, c, d); } /** * {@inheritDoc} * *

For lower limit \( a \), start of the density constant region \( b \), * end of the density constant region \( c \) and upper limit \( d \), the * mean is: * *

\[ \frac{1}{3(d+c-b-a)}\left(\frac{d^3-c^3}{d-c}-\frac{b^3-a^3}{b-a}\right) \] */ @Override public abstract double getMean(); /** * {@inheritDoc} * *

For lower limit \( a \), start of the density constant region \( b \), * end of the density constant region \( c \) and upper limit \( d \), the * variance is: * *

\[ \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 \] * *

where \( \mu \) is the mean. */ @Override public abstract double getVariance(); /** * Gets the start of the constant region of the density function. * *

This is the first shape parameter {@code b} of the distribution. * * @return the first shape parameter {@code b} */ public double getB() { return b; } /** * Gets the end of the constant region of the density function. * *

This is the second shape parameter {@code c} of the distribution. * * @return the second shape parameter {@code c} */ public double getC() { return c; } /** * {@inheritDoc} * *

The lower bound of the support is equal to the lower limit parameter * {@code a} of the distribution. */ @Override public double getSupportLowerBound() { return a; } /** * {@inheritDoc} * *

The upper bound of the support is equal to the upper limit parameter * {@code d} of the distribution. */ @Override public double getSupportUpperBound() { return d; } }





© 2015 - 2024 Weber Informatics LLC | Privacy Policy