SCOPE: Symmetric COvariance Population Estimator

SCOPE (Symmetric COvariance Population Estimator) is a hierarchical Bayesian regression framework implemented in R and Stan. It fits linear scaling relations between two observed quantities while fully accounting for measurement uncertainties in both coordinates, including asymmetric (skewed) uncertainties modeled via a skew-normal likelihood.

Companion Paper:
Alister W. Graham (2026), "Galaxy morphology dependent (black hole mass)-(velocity dispersion) relations: implications for gravitational wave forecasts and cosmological simulations"
https://arxiv.org/abs/2606.05808

Source Code & Repository:
https://github.com/A-Graham/SCOPE

Scientific & Statistical Overview:
Traditional ordinary least squares (OLS) and one-sided conditional regressions introduce significant mathematical bias when measurement errors exist in the independent variable (X). SCOPE models the intrinsic population as a bivariate normal distribution, treating both variables symmetrically during fitting. The reported regression slope is derived directly from the intrinsic population covariance:

    beta = rho * (sigma_Y / sigma_X)

with conditional intrinsic scatter:

    sigma_{Y|X} = sigma_Y * sqrt(1 - rho^2)

Key Features & Implementation:
• Symmetrical Error Treatment: Symmetrically fits joint intrinsic distributions, avoiding OLS regression dilution/bias.
• Asymmetric Error Handling: Automatically converts split (+/- 1 sigma) measurement uncertainties into skew-normal likelihood parameters, correctly evaluating p(true value | measurement).
• Robust MCMC Sampling: Implemented in Stan via rstan with non-centred parameterization and an LKJ prior on the correlation matrix.
• Automated Diagnostics & Publication Plots: Generates publication-ready figures (SCOPE_fit.pdf) and saves full posterior chains (Scout.dat) for downstream analysis.
• Broad Applicability: Tailored for astronomical scaling relations (e.g. M_BH-sigma, M_BH-M_bulge, luminosity-velocity dispersion), but generalizable to any two-variable dataset with two-dimensional errors.