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)
|