7
$\begingroup$

I have a dataset consisting of (t,y) Tuples, where each y has a known measurement error (standard deviation). I would like to fit a trapaze-function into the data and obtain the fitting paramaters and propagate the measurement errors onto those paramters. I was able to do the fitting but I can't managage to do the error propagation and I don't understand the methods used by NonLinearModelFit to calculate such errors, since the documentation is very vague. Since error propagation is a very standard procedure in physical sciences I think that there shoud be some way to accomplish this with mathematica (for fits).

An overview over the data: (the time values on the x-Axis are not displayed properly, the total time interval in which the data was aquired is about 0.5) data with errors

The fitted function and the parameters used: fitted function and the used parameters What I was able to do/find out:

I was able to do the fitting using NonLinearModelFit using the following model:

model=Piecewise[{{F+(FTransit-F)/(t2-t1)*(x-t1),t1<=x<t2}, {FTransit,t2<=x<t3},{FTransit+(F-FTransit)/(t4-t3)*(x-t3),t3<=x<t4}},F];

nlmVC=NonlinearModelFit[lightDataVCGraph,{model,{2.4597044744907`*^6<=t1<=2.4597046155903`*^6},{2.4597044744907`*^6<=t2<=2.4597046155903`*^6},{2.4597044744907`*^6<=t3<=2.4597046155903`*^6},{2.4597044744907`*^6<=t4<=2.4597046155903`*^6}},{{F,3.65},{FTransit,3.59},{t1,2.4597045370702376`*^6},{t2,2.4597045465046*10^6 +0.4},{t3,2.4597045465046*10^6 +0.2},{t4,2.45970457042957`*^6 +0.8}},x ];

In my understanding the "standard error" from the parameter table given by NonLinearModelFit via nlmVC[{"BestFit","ParameterTable"}]is the estimated standard deviation of the particular parameter from the model which I am interested in(correct?). Since I didn't give the fitting algorithm any measurement errors I assume that the "standard errors" are solely based on the goodness of the fit of the model and the sensitivity of the model considering the particular parameter.

Now I want to take into account that the y-part of each data touple has a measurement error/uncertainty, which should be propagated into the errors of the fitted parameters (and probably increase them according to th size of the errors).

What I found about how this is suggested to be done(optional if you know how it's done;D)

I found this guide: https://reference.wolfram.com/language/howto/FitModelsWithMeasurementErrors.html[FitModelsWithMeasurementErrors][3] which fits to my problem relatively well.

They suggest two things to incorporate the errors from the measurement:

  1. They input the measurement errors as weights into NonLinearModelFit by putting Weights -> 1/errors^2 into the options. In my understanding this weights the squared errors in the least squares-error function according to theire standard deviation/variance but does not change the way the "standard errors" of the parameters are computed. It minimally improves the fit(my errors are reletively constant) but does not change the obtained parameter errors from the fitting very much.

  2. They suggest setting VarianceEstimatorFunction->(1&) as an option in NonLinearModelFit.

To treat the weights as being computed from measurement errors, you can use the VarianceEstimatorFunction option in addition to Weights. VarianceEstimatorFunction explicitly defines the variance scale estimator that is used. For measurement errors, you want standard errors to be computed only from the weights and so the variance estimate should be the constant 1:<

According to this second guide: http://www.columbia.edu/~nas2173/tutorial.pdf this option "activates" the error propagation of the measurement errors (page 10 top) but doesnt explain how/why this works.

Now I dont understand what the variance estimator is in theory, why it should be set to constant 1 and how NonLinearModelFit should even get any information about the measurement errors since I only give it the weights which contain the errors only implicitely. I also dont understand why I would want the parameter errors to depend only on the measurement errors (which is apparently accomplished by this option). Shouldn't the parameter errors be influenced both by the measurement errors and the goodness of the fit of the model to this particular set of data?

The documentation of NonLinearModelFit and VarianceEstimatorFunction also do not explain how the error is propagated/what the VarianceEstimatorFunction does in this context. In the latter documentation the VarianceEstimatorFunction is defined:

With the setting VarianceEstimatorFunction->f, the variance scale is estimated by f[res,w] where res={y1-,y2-,…} is the list of residuals and w is the list of weights, as specified by the setting for the Weights option.` and they also suggest for data with measurement errors:

With VarianceEstimatorFunction->(1&) and Weights->{1/Δy12,1/Δy22,…}, Δyi is treated as the known uncertainty of measurement yi and parameter standard errors are effectively computed only from the weights.

but I still theoretically don't understand, why the second this is the option that makes NonLinearModelFit compute the parameter errors also from the measurement errors as I explained above.

I tried it out and it did not work: Also in practice the results are not convincing: I just tried out what the parameter errors are with not considering measurement errors, with weights, and with weights plus VarianceEstimatorFunction and they did not really change and even decreased in the last option, which is not plausile at all (they should increase due to added errors). In general the errors computed for "F","FTransit" should be way bigger since the data varies by about 0.2 on the y axis and the parameter error for FTransit is only 0.00326725 (have a look at the plot to see what I mean).

The parameters from the fit trying all 3 options (in the described order):

parameters with default options parameters with weighted fit parameters with weighted fit and using VarianceEstimatorFunction->(1&)

What I want to find out: So how is error propagation onto parameters done with mathematica for fitting nonlinear model functions? Why does the suggested method make sense (is it correct)? Why are the obtained errors so small and dont dont really change depending on if weights/VarianceEstimatorFunction is used? Is there any other way to solve this problem using Mathematica/another language?

TLDR: Managed to do a fit to measered data using NonLinearModelFit, now I want to propagate the measurement errors onto the parameters returned from the fit. As suggested by the documntation I tried using the options Weights -> 1/errors^2 and setting VarianceEstimatorFunction -> (1 &) which resulted in parameter errors that are propably too small. I also theoretically don't understand why this should do what I want.

How is error propagation onto fitting parameters done in mathematica, when using data with measurement errors?

Thanks in advance for helping me and all other people trying to do error analysis and fittings in mathematica.

I really appreciate your comments and the time you invest to help!

$\endgroup$
3
  • 1
    $\begingroup$ Your errors seems to be almost equal, since they are used as weights, it is somewhat expected for them to have small effect. WM estimates errors from residuals. You can also try to bootstrap your data (generate samples or use known sigmas to generate new data) to get parameter errors $\endgroup$
    – I.M.
    Apr 16, 2023 at 8:01
  • $\begingroup$ Whats "WM"? Do you have a resource that explains how the parameter errors are calculated from residuals? Which technique does mathemtica use and where is this documented? In general: Is there any additional resource that explains how mathematica functions work mathematically/which methods they use and not just how to apply them? Because the documentation doesn't really give information about the latter. $\endgroup$ Apr 16, 2023 at 16:31
  • $\begingroup$ Mathematica is closed source, exact details are not available. Assuming it does standard estimation based variance covariance matrix, you can look at scipy curve_fit source to get the idea. $\endgroup$
    – I.M.
    Apr 17, 2023 at 1:55

2 Answers 2

5
$\begingroup$

There are at least 2 issues:

  1. The predictor variable’s scale is so large that in the parameter table you’ve included, one sees the same value for the estimates of t1, t2, t3, and t4. The figure you present also has an “x-axis” labeling issue.
  2. There is some confusion as to what your data can estimate.

After digitizing the data and using 1, 2, 3, …, 95 as the predictor values I get the following:

data = {3.674, 3.670, 3.672, 3.675, 3.647, 3.659, 3.680, 3.666, 3.677,
    3.659, 3.653, 3.684, 3.670, 3.668, 3.674, 3.677, 3.656, 3.683, 
   3.657, 3.645, 3.663, 3.679, 3.638, 3.675, 3.668, 3.670, 3.650, 
   3.656, 3.663, 3.671, 3.653, 3.652, 3.644, 3.664, 3.633, 3.609, 
   3.621, 3.613, 3.605, 3.585, 3.597, 3.584, 3.581, 3.588, 3.568, 
   3.557, 3.570, 3.574, 3.563, 3.571, 3.597, 3.565, 3.587, 3.565, 
   3.573, 3.597, 3.574, 3.620, 3.630, 3.618, 3.637, 3.633, 3.649, 
   3.644, 3.666, 3.657, 3.650, 3.636, 3.666, 3.676, 3.641, 3.662, 
   3.652, 3.658, 3.644, 3.646, 3.656, 3.679, 3.653, 3.652, 3.646, 
   3.650, 3.648, 3.659, 3.657, 3.671, 3.665, 3.651, 3.654, 3.664, 
   3.659, 3.672, 3.664, 3.646, 3.659};
data = Transpose[{Range[Length[data]], data}];
ListPlot[data]

Data

Using NonlinearModelFit (and ignoring the measurement error for the moment) one gets the following:

model = Piecewise[{{F + (FTransit - F)/(t2 - t1)*(x - t1), t1 <= x < t2},
    {FTransit, t2 <= x < t3}, {FTransit + (F - FTransit)/(t4 - t3)*(x - t3), t3 <= x < t4}}, F];
nlm = NonlinearModelFit[data, model, 
  {{FTransit, 3.55}, {F, 3.66}, {t1, 30}, {t2, 45}, {t3, 55}, {t4, 65}}, x];

nlm["ParameterTable"]

Parameter table

The P-values are very small but that simply tells you that you're convinced the parameters aren't zero which in this case is not a surprise.

Are the standard errors for the parameters too small? I don't think so. Consider the estimate of $F$. We can take the response values below 25 and above 70 where the function has a mean value of $F$. The standard error for estimating $F$ from those values is found with

yF = Select[data, #[[1]] < 25 || #[[1]] > 70 &][[All, 2]];
StandardDeviation[yF]/Sqrt[Length[yF]]
(* 0.00172924 *)

That matches well with the value of 0.00151576 found in the parameter table.

Here is the data and fit:

Show[ListPlot[data], Plot[nlm[x], {x, 1, Length[data]}]]

Data and fit

Now because of the structure of your data (no repeat measurements at the same predictor value) and that the measurement error standard deviations are essentially the same value, you don't need to adjust for "measurement error" because it is already incorporated into the estimation procedure. Why do I say that?

A “model” consists of both a “fixed” effect (your Piecewise function) and an error structure. Mathematica’s NonlinearModelFit only allows the estimation of a single additive error term. A candidate model for your data and function is the following:

$$y_i=f(t_i)+\delta_i+\epsilon_i$$

where it is assumed (but needing to be checked) that $\delta_i \sim N(0,\sigma_\delta)$ and $\epsilon_i \sim N(0,\sigma)$. (I'm keeping with *Mathematica#'s habit of using the standard deviation rather than the usual variance in the definition of the normal distribution.) We have $\sigma_\delta$ representing the measurement error and $\sigma$ representing the lack-of-fit error from the proposed function. Estimating those two variance components needs to happen at the same time of fitting the model rather than adjusting sometime later.

Given that you only have a single measurement for each value of $t$ (and the measurement error standard deviation is essentially constant), if you don’t know $\sigma_\delta$ and $\sigma$, then all you can estimate is $\sigma_\delta^2+\sigma^2$ (i.e., the sum of the two variance terms). In other words, you’re fitting a model of the usual form

$$y_i=f(t_i)+e_i$$

where $e_i \sim N(0,\sqrt{\sigma_\delta^2+\sigma^2})$ so you're back to a standard model with a single error term. Therefore, when $\sigma_\delta$ and $\sigma$ are both unknown, Mathematica is giving you appropriate standard errors for the parameters and the measurement error is already incorporated.

You claim that you know $\sigma_\delta$ although you don’t give a value. (You might note a great deal of skepticism from me.) In that case you can estimate $\sigma$ and get appropriate and likely smaller standard errors for the parameters and predictions as you're better able to estimate $\sigma$ but I don't think that's necessary for this case.

$\endgroup$
4
  • $\begingroup$ Thank you for this answer! I understand what you say about the error structures in the models but I conceptionally dont understand how standard errors of parameters in non linear regression are determined in general/whats implemented in mathematica. Therefore I don't understand which results changing the error structure has. As I understand it nonlinear regression is done by minimizing an error function by an iterative algorithm optimizing the model parameters. How can this approach give estimates of the errors of the parameters? Do you have a good reccource how this is done? (+by mathemtica)? $\endgroup$ Apr 16, 2023 at 16:19
  • $\begingroup$ Comment on the data: This is from astronomical brightness observations over one night, the x axis is the Julian Date (days since −4712BC) (so t is very big and variies by about 0.5). The y axis is the relative britghtness of a star, the measurement error (std. dev.) of each datapoint is known from the shot noise of each pixel of the camera and then propagated into the total brightness of the star (y Axis). Therefore the measurement errors are known. $\endgroup$ Apr 16, 2023 at 16:46
  • $\begingroup$ I strongly recommend subtracting the minimum Julian Date, use that translated "date" in the analysis, and then use custom tick labels that represent hours (if by 0.5 you mean half of a day). If you had multiple observations at each time interval, then one could estimate both sources of error (assuming that having those two sources of error covers most of the deviation from the underlying function). $\endgroup$
    – JimB
    Apr 16, 2023 at 20:41
  • 1
    $\begingroup$ NonlinearModelFit uses by default (I think) maximum likelihood where the parameter values that maximize the log of the likelihood are found. The standard errors for the parameter estimators are found by taking the inverse of the second derivative of the log of the likelihood with respect to the parameters. My answer at mathematica.stackexchange.com/questions/164309/… gives an example of this for fitting a distribution but the general principle is the same for nonlinear regression. $\endgroup$
    – JimB
    Apr 16, 2023 at 20:47
0
$\begingroup$

Here is an illustration of VarianceEstimatorFunction & Weights options (since you do not provide data, I'm using a parabola).

Set model and objective:

(* model *)

SeedRandom[1] ;

ClearAll[model] ;
model[knobs_, x_] := 
  Block[{a, b, c}, {a, b, c} = knobs ; a*x^2 + b*x + c] ;

(* loss *)

ClearAll[objective] ;
objective[knobs_, x_, y_] := Total[(model[knobs, x] - y)^2 ];

(* knobs *)

{a, b, c} = {3.0, 2.0, 1.0} ;

(* data *)

x = Subdivide[-1.0, +1.0, 512] ;
y = model[{a, b, c}, x]  + RandomVariate[NormalDistribution[0.0, 0.25], Length[x]] ;

VarianceEstimatorFunction:

(* fit (VarianceEstimatorFunction) *)

fit = NonlinearModelFit[Transpose[{x, y}], 
   model[{pa, pb, pc}, t], {pa, pb, pc}, t, 
   VarianceEstimatorFunction -> (Mean[#^2] &)] ;
fit["ParameterErrors"]

(* standard errors (not corrected) *)

hessian = D[objective[{pa, pb, pc}, x, y], {{pa, pb, pc}, 2}] ;
residual = model[{pa, pb, pc} /. fit["BestFitParameters"] , x] - y ;
variance = Mean[residual^2] ;
Sqrt[Diagonal[Inverse[hessian] * 2 * variance]]
(* {0.03610732415967804`,0.01868201750546126`,0.01621072209504868`} *)
(* {0.036107324159678045`,0.01868201750546126`,0.01621072209504868`} *)

Weights:

(* fit (Weights) *)

fit = NonlinearModelFit[Transpose[{x, y}], 
   model[{pa, pb, pc}, t], {pa, pb, pc}, t, 
   VarianceEstimatorFunction -> (1 &), 
   Weights -> ConstantArray[1, Length[x]]/0.25^2] ;
fit["ParameterErrors"]

(* standard errors (not corrected) *)

hessian = D[objective[{pa, pb, pc}, x, y], {{pa, pb, pc}, 2}] ;
residual = model[{pa, pb, pc} /. fit["BestFitParameters"] , x] - y ;
variance =  0.25^2 ;
Sqrt[Diagonal[Inverse[hessian] * 2 * variance]]
(* {0.03687796134144208`,0.019080747061171525`,0.016556706891220286`} *)
(* {0.03687796134144208`,0.019080747061171525`,0.016556706891220286`} *)

Both:

(* fit (VarianceEstimatorFunction & Weights) *)

fit = NonlinearModelFit[Transpose[{x, y}], 
   model[{pa, pb, pc}, t], {pa, pb, pc}, t, 
   VarianceEstimatorFunction -> (Mean[#^2] &), 
   Weights -> ConstantArray[1, Length[x]]/0.25^2] ;
fit["ParameterErrors"]

(* standard errors *)

hessian = D[objective[{pa, pb, pc}, x, y], {{pa, pb, pc}, 2}] ;
residual = model[{pa, pb, pc} /. fit["BestFitParameters"] , x] - y ;
variance =  Mean[residual^2] / (1/0.25^2) ;
Sqrt[Diagonal[Inverse[hessian]  * 2 * variance]]
(* {0.00902683103991951`,0.004670504376365315`,0.00405268052376217`} *)
(* {0.009026831039919511`,0.004670504376365315`,0.00405268052376217`} *)

Bootstrap:

data = Table[{pa, pb, pc} /. With[{index = RandomChoice[Range[Length[x]], Length[x]]}, NonlinearModelFit[Transpose[{x[[index]], y[[index]]}], model[{pa, pb, pc}, t], {pa, pb, pc}, t] ["BestFitParameters"]], 1000] ;
Map[Mean, Transpose[data]]
Map[StandardDeviation, Transpose[data]]
(* {3.0000542529148317`,1.9912286713436478`,1.0066647051261475`} *)
(* {0.03548351260407338`,0.019119288196945958`,0.016017380719684282`} *)
$\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.