Astronomers often deal with a type of regression problems where the ordinate is luminosity, while the abscisca is a distance-independent observable. Such luminosity relations are some of the most useful correlations in astronomy because they enable the standard candle method to measure distances. Examples include the Leavitt (1912) law for Cepheid variables, the Tully-Fisher (1977) relation of disk galaxies, and the Phillips (1993) law of Type Ia Supernovae.

Generally speaking, the luminosity relations can be parametered as a linear relation in logarithmic scales:

\[y = \beta x + \gamma\]

where \(x\) is the luminosity predictor in logarithmic (e.g., pulsation period of a Cepheid, the rotation velocity of a galaxy, or the lightcurve decline rate of a supernova), and \(y \equiv \log L\) is the logarithmic luminosity.

To calibrate a luminosity relation is to determine the coefficients (\(\beta, \gamma\)) from a sample of standard candles where both \(x\) and \(y\) have been measured \(\{\tilde{x}_i,\tilde{y}_i\}_{i=1}^N\). However, this is not a typical linear regression problem, because the dependent variable is a composite quantity, \(y = f + d\), consisted of the logarithmic flux (\(f \equiv \log F\)) and the logarithmic distance squared (\(d \equiv \log 4\pi D^2\)). This distinction is important because observational selection affects the two quantities differently.

Once calibrated, one can use the coefficients and a measurement of \(x\) to predict the luminosity and compare it with the apparent flux to estimate its distance:

\[d = (\beta x + \gamma) - f\]

Evidently, the distance will be biased if the best estimated coefficients are biased in the calibration process. As the statistical precision of our distance measurements has reached percentage level, potential biases in these distance measurements begin to limit the accuracy of our measurements.

There are two inference biases in the calibration problem: (1) the bias introduced by sample selection such as the flux limit of a survey, and (2) the bias introduced by the measurement error of the independent variable.

The former bias is widely known among astronomers. It was first noticed by Sandage (1994), who wrote “the Hubble constant does not increase outward” when using the Tully-Fisher relation to estimate the Hubble constant. In the literature, it is called distance-dependent Malmquist bias (Butkevich et al. 2005) or simply incompleteness bias (Sandage 1988). The statistical mitigation method is to apply the observational selection function in the likelihood function (Willick 1994, Willick et al. 1997). Simply put, one needs to re-normalize the Gaussian probability density function of the deviate between the measured \(\tilde{y}\) and model prediction \(\beta x + \gamma\), because it has been truncated by the selection function.

stream data Above: Illustration of the distance-dependent Malmquist bias caused by flux limit. From left to right are the simulated data and the true correlation, the expected luminosity distribution vs. distance for sources with the same \(x\) (due to volume increase), and the distribution after the flux cut (due to observational selection).

The latter bias is widely known among statisticians, who have long realized that the regression coefficients are biased when the independent variable is measured with error. The regression coefficients are biased in opposite directions in the forward model (\(y\) vs. \(x\)) and the reverse model (\(x\) vs. \(y\)), and their means are only unbiased when \(x\) and \(y\) have comparable amount of error (i.e., \(\sigma_y \approx \beta \sigma_x\)). Mitigation methods include the bivariate correlated errors and intrinsic scatter estimator (BCES; Akritas and Bershady 1996) and the structural equation modeling with Gaussian mixture models (LINMIX; Kelly 2005). The latter work incorporated the Gaussian scatter of the independent variable and its distribution function in the likelihood function in a hierarchical structure similar to Willick (1994).

Recently, I realized that the bias of the regression coefficients is a generalized Eddington bias. The classic Eddington bias is the bias of a measurement relative to the mean of the true values that could have produced the measurement (Dyson 1926, Eddington 1940):

\[\tilde{x} - \langle x \rangle_{\tilde{x}} = -\sigma^2 \frac{d \ln p(\tilde{x})}{d\tilde{x}}\]

The bias depends on not only the measurement error but also the gradient of the distribution function \(p(\tilde{x})\). The measurement is biased because the distribution function is not uniform, causing an asymmetry between true values above the measured value and those below it. The measured value and the true value of the same quantity forms a correlation with a slope of unity, one can generalize the derivation to a correlation with a slope of \(\beta\) and obtain a similar analytical result (Fu 2025):

\[y(\tilde{x}) - \langle y \rangle_{\tilde{x}} = -\beta\sigma_x^2 \frac{d \ln p(\tilde{x})}{d\tilde{x}}\]

The generalized Eddington bias causes the mean \(y\) of the sample at fixed \(\tilde{x}\) to deviate from the correlation-predicted $y$, which consequently leads to the bias of the regression coefficients.

stream data Above: Illustration of the general Eddington bias caused by scatter and non-uniform distribution of \(x\). While the distribution from a horizontal cut is a Gaussian (top right), the one from a vertical cut is a skewed Gaussian with a mean shifted to lower luminosity (arrow in bottom right).

Although both biases are present at the same time, astronomers only deal with only one of the two biases. For example, Willick et al. (1997) provided the framework to deal with the Malmquist bias, but they provide likelihood functions for the forward model and the reverse model separately and the two models give different results when applied on the same data (meaning that the Eddington bias is untreated). On the other hand, Kelly (2007) provided the likelihood function to mitigate the Eddington bias, but they treated the dependent variable \(y\) as a whole entity instead of separating it into flux and distance; as a result, the distance-dependent Malmquist bias cannot be easily implemented. Fortunately, thanks to these previous work, it turns out to be easy to combine the two methods to mitigate both biases simultaneously. One just needs to decompose the joint probability into a series of conditional probabilities, marginalize over unknown true values, and properly incorporate selection function by renormalization (Fu 2025). In this manner, the forward and the reverse models of Willick et al. (1997) are unified into a bi-directional model that mitigate both biases by incorporating all of their contributing factors: (1) selection function of \(f\), (2) intrinsic disperison and measurement error in \(y\), (3) measurement error in \(x\), and (4) the distribution function of \(x\). The former two controls the Malmquist bias, while the latter two controls the Eddington bias. With simulated data set of galaxies following the Tully-Fisher relation, I showed that the bidirectional method can provide essentially unbiased inference of both coefficients of the luminosity relation even when some of these controlling factors are not known a priori. I hope this effort helps pave the way towards precise and accurate measurement of distance-related quantities in the future.

stream data Above: Mitigation of inference biases in regression coefficients demonstrated with simulated data. First, the posterior from simple \(\chi^2\) minimization is shown in gray (i.e., without bias correction). Next, we incorporate observational selection into the likelihood to mitigate the Malmquist bias (blue). Finally, we account for the scatter and the distribution of the luminosity predictor to correct both Malmquist and Eddington biases (red and green), yielding unbiased parameter estimates. To highlight the uncertainty in bias correction, we fix the parameters governing the Eddington bias [\(\sigma_x\), \(p(\tilde{x})\)] for the red contours, while allowing them to vary freely for the green contours. Bias uncertainty contributes minimally to the error in \(\beta\) but dominates the uncertainty in \(\gamma\).

PS. This is the longer version of the Opinion piece I wrote for the Astrostatistics Newsletter in Feb 2026.

References:

  • Akritas, M. G. and Bershady, M. A., “Linear Regression for Astronomical Data with Measurement Errors and Intrinsic Scatter”, The Astrophysical Journal, vol. 470, IOP, p. 706, 1996. doi:10.1086/177901.
  • Dyson, F., “A method for correcting series of parallax observations”, Monthly Notices of the Royal Astronomical Society, vol. 86, OUP, p. 686, 1926. doi:10.1093/mnras/86.9.686.
  • Eddington, A. S., “The correction of statistics for accidental error”, Monthly Notices of the Royal Astronomical Society, vol. 100, OUP, p. 354, 1940. doi:10.1093/mnras/100.5.354.
  • Fu, H., “Mitigating Eddington and Malmquist Biases in Latent-inclination Inference of the Tully─Fisher Relation”, The Astrophysical Journal, vol. 992, no. 2, Art. no. 214, IOP, 2025. doi:10.3847/1538-4357/adfa2c.
  • Kelly, B. C., “Some Aspects of Measurement Error in Linear Regression of Astronomical Data”, The Astrophysical Journal, vol. 665, no. 2, IOP, pp. 1489–1506, 2007. doi:10.1086/519947.
  • Sandage, A., “Cepheids as Distance Indicators When Used Near Their Detection Limit”, Publications of the Astronomical Society of the Pacific, vol. 100, IOP, p. 935, 1988. doi:10.1086/132255.
  • Sandage, A., “Bias Properties of Extragalactic Distance Indicators. I. The Hubble Constant Does Not Increase Outward”, The Astrophysical Journal, vol. 430, IOP, p. 1, 1994. doi:10.1086/174378.
  • Willick, J. A., Strauss, M. A., Dekel, A., and Kolatt, T., “Maximum Likelihood Comparisons of Tully-Fisher and Redshift Data: Constraints on Ω and Biasing”, The Astrophysical Journal, vol. 486, no. 2, IOP, pp. 629–664, 1997. doi:10.1086/304551.
  • Willick, J. A., “Statistical Bias in Distance and Peculiar Velocity Estimation. I. The ``Calibration’’ Problem”, The Astrophysical Journal Supplement Series, vol. 92, IOP, p. 1, 1994. doi:10.1086/191957.