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.linearalgebra.fitting; 22 23 import de.lmu.ifi.dbs.elki.math.MathUtil; 24 import net.jafama.FastMath; 25 26 /** 27 * Gaussian function for parameter fitting 28 * <p> 29 * Based loosely on fgauss in the book "Numerical Recipies".<br> 30 * We did not bother to implement all optimizations at the benefit of having 31 * easier to use parameters. Instead of position, amplitude and width used in 32 * the book, we use the traditional Gaussian parameters mean, standard deviation 33 * and a linear scaling factor (which is mostly useful when combining multiple 34 * distributions) The cost are some additional computations such as a square 35 * root and probably a slight loss in precision. This could of course have been 36 * handled by an appropriate wrapper instead. 37 * <p> 38 * Due to their license, we cannot use their code, but we have to implement the 39 * mathematics ourselves. We hope the loss in precision isn't big. 40 * <p> 41 * They are also arranged differently: the book uses 42 * <code>amplitude, position, width</code> whereas we use 43 * <code>mean, stddev, scaling</code>.<br> 44 * But we're obviously using essentially the same mathematics. 45 * <p> 46 * The function also can use a mixture of gaussians, just use an appropriate 47 * number of parameters (which obviously needs to be a multiple of 3) 48 * 49 * @author Erich Schubert 50 * @since 0.2 51 */ 52 public class GaussianFittingFunction implements FittingFunction { 53 /** 54 * Static instance 55 */ 56 public static final GaussianFittingFunction STATIC = new GaussianFittingFunction(); 57 58 /** 59 * Compute the mixture of Gaussians at the given position 60 */ 61 @Override eval(double x, double[] params)62 public FittingFunctionResult eval(double x, double[] params) { 63 final int len = params.length; 64 // We always need triples: (mean, stddev, scaling) 65 assert (len % 3) == 0; 66 67 double y = 0.0; 68 double[] gradients = new double[len]; 69 70 // Loosely based on the book: 71 // Numerical Recipes in C: The Art of Scientific Computing 72 // Due to their license, we cannot use their code, but we have to implement 73 // the mathematics ourselves. We hope the loss in precision is not too big. 74 for(int i = 2; i < params.length; i += 3) { 75 // Standardized Gaussian parameter (centered, scaled by stddev) 76 double stdpar = (x - params[i - 2]) / params[i - 1]; 77 double e = FastMath.exp(-.5 * stdpar * stdpar); 78 double localy = params[i] / (params[i - 1] * MathUtil.SQRTTWOPI) * e; 79 80 y += localy; 81 // mean gradient 82 gradients[i - 2] = localy * stdpar; 83 // stddev gradient 84 gradients[i - 1] = (stdpar * stdpar - 1.0) * localy; 85 // amplitude gradient 86 gradients[i] = e / (params[i - 1] * MathUtil.SQRTTWOPI); 87 } 88 89 return new FittingFunctionResult(y, gradients); 90 } 91 } 92