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
019import org.apache.commons.numbers.gamma.Gamma;
020import org.apache.commons.numbers.gamma.GammaRatio;
021import org.apache.commons.numbers.gamma.LogGamma;
022import org.apache.commons.numbers.gamma.RegularizedGamma;
023import org.apache.commons.rng.UniformRandomProvider;
024import org.apache.commons.rng.sampling.distribution.AhrensDieterMarsagliaTsangGammaSampler;
025import org.apache.commons.rng.sampling.distribution.SharedStateContinuousSampler;
026
027/**
028 * Implementation of the Nakagami distribution.
029 *
030 * <p>The probability density function of \( X \) is:
031 *
032 * <p>\[ f(x; \mu, \Omega) = \frac{2\mu^\mu}{\Gamma(\mu)\Omega^\mu}x^{2\mu-1}\exp\left(-\frac{\mu}{\Omega}x^2\right) \]
033 *
034 * <p>for \( \mu &gt; 0 \) the shape,
035 * \( \Omega &gt; 0 \) the scale, and
036 * \( x \in (0, \infty) \).
037 *
038 * @see <a href="https://en.wikipedia.org/wiki/Nakagami_distribution">Nakagami distribution (Wikipedia)</a>
039 */
040public final class NakagamiDistribution extends AbstractContinuousDistribution {
041    /** Support lower bound. */
042    private static final double SUPPORT_LO = 0;
043    /** Support upper bound. */
044    private static final double SUPPORT_HI = Double.POSITIVE_INFINITY;
045    /** Natural logarithm of 2. */
046    private static final double LN_2 = 0.6931471805599453094172321;
047
048    /** The shape parameter. */
049    private final double mu;
050    /** The scale parameter. */
051    private final double omega;
052    /** Density prefactor. */
053    private final double densityPrefactor;
054    /** Log density prefactor. */
055    private final double logDensityPrefactor;
056    /** Cached value for inverse probability function. */
057    private final double mean;
058    /** Cached value for inverse probability function. */
059    private final double variance;
060
061    /**
062     * @param mu Shape parameter (must be positive).
063     * @param omega Scale parameter (must be positive). Controls the spread of the distribution.
064     */
065    private NakagamiDistribution(double mu,
066                                 double omega) {
067        this.mu = mu;
068        this.omega = omega;
069        densityPrefactor = 2.0 * Math.pow(mu, mu) / (Gamma.value(mu) * Math.pow(omega, mu));
070        logDensityPrefactor = LN_2 + Math.log(mu) * mu - LogGamma.value(mu) - Math.log(omega) * mu;
071        final double v = GammaRatio.delta(mu, 0.5);
072        mean = Math.sqrt(omega / mu) / v;
073        variance = omega - (omega / mu) / v / v;
074    }
075
076    /**
077     * Creates a Nakagami distribution.
078     *
079     * @param mu Shape parameter (must be positive).
080     * @param omega Scale parameter (must be positive). Controls the spread of the distribution.
081     * @return the distribution
082     * @throws IllegalArgumentException  if {@code mu <= 0} or if
083     * {@code omega <= 0}.
084     */
085    public static NakagamiDistribution of(double mu,
086                                          double omega) {
087        if (mu <= 0) {
088            throw new DistributionException(DistributionException.NOT_STRICTLY_POSITIVE, mu);
089        }
090        if (omega <= 0) {
091            throw new DistributionException(DistributionException.NOT_STRICTLY_POSITIVE, omega);
092        }
093        return new NakagamiDistribution(mu, omega);
094    }
095
096    /**
097     * Gets the shape parameter of this distribution.
098     *
099     * @return the shape parameter.
100     */
101    public double getShape() {
102        return mu;
103    }
104
105    /**
106     * Gets the scale parameter of this distribution.
107     *
108     * @return the scale parameter.
109     */
110    public double getScale() {
111        return omega;
112    }
113
114    /** {@inheritDoc} */
115    @Override
116    public double density(double x) {
117        if (x <= SUPPORT_LO ||
118            x >= SUPPORT_HI) {
119            return 0;
120        }
121
122        return densityPrefactor * Math.pow(x, 2 * mu - 1) * Math.exp(-mu * x * x / omega);
123    }
124
125    /** {@inheritDoc} */
126    @Override
127    public double logDensity(double x) {
128        if (x <= SUPPORT_LO ||
129            x >= SUPPORT_HI) {
130            return Double.NEGATIVE_INFINITY;
131        }
132
133        return logDensityPrefactor + Math.log(x) * (2 * mu - 1) - (mu * x * x / omega);
134    }
135
136    /** {@inheritDoc} */
137    @Override
138    public double cumulativeProbability(double x) {
139        if (x <= SUPPORT_LO) {
140            return 0;
141        } else if (x >= SUPPORT_HI) {
142            return 1;
143        }
144
145        return RegularizedGamma.P.value(mu, mu * x * x / omega);
146    }
147
148    /** {@inheritDoc} */
149    @Override
150    public double survivalProbability(double x) {
151        if (x <= SUPPORT_LO) {
152            return 1;
153        } else if (x >= SUPPORT_HI) {
154            return 0;
155        }
156
157        return RegularizedGamma.Q.value(mu, mu * x * x / omega);
158    }
159
160    /**
161     * {@inheritDoc}
162     *
163     * <p>For shape parameter \( \mu \) and scale parameter \( \Omega \), the mean is:
164     *
165     * <p>\[ \frac{\Gamma(m+\frac{1}{2})}{\Gamma(m)}\left(\frac{\Omega}{m}\right)^{1/2} \]
166     */
167    @Override
168    public double getMean() {
169        return mean;
170    }
171
172    /**
173     * {@inheritDoc}
174     *
175     * <p>For shape parameter \( \mu \) and scale parameter \( \Omega \), the variance is:
176     *
177     * <p>\[ \Omega\left(1-\frac{1}{m}\left(\frac{\Gamma(m+\frac{1}{2})}{\Gamma(m)}\right)^2\right) \]
178     */
179    @Override
180    public double getVariance() {
181        return variance;
182    }
183
184    /**
185     * {@inheritDoc}
186     *
187     * <p>The lower bound of the support is always 0.
188     *
189     * @return 0.
190     */
191    @Override
192    public double getSupportLowerBound() {
193        return SUPPORT_LO;
194    }
195
196    /**
197     * {@inheritDoc}
198     *
199     * <p>The upper bound of the support is always positive infinity.
200     *
201     * @return {@link Double#POSITIVE_INFINITY positive infinity}.
202     */
203    @Override
204    public double getSupportUpperBound() {
205        return SUPPORT_HI;
206    }
207
208    @Override
209    public Sampler createSampler(UniformRandomProvider rng) {
210        // Generate using a related Gamma distribution
211        // See https://en.wikipedia.org/wiki/Nakagami_distribution#Generation
212        final double shape = mu;
213        final double scale = omega / mu;
214        final SharedStateContinuousSampler sampler =
215            AhrensDieterMarsagliaTsangGammaSampler.of(rng, shape, scale);
216        return () -> Math.sqrt(sampler.sample());
217    }
218}