Loading [MathJax]/extensions/TeX/AMSsymbols.js
LALPulsar 7.1.1.1-8a6b96f
All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
pulsar_parameter_estimation_nested.c File Reference

Parameter estimation code for known pulsar searches using the nested sampling algorithm. More...

Prototypes

INT4 main (INT4 argc, CHAR *argv[])
 

Detailed Description

Parameter estimation code for known pulsar searches using the nested sampling algorithm.

Author
Matthew Pitkin, John Veitch, Max Isi, Colin Gill

Description

This code is used to perform parameter estimation and evidence calculation in targeted/semi-targeted searches for gravitational waves from known pulsars. It may also be used to follow-up on signal candidates from semi-coherent all-sky searches for unknown sources.

It uses the Bayesian technique of 'Nested Sampling' to sample over a defined prior parameter space (unknown signal parameters such as the gravitational wave amplitude). These samples can then be used to create posterior probability density functions of the unknown parameters. The samples are also used to calculate the Bayesian evidence, also known as the marginal likelihood, for a given signal model. This can be compared with other models, in particular the model that the data is described by Gaussian noise alone.

As input the code requires time domain data that has been heterodyned using the known (or close to) phase evolution of the pulsar. The time domain input should consist of a three column text file containing the GPS time stamp of the data point, the real part of the heterodyned data and the imaginary part of the heterodyned data, e.g.

900000000.000000 1.867532e-24 -7.675329e-25
900000060.000000 2.783651e-24 3.654386e-25
...

Most commonly such data will have a sample rate of 1/60 Hz, giving a bandwidth of the same amount, but the code can accept any rate.

The code also requires that you specify which parameters are to be searched over, and the prior ranges over these. Any of the signal parameters can be searched over, including frequency, sky position and binary system parameters, although the bandwidth of the data and search efficiency need to be taken into account.

The 'Nested Sampling' algorithm (developed by [26]) used is that defined in LALinferenceNestedSampler.c (see [28]). It is essentially an efficient way to perform the integral

\[ Z = \int^{\mathbf{\theta}} p(d|\mathbf{\theta}) p(\mathbf{\theta}) \mathrm{d}\mathbf{\theta}, \]

where \( \mathbf{\theta} \) is a vector of parameters, \( p(d|\mathbf{\theta}) \) is the likelihood of the data given the parameters, and \( p(\mathbf{\theta}) \) is the prior on the parameters. It does this by changing the multi-dimensional integral over N parameters into a one-dimensional integral

\[ Z = \int^X L(X) \mathrm{d}X \approx \sum_i L(X_i) \Delta{}X_i, \]

where \( L(X) \) is the likelihood, and \( X \) is the prior mass. The algorithm will draw a number ( \( N \) ) of samples (live points) from the parameter priors, calculate the likelihood for each point and find the lowest likelihood value. The lowest likelihood value will be added to the summation in the above equation, with \( \log{\Delta{}X_i} \approx 1/N \) coming from the fact that the prior would be normalised to unity and therefore each point should occupy an equal fraction and at each iteration the prior volume will decrease geometrically (for \( \log{\Delta{}X_0} = 0 \) ). A new point is then drawn from the prior with the criterion that it has a higher likelihood than the previous lowest point and substitutes that point. To draw the new point a Markov Chain Monte Carlo (MCMC) procedure is used. The procedure is continued until a stopping criterion is reached, which in this case is that the remaining prior volume is less than the tolerance value set (see below). The implementation of this can be seen in [28] .

Usage

The usage format is given below and can also be found by running the code with

lalpulsar_parameter_estimation_nested --help

An example of running the code to search over the four unknown parameters \( h_0 \) , \( \phi_0 \) , \( \psi \) and \( \cos{\iota} \) , for pulsar J0534-2200, given heterodyned time domain data from the H1 detector in the file finehet_J0534-2200_H1, is:

lalpulsar_parameter_estimation_nested --detectors H1 --par-file J0534-2200.par --input-files finehet_J0534-2200_H1 --outfile ns_J0534-2200.hdf --prior-file prior_J0534-2200.txt --ephem-earth lscsoft/share/lalpulsar/earth05-09.dat --ephem-sun lscsoft/share/lalpulsar/sun05-09.dat --Nlive 1000 --Nmcmcinitial 0 --tolerance 0.25
LALDetector detectors[LAL_NUM_DETECTORS]

The par-file is a TEMPO(2)-style file containing the parameters of the pulsar used to perform the heterodyne (the frequency parameters are the rotation frequency and therefore not necessarily the gravitational wave frequency) e.g.

RA 12:54:31.87523895
DEC -54:12:43.6572033
PMRA 1.7
PMDEC 2.8
POSEPOCH 54320.8531
F0 123.7896438753
F1 4.592e-15
PEPOCH 54324.8753
DEC
RA
#define F0

The prior-file is a text file containing a list of the parameters to be searched over, the prior type ("uniform" or "gaussian") and their given lower/mean and upper/standard deviation ranges e.g.

h0 uniform 0 1e-21
phi0 uniform 0 6.283185307179586
cosiota uniform -1 1
psi uniform -0.785398163397448 0.785398163397448
double e

Note that if searching over frequency parameters the ranges specified in the prior-file should be given in terms of the pulsars rotation frequency and not necessarily the gravitational wave frequency e.g. for a triaxial star emitting gravitational waves at 100 Hz (which will be at twice the rotation frequency) if you wanted to search over 99.999 to 100.001 Hz then you should used

An example of running the code as above, but this time on fake data created using the Advanced LIGO design noise curves and with a signal injected into the data is:

lalpulsar_parameter_estimation_nested --fake-data AH1 --inject-file fake.par --par-file fake.par --outfile ns_fake.hdf --prior-file prior_fake.txt --ephem-earth lscsoft/share/lalpulsar/earth05-09.dat --ephem-sun lscsoft/share/lalpulsar/sun05-09.dat --Nlive 1000 --Nmcmcinitial 0 --tolerance 0.25
float data[BLOCKSIZE]
Definition: hwinject.c:360

In this case the inject-file parameter file must contain the values of h0, phi0, psi and cosiota, otherwise these will be set to zero by default. The parameter files given for inject-file and par-file do not have to be the same - the injection can be offset from the 'heterodyned' values, which will be reflected in the data. If an inject-output file is also specified then the fake data containing the signal, and a fake signal-only data set, will be output.

Definition in file pulsar_parameter_estimation_nested.c.

Go to the source code of this file.

Variables

LALStringVectorcorlist = NULL
 

Function Documentation

◆ main()

INT4 main ( INT4  argc,
CHAR argv[] 
)

Definition at line 146 of file pulsar_parameter_estimation_nested.c.

Variable Documentation

◆ corlist

LALStringVector* corlist = NULL

Definition at line 144 of file pulsar_parameter_estimation_nested.c.