Error in linear regression with predefined error in y

I'm fitting y=ax+b with polyfit. x has no errors, but every component y_i has an error equal to error_i = C_i*y_i. (So this is correlated right?) How do I determine the error in the slope a?

I've been thinking about not using polyfit and minimazing S = sum(w_i * ( y_i - fit_i)^2) myself. With w_i = 1/error_i^2. But I have no idea how this minimizing can be done.

Answer 1

[Tom Lane]

Take a look at the lscov function and see if it does what you need.

[Arjen]

it takes the weight factors and give the standard error i need. But only for Ax = b. I would like one which includes an intercept

[Tom Lane]

Just include an column of ones. Using the "y=X*b+error" notation, if you have a vector x and you want to fit a parabola, X = [ones(size(x)),x,x.^2].

[Arjen]

im not quiet sure if I understand this correctly. What's exactly the purpose of these ones? And is this intercept then also used in the calculation of the slope of the fit?

[Tom Lane]

You wanted an intercept. The equation a*1+b*x defines the intercept as "a." If every row of X has a 1 and an x value, you'll be fitting this equation with an intercept as the first element of the coefficient vector and the slope as the second element. The slope will then be computed for a general line, rather than one constrained to have an intercept equal to zero.

[Arjen]

ah yes, of course. Thanks a lot!!

Answer 2

[Wayne King]

polyfit returns a least-squares fit, but not with weights as you suggest. Do you have the Statistics Toolbox? If so consider, robustfit.m or LinearModel.fit, which has options for robust fitting.

Also, perhaps a simple first-order linear model is not adequate for your data?

[Arjen]

I've browsed a lot to find a simple linear fitting method. But I didnt find anything similar to a linear regression with errors on y... which actually suprises my a lot because this is a fairly easy problem... At least i thought so.

The robustfit is giving me what i need, but the method is rather frustrating... It asks for a function instead of a vector. I'm taking the stats.se as an error for the slope, but it 1e20 times bigger than the slope itself. And the fit is pretty good. Something is going wrong here...