diff --git a/docs/fully_coherent_search_using_MCMC_on_glitching_data.md b/docs/fully_coherent_search_using_MCMC_on_glitching_data.md
index 04ff8af78f0cbe7b25db64e0ced2296170f70277..7f2ea890e07a9e2cb37dcaef45bd5f2a9193800f 100644
--- a/docs/fully_coherent_search_using_MCMC_on_glitching_data.md
+++ b/docs/fully_coherent_search_using_MCMC_on_glitching_data.md
@@ -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
such a signal can have on a fully-coherent search. The complete script for this
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
-will modify the prior on `F0` to a flat uniform prior. The reason for this is
-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
+After importing `pyfstat`, We setup a flat prior on `F0` and `F1`, based on the
+values used to generate the signal:
```
from pyfstat import MCMCSearch
@@ -26,53 +23,57 @@ tref = 362750407.0
tstart = 1000000000
duration = 100*86400
-tend = tstart = duration
+tend = tstart + duration
-theta_prior = {'F0': {'type': 'unif', 'lower': F0-5e-5,
- 'upper': F0+5e-5},
- 'F1': {'type': 'norm', 'loc': F1, 'scale': abs(1e-6*F1)},
+theta_prior = {'F0': {'type': 'unif', 'lower': F0-1e-4, 'upper': F0+1e-4},
+ 'F1': {'type': 'unif', 'lower': F1*(1+1e-3), 'upper': F1*(1-1e-3)},
'F2': F2,
'Alpha': Alpha,
'Delta': Delta
}
```
-Next, we will use 10 temperatures and a larger number of walkers - these have
-been tuned to illustrate the bimodal nature of the posterior
-
+In this search, we will use paralllel tempering (to help the walkers move
+between the different peaks in the posterior).
```
-ntemps = 10
-betas = np.logspace(0, -30, ntemps)
-nwalkers = 500
-nsteps = [100, 100, 100]
+ntemps = 2
+log10temperature_min = -0.01
+nwalkers = 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',
theta_prior=theta_prior, tref=tref, tstart=tstart, tend=tend,
nsteps=nsteps, nwalkers=nwalkers, ntemps=ntemps,
- log10temperature_min=log10temperature_min, scatter_val=1e-6)
+ log10temperature_min=log10temperature_min)
mcmc.run()
mcmc.plot_corner(add_prior=True)
```
-Running this takes slightly longer than the basic example (a few minutes) and
-produces a multimodal posterior:
-
+Running this example, we obtain traces of the walkers like this:
+
-Clearly one central peak pertains to the original frequency `F0=30`. At
-`30+0.4e-5`, we find the second largest peak - this is the mode corresponding
-to the second half of the data. In reality, we would expect both peaks to be
-of equal size since the glitch occurs exactly half way through (see [how the
-data was made](make_fake_data.md)). We will confirm this later on by performing
-a grid search.
+Although it is not obvious at first, the large widths of these traces in fact
+show that the walkers are jumping between two bimodal peaks (for both `F0` and
+`F1): this is possible due to the tuning of the parallel tempering. To see this
+clearly, we also plot the corner plot:
+
-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()
-Max twoF: 411.595
+Max twoF: 1354.7
```
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.
diff --git a/docs/img/fully_coherent_search_using_MCMC_on_glitching_data_corner.png b/docs/img/fully_coherent_search_using_MCMC_on_glitching_data_corner.png
new file mode 100644
index 0000000000000000000000000000000000000000..c6002bd29400430f3bc46dbe13d66a095e36f059
Binary files /dev/null and b/docs/img/fully_coherent_search_using_MCMC_on_glitching_data_corner.png differ
diff --git a/docs/img/fully_coherent_search_using_MCMC_on_glitching_data_walkers.png b/docs/img/fully_coherent_search_using_MCMC_on_glitching_data_walkers.png
new file mode 100644
index 0000000000000000000000000000000000000000..e54b4a441f494db7ead5ade1793846d664a9e056
Binary files /dev/null and b/docs/img/fully_coherent_search_using_MCMC_on_glitching_data_walkers.png differ
diff --git a/examples/fully_coherent_search_using_MCMC_on_glitching_data.py b/examples/fully_coherent_search_using_MCMC_on_glitching_data.py
index 9e5992a01c759a8be0593a2b59459dbbaf774fdb..e7fbeea1251924f6416f8816a8dea8fff67d1103 100644
--- a/examples/fully_coherent_search_using_MCMC_on_glitching_data.py
+++ b/examples/fully_coherent_search_using_MCMC_on_glitching_data.py
@@ -18,10 +18,10 @@ theta_prior = {'F0': {'type': 'unif', 'lower': F0-1e-4, 'upper': F0+1e-4},
'Delta': Delta
}
-ntemps = 10
-log10temperature_min = -30
-nwalkers = 500
-nsteps = [100, 100, 100]
+ntemps = 2
+log10temperature_min = -0.01
+nwalkers = 100
+nsteps = [5000, 10000]
mcmc = MCMCSearch('fully_coherent_search_using_MCMC_on_glitching_data', 'data',
sftfilepath='data/*_glitch*.sft',
diff --git a/examples/make_fake_data.py b/examples/make_fake_data.py
index eeadca9108ea1cea366e786e0029b384ad48eab0..47072ce87a0ef622774436db581c0b676a065b47 100644
--- a/examples/make_fake_data.py
+++ b/examples/make_fake_data.py
@@ -24,7 +24,9 @@ data = Writer(
F2=F2, duration=duration, Alpha=Alpha, Delta=Delta, h0=h0, sqrtSX=sqrtSX)
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
dtglitch = duration/2.0
@@ -37,11 +39,6 @@ glitch_data = Writer(
dtglitch=dtglitch, delta_F0=delta_F0, delta_F1=delta_F1, detector='L1')
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
dtglitch = [duration/4.0, 4*duration/5.0]