  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.
