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:

`retflag`	Description
temp	Include the `tpars` temperature parameters in the retrieval
rad	Include the `rplanet` value as a retrieval parameter
press	Include log `refpressure` value as a retrieval parameter
mol	Include the `molpars` abundance parameters in the retrieval
ray	Include the `rpars` Rayleigh parameters in the retrieval
cloud	Include the `cpars` clouds parameters in the retrieval
patchy	Include the `fpatchy` value as a retrieval parameter
mass	Include the `mplanet` value as a retrieval parameter

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 H₂O 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:

`sampler`	Algorithm	References
snooker	Snooker Differential-Evolution MCMC	[Cubillos2017a]

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")