0
$\begingroup$

I am trying to fit some complicated function that depends on 18 parameters. The problematic is that I have a large correlated dataset (336 points).

When fitting, for simplicity, taking only the diagonal errors with NonlinearModelFit, the fit is extremely quick (0.15s). If, instead, I construct the chi-square function myself and use FindMinimum, it takes 8s and get the same result. If using NMinimize with good boundaries for the parameters, I get the same result in 28s.

Now, the problem is when I try to fit with the full covariance matrix. Now, I construct the chi squared function myself and try to use FindMinimum or NMinimize. The problem is that I get the following message:

No more memory available. Mathematica kernel has shut down. Try quitting other applications and then retry

and the fit of course stops and the kernel is empty.

Actually, if I try just to create the chi-square function itself, the same problem happens. I guess the problem is a sum of 336^2 terms with a complicated function resulting in a too long expression.

My naive workaround was to modify the function to be fitted f[x_,a_,b_,...] as f[x_?NumericQ,a_?NumericQ,b_?NumericQ,...] to prevent symbolic evaluation and to get ran out of memory. This wasn't actually enough, see here. Further, even for simple fits, setting the inputs as ?NumericQ will slow down the minimization considerably. I used then in FindMinimum the Hold[] attribute as in the previous thread, but then the fit is far from ideal (even if considering diagonal errors only) and gives the message

Encountered a gradient that is effectively zero. The result returned may not be a minimum; it may be a maximum or a saddle point.

Any ideas on how to handle this?

$\endgroup$
6
  • 1
    $\begingroup$ Your question is hard to answer in general. Perhaps if you showed your code so far and a sample of your dataset then some specific suggestions could be made. $\endgroup$
    – MarcoB
    Apr 25, 2022 at 12:43
  • $\begingroup$ Thanks for the comment. I will try to provide a MWE (the function is just too lengthy and with a simpler one no problem arises). By the way, I could try on a (better) different computer and it works without problem. However, I would prefer to use my own computer. $\endgroup$
    – pablo
    Apr 25, 2022 at 14:04
  • $\begingroup$ Giving the fixed and random effects (covariance structure) in LaTeX form rather than just Mathematica code would also be helpful so that the definition of "correlated data" could be understood. $\endgroup$
    – JimB
    Apr 25, 2022 at 14:17
  • $\begingroup$ Essentially, the experimentalists provide a set of data points (336) together with the covariance matrix (336*336 entries). I suspect the form its not essential, as the problem comes when creating such a long chi square. I obtain the same error if using a covariance matrix with the same size with random numbers as input. I believe the key points are (i) the length of the covariance matrix (ii) the complicated function. The latter is given in analytic form including logs, square roots, powers,... but is rather lengthy. I checked that, for a simple function, the problem would not exist. $\endgroup$
    – pablo
    Apr 26, 2022 at 7:53
  • $\begingroup$ From your description of the covariance matrix it appears that you are not estimating any covariance terms but rather those are given as known. (That seems more than a bit unbelievable that covariance terms would be known but the terms almost always easier to estimate are not. But let's ignore that for now.) Why you think in terms of a minimizing a chisquare is also a bit odd. Why not directly find the maximum likelihood estimates by maximizing the likelihood? Mathematica's LogLikelihood function would seem to be what you want to use. $\endgroup$
    – JimB
    Apr 26, 2022 at 14:34

1 Answer 1

2
$\begingroup$

I don't know what particular chisquare statistic you're trying to minimize (or even why you'd want a chisquare statistic to minimize - certainly not the Sum(O-E)^2/E type) but obtaining the maximum likelihood estimates of the parameters should be relatively straightforward.

I'm also assuming that you think you know the true covariance matrix as opposed to estimating it along with the other parameters.

(* Sample size *)
n = 336;

(* Construct covariance matrix and associated functions *)
(* This is an autoregressive covariance matrix *)
σ = 4;
ρ = 0.5;
Σ = Table[σ^2 ρ^Abs[i - j], {i, n}, {j, n}];

(* Generate a simple prediction equation and simulated data *)
(* Normally distributed error from covariance matrix *)
SeedRandom[12345];
ϵ = RandomVariate[MultinormalDistribution[Σ], 1][[1]];
(* Function to fit *)
f[x_, a_, b_, c_] := a + b x + c Log[x]
x = Range[n];
y = f[x, 1, 1/2, 2] + ϵ;

(* Plot data *)
ListPlot[Transpose[{x, y}]]

Scatter plot of data

Now define the log of the likelihood. I've removed the parts that remain constant but have also included the full likelihood in comments in case one needs to compare results with other software.

(* logAbsDet=Log[Abs[Det[Σ]]] *)
ΣInverse = Inverse[Σ];
(* logL[y_,x_,a_,b_,c_]:=-(n Log[2π]+logAbsDet+(y-f[x,a,b,c]).ΣInverse.(y-f[x,a,b,c]))/2 *)
logL[y_, x_, a_, b_, c_] := -(y - f[x, a, b, c]) . ΣInverse . (y - f[x, a, b, c])/2 // N

(* Get starting values for parameters *)
nlm = NonlinearModelFit[Transpose[{x, y}], f[z, a0, b0, c0], {a0, b0, c0}, z]["BestFitParameters"]
(* {a0 -> 1.25527, b0 -> 0.507967, c0 -> 1.45971} *)

(* Maximum likelihood estimates *)
mle = FindMaximum[logL[y, x, a, b, c], {{a, a0}, {b, b0}, {c, c0}} /. nlm]
(* {-166.9, {a -> 0.252364, b -> 0.505238, c -> 1.75736}} *)

(* Estimated parameter covariance matrix *)
(cov = -Inverse[(D[logL[y, x, a, b, c], {{a, b, c}, 2}]) /. mle[[2]]]) // MatrixForm 

Estimated parameter covariance matrix

(* Estimated parameter correlation matrix *)
(cor = Table[cov[[i, j]]/Sqrt[cov[[i, i]] cov[[j, j]]], {i, 3}, {j, 3}]) // MatrixForm

Estimated parameter correlation matrix

We see that the parameter estimators are highly correlated.

$\endgroup$
4
  • $\begingroup$ Hi, and sorry for my delay. Indeed, that's the way I proceeded (perhaps I was not quite explicit). Essentially, I created your logL function for my data (appart from the minus sign and the factor of 2, that's what I call the chi square). The problem is that, when I create the logL function, my computer runs out of memory and the kernel is restarted. I guess the problem is the length of my function (containing logs, square roots, etc.. and 18 parameters), as doing an example like yours works easily and better computers (I have 4GB ram on my laptop) proceeds easily. $\endgroup$
    – pablo
    May 4, 2022 at 10:15
  • $\begingroup$ If the objective is to get your function and estimation working on your 4GB laptop, then I would suggest determining where things stop working. Taking the inverse of the covariance matrix? If so, then it seems hopeless. If it's at the NonlinearModelFit stage, then it might also be hopeless. But if it's at the last stage (FindMaximum), then maybe some of your equation is linear in some of the parameters which means you can reduce the number of parameters in the FindMaximum. But without knowing the function you're trying to fit, that's only a guess. $\endgroup$
    – JimB
    May 4, 2022 at 16:01
  • $\begingroup$ Actually, the problem is at computing the logL. The inverse covariance matrix is computed without problems indeed. Apparently, Mathematica precomputes the expression for logL (note we have 336^2 long expressions). That's why I tried to append ?NumericQ to prevent evaluation before the model is given a point and produces a number. However, it turns out that fitting this way is terribly time consuming and couldn't produce results with FindMinimum (I don't have your overall - sign, so it is equivalent) or NMinimize. That was funny... Thanks a lot for your help anyway! $\endgroup$
    – pablo
    May 6, 2022 at 6:36
  • $\begingroup$ With 18 parameters (and based on your recent question), I would bet that your model is linear in the parameters. That means generalized least squares is what you want (as I mentioned in an answer to your most recent question). That is not an iterative approach which means there is an explicit formula for the estimates and the covariance matrix for the estimates. $\endgroup$
    – JimB
    Nov 9 at 16:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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