Parameter estimation code for known pulsar searches using the nested sampling algorithm. More...
Prototypes | |
INT4 | main (INT4 argc, CHAR *argv[]) |
Parameter estimation code for known pulsar searches using the nested sampling algorithm.
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.
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] .
The usage format is given below and can also be found by running the code with
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:
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.
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.
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:
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 | |
LALStringVector * | corlist = NULL |
Definition at line 146 of file pulsar_parameter_estimation_nested.c.
LALStringVector* corlist = NULL |
Definition at line 144 of file pulsar_parameter_estimation_nested.c.