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 /** */
26 private static final long serialVersionUID = 1L;
27
28 /** the shape parameter of the distribution, also often called α or k. */
29 private final double shape;
30
31 /** the scale parameter of the distribution, also often called θ. */
32 private final double scale;
33
34 /**
35 * constructs a new gamma distribution. The gamma distribution represents the time to complete some task, e.g. customer
36 * service or machine repair. The parameters are not rate-related, but average-related, so the mean is shape*scale (or
37 * αθ or kθ), and the variance is αθ<sup>2</sup>.
38 * @param stream StreamInterface; the random number stream
39 * @param shape double; is the shape parameter > 0, also known as α or k
40 * @param scale double; is the scale parameter> 0, also known as θ
41 * @throws IllegalArgumentException when shape <= 0.0 or scale <= 0.0
42 */
43 public DistGamma(final StreamInterface stream, final double shape, final double scale)
44 {
45 super(stream);
46 Throw.when(shape <= 0.0 || scale <= 0.0, IllegalArgumentException.class, "Error Gamma - shape <= 0.0 or scale <= 0.0");
47 this.shape = shape;
48 this.scale = scale;
49 }
50
51 @Override
52 public double draw()
53 {
54 // according to Law and Kelton, Simulation Modeling and Analysis, 1991
55 // pages 488-489
56 if (this.shape < 1.0)
57 {
58 double b = (Math.E + this.shape) / Math.E;
59 int counter = 0;
60 while (counter < 1000)
61 {
62 // step 1.
63 double p = b * this.stream.nextDouble();
64 if (p <= 1.0d)
65 {
66 // step 2.
67 double y = Math.pow(p, 1.0d / this.shape);
68 double u2 = this.stream.nextDouble();
69 if (u2 <= Math.exp(-y))
70 {
71 return this.scale * y;
72 }
73 }
74 else
75 {
76 // step 3.
77 double y = -Math.log((b - p) / this.shape);
78 double u2 = this.stream.nextDouble();
79 if (u2 <= Math.pow(y, this.shape - 1.0d))
80 {
81 return this.scale * y;
82 }
83 }
84 counter++;
85 }
86 CategoryLogger.always().info("Gamma distribution -- 1000 tries for alpha<1.0");
87 return 1.0d;
88 }
89 else if (this.shape > 1.0)
90 {
91 // according to Law and Kelton, Simulation Modeling and
92 // Analysis, 1991, pages 488-489
93 double a = 1.0d / Math.sqrt(2.0d * this.shape - 1.0d);
94 double b = this.shape - Math.log(4.0d);
95 double q = this.shape + (1.0d / a);
96 double theta = 4.5d;
97 double d = 1.0d + Math.log(theta);
98 int counter = 0;
99 while (counter < 1000)
100 {
101 // step 1.
102 double u1 = this.stream.nextDouble();
103 double u2 = this.stream.nextDouble();
104 // step 2.
105 double v = a * Math.log(u1 / (1.0d - u1));
106 double y = this.shape * Math.exp(v);
107 double z = u1 * u1 * u2;
108 double w = b + q * v - y;
109 // step 3.
110 if ((w + d - theta * z) >= 0.0d)
111 {
112 return this.scale * y;
113 }
114 // step 4.
115 if (w > Math.log(z))
116 {
117 return this.scale * y;
118 }
119 counter++;
120 }
121 CategoryLogger.always().info("Gamma distribution -- 1000 tries for alpha>1.0");
122 return 1.0d;
123 }
124 else
125 // shape == 1.0
126 {
127 // Gamma(1.0, scale) ~ exponential with mean = scale
128 return -this.scale * Math.log(this.stream.nextDouble());
129 }
130 }
131
132 @Override
133 public double getProbabilityDensity(final double x)
134 {
135 if (x <= 0)
136 {
137 return 0.0;
138 }
139 return (Math.pow(this.scale, -this.shape) * Math.pow(x, this.shape - 1) * Math.exp(-1 * x / this.scale))
140 / ProbMath.gamma(this.shape);
141 }
142
143 /**
144 * @return alpha
145 */
146 public double getShape()
147 {
148 return this.shape;
149 }
150
151 /**
152 * @return beta
153 */
154 public double getScale()
155 {
156 return this.scale;
157 }
158
159 @Override
160 public String toString()
161 {
162 return "Gamma(" + this.shape + "," + this.scale + ")";
163 }
164 }