Memory problems when fitting large (correlated) datasets

Question

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?

Answer 1

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

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 correlation matrix *)
(cor = Table[cov[[i, j]]/Sqrt[cov[[i, i]] cov[[j, j]]], {i, 3}, {j, 3}]) // MatrixForm

We see that the parameter estimators are highly correlated.