Bayesian Inference of Globular Cluster Properties Using Distribution Functions

In this study, we define the likelihood using a physical DF, f( θ ; d_i ), of the limepy lowered-isothermal model. Given a fixed set of data, d , of N stars, the likelihood is a function of the model parameters, θ , and the total mass, M_total, of the GC:

$$\begin{eqnarray}&&{ \mathcal L }({\boldsymbol{\theta }};{\boldsymbol{d}})=\prod _{i=1}^{N}\displaystyle \frac{f({\boldsymbol{\theta }};{d}_{i})}{{M}_{\mathrm{total}}}\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&=\prod _{i=1}^{N}\displaystyle \frac{f(g,{{\rm{\Phi }}}_{0},{M}_{\mathrm{total}},{r}_{h};{r}_{i},{v}_{i})}{{M}_{\mathrm{total}}},\end{eqnarray} \tag{ 4 }$$

where the stars are assumed to be independent.

For lowered-isothermal models, the DF f is calculated numerically via the limepy software (Gieles & Zocchi 2015), and thus the likelihood must be calculated numerically, too.

As mentioned in Section 2, we simulate the position and kinematic data of stars following a limepy model DF with parameters, θ , shown in Table 1, assume the likelihood defined in Equation (4), and define physically motivated informative priors on the model parameters.

Given that the likelihood is defined by the DF that was used to generate the data, we expect to obtain reasonable parameter estimates through inference made from the posterior distribution using Markov Chain Monte Carlo (MCMC) sampling. However, we are also going to impose prior distributions that are at least weakly informative, and so it is good practice to test whether the posterior can still be used to reliably infer the model parameters. Moreover, in the cases where the sampling of stars from the cluster is biased to inside the core or outside the core, we aim to understand how this sampling bias affects parameter inference.

Two advantages of Bayesian inference are the necessity to incorporate meaningful prior information and the requirement to state this explicitly. In order for the DF to correspond to a physically realistic collection of stars in a GC, all model parameters must be greater than zero. Negative parameter values are not allowed by the likelihood, but we also disallow negative parameter values via the priors (this increases efficiency and keeps the limepy model from returning errors).

One reason to use informative priors is that images and studies both within the the Milky Way and around other galaxies provide prior information on quantities like the mass and half-light radius of GCs. For example, GC masses span about an order of magnitude, and most astronomers would be comfortable setting the prior $p({\mathrm{log}}_{10}{M}_{\mathrm{total}})\sim N({\mu }_{M},{\sigma }_{M})$, where the hyperparameters μ_M and σ_M are defined in ${\mathrm{log}}_{10}{M}_{\mathrm{total}}$. This is the prior we choose, and it is also supported by the near-universal GC mass function (Brodie & Strader 2006; Harris 2010, 1996).

The limepy model works in M_total space, so we need to do a change of variables to obtain the prior, p(M_total). Using a change of variables, the prior on M_total is

$$\begin{eqnarray}&&p({M}_{\mathrm{total}})=\displaystyle \frac{N({\mu }_{M},{\sigma }_{M})}{{M}_{\mathrm{total}}\mathrm{ln}10}.\end{eqnarray} \tag{ 5 }$$

The half-light radius is another quantity of GCs for which we have considerable prior information. Images of GCs give an independent estimate of r_h , with a conservative measurement uncertainty of roughly 0.4 pc (e.g., de Boer et al. 2019). In this simulation study, we assume the observer has this prior information and set a truncated normal prior on r_h .

We have considerably less prior information on the values of g and Φ₀, aside from the physically allowable, positive values. For these parameters, we use truncated uniform distributions. In summary, we assume the parameters for the limepy model are distributed as

$$\begin{eqnarray}&&g\sim \mathrm{unif}(0.001,3.5),\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{0}\sim \mathrm{unif}(1.5,14),\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{M}_{\mathrm{total}}\sim \displaystyle \frac{N({\mu }_{M},{\sigma }_{M})}{{M}_{\mathrm{total}}\mathrm{ln}(10)},\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&\mathrm{and}\ \ {r}_{h}\sim N(a,b,{\mu }_{{r}_{h}},{\sigma }_{{r}_{h}}),\end{eqnarray} \tag{ 9 }$$

where μ_M = 5.85 and σ_M = 0.6 (defined in ${\mathrm{log}}_{10}{M}_{\mathrm{total}}$), and the hyperparameters for the lower and upper bounds of r_h are a = 0 and b = 30, respectively. The mean and standard deviation for the r_h parameter (${\mu }_{{r}_{h}}$ and ${\sigma }_{{r}_{h}}$) are chosen to reflect plausible information an observer would have for a given GC. Thus, for the average GCs in our analysis, we try different means, such as ${\mu }_{{r}_{h}}=3.4$, ${\mu }_{{r}_{h}}=3.1$, etc., with ${\sigma }_{{r}_{h}}=0.4$ pc. Our results are insensitive to the choice of the mean, as long as it is not too many standard deviations away from the true value.

Given the limepy model, we have a likelihood function ${ \mathcal L }(g,{{\rm{\Phi }}}_{0},M,{r}_{h};{\boldsymbol{d}})$ for the four unknown parameters, depending on the observed star data, d . Combining the above prior distributions with this limepy likelihood function leads to a posterior distribution, or target posterior density, via Bayes' theorem (Equation (2)),

$$\begin{eqnarray*}\begin{array}{rcl}p(g,{{\rm{\Phi }}}_{0},{M}_{\mathrm{total}},{r}_{h}| {\boldsymbol{d}}) & \propto & p({\boldsymbol{d}}| g,{{\rm{\Phi }}}_{0},{M}_{\mathrm{total}},{r}_{h})\\ & & \times \,p(g)p({{\rm{\Phi }}}_{0})p({M}_{\mathrm{total}})p({r}_{h}),\end{array}\end{eqnarray*}$$

where we assume independent priors. Our goal is to sample from the target distribution p(g, Φ₀, M_total, r_h ∣ d ), and perform inference of the parameter values, the CMP, and the mean-square velocity profile of the GC.

Ultimately, we explore and collect samples of this posterior density using a MCMC algorithm, specifically a version of the standard Metropolis algorithm (Metropolis et al. 1953) that includes automated, finite adaptive tuning (to be discussed later). First, however, we find optimal starting values; we use the differential-evolution optimizer function DEopt from the NMOF package (Schumann 2021; Gilli et al. 2019) in R (R Core Team 2019) to find modal (i.e., argmax) values of the four parameters, and then use these values as the initial state of our MCMC algorithm. Differential evolution was first introduced by Storn & Price (1997), and we refer the reader to this paper for details on the algorithm. This initial step allows an automated selection of good starting values, which helps to overcome the complicated structure of the posterior distribution, thereby making sampling more efficient. Once the starting values are obtained, we run an automated, finite adaptive-tuning method during the burn-in of the Markov chain. To describe the finite adaptive-tuning method, we first provide a brief review of proposal distributions and sampling efficiency.

Sampling a target or posterior distribution using the standard Metropolis algorithm requires a choice of proposal or "jumping" distribution. The latter is used to randomly suggest a new place in parameter space, θ^*, based on the current location, θ_i . Often, this suggestion is done using a normal distribution such that

$$\begin{eqnarray}&&{{\boldsymbol{\theta }}}^{* }={{\boldsymbol{\theta }}}_{i}+{\boldsymbol{Z}},\end{eqnarray} \tag{ 10 }$$

where Z ∼ N(0, Σ). Here, N(0, Σ) is the jumping distribution with a covariance matrix, Σ, set by the user. The value of Σ determines whether, on average, "big jumps" or "small jumps" are attempted from the current location of θ_i . These proposed jumps are either accepted or rejected according to the standard formula in the Metropolis algorithm. The efficiency of the sampling is dependent on the choice of this covariance matrix. For example, if the variance is too small then the algorithm will make jumps that are too small. If the variance is too large, then the algorithm will make jumps that are too large.

Finding a Σ that enables the most efficient sampling is sometimes accomplished through manual tuning: adjusting Σ until the appropriate acceptance rate is achieved. Obviously, this can be a tedious and time-consuming process, especially in the case of multiple parameters. Thankfully, there are methods which automate this task and that are founded in statistical theory.

In this paper, we use an automated, finite adaptive-tuning method during the burn-in of the Markov chain. This adaptive-tuning method is one in which the proposal step sizes are adjusted automatically and iteratively. We obtain a good covariance matrix for the proposal distribution using an adaptive Metropolis algorithm (Haario et al. 2001; Roberts & Rosenthal 2009), which repeatedly updates the Metropolis proposal distribution (i.e., the proposal covariance matrix) based on the empirical covariance of the run so far, in an effort to obtain a proposal covariance matrix equal to about (2.38)² times the target covariance matrix divided by the Markov chain's dimension, which has been shown to be optimal under appropriate assumptions (Roberts & Rosenthal 1997, 2001). Foundational works on the subject of adaptive Metropolis and convergence are found in the statistics literature (Roberts et al. 1997; Haario et al. 2001; Roberts & Rosenthal 2009).

The practice of using the adaptive Metropolis algorithm for an initial run and then fixing the proposal variance for the final run corresponds to "finite adaptation", as in Proposition 3 of Roberts & Rosenthal (2007). We require a minimum of five initial runs to update the proposal variance, but also automatically allow for further iterations as needed to achieve efficient sampling. Almost all of the GCs we analyze take no more than five iterations of the finite adaptive tuning, which takes one to five minutes per cluster on a simple laptop computer.

Once the finite adaptive step is complete, we run a standard Metropolis algorithm using the final (hopefully approximately optimal) proposal distribution found by the adaptive Metropolis step. The final sampling takes less than 15 minutes per cluster to complete. At the end, we discard an initial burn-in period, and take the remaining chain values as a sample from the posterior density.

The above procedure allows us to approximately sample from p(g, Φ₀, M_total, r_h ∣ d ), and hence (a) approximately compute the posterior means and other statistics of the four unknown parameters (g, Φ₀, M_total, r_h ), including Bayesian credible intervals; and (b) calculate a CMP of the GC for every sample from the target distribution.

Very generally, GCs may be classified as having an average, compact, or extended morphology based on their radius, r_h , or may be considered to have high or low concentration based on the value of Φ₀. Additionally, the spatial and kinematic data from stars may be a random sample from everywhere in the cluster, a random sample beyond some radius, or a random sample within some radius. We expect the ability of our method to recover the true mass, CMP, and mean-square velocity profile to depend on both GC morphology and the type of sampling of its stars. Understanding the bias in parameter inference that can occur as a result of biased sampling is important, since in reality we sometimes lack position and kinematic data from the inner or outer regions of the cluster. Thus, we investigate multiple combinations of the aforementioned cases to understand any possible bias.

Table 1 summarize the types of GCs we investigate. In our simulated GCs, all stars have the same brightness and mass, and so the half-light radius corresponds to the half-mass–radius. For each case, we simulate 50 GCs using the parameter values listed in Table 1, and subsample 500 stars either (1) randomly, (2) outside r_cut, or (3) inside r_cut. We choose an r_cut value of 1.5pc mostly for simplicity but also partly because recent work by the HSTPROMO Team indicates that proper motions are most often available for stars within the half-mass–radius (Watkins et al. 2013) but not beyond. Our conservative choice for r_cut is therefore half of the average effective radius of Galactic clusters (excluding very extended clusters with effective radii greater than 10 pc) (Baumgardt & Hilker 2018). We use this same cut-off radius when sampling outer stars (i.e., situation (2) above) as well. In this way, we investigate what happens when data are only available for outer stars (e.g., when Gaia kinematic data are used).

By repeating the analysis on 50 randomly generated GCs, we estimate and examine the coverage probabilities for the Bayesian credible regions for all sampling scenarios, for all GCs listed in Table 1 (Section 4).

For example, for the average cluster, we generate 50 simulated GCs with parameter values g = 1.5, Φ₀ = 5, M_total = 10⁵ M_⊙, and r_h = 3.0 pc, and randomly sample 500 stars from each GC. For each GC, we run the analysis on the subsample of stars, obtaining samples of the target distribution as described in the previous section. Next, we estimate the mean, interquartile range, and 95% credible interval of the posterior distribution using our MCMC samples from the target distribution. After doing this for all 50 average GCs, we count how many times the interquartile ranges and 95% credible intervals cover the true parameter value to estimate the coverage probability. If the Bayesian credible regions are reliable, then the interquartile ranges should cover the true parameter values 50% of the time, and the 95% credible intervals should cover the true parameter values 95% of the time.

The same procedure is repeated for all GC types listed in Table 1. For example, we look at GCs with different half-light radii, reflecting extended and compact clusters. For these clusters we use parameter values of g = 1.5, Φ₀ = 5, M_total = 10⁵ M_⊙, and r_h = 9.0 pc and g = 1.5, Φ₀ = 5, M_total = 10⁵ M_⊙, and r_h = 1.0 pc, respectively. To further explore the parameter space believed to be covered by Galactic GCs, and specifically to explore GCs that are more (less) concentrated, we also look at GCs with a high (low) Φ₀.

Using our estimate of the posterior distribution for a single GC, we can also estimate that GC's CMP. The CMP is an estimate of the mass contained within some distance, r, of the GC. To estimate the CMP, we follow the same procedure as described in Eadie & Jurić (2019), who used this approach to estimate the Milky Way's CMP. For every set of model parameters (g, Φ₀, M_total, r_h ) sampled by our algorithm (i.e., every row of parameter values in the Markov chain), we calculate the CMP determined by the limepy model. Because we have thousands of rows in our Markov chain, we obtain thousands of CMP estimates. These CMPs provide us with a visual and quantitative estimate that can be used to calculate Bayesian credible regions and that can be compared directly to the true CMP of the cluster.

Another quantity of interest that we can estimate using the posterior distribution for a single GC is the mean-square velocity. The mean-square velocity is equal to the sum of the velocity dispersion squared (or the variance) and the square of the mean velocity. The limepy code can calculate the mean-square velocity as a function of radius from the center of the cluster, given a specific set of model parameters. Thus, the estimate of the GC's mean-square velocity profile can be calculated in much the same way as the CMP, using the parameter samples from the posterior distribution.

In all of the GC examples, we assume that we know the complete position and velocity components of the stars. However, in reality we often have incomplete data. For example, we might only have projected measurements on the plane of the sky (i.e., projected distances in the x–y plane, and proper motions). This missing data might influence our mass and mass profile estimates in unexpected ways, and is important to study. In a Bayesian analysis one can treat the missing components as parameters in the model, but this also means that further prior distributions must be set. Given the complexity of the problem, we leave this to future work.

For the cases in which we randomly sample stars from everywhere in the cluster, we find the Bayesian credible regions to be reliable for all five GC types.

As an example, Figure 2 shows the 95% credible intervals (error bars) for each model parameter, for 50 realizations of an average cluster. The true parameter values are shown as vertical blue lines, and the number of times out of 50 that the 95% credible interval of the target distribution overlaps the true value is shown at the top of each panel. We can see that the credible intervals for each parameter reliably contains the true parameter approximately 95% of the time (Figure 2).

Zoom In Zoom Out Reset image size

As a second example, we show a similar plot for the case of the extended GCs (Figure 3). Here too, we find the 95% credible intervals to be reliable for the most part. The credible intervals for g and Φ₀ are slightly overconfident, since the true parameter value lies within the 95% credible intervals only 90% and 92% of the time, respectively.

Zoom In Zoom Out Reset image size

As a final and third example, Figure 4 shows the same type of plot for a more concentrated cluster with Φ₀ = 8. Again, the credible intervals are reliable, showing good coverage probabilities.

Zoom In Zoom Out Reset image size

Table 2 shows the estimated coverage probabilities for the M_total parameter in the case of random sampling, for all five types of clusters, found by calculating the fraction of times that the true M_total is contained within the Bayesian credible interval. We can see that both the 50% and 95% credible intervals for M_total are reliable when the stars are randomly sampled throughout the cluster, despite cluster type.

The MCMC samples can also be used to infer the CMP of the cluster under the limepy model. Figure 5 shows the CMP inferred for one example of an average, compact, extended, low-Φ₀, and high-Φ₀ cluster in the random sampling case. The posterior distribution samples of g, Φ₀, M_total, and r_h from the Markov chains are used to calculate the posterior estimate of the CMP, shown as transparent black curves. The red curve shows the true CMP given by the limepy model with the correct parameters.

Zoom In Zoom Out Reset image size

The CMPs provide not only a visual inspection of our method but also a quantitative one. The posterior curves for a given GC (e.g., the collection of black curves for the average GC in Figure 5) can be used to construct Bayesian credible intervals at all radii (e.g., the teal regions for the average GC in Figure 6). After constructing these credible regions for each GC, we can ask: "How often does the true CMP lie within these credible regions, at different radii?".

Zoom In Zoom Out Reset image size

As an example of this quantitative comparison, we use the results of all 50 realizations of average GCs to calculate the reliability of the CMP 95% credible regions. Table 3 shows how often the true M(r < R) fell within the 95% credible region at 10 logarithmically spaced distances, r, for the average GCs. The results show that the credible regions are reliable, with the true M(r < R) being recovered approximately 95% of the time at all radii.

In general, we find that the credible regions and CMPs are reliable for all types of GCs when the stars are sampled randomly throughout the cluster. It is reassuring that we can recover the true parameter values and the CMPs reliably from a random sample of only 500 stars.

In the case of real data, we will not know the true CMP of a GC. Thus, one might like to check whether the CMP inference is reasonable given the observed data. One way to compare the Bayesian-inferred CMP to the observed data is shown in Figure 6. Here, the 50%, 75%, and 95% credible regions are compared to the empirical cumulative distribution function (ECDF) of the 500 stars' positions, r. Another way to do this kind of comparison or check, which is not done here, would be to perform posterior predictive checks: to simulate data from the posterior distribution, and compare these simulated data to the real data (e.g., see Shen et al. 2022, where Bayesian posterior predictive checks are used to check inferences about the CMP of the Milky Way).

Additionally, we can inspect other physical quantities provided by the limepy model fit. For example, Figure 7 shows the mean-square velocity $\overline{{v}^{2}}$ profiles as a function of radius for one GC in each of the five morphologies. Under random sampling of the stars, we observe that the true mean-square velocity profile is well recovered by the MCMC samples. Similarly to the CMPs discussed above, Bayesian credible regions for velocity profiles could also be calculated, and posterior predictive checks could be performed to compare simulated data from the posterior to the observed data.

Zoom In Zoom Out Reset image size

In general, we find that biased sampling of stars from only inside or outside the cluster core results in model parameter estimates that are biased and in Bayesian credible intervals that are unreliable. While obtaining biased estimates from a biased data sample is not surprising, the reality is that this type of sampling mimics the data from some telescopes. Investigating these cases can illuminate the kind of biases we should expect and possibly correct for. Indeed, through our investigations of biased sampling, we find the success of the parameter inference and CMP inference is a combination of both the cluster's morphology and the type of biased sampling.

As an example, Figure 8 shows the 95% credible intervals for an average GC when only the outer stars' data are sampled. We can see that the credible intervals are unreliable, and that parameter estimates are biased. In particular, M_total, g, and r_h are consistently overestimated, while Φ₀ is underestimated.

Zoom In Zoom Out Reset image size

In contrast, biased sampling of outer stars of an extended cluster result in parameter estimates that are much more reliable (Figure 9). In this case, the extended GC's mass, M_total, and half-light radius, r_h , can actually be estimated reliably.

Zoom In Zoom Out Reset image size

In Table 4, we summarize how reliably we can recover M_total in the biased sampling cases. Only the extended cluster with sampling in the outer regions is reliable. Also note the *'s in the table, which indicate when the MCMC algorithm had trouble finding a stationary distribution with good mixing, leading to biased estimates of the total mass (Table 2). In these particular cases, the behavior of the Markov chain would be a clue to the observer that the model is having trouble describing the data.

For GCs with high Φ₀, biased sampling of stars in the inner regions also leads to poor parameter estimates and unreliable credible regions (Figure 10). For some of the GCs in the scenario, the Markov chains become stuck in one location. The estimates of the mean from these bad chains are shown as the open circles with a small dot in the middle (i.e., the "estimates" have a variance of zero because the chains became stuck at a single place in parameter space). The exact estimated parameter values in these bad cases are rather meaningless and random. Moreover, if a scientist were to see this behavior in a Markov chain from a real data analysis, then they would know not to trust the solution. However, in many cases of randomly generated GCs with high Φ₀ and biased sampling in the inner regions, the Markov chains do look reasonable even when their estimates are not. Thus, a scientist could mistakenly assume the convergence is giving reliable parameter estimates. We will return to this scenario shortly.

Zoom In Zoom Out Reset image size

As mentioned in the previous section,the CMPs provide more insight than simply looking at the parameter estimates and their credible intervals, both visually and quantitatively. Figure 11 shows example CMPs for each GC morphology when the stars in these GCs are sampled only in their outer or inner regions (first and second column, respectively). Looking at the first column in Figure 11, we see that when stars are sampled outside the core, the inner region of the cluster's profile tends to be underestimated, regardless of the GC morphology. The opposite is true for sampling inside the core (the second column). At the same time, sampling outside the core tends to lead to an overestimate of the total mass, while sampling inside the core leads to a (sometimes severe) underestimate.

Zoom In Zoom Out Reset image size

There are two exceptions to the observation that biased samples lead to biased CMPs, namely (i) when extended and low-Φ₀ clusters are sampled in the outer regions, and (ii) when compact clusters are sampled in the inner regions. For the extended and low-Φ₀ GC, our method is able to recover the true CMP reasonably well when stars outside the core are sampled, whereas this is certainly not the case when stars inside the core are sampled. For the compact GC, we see the opposite case: the CMP is reasonably well estimated when the sample contains stars inside the core versus outside the core.

These cases where biased samples still lead to unbiased estimates are not surprising: sampling stars in the outer region of an extended or less-concentrated cluster will provide a better representation of the true stellar distribution than sampling stars in its core, because these types of GCs are less dense in their inner regions (Figure 1). Likewise, sampling stars in the inner region of a compact cluster will be a better representation of the true stellar distribution than a sample from the outer region because compact GCs are more dense toward their centers.

Next, we test the reliability of the CMP Bayesian credible regions. As an example, we show the results for low-Φ₀ GCs when stars are sampled outside the core (Table 5). It is clear that the 95% c.r. at all radii are unreliable, with the inner regions being the most unreliable.

Next, we use the MCMC samples to estimate the mean-square velocity $\overline{{v}^{2}}$ profile as a function of radius. In Figure 12, each row corresponds to a specific GC type, and the columns indicate whether stars were sampled outside (left) or inside (right) the core of the GC. The light-blue dashed line shows the r_cut value, and along the bottom are semitransparent marks showing the exact positions of the stars in the sample.

Zoom In Zoom Out Reset image size

In the left-hand column of Figure 12, the estimated $\overline{{v}^{2}}$ profiles are reasonably well matched to the true profiles for three morphologies (average, extended, and low-Φ₀ GCs). Notably, the corresponding mass profile CMPs in Figure 11 are also some of better estimates of the entire set. For the other two types of GCs, it is the inner part of the profiles that do not match; the true mean-square velocity profile (red curve) in the center of the GC is much higher than the predicted profiles (black curves). Our findings suggest that a reasonable estimate of the $\overline{{v}^{2}}$ profile might be possible for outer regions of the GC when stars are sampled outside the core, but that it would be ill-advised to extrapolate the model fit to the inner regions when only stars outside of the core are available.

In the right-hand column of Figure 12, we see that for every type of GC the true mean-square velocity profile is poorly matched by the predictions at all radii, r. Within the r_cut value, the true profile is generally lower than the black curves, whereas it is much higher than the black curves outside r_cut. Thus, the kinematic information from inner-GC stars alone is not enough to constrain the model at any radii.

One aspect that we have not explored in the biased sampling cases is whether the r_cut value plays a significant role in determining parameter estimates, especially if that r_cut value was more directly linked to GC morphology. Here, we have used a fixed r_cut value mostly for simplicity, but in future work it would be worth exploring the impact of r_cut more fully. For example, the r_cut value for an extended or less-concentrated GC might be relatively smaller than that for the r_cut value for a compact or highly concentrated GC.

It is also worth mentioning that for the fits in the right-hand column of Figure 12, the Markov chains had trouble converging and/or the estimates of the parameter were at the lower or upper ends of the prior distributions (see Table 4). Both of these issues are red flags; the model has not been fit well to the data and any inference would be imprudent.