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.

Answer 1

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.

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 = μ

Answer 3

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 :)

Answer 4

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);