1 /* 2 * This file is part of ELKI: 3 * Environment for Developing KDD-Applications Supported by Index-Structures 4 * 5 * Copyright (C) 2018 6 * ELKI Development Team 7 * 8 * This program is free software: you can redistribute it and/or modify 9 * it under the terms of the GNU Affero General Public License as published by 10 * the Free Software Foundation, either version 3 of the License, or 11 * (at your option) any later version. 12 * 13 * This program is distributed in the hope that it will be useful, 14 * but WITHOUT ANY WARRANTY; without even the implied warranty of 15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 16 * GNU Affero General Public License for more details. 17 * 18 * You should have received a copy of the GNU Affero General Public License 19 * along with this program. If not, see <http://www.gnu.org/licenses/>. 20 */ 21 package de.lmu.ifi.dbs.elki.math.statistics.distribution; 22 23 import java.util.Random; 24 25 import de.lmu.ifi.dbs.elki.math.MathUtil; 26 import de.lmu.ifi.dbs.elki.utilities.documentation.Reference; 27 import de.lmu.ifi.dbs.elki.utilities.optionhandling.OptionID; 28 import de.lmu.ifi.dbs.elki.utilities.optionhandling.constraints.CommonConstraints; 29 import de.lmu.ifi.dbs.elki.utilities.optionhandling.parameterization.Parameterization; 30 import de.lmu.ifi.dbs.elki.utilities.optionhandling.parameters.DoubleParameter; 31 import de.lmu.ifi.dbs.elki.utilities.random.RandomFactory; 32 33 import net.jafama.FastMath; 34 35 /** 36 * Generalized normal distribution by adding a skew term, similar to lognormal 37 * distributions. 38 * <p> 39 * This is one kind of generalized normal distributions. Note that there are 40 * multiple that go by the name of a "Generalized Normal Distribution"; this is 41 * what is currently called "version 2" in English Wikipedia. 42 * <p> 43 * Reference: 44 * <p> 45 * J. R. M. Hosking, J. R. Wallis<br> 46 * Regional frequency analysis: an approach based on L-moments 47 * 48 * @author Erich Schubert 49 * @since 0.6.0 50 */ 51 @Reference(authors = "J. R. M. Hosking, J. R. Wallis", // 52 title = "Regional frequency analysis: an approach based on L-moments", // 53 booktitle = "Regional frequency analysis: an approach based on L-moments", // 54 url = "https://doi.org/10.1017/CBO9780511529443", // 55 bibkey = "doi:10.1017/CBO9780511529443") 56 public class SkewGeneralizedNormalDistribution extends AbstractDistribution { 57 /** 58 * Location. 59 */ 60 private double loc; 61 62 /** 63 * Scale. 64 */ 65 private double scale; 66 67 /** 68 * Skew. 69 */ 70 private double skew; 71 72 /** 73 * Constructor for Gaussian distribution 74 * 75 * @param loc Location 76 * @param scale Scale 77 * @param skew Skew 78 * @param random Random generator 79 */ SkewGeneralizedNormalDistribution(double loc, double scale, double skew, Random random)80 public SkewGeneralizedNormalDistribution(double loc, double scale, double skew, Random random) { 81 super(random); 82 this.loc = loc; 83 this.scale = scale; 84 this.skew = skew; 85 } 86 87 /** 88 * Constructor for Gaussian distribution 89 * 90 * @param loc Location 91 * @param scale Scale 92 * @param skew Skew 93 * @param random Random generator 94 */ SkewGeneralizedNormalDistribution(double loc, double scale, double skew, RandomFactory random)95 public SkewGeneralizedNormalDistribution(double loc, double scale, double skew, RandomFactory random) { 96 super(random); 97 this.loc = loc; 98 this.scale = scale; 99 this.skew = skew; 100 } 101 102 /** 103 * Constructor for Gaussian distribution 104 * 105 * @param loc Location 106 * @param scale Scale 107 * @param skew Skew 108 */ SkewGeneralizedNormalDistribution(double loc, double scale, double skew)109 public SkewGeneralizedNormalDistribution(double loc, double scale, double skew) { 110 this(loc, scale, skew, (Random) null); 111 } 112 113 /** 114 * Get the location parameter 115 * 116 * @return Location 117 */ getLocation()118 public double getLocation() { 119 return loc; 120 } 121 122 /** 123 * Get the scale parameter 124 * 125 * @return Scale 126 */ getScale()127 public double getScale() { 128 return scale; 129 } 130 131 /** 132 * Get the skew parameter. 133 * 134 * @return Skew parameter 135 */ getSkew()136 public double getSkew() { 137 return skew; 138 } 139 140 @Override pdf(double val)141 public double pdf(double val) { 142 return pdf(val, loc, scale, skew); 143 } 144 145 @Override logpdf(double val)146 public double logpdf(double val) { 147 return logpdf(val, loc, scale, skew); 148 } 149 150 @Override cdf(double val)151 public double cdf(double val) { 152 return cdf(val, loc, scale, skew); 153 } 154 155 @Override quantile(double q)156 public double quantile(double q) { 157 return quantile(q, loc, scale, skew); 158 } 159 160 @Override nextRandom()161 public double nextRandom() { 162 double y = random.nextGaussian(); 163 if(Math.abs(skew) > 0.) { 164 y = (1. - FastMath.exp(-skew * y)) / skew; 165 } 166 return loc + scale * y; 167 168 } 169 170 @Override toString()171 public String toString() { 172 return "SkewNormalDistribution(mean=" + loc + ", stddev=" + scale + ", skew=" + skew + ")"; 173 } 174 175 /** 176 * Probability density function of the skewed normal distribution. 177 * 178 * @param x The value. 179 * @param mu The mean. 180 * @param sigma The standard deviation. 181 * @return PDF of the given normal distribution at x. 182 */ pdf(double x, double mu, double sigma, double skew)183 public static double pdf(double x, double mu, double sigma, double skew) { 184 if(x != x) { 185 return Double.NaN; 186 } 187 x = (x - mu) / sigma; // Scale 188 if(skew == 0.) { 189 return MathUtil.ONE_BY_SQRTTWOPI / sigma * FastMath.exp(-.5 * x * x); 190 } 191 final double y = -FastMath.log1p(-skew * x) / skew; 192 if(y != y || y == Double.POSITIVE_INFINITY || y == Double.NEGATIVE_INFINITY) { // NaN 193 return 0.; 194 } 195 return MathUtil.ONE_BY_SQRTTWOPI / sigma * FastMath.exp(-.5 * y * y) / (1 - skew * x); 196 } 197 198 /** 199 * Probability density function of the skewed normal distribution. 200 * 201 * @param x The value. 202 * @param mu The mean. 203 * @param sigma The standard deviation. 204 * @return log PDF of the given normal distribution at x. 205 */ logpdf(double x, double mu, double sigma, double skew)206 public static double logpdf(double x, double mu, double sigma, double skew) { 207 if(x != x) { 208 return Double.NaN; 209 } 210 x = (x - mu) / sigma; 211 if(skew == 0.) { 212 return MathUtil.LOG_ONE_BY_SQRTTWOPI - FastMath.log(sigma) - .5 * x * x; 213 } 214 double y = -FastMath.log(1. - skew * x) / skew; 215 if(y != y || y == Double.NEGATIVE_INFINITY || y == Double.NEGATIVE_INFINITY) { 216 return Double.NEGATIVE_INFINITY; 217 } 218 return -.5 * y * y - FastMath.log(MathUtil.ONE_BY_SQRTTWOPI * sigma * (1 - skew * x)); 219 } 220 221 /** 222 * Cumulative probability density function (CDF) of a normal distribution. 223 * 224 * @param x value to evaluate CDF at 225 * @param mu Mean value 226 * @param sigma Standard deviation. 227 * @return The CDF of the given normal distribution at x. 228 */ cdf(double x, double mu, double sigma, double skew)229 public static double cdf(double x, double mu, double sigma, double skew) { 230 x = (x - mu) / sigma; 231 if(Math.abs(skew) > 0.) { 232 double tmp = 1 - skew * x; 233 if(tmp < 1e-15) { 234 return (skew < 0.) ? 0. : 1.; 235 } 236 x = -FastMath.log(tmp) / skew; 237 } 238 return .5 + .5 * NormalDistribution.erf(x * MathUtil.SQRTHALF); 239 } 240 241 /** 242 * Inverse cumulative probability density function (probit) of a normal 243 * distribution. 244 * 245 * @param x value to evaluate probit function at 246 * @param mu Mean value 247 * @param sigma Standard deviation. 248 * @return The probit of the given normal distribution at x. 249 */ quantile(double x, double mu, double sigma, double skew)250 public static double quantile(double x, double mu, double sigma, double skew) { 251 x = NormalDistribution.standardNormalQuantile(x); 252 if(Math.abs(skew) > 0.) { 253 x = (1. - FastMath.exp(-skew * x)) / skew; 254 } 255 return mu + sigma * x; 256 } 257 258 /** 259 * Parameterization class 260 * 261 * @author Erich Schubert 262 */ 263 public static class Parameterizer extends AbstractDistribution.Parameterizer { 264 /** 265 * Skew parameter 266 */ 267 public static final OptionID SKEW_ID = new OptionID("distribution.skewgnormal.skew", "Skew of the distribution."); 268 269 /** Parameters */ 270 double mean, sigma, skew; 271 272 @Override makeOptions(Parameterization config)273 protected void makeOptions(Parameterization config) { 274 super.makeOptions(config); 275 276 DoubleParameter meanP = new DoubleParameter(LOCATION_ID); 277 if(config.grab(meanP)) { 278 mean = meanP.doubleValue(); 279 } 280 281 DoubleParameter sigmaP = new DoubleParameter(SCALE_ID) // 282 .addConstraint(CommonConstraints.GREATER_THAN_ZERO_DOUBLE); 283 if(config.grab(sigmaP)) { 284 sigma = sigmaP.doubleValue(); 285 } 286 287 DoubleParameter skewP = new DoubleParameter(SKEW_ID); 288 if(config.grab(skewP)) { 289 skew = skewP.doubleValue(); 290 } 291 } 292 293 @Override makeInstance()294 protected SkewGeneralizedNormalDistribution makeInstance() { 295 return new SkewGeneralizedNormalDistribution(mean, sigma, skew, rnd); 296 } 297 } 298 } 299