What exactly is "ParameterErrors" of NonlinearModelFit?

Question

I am trying to get what I call a "$\chi^2$ fit". This means taking a function $f(x,p)$ of the variable $x$ with a parameter $p$ and finding the value of $p$ that minimizes the "distance" from a data set $\{x_i, d_i\}$, i.e. minimizes the quantity

$$\chi^2=\frac{\left(f(x_i,p) - d_i\right)^2}{d_i}$$

See this Yale course notes for a reference.

You can find the value of the parameter $p$ with NMinimize, of course, or use NonlinearModelFit. I figured out that to get the same result, i.e. perform the same minimization, you need to use

nlm = NonlinearModelFit[data, 1 + p x^2, {p}, x, Weights -> (1/#2 &)]

where I have taken f[x,p] = 1 + p*x^2 just to give an example

This procedure is supposed to give you the best estimate of the parameter $p$ and an error is associated to to this estimate. I thought that nlm["ParameterErrors"] would give this error, but I do not find any documentation for this.

Usually, for a one parameter estimation at 68% confidence level, the error comes from looking at the curve $\chi^2=1+\chi_\min^2$ where $\chi_\min^2$ is the value of the function $\chi^2$ for the parameter $p$ that minimizes $\chi^2$ (in other words the value at the best fit).

By direct computation and confronting with "ParameterErrors" it seems that "ParameterErrors" does not give this type of error estimate.

Now, is anybody aware of what exactly "ParameterErrors" does and how to instruct NonlinearModelFit to do what I need to do?

Answer 1

To fit a model and obtain a measure of precision for the parameter estimates one needs to consider the error structure in addition to matching the data with a function to be minimized or maximized. And finding the value of a parameter that minimizes some function does not necessarily result in a “best fit.” Again, it depends on what is a reasonable error structure for generating similar data.

Usually the $\chi^2$ fit involves count data (as shown in all of the examples of the notes you reference) but you don’t explicitly state that in your question. Using your example of

$$f(x,p)=1+p x^2$$

suppose $d_i \sim Poisson(1+p x_i^2)$. For the set of observed counts $d_i$ you want to calculate

$$\chi^2 = \sum_{i=1}^n {{(f(x_i,p)-d_i)^2}\over{d_i}}$$

But there are problems when $d_i=0$ and this isn't the definition in the referenced notes. The definition in the notes is

$$\chi^2 = \sum_{i=1}^n {{(observed_i-expected_i)^2}\over{expected_i}}$$

which translates to

$$\chi^2 = \sum_{i=1}^n {{(d_i-f(x_i,p))^2}\over{f(x_i,p)}}$$

where $f(x_i,p)>0$.

The most straightforward way to estimate $p$ and obtain an appropriate measure of precision is to use maximum likelihood. That can be done directly by constructing the likelihood function and maximizing it with respect to $p$ or using GeneralizedLinearModelFit. Your use of NonlinearModelFit results in an estimator with nearly the same optimal properties of the maximum likelihood estimator. (A Bayesian approach is also reasonable but I've left that out.)

To show the differences (and similarities) among the different estimation techniques first generate some data:

(* Generate some observed counts where p = 0.5 *)
n = 100;
x = Table[i/10, {i, n}];
e = 1 + p x^2;
SeedRandom[1234];
d = Flatten[Table[RandomVariate[PoissonDistribution[e[[i]] /. p -> 1/2], 1], {i, n}]];
data = Transpose[{x, d}];

Now estimate $p$ in a variety of ways:

(* Minimize sum of (observed-expected)^2/expected *)
pmin = p /. NMinimize[{Total[(d - e)^2/e], p > 0}, p, WorkingPrecision -> 100][[2]];
N[pmin, 20]
(* 0.52824973972704747790 *)

(* NonlinearModelFit with weights *)
pnlm = p /. 
   NonlinearModelFit[data, 1 + p xx^2, {p}, xx, Weights -> (1/#1 &),
     WorkingPrecision -> 100]["BestFitParameters"];
N[pnlm, 20]
(* 0.50953435937653171258 *)

(* Maximum likelihood estimate *)
(* First get log of likelihood *)
logL = -Sum[1 + p x[[i]]^2, {i, n}] + 
   Sum[o[[i]] Log[1 + p x[[i]]^2], {i, n}];
pmle = p /. NMaximize[{logL, p > 0}, p, WorkingPrecision -> 100][[2]];
N[pmle, 20]
(* 0.51267504624247847170 *)

(* GerneralizedLinearModelFit - logically the same as maximum likelihood estimate *)
glm = GeneralizedLinearModelFit[data, xx^2, xx,
   ExponentialFamily -> "Poisson", LinearOffsetFunction -> 1,
   IncludeConstantBasis -> False, LinkFunction -> "IdentityLink",
   WorkingPrecision -> 100];
N[glm["BestFitParameters"][[1]], 20]
(* 0.51267504624247847170 *)

An estimate of the standard error of the maximum likelihood estimator of $p$ is found with

N[glm["ParameterErrors"][[1]], 20]
(* 0.012605361614497436926 *)