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} & \text{for } a\le x \lt b \\ 029 * \frac{2}{d+c-a-b} & \text{for } b\le x \lt c \\ 030 * \frac{2}{d+c-a-b}\frac{d-x}{d-c} & \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}