Problem with Covariance Matrix Output in NonlinearModelFit

Question

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_] := {
      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];
   
   NLM["ParameterTable"]
   

But what this does is output the following:

Similarly, when I output the covariance matrix I get:

NLM["CovarianceMatrix"]
(* -> {{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

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.

NLM["ParameterTable"]