BAYESIAN ANALYSIS OF TWO STELLAR POPULATIONS IN GALACTIC GLOBULAR CLUSTERS. I. STATISTICAL AND COMPUTATIONAL METHODS

Bayesian methods offer a principled, probability-based approach for combining information from the current data and our prior knowledge. They require a likelihood function—the distribution of the data given the model parameters. The likelihood function is the primary statistical tool for assessing the viability of a parameter value vis-à-vis the observed data under a postulated statistical model. The knowledge we have about the model parameters before considering the current data is specified in a prior distribution. Past and current information are combined in the posterior distribution of the parameters, which is related to the likelihood function and the prior distribution through Bayes' theorem. With generic data and model parameters represented by Y and ψ, Bayes' theorem gives the posterior distribution as

$$\begin{eqnarray}&&P(\psi | Y)=\displaystyle \frac{P(Y| \psi )P(\psi )}{P(Y)},\end{eqnarray} \tag{ 1 }$$

where $P(Y| \psi )\equiv L(\psi | Y)$ is the likelihood function and $P(\psi )$ the prior distribution. The term P(Y), sometimes called the "evidence," is a normalizing constant which makes $P(\psi | Y)$ a proper probability distribution. The posterior distribution provides a summary of the combined information in the data and our prior knowledge, and can be used to derive parameter estimates and uncertainties.

To build a Bayesian model for two-population globular clusters, we start by defining necessary notation and terminology in Section 2.2. In Section 2.3 we construct a preliminary likelihood function for a simple stellar population that accounts for measurement error, the presence of field stars, and the possibility of binary star systems. We extend this model to allow for two-stellar populations in Section 2.4, and specify the full prior distribution in Section 2.5.

For each star in a data set we obtain calibrated photometric magnitudes using at least two filters. Following DeGennaro et al. (2009), van Dyk et al. (2009) and Stein et al. (2013), we refer to the observed photometric magnitude in filter j for star i as x_ij for $j=1,\ldots ,n$ and $i=1,\ldots ,N$, where N is the number of stars in the data set and n is the number of filters. The observed photometric magnitudes for star i are tabulated in the column vector ${{\boldsymbol{X}}}_{i}={({x}_{i1},\ldots ,{x}_{{in}})}^{\top }$, and the known (independent) Gaussian measurement errors in the (diagonal) variance-covariance matrix, ${{\boldsymbol{\Sigma }}}_{i}$.

As discussed in Section 1, our statistical model is based on a hierarchy of parameters. We refer to the parameters that are common to cluster stars—specifically age, metallicity, distance, and absorption—as cluster parameters. These parameters are collected in the vector ${\boldsymbol{\Theta }}=({\theta }_{\mathrm{age}},{\theta }_{[\mathrm{Fe}/{\rm{H}}]},{\theta }_{m-{M}_{V}},{\theta }_{{A}_{V}})$. We refer to the parameters that are common to all stars belonging to a stellar population, but that vary from population to population within a cluster, as population parameters. We assume that only helium abundance differs between the populations; helium abundance and the proportion of stars for population k, denoted ${\phi }_{{Yk}}$ and ${\phi }_{{pk}}$, respectively, are the population parameters. (When discussing simple stellar populations we denote the single helium abundance with ${\phi }_{Y}$.) We refer to the population with the lower helium abundance as "Population 1," and that with the higher abundance as "Population 2." (This should not be confused with the traditional use of Population I versus Population II stars.) As a result, assuming that the two stellar populations result from two epochs of star formation, Population 1 corresponds to the first generation of stars and Population 2 to the second generation. For now the only stellar parameter specific to star i is its initial mass, M_i. (Two more are specified below.) The computer-based stellar evolution model, ${\boldsymbol{G}}$, takes $({M}_{i},{\boldsymbol{\Theta }},{\phi }_{Y})$ and outputs a $1\times n$ vector of predicted photometric magnitudes for a star with those parameters. We express the vector of predicted magnitudes as ${\boldsymbol{G}}({M}_{i},{\boldsymbol{\Theta }},{\phi }_{Y})$. For this study, ${\boldsymbol{G}}$ are the updated Dotter et al. (2008) models that include HST UV magnitudes.

Before considering the likelihood function for a two-population globular cluster, we first consider a "preliminary" likelihood function for a simple stellar population. Following van Dyk et al. (2009), we account for unresolved binaries because the added luminosity of a binary companion shifts a star off the main sequence on the CMD, which can result in systematic errors if not properly handled. We thus treat every observed star as a possible binary system and fit its primary initial mass, M_i, and ratio of the secondary and primary initial masses, ${R}_{i}\leqslant 1;$ a unitary system is expected to have a mass ratio near zero. Because stellar luminosities sum, and magnitudes are on a log-luminosity scale, the predicted magnitudes for (binary) star i are

$$\begin{eqnarray}&&{{\boldsymbol{\mu }}}_{i}=-2.5\,{\mathrm{log}}_{10}({10}^{-{\boldsymbol{G}}({M}_{i},{\boldsymbol{\Theta }},{\phi }_{Y})/2.5}+{10}^{-{\boldsymbol{G}}({M}_{i}{R}_{i},{\boldsymbol{\Theta }},{\phi }_{Y})/2.5}).\end{eqnarray} \tag{ 2 }$$

Owing to the nature of the stellar evolution models tabulated in ${\boldsymbol{G}}$, ${{\boldsymbol{\mu }}}_{i}$ is a complex nonlinear function of the underlying parameters.

We account for field star contamination by introducing indicator variables ${\boldsymbol{Z}}=({Z}_{1},\ldots ,{Z}_{N})$, where Z_i = 1 if star i is a cluster star, and Z_i = 0 if star i is a field star, following the example of van Dyk et al. (2009). These variables allow us to specify a different statistical model for the photometric magnitudes of cluster stars versus those of field stars. We model the observed photometric magnitudes of cluster stars as n-dimensional multivariate Gaussian distributions, such that

$$\begin{eqnarray}&&P({{\boldsymbol{X}}}_{i}| {{\boldsymbol{\Sigma }}}_{i},{M}_{i},{R}_{i},{\boldsymbol{\Theta }},{\phi }_{Y},{Z}_{i}=1)\,\\ &&\displaystyle =\frac{1}{\sqrt{{(2\pi )}^{n}| {{\boldsymbol{\Sigma }}}_{i}| }}\exp \left(-\displaystyle \frac{1}{2}{({{\boldsymbol{X}}}_{i}-{{\boldsymbol{\mu }}}_{i})}^{\top }{{\boldsymbol{\Sigma }}}_{i}^{-1}({{\boldsymbol{X}}}_{i}-{{\boldsymbol{\mu }}}_{i})\right).\end{eqnarray} \tag{ 3 }$$

While this model appears simple at first glance, it is actually quite complex due to the dependence of ${{\boldsymbol{\mu }}}_{i}$ on the stellar evolution parameters, and the complex interdependencies therein. That ${\boldsymbol{G}}$ cannot be expressed in closed form yields challenges for inference and computation.

Following van Dyk et al. (2009), we specify a simple model for field stars that does not depend on any of the parameters of interest. Each field star may have its own values for ${\theta }_{\mathrm{age}}$, ${\theta }_{[\mathrm{Fe}/{\rm{H}}]}$, ${\theta }_{m-{M}_{V}}$, ${\theta }_{{A}_{V}}$, and ${\phi }_{Y}$, and we cannot fit these parameters. We therefore simply assume that each field star magnitude is uniformly distributed over the range of the data, such that

$$\begin{eqnarray*}&&P({{\boldsymbol{X}}}_{i}| {Z}_{i}=0)=c\quad \mathrm{if}\,{\min }_{j}\,\leqslant \,{x}_{{ij}}\leqslant {\max }_{j},\,j=1,\ldots ,n,\end{eqnarray*}$$

and zero elsewhere, where $({\min }_{j},\,{\max }_{j})$ is the range of magnitude values for filter j, and $c={[{\prod }_{j=1}^{n}({\max }_{j}-{\min }_{j})]}^{-1}$. We could instead incorporate a more complex and realistic model for field stars; properties of field stars for specific Galactic fields exist and may assist in tuning the model (e.g., Robin et al. 2012). However, our work to date has not necessitated the additional effort because the simple model adequately identifies field stars. This is illustrated using a simulation study in Section 4.2.

A preliminary likelihood function for a simple stellar population can now be written,

$$\begin{eqnarray}&&{L}_{{\rm{p}}}({\boldsymbol{M}},{\boldsymbol{R}},{\boldsymbol{Z}},{\boldsymbol{\Theta }},{\phi }_{Y}| {\boldsymbol{X}},{\boldsymbol{\Sigma }})\\ &&=\displaystyle \prod _{i=1}^{N}\left[{Z}_{i}\times \displaystyle \frac{1}{\sqrt{{(2\pi )}^{n}| {{\boldsymbol{\Sigma }}}_{i}| }}\exp \left(-\displaystyle \frac{1}{2}{({{\boldsymbol{X}}}_{i}-{{\boldsymbol{\mu }}}_{i})}^{\top }\right.\right.\\ &&\quad \times {{\boldsymbol{\Sigma }}}_{i}^{-1}({{\boldsymbol{X}}}_{i}-{{\boldsymbol{\mu }}}_{i})\phantom{\rule{0ex}{2.90ex}})+(1-{Z}_{i})\times P({{\boldsymbol{X}}}_{i}| {Z}_{i}=0)\phantom{\rule{0ex}{3.60ex}}]\\ &&=\displaystyle \prod _{i=1}^{N}[{Z}_{i}\times P({{\boldsymbol{X}}}_{i}| {{\boldsymbol{\Sigma }}}_{i},{M}_{i},{R}_{i},{\boldsymbol{\Theta }},{\phi }_{Y},{Z}_{i}=1)\\ &&\quad +(1-{Z}_{i})\times P({{\boldsymbol{X}}}_{i}| {Z}_{i}=0)],\end{eqnarray} \tag{ 4 }$$

where ${\boldsymbol{M}}=({M}_{1},\ldots ,{M}_{N})$, ${\boldsymbol{R}}=({R}_{1},\ldots ,{R}_{N})$, ${\boldsymbol{X}}=({{\boldsymbol{X}}}_{1},\ldots ,{{\boldsymbol{X}}}_{N})$, and ${\boldsymbol{\Sigma }}=({{\boldsymbol{\Sigma }}}_{1},\ldots ,{{\boldsymbol{\Sigma }}}_{N})$. The sum in Equation (4) represents the fact that the sample of stars is a mixture of two subgroups: cluster stars and field stars; such distributions are known as finite mixture distributions in the statistics literature. Interested readers are referred to Andreon & Weaver (2015) for a review of the application of mixture models in astronomy.

Rather than embedding ${\boldsymbol{G}}$ into a statistical likelihood function as we do in Equation (4), the computer model can be accounted for using a computational approach known as Approximate Bayesian Computation (ABC). ABC is typically used in situations where the likelihood function is either unavailable or computationally expensive to evaluate, but forward simulation of synthetic data under the statistical model is relatively fast (e.g., Ishida et al. 2015). While synthetic data can be easily generated under the model in Equation (4), constructing a distance measure for comparing observed and synthetic data that accounts for the known Gaussian measurement errors, binary star systems, and field star contamination would be a challenge.

We now extend $P({{\boldsymbol{X}}}_{i}| {{\boldsymbol{\Sigma }}}_{i},{M}_{i},{R}_{i},{\boldsymbol{\Theta }},{\phi }_{Y},{Z}_{i}=1)$ in Equations (3) and (4) to account for the fact that the sample of cluster stars is itself a mixture of two subgroups that are the two stellar populations. This results in a model with three subgroups: field stars and two cluster populations. The likelihood function for a two-population cluster is then

$$\begin{eqnarray}&&L({\boldsymbol{M}},{\boldsymbol{R}},{\boldsymbol{\Theta }},{\boldsymbol{\Phi }},{\boldsymbol{Z}}| {\boldsymbol{X}},{\boldsymbol{\Sigma }})\\ &&\,\,=\displaystyle \prod _{i=1}^{N}\left[{Z}_{i}\times \displaystyle \sum _{k=1}^{2}{\phi }_{{pk}}P({{\boldsymbol{X}}}_{i}| {{\boldsymbol{\Sigma }}}_{i},{M}_{i},{R}_{i},{\boldsymbol{\Theta }},{\phi }_{{Yk}},{Z}_{i}=1)\right.\\ &&\,\,\quad +(1-{Z}_{i})\times P({{\boldsymbol{X}}}_{i}| {Z}_{i}=0)\phantom{\rule{0ex}{3.60ex}}].\end{eqnarray} \tag{ 5 }$$

Evaluating Equation (5) involves computing the expected photometry for each star as if it were a member of each population, i.e.,

$$\begin{eqnarray}{{\boldsymbol{\mu }}}_{{ik}} & = & -2.5\,{\mathrm{log}}_{10}({10}^{-{\boldsymbol{G}}({M}_{i},{\boldsymbol{\Theta }},{\phi }_{{Yk}})/2.5}\\ & & +{10}^{-{\boldsymbol{G}}({M}_{i}{R}_{i},{\boldsymbol{\Theta }},{\phi }_{{Yk}})/2.5}),\end{eqnarray} \tag{ 6 }$$

for $i=1,...,N$, k = 1, 2. The population proportions in Equation (5) must sum to one: ${\phi }_{p1}+{\phi }_{p2}=1$.

A key advantage to adopting a Bayesian approach is its ability to directly incorporate previous (independent) results through the joint prior distribution, which we specify via a set of independent priors on each parameter. For example, we construct a prior distribution on initial mass that is derived from the Miller & Scalo (1979) initial mass function. In particular, we specify a Gaussian prior distribution on the ${\mathrm{log}}_{10}$ of primary initial masses:

$$\begin{eqnarray}&&P({\mathrm{log}}_{10}({M}_{i}))\propto \exp \left(-\displaystyle \frac{1}{2}{\left(\displaystyle \frac{{\mathrm{log}}_{10}({M}_{i})+1.02}{0.677}\right)}^{2}\right),\end{eqnarray} \tag{ 7 }$$

truncated to 0.1 M_⊙ to 8 M_⊙, where the numerical constants are taken from Miller & Scalo (1979). For the ratio of the secondary and primary masses we use a uniform prior distribution on $[0,1]$. We need not truncate the lower end of the secondary mass because low secondary masses indicate that the star is a unitary system (van Dyk et al. 2009).

For the cluster parameters, ${\boldsymbol{\Theta }}$, we incorporate ancillary information to specify informative (i.e., narrow) prior distributions when available, and use relatively diffuse prior distributions when such information is lacking. In particular, for ${\theta }_{[\mathrm{Fe}/{\rm{H}}]}$, ${\theta }_{m-{M}_{V}}$, and ${\theta }_{{A}_{V}}$, we use Gaussian prior distributions (truncated to be positive in the case of ${\theta }_{{A}_{V}}$), with means set according to previously published values from independent data sets and standard deviations chosen to be reasonably large. For age, ${\theta }_{\mathrm{age}}$, we use a uniform prior distribution truncated to the reasonable range of 1–15 Gyr, which includes all Galactic globular clusters.

Because the population parameters, ${\boldsymbol{\Phi }}=({\phi }_{Y1},{\phi }_{Y2},$ ${\phi }_{p1},{\phi }_{p2})$, are the primary parameters of scientific interest, we use uniform prior distributions subject to physical constraints on their ranges. A uniform prior distribution on the interval $[0.15,0.3]$ is used for ${\phi }_{Y1};$ this bounds the helium fraction between 15% and 30%. Similarly, a uniform prior distribution on the interval $[0.15,0.4]$ is used for ${\phi }_{Y2}$, and we impose the constraint ${\phi }_{Y2}\gt {\phi }_{Y1}$. Because we do not typically have prior knowledge for the proportion of stars in each population, ${\phi }_{p1}$ is given a uniform prior distribution on the interval $[0,1]$. When such prior knowledge is available we advocate using a more general beta prior distribution⁷ ; a uniform distribution on the interval $[0,1]$ is equivalent to a beta (1, 1) distribution. We do not need to specify a prior distribution for ${\phi }_{p2}$ because ${\phi }_{p2}=1-{\phi }_{p1}$.

Ancillary measurements (e.g., proper motions) can be used to probabilistically separate field stars from cluster stars. When such ancillary measurements are unavailable, we use $P({Z}_{i}=1)=\alpha$ for $i=1,...,N$, where α is based on the expected fraction of cluster stars in the data set. As we show in Section 5, our results are not sensitive to reasonable choices of α.

With the full joint posterior distribution denoted by $P({\boldsymbol{\Theta }},{\boldsymbol{\Phi }},{\boldsymbol{M}},{\boldsymbol{R}},{\boldsymbol{Z}}| {\boldsymbol{X}})$, the marginal posterior distribution of $({\boldsymbol{\Theta }},{\boldsymbol{\Phi }})$ is given by

$$\begin{eqnarray}P({\boldsymbol{\Theta }},{\boldsymbol{\Phi }}| {\boldsymbol{X}}) & = & \displaystyle \int \cdots \displaystyle \int \\ & & \times \left(\displaystyle \sum _{{Z}_{1}}\cdots \displaystyle \sum _{{Z}_{N}}P({\boldsymbol{\Theta }},{\boldsymbol{\Phi }},{\boldsymbol{M}},{\boldsymbol{R}},{\boldsymbol{Z}}| {\boldsymbol{X}})\right)d{\boldsymbol{M}}d{\boldsymbol{R}}\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&\propto P({\boldsymbol{\Theta }},{\boldsymbol{\Phi }})\displaystyle \prod _{i=1}^{N}\phantom{\rule{0ex}{2.60ex}}[c({\boldsymbol{\Theta }},{\boldsymbol{\Phi }})P({Z}_{i}=1)\\ &&\quad +P({{\boldsymbol{X}}}_{i}| {Z}_{i}=0)P({Z}_{i}=0)\phantom{\rule{0ex}{2.60ex}}],\end{eqnarray} \tag{ 9 }$$

where $c({\boldsymbol{\Theta }},{\boldsymbol{\Phi }})$

$$\begin{eqnarray}&&=\displaystyle \iint \displaystyle \sum _{k=1}^{2}{\phi }_{{pk}}P({{\boldsymbol{X}}}_{i}| {{\boldsymbol{\Sigma }}}_{i},{M}_{i},{R}_{i},{\boldsymbol{\Theta }},{\phi }_{{Yk}},{Z}_{i}=1)\\ &&\,\,\times P({M}_{i},{R}_{i}| {Z}_{i}=1){{dM}}_{i}{{dR}}_{i}.\end{eqnarray} \tag{ 10 }$$

When $P({Z}_{i}=1)=\alpha$ for $i=1,...,N$ (i.e., when all stars have the same probability of being cluster stars), Equations (8) and (9) reduce to $P({\boldsymbol{\Theta }},{\boldsymbol{\Phi }}| {\boldsymbol{X}})\propto$

$$\begin{eqnarray}&&P({\boldsymbol{\Theta }},{\boldsymbol{\Phi }})\displaystyle \prod _{i=1}^{N}\phantom{\rule{0ex}{2.60ex}}[(1-\alpha )P({{\boldsymbol{X}}}_{i}| {Z}_{i}=0)\\ &&\,\,+\displaystyle \iint \alpha \displaystyle \sum _{k=1}^{2}{\phi }_{{pk}}P({{\boldsymbol{X}}}_{i}| {{\boldsymbol{\Sigma }}}_{i},{M}_{i},{R}_{i},{\boldsymbol{\Theta }},{\phi }_{{Yk}},{Z}_{i}=1)\\ &&\,\,\times P({M}_{i},{R}_{i}| {Z}_{i}=1){{dM}}_{i}{{dR}}_{i}\phantom{\rule{0ex}{2.60ex}}].\end{eqnarray} \tag{ 11 }$$

This integral cannot be evaluated analytically because $P({{\boldsymbol{X}}}_{i}| {M}_{i},{R}_{i},{\boldsymbol{\Theta }},{\phi }_{{Yk}},{Z}_{i}=1)$ depends on M_i and R_i through ${\boldsymbol{G}}$ (Stein et al. 2013). Instead, we employ brute-force numerical integration via Riemann sums. By marginalizing out the $3N$ stellar parameters we reduce the dimension of the posterior distribution from typically thousands to just seven.

Because the remaining parameter vector $({\boldsymbol{\Theta }},{\boldsymbol{\Phi }})$ after marginalizing out ${\boldsymbol{M}}$, ${\boldsymbol{R}}$, and ${\boldsymbol{Z}}$ is just seven-dimensional, we initially implemented a standard Metropolis algorithm to sample from $P({\boldsymbol{\Theta }},{\boldsymbol{\Phi }}| {\boldsymbol{X}})$. However, we found the trial-and-error approach to tuning the (seven-dimensional) proposal distribution to be difficult; this is not surprising given the correlations among the components of $({\boldsymbol{\Theta }},{\boldsymbol{\Phi }})$ in ${\boldsymbol{G}}$. To avoid arduous fine-tuning and make BASE-9 more accessible to users less familiar with MCMC, we implement an Adaptive Metropolis (AM) algorithm (e.g., Haario et al. 2001; Roberts & Rosenthal 2009; Rosenthal et al. 2011, p. 93). Whereas an iteration of a standard Metropolis algorithm only depends on the most recent value in the MCMC chain, an AM algorithm uses the entire history of the chain to adapt the proposal distribution at each iteration. This, however, violates a defining property of a Markov chain: the distribution of a value in the chain can only depend on the history of the chain through its most recent value. Thus, care must be taken to guarantee an AM algorithm converges properly. As recounted in Rosenthal et al. (2011, p. 93), AM algorithms must satisfy the Diminishing Adaptation Condition: the amount of adaptation at iteration l must go to 0 as $l\to \infty$. The Diminishing Adaptation Condition is key; other technical conditions are almost always satisfied except in specially constructed examples (Rosenthal et al. 2011, p. 93). Readers interested in additional mathematical details are encouraged to consult Rosenthal et al. (2011, p. 93) and references therein.

After marginalizing over ${\boldsymbol{M}}$, ${\boldsymbol{R}}$, and ${\boldsymbol{Z}}$, the resulting marginal posterior distribution $P({\boldsymbol{\Theta }},{\boldsymbol{\Phi }}| {\boldsymbol{X}})$ given in Equation (11) appears roughly Gaussian. This is illustrated in Figure 2, which displays the matrix of two-dimensional scatter plots of 25,000 posterior draws from $P({\boldsymbol{\Theta }},{\boldsymbol{\Phi }}| {\boldsymbol{X}})$. (The data we used to construct these plots are photometric magnitudes from NGC 5272. Details are provided in Section 5.) Based on results in Gelman et al. (1996), the optimal proposal distribution for a Gaussian posterior distribution with a d-dimensional covariance matrix ϒ is itself a Gaussian distribution with covariance matrix $[{(2.38)}^{2}/d]{\rm{\Upsilon }}$. Because the actual form of the posterior distribution is unknown, for iteration $l+1$ we use a multivariate t proposal distribution with six degrees of freedom⁸ , centered at the current value of $({\boldsymbol{\Theta }},{\boldsymbol{\Phi }})$ and with scale equal to $[{(2.38)}^{2}/7]{{\boldsymbol{\Xi }}}^{(l)}$. Here, ${{\boldsymbol{\Xi }}}^{(l)}$ is the empirical variance-covariance matrix of $\{({{\boldsymbol{\Theta }}}^{(1)},{{\boldsymbol{\Phi }}}^{(1)}),...,({{\boldsymbol{\Theta }}}^{(l)},{{\boldsymbol{\Phi }}}^{(l)})\}$. Because we recalculate ${{\boldsymbol{\Xi }}}^{(l)}$ at every iteration, the proposal distribution adapts at every iteration based on the past history of the chain. As $l\to \infty$, the empirical distribution of $\{({{\boldsymbol{\Theta }}}^{(1)},{{\boldsymbol{\Phi }}}^{(1)}),...,({{\boldsymbol{\Theta }}}^{(l)},{{\boldsymbol{\Phi }}}^{(l)})\}$ approaches the marginal posterior distribution $P({\boldsymbol{\Theta }},{\boldsymbol{\Phi }}| {\boldsymbol{X}})$, improving efficiency. Furthermore, ${{\boldsymbol{\Xi }}}^{(l)}$ stabilizes and thus the adaptation diminishes as required.

Zoom In Zoom Out Reset image size

Alternative modern MCMC approaches include Hamiltonian Monte Carlo (HMC) (see, e.g., http://mc-stan.org) and Riemann manifold Monte Carlo (RMMC) methods (Girolami & Calderhead 2011). Such methods are particularly useful when the posterior distribution exhibits strong correlations and curving degeneracies. HMC, for example, borrows ideas from Hamiltonian dynamics to make sampling more efficient, but typically requires the likelihood to be available analytically, and that its derivatives with respect to the model parameters be available. Similar caveats apply to RMMC. While we could develop an analytical emulator of the function ${\boldsymbol{G}}$ to deploy HMC or RMMC, the additional effort is unnecessary due to the roughly Gaussian shape of $P({\boldsymbol{\Theta }},{\boldsymbol{\Phi }}| {\boldsymbol{X}});$ the AM algorithm automatically improves sampling efficiency by adapting the proposal distribution.

When implementing the AM algorithm, we first run the sampler in "tuning" mode. The goal of this tuning period is not to obtain an optimal proposal distribution, but rather to sufficiently explore the posterior distribution and generate a reasonable ${{\boldsymbol{\Xi }}}^{(1)}$ for the AM algorithm (see the Appendix for details). Once we have calculated ${{\boldsymbol{\Xi }}}^{(1)}$ from the tuning period, the first 1000 iterations of the AM algorithm use the multivariate t proposal distribution described above, with ${{\boldsymbol{\Xi }}}^{(l)}={{\boldsymbol{\Xi }}}^{(1)}$ for $l=1,...,1000$. This non-adaptive period is necessary to generate a sufficiently large sample to estimate posterior covariances before adapting the proposal distribution. At iteration 1001, and at every subsequent iteration, ${{\boldsymbol{\Xi }}}^{(l)}$ is the empirical covariance matrix of the previous l iterations.

The efficiency of our AM algorithm is demonstrated in Figure 3. There, we compare the performance of our AM algorithm to that of a standard Metropolis algorithm for sampling the same posterior distribution. The standard Metropolis sampler is identical to the AM sampler except ${{\boldsymbol{\Xi }}}^{(l)}$ is fixed at ${{\boldsymbol{\Xi }}}^{(1)}$ throughout. Both algorithms are implemented with the same starting values, the same tuning period, and use the same data as in Figure 2. The proposal distribution in the AM algorithm begins adapting at iteration 1001, after which there is an obvious difference in performance between AM and standard Metropolis. Standard Metropolis struggles as the MCMC chain repeatedly becomes stuck. While AM also sticks initially, adapting the proposal distribution quickly frees the chain and leads to increased efficiency. As expected, the AM algorithm becomes increasingly efficient as the number of iterations increases.

Zoom In Zoom Out Reset image size

As an initial test of our method, we simulate two-population globular clusters under three scenarios, with ten replicate clusters per scenario. The three scenarios differ in the percentage of stars belonging to Population 1: 50%, 80%, and 100% for Scenario 1, 2, and 3, respectively. Scenario 3 therefore contains ten replicates of a single-population cluster, which we intentionally fit with our (incorrect) two-population model to demonstrate how model misspecification can potentially be identified. Each cluster is simulated with ${\theta }_{\mathrm{age}}=10.08$, ${\theta }_{m-{M}_{V}}=15.375$, ${\theta }_{{A}_{V}}=0.372$, and ${\theta }_{[\mathrm{Fe}/{\rm{H}}]}=-1.5$, which are "average" published values across the clusters compiled in Harris (1996, and updated in 2010 at http://physwww.mcmaster.ca/~harris/mwgc.ref). In the simulations we set ${{\phi }_{Y}}_{1}=0.24$ and ${{\phi }_{Y}}_{2}=0.29$, so that the true difference in helium abundance is 0.05. We simulate 30,000 cluster stars and 1000 field stars per cluster, and every star is generated as a single-star system. (Future work will include binaries.) For each cluster we generate photometric magnitudes in five filters, corresponding to the filters in HST UVIS photometry (Piotto et al. 2015): $F275W$, $F336W$, $F438W$, $F606W$, and $F814W$. Details about these filters are provided in Section 5. The photometric magnitudes for each star are simulated with uncorrelated Gaussian measurement error that is a function of both the wavelength band and the magnitude, as depicted in Figure 4.

Zoom In Zoom Out Reset image size

Of the 31,000 stars generated per cluster, about 90% are dropped from the simulated cluster for one of two reasons. First, there is a threshold signal-to-noise ratio that eliminates stars too dim to be observed under realistic conditions. Second, we believe the stellar evolution models are inaccurate for fainter main sequence stars and we therefore impose a magnitude cutoff on real photometry (DeGennaro et al. 2009; van Dyk et al. 2009). We impose the same cutoff on simulated photometry so that our simulation results are as informative as possible. The exact cutoff we use depends on the assumed distance to the cluster (see Section 5). In the simulations, we discard stars with a photometric magnitude in the $F275W$ filter greater than 23. After losing stars due to low signal-to-noise and the $F275W$ magnitude cutoff, about 3000 simulated stars remain per cluster.

We use prior distributions for ${\theta }_{m-{M}_{V}}$, ${\theta }_{{A}_{V}}$, and ${\theta }_{[\mathrm{Fe}/{\rm{H}}]}$ with means equal to their true values under the simulation; the prior standard deviations were set to 0.05, 0.124, and 0.05, respectively. The prior distributions on the population parameters are as described in Section 2.5. We assign $P({Z}_{i}=1)=\alpha =0.95$ for $i=1,...,N$. This is the value for α we use when analyzing NGC 5272 in Section 5 after testing the sensitivity to the choice of α.

We use our AM algorithm to explore $P({\boldsymbol{\Theta }},{\boldsymbol{\Phi }}| {\boldsymbol{X}})$ for each of the 30 simulated clusters. We run one chain per cluster for 25,000 iterations after the tuning period. Inspection of the trace plot for each chain shows that all the chains reach their apparent stationary distributions within the first 5000 iterations. We discard the first 5000 iterations of each chain as burn-in and base inference on the remaining 20,000 iterations. Results for the three scenarios appear in Figure 5.

Zoom In Zoom Out Reset image size

The results for Scenarios 1 and 2 are presented in the top four rows of Figure 5. There, we observe that our method is performing reasonably well with respect to recovering the difference in helium abundance and the proportion of stars in each population. This is encouraging, as our main inferential goal is to recover the difference in helium abundance. Unfortunately, there is a systematic difference between the fitted parameters and the true values of the parameters under the simulation. The reasons for this discrepancy are examined and discussed in detail in Section 4.3 of Stenning (2015). It was discovered that the deviations increase with the size of the measurement errors, suggesting an influence of the prior distribution. The cause is that as the sample size increases, the influence of the prior distribution on primary initial mass does not diminish because there is only one observation (i.e., one star) per mass parameter. Future work will focus on fitting the distribution of the masses to hopefully eliminate the deviations. For now, we simply note that the systematic deviations are small relative to both the systematic errors stemming from the underlying stellar evolution model and to the best available statistical errors on these parameters using other methods. For example, minimum star-by-star [Fe/H] statistical and systematic errors are approximately 0.015 and 0.03 dex, respectively (Carretta et al. 2009). Typical statistical errors for distance moduli are $\sigma (m-{M}_{V})$ = 0.1 mag, and those for absorption are σ(${A}_{{\rm{V}}}$) = $0.1{A}_{{\rm{V}}}$, with a lower limit of 0.03 mag (Harris 1996, and as updated at http://physwww.mcmaster.ca/~harris/mwgc.ref). Furthermore, we can adequately recover the relative difference in helium abundance because the systematic differences are in the same direction and to a similar degree for both populations.

To check that the recovered helium abundance difference is due to the presence of two populations and is not an artifact of the method, in Scenario 3 we intentionally fit simulated single-population clusters with the (now incorrect) two-population model. The results for the cluster parameters are similar to those in Scenarios 1 and 2; see row 5 of Figure 5. This is expected because the cluster parameters are common to both populations. However, in the last row of Figure 5, the results for the proportion of stars in Population 1 indicate model misspecification. Specifically, we observe that the fitted value is close to either zero or one (replicates 2, 3, 8, and 10) and/or the 95% interval is very wide, spanning most of the range from $[0,1]$ (replicates 1, 2, 4, 5, 6, 7, 9, and 10). Both of these outcomes suggest that a second population may not be present in the data.

We caution that investigating intervals and estimates in this way does not provide a formal diagnostic for model misspecification but results such as those under Scenario 3 should be considered as "smoking-gun" evidence that the two-population model has been applied to a single-population cluster. For now, we intend our model to be used in cases where there are two prominent populations as viewed in CMDs; e.g., two populations for NGC 5272 can be seen in the rightmost CMD in Figure 6. Formal criteria to infer the number of populations in a cluster will be included when our model and algorithms are extended to accommodate three or more populations.

Zoom In Zoom Out Reset image size

To test the adequacy of using the simple uniform model described in Section 2.3 for field star magnitudes, we simulate five replicate single-population clusters with parameters equal to those reported for NGC 5272 in the updated Harris (1996) globular cluster catalog. Field stars are simulated from the Besançon model (Robin et al. 2003) with Galactic $l,b$ = 42.2170, +78.7069, though we increase the field size $100\times$ (from 0.0013 to 0.130 square degrees) relative to that for NGC 5272 to provide an ample sample in each replication. After removing stars due to low signal-to-noise and an imposed magnitude cutoff at $F275W=22.074$ (see Section 5 for a discussion regarding our choice of cutoff), there are approximately 2100 cluster stars and 110 field stars per replicate data set. Using BASE-9, we are able to infer the posterior probability that a star is a cluster member; see Stein et al. (2013). Those stars with greater than 50% posterior probability are classified as cluster stars. We can evaluate the resulting classification using the confusion matrices in Table 2. A confusion matrix is a table with columns representing true classifications and rows representing predicted classifications. For example, the confusion matrix for Replication 1 reveals that 111 field stars are correctly identified as such, while one field star is misclassified as a cluster star; all 2103 clusters stars in Replication 1 are correctly classified. In this simulation, the simple model for field star magnitudes misidentifies field stars as cluster stars $\lt 2 \%$ of the time, and never misidentifies cluster stars as field stars. Based on these results, a more complex model for field stars seems unnecessary.