Classes and functions

beansp.beans.prior_func(theta_in)

This function is the default prior and implements a simple box prior for all the parameters

Parameters:

theta_in – parameter vector, with X, Z, Q_b, d, xi_b, xi_p, and (optionally) mass, radius, f_E & f_a

Returns:

prior probability

beansp.beans.prior_1808(theta_in)

This function implements a simple box prior for all the parameters excluding mass, radius, and the metallicity, which come instead from :function:`beansp.beans.lnZprior`

This prior is explicitly intended for use with SAX J1808.4-3658, and should only be used with extreme caution in other cases!

Parameters:

theta_in – parameter vector, with X, Z, Q_b, d, xi_b, xi_p, mass, radius, and (optionally) f_E & f_a

Returns:

prior probability

beansp.beans.lnZprior(z)

This beta function for the metallicity prior is from Andy Casey and is an approximation of the metallicity of a mock galaxy at 2.5-4.5 kpc for the location of SAX J1808.4-3658. Assuming ZCNO = 0.01 is average value.

Parameters:

z – CNO metallicity (mass fraction)

Returns:

prior probability

beansp.beans.corr_goodwin19(burst, **kwargs)

This is an example Settle correction function that applies the correction of Goodwin et al. (2019), multiplying the recurrence time and burst energy by 0.65. With this correction there is no need to modify the alpha value. The **kwargs can be used to incorporate other parameters into the correction; from settle we pass the base flux F, mdot M, H-fraction X, metallicity Z, and neutron star radius R and mass M (none used for this example)

Parameters:

burst – 3-element tuple output from settle, with alpha, trec [hr], and burst energy [1e39 erg], all values in the observer frame

Returns:

tuple with any necessary corrections performed

class beansp.Beans(prior=<function prior_func>, corr=None, config_file=None, run_id='test', nwalkers=200, nsteps=100, obsname=None, burstname=None, gtiname=None, interp='linear', smooth=0.02, theta=(0.58, 0.013, 0.4, 3.5, 1.0, 1.0, 1.5, 11.8), stretch_a=2.0, fluen=True, alpha=True, numburstssim=3, bc=2.21, ref_ind=1, threads=4, test_model=True, restart=False, **kwargs)

The main object class that includes the basic functionality required for beans. The code will read in burst (and observation) data and attempt to simulate bursts to match the observed burst properties. There are two principle modes; the original function generates a “train” of individual bursts, observed (for example) during a transient outburst, as for the original application to the 2002 outburst of SAX J1808.4-3658, observed with RXTE/PCA. The alternative is to match to a set of non-contiguous bursts (“ensemble” mode)

__init__(prior=<function prior_func>, corr=None, config_file=None, run_id='test', nwalkers=200, nsteps=100, obsname=None, burstname=None, gtiname=None, interp='linear', smooth=0.02, theta=(0.58, 0.013, 0.4, 3.5, 1.0, 1.0, 1.5, 11.8), stretch_a=2.0, fluen=True, alpha=True, numburstssim=3, bc=2.21, ref_ind=1, threads=4, test_model=True, restart=False, **kwargs)

Initialise a Beans object

Parameters:

• prior – prior function to use

• corr – correction function for bursts, or None

• config_file – file to read in configuration parameters (in which case the keyword params below are ignored)

• run_id – string identifier for the run, used to label all the result files, and where you want output to be saved

• nwalkers – number of walkers for the emcee run

• nsteps – number of MCMC steps to run

• obsname – name of the file including the flux history, from which the mdot is estimated to generate to generate the burst train (set obsname=None for a non-contiguous, or “ensemble” mode run)

• burstname – name of the burst data file, listing the bursts

• gtiname – name of the GTI file, set to None to turn off checking

• interp – interpolation mode for the flux; possible values are ‘linear’, or ‘spline’

• smooth – smoothing factor for spline interpolation

• theta – initial centroid values for walker model parameters, with X, Z, Q_b, d, xi_b, xi_p, and (optionally) mass, radius, f_E & f_a

• stretch_a – the Goodman & Weare (2010) stretch move scale parameter, passed to emcee

• fluen – set to True (default) to include the fluences in the data for comparison, or False to omit

• alpha – set to True (default) to include the alphas in the data for comparison, or False to omit

• numburstssim – number of bursts to simulate, for the “train” mode, both earlier and later than the reference burst; i.e. set to half the total number of bursts you want to simulate. Don’t forget to account for missed bursts!

• bc – bolometric correction to adopt for the flux history (set to 1.0 if the fluxes are already bolometric):

• ref_ind – rank of “index” burst, against which the other burst times are relative to. For the “train” mode, should be around the middle of the predicted burst train. This burst will not be simulated but will be used as a reference to predict the other bursts.

• threads – number of threads for emcee to use (e.g. number of cores your computer has). Set to None to use all available

• test_model – flag to test the model during the setup process

• restart – set to True to continue a previously interrupted run

Result:

Beans object including all the required data

burst_table(show=True, predicted=False, key=None)

This method creates and, optionally, displays a table from the model and observed bursts.

The table is returned as an astropy Table object that can be saved as an MRT file, like so: tab.write('test.dat',format='mrt') although you will have to then edit the file afterwards to complete the metadata

Parameters:

• show – set to True to display each row (in LaTeX format)

• predicted – set to True to show predicted bursts also (not yet implemented)

• key – for the case where we have more than one set of predictions (e.g. for different numbers of predicted bursts), the key will identify which to return

Returns:

astropy table of burst properties

do_analysis(options=['autocor', 'posteriors'], part=None, truths=None, burnin=2000, savefig=False)

This method is for running standard analysis and displaying the results. Nothing is returned, but various options for analysis are available, which will (optionally) populate some of the beansp.Beans attributes, including

reader - sampler object read in from HDF file

sampler - chain data

nsteps_completed - number of steps completed

samples - flattened samples, omitting first samples_burnin

last - last walker positions

probs - likelihoods (total, prior, and contributions from each parameter) for the last walker positions

cc - ChainConsumer object with samples and derived cosi, gravity, redshift

cc_parameters - plot labels for cc object

model_pred - dictionary with model realisations read in from the “blobs”

By default the method will also create several files, labeled by the run_id; drawn from

{}_autocorrelationtimes.pdf (via plot_autocorr)

{}_posteriors.pdf

{}_parameterconstraints_pred.txt

TODO: need to reorganise a bit, and add more options

Parameters:

• options – array of strings corresponding to various analysis options, listed in the analyses dict below

• part – string or array “partition” dividing the set of samples into two or more separate groups for analysis

• truths – parameter vector to overplot on (one of the) corner plots (TODO: need to check if >1 corner plot are selected +truths, which will likely result in an error due to incompatible number of parameters)

• burnin – number of steps to discard when plotting the posteriors

• savefig – set to True to save figures to .pdf files, False to skip

Returns:

none

do_run(plot=False, analyse=True, burnin=2000, **kwargs)

This routine runs the chain for as many steps as is specified in the init call. Emcee parameters are defined in runemcee module. Have previously used multiprocessing to speed things up, but that’s not currently active

Parameters:

• plot – set to True to do the plot of the initial guess. This seems redundant since it’s also plotted at the __init__ stage

• analyse – set to True to call do_analysis automatically once the chains finish

• burnin – number of steps to ignore at the start of the run, passed to do_analysis

Returns:

lnprob(theta_in, x, y, yerr)

The full log-probability function incorporating the priors (via the Beans.lnprior attribute), and model likelihood (via Beans.lnlike()), that is passed to runemcee when creating the sampler (in the Beans.do_run() method).

Parameters:

• theta_in – parameter vector, with X, Z, Q_b, d, xi_b, xi_p, mass, radius, and (optionally) f_E & f_a

• x – the “independent” variable, passed to Beans.lnlike()

• y – the “dependent” variable (i.e. measurements), passed to Beans.lnlike()

• yerr – erorr estimates on y

Returns:

total (prior+model) likelihood, prior likelihood, model array (from Beans.lnlike())

plot(show_model=True, model=None, mdot=True, title=None, savefig=False)

Display a plot of the data and model results, for a burst train calculated with generate_burst_train Adapted from the example at https://matplotlib.org/gallery/api/two_scales.html

Parameters:

• show_model – set to False to skip the model generation step, in which case only the observed data will be plotted

• model – array of packed model prediction, OR dict giving full model results

• mdot – flag to show mdot rather than flux (only possible if you’re passing the full model)

• title – add a title, if required

• savefig – set to True to save the figure to PDF

sim_data(file=None)

This method generates a table of simulated data based on the current theta vector, that can be used for tests of the MCMC method

Errors are taken as the mean of the observational errors; bursts matching observed bursts (according to the matching algorithm) are flagged as True in the 6th column

In lieu of an output file, you can just paste the results into a text file and use as the input (specify as burstname)

Parameters:

file – file to save the results to (not yet implemented)

Returns: