I am running NonlinearModelFit based off of some simulated data and trying to fit to a function with more than one parameter. Eventually, I would like to fit to 5 parameters (Right now I'm just trying to get it working with 2 parameters), but the output covariance matrix is almost all 0's. I believe this stems from the fact that the estimated error for some of my parameters is 0, which I know cannot be the case.

I will outline as best as possible my methodology below:

  1. Data creation. (This data takes the model and adds a Gaussian error for each data point.) Note that I had to set the precision higher because of some numerical accuracy issues.

    AU = 149597871000;
    G = 6.67428`20*^-11;
    GMsun = 1.32712442099`20*^20;
    GMjup = GMsun/1047.348644`20;
    rJup = 5.2`20 AU;
    lambda = 20 AU;
    dela = 10^-11;
    Data[dist_] := {
       GMsun/dist^2 + (GMjup dist)/(dist^2 + rJup^2)^(3/2) + 
        RandomVariate[NormalDistribution[], WorkingPrecision -> 20] dela
    Model[dist_, alpha_, jupiter_] := 
      GMsun/dist^2 (1 + alpha Exp[-dist/lambda]) +
       (G jupiter dist)/(dist^2 + rJup^2)^(3/2);
    Dat = Table[Data[x], {x, AU, 100 AU, AU}];
  2. Then for a simple case I tried to fit to my model:

    NLM = NonlinearModelFit[Dat, Model[dist, alpha, jupiter], {{alpha, 10^-4}, {jupiter, 10^27}}, {dist}, WorkingPrecision -> 20];

But what this does is output the following:

Table of parameter values

Similarly, when I output the covariance matrix I get:

(* -> {{2.762131253783*10^-18, 0},
       {0,                     0}} *)

My question is this: Why do I get 0's in my covariance matrix? Clearly these 0's stem from the 0 standard error for the jupiter parameter, but shouldn't there be some associated standard error with jupiter just as there is with alpha?

Thank you so much

  • $\begingroup$ I edited your code to deal with the numerical issues in a more robust way so that this doesn't distract people's attention from the substance of your question. Hope this is okay. $\endgroup$ Commented Jul 17, 2013 at 2:05
  • $\begingroup$ I noticed, thank you! $\endgroup$
    – diracula
    Commented Jul 17, 2013 at 2:37

1 Answer 1


I suspect a singular value decomposition being computed under the hood is returning zeros where it perhaps shouldn't. Note that NonlinearModelFit has an undocumented Tolerance option.

FilterRules[Options[NonlinearModelFit], Tolerance]

(*{Tolerance -> Automatic}*)

For some reason it doesn't allow one to set it to exactly zero but we can choose something really small, here 10^-50.

NLM = NonlinearModelFit[Dat, Model[dist, alpha, jupiter], 
        {{alpha, 10^-4}, {jupiter, 10^27}}, {dist}, 
           Tolerance -> 10^-50, WorkingPrecision -> 20];

This seems to alleviate the problem.


enter image description here

  • $\begingroup$ Oh, interesting, this is very helpful. And if you compute the covariance matrix, you now get the correct values along the diagonal. However, there are still zeroes everywhere else in the matrix, which shouldn't be the case. Do you have any idea why this would happen? $\endgroup$
    – diracula
    Commented Jul 18, 2013 at 17:55
  • 1
    $\begingroup$ I suspect it has something to do with the scale of the values being so radically different that 20 digits precision isn't enough to capture it. Try working at much higher precision or rescaling your data and model somehow so that things are measured at the same order of magnitude. $\endgroup$
    – Andy Ross
    Commented Jul 18, 2013 at 18:25
  • 1
    $\begingroup$ Oh, so working with approximately 40 digits makes it work fine, that's very interesting, thanks a bunch! $\endgroup$
    – diracula
    Commented Jul 19, 2013 at 7:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.