Monte-Carlo Sampling Higher Dimensional Spaces¶
For one- or two-dimensional parameter studies, a gridded search of parameter space (as in Simple Parameter Study: 2-D Model Grid) is a reasonable approach. However, as the dimensionality grows, things quickly get out of hand. As a result, it can sometimes be advantageous to run Monte Carlo simulations instead, sampling more sparsely (but more efficiently) a high-dimensional space.
You can do this in ARES using the ModelSample
class, which is just a wrapper around the ModelGrid
class. As a result, the problem setup is very similar to that in Simple Parameter Study: 2-D Model Grid, and that structure of the output data are identical, which means the routines documented in Analyzing Model Grids / Monte Carlo Simulations translate as well.
Before we start, the few usual imports:
import ares
import numpy as np
Our “go-to” Efficient Example: \(tanh\) model for the global 21-cm signal¶
To facilitate a comparison between the model grid results, let’s start by choosing the same blobs as in Simple Parameter Study: 2-D Model Grid:
blobs_scalar = ['z_D', 'dTb_D', 'tau_e']
blobs_1d = ['cgm_h_2', 'igm_Tk', 'dTb']
blobs_1d_z = np.arange(5, 21)
base_pars = \
{
'problem_type': 101,
'tanh_model': True,
'blob_names': [blobs_scalar, blobs_1d],
'blob_ivars': [None, [('z', blobs_1d_z)]],
'blob_funcs': None,
}
Now, instead of creating a ModelGrid
instance, we make a ModelSample
instance:
mc = ares.inference.ModelSample(**base_pars)
At this point we have yet to specify which parameters will to sample. Because we are now doing Monte Carlo simulations, we must define the distributions from which to draw samples in each parameter of interest, rather than the grid of values to sample. To do this we need Keith Tauscher’s distpy <https://bitbucket.org/ktausch/distpy> package, which is used in MCMC calculations (e.g., Fitting the Global 21-cm Signal) as well:
from distpy import DistributionSet
ps = DistributionSet()
Now, let’s study the same parameters as Simple Parameter Study: 2-D Model Grid with one addition: the duration of “reheating”:
from distpy import UniformDistribution
# Draw samples from a uniform distribution between supplied (min, max) values for each parameter
ps.add_distribution(UniformDistribution(6, 12), 'tanh_xz0')
ps.add_distribution(UniformDistribution(0.1, 8), 'tanh_xdz')
ps.add_distribution(UniformDistribution(0.1, 8), 'tanh_Tdz')
# Give distributions to the ModelSample instance
mc.prior_set = ps
Note
You can also draw samples from a Gaussian (via GaussianPrior
), a truncated Gaussian (TruncatedGaussianPrior
), and many more. See ares.inference.Priors for a complete listing.
One last thing: we must specify how many random samples to draw:
mc.N = 2e3 # Number of models to run
Finally, to run it:
mc.run('test_3d_mc', clobber=True, save_freq=100)
To analyze the results, create an analysis instance,
anl = ares.analysis.ModelSet('test_3d_mc')
and, for example, plot the 2-d parameter space with points color-coded by tau_e
anl.Scatter(['tanh_xz0', 'tanh_xdz'], c='tau_e', edgecolors='none')
Now that we have also varied the thermal history through tanh_Tdz
, we can look at the interplay between reionization and reheating in setting the emission maximum of the global signal, e.g.,
anl.Scatter(['tanh_xdz', 'tanh_Tdz'], c='dTb_D', edgecolors='none', fig=2)
See Analyzing Model Grids / Monte Carlo Simulations for more information.