6
$\begingroup$

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}$).

$\endgroup$
2
  • $\begingroup$ In this case minimizing $\chi^2$ is the same as maximizing the likelihood (assuming multivariate Gaussian errors). Is this minimization of $\chi^2$ statistics for regression models something that's taught in physics classes as opposed to explicitly approaching this as a maximum likelihood problem? (I probably haven't asked the question appropriately. It is not meant as a judgement. I'm just curious as where that approach comes from and if it's commonly taught that way.) $\endgroup$
    – JimB
    Mar 26, 2020 at 3:20
  • $\begingroup$ Both pictures, $\chi^2$ and maximum likelihood, are used pretty widely in physics, depending on the situation. I can't say that one is taught in preference to the other. $\endgroup$ Apr 9, 2020 at 14:08

2 Answers 2

5
$\begingroup$

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

Example code result

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

Plot of data, original function and fit function

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

$\endgroup$
4
  • 1
    $\begingroup$ Would you add a specific example of the use of covariantFit1 ? That would be great. $\endgroup$
    – JimB
    Mar 23, 2020 at 23:17
  • $\begingroup$ What version of Mathematica are you using? I get multiple errors with version 12.0. For example, the statement RandomVariate[MultinormalDistribution[cov]] is missing the sample size parameter and the dimensions of cov is {7,1,7,1} rather than {7,7}. It appears that how you pasted in the code didn't always work. For example when defining cov you have basis 1) which probably should be basis, 1). $\endgroup$
    – JimB
    Mar 24, 2020 at 16:53
  • $\begingroup$ Oops, that was a typo. Fixed. Does it work for you now? $\endgroup$ Mar 24, 2020 at 17:03
  • $\begingroup$ Yes, thanks. That fixed it. $\endgroup$
    – JimB
    Mar 24, 2020 at 17:11
1
$\begingroup$

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

Covariance matrix

(* 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.)

$\endgroup$

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.