Example 2: Normal linear model for accessing the relationship between central black hole mass and bulge velocity dispersion

Our first practical application illustrates the use of a Bayesian Gaussian model for describing the relation between the mass \(M_\bullet\) ofthe supermassive black hole present in the center of a galaxy and the velocity dispersion \(\sigma_\epsilon\) of the stars in its bulge.

The study of the global properties of individual galaxies has attracted attention in the astronomical community.
Merritt (2000) the so-called \(M_\bullet-\sigma_\epsilon\) relation and took the form
\[ \frac{M_\bullet}{10^8 M_{\odot}} \approx 3.1 \times \left( \frac{\sigma_\epsilon}{200 km s^{-1}} \right)^4 \] originally known as the Faber-Jackson law for black holes. Here \(M_{\odot}\) is the mass of the Sun. In this example we follow Harris et al (2013) who model the relation as: \[ \log \left( \frac{M_\bullet}{M_{\odot}} \right) = \alpha + \beta \log\left( \frac{\sigma_\epsilon}{\sigma_0} \right) \] where \(\alpha_0\) is a reference value, usually chosen to be 200 km/s, and \(\alpha\) and \(\beta\) are the linear predictor coefficients to be determined.

The data set used in this example was presented by Harris et al (2013). This is a compilation of literature data from a variety of sources obtained with the Hubble Space Telescope as well as with a wide range of other ground-based facilities. The original data set was composed of data from 422 galaxies, 46 of which have available measurements of \(M_\bullet\) and \(\sigma_\epsilon\).

The logarithm of the observed central black hole mass, \(M \equiv \log(M_\bullet/M_\odot)\), represents our response variable and the logarithm of the observed velocity dispersion is our explanatory variable, \(\sigma \equiv \log(\sigma_\epsilon/\sigma_0)\). Their true values are denoted by \(M^{\mbox{true}} \equiv \log(M_\bullet^{\mbox{true}}/M_\cdot)\) and \(\sigma^{\mbox{true}} \equiv \log(\sigma_\epsilon^{\mbox{true}}/\sigma_0)\) respectively.

We want to take into account the measurement errors in both variables \((\varepsilon_M, \varepsilon_\sigma)\). In the following statistical model, \(N = 46\) is the total number of galaxies, \(\eta = \mu = \alpha + \beta \sigma\) is the linear predictor, \(\alpha\) is the intercept, \(\beta\) is the slope, and \(\varepsilon\) is the intrinsic scatter around \(\mu\). Our model is

\[ \begin{array}{rl} M_i & \sim N(M_i^{\mbox{true}},\varepsilon_{M,i}^2) \\ M_i^{\mbox{true}} & \sim N(\mu_i,\varepsilon^2) \\ \mu_i & = \alpha + \beta \dot \sigma_i \\ \sigma_i & \sim N(\sigma_i^{\mbox{true}},\varepsilon_{\sigma,i}^2) \\ \end{array} \] for \(i = 1,\ldots, N\). We assign non-informative priors to \(\alpha\), \(\beta\), \(\sigma\), and \(\varepsilon\) with \[ \begin{array}{rl} \alpha & \sim N(0,10^3) \\ \beta & \sim N(0,10^3) \\ \varepsilon^2 & \sim \mbox{Gamma}(10^{-3},10^{-3}) \end{array} \]

Load rstan and set options.

library(rstan)
## Warning: package 'rstan' was built under R version 3.4.3
## Loading required package: ggplot2
## Loading required package: StanHeaders
## Warning: package 'StanHeaders' was built under R version 3.4.3
## rstan (Version 2.16.2, packaged: 2017-07-03 09:24:58 UTC, GitRev: 2e1f913d3ca3)
## For execution on a local, multicore CPU with excess RAM we recommend calling
## rstan_options(auto_write = TRUE)
## options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

Read and extract the variables into R.

# Data
path_to_data = "./M_sigma.csv"

# Read data
MS<-read.csv(path_to_data,header = T)

# Identify variables
N <- nrow(MS) # number of data points
obsx <- MS$obsx # log black hole mass
errx <- MS$errx # error on log black hole mass
obsy <- MS$obsy # log velocity dispersion
erry <- MS$erry # error on log velocity dispersion

# Prepare data  
MS_data <- list(
  obsx = obsx,
  obsy = obsy,
  errx = errx,
  erry = erry,
  N = N
)

Specify the rstan model.

# Stan Gaussian model
stan_code="
data {
    int<lower=0> N;          // number of data points
    vector[N] obsx;          // obs velocity dispersion
    vector<lower=0>[N] errx; // errors in velocity dispersion measurements
    vector[N] obsy;          // obs black hole mass
    vector<lower=0>[N] erry; // errors in black hole mass measurements
}
parameters {
    real alpha;            // intercept
    real beta;             // angular coefficient
    real<lower=0> epsilon; // scatter around true black hole mass
    vector[N] x;           // true velocity dispersion
    vector[N] y;           //  true black hole mass
}
model {
    // likelihood and priors
    alpha ~ normal(0, 1000);
    beta ~ normal(0, 1000);
    epsilon ~ gamma(0.001, 0.001);
    for (i in 1:N){
        x[i] ~ normal(0, 1000);
        y[i] ~ normal(0, 1000);
    }
    obsx ~ normal(x, errx);
    y ~ normal(alpha + beta*x, epsilon);
    obsy ~ normal(y, erry);
}
"

Fit the model using rstan

oldcat <- cat()
cat <- function(...) return(invisible(NULL))

# Run mcmc
fit = stan(model_code=stan_code, data=MS_data, iter=15000, chains=4, warmup=5000)
## In file included from C:/Users/John/Documents/R/win-library/3.4/BH/include/boost/config.hpp:39:0,
##                  from C:/Users/John/Documents/R/win-library/3.4/BH/include/boost/math/tools/config.hpp:13,
##                  from C:/Users/John/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/var.hpp:7,
##                  from C:/Users/John/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
##                  from C:/Users/John/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core.hpp:12,
##                  from C:/Users/John/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/mat.hpp:4,
##                  from C:/Users/John/Documents/R/win-library/3.4/StanHeaders/include/stan/math.hpp:4,
##                  from C:/Users/John/Documents/R/win-library/3.4/StanHeaders/include/src/stan/model/model_header.hpp:4,
##                  from file22e8597c51b1.cpp:8:
## C:/Users/John/Documents/R/win-library/3.4/BH/include/boost/config/compiler/gcc.hpp:186:0: warning: "BOOST_NO_CXX11_RVALUE_REFERENCES" redefined
##  #  define BOOST_NO_CXX11_RVALUE_REFERENCES
##  ^
## <command-line>:0:0: note: this is the location of the previous definition
## In file included from C:/Users/John/Documents/R/win-library/3.4/BH/include/boost/multi_array/base.hpp:28:0,
##                  from C:/Users/John/Documents/R/win-library/3.4/BH/include/boost/multi_array.hpp:21,
##                  from C:/Users/John/Documents/R/win-library/3.4/BH/include/boost/numeric/odeint/util/multi_array_adaption.hpp:29,
##                  from C:/Users/John/Documents/R/win-library/3.4/BH/include/boost/numeric/odeint.hpp:61,
##                  from C:/Users/John/Documents/R/win-library/3.4/StanHeaders/include/stan/math/prim/arr/functor/integrate_ode_rk45.hpp:13,
##                  from C:/Users/John/Documents/R/win-library/3.4/StanHeaders/include/stan/math/prim/arr.hpp:38,
##                  from C:/Users/John/Documents/R/win-library/3.4/StanHeaders/include/stan/math/prim/mat.hpp:298,
##                  from C:/Users/John/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/mat.hpp:12,
##                  from C:/Users/John/Documents/R/win-library/3.4/StanHeaders/include/stan/math.hpp:4,
##                  from C:/Users/John/Documents/R/win-library/3.4/StanHeaders/include/src/stan/model/model_header.hpp:4,
##                  from file22e8597c51b1.cpp:8:
## C:/Users/John/Documents/R/win-library/3.4/BH/include/boost/multi_array/concept_checks.hpp: In static member function 'static void boost::multi_array_concepts::detail::idgen_helper<N>::call(Array&, const IdxGen&, Call_Type)':
## C:/Users/John/Documents/R/win-library/3.4/BH/include/boost/multi_array/concept_checks.hpp:42:43: warning: typedef 'index_range' locally defined but not used [-Wunused-local-typedefs]
##        typedef typename Array::index_range index_range;
##                                            ^
## C:/Users/John/Documents/R/win-library/3.4/BH/include/boost/multi_array/concept_checks.hpp:43:37: warning: typedef 'index' locally defined but not used [-Wunused-local-typedefs]
##        typedef typename Array::index index;
##                                      ^
## C:/Users/John/Documents/R/win-library/3.4/BH/include/boost/multi_array/concept_checks.hpp: In static member function 'static void boost::multi_array_concepts::detail::idgen_helper<0ull>::call(Array&, const IdxGen&, Call_Type)':
## C:/Users/John/Documents/R/win-library/3.4/BH/include/boost/multi_array/concept_checks.hpp:53:43: warning: typedef 'index_range' locally defined but not used [-Wunused-local-typedefs]
##        typedef typename Array::index_range index_range;
##                                            ^
## C:/Users/John/Documents/R/win-library/3.4/BH/include/boost/multi_array/concept_checks.hpp:54:37: warning: typedef 'index' locally defined but not used [-Wunused-local-typedefs]
##        typedef typename Array::index index;
##                                      ^
## In file included from C:/Users/John/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core.hpp:42:0,
##                  from C:/Users/John/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/mat.hpp:4,
##                  from C:/Users/John/Documents/R/win-library/3.4/StanHeaders/include/stan/math.hpp:4,
##                  from C:/Users/John/Documents/R/win-library/3.4/StanHeaders/include/src/stan/model/model_header.hpp:4,
##                  from file22e8597c51b1.cpp:8:
## C:/Users/John/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/set_zero_all_adjoints.hpp: At global scope:
## C:/Users/John/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/set_zero_all_adjoints.hpp:14:17: warning: 'void stan::math::set_zero_all_adjoints()' defined but not used [-Wunused-function]
##      static void set_zero_all_adjoints() {
##                  ^
cat <- oldcat

Print the summary statistics of the parameter posterior distributions and plot the fitted distributions.

print(fit, pars = c("alpha", "beta", "epsilon", "lp__"))
## Inference for Stan model: 3427dba4a84c9dfe554684bec1385ea9.
## 4 chains, each with iter=15000; warmup=5000; thin=1; 
## post-warmup draws per chain=10000, total post-warmup draws=40000.
## 
##          mean se_mean   sd   2.5%    25%   50%   75% 97.5% n_eff Rhat
## alpha    8.35     0.0 0.06   8.24   8.31  8.35  8.39  8.46 40000    1
## beta     4.46     0.0 0.33   3.80   4.24  4.46  4.68  5.09 28058    1
## epsilon  0.27     0.0 0.06   0.16   0.23  0.27  0.31  0.40 10101    1
## lp__    -9.29     0.1 9.03 -27.00 -15.20 -9.32 -3.49  8.71  7671    1
## 
## Samples were drawn using NUTS(diag_e) at Thu Dec 07 07:17:23 2017.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).
plot(fit, pars = c("alpha", "beta", "epsilon", "lp__"))
## ci_level: 0.8 (80% intervals)
## outer_level: 0.95 (95% intervals)

pairs(fit, pars = c("alpha", "beta", "epsilon", "lp__"))

Using another method Harris et al (2013) reported parameter values are \(\alpha = 8.412 \pm 0.067\) and \(\beta = 4.610\pm 0.403\).

References

From: Bayesian Models for Astrophysical Data, Cambridge Univ. Press (c) 2017, Joseph M. Hilbe, Rafael S. de Souza and Emille E. O. Ishida

Code 10.1 - Normal linear model in R using JAGS for accessing the relationship between central black hole mass and bulge velocity dispersion