Hi. > >However, as I've raised in > > https://issues.apache.org/jira/browse/MATH-512 > >I think that "GaussianFitter" needs some refactoring. > >Maybe you could try to see what changes would be necessary to deal with > >MATH-512 and such that the upgraded "GaussianFitter" will meet with your > >requirements at the same time. > > I went ahead and implemented a "Normal" GaussianFitter. Turns out there's > not that much difference in the resulting best fit parameters, and it seems > four parameters are better than three.
Maybe the data you are fitting are not from a "Normal" Gaussian... [You add one more parameter and get a better fit but I think that it does not necessarily mean that your data is indeed best represented by the sum of a constant and a Gaussian.] > I agree with your comments in 512. I can refactor the classes as suggested, > just let me know... Yes, please give it a try; however before you do it, I'm thinking of another improvement: Why not move the defintion of "ParametricGaussianFunction" to the "Gaussian" class (in package "function")? [I've done just that in the attached copy of "Gaussian.java".] Thinking about it, it seems that the "ParametricRealFunction" interface might be of a more general use than just in fitting, so I'd move it over to the package "analysis" (where other function interfaces are defined). > Also if commons-math is interested I can submit the classes for the "Normal" > GaussianFitter. I thought about combining the two and I think the > implementation will be much cleaner if there are two separate implementations. At first sight, I don't think so. Once refactored, the "GaussianFitter" could have 2 constructors, one would take a "Gaussian.Parametric" function (that would be fitting the "Normal" case) and the other would take a "FourParameterGaussianParametricFunction". But I still wonder if the use of the word "Gaussian" in the latter is really appropriate. I'd even say that such a function shouldn't be in CM at all, unless one is willing to accept the implementation of a "FiveParameterGaussianParametricFunction" (where we would add, say, a quadratic function) and a "SixParameterGaussian...", etc. > Maybe the one that reflects the Wikipedia definition should be called > "GaussianFitter" and the current one should be Called > "GaussianFitterWithHeightFactor"... > > or perhaps > > ThreeParameterGaussianFitter > FourParameterGaussianFitter... As explained above, I'd rather leave such fuzzy names to the user-code layer. Regards, Gilles
/* * Licensed to the Apache Software Foundation (ASF) under one or more * contributor license agreements. See the NOTICE file distributed with * this work for additional information regarding copyright ownership. * The ASF licenses this file to You under the Apache License, Version 2.0 * (the "License"); you may not use this file except in compliance with * the License. You may obtain a copy of the License at * * http://www.apache.org/licenses/LICENSE-2.0 * * Unless required by applicable law or agreed to in writing, software * distributed under the License is distributed on an "AS IS" BASIS, * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. * See the License for the specific language governing permissions and * limitations under the License. */ package org.apache.commons.math.analysis.function; import org.apache.commons.math.analysis.UnivariateRealFunction; import org.apache.commons.math.analysis.DifferentiableUnivariateRealFunction; import org.apache.commons.math.exception.NotStrictlyPositiveException; import org.apache.commons.math.exception.NullArgumentException; import org.apache.commons.math.exception.DimensionMismatchException; import org.apache.commons.math.util.FastMath; import org.apache.commons.math.optimization.fitting.ParametricRealFunction; /** * <a href="http://en.wikipedia.org/wiki/Gaussian_function";> * Gaussian</a> function. * * @version $Revision$ $Date$ * @since 3.0 */ public class Gaussian implements DifferentiableUnivariateRealFunction { /** Mean. */ private final double mean; /** Inverse of twice the square of the standard deviation. */ private final double i2s2; /** Normalization factor. */ private final double norm; /** * Gaussian with given normalization factor, mean and standard deviation. * * @param norm Normalization factor. * @param mean Mean. * @param sigma Standard deviation. * @throws NotStrictlyPositiveException if {@code sigma <= 0}. */ public Gaussian(double norm, double mean, double sigma) { if (sigma <= 0) { throw new NotStrictlyPositiveException(sigma); } this.norm = norm; this.mean = mean; this.i2s2 = 1 / (2 * sigma * sigma); } /** * Normalized gaussian with given mean and standard deviation. * * @param mean Mean. * @param sigma Standard deviation. * @throws NotStrictlyPositiveException if {@code sigma <= 0}. */ public Gaussian(double mean, double sigma) { this(1 / (sigma * FastMath.sqrt(2 * Math.PI)), mean, sigma); } /** * Normalized gaussian with zero mean and unit standard deviation. */ public Gaussian() { this(0, 1); } /** {@inheritDoc} */ public double value(double x) { return value(x - mean, norm, i2s2); } /** {@inheritDoc} */ public UnivariateRealFunction derivative() { return new UnivariateRealFunction() { /** {@inheritDoc} */ public double value(double x) { final double diff = x - mean; final double g = Gaussian.value(diff, norm, i2s2); if (g == 0) { // Avoid returning NaN in case of overflow. return 0; } else { return -2 * diff * i2s2 * g; } } }; } /** * Parametric function where the parameters array ordered as * follows: * <ul> * <li>Norm</li> * <li>Mean</li> * <li>Standard deviation</li> * </ul> */ public static class Parametric implements ParametricRealFunction { /** * Computes the value of the Gaussian at {@code x}. * * @param x Value for which the function must be computed. * @param param Values of norm, mean and standard deviation. * @return the value of the function. * @throws NullArgumentException if {@code param} is {@code null}. * @throws DimensionMismatchException if the size of {@code param} is * not 3. * @throws NotStrictlyPositiveException if {@code param[2]} is negative. */ public double value(double x, double[] param) { validateParameters(param); final double diff = x - param[1]; final double i2s2 = 1 / (2 * param[2] * param[2]); return Gaussian.value(diff, param[0], i2s2); } /** * Computes the value of the gradient at {@code x}. * The components of the gradient vector are the partial * derivatives of the function with respect to each of the * <em>parameters</em> (norm, mean and standard deviation). * * @param x Value at which the gradient must be computed. * @param param Values of norm, mean and standard deviation. * @return the gradient vector at {@code x}. * @throws NullArgumentException if {@code param} is {@code null}. * @throws DimensionMismatchException if the size of {@code param} is * not 3. * @throws NotStrictlyPositiveException if {@code param[2]} is negative. */ public double[] gradient(double x, double[] param) { validateParameters(param); final double diff = x - param[1]; final double norm = param[0]; final double sigma = param[2]; final double i2s2 = 1 / (2 * sigma * sigma); final double n = Gaussian.value(diff, 1, i2s2); final double m = norm * n * 2 * i2s2 * diff; final double s = m * diff / sigma; return new double[] { n, m, s }; } /** * Validates parameters to ensure they are appropriate for the evaluation of * the {@link #value(double,double[])} and {@link #gradient(double,double[])} * methods. * * @param param Values of norm, mean and standard deviation. * @throws NullArgumentException if {@code param} is {@code null}. * @throws DimensionMismatchException if the size of {@code param} is * not 3. * @throws NotStrictlyPositiveException if {@code param[2]} is negative. */ private void validateParameters(double[] param) { if (param == null) { throw new NullArgumentException(); } if (param.length != 3) { throw new DimensionMismatchException(param.length, 3); } if (param[2] <= 0) { throw new NotStrictlyPositiveException(param[2]); } } } /** * @param xMinusMean {@code x - mean}. * @param norm Normalization factor. * @param i2s2 Inverse of twice the square of the standard deviation. * @return the value of the Gaussian at {@code x}. */ private static double value(double xMinusMean, double norm, double i2s2) { return norm * FastMath.exp(-xMinusMean * xMinusMean * i2s2); } }
--------------------------------------------------------------------- To unsubscribe, e-mail: dev-unsubscr...@commons.apache.org For additional commands, e-mail: dev-h...@commons.apache.org
