Gaussian fit in C#

Question

In a project I'm working on I need to obtain a Gaussian fit from a set of points - needing mean and variance for some processing, and possibly an error degree (or accuracy level) to let me figure out if the set of points really have a normal distribution.

I've found this question

but it is limited to 3 points only - whereas I need a fit that can work with any number of points.

What I need is similar to the labview Gaussian Peak Fit

I have looked at mathdotnet and aforge.net (using both in the same project), but I haven't found anything.

Does anybody know any C# or (easily convertible) C/C++ or Java solutions?

Alternatively, I've been told that an iterative algorithm should be used - I could implement it by myself (if not too much math-complicated). Any idea about what I can use? I've read a lot of articles (on Wikipedia and others found via Google) but I haven't found any clear indication of a solution.

Wherewith

see Sinusoidal fitting classes for c# and replace "sinusoidal" by "gaussian" in it. – Instep

Answer 1

I've found a good implementation in ImageJ, a public domain image processing program, whose source code is available here

Converted to C# and adapted to my needs.

Thanks to you guys for your answers... not strictly related to the solution I found, but somehow I got there with your help :)

Wherewith

@ephraim sorry I don't think I have that code anymore, it was long ago – Wherewith

Answer 2

In Math.Net (nuget), you can do:

var real_σ = 0.5;
var real_μ = 0;

//Define gaussian function
var gaussian = new Func<double, double, double, double>((σ, μ, x) =>
{
    return Normal.PDF(μ, σ, x);
});

//Generate sample gaussian data
var data = Enumerable.Range(0, 41)
    .Select(x => -2 + x * 0.1)
    .Select(x => new { x, y = gaussian(real_σ, real_μ, x) });

var xs = data.Select(d => d.x).ToArray();
var ys = data.Select(d => d.y).ToArray();
var initialGuess_σ = 1;
var initialGuess_μ = 0;

var fit = Fit.Curve(xs, ys, gaussian, initialGuess_σ, initialGuess_μ);
//fit.Item1 = σ, fit.Item2 = μ

Exogenous

Answer 3

Here I show an example of how you can fit an arbitrary function with an arbitrary number of parameters with upper/lower bounds for each individual parameter. Just as RexCardan's example it is done using the MathNet library.

It is not very fast, but it works.

In order to fit a different function change the double gaussian(Vector<double> vectorArg) method. If you change the number of vectorArgs you also need to adjust:

1. The number of elements in lowerBound, upperBound and initialGuess in CurveFit.
2. Change the number of parameters of return z => f(new DenseVector(new[] { parameters[0], parameters[1], parameters[2], parameters[3], parameters[4], parameters[5], z }));
3. Change the number of parameters of t => f(new DenseVector(new[] { p[0], p[1], p[2], p[3], p[4], p[5], t }))

Example code for a double gaussian

using MathNet.Numerics;
using MathNet.Numerics.Distributions;
using MathNet.Numerics.LinearAlgebra;
using MathNet.Numerics.LinearAlgebra.Double;
using System;
using System.Linq;

static class GaussianFit
{
    /// <summary>
    /// Non-linear least square Gaussian curve fit to data.
    /// </summary>
    /// <param name="mu1">x position of the first Gaussian.</param>
    /// <param name="mu2">x position of the second Gaussian.</param>
    /// <returns>Array of the Gaussian profile.</returns>
    public Func<double, double> CurveFit(double[] xData, double[] yData, double mu1, double mu2)
    {
        //Define gaussian function
        double gaussian(Vector<double> vectorArg)
        {
            double x = vectorArg.Last();
            double y = 
                vectorArg[0] * Normal.PDF(vectorArg[1], vectorArg[2], x)
                + vectorArg[3] * Normal.PDF(vectorArg[4], vectorArg[5], x);
            return y;
        }

        var lowerBound = new DenseVector(new[] { 0, mu1 * 0.98, 0.05, 0, mu2 * 0.98, 0.05 });
        var upperBound = new DenseVector(new[] { 1e10, mu1 * 1.02, 0.3, 1e10, mu2 * 1.02, 0.3 });
        var initialGuess = new DenseVector(new[] { 1000, mu1, 0.2, 1000, mu2, 0.2 });

        Func<double, double> fit = CurveFuncConstrained(
            xData, yData, gaussian, lowerBound, upperBound, initialGuess
        );

        return fit;
    }

    /// <summary>
    /// Non-linear least-squares fitting the points (x,y) to an arbitrary function y : x -> f(p0, p1, p2, x),
    /// returning a function y' for the best fitting curve.
    /// </summary>
    public static Func<double, double> CurveFuncConstrained(
        double[] x,
        double[] y,
        Func<Vector<double>, double> f,
        Vector<double> lowerBound,
        Vector<double> upperBound,
        Vector<double> initialGuess,
        double gradientTolerance = 1e-5,
        double parameterTolerance = 1e-5,
        double functionProgressTolerance = 1e-5,
        int maxIterations = 1000
    )
    {
        var parameters = CurveConstrained(
            x, y, f,
            lowerBound, upperBound, initialGuess,
            gradientTolerance, parameterTolerance, functionProgressTolerance,
            maxIterations
        );
        return z => f(new DenseVector(new[] { parameters[0], parameters[1], parameters[2], parameters[3], parameters[4], parameters[5], z }));
    }

    /// <summary>
    /// Non-linear least-squares fitting the points (x,y) to an arbitrary function y : x -> f(p0, p1, p2, x),
    /// returning its best fitting parameter p0, p1 and p2.
    /// </summary>
    public static Vector<double> CurveConstrained(
        double[] x,
        double[] y,
        Func<Vector<double>, double> f,
        Vector<double> lowerBound,
        Vector<double> upperBound,
        Vector<double> initialGuess,
        double gradientTolerance = 1e-5,
        double parameterTolerance = 1e-5,
        double functionProgressTolerance = 1e-5,
        int maxIterations = 1000
    )
    {
        return FindMinimum.OfFunctionConstrained(
            (p) => Distance.Euclidean(
                Generate.Map(
                    x, 
                    t => f(new DenseVector(new[] { p[0], p[1], p[2], p[3], p[4], p[5], t }))
                    ), 
                y),
            lowerBound,
            upperBound,
            initialGuess,
            gradientTolerance,
            parameterTolerance,
            functionProgressTolerance,
            maxIterations
        );
    }

Example

To fit two Gaussians at x positions 10 and 20:

Func<double, double> fit = GaussianFit.Curvefit(x_data, y_data, 10, 20);

Strobe

Thanks man, old question, it's nice to see that they never expire :-) Nice work, I can't test it (no longer working on that app, not even on Windows) but regardless I think you deserve a thumb up for the effort ;-). By looking at the code, it looks like it's what I was looking for. – Wherewith

Can confirm it's working with the current version of MathNet.Numerics. With Bounded BFGS solver you don't really need to provide means for gaussians like in the given example, you can just init with the center X location of your data or anything else. This is data dependent, of course, maybe my data didn't have enough noise for the initialization to make a difference. – Bulk

I modified this so that it should work for one Gaussian and then tried to fit some of my data. Even when I give it a good starting point the result makes no sense. Is there a way to use BfgsMinimizer directly and see what's going on? – Bannockburn

@Bannockburn Did you also adjust the bounds/initial guesses for the standard deviation and the amplitude? I'm not sure regarding your question about BfgsMinimizer. – Strobe

Answer 4

Just calculate the mean and standard deviation of your sample, those are the only two parameters of a Gaussian distribution.

For "goodness of fit", you can do something like mean-square error of the CDF.

Sociality

Answer 5

Concerning answer 1 by Antonio Oct 18 '11 at 18:03:

I've found a good implementation in ImageJ, a public domain image processing program, whose source code is available here.

Unfortunately, the link is broken, however you can still find a snapshot on archive.org:

https://imagej.nih.gov/ij/developer/source/ij/measure/CurveFitter.java.html

(I would have posted this as a comment to answer 1, but as a new user I am apparently not allowed to do that.)

Ukko

Hedvig