BAYESIAN TECHNIQUES FOR COMPARING TIME-DEPENDENT GRMHD SIMULATIONS TO VARIABLE EVENT HORIZON TELESCOPE OBSERVATIONS

We follow the standard Bayesian approach to calculate the posterior likelihood $P({\boldsymbol{w}}| \mathrm{data}$) that a theoretical model described by a vector of parameters $w$ is in agreement with a set of data, i.e.,

$$\begin{eqnarray}&&P({\boldsymbol{w}}| \mathrm{data})={{CP}}_{\mathrm{pri}}({\boldsymbol{w}})P(\mathrm{data}| {\boldsymbol{w}}).\end{eqnarray} \tag{ 1 }$$

Here, P_pri($w$) is the prior likelihood over the model parameters and we calculate the proportionality constant C such that the integral of the posterior likelihood over the parameter space is equal to unity.

The vector of parameters for the GRMHD models that will be compared to EHT observations includes the properties of the black hole spacetime (its mass, spin, and any deviations from the Kerr solution), geometric properties of the orientation of the observer (the orientation and inclination of the spin axis on the observer's sky), and parameters related to the microphysical properties of the plasma (e.g., the ratio of electron-to-ion temperatures in phenomenological models). Moreover, the results of the GRMHD simulations will need to be convolved with a model that describes the image blurring caused by interstellar scattering, which itself will have a number of parameters. Even though we explore below the characteristics of this Bayesian framework using a simplified two-parameter comparison between simulations and data, our methods are general and the full parameter space will need to be considered and explored when fitting actual EHT data.

Hereafter, we will consider a small set for the black hole spin, which has little effect on the overall properties of the black hole image (see Johannsen & Psaltis 2010; Broderick et al. 2011), and fix its inclination with respect to the observer, which can be constrained significantly by the overall size of the 1.3 mm image of the black hole (Psaltis et al. 2015). We will also consider discrete models for the plasma thermodynamics (see Chan et al. 2015b), with parameters chosen such that the simulations reproduce the broadband spectral properties of Sgr A*. This set of choices leaves us with two parameters that need to be obtained by comparing models to interferometric data, i.e., an overall multiplication factor for the flux density, F₀, and the angle ξ between the east–west axis and the projection of the black hole spin angular momentum onto the sky plane, measured in degrees east of north (e.g., Broderick et al. 2011).⁴ Under these assumptions, Bayes' theorem becomes

$$\begin{eqnarray}&&P(\xi ,{F}_{0}| \mathrm{data})={{CP}}_{\mathrm{pri}}(\xi ){P}_{\mathrm{pri}}({F}_{0})P(\mathrm{data}| \xi ,{F}_{0}).\end{eqnarray} \tag{ 2 }$$

We assume flat prior likelihoods over the two model parameters P_pri(ξ) and P_pri(F₀), with −180° ≤ ξ < 180° and F_min ≤ F₀ ≤ F_max, with the range of the overall flux normalization to be specified later.

The last term in Equation (2), $P(\mathrm{data}| \xi ,{F}_{0})$, measures the likelihood that a particular set of data is consistent with a set of model parameters. For the case of EHT, data and simulations need to be compared directly as interferometric observables. A fundamental interferometric observable is the visibility, which is a sample of the Fourier transform of the sky image, I(α, β), at discrete spatial frequencies. The visibility, ${ \mathcal V }(u,v)$, is defined as

$$\begin{eqnarray}&&{ \mathcal V }(u,v)=\iint I(\alpha ,\beta )\exp [-2\pi i(u\alpha +v\beta )]d\alpha d\beta \end{eqnarray} \tag{ 3 }$$

for a projected antenna separation of (u, v) wavelengths within the plane perpendicular to the line of sight.

At any given observing epoch, the EHT will generate a set of simultaneous visibility amplitudes at many baselines as well as sets of closure phases on baseline triangles and closure amplitudes. For the purposes of this initial investigation, and to simplify our notation, we will assume that we only have observations of visibility amplitudes and that each measurement has independent uncertainties from the others. In other words, we will assume that

$$\begin{eqnarray}&&P(\mathrm{data}| \xi ,{F}_{0})\equiv \displaystyle \prod _{i=1}^{M}\displaystyle \prod _{j=1}^{N}{P}_{{ij}}(\mathrm{data}| \xi ,{F}_{0}),\end{eqnarray} \tag{ 4 }$$

where i = 1, 2,..., M and j = 1, 2,..., N denote baselines and epochs, respectively.

Each observation is characterized by a likelihood ${P}^{\mathrm{obs}}({ \mathcal V };{{ \mathcal V }}_{{ij}},{\sigma }_{{ij}})$, centered at ${{ \mathcal V }}_{{ij}}$ and with a dispersion of σ_ij. If, for any pair of model parameters (ξ, F₀), the model made a single prediction for the visibility amplitude for baseline i and epoch j, say ${{ \mathcal V }}_{{ij}}^{\mathrm{sim}}(\xi ,{F}_{0})$, then the single-baseline, single-epoch likelihood in Equation (4) would simply be equal to

$$\begin{eqnarray}&&{P}_{{ij}}(\mathrm{data}| \xi ,{F}_{0})={P}^{\mathrm{obs}}[{{ \mathcal V }}_{{ij}}^{\mathrm{sim}}(\xi ,{F}_{0});{{ \mathcal V }}_{{ij}},{\sigma }_{{ij}}].\end{eqnarray} \tag{ 5 }$$

However, GRMHD models are highly variable (e.g., Chan et al. 2015a) and, for each baseline, they predict a range of visibility amplitudes, with a likelihood that we denote by ${P}_{{ij}}^{\mathrm{sim}}({ \mathcal V };\xi ,{F}_{0})$. In this case, we write the posterior likelihood that a particular observation is consistent with the model predictions as

$$\begin{eqnarray}&&{P}_{{ij}}(\mathrm{data}| \xi ,{F}_{0})=\int {P}^{\mathrm{obs}}({ \mathcal V };{{ \mathcal V }}_{{ij}},{\sigma }_{{ij}}){P}_{{ij}}^{\mathrm{sim}}({ \mathcal V };\xi ,{F}_{0})d{ \mathcal V }.\end{eqnarray} \tag{ 6 }$$

Our goal is to find the set of model parameters (ξ and F₀ in the example here) that maximizes the posterior likelihood in Equation (4), given a set of observations. This Bayesian analysis process, presented in Figure 1 as a flow chart, selects models for which the distribution of predicted visibilities in each baseline matches the distribution of observed visibilites.

Zoom In Zoom Out Reset image size

The first factor in Equation (6), ${P}^{\mathrm{obs}}({ \mathcal V };{{ \mathcal V }}_{{ij}},{\sigma }_{{ij}})$, is well understood. The real and imaginary components of the complex visibility individually follow Gaussian distributions. The visibility amplitude is the magnitude of the complex vector described by those variables, and its posterior likelihood is not Gaussian. The appropriate posterior likelihood for the visibility amplitude, known as a Rice distribution, is

$$\begin{eqnarray}&&{P}^{\mathrm{obs}}({ \mathcal V };{{ \mathcal V }}_{{ij}},{\sigma }_{{ij}})=\displaystyle \frac{{ \mathcal V }}{{\sigma }_{{ij}}^{2}}\exp \left[-\displaystyle \frac{({{ \mathcal V }}^{2}+{{ \mathcal V }}_{{ij}}^{2})}{2{\sigma }_{{ij}}^{2}}\right]{I}_{0}\left(\displaystyle \frac{{ \mathcal V }{{ \mathcal V }}_{{ij}}}{{\sigma }_{{ij}}^{2}}\right),\end{eqnarray} \tag{ 7 }$$

where I₀ is the zeroth order modified Bessel function of the first kind. This likelihood approaches a Gaussian for large signal-to-noise ratios (see Chapter 6 of Thompson et al. 2001).

In order to simplify the calculation of integral (6), we make use of the fact that the distribution of simulated visibility amplitudes ${P}_{{ij}}^{\mathrm{sim}}({ \mathcal V };\xi ,{F}_{0})$ has discrete values rather than a continuous distribution, as it is obtained from snapshot images. For T snapshots,

$$\begin{eqnarray}&&{P}_{{ij}}^{\mathrm{sim}}({ \mathcal V };\xi ,{F}_{0})=\displaystyle \frac{1}{T}\displaystyle \sum _{t=1}^{T}\delta [{ \mathcal V }-{{ \mathcal V }}_{{ij}}^{\mathrm{sim}}(t;\xi ,{F}_{0})],\end{eqnarray} \tag{ 8 }$$

where ${{ \mathcal V }}_{{ij}}^{\mathrm{sim}}(t;\xi ,{F}_{0})$ is the visibility amplitude sampled at (u_ij, v_ij) coordinates from the tth snapshot with the given set of parameters (ξ, F₀). Therefore, Equation (6) simplifies to

$$\begin{eqnarray}&&{P}_{{ij}}(\mathrm{data}| \xi ,{F}_{0})\\ &&\quad =\displaystyle \frac{1}{T}\displaystyle \sum _{t=1}^{T}\displaystyle \frac{{{ \mathcal V }}_{{ij}}^{\mathrm{sim}}(t;\xi ,{F}_{0})}{{{\sigma }_{{ij}}}^{2}}\exp \left[-\displaystyle \frac{{({{ \mathcal V }}_{{ij}}^{\mathrm{sim}}(t;\xi ,{F}_{0})+{{ \mathcal V }}_{{ij}})}^{2}}{2{{\sigma }_{{ij}}}^{2}}\right]\\ &&\qquad \times {I}_{0}\left(\displaystyle \frac{{{ \mathcal V }}_{{ij}}^{\mathrm{sim}}(t;\xi ,{F}_{0}){{ \mathcal V }}_{{ij}}}{{{\sigma }_{{ij}}}^{2}}\right)\end{eqnarray} \tag{ 9 }$$

for the likelihood ${P}_{{ij}}^{\mathrm{sim}}({ \mathcal V };\xi ,{F}_{0})$, with given parameters.

In the following set of examples, we calculate the posterior likelihoods over the model parameters using Equations (2), (4), and (9). In real applications, however, observations on different baselines are performed simultaneously, so we must consider a joint probability distribution in visibility space. For example, when an observation is carried out with three stations (or equivalently, three baselines), the posterior likelihood for the observations, ${P}^{\mathrm{obs}}({ \mathcal V })$, becomes a multivariate Rician. Moreover, we need to consider the posterior likelihood, ${P}^{\mathrm{sim}}({ \mathcal V })$, that the model predicts a particular combination of simultaneous visibilities on the same baselines.

Finally, we can trivially incorporate to our Bayesian inference other interferometric observables, such as closure phases and closure amplitudes. The closure phase is the sum of visibility phases around the triangle of baselines formed by any set of three stations. The closure phase is independent of instrumental and atmospheric phase fluctuations and depends only on the visibility phase of the source. If we have ${ \mathcal N }$ stations where ${ \mathcal N }\geqslant 3$, ${}_{{ \mathcal N }-1}{C}_{2}=({ \mathcal N }-1)({ \mathcal N }-2)/2$ independent closure phases can be measured. The posterior likelihood for closure phase data becomes

$$\begin{eqnarray}&&{P}_{{lj}}(\mathrm{data}| \xi ,{F}_{0})=\int {P}^{\mathrm{obs}}({\rm{\Phi }};{{\rm{\Phi }}}_{{lj}},{\sigma }_{{lj}}){P}_{{lj}}^{\mathrm{sim}}({\rm{\Phi }};\xi ,{F}_{0})d{\rm{\Phi }}\end{eqnarray} \tag{ 10 }$$

for $l=1,2,\ldots {,}_{{ \mathcal N }-1}{C}_{2}$ sets of telescopes (closure phase "triangles"), and j = 1, 2,..., N observations. Here, P^obs(Φ; Φ_lj, σ_lj) and P^sim_lj (Φ; ξ, F₀) are the likelihoods for closure phases from the observation and simulations, similar to ${P}^{\mathrm{obs}}({ \mathcal V };{{ \mathcal V }}_{{ij}},{\sigma }_{{ij}})$ and ${P}_{{ij}}^{\mathrm{sim}}({ \mathcal V };\xi ,{F}_{0})$ in Equation (6).

In addition to inferring the most likely parameters of a given model from observations, we are also interested in comparing the ability of different models to describe the observed distributions of visibility amplitudes and closure quantities.

Information criteria approaches such as the Akaike Information Criterion and the Bayesian Information Criterion (BIC) have been adopted elsewhere to evaluate the quality of astrophysical models (e.g., Liddle 2007). Broderick et al. (2011) applied these criteria in evaluating the relative success of time-independent models of Sgr A* in fitting EHT data. The BIC is defined as

$$\begin{eqnarray}&&\mathrm{BIC}\equiv -2\mathrm{log}{{ \mathcal L }}_{\max }+k\mathrm{log}N,\end{eqnarray} \tag{ 11 }$$

where ${{ \mathcal L }}_{\max }$ is the maximum likelihood using the parameters within the model, k is the number of free parameters, and N is the number of data points. The definition implies that the smaller BIC value suggests a better fit to the data. However, criteria such as the BIC consider only the maximum posterior likelihoods and do not account for the relative extent of the volumes in parameter space within which each model is consistent with the data.

In the context of our work, we will perform model comparisons in terms of an appropriately defined Bayesian evidence. If we consider each model as a single point in a discrete parameter space, then, given the data, we can use Bayes' theorem to calculate the posterior likelihood for each model as

$$\begin{eqnarray}&&P({{ \mathcal M }}_{m}| \mathrm{data})={{CP}}_{\mathrm{pri}}({{ \mathcal M }}_{m})P(\mathrm{data}| {{ \mathcal M }}_{m}).\end{eqnarray} \tag{ 12 }$$

Here m is the discrete index specifying each model {${{ \mathcal M }}_{m}$}, which is governed by a set of parameters ${{\boldsymbol{w}}}_{m}$. If we denote by P_pri(m) the prior likelihood for each model and by P_pri,m(${{\boldsymbol{w}}}_{m}$) the prior likelihoods for each set of parameters within model m, then we can write

$$\begin{eqnarray}&&P({{ \mathcal M }}_{m}| \mathrm{data})={{CP}}_{\mathrm{pri}}(m){P}_{\mathrm{pri},m}({{\boldsymbol{w}}}_{m}){P}_{m}(\mathrm{data}| {{\boldsymbol{w}}}_{m}).\end{eqnarray} \tag{ 13 }$$

The last term in Equation (13) is the posterior Bayesian likelihood of the data, given a model m and a set of model parameters, e.g., Equation (4).

In order to calculate the Bayesian evidence for model m, we marginalize Equation (13) over the space of model parameters, i.e.,

$$\begin{eqnarray}&&{ \mathcal L }(m| \mathrm{data})={{CP}}_{\mathrm{pri}}(m)\int {P}_{\mathrm{pri},m}({{\boldsymbol{w}}}_{m}){P}_{m}(\mathrm{data}| {{\boldsymbol{w}}}_{m})d{{\boldsymbol{w}}}_{m}.\end{eqnarray} \tag{ 14 }$$

Assuming a flat prior of models and the same prior likelihood over all model parameters for each model, this intergral simplifies to

$$\begin{eqnarray}&&{ \mathcal L }(m| \mathrm{data})={C}^{\prime }\int {P}_{m}(\mathrm{data}| {{\boldsymbol{w}}}_{m})d{{\boldsymbol{w}}}_{m},\end{eqnarray} \tag{ 15 }$$

where C' is an appropriate normalization constant.

The integrated likelihood ${ \mathcal L }(m| \mathrm{data})$ provides a measure of model selection statistics. As is the standard practice in Bayesian inference, we use the Bayesian evidence only for comparison among models and not to assign a posterior likelihood to any model individually.

In a recent study, Chan et al. (2015a) investigated the variability and millimeter/infrared flare characteristics of GRMHD simulations. The simulation of the accretion flow was performed with the three-dimensional GRMHD code HARM (Gammie et al. 2003; McKinney 2006; McKinney & Blandford 2009; Narayan et al. 2012; Saḑowski et al. 2013) and the ray-tracing of photon trajectories was performed with the fast GPU-based algorithm GRay (Chan et al. 2013). The simulated Sgr A* images with horizon-scale resolution and their spectra were fit to both the observed broadband spectra and the overall 1.3 mm source size determined by EHT observations (Chan et al. 2015b).

Chan et al. (2015a, 2015b) simulated high cadence images of the accretion flow around Sgr A* for each model. The simulation results were recorded with a timestep of 10 GMc⁻³, which corresponds to 212 s for the mass of Sgr A*. Each snapshot contains four levels of field of view: 16M × 16M, 64M × 64M, 256M × 256M, and 1024M × 1024M, containing 512 × 512 pixels. Here 1 GMc⁻² corresponds to 5.1 μas for the distance to Sgr A*. We use a composite of the multi-level resolutions for our analysis in order to encompass all the emission near and away from the black hole and to generate interferometric visibilities. It is reasonable to limit our images to 1024M × 1024M, because the corresponding 10.2 μas pixel scale is considerably finer than the highest spatial frequency (58.9 μas fringe spacing) sampled by the EHT data we use as a comparison.

The GRMHD models of Chan et al. (2015a, 2015b, see also Narayan et al. 2012; Saḑowski et al. 2013) were categorized into two classes based on initial magnetic field configurations: Standard And Normal Evolution (SANE, disk-dominated) and Magnetically Arrested Disk (MAD, jet-dominated) models, which use multi-loop and single-loop initial magnetic fields, respectively. The simulation images at 1.3 mm for these two classes of models show dissimilar features from each other: in SANE models, the emission originates from a crescent-like shape in the disk region, whereas it is concentrated at the footprints of the outflowing material in the MAD models. Chan et al. (2015a, 2015b) also explored the dependence of the predictions of the GRMHD models on the black hole spin parameter a and on the way in which the electron temperature was prescribed.

In some of their models, the electron temperature in the funnel region of the accretion flow was fixed to be constant (constant T_e,funnel). In other models, the ratio of the electron and ion temperatures in the funnel was fixed (constant θ_funnel). In both cases, the electron and ion temperatures have a constant ratio in the disk. Among all the configurations they explored, they narrowed the possibilities down to five models that account for both the broadband spectrum of Sgr A* and its overall image size at 1.3 mm:

• (1)  
  Model A: a = 0.7, SANE, constant T_e,funnel.
• (2)  
  Model B: a = 0.9, SANE, constant T_e,funnel.
• (3)  
  Model C: a = 0.0, MAD, constant θ_funnel.
• (4)  
  Model D: a = 0.9, MAD, constant T_e,funnel.
• (5)  
  Model E: a = 0.9, MAD, constant θ_funnel.

All models show variability over the entire length of the simulation, i.e., over ∼10,000 GMc⁻³ (∼60 hrs). Figure 2 shows the 1.3 mm light curves for the simulations. Overall, the SANE models are observed to produce more variability than the MAD models.

Zoom In Zoom Out Reset image size

We generated visibilities from the simulations using our own VLBI simulation code. The algorithm calculates (u, v) coordinates from the observing time and celestial coordinates of the source, performs Fourier transform of the input image, and samples complex visibilities according to the telescope array setup. We tested the results of this code against the MIT Array Performance Simulator⁵ and found excellent agreement between the simulated visibilities.

Figure 3 shows the range of predicted visibility amplitudes for the five GRMHD models. The SANE models show more variability in the visibility space than the MAD models, just as they do in total intensity (Figure 2). This variability is explored in more detail in the companion paper, Medeiros et al. (2016).

Zoom In Zoom Out Reset image size

In the following sections, we describe tests of our Bayesian analysis using mock data from these simulations. These tests demonstrate that our statistical method successfully recovers model parameters (Section 3.2) and can faithfully distinguish between models (Section 3.5).

In our first test we investigated the performance of our statistical method in inferring model parameters. We constructed the test as follows:

• (1)  
  We assumed a fiducial set of model parameters (ξ, F₀) = (−30°, 1).
• (2)  
  We generated data using three baselines with three stations: Hawaii—Arizona, Hawaii—California, and Arizona—California.
• (3)  
  We sampled (u, v) coordinates during the period when three stations can perform simultaneous observations of Sgr A*.
• (4)  
  We generated visibilities using model D, which is an MAD model.
• (5)  
  We randomly selected visibility amplitudes from the distribution of simulated visibilities, ${P}_{{ij}}^{\mathrm{sim}}({ \mathcal V };-30^\circ ,1)$, for each baseline. Visibilities were taken from the same simulation frame for all three baselines. We then added random Gaussian noise to the complex visibilities to produce mock observational data from the visibility amplitude samples.
• (6)  
  We computed the posterior likelihood $P(\mathrm{data}| \xi ,{F}_{0})$ using simultaneous observations on three baselines, as described in Equation (4). We considered orientation ξ in the range from −90° to 90°, as visibility amplitudes without phase information are degenerate under rotations of 180°. We repeated the calculation for 5000 different realizations of the mock data.

We averaged the posterior likelihoods to form 500 sets of data, where each set contains N = 10 observations. For every set of data, we found the orientation and flux scale where the maximum likelihood occurs. The distributions of these maximum likelihood parameters for our input simulation are shown in Figures 4(a) and (b). The red vertical lines mark the model parameters (ξ, F₀) that were assumed when generating the mock visibility data. The histograms indicate that the orientation and flux scale of the highest likelihoods of all data sets are distributed around the assumed parameters. We do not see any bias in the inferred parameters.

Zoom In Zoom Out Reset image size

Figure 4(c) shows the averaged posterior likelihoods of 500 sets of data. The red dotted lines in the parameter space again mark the combination of the assumed model parameters (ξ, F₀) = (−30°, 1). The lines cross very near the peak of the posterior likelihood, and their intersection is located inside the 68% credible region of the likelihood. Figure 4(d) shows the averaged likelihoods when the 5000 mock visibilities are grouped into 200 sets, each including N = 25 observations. Compiling a larger number of observations causes the allowed region of the parameter space to contract, as expected.

In the near future, more stations will be incorporated into the EHT, such as the Atacama Large Millimeter/submillimeter Array (ALMA) and the South Pole Telescope (SPT). We performed a mock observation test to explore the effect of future EHT stations. The test is similar to that of the previous section, but assuming we have 15 baselines with 6 stations. The stations are the Submillimeter Array (SMA) in Hawaii, the Combined Array for Research in Millimeter-wave Astronomy (CARMA) in California, the Submillimeter Telescope (SMT) in Arizona, the Large Millimeter Telescope (LMT) in Mexico, ALMA in Chile, and the SPT at the South Pole. The sampled (u, v) coordinate sets include non-simultaneous observations because considering only simultaneous detections for six stations greatly limits the observing period. Visibility amplitude errors are estimated from the baseline sensitivity, which is proportional to the geometric mean of the System Equivalent Flux Densities (SEFDs) of antennas.⁶ We assumed 10 s coherent integration time and 2 GHz bandwidth for the sensitivity.

There are three factors that affect the allowed parameter ranges in the Bayesian analysis: the overall sensitivity of the interferometric array, (u, v) coverage in the visibility space, and sampling of the intrinsic variability of the source. Figure 4(e) shows the combined effect of these: the increased sensitivity and (u, v) coverage of a six-station EHT will narrow the allowed parameter space compared to Figure 4(d), which has the same number of observations but fewer stations and lower sensitivity. To assess the importance of source variability in parameter estimation, we constructed three array setups as follows:

• (a)  
  The early EHT array of three stations at Hawaii, Arizona, and California, with 920 observing epochs.
• (b)  
  The full EHT array with six stations when the sensitivity of ALMA is replaced with that of the co-located Atacama Pathfinder Experiment (APEX) telescope,⁷ with 61 observing epochs.
• (c)  
  The full EHT array with six stations, including ALMA, with four observing epochs.

In these three scenarios, we have varied the number of observing epochs such that the accumulated point source sensitivity of the array is constant. That is, the total thermal noise is the same in all the arrays, after accounting for the varying number of observing epochs. The comparison of the six-station array with APEX is in interesting counterpoint to the case with ALMA because it holds the (u, v) coverage fixed (and the sensitivity), but varies the number of observing epochs. Figure 5 shows that the parameter constraints are much tighter when we average over more observing epochs, emphasizing the importance of sampling the intrinsic source variability when attempting to constrain the model parameters with EHT observations. Figure 5(a) has the poorest (u, v) coverage by far, but achieves the tightest constraint, while fixing both (u, v) coverage and sensitivity in Figures 5(b) and (c) again show that simply sampling more epochs will tighten the parameter constraints significantly.

Zoom In Zoom Out Reset image size

In the parameter estimation test, we calculated likelihoods over parameter space with both time-dependent and averaged GRMHD images to compare the effect of time averaging. Figure 6 contrasts the likelihoods of time-resolved (black contours) and time-independent (red contours) inputs to our Bayesian method when the same number of mock observations are analyzed. Figures 6(a) and (b) are the likelihoods assuming the early three-station baselines and the six-station EHT with ALMA, respectively. The figure shows that a much smaller region of parameter space is allowed when comparing the mock data to a time-averaged image. However, the best-fit parameters in the time-averaged case are offset from the true parameters, as is most clearly evident in Figure 6(b). This demonstrates the false certainty implied by the small allowed parameter region. The time-resolved analysis, by contrast, delivers allowed parameter ranges that are consistent with the input simulation.

Zoom In Zoom Out Reset image size

We now use the Bayesian method to identify the relative posterior likelihood among several models, as discussed in Section 2.2. For this test we generated mock visibilities following the same procedure as the test of the previous section. We used (u, v) coordinates identical to those in the EHT observations of 2007 and 2009. As discussed earlier, SANE models present greater variability than MAD models (Figure 3). Therefore, we generated two mock data sets, one from model A (SANE) and one from model E (MAD), and performed our model selection test twice (labeled test A and test E, respectively).

We present in Table 1 the marginal likelihoods in Equation (15) and the $\bigtriangleup \mathrm{BIC}$ values in Equation (11). We obtained these values by fitting the mock model A data and mock model E data with all five models. The ratio of the marginal likelihoods for two models is the Bayes factor, and it measures the relative strength of the evidence against another one. Jeffreys's scale provides guidance for the interpretation of the Bayes factor (Kass & Raftery 1995). A ratio between 10 and 100 is translated as strong, and a ratio greater than 100 is regarded as decisive evidence. A similar approach using the difference of information criteria exists and $\bigtriangleup \mathrm{IC}$ greater than 5 and 10 are understood as strong and decisive, respectively.

Table 1.  Marginal Likelihoods of Model Estimation Tests

	Test A		Test E
	${ \mathcal L }(m| \mathrm{data})$	△BIC	${ \mathcal L }(m| \mathrm{data})$	△BIC
Model A	1.00	0.0	5.36 × 10−7	28.5
Model B	0.799	1.5	5.31 × 10−10	42.8
Model C	7.20 × 10−9	37.5	9.90 × 10−4	12.7
Model D	1.81 × 10−4	17.3	3.47 × 10−5	20.8
Model E	2.98 × 10−10	41.4	1.00	0.0

Note. Test A uses mock data using model A, and test E uses mock data using model E. The marginalized likelihoods, ${ \mathcal L }(m| \mathrm{data})$, are normalized to the best model for both tests, as are the BIC values.

Download table as:  ASCIITypeset image

In each test, our analysis prefers the model from which the data were generated. In Test A, model B is only slightly disfavored compared to model A. This results from the similar shapes of the visibility amplitude distributions in the SANE simulations. The fact that model B has a different spin than model A is also an indication that black hole spin has a very small effect on the overall appearance of SANE models for Sgr A*. In this test, the MAD models are decisively rejected, as expected, because the smaller variability of the model visibilities in the MAD simulations is inconsistent with the mock data. Conversly, Test E decisively favors model E but rejects all four other models.