Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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:

Table of parameter values

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

share|improve this question
    
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
add comment

1 Answer

up vote 1 down vote accepted

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

enter image description here

share|improve this answer
    
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
1  
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
1  
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
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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