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 */
017package org.apache.commons.statistics.distribution;
018
019/**
020 * Implementation of the Gumbel distribution.
021 *
022 * <p>The probability density function of \( X \) is:
023 *
024 * <p>\[ f(x; \mu, \beta) =  \frac{1}{\beta} e^{-(z+e^{-z})} \]
025 *
026 * <p>where \[ z = \frac{x - \mu}{\beta} \]
027 *
028 * <p>for \( \mu \) the location,
029 * \( \beta &gt; 0 \) the scale, and
030 * \( x \in (-\infty, \infty) \).
031 *
032 * @see <a href="https://en.wikipedia.org/wiki/Gumbel_distribution">Gumbel distribution (Wikipedia)</a>
033 * @see <a href="https://mathworld.wolfram.com/GumbelDistribution.html">Gumbel distribution (MathWorld)</a>
034 */
035public final class GumbelDistribution extends AbstractContinuousDistribution {
036    /** Support lower bound. */
037    private static final double SUPPORT_LO = Double.NEGATIVE_INFINITY;
038    /** Support upper bound. */
039    private static final double SUPPORT_HI = Double.POSITIVE_INFINITY;
040    /** &pi;<sup>2</sup>/6. */
041    private static final double PI_SQUARED_OVER_SIX = Math.PI * Math.PI / 6;
042    /**
043     * <a href="https://en.wikipedia.org/wiki/Euler%27s_constant">
044     * Approximation of Euler's constant</a>.
045     */
046    private static final double EULER = 0.57721566490153286060;
047    /** ln(ln(2)). */
048    private static final double LN_LN_2 = -0.3665129205816643270124;
049    /** Location parameter. */
050    private final double mu;
051    /** Scale parameter. */
052    private final double beta;
053
054    /**
055     * @param mu Location parameter.
056     * @param beta Scale parameter (must be positive).
057     */
058    private GumbelDistribution(double mu,
059                               double beta) {
060        this.beta = beta;
061        this.mu = mu;
062    }
063
064    /**
065     * Creates a Gumbel distribution.
066     *
067     * @param mu Location parameter.
068     * @param beta Scale parameter (must be positive).
069     * @return the distribution
070     * @throws IllegalArgumentException if {@code beta <= 0}
071     */
072    public static GumbelDistribution of(double mu,
073                                        double beta) {
074        if (beta <= 0) {
075            throw new DistributionException(DistributionException.NOT_STRICTLY_POSITIVE, beta);
076        }
077        return new GumbelDistribution(mu, beta);
078    }
079
080    /**
081     * Gets the location parameter of this distribution.
082     *
083     * @return the location parameter.
084     */
085    public double getLocation() {
086        return mu;
087    }
088
089    /**
090     * Gets the scale parameter of this distribution.
091     *
092     * @return the scale parameter.
093     */
094    public double getScale() {
095        return beta;
096    }
097
098    /** {@inheritDoc} */
099    @Override
100    public double density(double x) {
101        if (x <= SUPPORT_LO) {
102            return 0;
103        }
104
105        final double z = (x - mu) / beta;
106        final double t = Math.exp(-z);
107        return Math.exp(-z - t) / beta;
108    }
109
110    /** {@inheritDoc} */
111    @Override
112    public double logDensity(double x) {
113        if (x <= SUPPORT_LO) {
114            return Double.NEGATIVE_INFINITY;
115        }
116
117        final double z = (x - mu) / beta;
118        final double t = Math.exp(-z);
119        return -z - t - Math.log(beta);
120    }
121
122    /** {@inheritDoc} */
123    @Override
124    public double cumulativeProbability(double x) {
125        final double z = (x - mu) / beta;
126        return Math.exp(-Math.exp(-z));
127    }
128
129    /** {@inheritDoc} */
130    @Override
131    public double survivalProbability(double x) {
132        final double z = (x - mu) / beta;
133        return -Math.expm1(-Math.exp(-z));
134    }
135
136    /** {@inheritDoc} */
137    @Override
138    public double inverseCumulativeProbability(double p) {
139        ArgumentUtils.checkProbability(p);
140        if (p == 0) {
141            return Double.NEGATIVE_INFINITY;
142        } else if (p == 1) {
143            return Double.POSITIVE_INFINITY;
144        }
145        return mu - Math.log(-Math.log(p)) * beta;
146    }
147
148    /** {@inheritDoc} */
149    @Override
150    public double inverseSurvivalProbability(double p) {
151        ArgumentUtils.checkProbability(p);
152        if (p == 1) {
153            return Double.NEGATIVE_INFINITY;
154        } else if (p == 0) {
155            return Double.POSITIVE_INFINITY;
156        }
157        return mu - Math.log(-Math.log1p(-p)) * beta;
158    }
159
160    /**
161     * {@inheritDoc}
162     *
163     * <p>For location parameter \( \mu \) and scale parameter \( \beta \), the mean is:
164     *
165     * <p>\[ \mu + \beta \gamma \]
166     *
167     * <p>where \( \gamma \) is the
168     * <a href="https://mathworld.wolfram.com/Euler-MascheroniConstantApproximations.html">
169     * Euler-Mascheroni constant</a>.
170     */
171    @Override
172    public double getMean() {
173        return mu + EULER * beta;
174    }
175
176    /**
177     * {@inheritDoc}
178     *
179     * <p>For scale parameter \( \beta \), the variance is:
180     *
181     * <p>\[ \frac{\pi^2}{6} \beta^2 \]
182     */
183    @Override
184    public double getVariance() {
185        return PI_SQUARED_OVER_SIX * beta * beta;
186    }
187
188    /**
189     * {@inheritDoc}
190     *
191     * <p>The lower bound of the support is always negative infinity.
192     *
193     * @return {@link Double#NEGATIVE_INFINITY negative infinity}.
194     */
195    @Override
196    public double getSupportLowerBound() {
197        return SUPPORT_LO;
198    }
199
200    /**
201     * {@inheritDoc}
202     *
203     * <p>The upper bound of the support is always positive infinity.
204     *
205     * @return {@link Double#POSITIVE_INFINITY positive infinity}.
206     */
207    @Override
208    public double getSupportUpperBound() {
209        return SUPPORT_HI;
210    }
211
212    /** {@inheritDoc} */
213    @Override
214    double getMedian() {
215        // Overridden for the probability(double, double) method.
216        // This is intentionally not a public method.
217        // u - beta * ln(ln(2))
218        return mu - beta * LN_LN_2;
219    }
220}