Chi^2 fitting for correlated data

Question

Suppose you have $N$ correlated data points $\vec{y}_\mathrm{data}$ and a model that is a function of $M$ parameters $\vec{x}$. The associated $\chi^2$ statistic is

$\chi^2 = (\vec{y}_\mathrm{data} - \vec{y}_\mathrm{theo}(\vec{x})) \cdot C^{-1} \cdot (\vec{y}_\mathrm{data} - \vec{y}_\mathrm{theo}(\vec{x}))$ ,

where $\vec{y}_\mathrm{data}$ is the vector with the data points, $\vec{y}_\mathrm{theo}(\vec{x})$ is the fit function and $C$ is the covariance matrix associated the data points. To best fit the model, one minimizes $\chi^2$ with respect to $\vec{x}$.

What is the best function on Mathematica to do this? NonlinearModelFit doesn't handle correlated data and FindMinimum doesn't provide useful statistics (like the covariance matrix associated with $\vec{x}$).

Answer 1

Here is a wrapper for NonlinearModelFit that handles data with a known correlation matrix:

covariantFit1::usage = "Wrapper for NonlinearModelFit that \
includes a covariance matrix.  It does not handle the simple \
data format or handle constraints.  It is computationally \
inefficient since it evaluates \"form\" n^2 times for every \
step, where n is the number of data points.  Some methods \
associated with the fit function itself will not work."; 
Options[covariantFit1] = FilterRules[Options[NonlinearModelFit],
            {Except[Weights],Except[VarianceEstimatorFunction]}]; 
covariantFit1[cov_, data_?MatrixQ, Except[_List, form_], pars_, 
   vars_, opts : OptionsPattern[]] := 
 Block[{err2, oo, ii, transform, ff}, 
  {err2, oo} = Eigensystem[cov]; 
  transform = 
   oo.Map[(form/.MapThread[Rule, {vars, Drop[#, -1]}]) &, data];
  ff = NonlinearModelFit[oo.Map[Last, data] ,
    (* Expand the variable transform, 
    but don't take the part unless it is numeric. *)
    If[ii > 0, #[[ii]]] &[transform], pars, ii, 
       VarianceEstimatorFunction->(1&), Weights->1/err2, opts];
  Unprotect[FittedModel];
  ff["chiSquared"] =  Block[{x=ff["FitResiduals"]},x.(x/err2)];
  ff["BestFit"] = form/.ff["BestFitParameters"];
  ff["Function"] = Function[Evaluate[
     form/.Join[ff["BestFitParameters"],
     MapIndexed[Rule[#1, Apply[Slot, {#2[[1]]}]]&, vars]]]];
  Protect[FittedModel];
  ff];

Here is some example code to try it out:

func = Function[x, c0 Exp[-c1*Norm[x] - c2*x.x]]; 
values = {c0 -> 1.0, c1 -> -2.5, c2 -> 0.8}; 
basis = Table[{x}, {x, 1, 4, 0.5}]; 
cov = Outer[(0.2*(1.01*Exp[-Norm[#1 - #2]^2/4.0] - 0.01)) &, 
      basis, basis, 1]; 
data = MapThread[Append, {basis,
       RandomVariate[MultinormalDistribution[cov]] + 
       Map[func, basis] /. values}]; 
ff1 = covariantFit1[cov, data, func[{x}], 
      {{c0, 1.0}, {c1, -2}, {c2, 0.75}}, {x}]; 
Print[Row[{"chi^2:", ff1["chiSquared"], "for", Length[data] -
           Length[ff1["BestFitParameters"]], "d.o.f."}, " "]];
ff1["ParameterTable"]

And plot the result (note that ff1[x] won't work):

Needs["ErrorBarPlots`"];
Show[ErrorListPlot[MapThread[{#1, ErrorBar[Sqrt[#2]]} &, 
           {data, Diagonal[cov]}], PlotRange -> All], 
  Plot[{ff1["Function"][x], func[{x}]/.values},
      {x, 0, 5}]]

Here is a fancier version on GitHub: covariant-fit.m. It also includes a more efficient version covariantFit2 which calls FindMinimum directly.

Answer 2

An alternative approach is to find the maximum likelihood estimators of the parameters given a known covariance structure. (In addition, this approach is amenable to estimating the parameters of the covariance matrix if just the structure of that matrix is known.) One of the advantages of the other answer is an easily obtained nicely formatted table of the estimates and standard error.

(* Generate data in the same way as the previous answer *)
func = Function[x, c0 Exp[-c1*Norm[x] - c2*x.x]];
values = {c0 -> 1.0, c1 -> -2.5, c2 -> 0.8};
basis = Table[{x}, {x, 1, 4, 0.5}];
cov = Outer[(0.2*(1.01*Exp[-Norm[#1 - #2]^2/4.0] - 0.01)) &, basis, basis, 1];
SeedRandom[321];
data = MapThread[Append, {basis, RandomVariate[MultinormalDistribution[cov]] + 
  Map[func, basis] /. values}];

(* Log of the likelihood *)
logL = LogLikelihood[
   MultinormalDistribution[cov], {data[[All, 2]] - (func[{#}] & /@ data[[All, 1]])}];

(* Get initial estimates of parameters by ignoring the covariance structure *)
nlm = NonlinearModelFit[data, func[{x}], {c0, c1, c2}, x];
{c00, c10, c20} = {c0, c1, c2} /. nlm["BestFitParameters"];

(* Maximum likelihood estimates *)
mle = FindMaximum[{logL, c0 > 0}, {{c0, c00}, {c1, c10}, {c2, c20}}];
(* {10.9581, {c0 -> 0.994929, c1 -> -2.51037, c2 -> 0.802505}} *)

(* Covariance matrix for parameter estimators *)
(pcov = -Inverse[(D[logL, {{c0, c1, c2}, 2}]) /. mle[[2]]]) // MatrixForm

(* Standard errors of parameter estimators *)
Diagonal[pcov]^0.5
(* {0.0818945, 0.0434231, 0.0125824} *)

The results match with the other given answer and (with enough data) allows for the estimation of the covariance parameters. Maybe in Physics and Chemistry might there be cases where the parameters for the fixed effects of the model are unknown but the associated variance/covariance terms are known. But that is nearly unheard of in the biological sciences. (In other words, how many times is the mean unknown but the variance is known? The answer is: almost never.)