regression with error bars

  • Follow


Hello

I have two columns of data y (with errors Err+ and Err- for each data
point) and x (with Err). I would like to find the linear regression
coefficients A and B (in y=Ax+B) considering the errors in my data
point. Is their a function in IDL which would do the trick? and the
regression has to do with the errors in y only but not in x, right?

Thank you for your help.
0
Reply Meagan 5/4/2010 7:57:04 PM

On May 4, 3:57=A0pm, Meagan Adams <meagan....@gmail.com> wrote:
> Hello
>
> I have two columns of data y (with errors Err+ and Err- for each data
> point) and x (with Err). I would like to find the linear regression
> coefficients A and B (in y=3DAx+B) considering the errors in my data
> point. Is their a function in IDL which would do the trick? and the
> regression has to do with the errors in y only but not in x, right?

I don't believe there is a built-in function in IDL which can do this.

You can use MPFIT.  The trick is that you will need to be able to form
your own scaled residual.  Normally you would have a user function
something like this,
FUNCTION MYFUNCT, P, X=3DX, Y=3DY, ERROR=3DERROR
  MODEL =3D F(X,P)  ;; User model
  RETURN,  (Y - MODEL)/ERROR
END

but now you will need to reformulate as,
FUNCTION MYFUNCT, P, X=3DX, Y=3DY, ERROR_PLUS=3DERROR_PLUS,
ERROR_MINUS=3DERROR_MINUS
  MODEL =3D F(X,P)  ;; User model
  RESID_PLUS  =3D (Y - MODEL) / ERROR_PLUS   ;; Residual for positive
values
  RESID_MINUS =3D (Y - MODEL) / ERROR_MINUS  ;; Residual for negative
values
  ;; Combine the two kinds of residuals
  RESID =3D (RESID_PLUS GT 0)*RESID_PLUS + (RESID_MINUS LT
0)*RESID_MINUS
  RETURN, RESID
END

Happy fitting,
Craig



0
Reply Craig 5/5/2010 2:19:26 AM


Hello,

Thank you for your reply and examples. I am a newbie to IDL and I've never used MPFIT. I read the tutorial from this website:
http://cow.physics.wisc.edu/~craigm/idl/mpfittut.html

and I wrote a pro to fit a line to x and y if we had errors in y only (+/-):

pro test

xval = [1,2,3,4,5]
yval = [1.2,1.9, 3.1, 4.5, 5]
errval = [0.1,0.001,0.1,0.05,0.2]

expr = 'P[0]+P[1]*X'  ;line y=A+Bx

result = MPFITEXPR(expr, xval, yval,errval)

plot, xval, yval, psym=2
oplot, xval, result[0]+result[1]*xval

l=linfit(xval, yval) ;no errors considered
oplot, xval, l[0]+l[1]*xval, linestyle=1

end


Are the error values (errval) used in this way in the function MPFITEXPR? And if not, after I calculate the residuals using your method, how do I feed it into the function MPFITFUN ?

Thank you!!
0
Reply Meagan 5/22/2010 1:34:26 PM

On May 22, 9:34=A0am, Meagan A. <u...@compgroups.net/> wrote:
> Hello,
>
> Thank you for your reply and examples. I am a newbie to IDL and I've neve=
r used MPFIT. I read the tutorial from this website:http://cow.physics.wisc=
..edu/~craigm/idl/mpfittut.html
>
> and I wrote a pro to fit a line to x and y if we had errors in y only (+/=
-):
>
> pro test
>
> xval =3D [1,2,3,4,5]
> yval =3D [1.2,1.9, 3.1, 4.5, 5]
> errval =3D [0.1,0.001,0.1,0.05,0.2]
>
> expr =3D 'P[0]+P[1]*X' =A0;line y=3DA+Bx
>
> result =3D MPFITEXPR(expr, xval, yval,errval)
>
> plot, xval, yval, psym=3D2
> oplot, xval, result[0]+result[1]*xval
>
> l=3Dlinfit(xval, yval) ;no errors considered
> oplot, xval, l[0]+l[1]*xval, linestyle=3D1
>
> end
>
> Are the error values (errval) used in this way in the function MPFITEXPR?=
 And if not, after I calculate the residuals using your method, how do I fe=
ed it into the function MPFITFUN ?

There are several MPFIT* functions.  You can't use MPFITFUN or
MPFITEXPR because you are changing the standard definition of
"residual" because you have different + and - error bars.  Instead you
will need to use the core engine MPFIT(), and you will need to write a
user function like I described in my previous post.

Craig
0
Reply Craig 5/22/2010 3:02:16 PM

Perfect!! Thank you!!!!
0
Reply Meagan 5/22/2010 4:29:27 PM

One more question, please.

In:
RESID =3D (RESID_PLUS GT 0)*RESID_PLUS + (RESID_MINUS LT 
0)*RESID_MINUS 

Shouldn't we subtract the two resids instead of adding them? Thank you again!
0
Reply Meagan 5/22/2010 4:44:56 PM

On May 22, 12:44=A0pm, Meagan A. <u...@compgroups.net/> wrote:
> One more question, please.
>
> In:
> RESID =3D (RESID_PLUS GT 0)*RESID_PLUS + (RESID_MINUS LT
> 0)*RESID_MINUS
>
> Shouldn't we subtract the two resids instead of adding them? Thank you ag=
ain!

Not really.  The "addition" there is an IDL vectorized short-hand
notation for, "if residuals are greater than zero use RESID_PLUS;
otherwise use RESID_MINUS".   Only one clause of expression is non-
zero, so you are never really adding or subtracting residual values.

Craig
0
Reply Craig 5/22/2010 5:53:34 PM

9 Replies
717 Views

(page loaded in 0.074 seconds)

Similiar Articles:













7/22/2012 10:55:46 AM


Reply: