A geometric approach to non-linear correlations with intrinsic scatter

TL;DR

This paper introduces a Riemannian geometric model for non-linear correlations with intrinsic scatter, enabling more accurate Bayesian estimation and application to astrophysical data, notably revising the slope of the black hole mass-velocity dispersion relation.

Contribution

It presents a novel geometric framework for modeling non-linear correlations with intrinsic scatter, incorporating Bayesian methods and analytic likelihoods, applicable to complex datasets.

Findings

The model is symmetric and invariant under coordinate transformations.

It effectively incorporates censored and truncated data with measurement errors.

Applied to astrophysical data, it suggests a higher slope (~6) for the black hole mass-velocity dispersion relation.

Abstract

We propose a new mathematical model for $n-k$-dimensional non-linear correlations with intrinsic scatter in $n$-dimensional data. The model is based on Riemannian geometry, and is naturally symmetric with respect to the measured variables and invariant under coordinate transformations. We combine the model with a Bayesian approach for estimating the parameters of the correlation relation and the intrinsic scatter. A side benefit of the approach is that censored and truncated datasets and independent, arbitrary measurement errors can be incorporated. We also derive analytic likelihoods for the typical astrophysical use case of linear relations in $n$-dimensional Euclidean space. We pay particular attention to the case of linear regression in two dimensions, and compare our results to existing methods. Finally, we apply our methodology to the well-known $M_\text{BH}$-$\sigma$ correlation between the mass of a supermassive black hole in the centre of a galactic bulge and the corresponding bulge velocity dispersion. The main result of our analysis is that the most likely slope of this correlation is $\sim 6$ for the datasets used, rather than the values in the range $\sim 4$-$5$ typically quoted in the literature for these data.

Tables

Table 1. Table 1: Slopes of linear relations fit with the maximum a posteriori estimate of equation ( 37 ) (MAP) and the linmix_err and FITEXY methods.

Slope of data	$0.1$	$1$	$10$
Method	Slope of fit $\pm$ 1-$\sigma$ uncertainties
	$\alpha =0$
linmix_err fwd	$0.09\pm 0.02$	$0.99\pm 0.04$	$1.67\pm 0.70$
linmix_err inv	$0.42\pm 0.12$	$1.08\pm 0.05$	$15.4\pm 5.99$
FITEXY fwd	$0.09\pm 0.02$	$0.99\pm 0.04$	$1.68\pm 0.39$
FITEXY inv	$0.41\pm 0.06$	$1.08\pm 0.04$	$15.5\pm 6.10$
Geo-MAP fwd	$0.10\pm 0.03$	$1.04\pm 0.05$	$11.8\pm 2.80$
Geo-MAP inv	$0.10\pm 0.03$	$1.04\pm 0.04$	$11.8\pm 2.79$
	$\alpha =10$
linmix_err fwd	$0.04\pm 0.02$	$0.93\pm 0.04$	$1.79\pm 0.57$
linmix_err inv	$0.51\pm 0.27$	$1.01\pm 0.04$	$10.9\pm 3.31$
FITEXY fwd	$0.04\pm 0.02$	$0.93\pm 0.04$	$1.80\pm 0.35$
FITEXY inv	$0.52\pm 0.13$	$1.01\pm 0.04$	$10.9\pm 3.24$
Geo-MAP fwd	$0.04\pm 0.02$	$0.98\pm 0.04$	$8.66\pm 1.68$
Geo-MAP inv	$0.04\pm 0.02$	$0.98\pm 0.04$	$8.66\pm 1.68$

Table 2. Table 2: M BH subscript 𝑀 BH M_{\text{BH}} – σ 𝜎 \sigma relations log 10 ⁡ ( M BH / M ⊙ ) = α + β ​ log 10 ⁡ ( σ / 200 ​ km ​ s − 1 ) subscript 10 subscript 𝑀 BH subscript M direct-product 𝛼 𝛽 subscript 10 𝜎 200 km superscript s 1 \log_{10}(M_{\text{BH}}/\text{M}_{\odot})=\alpha+\beta\log_{10}(\sigma/200\,$\mathrm{k}\mathrm{m}\,\mathrm{s}^{-1}$) derived using the datasets in Gültekin et al. ( 2009 ) , McConnell et al. ( 2011 ) , Graham et al. ( 2011 ) and van den Bosch ( 2016 ) , using methods FITEXY , linmix_err and the Geometric method of this paper. For the Geometric method, results with both bootstrapped (Geo-MAP) and MCMC-sampled (Geo-MCMC) maximum a posteriori estimate are shown. The results for the FITEXY and linmix_err methods for the first three datasets are from Park et al. ( 2012 ) .

Method	Intercept $\alpha$	Slope $\beta$	${\sigma }_{\text{int},M_{\text{BH}}}$a	${\sigma }_{\text{int}}$a
Gültekin et al. (2009) data
FITEXY fwd	$8.19\pm 0.06$	$4.06\pm 0.32$	$0.39\pm 0.06$	$0.093$
FITEXY inv	$8.21\pm 0.07$	$5.35\pm 0.66$	$0.45\pm 0.09$	$0.083$
linmix_err fwd	$8.19\pm 0.07$	$4.04\pm 0.40$	$0.42\pm 0.05$	$0.101$
linmix_err inv	$8.21\pm 0.08$	$5.44\pm 0.56$	$0.49\pm 0.09$	$0.089$
Geo-MAP fwd	$8.22\pm 0.07$	$5.42\pm 0.72$	$0.43\pm 0.08$	$0.076\pm 0.008$
Geo-MCMC fwd	${8.21}_{-0.14}^{+0.14}$	${5.41}_{-0.79}^{+1.10}$	${0.44}_{-0.09}^{+0.12}$	${0.080}_{-0.016}^{+0.022}$
McConnell et al. (2011) data
FITEXY fwd	$8.28\pm 0.06$	$5.07\pm 0.36$	$0.43\pm 0.05$	$0.083$
FITEXY inv	$8.32\pm 0.06$	$6.29\pm 0.49$	$0.47\pm 0.06$	$0.073$
linmix_err fwd	$8.27\pm 0.06$	$5.06\pm 0.36$	$0.44\pm 0.05$	$0.085$
linmix_err inv	$8.32\pm 0.07$	$6.31\pm 0.46$	$0.49\pm 0.07$	$0.077$
Geo-MAP fwd	$8.32\pm 0.05$	$6.26\pm 0.51$	$0.45\pm 0.06$	$0.070\pm 0.008$
Geo-MCMC fwd	${8.32}_{-0.12}^{+0.13}$	${6.31}_{-0.71}^{+0.90}$	${0.46}_{-0.08}^{+0.11}$	${0.071}_{-0.013}^{+0.017}$
Graham et al. (2011) data
FITEXY fwd	$8.15\pm 0.05$	$5.08\pm 0.34$	$0.31\pm 0.04$	$0.060$
FITEXY inv	$8.16\pm 0.05$	$5.84\pm 0.42$	$0.33\pm 0.05$	$0.056$
linmix_err fwd	$8.15\pm 0.05$	$5.08\pm 0.36$	$0.31\pm 0.05$	$0.060$
linmix_err inv	$8.17\pm 0.06$	$5.85\pm 0.42$	$0.34\pm 0.06$	$0.057$
Geo-MAP fwd	$8.17\pm 0.04$	$6.00\pm 0.41$	$0.31\pm 0.05$	$0.051\pm 0.006$
Geo-MCMC fwd	${8.17}_{-0.11}^{+0.11}$	${5.98}_{-0.67}^{+0.87}$	${0.32}_{-0.10}^{+0.11}$	${0.052}_{-0.016}^{+0.018}$
van den Bosch (2016) data
linmix_err fwdb	$8.32\pm 0.04$	$5.30\pm 0.22$	$0.49\pm 0.03$	$0.091$
Geo-MAP fwd	$8.42\pm 0.04$	$5.90\pm 0.26$	$0.55\pm 0.06$	$0.091\pm 0.008$
Geo-MCMC fwd	${8.44}_{-0.09}^{+0.09}$	${6.15}_{-0.40}^{+0.42}$	${0.54}_{-0.06}^{+0.07}$	${0.087}_{-0.009}^{+0.011}$
van den Bosch (2016) data, no upper limits
FITEXY fwd	$8.35\pm 0.04$	$4.91\pm 0.23$	$0.48$c	$0.096$
FITEXY inv	$8.42\pm 0.05$	$6.51\pm 0.30$	$0.55$c	$0.084$
linmix_err fwd	$8.34\pm 0.04$	$4.92\pm 0.25$	$0.49\pm 0.03$	$0.098$
linmix_err inv	$8.43\pm 0.05$	$6.57\pm 0.32$	$0.56\pm 0.05$	$0.084$
Geo-MAP fwd	$8.43\pm 0.04$	$6.70\pm 0.40$	$0.54\pm 0.07$	$0.080\pm 0.008$
Geo-MCMC fwd	${8.43}_{-0.09}^{+0.10}$	${6.68}_{-0.55}^{+0.67}$	${0.55}_{-0.06}^{+0.08}$	${0.082}_{-0.010}^{+0.011}$

Equations

$$\begin{array}{rl}S& =\{m\in M|f(m)=0\}\\ & =\bigcap _{i=1,\dots ,k}{S}_i=\bigcap _{i=1,\dots ,k}\{m\in M|f_i(m)=0\}\end{array}$$

$$NS=\bigcup _{p\in S}{N}_{p}S=\bigcup _{p\in S}\{z\in T_{p}M|g(z,w)=0,\forall w\in T_{p}S\}$$

$$m_x$$

$$m_y$$

$$p_{\text{int}}(\mathbf{x})=\sum _{i\in I}{p}_{\text{int}}({\mathbf{\tau }}_i(\mathbf{x}),{\mathbf{\nu }}_i(\mathbf{x})),$$

$$p_{\text{int}}(\mathbf{\nu };\mathbf{\Sigma })=\frac{\exp \left[-\frac{1}{2}{\mathbf{\nu }}^T{\mathbf{\Sigma }}^{-1}\mathbf{\nu }\right]}{\sqrt{(2\pi {)}^k\det (\mathbf{\Sigma })}},$$

$$p_{\text{int}}(\nu ;\sigma )=\frac{1}{\sqrt{2\pi {\sigma }^2}}\exp \left(-\frac{{\nu }^2}{2{\sigma }^2}\right),$$

$$p_{\text{int}}(x,y)=\frac{1}{\sqrt{2\pi {\sigma }^2}}\exp \left[-\frac{(y-\alpha -\beta x{)}^2}{2{\sigma }^2(1+{\beta }^2)}\right],$$

$$p_x(x)={\int }_{\mathbb{R}}{p}_{\text{int}}(x,y)\mathrm{d}y=\sqrt{1+{\beta }^2}$$

$$p_y(y)={\int }_{\mathbb{R}}{p}_{\text{int}}(x,y)\mathrm{d}x=\sqrt{1+{\beta }^2}/|\beta |,$$

$$\begin{array}{rl}& \mathcal{L}(\mathbf{\theta },\mathbf{\phi }|\mathcal{D})=\\ & \prod _{i=1}^{n_d}{\int }_M{h}_i(\mathbf{x}(\mathbf{\tau },\nu ;\mathbf{\theta }))p_{\text{det}}(\mathbf{x}(\mathbf{\tau },\nu ;\mathbf{\theta }))p_{\text{int}}(\mathbf{\tau },\nu ;\mathbf{\phi }){\text{vol}}_{{\varphi }_f}=\\ & \prod _{i=1}^{n_d}{\int }_M{h}_i(\mathbf{x})p_{\text{det}}(\mathbf{x})\sum _{j\in I}{p}_{\text{int}}({\mathbf{\tau }}_j(\mathbf{x};\mathbf{\theta }),{\mathbf{\nu }}_j(\mathbf{x};\mathbf{\theta });\mathbf{\phi }){\text{vol}}_{{\varphi }_x},\end{array}$$

$$p_{\text{post}}(\mathbf{\theta },\mathbf{\phi })=\frac{\mathcal{L}(\mathbf{\theta },\mathbf{\phi }|\mathcal{D}){\pi }_{\mathbf{\theta },\mathbf{\phi }}(\mathbf{\theta },\mathbf{\phi })}{\iint \mathcal{L}({\mathbf{\theta }}^{\prime },{\mathbf{\phi }}^{\prime }|\mathcal{D}){\pi }_{\mathbf{\theta },\mathbf{\phi }}({\mathbf{\theta }}^{\prime },{\mathbf{\phi }}^{\prime })\mathrm{d}{\mathbf{\theta }}^{\prime }\mathrm{d}{\mathbf{\phi }}^{\prime }}.$$

$$f(\mathbf{x};\mathbf{\theta })=(\mathbf{x}-\mathbf{p}{)}^T\mathbf{p}={\mathbf{x}}^T\mathbf{p}-{\Vert \mathbf{p}\Vert }^2=0,$$

$$x_1=a_1+\sum _{j=2}^n{a}_i{x}_i.$$

$$a_1$$

$$p_1$$

$$h(\mathbf{x};\mathbf{m},\mathbf{\Sigma })=\frac{\exp \left[-\frac{1}{2}(\mathbf{x}-\mathbf{m}{)}^T{\mathbf{\Sigma }}^{-1}(\mathbf{x}-\mathbf{m})\right]}{\sqrt{(2\pi {)}^n\det (\mathbf{\Sigma })}},$$

$$\text{R}(\mathbf{p})=I-\hat{\mathbf{w}}{\hat{\mathbf{w}}}^T-{\hat{\mathbf{e}}}_n{\hat{\mathbf{e}}}_n^T+\left(\hat{\mathbf{w}}{\hat{\mathbf{e}}}_n\right){\text{R}}_2(\theta ){\left(\hat{\mathbf{w}}{\hat{\mathbf{e}}}_n\right)}^T,$$

$${\text{R}}_2(\theta )=\left(\begin{array}{cc}\cos \theta & \sin \theta \\ -\sin \theta & \cos \theta \end{array}\right)$$

$$\mathcal{L}(\mathbf{p},\sigma |\mathbf{m},\mathbf{\Sigma })=\frac{1}{\sqrt{2\pi ({\tilde{\Sigma }}_{nn}+{\sigma }^2)}}\exp \left[\frac{-{\tilde{m}}_n^2}{2({\tilde{\Sigma }}_{nn}+{\sigma }^2)}\right],$$

$$\tilde{\mathbf{m}}$$

$$\tilde{\mathbf{\Sigma }}$$

$$\tilde{\mathbf{m}}$$

$$\tilde{\Sigma }$$

$$\begin{array}{rl}({\varphi }_f\circ {\varphi }_x^{-1})(\mathbf{x})& ={\text{R}}_{n-k+1}^T\cdots {\text{R}}_n^T(\mathbf{x}-{\mathbf{p}}_n)\\ & -{\text{R}}_{n-k+1}^T\cdots {\text{R}}_{n-1}^T{\mathbf{p}}_{n-1}-\cdots \\ & -{\text{R}}_{n-k+1}^T{\mathbf{p}}_{n-k+1}\\ & ={\text{R}}^T(\mathbf{x}-\mathbf{p})\\ & =(\mathbf{\tau },\mathbf{\nu })\in {\mathbb{R}}^{n-k}\times {\mathbb{R}}^k,\end{array}$$

$${\text{R}}^T={\text{R}}_{n-k+1}^T\cdots {\text{R}}_n^T$$

$$\mathbf{p}={\mathbf{p}}_n+{\text{R}}_n{\mathbf{p}}_{n-1}+\cdots +{\text{R}}_n\cdots {\text{R}}_{n-k+2}{\mathbf{p}}_{n-k+1}.$$

$$f(\mathbf{x};\mathbf{\theta })=[{\text{R}}^T(\mathbf{x}-\mathbf{p}){]}_{n-k+1,\dots ,n}=0,$$

$$\mathcal{L}(\mathbf{\theta },{\mathbf{\Sigma }}_{\text{int}}|\mathbf{m},\mathbf{\Sigma })=\frac{\exp \left[-\frac{1}{2}{\tilde{\mathbf{m}}}_{\nu }^T({\tilde{\mathbf{\Sigma }}}_{\nu \nu }+{\mathbf{\Sigma }}_{\text{int}}{)}^{-1}{\tilde{\mathbf{m}}}_{\nu }\right]}{\sqrt{(2\pi {)}^k\det \left({\tilde{\mathbf{\Sigma }}}_{\nu \nu }+{\mathbf{\Sigma }}_{\text{int}}\right)}},$$

$$\tilde{\mathbf{m}}$$

$$\tilde{\mathbf{\Sigma }}$$

$$f(x,y;\mathbf{\theta })=y-\beta x-\alpha =0,$$

$$\begin{array}{rl}& h(x,y;x_0,y_0,{\sigma }_x,{\sigma }_y,\rho )=\\ & \frac{\exp \left\{-\frac{1}{2(1-{\rho }^2)}\left[\frac{(x-x_0{)}^2}{{\sigma }_x^2}+\frac{(y-y_0{)}^2}{{\sigma }_y^2}-\frac{2\rho (x-x_0)(y-y_0)}{{\sigma }_x{\sigma }_y}\right]\right\}}{2\pi {\sigma }_x{\sigma }_y\sqrt{1-{\rho }^2}},\end{array}$$

Full text

††pubyear: 2017††pagerange: A geometric approach to non-linear correlations with intrinsic scatter–B

1 Introduction

An important question in all the sciences is whether different measured observables seem to be connected by some form of mathematical relation, or whether they seem to be completely independent. If a relation is suspected, it then becomes important to estimate the type of relation connecting the measurements and to estimate its parameters.

Historically, these problems have been most often approached by the use of various correlation coefficients and linear and non-linear regression. Linear regression in particular has been the subject of lively debates and much research over the past two centuries or so, ever since first derived in the guise of least-squares regression by Gauss (or possibly Legendre, see Stigler 1981). The discussion has included such points of contention as how to treat the variables symmetrically – that is how to avoid dividing the variables into dependent and independent variables – and whether this is necessary (Pearson, 1901; Boggs et al., 1987; Isobe et al., 1990; Feigelson & Babu, 1992; Robotham & Obreschkow, 2015).

This issue is particularly important when investigating the hypothesis that the observables obey some mathematical relation. By definition, such a relation as a geometric structure (e.g. a line, a plane, or some more complicated subset) is unique, even though it may have many equivalent formulations. As such, any method used for estimating such a relation should eventually yield the same structure even if the variables are permuted or a different variable is chosen as the dependent variable. This requires that the method must involve all observed variables in an essentially symmetric way. However, this is not true for least squares regression and many other similar approaches as well (see Section 3.3.1). This fact has been well known for a long time (see e.g. Pearson, 1901; Isobe et al., 1990), but it appears somewhat underappreciated. Since the focus of this paper is on formulating and estimating relations between observables, we will pay close attention to this symmetry invariance throughout.

Another recurrent theme is the question of how to treat data subject to censoring (i.e. lower and upper limits), truncation (non-detections) and heteroscedastic, independent and correlated errors in general, and finally, how to incorporate intrinsic scatter, or intrinsic uncertainty in the regression hyperplane (see e.g. Kelly 2007, Hogg et al. 2010 and Robotham & Obreschkow 2015 and the references therein).

Intrinsic scatter can be intuitively understood to mean that the data probability distribution can be non-zero also in the neighbourhood of the subset defined by the relation, that is we also allow for data that may be somewhat off the relation and thus not exactly satisfying the equation(s) defining the relation. For example, in the case of a line in two dimensions, intrinsic scatter would manifest as a data distribution that is not a linear delta function ridge, but a ‘fuzzy’ line instead. Despite this intuitive obviousness, intrinsic scatter has a history of various somewhat non-rigorous definitions, such as through an additional additive component in measurement errors. In this paper, we will present a more rigorous definition of intrinsic scatter which is also usable in non-Euclidean contexts.

Specific examples of papers addressing some of the problems above include Pearson (1901), who introduced a least squares method to fit lines and hyperplanes based on minimizing residuals orthogonal to the regression plane (OR, orthogonal regression). Later, algorithms were developed to solve this problem for non-linear relations as well (see e.g. Boggs et al. 1987 and Boggs et al. 1988 and the references therein). In Kelly (2007) and Hogg et al. (2010), problems related to outliers, truncation, censoring and intrinsic scatter are solved with a fully Bayesian approach, with some limitations. Namely, Kelly (2007) requires specifying a single dependent variable, and it is for this variable only that censoring is supported, while Hogg et al. (2010) only consider the two-dimensional case. Finally, Robotham & Obreschkow (2015) extends the results in Hogg et al. (2010) to $n$ dimensions, for arbitrary $n$, and presents analytic likelihoods for $n-1$-dimensional hyperplanes with intrinsic scatter (i.e. lines in a plane, planes in three-dimensional space and so on), assuming data with Gaussian errors.

However, there still exist some remaining issues related to fitting non-linear relations to data that have been discussed to a much lesser extent. One of these is the notion of pre-existing geometry in the measured quantities, such as for angular quantities, or measurements of points on curved surfaces (however, see e.g. Pennec 2006 and Calin & Udriste 2014). Another is the question of a proper characterisation of intrinsic scatter for non-linear relations, and how to incorporate intrinsic scatter when the $n$-dimensional data is not well described by an $n-1$-dimensional subspace (i.e. of codimension one), but an $n-k$-dimensional subspace (codimension $k$), for an arbitrary $k>1$. A suitable resolution of these issues is of interest, since it would enable powerful hypothesis testing via finding the most likely value of $k$ and distribution of intrinsic scatter for each proposed linear or non-linear relation simultaneously.

In this paper, we propose a solution to these issues by formulating the concept of non-linear relations with intrinsic scatter of general codimension $k$ through Riemannian geometry. The novelty of our approach is in extending the idea of intrinsic scatter to curved spaces and correlations that are non-linear and have codimension greater than one, that is, not restricted to linear hyperplanes of dimension $n-1$. The approach also accommodates arbitrary measurement errors, censoring and truncation. Furthermore, we give analytic results for the likelihood and posterior probability for linear $n-k$-dimensional relations in Euclidean spaces. These are easy to implement in fitting codes, and useful for hypothesis testing by enabling quick determination of the most likely codimension of a potential correlation in $n$-dimensional data, along with the parameters of the intrinsic scatter.

In Section 2, we give definitions of correlation relations and intrinsic scatter distributions through the use of Riemannian submanifolds and an intrinsic coordinate system. The definitions allow the concept of intrinsic scatter distributions to be smoothly extended to non-linear relations and correlations of codimension greater than one. We also give a Bayesian solution to the problem of finding the most likely parameters for a given relation and intrinsic scatter, incorporating censoring, truncation and independent, general distributions of measurement errors. In Section 3 we apply our formulation to the special case of linear relations in Euclidean spaces. In particular, we extend the results in Robotham & Obreschkow (2015) to $n-k$-dimensional subspaces with intrinsic scatter,

which correspondingly extends their usefulness to a much wider range of applications. We also extensively discuss the implications of the formalism in the case of linear regression in two dimensions, with comparisons to well-established existing methods. Finally, in Section 4, we apply our methodology to the well-known $M_{\text{BH}}$–$\sigma$ relation between the mass $M_{\text{BH}}$ of a supermassive black hole in a galactic bulge and the bulge velocity dispersion $\sigma$. Finally, we present a summary of the paper in Section 5.

We note that Section 2 is by nature somewhat technical, and relies heavily on concepts in Riemannian geometry. However, the salient points of the formulation we propose are illustrated in Figures 1 through 3. For the reader with immediate applications in mind, we suggest looking at the likelihood functions derived for the special cases of an $n-1$-dimensional linear relation, equation (19), $n-k$-dimensional linear relation, equation (28) and the results for the two-dimensional case, equations (33) and (81)–(84).

2 Relations and intrinsic scatter

2.1 Relations as submanifolds

The dual aim of this paper is firstly to give a geometric formulation of relations with intrinsic scatter that is intuitive and adapted to the geometry of the relation and the space of measurements. Secondly, to fit these relations to data, we aim to provide a method that is fully Bayesian and places no observable in a privileged position, i.e. there is no division into independent and dependent variables. In addition, we wish to accommodate truncated and censored datasets and independent and general distributions of measurement errors. To this end, we need to carefully define what we mean by a relation and its intrinsic scatter. We also need equal care in defining what we mean by being off the relation to be able to define intrinsic scatter in the first place.

In this paper, we will consider relations as submanifolds of the manifold of all possible combinations of measurement values, which we take to have a Riemannian geometry. Being off the relation will be related to the geodesic (shortest path) distance between a point in the measurement manifold, representing a single possible combination of measurements, and the relation submanifold. In this approach, the joint distribution of the measured observables is the primary object of interest. It is this joint distribution that we wish to model with the relation submanifold and the intrinsic scatter. Consequently, the distributions of the individual observables, obtained as marginal distributions, are of secondary interest only.

We start by defining a relation as a smooth map $f:M\rightarrow\mathbb{R}^{k}$, where $M$ is an $n$-dimensional manifold with a Riemannian metric $g_{ab}$ and $k\leq n$. The relation $f$ is fully specified by the $k$ component functions $f_{i}$, so that $f(m)=(f_{1}(m),\ldots,f_{k}(m))$, for $m\in M$. Furthermore, the relation $f$ may depend on $n_{p}$ parameters $\bm{\bm{\theta}}\in\mathbb{R}^{n_{p}}$. In the following, we assume that we are always working with some local coordinate chart $\phi_{x}:M\rightarrow\mathbb{R}^{n}$ and at times use a shorthand $\bm{x}$ for the point $m\in M$ for which $\phi_{x}(m)=(x^{1},\ldots,x^{n})=\bm{x}$. These coordinates represent the observable quantities (e.g. angle, charge, distance, velocity, mass and so on) in the chosen units, or as scale-free logarithmic measurements. The manifold $M$ then represents all possible combinations of measured values for these observables.

Given $M$ and $f$, the level set

[TABLE]

defines an $n-k$-dimensional (or of codimension $k$) subset $S\subseteq M$.111A general level set would be $f(m)=\bm{c}$, $\bm{c}\in\mathbb{R}^{k}$, but we assume that the constant $\bm{c}$ has been subsumed in the definition of $f$.

We further require that $S$ is a regular level set, that is, there are no $m\in M$ for which $f(m)=0$ and the pushforward $\mathrm{d}f_{m}$ fails to be surjective, so that $S$ is additionally an embedded (or regular) submanifold of $M$, as is each $S_{i}$ (see Lee, 2013, for a complete discussion). $S$ then represents the locus of our relation.

The most geometrically obvious way to define being off the relation is to relate it to the geodesic distance of a point from the submanifold $S$. To this end, it is useful to construct a new coordinate system defined with the help of the normal bundle

[TABLE]

of $S$, essentially consisting of all the tangent vectors of $M$ that are orthogonal to $S$ at each point $p$ on $S$. We now define a coordinate system $\phi_{f}$ on $M$ by defining $\phi_{f}(m)=(\bm{\tau}(p),\bm{\nu}(z))\in\mathbb{R}^{n-k}\times\mathbb{R}^{k}$ where $m=\exp_{p}(z)$, $(p,z)\in NS$, and $\bm{\tau}$ is an arbitrary coordinate system on $S$ and $\bm{\nu}$ is an orthonormal coordinate system on $N_{p}S$. Here $\exp_{p}:TM\rightarrow M$ is the exponential mapping around point $p$. In these coordinates, the geodesic distance $\delta$ between $m$ and the submanifold $S$ is just $\delta=\sqrt{g(z,z)}$. Intuitively, the coordinates $\phi_{f}$ correspond to covering the original manifold $M$ using normal coordinates around each point $p\in S$, or $M\subseteq\cup_{p\in S}\{p\}\times EN_{p}$, where we have used a shorthand $EN_{p}=\exp_{p}(N_{p}S)$. See the explanatory Figure 1. The spaces $EN_{p}$ contain all the points that can be reached from $p$ by geodesics orthogonal to $S$ at $p$. However, the subspaces $EN_{p}$ do not necessarily fit together neatly. In general, $\phi_{f}^{-1}$ fails to be injective. This happens when two or more geodesics $\gamma_{i}(t)=\exp_{p_{i}}(tz_{i})$, $i\in I\subset\mathbb{Z}_{+}$, intersect each other at some $m\in M$. Usually normal coordinates are only defined up to these points, which define the boundary of the region where $\exp_{p}$ is a diffeomorphism. In this case, however, it is advantageous to let $EN_{p}$ contain the points reachable from $p$ by geodesics of unlimited length. In this case, $\phi_{f}$ may be multivalued, such that $\phi_{f}(m)=\{(\bm{\tau}(p_{i}),\bm{\nu}(v_{i})\}$, $i\in I$. This is discussed in the Sections 2.2 and 2.3. When $M$ is complete and connected and $S$ is closed, it is known that $\phi_{f}^{-1}$ is surjective (Wolf & Zierau, 1996). In some pathological cases $\phi_{f}^{-1}$ may fail to be surjective, such as when the relation submanifold $S$ has sharp corners, but in this paper we consider only cases where $S$ is sufficiently smooth.

2.1.1 A concrete example

A prototypical astrophysical case is that of a relation between two real-valued observables. In this case, the manifold $M$ is $\mathbb{R}^{2}$, the local coordinates are the identity map $\phi_{x}(x,y)=(x,y)$ and the relation is $f(x,y;\bm{\bm{\theta}})$, parameterized by $n_{p}$ parameters $\bm{\bm{\theta}}\in\mathbb{R}^{n_{p}}$. The normal spaces are $N_{p}S=\{\lambda\cdot(\nabla f)(p)|\lambda\in\mathbb{R}\}$ and we can identify $\exp_{p}(N_{p}S)$ with $N_{p}S$ itself. If we can put our relation $f$ in the form $y-g(x)=0$, we can find the coordinate transformation $\phi_{f}\circ\phi_{x}^{-1}$ as follows. If $m\in M$ is the point under consideration, and $p\in S$ is the point on $S$ geodesically closest to $m$, with $\phi_{x}(m)=(m_{x},m_{y})$ and $\phi_{x}(p)=(p_{x},g(p_{x}))$, then we can take for example $\tau(p_{x},g(p_{x}))=p_{x}$. With this definition, we can in principle solve $\nu$ and $\tau$ as a function of $m_{x}$ and $m_{y}$ from

[TABLE]

See Figure 2. Unfortunately, analytic solutions for these equations are not straightforward to derive except in the case where $g(x)$ is linear.

2.2 Intrinsic distributions

We now have a coordinate system $\phi_{f}$ with which we can define intrinsic scatter in a satisfying way. In fact, we can go beyond simple intrinsic scatter, and define a general unnormalized probability distribution $p_{\text{int}}(\bm{\tau},\bm{\nu})$ on $M$, adapted to the relation $f$ through the coordinates $\phi_{f}$. This amounts to setting up a $k$-dimensional probability distribution on $N_{p}S$ for each point $p\in S$ and using the exponential mapping to create an $n$-dimensional distribution on $M$. However, as mentioned above, $\phi_{f}$ is not injective when the geodesics normal to $S$ intersect each other. In these cases, all the contributing points should be included, with corresponding probabilities added, so that for a point $m\in M$ given in the original coordinates, $\phi_{x}(m)=\bm{x}$, we have

[TABLE]

where $I\subset\mathbb{Z}_{+}$ indexes all the points on $S$ corresponding to $m$. This definition is necessary for cases such as the wrapped normal distribution on a circle $S^{1}$.

The distribution $p_{\text{int}}$ can represent e.g. a physical process producing objects with correlated values of some physical parameters, but with intrinsic scatter given by some probability distribution. An example would be the process producing the correlation between the central black hole mass and velocity dispersion in galaxies (see Section 4). The distributions of the possible physical measurements, as observed in the chart $\phi_{x}$, are implicitly defined as the marginal distributions of the distribution $p_{\text{int}}$. These can be improper, as can the distribution $p_{\text{int}}$ itself. That is, we may have $\int_{M}p_{\text{int}}\,\text{vol}_{\phi_{x}}=\infty$, where $\text{vol}_{\phi_{x}}=\sqrt{\left|g\right|}\mathrm{d}x^{1}\wedge\cdots\wedge\mathrm{d}x^{n}$ is the volume form on $M$ in the chart $\phi_{x}$, and $g=\det(g_{ab})$. A concrete example of an improper $p_{\text{int}}$ is given below.

For applications, a useful distribution is the $k$-dimensional normal distribution, independent of $\bm{\tau}$, or the position on $S$, defined by

[TABLE]

where $\bm{\Sigma}$ is the covariance matrix. Of particular importance is the case where $k=1$, in which case we have a one-dimensional normal distribution depending only on $\nu$, for which

[TABLE]

where the parameter $\sigma$, standard deviation, is now an accurate representation of what is typically called ‘intrinsic scatter’ in astrophysical literature.

It should be noted that the intrinsic scatter distributions (5) and (6) can yield a probability distribution that is not proper. An example is a line $y=\alpha+\beta x$ in flat two-dimensional space, for which, in the coordinate chart, we have (for more details, see Section 3.3)

[TABLE]

for which $\int_{\mathbb{R}^{2}}p_{\text{int}}(x,y)\,\mathrm{d}x\mathrm{d}y$ diverges. Similarly, in this case the marginal distributions for $x$ and $y$,

[TABLE]

are similarly improper.

If necessary, improper distributions $p_{\text{int}}$ arising this way can be made proper e.g. with a suitable truncation in the coordinates $\bm{\tau}$ followed by normalization.

Figure 3 depicts samples of the distribution (6) for three non-linear relations. From the figure it is easy to appreciate the fact that intrinsic scatter can lead to distinct patterns of amplification and attenuation in the probability density near the regions where the underlying relation is strongly curved. As such, a misleading result would be obtained from a fit to such data, if the intrinsic scatter was not modelled properly, or at all.

Finally, we note that intrinsic scatter could also have been defined in two other ways, either along coordinate directions or the scalar geodesic distance to $S$. This seems to exhaust the possibilities, since in general the only directions we have available are the coordinate directions, or alternatively directions along the relation and orthogonal to the direction. In addition, only geodesic distance is a sensible natural choice in a Riemannian context. It should be emphasized that the definitions of intrinsic scatter using the intrinsic coordinates $\phi_{f}$ or the coordinate chart $\phi_{x}$ are in principle completely equivalent, and merely formulated in different coordinates. Both are more general than a definition based only on geodesic distance from the relation $S$.

A definition based on coordinate directions is often used at least in astronomical literature, where in the case of a line in two dimensions $(x,y)$, the intrinsic scatter is typically taken to lie in the direction of the $y$-coordinate (e.g. Press et al., 1992; Tremaine et al., 2002; Gültekin et al., 2009). This gives a relation that is ‘puffed up’ in the $y$-direction. This approach has the obvious problem of not being invariant with respect to a change of coordinates. However, it may be reasonable in some cases, where we are certain that whatever processes are producing the intrinsic scatter only operates strictly along a particular measured observable. Nevertheless, considering the often very correlated nature of astronomical observables, establishing this fact seems to be a difficult proposition. Consider for example the $M_{\text{BH}}$–$\sigma$relation between the mass of a supermassive black hole and the central velocity dispersion of its host galaxy, where the intrinsic scatter is typically taken to lie entirely along the $M_{\text{BH}}$ axis. As galaxies together with their SMBHs evolve, they move on the $M_{\text{BH}}$–$\sigma$plane, but with different galaxies starting at (nearly) the same point, not necessarily ending up near each other later due to differences in composition, environment, merger history and so on. This produces intrinsic scatter. However, it seems very difficult to justify that these processes should operate only along the black hole mass axis. Rather, it seems that they might operate in any which direction indeterminately, but of these, the only one evident to us is the direction away from the underlying relation. Indeed, if the processes producing intrinsic scatter would only operate along the relation, the resulting data would have no scatter around the relation at all. Finally, comparing the magnitude of intrinsic scatter between datasets becomes difficult, since the value will end up depending on the orientation of the line (see Section 3.3.1). In this sense, it seems reasonable to formulate intrinsic scatter problems in terms of the intrinsic coordinates $\phi_{f}$, separately considering the directions orthogonal and parallel to the proposed relation.

As the last option, the intrinsic scatter could also have been defined as a probability density in the original measurement coordinate system $\phi_{x}$, depending only on the scalar geodesic distance to $S$. This approach would give relations that are ‘puffed up’ symmetrically in all directions, with density falling off as a given function of the geodesic distance. Probability densities using this approach would not present the enhancements shown in Figure 3. However, this approach is strictly less general than the one outlined above, and does not allow for e.g. wrapped distributions on spheres, toroids and so on, and cannot cope with directionally varying scatter for submanifolds $S$ with codimension higher than one.

2.3 Fitting a relation with intrinsic scatter

In a typical case, we have obtained $n_{d}$ measurements $\mathcal{D}=\{(m_{1},h_{1}),\ldots,(m_{n_{d}},h_{n_{d}})\}$, where $m_{i}\in M$ are the measured values and $h_{i}:M\rightarrow\mathbb{R}$ are the probability distributions $p(\eta_{i}|m_{i})$ for the true value $\eta_{i}$, given the measurement $m_{i}$, typically representing the estimated measurement errors. The $h_{i}$ are otherwise unspecified, and arbitrary distributions of measurement errors, including upper or lower limits (left or right censoring, respectively) can be naturally accommodated. For optimal results, the distributions $h_{i}$ should be Bayesian posterior probability densities obtained in the measurement process, fully representing the available information. However, in some cases $h_{i}$ are completely unknown, in which case they must be estimated, typically as Gaussians with unknown covariance matrices, which then must taken as additional nuisance parameters, in addition to the true values $\eta$. In addition to the measurements, we have a postulated relation $f$, parameterized by $n_{p}$ parameters $\bm{\theta}\in\mathbb{R}^{n_{p}}$, from which we derive coordinates $\phi_{f}(m;\bm{\theta})$, also parameterized by $\bm{\theta}$. We believe the data to fulfil the relation $f$, up to some intrinsic scatter $p_{\text{int}}$, parameterized by $n_{s}$ parameters $\bm{\varphi}\in\mathbb{R}^{n_{s}}$. We would now like to find the most probable values of $\bm{\theta}$ and $\bm{\varphi}$, constrained by the data $\mathcal{D}$.

The Bayesian way to proceed is to note that for a single datum $(m,h)$, the conditional probability of the measurement factors as $p(m,h|\bm{\theta},\bm{\varphi})=\int p(m,h|\eta)p_{\text{int}}(\eta|\bm{\theta},\bm{\varphi})\,\mathrm{d}\eta$, where $\eta$ is the true value, drawn from the intrinsic distribution, $p_{\text{int}}$. As such, the true value $\eta$ is taken as a nuisance parameter, and integrated out. In addition, it may be that the measurement process is not fully sensitive across the range of possible values, leading to truncation, or data points that are not seen in the sample. If the probability that a given value $\eta$ leads to a detection is given by $p_{\text{det}}(\eta)$, we have for the measurement probability that $p(m,h|\bm{\theta},\bm{\varphi})=\int p(m,h|\eta)p_{\text{det}}(\eta)p_{\text{int}}(\eta|\bm{\theta},\bm{\varphi})\,\mathrm{d}\eta$. See Appendix A for full derivation.

Now, taken as a function of the parameters, the probability $p(\mathcal{D}|\bm{\theta},\bm{\varphi})$ defines the likelihood $\mathcal{L}$ of $\bm{\theta}$ and $\bm{\varphi}$. Assuming that the data $\mathcal{D}$ are independent, we get

[TABLE]

where the integration can be done either in the coordinate chart $\phi_{x}$ or the intrinsic chart $\phi_{f}$. The integration in $\phi_{f}$ automatically takes care of points $m\in M$ that have multiple representations $(\bm{\tau}_{j},\bm{\nu}_{j})$, but if the integration is done over $\phi_{x}$, these have to be explicitly summed over.

If $\bm{\theta}$ and $\bm{\varphi}$ have a joint prior distribution $\pi_{\bm{\theta},\bm{\varphi}}$ (possibly uninformative), we can define the joint posterior probability distribution of $\bm{\theta}$ and $\bm{\varphi}$,

[TABLE]

The posterior distribution, equation (11), contains all the information of how the data constrains the parameters $\bm{\theta}$ of the relation itself, as well as the parameters $\bm{\varphi}$ of the intrinsic scatter model, given our prior information $\pi_{\bm{\theta},\bm{\varphi}}$. In some special cases, the posterior distribution can be computed analytically, but in general it is necessary to employ numerical techniques, such as a suitable Markov Chain Monte Carlo (MCMC) method

(see e.g. Gelman et al., 2013). Finally, point estimates such as the Maximum (Log-)Likelihood Estimate (MLE) and Maximum A Posteriori Estimate (MAP) can be obtained using the likelihood via standard procedures. Full derivation of the likelihood (10) and the posterior distribution (11) following the style of Jaynes (2003) is found in Appendix A.

2.4 Some considerations

For the formulation presented above to make sense, a reasonable geometry must exist between the measured quantities. This is not a problem e.g. for measurements of a 3D angle, where we have the natural geometry of the sphere $S^{2}$. In other cases, such as for joint measurements of galaxy luminosity and velocity dispersion it is reasonable to question whether a distance measured ‘across’ the measured quantities makes any sense. This observation has been discussed numerous times since Pearson in 1901 originally presented a least squares method using orthogonal distances to the regression line (or plane), typically called Orthogonal Regression (OR) (Pearson, 1901). See for example Isobe et al. (1990) for a discussion.

The Riemannian approach presented in this paper makes it possible to specify the units and any possible geometry between the measured quantities beforehand. After the geometry is set, the metric character of the approach ensures that results are invariant with respect to all coordinate transformations, which in this case would typically represent scaling of the measured quantities (i.e. a change of units). Another natural possibility is to introduce the measurements in an entirely scale-free manner by using logarithms.222 Linear measurements $x^{i}$ with a metric $g=\operatorname*{diag}(1,\ldots,1)$ are equivalent to logarithmic measurements $\log_{10}x^{i}$ and a metric with an exponential dependence on the measured values, $g=\log(10)\operatorname*{diag}(10^{x^{1}},\ldots,10^{x^{n}})$.

A practical problem is the fact that the coordinate transformation between $\phi_{x}$ and $\phi_{f}$ may prove difficult or impossible to derive analytically. Finding the transformation and computing the posterior probabilities numerically necessitates computing numerous geodesic distances in a curved space. In general this is highly computationally demanding, although efficient algorithms exist for some special cases (Crane et al., 2013).

Finally, we note that there are to be some consequences from and natural restrictions to the choice of parameterization and the prior distributions of the parameters. Firstly, it is clear that an arbitrary transformation of parameters $\bm{\theta}=g(\bm{\theta}^{\prime})$ will result in the change of the prior probability measure so that $\pi(\bm{\theta})\mathrm{d}\bm{\theta}=\pi(g(\bm{\theta}^{\prime}))\det(J_{g})\mathrm{d}\bm{\theta}^{\prime}$, where $J_{g}$ is the Jacobian of the transformation. This can be understood also through the fact that a probability distribution $\pi$ defined in the parameter space also yields through the relation $f$ and $p_{\text{int}}$ a probability distribution on $M$. This latter distribution, once defined, should naturally be invariant under a change of parameterization, which is guaranteed by the transformation formula above.

Secondly, prior distributions are usually considered to be somewhat arbitrary and application-specific, often chosen for convenience to yield analytic results for the posterior probability. However, as noted above, this freedom can be restricted by what we know of how the data is distributed on $M$ combined with the choice of a relation $f$ and $p_{\text{int}}$. An important special case is if we have no prior knowledge of how the data is distributed, but only the possible values that it may take, through the choice of $M$. In this case the maximum entropy distribution on $M$ is uniform. This distribution places restrictions on the choice of $f(\bm{x};\bm{\theta})$, and after $f$ is chosen, forces a particular choice of $\pi(\bm{\theta})$, which can then in this context be called the uninformative prior. The general argument is outside the scope of this paper, but see e.g. George & McCulloch (1993), Broemeling & Broemeling (2003) and Fraser et al. (2010) for historical and modern discussions on vague and uninformative priors from different points of view. However, a special case of an uninformative prior defined in the way described above is considered in detail in Section 3.3 and later used in Sections 3.3.1 and 4.

3 Worked out examples

The following sections contain analytical results for particular choices of the measurement manifold $M$ and the relation $f$. In each case, we show how the framework presented above is used to obtain the likelihood, which can then be used for parameter estimation. The outline of the process, or a ‘practical how-to’, is as follows:

Determine the manifold $M$, along with the metric $g$. These are typically fixed by the nature of the measured quantities. 2. 2.

Fix the functional form $f(x;\bm{\theta})$ of the relation to be fit, along with the parameters $\bm{\theta}$ of the relation. 3. 3.

Define the intrinsic coordinates $(\tau,\nu)$. A well-motivated choice for the coordinates $\tau$ along the relation is using geodesic distances from a given fixed origin. The coordinates $\nu$ are then defined up to a rotation of the $\nu$-coordinate axes. 4. 4.

Determine the coordinate transformation $\phi_{x}\circ\phi_{f}^{-1}$ from $(\tau,\nu)$ to $x$ or the inverse transformation $\phi_{f}\circ\phi_{x}^{-1}$. This may only be possible analytically in one direction only, or may not be possible analytically at all. 5. 5.

Choose a suitable intrinsic scatter model $p_{\text{int}}(\nu,\tau)$ as a function of the intrinsic coordinates. 6. 6.

Evaluate the likelihood integral, equation (10). For analytical results, this typically needs to be done only for a single arbitrary datapoint, if all the datapoints have a similar error distribution model. 7. 7.

Compute the likelihood over all data as a product of the likelihood of the individual likelihoods for each data point. Equivalently, choose a suitable prior and compute the posterior distribution, equation (11), over all the data.

We will refer to these steps in the following derivations so that the process can be more easily followed.

3.1 Linear $n-1$-dimensional case

A highly useful special case is the case of normally distributed intrinsic scatter, equation (6), in a linear relation $f:\mathbb{R}^{n}\rightarrow\mathbb{R}$ between $n$ observables, defining an $n-1$-dimensional (codimension $1$) hyperplane $S$, with measurements drawn from an Euclidean space $M=\mathbb{R}^{n}$. For now, we assume that the metric $g$ is given by the $n\times n$ identity matrix I. These definitions complete step 1 from above. This special case represents well the typical astrophysical problem of modelling a linear relation between multiple physical variables that have no special geometry, such as mass, luminosity or velocity dispersion.

To address step 2, we fix the parametrization. All the planes $S$ can be parameterized with $\bm{\bm{\theta}}=\bm{p}\in\mathbb{R}^{n}$ so that

[TABLE]

where $\bm{x}\in\mathbb{R}^{n}$. With this parameterization, the vector $\bm{p}$ corresponds to the point of $S$ where it is closest to the origin, and as such $\bm{p}$ is also normal to the plane $S$. A different parameterization, more often used in the astrophysical literature, emphasizes one particular coordinate dimension, written as

[TABLE]

Converting between these two parameterizations is accomplished through

[TABLE]

We further assume that for each measurement $(\bm{m},h)$, where now $\bm{m}\in\mathbb{R}^{n}$, the measurement error can be modelled with an $n$-dimensional normal distribution

[TABLE]

where $\bm{\Sigma}$ is the measurement error covariance matrix. This assumption is necessary to obtain the analytic results below, but it is also well suited to published astrophysical data, for which the complete posterior distributions for each datum are typically not available. Upper limits can be incorporated by substituting one or more degrees of freedom $x_{i}$ in equation (16) with suitable one-dimensional distributions representing the limit distributions.

To complete steps 3 and 4, we need to specify the intrinsic coordinates and the coordinate transformations. One method to construct the coordinates $\phi_{f}$ is to transform the natural coordinate frame of the measurement space $M$ with a translation by $\bm{p}$ and a rotation $\textbf{{R}}\in SO(n)$, which takes the $n$’th basis vector $\hat{\bm{e}}_{n}$ to the unit vector $\hat{\bm{p}}$ (hereafter, we use a hat to signify a unit vector). The new components of a vector $\bm{x}$ are then $\phi_{f}(\bm{x})=\textbf{{R}}^{T}(\bm{x}-\bm{p})=(\bm{\tau},\nu)\in\mathbb{R}^{n-1}\times\mathbb{R}$. Defining $\bm{w}=(p_{0},\ldots,p_{n-1},0)$, $\cos\theta=\hat{\bm{p}}^{T}\hat{\bm{e}}_{n}$ and $\sin\theta=\hat{\bm{p}}^{T}\hat{\bm{w}}$, we can compute R with

[TABLE]

where

[TABLE]

is the usual $2$-dimensional rotation matrix. The origin of the $\bm{\tau}$ coordinates is then at $\bm{p}$ (in chart $\phi_{x}$), with the orientation of the coordinate axes $\hat{\bm{\tau}}_{i}$ given by $\hat{\bm{\tau}}_{i}=\textbf{{R}}\hat{\bm{e}_{i}}$, for $i=1,\ldots,n-1$. The orthogonal intrinsic coordinate axis is $\hat{\bm{\nu}}_{1}=\textbf{{R}}\hat{\bm{e}}_{n}$.

Step 5 requires specification of the intrinsic scatter model $p_{\text{int}}$. We assume the normal intrinsic distribution, equation (6). With this definition, combined with the assumed measurement error distribution, and assuming no censoring, that is no upper or lower limits, we can complete step 6. We evaluate the equation (10) for the likelihood of a single measurement datapoint. The result is

[TABLE]

where

[TABLE]

We can also write $\tilde{m}_{n}=\hat{\bm{p}}^{T}\bm{m}-\left\|\bm{p}\right\|$ and $\tilde{\Sigma}_{nn}=\hat{\bm{p}}^{T}\bm{\Sigma}\hat{\bm{p}}$, which transforms equation (19) to the form used in Robotham & Obreschkow (2015).

We can also address the case where the measurements are given in coordinates where the metric is represented by a constant matrix $\textbf{{G}}\neq\textbf{{I}}$. This can be the case when the metric G encodes the choice of units for the variables, or when it has off-diagonal terms. A situation that would lead to a non-diagonal metric is, for example, a case where we are interested in two observables $A_{1}$ and $A_{2}$, but can only measure $B_{1}=A_{1}$ and $B_{2}=A_{1}+A_{2}$. The observables may be e.g. inflows of current or liquid from two independent sources, but we can only measure one source directly, and the sum of the flows somewhere downstream. The flows $A_{i}$ can by assumption have any values, so the metric in coordinates $A_{i}$ should be a diagonal product metric, whereas in coordinates $B_{i}$ off-diagonal terms appear. For the case of a non-identity constant metric $G$, the equation (19) is still valid, but we have

[TABLE]

where W is a diagonal matrix of the square roots of the eigenvalues of G, and $\textbf{{P}}\in O(n)$, so that $\textbf{{W}}^{-1}\textbf{{P}}^{T}\textbf{{G}}\textbf{{P}}\textbf{{W}}^{-1}=\textbf{{I}}$. This is possible, since the metric is assumed to be Riemannian, in which case all the eigenvalues of G must be positive. Note that if the metric is not constant, the shortcut presented above does not work, and all the steps 1-7 have to be followed using the metric explicitly.

3.2 Linear $n-k$-dimensional case

The approach in the previous section can be readily generalized to the case of $k>1$ simultaneous linear relations, defining an $n-k$-dimensional (codimension $k$) affine subspace $S_{n-k}$. In this case, it is easiest to define the relation through defining the intrinsic coordinates first. As such, we go through steps 3 and 4 first. This can be accomplished by starting with the $n-1$-dimensional subspace $S_{n-1}$ defined by the first relation, through the $n$ parameters $\bm{p}_{n}\in\mathbb{R}^{n}$, as above. Now, the subspace $S_{n-k}$ must lie in the intersection of $S_{n-1}$ and an $n-2$-dimensional subspace $S_{n-2}$, which itself must also lie entirely in $S_{n-1}$. We may parameterize $S_{n-2}$ within $S_{n-1}$ with the $n-1$ parameters $\bm{p}_{n-1}=(p_{n-1,1},\ldots,p_{n-1,n-1},0)\in\mathbb{R}^{n}$. This process is then continued until we have given $p_{n-k+1}$, parameterizing $S_{n-k}$ in $S_{n-k+1}$. The number of parameters is $(2nk-k^{2}+k)/2$ in total. Each parameter vector $\bm{p}_{i}$ yields a rotation matrix $\textbf{{R}}_{i}$ as in the section above, which together with $\bm{p}_{i}$ defines the coordinate transformation from the intrinsic coordinates $\phi_{f}$ on $S_{i}$ to the intrinsic coordinates on $S_{i-1}$, with the understanding that $S_{n}=M$. Following this process through, we find that the coordinate transformation from $\phi_{x}$ to the intrinsic coordinates on $S_{n-k}$ is given by

[TABLE]

where

[TABLE]

and

[TABLE]

Again, the origin of the intrinsic coordinate axes is at $\bm{p}$. The orientation is similarily given by $\hat{\bm{\tau}}_{i}=\textbf{{R}}\hat{\bm{e}}_{i}$ for $i=1,\ldots,n-k$ and $\hat{\bm{\nu}}_{i}=\textbf{{R}}\hat{\bm{e}}_{n-k+i}$ for $i=1,\ldots,k$. This process defines the relation $f$ in a roundabout way as the equation $\bm{\nu}=\bm{0}$, which in the chart $\phi_{x}$ is then given by the system of $k$ equations

[TABLE]

where $\bm{\theta}=(\bm{p}_{n},\ldots,\bm{p}_{n-k+1})$. This completes step 2.

We again assume that measurement errors are normally distributed, given by equation (16). For step 5 we now assume that the intrinsic scatter is normally distributed in $\bm{\nu}$ with a covariance matrix $\bm{\Sigma}_{\text{int}}$, as in equation (5). With this definition, we can complete step 6 and compute the likelihood given by a single measurement $\bm{m}$. The result is

[TABLE]

where

[TABLE]

A possible constant non-Euclidean metric can be accommodated as in the previous section.

The likelihood, equation (28) is readily generalized for intrinsic distributions other than the multivariate normal distribution through mixture models (see e.g. the approach in Kelly, 2007). In addition, the result in this section is useful for hypothesis testing in the sense of finding the most likely value of $k$ for a given $n$-dimensional dataset.

3.3 The line in two dimensions

It is useful to work out the two-dimensional special case of a line with intrinsic scatter in detail, considering the amount of literature focusing on symmetric fitting of linear relations with and without intrinsic scatter (e.g. Pearson, 1901; Boggs et al., 1987; Isobe et al., 1990; Feigelson & Babu, 1992; Robotham & Obreschkow, 2015, and many others). As such, step 1 consists of setting $M=\mathbb{R}^{2}$ with an Euclidean metric. In step 2, for ease of comparison with existing methods, instead of the parameterization used above, we define the relation $f$ and the line $S$ with

[TABLE]

where now $\bm{\bm{\theta}}=(\alpha,\beta)\in\mathbb{R}^{2}$. The measurement error in the case of no censoring can now be represented with a two-dimensional normal distribution, given by

[TABLE]

where $x_{0}$ and $y_{0}$ specify the measured values, $\sigma_{x}$, $\sigma_{y}$ represent the measurement uncertainties in the $x$ and $y$ directions, and $\rho\in[-1,1]$ specifies the correlation between the measurement errors. Upper and lower limits can be introduced as in Section 3.1. While the parameterization through $\alpha$ and $\beta$ is convenient and intuitive, it is not manifestly symmetric with respect to the coordinates. This will be investigated further below.

Keeping the assumption that the internal scatter is normally distributed, via equation (6), we can use the results for steps 3 to 6 from Section 3.1. The likelihood of a single measurement in the case of no censoring is then obtained from equation (19), yielding

[TABLE]

where now

[TABLE]

so that $\nu^{2}$ is the squared orthogonal distance from the relation, and $\tilde{\sigma}$ is an extended uncertainty incorporating both intrinsic scatter and the measurement errors. Analytic results for particular censored error distributions can also be found. See Appendix B.

We note that equation (35) has the correct asymptotic behaviour with respect to $\beta$, in the sense that when the regression line tends towards the vertical, or $\beta\rightarrow\infty$, we have $\tilde{\sigma}\rightarrow\sigma_{x}^{2}+\sigma^{2}$, agreeing with intuitive result that all of the uncertainty should in this case be a combination of the horizontal (along $x$-axis) and intrinsic scatter. Likewise, when the regression line tends towards the horizontal, or $\beta\rightarrow 0$, we have $\tilde{\sigma}\rightarrow\sigma_{y}^{2}+\sigma^{2}$, similarly agreeing with geometric intuition.

We now apply the discussion in Section 2.4 and derive a strictly unique uninformative prior probability density for the parameters $\bm{\theta}=(\alpha,\beta)$. This is appropriate for the case where no data has yet been collected, and we have no information on where on the plane the relation is and in what orientation. Thus, we need to find a prior density for $\bm{\theta}$ which yields a uniform distribution on $M$ when integrated over. That is, the lines must cover the plane evenly. A direct approach seems to necessitate set-based analysis, but in this special case we can apply existing results from the literature, derived via other means.

For example, demanding invariance under Euclidean coordinate transformations (simultaneous rotation and translation) yields an invariant prior probability measure $\pi(\alpha,\theta)\mathrm{d}\alpha\mathrm{d}\theta=\mathrm{d}\alpha\mathrm{d}\sin\theta$, where $\theta=\arctan\beta$ is the angle between the line and the $x$-axis. This is equivalent to $\pi(\alpha,\beta)=(1+\beta^{2})^{-3/2}$. The same result can be obtained by demanding that the prior density be invariant under a switch of the coordinates, that is under $x\leftrightarrow y$, $\sigma_{x}\leftrightarrow\sigma_{y}$, $\beta\leftrightarrow 1/\beta$ and $\alpha\leftrightarrow-\alpha/\beta$. This result was apparently originally derived by E.T. Jaynes in 1976 (reprinted in Jaynes 1983).

Finally, for $\sigma$ there is considerably more leeway in the literature with regards to the choice of an uninformative prior. A typical choice is the Jeffreys prior, $\pi(\sigma)\propto 1/\sigma$. The complete uninformative prior for this special case is then

[TABLE]

Equation (36) and (33) can now be combined to yield the unnormalized posterior distribution

[TABLE]

This gives a Bayesian solution to the problem of symmetric fitting of a linear relation with intrinsic scatter to two-dimensional data with heteroscedastic errors in both measured variables, including possible upper or lower limits and truncation. This result thus extends the earlier Bayesian results in Zellner (1971), Gull (1989), Jaynes (1991, unpublished) and Kelly (2007). It should be noted that while the likelihood, equation (33) is invariant under the switch of coordinates, the posterior distribution in itself is not, if the Jacobian of the transformation of parameters is not included, as discussed in Section 2.4. In practice this means that in the usual case of $f=y-\alpha-\beta x$ we wish to find the maximum of

[TABLE]

In the inverse case, where $f=x-\alpha^{\prime}-\beta^{\prime}y$, and $\alpha^{\prime}=-\alpha/\beta$, $\beta^{\prime}=1/\beta$, we should maximize

[TABLE]

where $\beta^{\prime 3}$ is the Jacobian of the transformation $(\alpha,\beta)\mapsto(\alpha^{\prime},\beta^{\prime})$. This works since for the priors we have $p(\alpha^{\prime},\beta^{\prime},\sigma)=\left|\beta^{3}\right|p(\alpha,\beta,\sigma)$ and for the likelihood we have

[TABLE]

as found above.

The result in this section can be easily extended to cases where the intrinsic scatter is not normally distributed or the measurement errors are not distributed with a bivariate normal distribution, by approximating the distributions with a weighted sum of Gaussians, as used e.g. in the linmix_err-method of Kelly (2007). However, the result is not applicable for situations more general than the linear case considered here, such as when the $M$ itself is not trivially Euclidean, but includes e.g. angular or directional measurements or when the relation $f$ is non-linear. In these cases, the likelihood function and consequently the posterior distribution may have to be evaluated numerically.

3.3.1 Comparison to some existing approaches

The equations (33), (34) and (35) are a fundamentally symmetric way to describe the likelihood of parameters $(\alpha,\beta,\sigma)$ specifying a line with orthogonal intrinsic scatter. Several earlier works have incorporated intrinsic scatter as an error parameter, added in quadrature to the measurement errors along some specific measurement axis, as in e.g. the least likelihood method in Gültekin et al. (2009) and the FITEXY method (Press et al., 1992) as modified in Tremaine et al. (2002). These methods introduce intrinsic scatter into the measurement errors of the dependent variable (for now taken to be $y$), yielding a total variance of the form

[TABLE]

in the case of normally distributed intrinsic scatter in the $y$-direction, where $\sigma^{2}_{\text{int},y}$ is the variance of the intrinsic scatter and $\rho=0$ is assumed. This is not equivalent to equation (35), and in particular, the missing normalization factor $1+\beta^{2}$ leads to $\bar{\sigma}\rightarrow\infty$ as $\beta\rightarrow\infty$. We will now discuss how this difference arises.

Firstly, assume that the underlying distribution of the data is a normal distribution orthogonal to a line $y=\alpha+\beta x$, i.e. given by equation (6) with $\nu=(y-\alpha-\beta x)/\sqrt{1+\beta^{2}}$. In this case the conditional distributions of the intrinsic scatter in the $x$- and $y$-directions are also normal, with variances

[TABLE]

However, the marginal distributions of $x$ and $y$ are not normal. Now, the result in equation (41) can be obtained from the likelihood produced by a single measurement,

[TABLE]

where $h$ is the 2-dimensional Gaussian given by equation (32) and $p_{\text{int}}$ is the normal distribution, equation (6). The assumptions yielding equation (41) are to take $\nu=y-\alpha-\beta x$, i.e. distance from the regression line in the $y$-direction, and to set $\sigma=\sigma_{\text{int},y}$. This yields a normal distribution, as in equation (33), but with

[TABLE]

where the subscript $a$ refers to asymmetric, which is a point we will discuss below. If we set $\rho=0$, we have the likelihood used in Gültekin et al. (2009),

[TABLE]

as well as the $\chi^{2}$ value used in Tremaine et al. (2002). Likewise, if we instead use $\nu=(y-\alpha-\beta x)/\sqrt{1+\beta^{2}}$, i.e. orthogonal distance from the regression line, in the equation (44), the result is the set of equations (33)-(35). However, at this point we find that making the obvious substitution from equation (43) does not make the set of equations (33)-(35) equivalent to equation (33) combined with the substitutions (45) and (46). Furthermore the $\chi^{2}$ in equation (47) is invariant under the switch of independent and dependent coordinates, together with $\sigma_{\text{int},y}\leftrightarrow\sigma_{\text{int},x}$. However, the variance $\tilde{\sigma}^{2}_{a}$ and consequently the likelihood are manifestly not invariant. The equation (47) will give smaller likelihoods than equation (33) for $\beta\gg 1$ or $\beta\sim 0$, depending on which variable is taken as the dependent one.

The fundamental reason for this state of affairs is the fact that if the intrinsic scatter is modelled as a proper probability distribution in any particular coordinate direction, the normalization of the intrinsic scatter as a distribution orthogonal to the relation will change in normalization with change in $\beta$. This reflects the breaking of the rotational symmetry of the problem. For example, assume that the intrinsic scatter is modelled as a normal distribution in the $y$-direction, with variance $\sigma^{2}_{\text{int},y}$. In this case the distribution of the data, being the intrinsic scatter in a direction orthogonal to the regression line, will need to have a normalization constant that goes down with increasing slope $\beta$. Indeed, as $\beta\rightarrow\infty$, the normalization constant will need to tend towards zero, to keep the conditional probability distribution in the $y$-direction normalized to unity. This leads to a vanishing likelihood for a vertical line $\beta\rightarrow\infty$, which results from the fact that in this case $\tilde{\sigma}_{a}\rightarrow\infty$ and consequently the resulting likelihood $\mathcal{L}(\alpha,\beta\rightarrow\infty|\mathcal{D})\rightarrow 0$. In mathematical terms, if we have a distribution of the data depending only on the orthogonal distance, $p_{\text{int}}(\nu)$, and we have demanded that in $y$-direction we should have a proper probability density function $p_{\text{int},y}(y)$, that is

[TABLE]

then necessarily

[TABLE]

which goes to zero as $\beta\rightarrow\infty$. This is inconsistent, if we expect the distribution of the data, $p_{\text{int}}$, to have the same normalization no matter which way the regression line might point. Armed with this knowledge, if we now multiply the right side of equation (47) by $\sqrt{1+\beta^{2}}$ and substitute $\sigma^{2}_{\text{int},y}$ from equation (43), we do arrive at the equations (33)-(35). Based on this analysis, the likelihood of equation (33) is recommended for the purpose of fitting a regression line with normally distributed intrinsic scatter to data with normally distributed measurement errors. This point is also raised in Robotham & Obreschkow (2015).

To investigate this conclusion numerically, we generated sets of simulated observations from linear relations with intrinsic scatter. The data were then fit using linmix_err, FITEXY,333We used the implementation available in the MPFIT package (Markwardt, 2009) through the MPFITEXY wrapper routine (Williams et al., 2010). and by maximizing the posterior probability, equation (37). This last method is referred to as Geo-MAP (Geometric Maximum A Posteriori) hereafter. The parameter estimates and 1-$\sigma$ uncertainties for the linmix_err method were taken using the median and standard deviation of the Markov chain outputs. For the FITEXY method the estimate and uncertainty provided by the algorithm were used. For the Geo-MAP method a bootstrapping method was used together with a numerical maximization of the posterior distribution to yield the parameter estimates and 1-$\sigma$ uncertainties. The simulated datasets were generated assuming a relation $y=\alpha+\beta x$, with a normally distributed orthogonal intrinsic scatter $\sigma=0.1$ and normally distributed measurement error with $\sigma_{x}=\sigma_{y}=0.1$ and $\rho=0$. In total six datasets were generated, with $n_{d}=100$ data points each and with $\alpha\in\{0,10\}$ and $\beta\in\{0.1,1,10\}$. The datasets represent a situation where the intrinsic scatter and measurement errors contribute equally to the observed scatter, and as such approximations assuming the relative smallness of either are maximally violated. The data were generated with a uniform distribution $U(-1,1)$ along the relation, and then shifted first according to $\sigma$ and then according to $\sigma_{x}$ and $\sigma_{y}$. The data was then fit both in the forward direction, with $y=\alpha+\beta x$ as well as the inverse direction, with $x=\alpha^{\prime}+\beta^{\prime}y$. For the inverse direction, the slope $\beta=1/\beta^{\prime}$ was computed after the fit. The slopes of the resulting fits are listed in Table 1. Figure 4 illustrates the differences between the methods for the $\alpha=0$, $\beta=10$ case.

The results show that for values of $\beta\gg 1$ the inverse fits with the linmix_err and FITEXY methods agree better with the data, albeit with high uncertainties. Similarly, for $\beta\ll 1$ the forward fits give a better agreement with the data. This is in line with the observations above, showing that these methods give results tending towards lower values of $\beta$ for $\beta\gg 1$ when used in the forward direction, and correspondingly towards higher values of $\beta$ for $\beta\sim 0$, when used in the inverse direction. It should be noted that in addition to yielding estimates of $\beta$ that are much too high or low in these cases, both linmix_err and FITEXY methods give corresponding uncertainty estimates that are small, so that the true value of $\beta$ ends up as a multiple-$\sigma$ outlier. Similarly as expected, the Geo-MAP method gives nearly identical results for forward and inverse fits, up to the noise caused by the bootstrapping procedure. In addition, the true value of $\beta$ is always contained within the 1-$\sigma$ bounds except for the case of $(\alpha,\beta)=(10,0.1)$. Finally, we note that none of the methods appear invariant with respect to a shift with $\alpha$.

We conclude that methods based on likelihood examining distances from the regression line in the direction of a particular coordinate may give estimates of slope that are far from the true value while simultaneously providing tight error bounds. In particular, if the distance is measured along the dependent variable, the resulting method will be biased towards $\beta=0$. For methods, like the maximum-likelihood method in (Gültekin et al., 2009), this effect remains even when the intrinsic scatter $\sigma\rightarrow 0$, in which case the distribution $p_{\text{int}}$ becomes a delta-function ridge, and the problem reduces to Ordinary Least Squares (OLS) regression. This tendency towards lower slopes was reported in Park et al. (2012) for all the methods used in the paper, namely OLS, BCES (Akritas & Bershady, 1996), FITEXY, linmix_err, and the Gültekin et al. (2009) maximum-likelihood method. Our results also agree with Isobe et al. (1990), who studied different regression methods in the limit of vanishing measurement error.

3.4 Ellipses in Euclidean space

Next, we consider a case where the relation $f$ is properly non-linear. With astrophysical applications in mind, we seek analytic results for ellipses embedded in a Euclidean space. A possible application could be, for example, to model the orbit of a stream of low mass particles. As such, for step 1 we consider our measurements to be points in $M=\mathbb{R}^{3}$ with the Euclidean metric. For step 2, we define the ellipse with

[TABLE]

where $\bm{\bm{\theta}}=(\bm{p},a,e,i,\omega,\Omega)$. Here the parameters $\bm{\bm{\theta}}$ constitute $\bm{p}\in\mathbb{R}^{3}$, which specifies the position of the centre of the ellipse, and $a$, $e$, $i$, $\omega$ and $\Omega$, which are the orbital elements, namely semimajor axis, eccentricity, inclination, argument of pericentre and longitude of the ascending node, respectively. The tilde signifies a coordinate transformation to a vector $\tilde{\bm{x}}=\textbf{{R}}^{T}(\bm{x}-\bm{p})$, where the matrix R rotates the Euclidean coordinate triad $\{\hat{\bm{e}}_{1},\hat{\bm{e}}_{2},\hat{\bm{e}}_{3}\}$ into a coordinate triad $\{\tilde{\bm{e}}_{1},\tilde{\bm{e}}_{2},\tilde{\bm{e}}_{3}\}$ aligned with the ellipse. That is, $\tilde{\bm{e}}_{1}$ points along the semimajor and $\tilde{\bm{e}}_{2}$ along the semiminor axis. The matrix R can be written explicitly as the product

[TABLE]

In the ellipse-aligned frame, the ellipse can now be characterized also as the curve

[TABLE]

where $E\in[0,2\pi)$ is typically called the eccentric anomaly. To use this convenient form, we will work in the coordinate frame aligned with ellipse in the following.

We can now define the intrinsic coordinates (step 3) by defining $\tau$ to be the arc length of the ellipse, with $\tau=0$ in the direction of $\tilde{\bm{e}}_{1}$ and increasing towards $\tilde{\bm{e}}_{2}$. We will need the relation between the eccentric anomaly $E$ and the arc length $\tau$, given by

[TABLE]

where

[TABLE]

is the incomplete elliptic integral of the second kind. The orthogonal coordinate axes can be defined so that $\hat{\bm{\nu}}_{1}$ is parallel to $\tilde{\bm{e}}_{3}$ and $\hat{\bm{\nu}}_{2}$ points orthogonally away from the ellipse, and thus can be given as

[TABLE]

with normalization omitted for clarity.

For step 4, we can now construct the coordinate transformation $\phi_{x}\circ\phi_{f}^{-1}$ to find

[TABLE]

However, the inverse transformation is already very difficult to find analytically, and further problems caused by the nonlinearity will now also manifest.

At this point, step 5, we should specify the intrinsic scatter model. To mitigate analytic problems, we would again like to use the simplest normal distribution, so that

[TABLE]

However, we find that even if we consider a similarly simplified measurement error distribution,

[TABLE]

it is very difficult to obtain a closed form result for the likelihood integral, equation (10). Neither is it obvious how to compute the form of $p_{\text{int}}$ in the measurement coordinates $\phi_{x}$, whether in the original or ellipse-aligned form.

To make progress, we have to make some assumptions. The main assumption we need is that the intrinsic scatter $\sigma$ is small compared to the minimum radius of curvature of the ellipse, or $\sigma\ll(1-e^{2})b$. If this is the case, we can approximate the intrinsic scatter distribution with a convolution of the ellipse and a normal distribution, which in the ellipse-aligned coordinate chart yields

[TABLE]

where $\delta$ is the Dirac delta distribution. While this integral also resists evaluation in closed form, we can approximate it by a mixture of normal distributions, by distributing individual normal distributions on points equidistantly spaced in arc length. The number $N$ of component distributions required to make the distribution smooth depends on the intrinsic scatter $\sigma$. A Euclidean separation of less than $\sigma$ is sufficient, and we can thus set

[TABLE]

where the brackets denote the ceiling function. We can now complete step 6, and evaluate the likelihood for a single measurement, yielding

[TABLE]

where the eccentric anomalies $E_{i}$ need to be numerically solved from

[TABLE]

The equation (61) is then straightforward to evaluate numerically for each datapoint, completing step 7.

3.5 Circles on a spherical surface

Finally, we consider an example of a relation defined in a non-Euclidean space. Specifically, we assume the measurements lie on $M=S^{2}$, the two-dimensional spherical shell. We use the standard spherical coordinate chart $(\theta,\phi)\in\mathbb{R}^{2}$, for which the metric is given by the matrix

[TABLE]

This choice is appropriate for positions on the plane of the sky, for example. This definition completes step 1. For step 2, we choose a relation that defines a circle, by setting

[TABLE]

where $\bm{\bm{\theta}}=(\theta_{c},\phi_{c},\alpha_{c})$. Here the parameters $(\theta_{c},\phi_{c})\in\mathbb{R}^{2}$ specify the centre of the circle, and $\alpha_{c}$ is equal to one-half of the angular size of the circle. The tilde denotes a coordinate transformation which takes the centre of the circle to the north pole, or $\tilde{\theta}_{c}=0$. This transformation is most conveniently realized by first transforming the points $(\theta,\phi)$ to Cartesian three-vectors, using the rotation matrix constructed in Section 3.1, and then converting back to spherical coordinates. The analytic form of the transformation $(\theta,\phi)\mapsto(\tilde{\theta},\tilde{\phi})$ can be given directly in terms of $(\theta,\phi)$ and $(\theta_{c},\phi_{c})$ but the result is unwieldy and omitted.

For steps 3 and 4 we now define the intrinsic coordinates and the coordinate transformation between them and the tilde-transformed coordinates. A convenient choice is to set $\tilde{\phi}=\tau$, so that the $\tau$ coordinate increases with longitude as we go along the circle. We likewise set $\tilde{\theta}=\alpha_{c}+\nu$. This coordinate transformation is multivalued in the sense that a given point $(\tilde{\theta},\tilde{\phi})$ corresponds to a infinite number of intrinsic coordinate pairs, which can be written as

[TABLE]

The second set of intrinsic coordinate pairs above corresponds to geodesics originating from the point of the circle on the other side of the sphere, at $\tilde{\phi}+\pi$ (modulo $2\pi$).

For the definition of the intrinsic scatter distribution in step 5 we again use the one-dimensional normal distribution, equation (6). At this point, it is advantageous to express the intrinsic scatter distribution in the original coordinate frame as well. Since the coordinate transformation is now multivalued, we must proceed as explained in Section 2.2, to get

[TABLE]

where $\vartheta$ is a Jacobi theta function corresponding to the definition

[TABLE]

with $q,z\in\mathbb{C}$, $|q|<1$. Note that there is no dependence on $\tilde{\phi}$. From this form of the $p_{\text{int}}$ distribution it is easy to appreciate the close relation to the wrapped normal distribution defined on the circle, also defined through a Jacobi theta function. As such, this example shows how the wrapped distributions naturally arise from the intrinsic coordinate formalism.

To proceed further, we need to make some assumptions of the error distribution of the measurements, which we denote by pairs $(\tilde{\theta}_{m},\tilde{\phi}_{m})$. An attractive choice would be one of the spherical generalizations of the normal distribution, such as the Kent distribution or the Von Mises–Fisher distribution. However, these combined with the intrinsic scatter distribution (66) and the metric of $S^{2}$ do not seem to easily yield closed form results for the likelihood in step 6. Instead, we have to assume that the measurement error is very small, and approximately given by a uniform distribution within a small cell $[\tilde{\theta}_{m}\pm\Delta_{m}/2,\tilde{\phi}_{m}\pm\Delta_{m}/(2\sin\tilde{\theta}_{m})]\subset S^{2}$, where $\Delta_{m}$ specifies the (small) magnitude of the error.

With the help of this admittedly severe assumption, we can proceed to step 6 and compute the likelihood for a single measurement, finding

[TABLE]

to second order in $\Delta_{m}$. While the likelihood (68) is complicated by the presence of the theta functions and their derivatives, these can be efficiently evaluated numerically, and it is again straightforward to numerically compute the value of the likelihood for each datapoint to complete step 7.

We can compare the result to some existing approaches, such as presented in Jupp & Kent (1987) and Fujiki & Akaho (2009). The former method is based on splines, and the latter is based on Euclidean approximation of small spherical distances and least squares fitting. The method in Jupp & Kent (1987) can in principle cope with arbitrarily large normally distributed measurement errors, but cannot incorporate intrinsic scatter, limits or truncation. In addition, since the result is always a spline, no parameter estimation for a predetermined curve is possible. The method in Fujiki & Akaho (2009) is suitable for parameter estimation, but otherwise shares the disadvantages of the Jupp & Kent (1987) method in addition to not incorporating measurement errors at all. In this sense, the geometric approach presented here compares favourably, although it presents mathematical difficulties if curves other than circles are to be fit.

3.6 What is learned from the examples

The examples clearly demonstrate that analytic results are only easy to obtain in the case where $M$ is Euclidean and the relation to be fit is linear, or more exactly a geodesic. In the last two examples, simplifying assumptions are necessary to obtain a closed form for the likelihood. As such, while the linear Euclidean case of normally distributed orthogonal scatter is essentially solved here completely, much future work is required in the context on non-linear relations and non-Euclidean spaces.

In addition, the choice of parameterization of the relation to be fit has some important consequences. Firstly, the chosen parameterization directly affects the form of the prior probability distribution, as discussed in Section 2.2. In addition, different parameterizations may not be numerically well-behaved everywhere. For example, the $\bm{p}$-parameterization used in Sections 3.1 and 3.2 will likely exhibit numerical instability if any of the relation hyperplanes passes close to the coordinate origin, in which case $\left\|\bm{p}\right\|\rightarrow 0$ and the orientation of the plane becomes indeterminate. The $(\alpha,\beta)$-parameterization in Section 3.3 suffers from similar difficulties for nearly vertical lines, as then the values $\beta\rightarrow\pm\infty$ become degenerate.

4 The $M_{\text{BH}}$–$\sigma$ relation

As a typical application, the methodology described above can be applied to the $M_{\text{BH}}$–$\sigma$ relation, typically written in the form

[TABLE]

The $M_{\text{BH}}$–$\sigma$ relation is an important correlation between the mass $M_{\text{BH}}$ of the supermassive black hole in the centre of a galaxy and the velocity dispersion $\sigma$ of the galactic bulge.444To avoid confusion, the intrinsic scatter will be denoted with $\sigma_{\text{int}}$ in this section.

Since the initial discovery of this correlation (Ferrarese & Merritt 2000; Gebhardt et al. 2000, but see also Magorrian et al. 1998), it has been re-established several times, using both larger datasets and different statistical procedures (e.g. Tremaine et al., 2002; Novak et al., 2006; Gültekin et al., 2009; McConnell et al., 2011; Graham et al., 2011; Beifiori et al., 2012; McConnell & Ma, 2013; Saglia et al., 2016; van den Bosch, 2016). The existence of the $M_{\text{BH}}$–$\sigma$ relation and analogous relations, such as the correlation $M_{\text{BH}}\text{--}M_{\text{bulge}}$ with the galaxy bulge mass, have also been confirmed in numerical simulations, both in galaxy merger simulations (e.g. Di Matteo et al., 2005; Johansson et al., 2009a, b; Choi et al., 2014) and in cosmological simulations (e.g. Sijacki et al., 2007; Di Matteo et al., 2008; Booth & Schaye, 2009; Sijacki et al., 2015). However, all of these studies have used statistical methods which treat one observable as independent and the other observable as dependent. Redoing the fit with the independent observable as the dependent and vice versa produces a significantly different correlation slope, and in some cases also affects the estimated value of the intrinsic scatter (Park et al., 2012). Consequently, there is some controversy in the astronomical literature as to whether it is more suitable to fit $M_{\text{BH}}$ as a function of $\sigma$ (forward regression, using the definition in Park et al. 2012), or the other way around (inverse regression, respectively) in the presence of intrinsic scatter. For a review of the debate, see Graham (2016) and the references therein. That different slopes are produced by switching dependent and independent variables has been known for a long time (e.g. Pearson, 1901), and amounts to asking two different questions: what is the most likely value for $M_{\text{BH}}$ given $\sigma$ or vice versa. However, it is equally well known that if one is interested in the functional relation between the observables, then a symmetric method should be used (see e.g. Isobe et al., 1990).

The approach in Section 3.3 presents a symmetric Bayesian solution to fitting a linear relation between two observables, incorporating heteroscedastic errors, upper limits and intrinsic scatter. As such, we would expect it to yield a non-biased estimate of the intrinsic slope $\beta$ of the relation (69), regardless of which way the relation is fit. To investigate this, the three datasets used in Table 1 of Park et al. (2012) were analysed. The data are originally from Gültekin et al. (2009), McConnell et al. (2011) and Graham et al. (2011). In order to compare with the results in Park et al. (2012), the measurement errors were modelled as uncorrelated bivariate Gaussians in the logarithmic space, leaving out all upper limits. Following Park et al. (2012), the standard deviations were set equal to mean errors, i.e. $\sigma_{\log_{10}M_{\text{BH}}}=(\log_{10}M_{\text{BH,high}}-\log_{10}M_{\text{BH,low}})/2$ and similarly for the velocity dispersions. In addition, for the Graham et al. (2011) sample, which lacks velocity dispersion error data, the velocity dispersion errors were set to 10% and then propagated to averaged logarithmic errors. Finally, a dataset compiled from a multitude of sources used in van den Bosch (2016) was used. For this dataset, the errors were used as given, and upper limits were also incorporated in the fit.

For each dataset, the values for the same set of parameters as in Park et al. (2012) were computed: the intercept $\alpha$, slope $\beta$ and intrinsic scatter along the $M_{\text{BH}}$-axis $\sigma_{\text{int},M_{\text{BH}}}$. The results were computed for forward regression, equation (69), as well as the inverse regression, converting back to equivalent forward values in the end. The parameter values and uncertainties were estimated using the posterior distribution, equation (37), in two complementary ways.

The first set of estimates (Geo-MAP, hereafter), were derived by numerically maximizing the posterior distribution combined with a bootstrap resampling procedure. For each dataset, a total of $n_{d}(\log n_{d})^{2}$ bootstrap samples were constructed (Babu & Singh, 1983; Feigelson & Babu, 2012), where $n_{d}$ is the number of data points. The parameters were then estimated using the median parameter values. Parameter uncertainties at 1-$\sigma$ level were estimated using median absolute deviations scaled to correspond to standard deviations for a normal distribution.

For the second set of estimates (Geo-MCMC, hereafter), the posterior distribution was sampled with the Markov chain Monte Carlo (MCMC) sampler emcee (Foreman-Mackey et al., 2013). Convergence during sampling was monitored with the potential scale reduction factor $R$ (Gelman et al., 2013), and sampling was continued until $R<1.01$ was achieved. The parameter estimates were obtained as the values corresponding to the sample with maximum posterior probability. The parameter uncertainties at 1-$\sigma$ level were then estimated by constructing credible regions containing $0.6827$ of the posterior probability mass. This was done by starting from the maximum posterior probability sample, and descending in posterior probability until the limit was exceeded. The extent of the credible region in each parameter direction was then used to compute the upper and lower 1-$\sigma$ limits.

The results of this analysis are displayed in Table 2. Shown also are the results from Park et al. (2012) containing both forward and inverse regressions with the FITEXY (Press et al., 1992; Tremaine et al., 2002) and ‘Bayesian’ (i.e. linmix_err) methods. In addition, fits for the van den Bosch (2016) dataset without upper limits were computed separately for all methods. A graphical representation of the datasets and the relations obtained with the different methods is shown in Figure 5.

Table 2 indicates that the results computed using the geometrical approach presented in this paper are more consistent with the inverse fits done using the linmix_err and FITEXY algorithms. This is not surprising in light of the discussion in the previous section, since the slope of the $M_{\text{BH}}$–$\sigma$ relation is high and we expect the forward fits to be biased towards low slopes in this situation. The effects of this bias can also be seen from the best-fitting values of the orthogonal intrinsic scatter, $\sigma_{\text{int}}$, for which the forward fits yield consistently higher values than the inverse fits for linmix_err and FITEXY. It seems that in these methods the lower value for the slope is compensated by a higher estimate for the intrinsic scatter. Note that this behaviour cannot be appreciated by looking at the values of the scatter in $M_{\text{BH}}$-direction, $\sigma_{\text{int},M_{\text{BH}}}$, since these are not truly intrinsic, but depend on $\beta$ and indeed show opposite behaviour. In contrast, the geometric methods yield identical results (up to sampling noise) in forward and inverse directions, so only forward results are shown in Table 2. The geometric methods also give consistently smaller estimates for the orthogonal intrinsic scatter. The 1-$\sigma$ errors for all parameters are in general comparable between the methods, with the exception of the fully Bayesian MCMC approach, which consistently gives more conservative 1-$\sigma$ errors.

Finally, it seems that if there indeed is a fundamental approximately linear relation between the logarithms of the mass of a supermassive black hole mass and the velocity dispersion of its host galaxy, the slope of the relation is likely to be $\gtrsim 6$, at least based on the van den Bosch (2016) data, which is the most comprehensive dataset used here. This is in contrast to the values around $4$–$5$ often obtained in the literature (e.g. $4.8$, Ferrarese & Merritt 2000; $3.75$, Gebhardt et al. 2000; $4.02$, Tremaine et al. 2002; $4.24$, Gültekin et al. 2009; $5.64$, McConnell & Ma 2013), including the value $5.35$ derived in van den Bosch (2016). These values were all derived using a method or a variation of a method described in Section 3.3.1, fitting the $M_{\text{BH}}$–$\sigma$ relation in the ‘forward’ direction ($M_{\text{BH}}$ as a function of $\sigma$), in which it has a high numerical value for the slope. This then causes the fitted value of the slope to be biased towards lower values. However, as can be seen from Figure 5, the numerically rather different estimates of the slope are not visually at all that obvious, with a change in slope of $5$ to $6$ corresponding only to a \sim$$$ change in angle of the regression line with respect to a fixed direction, such as the x$-axis.

The obtained slope also crucially depends on the observational and other biases the data may have. Indeed, it is evident from Table 2 that the choice of dataset has an effect that is roughly comparable to the effect of the choice of method. It should also be noted that the fundamental relation, if there is one, may involve more than two physical quantities or their measurable proxies, as suggested in van den Bosch (2016). Fitting relations to a sampled projection of this hyperplane would then in general yield a higher value for the intrinsic scatter and slope that is offset from the ‘true’ value. As such, there is an urgent need for more high-quality data in order to say much with any certainty regarding the slope of the $M_{\text{BH}}$–$\sigma$ relation (or the possible multi-variable generalizations), its possible evolution with redshift, or whether it truly is linear across the entire range of black hole masses.

5 Summary

We have presented a mathematical formalism for representing physical relations as submanifolds $S$ of a Riemannian manifold of observables $M$. In this geometric approach, intrinsic scatter in the relation can be accommodated with probability distributions defined on the normal spaces of $S$. Given a probability distribution of data, and parameterizations of the relation and the intrinsic scatter distribution, the formalism then yields a Bayesian posterior probability for the parameters of the relation and the intrinsic scatter distribution, equation (11). The novelty of our formulation is that it fully accommodates arbitrary measurement errors, both left and right censored data (upper and lower limits, respectively), truncation (non-detections) and extends the concept of intrinsic scatter both to non-linear relations and relations that define a submanifold of codimension greater than one.

We have derived explicit analytic results for the likelihood and the posterior distribution first in the case where the postulated relation defines a linear $n-1$-dimensional hyperplane. We then extended this result to the case where the relation defines an $n-k$-dimensional affine subspace, for an arbitrary $k$, a result we believe to be potentially highly useful for seeking out the most likely codimension of a correlation within a set of $n$-dimensional data. Finally, we have derived the likelihood and posterior distribution in the case of a line in two dimensions, and discussed its implications at length. We also compared the results given by our method with two established methods widely used in astronomical literature, namely FITEXY and linmix_err. We demonstrated that our inherently symmetrical geometrical approach is preferable in situations where the data obeys a relation with a slope much larger or smaller than one, and measurement errors and intrinsic scatter are severe and equally important.

Finally, we used our method to fit the $M_{\text{BH}}$–$\sigma$ relation, between the mass $M_{\text{BH}}$ of a supermassive black hole in a galactic bulge, and the stellar velocity dispersion $\sigma$ of the bulge, using several published datasets. We compared our results to the fits in the literature, and find that our results support a slope of $\sim 6$, clearly higher than the slopes $\sim 4\text{--}5$ derived in the literature. We note that this difference is mainly due to the methods used to derive the literature results. We show that if these methods are used ‘in reverse’, to fit $\sigma$ as a function of $M_{\text{BH}}$, and then inverting the slope, the results are in much better agreement with ours. This is due to the tendency of standard methods, such as FITEXY or linmix_err, which do not respect the geometric symmetry of the problem, to misestimate steep slopes in the presence of intrinsic scatter (see Section 3.3.1 for the discussion).

Acknowledgements

The author is most grateful to the anonymous referee for comments and ideas that have contributed to a much improved manuscript, as well as to Peter H. Johansson for extensive comments on the manuscript drafts. The numerical computations in this work have benefited from the Python libraries NumPy, SciPy and SymPy (van der Walt et al., 2011; Jones et al., 01; Meurer et al., 2017). The Figures have been rendered with the help of the Python library Matplotlib (Hunter, 2007). This work has made use of NASA’s Astrophysics Data System Bibliographic Services. The research for this publication was supported by the Academy of Finland grant no. 1274931.

Appendix A Derivation of the posterior distribution

We present a derivation of the likelihood, equation (10), and the posterior distribution, equation (11), following the Bayesian style promoted in Jaynes (2003). We consider the following propositions:

• •

$\eta$ = ‘the true value of the observable is $\eta\in M$’

• •

$X$ = ‘a value $x\in M$ was measured for the observable’

• •

$V$ = ‘a detection was made’

• •

$H$ = ‘the true value of the observable is drawn from the distribution $p_{\text{int}}$ defined on a relation $S$, with parameters $\bm{\varphi}$ and $\bm{\theta}$, respectively’

In addition, we will use $I$ to specify all the other relevant prior information. This includes the prior distribution $\pi_{\bm{\theta},\bm{\varphi}}$ of $\bm{\theta}$ and $\bm{\varphi}$, the measurement error distribution $h(\eta;x)$ and the fact that the probability of a detection for a value $\eta$ of the observable is given by $p_{\text{det}}(\eta)$.

At the outset, we then know the following probabilities

[TABLE]

We will need the probability $P(X|\eta VI)$ of a measured value, given a true value and the fact that there is a detection. This is obtained with Bayes’ theorem

[TABLE]

as a function of $x$ (which we write as $h(x;\eta)$ in the following), since the prior probabilities $P(X|VI)$ and $P(\eta|VI)$ must be uninformative and equal everywhere, since nothing in $I$ tells us where the true and measured values are a priori.

We can now compute the probability of a measured value given the true value, yielding

[TABLE]

where $\bar{V}$ is the negation of $V$ (i.e. there was no detection). Here we used the fact that $V$ and $\bar{V}$ form a complete ($P(V+\bar{V}|I)=1$) and independent set of propositions, together with the identity $P(AB|C)=P(A|C)P(B|AC)=P(B|C)P(A|BC)$.

Using Bayes’ theorem again we can now obtain the probability $P(H|XI)$ that the data obey the relation and are drawn from $p_{\text{int}}$, given the measured value. The theorem gives

[TABLE]

Since the $\eta$ also form a complete and independent set of propositions (the true value must be somewhere, and the possible positions are independent), we can write

[TABLE]

which gives the likelihood, equation (10). The factor $P(H|I)=\pi_{\bm{\theta},\bm{\varphi}}(\bm{\theta},\bm{\varphi})$ is the parameter prior probability, and the denominator is a normalizing constant, given formally by

[TABLE]

We then have the posterior probability, equation (11),

[TABLE]

Appendix B An example of limit distributions

Here we present a simple example of how the framework presented in the paper works with upper limits, and how analytic results can be obtained for upper limits as well. This is desirable since typically the analytic expression for the likelihood is much faster to compute than performing a numerical integration to obtain the required value. As mentioned in Section 2.3, limits are taken into account by including them into the measurement error distribution. This can be done for the examples in Sections 3.1–3.3 by substituting one or more degrees of freedom $x_{i}$ in equations (16) and (32) with suitable one-dimensional distributions. Conceptually the simplest possibility is the uniform distribution, so that

[TABLE]

where $\chi_{[a,b]}(x)=1$ if $x\in[a,b]$ and [math] otherwise, and $x_{i}^{u}$ is the limiting value. Other choices with less pronounced cutoffs are naturally also possible. If the measurements are already given in logarithmic units with base $k$, the form

[TABLE]

must be used instead. Lower limits can be introduced in an analogous manner.

Analytic results of reasonable complexity based on these limit distributions can be obtained for the example in Section 3.3. For the linear case we have

[TABLE]

[TABLE]

for the upper limits in $x$ and $y$ directions respectively, where $x^{u}$ and $y^{u}$ are the corresponding limiting values. If the measurements have been given in logarithmic scale with base $k$, we have instead

[TABLE]

[TABLE]

where $\operatorname{erfc}(x)=1-\operatorname{erf}(x)$. Note that here the limiting values $x^{u}$ and $y^{u}$ are given in linear scale. For numerical applications, it should be noted that the arguments of the $k$-exponential and the $\operatorname{erfc}$ function may have large numerical values, and asymptotic expansions should be used when necessary.