1
$\begingroup$

I am trying to fit data that is correlated, such that the $\chi^2 = (f(x_i,pars)-y_i)(Cov^{-1})_{ij}(f(x_j,pars)-y_j)$. However, NonlinearModelFit only accepts diagonal weights. Of course, the obvious approach would be to use FindMinimum or NMinimize.

However, I find that either method has an inferior performance as compared to NonlinearModelFit. In particular, and taking a diagonal covariance matrix to compare NonlinearModelFit vs. FindMinimum or NMinimize, I find that the minima found by the latter two methods are not as small as those found by NonlinearModelFit. Also, the time it takes is, at least, 10 times bigger.

The time performance of course gets worse when building the $\chi^2$ for the full covariance matrix, that I can handle so far only via FindMinimum or NMinimize.

In principle, this might not be a problem if I only have to wait for a couple of minutes. However, since I have to repeat this fit around 8000 times for a MonteCarlo procedure, this becomes extemely important.

Any ideas how to use NonlinearModelFit to minimize a given $\chi^2$?

(Edited below)

As a MWE, the covariance matrix MyCovSqMat

MyCorr = {{1, 0.7, 0.4, 0.1}, {0.7, 1, 0.7, 0.4}, {0.4, 0.7, 1, 
    0.7}, {0.1, 0.4, 0.7, 1}};
MyErr = {0.8, 0.5, 0.6, 0.9};
MyCovSqMat = 
  Table[MyCorr[[i, j]] MyErr[[i]] MyErr[[j]], {i, 1, 4}, {j, 1, 4}];

Then data with noise (MyData) is generated for a linear function 2+3x

MyNoise = 
  RandomReal[MultinormalDistribution[{0, 0, 0, 0}, MyCovMat], 1];
MyData = Table[{i, 2 + 3 i + MyNoise[[1, i]]}, {i, 1, 4}]

Then we could fit taking only the diagonal elements of the covariance matrix either through NonlinearModelFit or minimization methods:

NonlinearModelFit[MyData, a + b x, {a, b}, x, Weights -> MyErr^-2]
MyChiSq = 
  Sum[((a + b MyData[[i, 1]] - MyData[[i, 2]])/MyErr[[i]])^2, {i, 1, 
    4}];
FindMinimum[MyChiSq, {a, b}]

However, the full covariance matrix can only be handled through FindMinimum or NMinimize:

MyInvCovSqMat = Inverse[MyCovSqMat];
MyFullChiSq = 
  Sum[(a + b MyData[[i, 1]] - MyData[[i, 2]]) MyInvCovSqMat[[i, 
      j]] (a + b MyData[[j, 1]] - MyData[[j, 2]]), {i, 1, 4}, {j, 1, 
    4}];
FindMinimum[MyFullChiSq, {a, b}]

However, what I would like is to minimize MyFullChiSq above via NonlinearModelFit due to its apparently superior performance in more complicated fits

Solution

Finally there was an easy fix following @Ulrich Neuman:

NonlinearModelFit[{1, 0}, MyFullChiSq x, {a, b}, x]
$\endgroup$
5
  • 1
    $\begingroup$ Please, present a data for your $\chi^2$. In other case your question is too general to answer it. $\endgroup$
    – user64494
    Nov 8, 2023 at 9:15
  • $\begingroup$ I understand, my particular case is complicated. I will try to find a MWE, that might also be useful for myself anyhow. In any case, essentially I would need to have some function that depends on n variables and that should be minimized bye means of NonlinearModelFit $\endgroup$
    – pablo
    Nov 8, 2023 at 9:33
  • $\begingroup$ I have edited the post and included a MWE hoping my question is clearer $\endgroup$
    – pablo
    Nov 8, 2023 at 10:56
  • $\begingroup$ At least for the example you give, NonlinearModelFit takes twice the amount of time as FindMinimum. Also, you don't get any measures of precision. To assume that the covariance matrix is known but the regression parameters are not, is hard to believe. I think a more realistic example is needed. $\endgroup$
    – JimB
    Nov 9, 2023 at 0:43
  • $\begingroup$ I understand, but that's not necessarily the case in more complicated examples, as I commented. Unfortunately, the inpout required is too long. Regarding the covariance matrix I mean that from the data, not the one for the parameters. $\endgroup$
    – pablo
    Nov 9, 2023 at 21:15

2 Answers 2

3
$\begingroup$

As MWE we could examplary try to fit a circle to random data xy

xy = Table[(1 + RandomReal[{-.1, .1}]) {Cos[phi], Sin[phi ]}, {phi, 
RandomReal[{-Pi, Pi}, 100]}];

Add 0 to the data

xy0 = Map[Join[#, {0}] &, xy];

and fit for unknown radius R

fit = NonlinearModelFit[xy0, (x^2 + y^2 - R^2)^2, R, {x, y}]
fit["BestFitParameters"](*{R -> 1.0008}*)

Hope it helps!

$\endgroup$
4
  • $\begingroup$ That's not exactly what I need (it seems a standard NonlinearModelFit). Thanks anyhow! I will try to provide now a MWE $\endgroup$
    – pablo
    Nov 8, 2023 at 10:48
  • $\begingroup$ @paplo You try to minimize a functional>=0 depending on xi,yi,pars , that's what I show in my answer! $\endgroup$ Nov 8, 2023 at 10:54
  • $\begingroup$ I just edited my entry, hope it makes it clearer. The trick to add 0 and minimize it is a nice one! but, still, you have many data points, while I have a single function (say, MyFullChiSq above) that depends on some parameters and should be minimize (the 0 in this case might help, but I don't have a data set) $\endgroup$
    – pablo
    Nov 8, 2023 at 10:56
  • $\begingroup$ Ok I see how to fix it, that helped a lot in order to figure out the solution! :D $\endgroup$
    – pablo
    Nov 8, 2023 at 11:04
0
$\begingroup$

The standard approach (going back to at least 1935) for your linear model ($a+bx$) with a known covariance matrix is Generalized Least Squares (https://en.wikipedia.org/wiki/Generalized_least_squareshttps://en.wikipedia.org/wiki/Generalized_least_squares).

Directly programming this approach is far faster than using NonlinearModelFit and just as important you can obtain the covariance matrix for the estimates of the parameters.

Using SeedRandom[12345]; and using your data generation approach we have

(x = Transpose[{ConstantArray[1, Length[MyData]], MyData[[All, 1]]}]) // MatrixForm

Design matrix

y = MyData[[All, 2]]
(* {4.209, 8.29689, 11.1309, 14.5286} *)
ΣInv = MyInvCovSqMat;
AbsoluteTiming[estimates = Inverse[(Transpose[x] . ΣInv . x)] . Transpose[x] . ΣInv . y]
(* {0.000058, {1.75184, 3.2576}} *)
(cov = Inverse[Transpose[x] . ΣInv . x]) // MatrixForm

Covariance matrix

se = Sqrt[Diagonal[cov]]
(* {0.949452, 0.356867} *)

This gives the same estimates of the coefficients as in your approach and as in @UlrichNeumann 's answer but with the added bonus of the essential covariance matrix for having some idea of the precision of the estimates.

In real life "knowing" the covariance matrix without error is at best wishful thinking.

$\endgroup$
3
  • $\begingroup$ However, my problem is far from linear, nor it can be linearized, and cannot be solved analytically as in linear model cases, that was just a MWE to illustrate my needs. Somehow, the simpler version (diagonal covariance matrix) is handled much faster by nonlinearmodelfit, that the FindMinimum, that a nonlinear problem requires. Of course, the covariance matrix for the parameters is important, but in mycase will be handled through a MC procedure. $\endgroup$
    – pablo
    Nov 9, 2023 at 21:14
  • 1
    $\begingroup$ Then your example should be nonlinear as it didn't illustrate your needs. You have a single solution in mind. This forum can be great especially when someone isn't as restricted. This is not to be critical but there's more out there than any of us thinks we know. Finally, you currently don't have measures of precision which are essential to convinces others of your results (unless maybe you're doing some sort of bootstrap process). $\endgroup$
    – JimB
    Nov 9, 2023 at 21:35
  • $\begingroup$ Just as another check to make sure: Should I assume that your model is not linear in the coefficients? I ask because $y=a_0+a_2 e^x+a_3 x^3+a_4 \sin(x)$ is a linear model with respect to $a_0$, $a_1$, $a_2$, $a_3$, and $a_4$ and therefore generalized least squares would be appropriate for such a model. $\endgroup$
    – JimB
    Nov 12, 2023 at 1:22

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.