-
Notifications
You must be signed in to change notification settings - Fork 25
QUfitting
The QUfitting module can be used to fit models to Stokes Q and U, extract model parameters, and compare goodness of fit between different models. This takes the form of a single tool that fits a user supplied model and reports the best fit parameters and several goodness of fit metrics, allowing the user to re-run with different models and compare the resulting metrics.
QUfit is suitable for individual polarization spectra; no 3D version has been implemented or planned, in part due to the computational cost of the fitting process.
QUfit using multinest sampling to explore the parameter space, find the best fit, and estimate the uncertainties in the fitted parameters. To do this, QUfit uses the bilby module, which provides an interface to a variety of multinest sampling packages. The default sampler used by Bilby is dynesty, but other samplers can be installed and used seamlessly. Based on experiments on the performance of different samplers, we have found dynesty and nestle work well in the cases we've tested.
If RM-Tools has been installed with terminal commands, then QUfitting can be invoked from the command line as qufit <datafile> <options>, otherwise it must be invoked as a Python script as python <RMtools directory>/RMtools_1D/do_QUfit_1D_mnset.py <datafile> <options>.
Invoking with the -h flag set will produce the help text with the list of option flags.
QUfit is not currently designed to be runable/importable in a script. This functionality could be developed if requested.
QUfit expect an ASCII file in the same format as for RMsynth1D, with 7 columns: frequency (in Hz), Stokes I, Stokes Q, Stokes U, error in I, error in Q, error in U. Each row is a single channel. No header is expected.
If no Stokes I data is available, RMsynth1D will also accept a file with 5 columns: frequency (in Hz), Stokes Q, Stokes U, error in Q, error in U.
Examples of these files can be generated by running mk_test_ascii_data.py in the RMtools_1D/ directory.
QUfit creates a directory with the same base name as the input file supplied to it. Inside this directory several files are saved. Some are associated with the multinest code and are not needed normally: ev.dat, live.points, phys_live.points, post_equal_weights.dat, and resume.dat.
The files that present the results are: stats.dat, fig_m1_corner.pdf, and fig_m1_specfit.pdf.
- stats.dat gives the fitted parameter values and uncertainties.
- fig_m1_corner.pdf shows the marginal probability distributions for all the model parameters
- fig_m1_specfit.pdf shows the frequency spectra and the best fit
Setting the -v flag activates verbose output, which causes QUfit to print the results to the terminal. In the current version of QUfit, this includes some additional outputs that aren't currently saved to files (this will be rectified in a future update). The outputs include the best-fit parameter values with uncertainties, and various quality of fit metrics (including the Akaike Information Criterion, AIC, and Bayesian Information Criterion, BIC).
-h Displays the help documentation to the terminal.
-m MODELNUM Supply the number of the model that you want to fit (default 1). The default models are stored in <RMtools directory>/RM-Tools/RMtools_1D/models_ns. QUfit will first look for as models_ns/ directory in the working directly when it is being run, and if it does not find the model there it will look in the default directly.
-f FUNCTION Supply a string to specify the fitting function for Stokes I: 'linear' polynomial or 'log' polynomial (power law)
-o POLYORD Supply an integer to set the order of polynomial with which to fit the Stokes I spectrum. Values from -5 to 5 are accepted: positive values will fit that order of polynomial, negative values will use dynamic order selection starting with zero and increasing up to a maximum of the equivalent positive order. Default is 2 (quadratic fit).
-i Ignore the Stokes I spectrum (do not divide Stokes Q and U by the fitted Stokes I spectrum). Default behaviour is to fit and divide the Stokes I spectrum.
-d Turn on debugging messages.
-p Output the figures to screen. Default is to save the plots without displaying them.
-v Verbose mode: output the fitting results to the terminal.
--sampler SAMPLER Which sampler to use with Bilby (see https://lscsoft.docs.ligo.org/bilby/samplers.html) Default is dynesty.
--ncores NCORES Number of cores to use for sampling [1].
--nlive NLIVE Number of live points to use for sampling [300].
A full description of polarization models is beyond the scope of this wiki; such descriptions can be found in papers like Burn (1966), Sokoloff et al. (1998), and Tribble (1991).
There are currently 10 models implemented in RM-Tools:
- m1: Simple Faraday thin source
- m11: Two Faraday thin components
- m111: Three Faraday thin components
- m2: Single Faraday component with Burn depolarization (Gaussian model)
- m3: Two Faraday components with a common Burn depolarization term (single Gaussian sigma)
- m4: Two Faraday components with different Burn depolarization terms (individual sigma values)
- m5: Single component uniform slab (Burn slab)
- m6: Two components of uniform slab
- m7: Single component of turbulent mixed medium (i.e. both large-scale and turbulent magnetic components in the emitting volume)
- m12: Single component of turbulent mixed medium with foreground Burn (turbulent) screen
These models have default limits set on the parameter space they explore, particularly for the RM and sigmaRM parameters. Make sure to check that these limits are appropriate for your data, and modify as necessary before running! The limits are set in the model files as the bounds keyword for the relevant parameters.
Creating new models for QUfitting is straightforward; copying and modifying an existing model is the easiest way. Each model has two parts: a model equation (given by the model function at the top of the file), and parameter definitions.
The model equation takes in an array of parameters and the array of channel lambda^2 values, and returns an array of complex polarization values (q+iu). The order of the parameters in the array is the same as the order below in the parameter definitions.
The parameter definitions list the parameters, their short names and labels (when used as plot axes), an initial value for the fit, the bounds of the acceptable range of values, the type of prior on the distribution of values, and whether the values should be allowed to wrap from the upper bound back to the lower bound (like how angles wrap back to zero after a full revolution).
New models are welcome for addition to the default list, so please contact the package maintainers about incorporating your models.