This project has moved and is read-only. For the latest updates, please go here.

Problem to define the user-defined function (FitToFunction)

Aug 3, 2012 at 1:48 PM
Edited Aug 3, 2012 at 2:00 PM

I'm trying to fitting uncertain mesurements with an user-defined function (FitToFunction method). I'm using C# with VisualStudio11 Beta (Any CPU).

To start, i would like to check if the exemple found in this website works :

UncertainMeasurementSample data = new UncertainMeasurementSample();

data.Add(new UncertainMeasurement<double>(1.0, new UncertainValue(2.0,3.0)));

data.Add(new UncertainMeasurement<double>(2.0, 1.0, 3.0));

data.Add(3.0, 2.0, 1.0);

Function<double[], double, double> model = delegate (double[] a, double x) { double y = a[0] * Math.Sin(2.0 * Math.PI / a[1] + a[2]); return(y); };

FitResult oscilation = data.FitToFunction(model, new double[] { 1.0, 1.0, 0.0 });

First I got a problem with the declaration of Function : I have to replace "Function" by "Func"

After compilation and execution of code I got this error (english translation) :

System.NullReferenceException: The reference of this object is not defined to one instance of one object from Meta.Numerics.Statistics.UncertainMeasurementSample`1.FitToFunction(Func`3 function, Double[] start)

Any ideas ?

Collins

Aug 4, 2012 at 9:27 AM

Hi Collins,

With release 2.0, Meta.Numerics changed from using its own Function delegate to the .NET-standard Func delegate. Clearly we need to update the tutorial documentation. That is, as you pointed out, an easy fix.

The NullReferenceException is more troublesome. I can reproduce your result and see what's going on in the code: once we have found the parameter values that minimize the chi-squared function and determined the curvature matrix of the chi-squared function at the minimum, we try to invert the curvature matrix in order to determine the covariance matrix of the fit parameters. That inversion is failing, indicating that the curvature matrix is not positive definite, even though it is supposed to be so at a minimum point. I don't yet know why that is happening (a bug in the inversion code? a bug in the minimum finder? a bug in the curvature determination?) I do know the bug is new since that documentaiton was written, because it did work for us at the time.

In any case, thanks for the report. I will reply again when I know more.

Aug 4, 2012 at 10:11 AM

Okay, I figured out what's going on here. There is a second problem with the tutorial documentation, and that problem has "infected" your example code. This problem makes the whole fit procedure ill-defined. When this problem is corrected, the fit works fine.

Let's just start again with revised example code:

  // Create a new sample and add some vaguely oscilatory data
  UncertainMeasurementSample data = new UncertainMeasurementSample();
  data.Add(2.0,2.5,0.5);
  data.Add(3.0,2.0,1.0);
  data.Add(5.0,-1.0,0.5);
  data.Add(8.0,-4.0,0.5);
  data.Add(10.0,-1.5,0.5);
  data.Add(12.0, 0.0, 1.0);

  // Define an oscilatory model
  // a[0] ~ amplitude, a[1] ~ period, a[2] ~ phase
  Func<double[], double, double> model = (double[] a, double x) => a[0] * Math.Sin(2.0 * Math.PI * x / a[1] + a[2]);

  // Fit the model to the data, providing an initial guess for the parameters
  FitResult oscilation = data.FitToFunction(model, new double[] { 2.0, 10.0, 0.0 });

  // Print the results
  for (int i = 0; i < oscilation.Dimension; i++) {
    Console.WriteLine(oscilation.Parameter(i));
  }

This code example is better in lots of ways: it is better commented, fixes the Func issue, has more data points than parameters, and uses modern notation for the model delagate. But the key change is this: the model function actually depends on x. Note that, in the original example, x never appeared in the delegate definition. (I suspect this was simply a transcription typo.) This meant that the function was, from the point of view of the model, just a complicated constant computed from the parameters. The best the fit function could do is to adjust the parameters so that the constant equaled the average value of the data. But since there are multiple ways to do this, there are zero directions at the "minimum" and the curvature matrix is indeed not positive definite and therefore indeed not invertible. In any real application, your model function will actually depend on your independent variable, so in practice you should not encounter this problem.

Of course, our code should still behave better when given a non-sensicle model, but that is a considerably less worrying problem than I feared we had when I last wrote. Thanks for pointing out our doc bug. Please let me know if you have any more problems after studying the revised example code.

Aug 6, 2012 at 8:32 AM

I copy past your exemple, but it dosn't work. It seems no convergence : "Meta.Numerics.NonconvergenceException' arrived to Meta.Numerics.Functions.FunctionMath.FindMinimum(Func`2 f, Double[] x) for Meta.Numerics.Statistics.UncertainMeasurementSample`1.FitToFunction(Func`3 function, Double[] start)".

So i try an easier exemple for fitting " y=2x^2+3x+7 " :

var data3 = new Meta.Numerics.Statistics.UncertainMeasurementSample();

            data3.Add(1.0, 12.0, 0.01);
            data3.Add(2.0, 21.0, 0.01);
            data3.Add(3.0, 34.0, 0.01);
            data3.Add(4.0, 51.0, 0.01);
            data3.Add(5.0, 72.0, 0.01);
            data3.Add(6.0, 97.0, 0.01);
            data3.Add(7.0, 126.0, 0.01);
            data3.Add(8.0, 159.0, 0.01);
            data3.Add(9.0, 196.0, 0.01);
            data3.Add(10.0, 237.0, 0.01);

Func<double[], double, double> model = (double[] a, double x) => a[0]*x*x + a[1] * x + a[2]; //V1

Meta.Numerics.Statistics.FitResult polyX = data.FitToFunction(model, new double[] { 2.0, 3.0, 7.0 }); //V1

 

With this exemple, no error but the result is not good :
 -Param n°0:  -3,9657908577218E-18 ± 0,00435194139889244
 -Param n°1:  3 ± 0,0491210625687066
 -Param n°2:  1,27879034232044E-16 ± 0,117615191762516

Just to confirm, i try with FitToPolynomial :

Meta.Numerics.Statistics.FitResult polyX = data3.FitToPolynomial(2);

I got good result :

 ---X^0: 6,99999999999902 ± 0,0194079021706795
 ---X^1: 3,00000000000081 ± 0,0145471440193393
 ---X^2: 1,99999999999983 ± 0,00300058269399403

Aug 7, 2012 at 12:48 AM

Hi Collins,

I think you might have made an error in your comparison of FitToFunction and FitToPolynomial. I notice that you call FitToPolynomial on the data3 object but FitToFunction on the data object (no 3), which presumably contains different data points. When I call FitToFunction on the data3 object with the content as defined in your last post, I get

a[0] = 2 ± 0.000435194139889244
a[1] = 3 ± 0.00491210625687065
a[2] = 7 ± 0.0117615191762515

precisely the same coefficents as for FitToPolynomial. (By the way, one of our unit tests, DataSetTest.FitToFunctionPolynomialCompatibilityTest, which you can find in our source, tests precisely this relationship.) 

With regard to your other problem, the FitToFunction call that gave you a NonconvergenceException, I'm afraid I can't reproduce this. Are you sure you are running precisely the code above, with the same data, model, and starting point? I wonder if perhaps there was a similiar error there, in which a different or changed object that happened to be lying around the same scope was used instead. I am certainly interested in chasing down any nonconvergence behavior you might encounter, but I would need to be able to reproduce it.

Cheers, David

Aug 7, 2012 at 10:52 AM

So sorry !

Effectively, not a full copy past : data instead data3 ! Now Everything is working fine for polynominal fitting.

 

For your exemple with sin() i got :

 A[0]: 3,43481931601377 ± 0,36174906992542
 A[1]: 14,5075811130185 ± 0,989537309537895
 A[2]: 1,33983540139721 ± 0,222553876989726

After i choosed more closer data :

             data3.Add(2.0, 2, 0.5);
             data3.Add(3.0, 1.9, 0.5);
             data3.Add(5.0, 0, 0.5);
             data3.Add(8.0, -2.0, 0.5);
             data3.Add(10.0, 0, 0.5);
             data3.Add(12.0, 2.0, 0.5);

And i got :

 A[0]: 2,07832127074118 ± 0,265143580229382
 A[1]: 10,0458902158885 ± 0,893527889071437
 A[2]: 0,0274078446223988 ± 0,439562043057793

So as a conclusion, everything is working fine :)

Now i can start to play with my own experiments results (i want to fit something like a*sin²(b/x) )...

Thanks,