The paper presents a Bayesian method to account for measurement errors in linear regression of astronomical data, allowing for heteroscedastic and correlated measurement errors, as well as intrinsic scatter in the regression relationship. The method is based on deriving a likelihood function for the measured data, focusing on the case where the intrinsic distribution of the independent variables can be approximated using a mixture of Gaussians. The method is extended to incorporate multiple independent variables, non-detections, and selection effects (e.g., Malmquist bias). A Gibbs sampler is described for simulating random draws from the probability distribution of the parameters given the observed data. The effectiveness of the method is demonstrated through simulations, showing that it outperforms other common estimators, 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.8 radio-quiet quasars, confirming the correlation seen by other authors. IDL routines are provided for performing the regression.The paper presents a Bayesian method to account for measurement errors in linear regression of astronomical data, allowing for heteroscedastic and correlated measurement errors, as well as intrinsic scatter in the regression relationship. The method is based on deriving a likelihood function for the measured data, focusing on the case where the intrinsic distribution of the independent variables can be approximated using a mixture of Gaussians. The method is extended to incorporate multiple independent variables, non-detections, and selection effects (e.g., Malmquist bias). A Gibbs sampler is described for simulating random draws from the probability distribution of the parameters given the observed data. The effectiveness of the method is demonstrated through simulations, showing that it outperforms other common estimators, 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.8 radio-quiet quasars, confirming the correlation seen by other authors. IDL routines are provided for performing the regression.