1
$\begingroup$

Goodmorning everyone, I'm trying to fit the following data

data = {{5000, -1.4639821765728436`}, {6000, -1.4325662500369458`}, 
{7000, -1.4162299682382788`}, {8000, -1.3772742193337653`}, {9000, 
-1.3402034260214057`}, {10000, -1.2943361732789949`}, {11000, 
-1.230247683145763`}, {12000, -1.138513177660941`}, {13000, 
-1.0291857533160162`}, {14000, -0.8532565647149878`}, {15000, 
-0.5956459671206248`}, {16000, -0.24730617369058852`}, {17000, 
    0.14740352730643308`}, {18000, 0.4840565960651153`}, {19000, 
    0.7306087875188423`}, {20000, 0.8947255877423731`}, {21000, 
    1.0133521263419236`}, {22000, 1.1003114109932892`}, {23000, 
    1.1445450355558333`}, {24000, 1.2003397210835882`}, {25000, 
    1.2440706908215582`}, {26000, 1.2546264421376196`}, {27000, 
    1.289309625033251`}, {28000, 1.3089131631916515`}, {29000, 
    1.3265060820517545`}, {30000, 1.3420883816135596`}, {31000, 
    1.3634512116579702`}, {32000, 1.367221122842278`}, {33000, 
    1.376771564509191`}, {34000, 1.3843113868778067`}, {35000, 
    1.4074335088082275`}};

With the following errors:

errors = {0.0027834661607562166`, 0.0033401593929074597`, 
   0.0038968526250587033`, 0.004453545857209947`, 
   0.0025051195446805946`, 0.0027834661607562166`, 
   0.0030618127768318377`, 0.0033401593929074597`, 
   0.0036185060089830813`, 0.0038968526250587033`, 
   0.0016700796964537299`, 0.0008907091714419893`, 
   0.0004731892473285568`, 0.0020040956357444757`, 
   0.0021154342821747244`, 0.0022267729286049734`, 
   0.0023381115750352216`, 0.0024494502214654702`, 
   0.002560788867895719`, 0.002672127514325968`, 
   0.0027834661607562166`, 0.0028948048071864648`, 
   0.003006143453616714`, 0.003117482100046963`, 
   0.0032288207464772115`, 0.0033401593929074597`, 
   0.0034514980393377083`, 0.0035628366857679574`, 
   0.003674175332198206`, 0.003785513978628454`, 
   0.0038968526250587033`};

Using NonlinearModelFit (thanks to @MarcoB for the numerical precision fix)

model[r_, l_, c_, f_] := -ArcTan[(1/(2 \[Pi] f c)  - (2 \[Pi] f l))/r] ;

With[{prec = $MachinePrecision}, 
 nlmf = NonlinearModelFit[SetPrecision[data, prec], 
   model[r, l, c, f], {{r, 28}, {l, 17/20000}, {c, 99 10^-9}}, f, 
   Weights -> 
    1/SetPrecision[errors, prec]^2 VarianceEstimatorFunction -> (1 &),
    WorkingPrecision -> prec]]
nlmf["ParameterTable"]
Show[ListPlot[data], Plot[nlmf["BestFit"], {f, 5000, 35000}]]

The result clearly shows an underestimated error on the parameter r. The relative error is of the magnitude of 0.00000000001%. enter image description here

I have tryed to fit the same function with the same data and the same starting parameters both with Root and Python and they provide a reasonable error. For instance, this is the result I get using Root:

root fit results

I have tried using different Methods in Mathematica, such as Newton, LevenbergMarquardt, Gradient and ConjugateGradient with no improvement. However, I'm not that confortable with Mathematica so my question is if mathematica simply can't give a more reliable estimate of the parameters error for this particular fit or I'm doing something wrong.

Hopefully I've provided enough information for someone to enlight me on how to get a decent fit result from mathematica.

EDIT: As suggested by @JimB, the correlation matrix shows that there's really only 2 parameters that can be estimated. The new fitting function should be something like this:

model[a_, b_, f_] := -ArcTan[(1 - 4 \[Pi]^2 a f^2)/(2 \[Pi] b f)];
nlmf = NonlinearModelFit[data, 
   model[a, b, f], {{a, 8.3 10^-11}, {b, 2.772 10^-6}}, f, 
   Weights -> 1/errors^2, VarianceEstimatorFunction -> (1 &)];
nlmf["ParameterTable"]

where $a = c \times l$ and $b = c \times r$.

The new results are not perfect but definitely an improvement new fit results

However, this defeats the purpose of the fit, which was to get simultaneously an estimate of the 3 elements of an RLC circuit. The approach was flawed since the beginning I guess.

Any tip on how to get better errors on the model with 3 parameters is still welcome, since I still find it weird that both Root and Python manage to do that.

$\endgroup$
  • 4
    $\begingroup$ First off, you want to heed the warnings about numerical precision rather than ignore them. As you show, your result in MMA does not match the other one in Root. Try using With[{prec = $MachinePrecision}, nlmf = NonlinearModelFit[ SetPrecision[data, prec], model[r, l, c, f], {{r, 28}, {l, 17/20000}, {c, 99 10^-9}}, f, Weights -> 1/SetPrecision[errors, prec]^2, (*VarianceEstimatorFunction\[Rule](1&),*) WorkingPrecision -> prec ] ] instead. The values of the parameters are closer. The errors are not yet the same though. $\endgroup$ – MarcoB May 21 at 17:50
  • $\begingroup$ Thanks! I have updated the question. $\endgroup$ – Simone May 22 at 18:33
  • $\begingroup$ Thanks for the edit. But there are just 2 parameters that can be estimated. You just can't get 3 parameters out of 2. Just because you can get 3 estimates out of Root and Python doesn't mean one obtains a correct analysis. It seems that you still think that there ought to be a way to obtain (appropriate) estimates of all 3 parameters. $\endgroup$ – JimB May 22 at 18:51
  • $\begingroup$ Have you obtained similar results in Root and Python for the more parsimonious model? (The model with just 2 parameters.) $\endgroup$ – JimB May 22 at 19:01
  • $\begingroup$ I don't know why but I forgot to try that before the edit. Now (with just 2 parameters) the results obtained with Root and Mathematica are identical. Once again, thanks. $\endgroup$ – Simone May 23 at 12:28
5
$\begingroup$

When I have very little experience with a software function, I find it incumbent on me to figure out what I'm probably doing wrong rather than blaming the software function. (I've learned that the hard way.)

The model presented with 3 parameters ($c$, $l$, and $r$) to fit is

$$-\tan ^{-1}\left(\frac{\frac{1}{2 \pi c f}-2 \pi f l}{r}\right)=-\tan ^{-1}\left(\frac{1-4 \pi ^2 c f^2 l}{2 \pi c f r}\right)$$

But the re-arrangement of the terms on the right show that there's really only 2 parameters that can be estimated. Those two parameters can be constructed in a variety of ways but a straightforward way is to say we can estimate $c\times l$ and $c\times r$.

Had the OP looked at nlmf["CorrelationMatrix"]//MatrixForm, this would have been obvious.

While there might still be issues with obtaining unexpected (as opposed to "wrong") variances, the initial model needs to be fixed and the question edited appropriately. All of the "least chisquared mininum" items seem to be maximum likelihood in disguise and probably should be removed from the question.

| improve this answer | |
$\endgroup$
  • $\begingroup$ Thank you. I have edited the question as you suggested. However, and I hope I'm not going too off topic, isn't the maximum likelihood and the least chi square the same thing for gaussian distributions? Anyway I have removed them since they weren't really usefull. $\endgroup$ – Simone May 22 at 18:41
  • $\begingroup$ Generally that is true. But many folks seem to mix up minimizing a chisquare statistic built for count data (i.e., (o-e)^2/e) and minimizing a chisquare for continuous data. And it just seemed to not be necessary in your question. $\endgroup$ – JimB May 22 at 18:47
1
$\begingroup$

The documentation about variance structures in NonlinearModelFit is not in standard terms (at least in terms familiar to me, a statistician).

When one uses the options Weights -> 1/errors^2, VarianceEstimatorFunction -> (1 &), the following model is assumed:

$$y=f(x) + e \epsilon$$

where $y$ is the response, $f(x)$ is the nonlinear function of $x$, $e$ is the "known" measurement standard deviation, and $\epsilon\sim N(0,1)$. Note that there are no variance parameters to estimate, only the parameters in the nonlinear function of $x$. (The documentation does state that Weights -> 1/errors^2, VarianceEstimatorFunction -> (1 &) means that errors are "treated as the known uncertainty of measurement...and parameter standard errors are effectively computed only from the weights." But that just is not clear enough for me. I've only figured this out by specifying known models with lots of observations with the above error structure. Yes, I'm slow.

Such a variance model assumes there is only measurement error and absolutely no error associated with lack-of-fit to the nonlinear function. While such a situation could occur, I deal in biological systems and have never encountered a time where that would be appropriate. Maybe that is why you felt that the resulting standard errors were too small.

What might be better? Keeping the weights and dropping the VarianceEstimatorFunction one is fitting the following model:

$$y=f(x) + e \sigma \epsilon$$

which looks similar as before except now we have another variance parameter to estimate: $\sigma$. That gives larger standard errors because a specific kind of lack-of-fit can be accommodated. Only you know if that model is appropriate. (And, of course, it shouldn't be chosen if the actual standard errors of the parameters just "feel" better.)

Another model to incorporate some semblance of lack-of-fit is the following:

$$y=f(x) + e \epsilon + \delta$$

where $y$, $f(x)$, $e$, $\epsilon$, are as before but now we add in an independent random deviation from the model ($\delta$) that has a normal distribution with mean 0 and variance $\sigma^2$. So there is a measurement error and an independent lack-of-fit error. I don't know of a direct way that NonlinearModelFit can estimate $\sigma$ but the calculations aren't too difficult using the LogLikelihood function.

Note that to keep the calculations more numerically stable, I've scaled the three parameters ($a$, $b$, and $\sigma$) so that they are on approximately the same scale. That's generally good practice for fitting all models.

(* Log of the likelihood *)
logL = Total[Table[LogLikelihood[NormalDistribution[0, Sqrt[σ^2 + errors[[i]]^2]],
     {data[[i, 2]] - model[a/10^10, b/10^6, data[[i, 1]]]}],
    {i, Length[data]}]];

(* Maximum likelihood estimates *)
mle = FindMaximum[{logL, a > 0 && b > 0 && σ > 0}, {{a, 1}, {b, 3}, {σ, 8}}]
(* {107.561, {a -> 0.917771, b -> 2.85538, σ -> 6.85457}} *)

(* Approximate standard errors *)
cov = -Inverse[D[logL, {{a, b}, 2}] /. mle[[2]]];
n = Length[data];
{aSE, bSE} = Diagonal[cov]^0.5*Sqrt[n/(n - 2)]
(* {0.00101829, 0.0128674} *)

So reasonable standard errors for the parameters depends on the variance structure of the model. Only you might know what that is for your particular application. You could try various variance structures and choose the one with the smallest AIC value.

I'd certainly like to see Mathematica use explicit variance/covariance structures rather than "weights".

| improve this answer | |
$\endgroup$
  • $\begingroup$ Thank you for your detailed and enlightening answer! However I think your other answer might be more useful for anyone with the same problem as mine so I'm choosing that as "most helpful". $\endgroup$ – Simone May 23 at 12:37

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