Computation of Fitted Parameter Uncertainties

[amyszka]

Hello. I am having trouble understanding how uncertainties of the fitted parameters from a Gaussian model fit are determined. I am getting unrealistic values for parameters (such as Δheight, Δsigma, etc.) when using lmfit.models.GaussianModel() on a dataset with y errorbars.

As an example, I've generated a pure Gaussian with a range of x-values [0, 100] and y-values as the Gaussian of height 5, center 50 and sigma 1.5. I've also generated some normally-distributed error values for each y-value simply as dy = np.abs(np.random.normal(0, 1, len(y)). I've brought this through LMFIT's process using:

gmodel = lmfit.models.GaussianModel()
params = gmodel.make_params(height=5, center=50, sigma=1.5)
result = gmodel.fit(y, params, x=x, weights=1.0/dy)

which then reports parameter errors of basically zero, which should not be expected considering the inclusion of the Δy array as weights.

center: 50.0000000 +/- 1.7798e-14 (0.00%)
sigma: 1.50000000 +/- 1.4858e-14 (0.00%)
height: 5.00000025 +/- 3.7783e-14 (0.00%)

This also occurs when I increase the amount of uncertainty in y (such as using dy = np.abs(np.random.normal(0, 1000, len(y))). Is this the correct way to incorporate the errors of y-values into LMFIT's fitting process?

Here's a plot of the modeled Gaussian.

Thank you.

[newville]

@amyszka Please give a complete example, and always print out and read the fit report. If you are using the built-in lmfit.models.GaussianModel(), then it looks like you have not given an initial value for amplitude. That might explain why the values are stuck at the initial values.

[amyszka]

Hi @newville. I've included an initial value for amplitude and I am still recieving a fit with ~0 uncertainty in the fitted parameters, even though the error on each y-value is very large. Here's my example and fit report:

Example: https://pastebin.com/hT9wRv4d
Fit report: https://pastebin.com/6eGTWUG6

Thank you.

[newville]

@amyszka well, you are saying that the uncertainties are pretty high, by setting the weight to your fake noise weight=1/dy_values. But the array that you are fitting is y_values, which is computed by your own gauss function, and without any added noise. So the fitting method finds a very good fit with the right answer parameters (your gauss function is not normalized to 1, but that's a simple scale factor). And because you explicitly told it that the uncertainties were high when they were, in fact, very low, the uncertainties in the parameters are estimated to be very, very low, like near machine precision.

So, that looks reasonable to me. Perhaps you intended to add some noise to the data to be fitted? It is really hard to over-emphasize the importance of data visualization and exploratory data analysis.

[amyszka]

@newville Yes I see my confusion – there needs to be noise directly present in the y_values for there to be resulting uncertainties on the fitted parameters. I had assumed that having a clean and noise-less y_values but with a noise array incorporated into the weights would return suitable parameter uncertainties; adding the line y_values = y_values + dy_values does indeed result in appropriate parameter uncertainties.

In the case where there is an given array of y_values and a separately-given array of dy_values not related together by that simple addition (for example when analysing real-world data, not generated), how would you suggest making sure the output parameter uncertainties are reliable given the differently-defined inputs? Is just using weight=1/dy_values enough? (Example)

Thank you again for your time.

[newville]

@amyszka

Yes I see my confusion – there needs to be noise directly present in the y_values for there to be resulting uncertainties on the fitted parameters. I had assumed that having a clean and noise-less y_values but with a noise array incorporated into the weights would return suitable parameter uncertainties; adding the line y_values = y_values + dy_values does indeed result in appropriate parameter uncertainties.

I continue to be alarmed (and somewhat suspicious), about the confusion about the use weight and uncertainties. For the very clean data representing a Gaussian, the fit to a Gaussian was very, very good, and the uncertainties in the fitted parameters were very low -- that is correct. Those ARE the appropriate uncertainties.

In the case where there is an given array of y_values and a separately-given array of dy_values not related together by that simple addition (for example when analysing real-world data, not generated),

Sorry, I don't understand this. What is "not related together by that simple addition"?

how would you suggest making sure the output parameter uncertainties are reliable given the differently-defined inputs? Is just using weight=1/dy_values enough? (Example)

well, with "real world" data, the data would generally have some noise. And then one way (and perhaps "the normal way") to use weight is to have that be the inverse of the estimated uncertainties in the data.

In your example dy_values are negative. Uncertainties are strictly positive. So, that's kind of weird.

[amyszka]

@newville Sorry if I was unclear. My main focus is fitting a Gaussian to some scientific data (which should approximate a Gaussian but it is noisy to begin with) as y_values, as well as a separate dy_values array given by the instruments used to obtain the data. This is in contrast with the earlier examples above where the y_values array was generated using an initial Gaussian and then adding on the dy_values to alter the shape (the "simple addition"). I also forgot to abs() the dy_values before copying them over to the example, whoops!

Thank you for your time and patience, I'll see how this goes with the fitting!

[amyszka]

Hello,

My supervisors and I have come up with a solution for our particular examples, in case anyone else comes across this issue. We found that applying the 1/dy weighting rather than 1/dy*2 (which is usual in least-squares methods) did indeed return better results, with smaller final uncertainty values when compared, and are very similar to the values obtained when manually integrating across the error array. We also found that specifying scale_covar=0 in the fitting step (rather than leaving it on by default) solves the issue where near-perfect Gaussian data - but with significant uncertainties on its y-values - returns a fit with parameter uncertainties of practically ±0.

Thank you.

[newville]

@amyszka

My supervisors and I have come up with a solution for our particular examples, in case anyone else comes across this issue. We found that applying the 1/dy weighting rather than 1/dy*2 (which is usual in least-squares methods) did indeed return better results,

Um, that "rather than 1/dy*2 (which is usual in least-squares methods)" is completely wrong. Really, not slightly wrong, not "a simple misunderstanding", not a difference in terminology. It is entirely wrong.

Somehow many people here over the past several months made up some very bad advice that one should weight the residual with (1/data_uncertainty)**2 and then insisted that this point is somehow open to interpretation, confusing, or potentially confusing.
It is all complete, utter nonsense. The more times I see it, the more suspicious I am that someone is trying to be confusing.

The Lmfit authors never say, endorse, or give any credence to the utterly incorrect and frankly stupid and dangerous idea that one should, under any circumstances weight the residual by "(1/data_uncertainty)**2". Really, I do not know where this idea comes from. It is profoundly wrong and does not hold up to dimensional analysis. It's so wrong that it is a useful test for competence. If anyone tells that you should weight by "(1/data_uncertainty)**2" you can be certain that they are completely unsuited for doing numerical analysis and you should never trust another word they say.

:

[karlglazebrook]

This is a rather harsh reply considering. Perhaps I can help clarify what I think is just mismatched terminology.

Normally in model fitting one minimises a $\chi^2$ that looks something like:

$\chi^2 = \sum_i \left(\frac{y_i - y_{model}(x_i)}{\sigma_i}\right)^2 $

where $\sigma_i$ is the noise on the measurements. This is indeed dimensionless. It is very common to refer to this as saying the 'weights are $1/\sigma_i^2$'. It just depends on whether the weights are inside the square or not. Conventions vary. I am sure you have heard of the term 'inverse variance weighting' https://en.wikipedia.org/wiki/Weighted_least_squares which literally means 'square the $1/\sigma$'s' - it is not unreasonable to think that lmfit would need to be fed squared weight values to achieve this.

Regarding lmfit it does say "that weights*(data - fit) is minimized in the least-squares sense" (3/4 down https://lmfit.github.io/lmfit-py/model.html) which does indeed imply it does its own squaring, but takes some digging to find. It's a complex package with complex documentation.

I expect this is tripping up others, and perhaps including an explicit $\chi^2$ example higher up in the documentation would help as this is very common. (And also highlight one wants scale_covar off for this to run correctly if the weights are 1/proper errors). This might help why this keeps coming up.