(This is a note personal only! )

There is a powerful tool dealing with data containing both “X” and “Y” errors to make a lease-square fitting in IDL, that is mpfit. This routine allows you to use a self-defined function to do the fitting. Also, you can set up some conditions such as the parameters’ limit or fixed parameters etc. This is fabulous! So beginning with the actual astronomical data have errors both in “X” and “Y”, you need to use a function called “linfitex” (you can see more info from the mpfit website) to both define the function you want to fit and calculate the least-square needed in the fitting. Here is an example of the function:

FUCNTION myfunct, p, x=x, y=y \$
     , sigma_x=sigma_x, sigma_y=sigma_y, _EXTRA=extra
  a = p[0]
  b = p[1]
  f = a*alog10(X)+b
;#########calculate the derivative of the function#########
;##put into resid as (y-f)/sqrt(sigma_y^2+(y'*sigma_x)^2)## 
  resid = (y - f)/sqrt(sigma_y^2 + (a/(x*alog(10))*sigma_x)^2) 
RETURN, resid

In the code above I just defined a function as

Then you can read your data. and when you want to fit the data, just using mpfit as:

p      = [1,1,1]
p[1]   = 1         
p[2]   = 2
p[3]   = 3 ;initial guess of the parameters
param  = MPFIT('MYFUNCT', params,\$
      FUNCTARGS={X:x,Y:y,SIGMA_X:x_er,SIGMA_Y:y_er}, PERROR=p_er)

The FUNCTARGS is the place where your data go into the function. And then you can find the results of the parameters from the fitting by print “param” and also the 1σ error from printing p_er. If you only have the error on “Y”, you can use a function like this:

FUNCTION myfunct, p, x=x, y=y, er=er
  a = p[0]     
  b = p[1]
  c = p[2]
  f = a*alog10(b)+c
RETURN, (y-f)/er
END

And then using the same way described above to call “mpfit” and get the results and errors.

Here is also a PYTHON version of this procedure.