This paper presents a Bayesian method for handling measurement errors in linear regression of astronomical data. The method accounts for heteroscedastic and possibly correlated measurement errors, as well as intrinsic scatter in the regression relationship. It is based on deriving a likelihood function for the measured data, with a focus on cases where the intrinsic distribution of the independent variables can be approximated using a mixture of Gaussians. The method is generalized to incorporate multiple independent variables, non-detections, and selection effects. A Gibbs sampler is described for simulating random draws from the probability distribution of the parameters, given the observed data. Simulations show that the Gaussian mixture model outperforms other common estimators and can effectively give constraints on the regression parameters, even when measurement errors dominate the observed scatter, source detection fraction is low, or the intrinsic distribution of the independent variables is not a mixture of Gaussians. The method is applied to fit the X-ray spectral slope as a function of Eddington ratio using a sample of 39 z < 0.83 radio-quiet quasars. The results confirm the correlation between the X-ray spectral slope and the Eddington ratio, where the X-ray spectral slope softens as the Eddington ratio increases. IDL routines are provided for performing the regression. The paper discusses the effects of measurement error on the estimates of the regression slope and correlation coefficient, and describes the statistical model and methods for incorporating selection effects and non-detections. The computational methods include a Bayesian approach for computing estimates of the regression parameters and their uncertainties, using Markov Chain Monte Carlo (MCMC) methods. The paper also describes the prior density and Gibbs sampler for sampling from the posterior distribution.This paper presents a Bayesian method for handling measurement errors in linear regression of astronomical data. The method accounts for heteroscedastic and possibly correlated measurement errors, as well as intrinsic scatter in the regression relationship. It is based on deriving a likelihood function for the measured data, with a focus on cases where the intrinsic distribution of the independent variables can be approximated using a mixture of Gaussians. The method is generalized to incorporate multiple independent variables, non-detections, and selection effects. A Gibbs sampler is described for simulating random draws from the probability distribution of the parameters, given the observed data. Simulations show that the Gaussian mixture model outperforms other common estimators and can effectively give constraints on the regression parameters, even when measurement errors dominate the observed scatter, source detection fraction is low, or the intrinsic distribution of the independent variables is not a mixture of Gaussians. The method is applied to fit the X-ray spectral slope as a function of Eddington ratio using a sample of 39 z < 0.83 radio-quiet quasars. The results confirm the correlation between the X-ray spectral slope and the Eddington ratio, where the X-ray spectral slope softens as the Eddington ratio increases. IDL routines are provided for performing the regression. The paper discusses the effects of measurement error on the estimates of the regression slope and correlation coefficient, and describes the statistical model and methods for incorporating selection effects and non-detections. The computational methods include a Bayesian approach for computing estimates of the regression parameters and their uncertainties, using Markov Chain Monte Carlo (MCMC) methods. The paper also describes the prior density and Gibbs sampler for sampling from the posterior distribution.