1 /*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * https://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17
18 /*
19 * This is not the original file distributed by the Apache Software Foundation
20 * It has been modified by the Hipparchus project
21 */
22 package org.hipparchus.distribution.continuous;
23
24 import org.hipparchus.exception.LocalizedCoreFormats;
25 import org.hipparchus.exception.MathIllegalArgumentException;
26 import org.hipparchus.special.Gamma;
27 import org.hipparchus.util.FastMath;
28
29 /**
30 * Implementation of the Gamma distribution.
31 *
32 * @see <a href="http://en.wikipedia.org/wiki/Gamma_distribution">Gamma distribution (Wikipedia)</a>
33 * @see <a href="http://mathworld.wolfram.com/GammaDistribution.html">Gamma distribution (MathWorld)</a>
34 */
35 public class GammaDistribution extends AbstractRealDistribution {
36 /** Serializable version identifier. */
37 private static final long serialVersionUID = 20120524L;
38 /** The shape parameter. */
39 private final double shape;
40 /** The scale parameter. */
41 private final double scale;
42 /**
43 * The constant value of {@code shape + g + 0.5}, where {@code g} is the
44 * Lanczos constant {@link Gamma#LANCZOS_G}.
45 */
46 private final double shiftedShape;
47 /**
48 * The constant value of
49 * {@code shape / scale * sqrt(e / (2 * pi * (shape + g + 0.5))) / L(shape)},
50 * where {@code L(shape)} is the Lanczos approximation returned by
51 * {@link Gamma#lanczos(double)}. This prefactor is used in
52 * {@link #density(double)}, when no overflow occurs with the natural
53 * calculation.
54 */
55 private final double densityPrefactor1;
56 /**
57 * The constant value of
58 * {@code log(shape / scale * sqrt(e / (2 * pi * (shape + g + 0.5))) / L(shape))},
59 * where {@code L(shape)} is the Lanczos approximation returned by
60 * {@link Gamma#lanczos(double)}. This prefactor is used in
61 * {@link #logDensity(double)}, when no overflow occurs with the natural
62 * calculation.
63 */
64 private final double logDensityPrefactor1;
65 /**
66 * The constant value of
67 * {@code shape * sqrt(e / (2 * pi * (shape + g + 0.5))) / L(shape)},
68 * where {@code L(shape)} is the Lanczos approximation returned by
69 * {@link Gamma#lanczos(double)}. This prefactor is used in
70 * {@link #density(double)}, when overflow occurs with the natural
71 * calculation.
72 */
73 private final double densityPrefactor2;
74 /**
75 * The constant value of
76 * {@code log(shape * sqrt(e / (2 * pi * (shape + g + 0.5))) / L(shape))},
77 * where {@code L(shape)} is the Lanczos approximation returned by
78 * {@link Gamma#lanczos(double)}. This prefactor is used in
79 * {@link #logDensity(double)}, when overflow occurs with the natural
80 * calculation.
81 */
82 private final double logDensityPrefactor2;
83 /**
84 * Lower bound on {@code y = x / scale} for the selection of the computation
85 * method in {@link #density(double)}. For {@code y <= minY}, the natural
86 * calculation overflows.
87 */
88 private final double minY;
89 /**
90 * Upper bound on {@code log(y)} ({@code y = x / scale}) for the selection
91 * of the computation method in {@link #density(double)}. For
92 * {@code log(y) >= maxLogY}, the natural calculation overflows.
93 */
94 private final double maxLogY;
95
96 /**
97 * Creates a new gamma distribution with specified values of the shape and
98 * scale parameters.
99 *
100 * @param shape the shape parameter
101 * @param scale the scale parameter
102 * @throws MathIllegalArgumentException if {@code shape <= 0} or
103 * {@code scale <= 0}.
104 */
105 public GammaDistribution(double shape, double scale) throws MathIllegalArgumentException {
106 this(shape, scale, DEFAULT_SOLVER_ABSOLUTE_ACCURACY);
107 }
108
109
110 /**
111 * Creates a Gamma distribution.
112 *
113 * @param shape the shape parameter
114 * @param scale the scale parameter
115 * @param inverseCumAccuracy the maximum absolute error in inverse
116 * cumulative probability estimates (defaults to
117 * {@link #DEFAULT_SOLVER_ABSOLUTE_ACCURACY}).
118 * @throws MathIllegalArgumentException if {@code shape <= 0} or
119 * {@code scale <= 0}.
120 */
121 public GammaDistribution(final double shape,
122 final double scale,
123 final double inverseCumAccuracy)
124 throws MathIllegalArgumentException {
125 super(inverseCumAccuracy);
126
127 if (shape <= 0) {
128 throw new MathIllegalArgumentException(LocalizedCoreFormats.SHAPE, shape);
129 }
130 if (scale <= 0) {
131 throw new MathIllegalArgumentException(LocalizedCoreFormats.SCALE, scale);
132 }
133
134 this.shape = shape;
135 this.scale = scale;
136 this.shiftedShape = shape + Gamma.LANCZOS_G + 0.5;
137 final double aux = FastMath.E / (2.0 * FastMath.PI * shiftedShape);
138 this.densityPrefactor2 = shape * FastMath.sqrt(aux) / Gamma.lanczos(shape);
139 this.logDensityPrefactor2 = FastMath.log(shape) + 0.5 * FastMath.log(aux) -
140 FastMath.log(Gamma.lanczos(shape));
141 this.densityPrefactor1 = this.densityPrefactor2 / scale *
142 FastMath.pow(shiftedShape, -shape) *
143 FastMath.exp(shape + Gamma.LANCZOS_G);
144 this.logDensityPrefactor1 = this.logDensityPrefactor2 - FastMath.log(scale) -
145 FastMath.log(shiftedShape) * shape +
146 shape + Gamma.LANCZOS_G;
147 this.minY = shape + Gamma.LANCZOS_G - FastMath.log(Double.MAX_VALUE);
148 this.maxLogY = FastMath.log(Double.MAX_VALUE) / (shape - 1.0);
149 }
150
151 /**
152 * Returns the shape parameter of {@code this} distribution.
153 *
154 * @return the shape parameter
155 */
156 public double getShape() {
157 return shape;
158 }
159
160 /**
161 * Returns the scale parameter of {@code this} distribution.
162 *
163 * @return the scale parameter
164 */
165 public double getScale() {
166 return scale;
167 }
168
169 /** {@inheritDoc} */
170 @Override
171 public double density(double x) {
172 /* The present method must return the value of
173 *
174 * 1 x a - x
175 * ---------- (-) exp(---)
176 * x Gamma(a) b b
177 *
178 * where a is the shape parameter, and b the scale parameter.
179 * Substituting the Lanczos approximation of Gamma(a) leads to the
180 * following expression of the density
181 *
182 * a e 1 y a
183 * - sqrt(------------------) ---- (-----------) exp(a - y + g),
184 * x 2 pi (a + g + 0.5) L(a) a + g + 0.5
185 *
186 * where y = x / b. The above formula is the "natural" computation, which
187 * is implemented when no overflow is likely to occur. If overflow occurs
188 * with the natural computation, the following identity is used. It is
189 * based on the BOOST library
190 * http://www.boost.org/doc/libs/1_35_0/libs/math/doc/sf_and_dist/html/math_toolkit/special/sf_gamma/igamma.html
191 * Formula (15) needs adaptations, which are detailed below.
192 *
193 * y a
194 * (-----------) exp(a - y + g)
195 * a + g + 0.5
196 * y - a - g - 0.5 y (g + 0.5)
197 * = exp(a log1pm(---------------) - ----------- + g),
198 * a + g + 0.5 a + g + 0.5
199 *
200 * where log1pm(z) = log(1 + z) - z. Therefore, the value to be
201 * returned is
202 *
203 * a e 1
204 * - sqrt(------------------) ----
205 * x 2 pi (a + g + 0.5) L(a)
206 * y - a - g - 0.5 y (g + 0.5)
207 * * exp(a log1pm(---------------) - ----------- + g).
208 * a + g + 0.5 a + g + 0.5
209 */
210 if (x < 0) {
211 return 0;
212 }
213 final double y = x / scale;
214 if ((y <= minY) || (FastMath.log(y) >= maxLogY)) {
215 /*
216 * Overflow.
217 */
218 final double aux1 = (y - shiftedShape) / shiftedShape;
219 final double aux2 = shape * (FastMath.log1p(aux1) - aux1);
220 final double aux3 = -y * (Gamma.LANCZOS_G + 0.5) / shiftedShape +
221 Gamma.LANCZOS_G + aux2;
222 return densityPrefactor2 / x * FastMath.exp(aux3);
223 }
224 /*
225 * Natural calculation.
226 */
227 return densityPrefactor1 * FastMath.exp(-y) * FastMath.pow(y, shape - 1);
228 }
229
230 /** {@inheritDoc} **/
231 @Override
232 public double logDensity(double x) {
233 /*
234 * see the comment in {@link #density(double)} for computation details
235 */
236 if (x < 0) {
237 return Double.NEGATIVE_INFINITY;
238 }
239 final double y = x / scale;
240 if ((y <= minY) || (FastMath.log(y) >= maxLogY)) {
241 /*
242 * Overflow.
243 */
244 final double aux1 = (y - shiftedShape) / shiftedShape;
245 final double aux2 = shape * (FastMath.log1p(aux1) - aux1);
246 final double aux3 = -y * (Gamma.LANCZOS_G + 0.5) / shiftedShape +
247 Gamma.LANCZOS_G + aux2;
248 return logDensityPrefactor2 - FastMath.log(x) + aux3;
249 }
250 /*
251 * Natural calculation.
252 */
253 return logDensityPrefactor1 - y + FastMath.log(y) * (shape - 1);
254 }
255
256 /**
257 * {@inheritDoc}
258 *
259 * The implementation of this method is based on:
260 * <ul>
261 * <li>
262 * <a href="http://mathworld.wolfram.com/Chi-SquaredDistribution.html">
263 * Chi-Squared Distribution</a>, equation (9).
264 * </li>
265 * <li>Casella, G., & Berger, R. (1990). <i>Statistical Inference</i>.
266 * Belmont, CA: Duxbury Press.
267 * </li>
268 * </ul>
269 */
270 @Override
271 public double cumulativeProbability(double x) {
272 double ret;
273
274 if (x <= 0) {
275 ret = 0;
276 } else {
277 ret = Gamma.regularizedGammaP(shape, x / scale);
278 }
279
280 return ret;
281 }
282
283 /**
284 * {@inheritDoc}
285 *
286 * For shape parameter {@code alpha} and scale parameter {@code beta}, the
287 * mean is {@code alpha * beta}.
288 */
289 @Override
290 public double getNumericalMean() {
291 return shape * scale;
292 }
293
294 /**
295 * {@inheritDoc}
296 *
297 * For shape parameter {@code alpha} and scale parameter {@code beta}, the
298 * variance is {@code alpha * beta^2}.
299 *
300 * @return {@inheritDoc}
301 */
302 @Override
303 public double getNumericalVariance() {
304 return shape * scale * scale;
305 }
306
307 /**
308 * {@inheritDoc}
309 *
310 * The lower bound of the support is always 0 no matter the parameters.
311 *
312 * @return lower bound of the support (always 0)
313 */
314 @Override
315 public double getSupportLowerBound() {
316 return 0;
317 }
318
319 /**
320 * {@inheritDoc}
321 *
322 * The upper bound of the support is always positive infinity
323 * no matter the parameters.
324 *
325 * @return upper bound of the support (always Double.POSITIVE_INFINITY)
326 */
327 @Override
328 public double getSupportUpperBound() {
329 return Double.POSITIVE_INFINITY;
330 }
331
332 /**
333 * {@inheritDoc}
334 *
335 * The support of this distribution is connected.
336 *
337 * @return {@code true}
338 */
339 @Override
340 public boolean isSupportConnected() {
341 return true;
342 }
343 }