Commit 27e12115 authored by Gregory Ashton's avatar Gregory Ashton
Browse files

Updates docs fixing errors in fc on glitching data

parent 115bcbac
...@@ -5,14 +5,11 @@ search using MCMC](fully_coherent_search_using_MCMC.md), to the glitching signal ...@@ -5,14 +5,11 @@ search using MCMC](fully_coherent_search_using_MCMC.md), to the glitching signal
[make fake data](make_fake_data.md]). The aim here is to illustrate the effect [make fake data](make_fake_data.md]). The aim here is to illustrate the effect
such a signal can have on a fully-coherent search. The complete script for this such a signal can have on a fully-coherent search. The complete script for this
example canbe found example canbe found
[here](../example/fully_cohrent_search_on_glitching_data.py). [here](../example/fully_coherent_search_using_MCMC_on_glitching_data.py).
We use the same prior as in the basic fully-coherent search, except that we After importing `pyfstat`, We setup a flat prior on `F0` and `F1`, based on the
will modify the prior on `F0` to a flat uniform prior. The reason for this is values used to generate the signal:
to highlight the multimodal nature of the posterior which results from the
glitch (a normal prior centered on one of the modes will bias one mode over
the other). So our initial set up is
``` ```
from pyfstat import MCMCSearch from pyfstat import MCMCSearch
...@@ -26,53 +23,57 @@ tref = 362750407.0 ...@@ -26,53 +23,57 @@ tref = 362750407.0
tstart = 1000000000 tstart = 1000000000
duration = 100*86400 duration = 100*86400
tend = tstart = duration tend = tstart + duration
theta_prior = {'F0': {'type': 'unif', 'lower': F0-5e-5, theta_prior = {'F0': {'type': 'unif', 'lower': F0-1e-4, 'upper': F0+1e-4},
'upper': F0+5e-5}, 'F1': {'type': 'unif', 'lower': F1*(1+1e-3), 'upper': F1*(1-1e-3)},
'F1': {'type': 'norm', 'loc': F1, 'scale': abs(1e-6*F1)},
'F2': F2, 'F2': F2,
'Alpha': Alpha, 'Alpha': Alpha,
'Delta': Delta 'Delta': Delta
} }
``` ```
Next, we will use 10 temperatures and a larger number of walkers - these have In this search, we will use paralllel tempering (to help the walkers move
been tuned to illustrate the bimodal nature of the posterior between the different peaks in the posterior).
``` ```
ntemps = 10 ntemps = 2
betas = np.logspace(0, -30, ntemps) log10temperature_min = -0.01
nwalkers = 500 nwalkers = 100
nsteps = [100, 100, 100] nsteps = [5000, 10000]
mcmc = MCMCSearch('fully_coherent_on_glitching_data', 'data', mcmc = MCMCSearch('fully_coherent_search_using_MCMC_on_glitching_data', 'data',
sftfilepath='data/*_glitch*.sft', sftfilepath='data/*_glitch*.sft',
theta_prior=theta_prior, tref=tref, tstart=tstart, tend=tend, theta_prior=theta_prior, tref=tref, tstart=tstart, tend=tend,
nsteps=nsteps, nwalkers=nwalkers, ntemps=ntemps, nsteps=nsteps, nwalkers=nwalkers, ntemps=ntemps,
log10temperature_min=log10temperature_min, scatter_val=1e-6) log10temperature_min=log10temperature_min)
mcmc.run() mcmc.run()
mcmc.plot_corner(add_prior=True) mcmc.plot_corner(add_prior=True)
``` ```
Running this takes slightly longer than the basic example (a few minutes) and Running this example, we obtain traces of the walkers like this:
produces a multimodal posterior: ![](img/fully_coherent_search_using_MCMC_on_glitching_data_walkers.png)
![](img/fully_coherent_on_glitching_data_corner.png)
Clearly one central peak pertains to the original frequency `F0=30`. At Although it is not obvious at first, the large widths of these traces in fact
`30+0.4e-5`, we find the second largest peak - this is the mode corresponding show that the walkers are jumping between two bimodal peaks (for both `F0` and
to the second half of the data. In reality, we would expect both peaks to be `F1): this is possible due to the tuning of the parallel tempering. To see this
of equal size since the glitch occurs exactly half way through (see [how the clearly, we also plot the corner plot:
data was made](make_fake_data.md)). We will confirm this later on by performing ![](img/fully_coherent_search_using_MCMC_on_glitching_data_corner.png)
a grid search.
Finally, the maximum twoF value found is From this corner plot, we that unlike the in the [single glitch fully-coherent
search](full_coherent_search_using_MCMC.md), the posterior, even after a large
number of steps, is multimodal. However, these two peaks do **not** correspond
exactly to the two frequencies before and after the glitch, which would be
`30` and `30+4e5` (to see this, see how the data is
[generated](../examples/make_dake_data.py)). This is partly due to the noise
and partly due to the fact that the maximum detection statistic in the case
of glitches can occur at point *in between* the two frequencies. Moreover, we
see bimodality in `F1`, which did does not change during the glitch.
``` ```
>>> mcmc.print_summary() >>> mcmc.print_summary()
Max twoF: 411.595 Max twoF: 1354.7
``` ```
That is, compared to the basic search (on a smooth signal) which had a twoF of That is, compared to the basic search (on a smooth signal) which had a twoF of
`1756.44177246` (in agreement with the predicted twoF), we have lost a large `~1764` (in agreement with the predicted twoF), we have lost a large
fraction of the SNR due to the glitch. fraction of the SNR due to the glitch.
...@@ -18,10 +18,10 @@ theta_prior = {'F0': {'type': 'unif', 'lower': F0-1e-4, 'upper': F0+1e-4}, ...@@ -18,10 +18,10 @@ theta_prior = {'F0': {'type': 'unif', 'lower': F0-1e-4, 'upper': F0+1e-4},
'Delta': Delta 'Delta': Delta
} }
ntemps = 10 ntemps = 2
log10temperature_min = -30 log10temperature_min = -0.01
nwalkers = 500 nwalkers = 100
nsteps = [100, 100, 100] nsteps = [5000, 10000]
mcmc = MCMCSearch('fully_coherent_search_using_MCMC_on_glitching_data', 'data', mcmc = MCMCSearch('fully_coherent_search_using_MCMC_on_glitching_data', 'data',
sftfilepath='data/*_glitch*.sft', sftfilepath='data/*_glitch*.sft',
......
...@@ -24,7 +24,9 @@ data = Writer( ...@@ -24,7 +24,9 @@ data = Writer(
F2=F2, duration=duration, Alpha=Alpha, Delta=Delta, h0=h0, sqrtSX=sqrtSX) F2=F2, duration=duration, Alpha=Alpha, Delta=Delta, h0=h0, sqrtSX=sqrtSX)
data.make_data() data.make_data()
print 'Predicted fstat value:', data.predict_fstat() # The predicted twoF, given by lalapps_predictFstat can be accessed by
twoF = data.predict_fstat()
print 'Predicted twoF value: {}\n'.format(twoF)
# Next, taking the same signal parameters, we include a glitch half way through # Next, taking the same signal parameters, we include a glitch half way through
dtglitch = duration/2.0 dtglitch = duration/2.0
...@@ -37,11 +39,6 @@ glitch_data = Writer( ...@@ -37,11 +39,6 @@ glitch_data = Writer(
dtglitch=dtglitch, delta_F0=delta_F0, delta_F1=delta_F1, detector='L1') dtglitch=dtglitch, delta_F0=delta_F0, delta_F1=delta_F1, detector='L1')
glitch_data.make_data() glitch_data.make_data()
# The predicted twoF, given by lalapps_predictFstat can be accessed by
print 'Predicted fstat value:', data.predict_fstat()
# Making data with two glitches # Making data with two glitches
dtglitch = [duration/4.0, 4*duration/5.0] dtglitch = [duration/4.0, 4*duration/5.0]
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment