# Problem with Covariance Matrix Output in NonlinearModelFit

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.6742820*^-11;
GMsun = 1.3271244209920*^20;
GMjup = GMsun/1047.34864420;
rJup = 5.220 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

-
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. – Oleksandr R. Jul 17 '13 at 2:05
I noticed, thank you! – diracula Jul 17 '13 at 2:37

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"]


-
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? – diracula Jul 18 '13 at 17:55
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. – Andy Ross Jul 18 '13 at 18:25
Oh, so working with approximately 40 digits makes it work fine, that's very interesting, thanks a bunch! – diracula Jul 19 '13 at 7:44