1 package nl.tudelft.simulation.jstats.distributions;
2
3 import org.djutils.exceptions.Throw;
4 import org.djutils.logger.CategoryLogger;
5
6 import nl.tudelft.simulation.jstats.math.ProbMath;
7 import nl.tudelft.simulation.jstats.streams.StreamInterface;
8
9 /**
10 * The Gamma distribution. For more information on this distribution see
11 * <a href="https://mathworld.wolfram.com/GammaDistribution.html"> https://mathworld.wolfram.com/GammaDistribution.html </a><br>
12 * The parameters are not rate-related, but average-related, so the mean is shape*scale (or αθ), and the variance is
13 * αθ<sup>2</sup>.
14 * <p>
15 * Copyright (c) 2002-2025 Delft University of Technology, Jaffalaan 5, 2628 BX Delft, the Netherlands. All rights reserved. See
16 * for project information <a href="https://simulation.tudelft.nl/dsol/manual/" target="_blank">DSOL Manual</a>. The DSOL
17 * project is distributed under a three-clause BSD-style license, which can be found at
18 * <a href="https://simulation.tudelft.nl/dsol/docs/latest/license.html" target="_blank">DSOL License</a>.
19 * </p>
20 * @author <a href="https://www.linkedin.com/in/peterhmjacobs">Peter Jacobs </a>
21 * @author <a href="https://github.com/averbraeck">Alexander Verbraeck</a>
22 */
23 public class DistGamma extends DistContinuous
24 {
25 /** the shape parameter of the distribution, also often called α or k. */
26 private final double shape;
27
28 /** the scale parameter of the distribution, also often called θ. */
29 private final double scale;
30
31 /**
32 * constructs a new gamma distribution. The gamma distribution represents the time to complete some task, e.g. customer
33 * service or machine repair. The parameters are not rate-related, but average-related, so the mean is shape*scale (or
34 * αθ or kθ), and the variance is αθ<sup>2</sup>.
35 * @param stream the random number stream
36 * @param shape is the shape parameter > 0, also known as α or k
37 * @param scale is the scale parameter> 0, also known as θ
38 * @throws IllegalArgumentException when shape <= 0.0 or scale <= 0.0
39 */
40 public DistGamma(final StreamInterface stream, final double shape, final double scale)
41 {
42 super(stream);
43 Throw.when(shape <= 0.0 || scale <= 0.0, IllegalArgumentException.class, "Error Gamma - shape <= 0.0 or scale <= 0.0");
44 this.shape = shape;
45 this.scale = scale;
46 }
47
48 @Override
49 public double draw()
50 {
51 // according to Law and Kelton, Simulation Modeling and Analysis, 1991
52 // pages 488-489
53 if (this.shape < 1.0)
54 {
55 double b = (Math.E + this.shape) / Math.E;
56 int counter = 0;
57 while (counter < 1000)
58 {
59 // step 1.
60 double p = b * this.stream.nextDouble();
61 if (p <= 1.0d)
62 {
63 // step 2.
64 double y = Math.pow(p, 1.0d / this.shape);
65 double u2 = this.stream.nextDouble();
66 if (u2 <= Math.exp(-y))
67 {
68 return this.scale * y;
69 }
70 }
71 else
72 {
73 // step 3.
74 double y = -Math.log((b - p) / this.shape);
75 double u2 = this.stream.nextDouble();
76 if (u2 <= Math.pow(y, this.shape - 1.0d))
77 {
78 return this.scale * y;
79 }
80 }
81 counter++;
82 }
83 CategoryLogger.always().info("Gamma distribution -- 1000 tries for alpha<1.0");
84 return 1.0d;
85 }
86 else if (this.shape > 1.0)
87 {
88 // according to Law and Kelton, Simulation Modeling and
89 // Analysis, 1991, pages 488-489
90 double a = 1.0d / Math.sqrt(2.0d * this.shape - 1.0d);
91 double b = this.shape - Math.log(4.0d);
92 double q = this.shape + (1.0d / a);
93 double theta = 4.5d;
94 double d = 1.0d + Math.log(theta);
95 int counter = 0;
96 while (counter < 1000)
97 {
98 // step 1.
99 double u1 = this.stream.nextDouble();
100 double u2 = this.stream.nextDouble();
101 // step 2.
102 double v = a * Math.log(u1 / (1.0d - u1));
103 double y = this.shape * Math.exp(v);
104 double z = u1 * u1 * u2;
105 double w = b + q * v - y;
106 // step 3.
107 if ((w + d - theta * z) >= 0.0d)
108 {
109 return this.scale * y;
110 }
111 // step 4.
112 if (w > Math.log(z))
113 {
114 return this.scale * y;
115 }
116 counter++;
117 }
118 CategoryLogger.always().info("Gamma distribution -- 1000 tries for alpha>1.0");
119 return 1.0d;
120 }
121 else
122 // shape == 1.0
123 {
124 // Gamma(1.0, scale) ~ exponential with mean = scale
125 return -this.scale * Math.log(this.stream.nextDouble());
126 }
127 }
128
129 @Override
130 public double getProbabilityDensity(final double x)
131 {
132 if (x <= 0)
133 {
134 return 0.0;
135 }
136 return (Math.pow(this.scale, -this.shape) * Math.pow(x, this.shape - 1) * Math.exp(-1 * x / this.scale))
137 / ProbMath.gamma(this.shape);
138 }
139
140 /**
141 * @return alpha
142 */
143 public double getShape()
144 {
145 return this.shape;
146 }
147
148 /**
149 * @return beta
150 */
151 public double getScale()
152 {
153 return this.scale;
154 }
155
156 @Override
157 public String toString()
158 {
159 return "Gamma(" + this.shape + "," + this.scale + ")";
160 }
161 }