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\).
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
Harris, G.L.H. et al (2013). “Globular clusters and supermassive black holes in galaxies: further analysis and a larger sample.” Monthly Notices of the Royal Astronomical Society, 438, 2117-2130.
Merritt, D. (2000). “Black holes and galaxy evolution. Dynamics of Galaxies: from the Early Universe to the Present”, eds. F. Combes, G. A. Mamon, and V. Charmandaris Vol. 197. Astronomical Society of the Pacific Conference Series, p. 221.