How to do curve fitting on a data set with errors associated to each data point (y ± Δ y)

I have a data set in the form x (independent variable), y ± Δy (dependent variable). As you can see there's an error associated with each data point in dependent variable. I want to do curve fitting for this data set and can you help me with how to do that.

I know how to use curve fitting tool box to fit a curve for a data set in the form x,y. But I don't know how to use matlab to do curve fitting when there's an error associated with each data point.

[Star Strider]

Use the errors as weights for each data point. The typical strategy is to use inverse-variance weighting, so you will need to convert your errors to variances, then use 1/variance for each data point as the weight for that point.

I do not have the Curve Fitting Toolbox, so you will need to discover how to do this yourself.

[SWATI BHARGAVA]

I have exactly same question as posted and I'm also using curve fitting toolbox and need to fit each data point which has error bars associated with it. If I do NOT use error bars with data, I'm able to fit the data using Smoothing Spline fit in curve fitting tootlbox. But I need to know how to do this when error bars are added to each y data point.

Is the approach mentioned by Star Strider above correct? Since curve fitting toolbox allows us to choose the weights, so this approach is pretty doable , if correct. Anyone , any leads?

[SWATI BHARGAVA]

• Once above is done, I finally want to find the standard deviation of the data in fitted curve as well.

Answer 1

[John D'Errico]

Since you have explicit bounds on y +/- deltay, you cannot just use simple weights. As well, you cannot use the curve fitting toolbox. Instead, you can set this up as an fmincon problem, with nonlinear constraints. Each data point will provide two nonlinear constraints, an upper and lower bound for the curve at each given point. This is much simpler of course if the problem is a linear one, since then you could use lsqlin instead with simple linear inequality constraints.

So, is your model linear in the parameters? I.e., any polynomial model, is linear in the parameters. So models like y = a0+a1*x, a0+a1*x+a2*x^2, etc. Each of those models are linear in the parameters, though they need not be linear in x.

Alternatively, a model like y=a0+a1*sin(a2*x) is linear in two of the parameters, but nonlinear in a2. This would take somewhat more effort to solve than the simple polynomial models I described, but it is still solvable. Other models, thus perhaps a1*exp(a2*x) wuld also have a solution method that would be pehaps most appropriate.

Now, perhaps you have no a-priori choice of model? In that case, a regression spline model probably makes sense. And I would point out that my SLM toolbox actually allows the suer to provide an explicit set of upper and lower bounds that can be used to constrain the curve at each data point.

The point is, the most appropriate solution method will depend on the model, but I don't know what model you have.

Answer 2

[Ankit Labh]

I have exactly the same problem. I wanted to fit gaussian with all y values as y +- Delta_y. Considering the errors associated to y value, I still do not get the answer to do the fitting. If you get any answers, do post here.

Thank you very much

Answer 3

[Boaz Zick]

I just found something that can maybe help you.

https://www.youtube.com/watch?v=99hymvgpkiI&ab_channel=SteveSpicklemire

watch from 10:10.

Monte Carlo method

the idea is to add noise to your data to simulate alot of measurments that could have been made with your error in measurment, and then curve fit each one of those measurments, and check the mean and standart deviation of your coefficients.

here is an example for a linear fit, in order to fit other curves change fitlm and add coefficients as necessary.

%%

rep=500;
a_lst=[a_coefficient]; %first value can be your curve fitting coefficients from ideal data points
b_lst=[]; %or left blank
for i = 1:rep
    xd=X+X_error.*normrnd(0,1,1,length(x));
    yd=Y+Y_error.*normrnd(0,1,1,length(y));
    d_linear=fitlm(xd,yd);
    bd_coefficient = d_linear.Coefficients.Estimate(1);
    b_lst=[b_lst,bd_coefficient];
    ad_coefficient = d_linear.Coefficients.Estimate(2);
    a_lst=[a_lst,ad_coefficient];
    i
end
histogram (a_lst,100);
a=mean(a_lst);
a_delta=std(a_lst);

b=mean(b_lst);
b_delta=std(b_lst);