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 > 0 \) the shape, 035 * \( \Omega > 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}