MCMC Tutorial
This mode allows you to fit spectra to observed exoplanet data.
Pyrat Bay
uses the MC3
package (https://github.com/pcubillos/mc3) to retrieve best-fitting
parameters and credible regions for the atmospheric parameters in a
Bayesian (MCMC) framework.
Sample Configuration File
Here is an example of an opacity-table configuration file (mcmc_eclipse.cfg):
[pyrat]
# Pyrat Bay run mode, select from: [tli atmosphere spectrum opacity mcmc]
runmode = mcmc
# Radiative-transer observing geometry, select from: [transit emission]
rt_path = emission
# Atmospheric model:
atmfile = WASP-00b.atm
# Opacity file:
extfile = ./exttable_100-3000K_1.0-5.0um.npz
# Cross-section opacity files:
csfile = {ROOT}/pyratbay/data/CIA/CIA_Borysow_H2H2_0060-7000K_0.6-500um.dat
# Spectrum boundaries and sampling rate:
wllow = 1.0 um
wlhigh = 5.0 um
wnstep = 1.0
wnosamp = 2160
# System parameters:
rstar = 1.27 rsun
tstar = 5800.0
gstar = 22908.0
smaxis = 0.045 au
mplanet = 0.6 mjup
rplanet = 1.0 rjup
tint = 100.0
# Reference pressure level at rplanet:
refpressure = 0.1 bar
# Filter bandpasses:
filters =
filters/tutorial_band00.dat
filters/tutorial_band01.dat
filters/tutorial_band02.dat
filters/tutorial_band03.dat
filters/tutorial_band04.dat
filters/tutorial_band05.dat
filters/tutorial_band06.dat
filters/tutorial_band07.dat
filters/tutorial_band08.dat
filters/tutorial_band09.dat
# Eclipse-depth data:
data =
0.0003296 0.0003445 0.0003014 0.0002637 0.0003159
0.0004956 0.0006660 0.0007359 0.0009185 0.0009853
uncert =
0.0000204 0.0000219 0.0000235 0.0000250 0.0000264
0.0000280 0.0000295 0.0000311 0.0000327 0.0000342
# Kurucz stellar model:
kurucz = ./fp00k2odfnew.pck
# Abundance models:
molmodel = vert
molfree = H2O
bulk = H2 He
# Radius-profile model, select from: [hydro_m hydro_g]
radmodel = hydro_m
# Temperature-profile model, select from [isothermal guillot madhu]
tmodel = guillot
# Rayleigh models, select from: [lecavelier dalgarno_H dalgarno_He dalgarno_H2]
rayleigh = lecavelier
rpars = 0.0 -4.0
# Retrieval models, select from: [temp rad mol ray cloud patchy mass]
retflag = temp mol
# Fitting parameters:
# log(k') log(g1) log(g2) alpha Tirr Tint log(fH2O)
params = -3.0 -1.0 0.0 0.0 1500 0 -3.0
pmin = -6.0 -3.0 -3.0 0.0 10 0 -9.0
pmax = 1.0 3.0 3.0 1.0 3000 100 0.0
pstep = 0.3 0.3 0.0 0.0 10 0 0.5
# MCMC temperature boundaries (K):
tlow = 100
thigh = 3000
# MCMC parameters:
sampler = snooker
nsamples = 100000
burnin = 3000
nchains = 14
ncpu = 7
# Log-scale X-axis:
logxticks = 1.0 1.4 2.0 2.3 3.0 4.0 5.0
# Maximum optical depth to calculate:
maxdepth = 10.0
# Output file names:
mcmcfile = MCMC_emission_tutorial.npz
# Verbosity level (<0:errors, 0:warnings, 1:headlines, 2:details, 3:debug):
verb = 2
Note
Note that an ‘mcmc’ run requires the user to define an
opacity table (extfile
) to allow the code to finish
within a Hubble time.
Observed Data
A retrieval run must include the data
, uncert
, and filters
keys to define the observed data that will guide the retrieval
sampler.
The data
key sets the observed transit depths for transmission
observations, or eclipse depths for emission observations. The
uncert
key sets the data uncertainties, and the filters
key
set the path to the filter pass-band files corresponding to each data
value. Use the dunits
key to specify the units of the data
and
uncert
values (default: dunits = none
).
For eclipse geometry, the user needs to input a stellar flux model to compute eclipse depths (see Stellar Spectrum).
Retrieval Models and Parameters
Use the retflag
key to set which model parameters will vary in the
retrieval. The following table list the available options:
|
Description |
---|---|
temp |
Include the |
rad |
Include the |
press |
Include log |
mol |
Include the |
ray |
Include the |
cloud |
Include the |
patchy |
Include the |
mass |
Include the |
Use the params
key to set initial values for the retrieval
parameters. Note that the user need to set these parameters in the
same order as listed above, i.e.: first the temperature parameters,
then radius, then the mass, the abundance parameters, and so on.
Details on the available models and their parameters are described in
the Atmosphere Tutorial and Spectrum Tutorial. The number of
rad
, press
, and mass
parameters is always one each.
The units for log refpressure
are always bars.
The units for rplanet
and mplanet
are set by the runits
and mpunits
keys, or else, will adopted from the units
specified in the rplanet
and mplanet
inputs.
Note
Pro-tip: You can use the params
key to run spectrum
forward models as well.
Use the pmin
and pmax
keys to set minimum and maximum
retrieval boundaries for each parameter listed in params
,
respectively.
Use the pstep
key to handle the behavior of the retrieval
parameters. A pstep
value of zero keeps a parameter fixed during
the retrieval. A negative-integer pstep
value indicates that such
parameter is shared with the indexed value (i.e., a parameter with a
pstep
value of -2 will take the value from the second parameter).
Parameter Priors
Currently, Pyrat Bay
enables uniform (default) or Gaussian priors
for the retrieval parameters. Use the prior
, priorlow
, and
priorup
keys to set the mean, upper- and lower-standard deviation
value of the Gaussian priors, respectively.
If both the priorlow
and priorup
values are zero for a given
parameter, the prior will remain uniform between the pmin
and
pmax
values.
For example, for the configuration file above, the following config sets uniform priors for the temperature parameters, and a Gaussian prior for the H2O abundance parameter of \(\log({\rm H2O}) = -3.5 \pm 0.3\):
# log(kappa) log(g1) log(g2) alpha beta log(H2O)
prior = 0.0 0.0 0.0 0.0 0.0 -3.5
priorlow = 0.0 0.0 0.0 0.0 0.0 0.3
priorup = 0.0 0.0 0.0 0.0 0.0 0.3
Retrieval Setup
Pyrat Bay
enables posterior sampling via the Markov-chain Monte
Carlo (MCMC) or Nested Sampling technique. Use the sampler
key to set
the sampling algorithm. The following table list the available
options and references of their implementation:
|
Algorithm |
References |
---|---|---|
snooker |
Snooker Differential-Evolution MCMC |
The ‘snooker’ option implements the DEMC-z algorithm with snooker proposals, described in [terBraak2008].
MCMC Retrieval
The nsamples
key sets the total number of MCMC samples.
The nchains
key sets the number of parallel MCMC chains to use.
The burnin
key sets the number of iterations to remove (burn) at
the beginning of each chain.
Examples
Since this is an eclipse retrieval, the code requires a stellar spectrum model to compute the planet-to-star flux ratios. You will also to download the configuration file shown above and the observing filter files, e.g., with these shell commands:
# Download Kurucz stellar model:
wget http://kurucz.harvard.edu/grids/gridp00odfnew/fp00k2odfnew.pck
tutorial_path=https://raw.githubusercontent.com/pcubillos/pyratbay/master/examples/tutorial
# Download the configuration file:
wget $tutorial_path/mcmc_eclipse.cfg
# Download the filter files:
for i in {0..9}
do
wget $tutorial_path/filters/tutorial_band0${i}.dat
done
mkdir filters/
mv filters/tutorial_band0?.dat filters/
Also, be sure to have create the opacity file as described in Opacity-grid Tutorial.
As in a spectrum run, an MCMC run returns a ‘pyrat’ object in an interactive run. The following Python script computes an opacity file using the configuration file found at the top of this tutorial. Just like before, to run the MCMC modeling, simply execute this command:
import pyratbay as pb
pyrat = pb.run("mcmc_eclipse.cfg")