Why am I getting a System.ArithmeticException in this example doing a non-linear regression to a Gaussian

Question

I am porting some code that I prototyped in Matlab to C# and need perform do a non-linear regression of a Gaussian to my sample data. I am currently trying Accord.Net and came up with the code below following the example here.

Unfortunately I am getting this exceptions and do not understand why:

System.ArithmeticException: Error calculation has produced a non-finite number. Please make sure that there are no constant columns in the input data."

with error message:

Error calculation has produced a non-finite number. Please make sure that there are no constant columns in the input data.

This is the code that I came up with:

    public static double[] FitGaussian()
    {
        var nls = CreateGaussianFitObject(new[] { 1.43e+04, 0.05093, 0.8098 }, 100, 0);

        double[,] data =
        {
            { -0.2,    12973.3071 },
            { -0.1,    13846.1569 },
            { -0,      14243.9094 },
            {  0.1,    14215.6044 },
            {  0.2,    13840.9077 },
        };

        // Extract inputs and outputs
        double[][] inputs = data.GetColumn(0).ToJagged();
        double[] outputs = data.GetColumn(1);

        var regression = nls.Learn(inputs, outputs);

        var result = regression.Coefficients;
        return result;
    }

    public static NonlinearLeastSquares CreateGaussianFitObject(double[] StartValues, int maxIterations, double tolerance)
    {
        Func<double, double, double, double> expFnc = (mu, sig, x) => m.Exp(-m.Pow((x - mu) / sig, 2));

        var nls = new NonlinearLeastSquares()
        {
            NumberOfParameters = 3,

            // Reference Gaussian Distributioon used in Matlab: f(x) =  a1*exp(-((x-b1)/c1)^2)
            //Function = (w, x) => w[0] * System.Math.Exp(-System.Math.Pow((x[0] - w[1]) / w[2], 2)),
            Function = (w, x) => w[0] * expFnc(w[1],w[2],x[0]),

            // Derivative in respect to the weights:
            Gradient = (w, x, r) =>
            {
                r[0] = expFnc(w[1], w[2], x[0]); // e^(-(x - b1)^2/c1^2)   <=>   diff a1*exp(-((x-b1)/c1)^2) w.r.t. a1
                r[1] = (2 * w[0] * (x[0] - w[1]) * expFnc(w[1], w[2], x[0])) / m.Pow(w[2],2); // (2 a1 (x - b1) e^(-(x - b1)^2/c1^2))/c1^2   <=>   diff a1*exp(-((x-b1)/c1)^2) w.r.t. b1
                r[2] = (2 * w[0] * m.Pow((x[0] - w[1]),2) * expFnc(w[1], w[2], x[0])) / m.Pow(w[2],3); // (2 a1 (x - b1)^2 e^(-(x - b1)^2/c1^2))/c1^3   <=>   diff a1*exp(-((x-b1)/c1)^2) w.r.t. c1
            },

            Algorithm = new LevenbergMarquardt()
            {
                MaxIterations = maxIterations,
                Tolerance = tolerance
            }
        };

        return nls;
    }

I have double-checked the gradient function (which I suspect as being the culprit) and it seems OK.

The corresponding Matlab code is this:

x = [-0.2,-0.1,0,0.1,0.2];
y = [12973.3071,13846.1569,14243.9094,14215.6044,13840.9077];
f = fit(x.',y.','gauss1');
fitObject = f;

which gives this output:

f = 

     General model Gauss1:
     f(x) =  a1*exp(-((x-b1)/c1)^2)
     Coefficients (with 95% confidence bounds):
       a1 =    1.43e+04  (1.419e+04, 1.441e+04)
       b1 =     0.05093  (0.03501, 0.06686)
       c1 =      0.8098  (0.7263, 0.8932)

As you can see I even tried giving the expected coefficients from Matlab as starting value to the C# implementation.

I am furthermore confused by the example that I cited above, because it seems to me they switched the coefficients and the inputs(?).

Answer 1

I forgot to use the argument StartValues in the constructor NonlinearLeastSquares(), which caused the exception. It should read:

    var nls = new NonlinearLeastSquares()
    {
        NumberOfParameters = 3,
        StartValues = StartValues,
        // Reference Gaussian Distributioon used in Matlab: f(x) =  a1*exp(-((x-b1)/c1)^2)
        //Function = (w, x) => w[0] * System.Math.Exp(-System.Math.Pow((x[0] - w[1]) / w[2], 2)),
        Function = (w, x) => w[0] * expFnc(w[1],w[2],x[0]),

        // Derivative in respect to the weights:
        Gradient = (w, x, r) =>
        {
            r[0] = expFnc(w[1], w[2], x[0]); // e^(-(x - b1)^2/c1^2)   <=>   diff a1*exp(-((x-b1)/c1)^2) w.r.t. a1
            r[1] = (2 * w[0] * (x[0] - w[1]) * expFnc(w[1], w[2], x[0])) / m.Pow(w[2],2); // (2 a1 (x - b1) e^(-(x - b1)^2/c1^2))/c1^2   <=>   diff a1*exp(-((x-b1)/c1)^2) w.r.t. b1
            r[2] = (2 * w[0] * m.Pow((x[0] - w[1]),2) * expFnc(w[1], w[2], x[0])) / m.Pow(w[2],3); // (2 a1 (x - b1)^2 e^(-(x - b1)^2/c1^2))/c1^3   <=>   diff a1*exp(-((x-b1)/c1)^2) w.r.t. c1
        },

        Algorithm = new LevenbergMarquardt()
        {
            MaxIterations = maxIterations,
            Tolerance = tolerance
        }
    };

Apparently StartValues is set to an array of zeros, if it is not explicitly set. This caused a division by zero and thus the exception. I found this out by adding

Debug.WriteLine($"Weights: w[0]: {w[0]}, w[1]: {w[1]}, w[2]: {w[2]}");
Debug.WriteLine($"Result: r[0]: {r[0]}, r[1]: {r[1]}, r[2]: {r[2]}");

to the lambda-expression for the Gradient property.