Commit 92956a08 authored by Gregory Ashton's avatar Gregory Ashton
Browse files

Initial commit adding documentation for semi-coherent glitch search

parent 45faf9c1
...@@ -54,8 +54,8 @@ Running this example, we obtain traces of the walkers like this: ...@@ -54,8 +54,8 @@ Running this example, we obtain traces of the walkers like this:
![](img/fully_coherent_search_using_MCMC_on_glitching_data_walkers.png) ![](img/fully_coherent_search_using_MCMC_on_glitching_data_walkers.png)
Although it is not obvious at first, the large widths of these traces in fact 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 show that the walkers are jumping between multiple peaks (for both `F0` and
`F1): this is possible due to the tuning of the parallel tempering. To see this `F1): this is possible, due to the tuning of the parallel tempering. To see this
clearly, we also plot the corner plot: clearly, we also plot the corner plot:
![](img/fully_coherent_search_using_MCMC_on_glitching_data_corner.png) ![](img/fully_coherent_search_using_MCMC_on_glitching_data_corner.png)
...@@ -71,7 +71,7 @@ see bimodality in `F1`, which did does not change during the glitch. ...@@ -71,7 +71,7 @@ see bimodality in `F1`, which did does not change during the glitch.
``` ```
>>> mcmc.print_summary() >>> mcmc.print_summary()
Max twoF: 1354.7 Max twoF: 422.97
``` ```
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
`~1764` (in agreement with the predicted twoF), we have lost a large `~1764` (in agreement with the predicted twoF), we have lost a large
......
# Semi-coherent glitch search on data with a single glitch using MCMC
In this example, based on [this
script](../examples/semi_coherent_glitch_search_using_MCMC.py), we show the
basic setup for a single-glitch search. We begin, in the usual way, with
defining some the prior
```python
import pyfstat
F0 = 30.0
F1 = -1e-10
F2 = 0
Alpha = 5e-3
Delta = 6e-2
tref = 362750407.0
tstart = 1000000000
duration = 100*86400
tend = tstart + duration
theta_prior = {'F0': {'type': 'norm', 'loc': F0, 'scale': abs(1e-6*F0)},
'F1': {'type': 'norm', 'loc': F1, 'scale': abs(1e-6*F1)},
'F2': F2,
'Alpha': Alpha,
'Delta': Delta,
'delta_F0': {'type': 'halfnorm', 'loc': 0,
'scale': 1e-5*F0},
'delta_F1': 0,
'tglitch': {'type': 'unif',
'lower': tstart+0.1*duration,
'upper': tstart+0.9*duration},
}
```
For simplicity, we have chosen a prior based on the known inputs. The important
steps here are the definition of `delta_F0`, `delta_F1` and `tglitch`, the
prior densities for the glitch-parameters. We then use a parallel-tempered
set-up, in addition to an initialisation step and run the search:
```python
ntemps = 4
log10temperature_min = -1
nwalkers = 100
nsteps = [5000, 1000, 1000]
mcmc = pyfstat.MCMCGlitchSearch(
'semi_coherent_glitch_search_using_MCMC', 'data',
sftfilepath='data/*_glitch*sft', theta_prior=theta_prior, tref=tref,
tstart=tstart, tend=tend, nsteps=nsteps, nwalkers=nwalkers,
scatter_val=1e-10, nglitch=1, ntemps=ntemps,
log10temperature_min=log10temperature_min)
mcmc.run()
mcmc.plot_corner(add_prior=True)
mcmc.print_summary()
```
The posterior for this search demonstrates that we recover the input parameters:
![](img/semi_coherent_search_using_MCMC_corner.png)
# Semi-coherent glitch search on data with two-glitches using MCMC
![](img/semi_coherent_twoglitch_search_walkers.png)
![](img/semi_coherent_twoglitch_search_corner.png)
...@@ -36,7 +36,7 @@ delta_F1 = 0 ...@@ -36,7 +36,7 @@ delta_F1 = 0
glitch_data = Writer( glitch_data = Writer(
label='glitch', outdir='data', tref=tref, tstart=tstart, F0=F0, F1=F1, label='glitch', outdir='data', tref=tref, tstart=tstart, F0=F0, F1=F1,
F2=F2, duration=duration, Alpha=Alpha, Delta=Delta, h0=h0, sqrtSX=sqrtSX, F2=F2, duration=duration, Alpha=Alpha, Delta=Delta, h0=h0, sqrtSX=sqrtSX,
dtglitch=dtglitch, delta_F0=delta_F0, delta_F1=delta_F1, detector='L1') dtglitch=dtglitch, delta_F0=delta_F0, delta_F1=delta_F1)
glitch_data.make_data() glitch_data.make_data()
# Making data with two glitches # Making data with two glitches
......
...@@ -24,13 +24,17 @@ theta_prior = {'F0': {'type': 'norm', 'loc': F0, 'scale': abs(1e-6*F0)}, ...@@ -24,13 +24,17 @@ theta_prior = {'F0': {'type': 'norm', 'loc': F0, 'scale': abs(1e-6*F0)},
'upper': tstart+0.9*duration}, 'upper': tstart+0.9*duration},
} }
nwalkers = 500 ntemps = 4
nsteps = [1000, 1000, 1000] log10temperature_min = -1
nwalkers = 100
nsteps = [5000, 1000, 1000]
mcmc = pyfstat.MCMCGlitchSearch( mcmc = pyfstat.MCMCGlitchSearch(
'semi_coherent_glitch_search', 'data', sftfilepath='data/*_glitch*sft', 'semi_coherent_glitch_search_using_MCMC', 'data',
theta_prior=theta_prior, tref=tref, tstart=tstart, tend=tend, sftfilepath='data/*_glitch*sft', theta_prior=theta_prior, tref=tref,
nsteps=nsteps, nwalkers=nwalkers, scatter_val=1e-10, nglitch=1) tstart=tstart, tend=tend, nsteps=nsteps, nwalkers=nwalkers,
scatter_val=1e-10, nglitch=1, ntemps=ntemps,
log10temperature_min=log10temperature_min)
mcmc.run() mcmc.run()
mcmc.plot_corner(add_prior=True) mcmc.plot_corner(add_prior=True)
......
...@@ -16,9 +16,12 @@ theta_prior = {'F0': {'type': 'norm', 'loc': F0, 'scale': abs(1e-6*F0)}, ...@@ -16,9 +16,12 @@ theta_prior = {'F0': {'type': 'norm', 'loc': F0, 'scale': abs(1e-6*F0)},
'F2': F2, 'F2': F2,
'Alpha': Alpha, 'Alpha': Alpha,
'Delta': Delta, 'Delta': Delta,
'delta_F0': {'type': 'halfnorm', 'loc': 0, 'delta_F0_0': {'type': 'halfnorm', 'loc': 0,
'scale': 1e-7*F0}, 'scale': 1e-7*F0},
'delta_F1': 0, 'delta_F0_1': {'type': 'halfnorm', 'loc': 0,
'scale': 1e-7*F0},
'delta_F1_0': 0,
'delta_F1_1': 0,
'tglitch_0': {'type': 'unif', 'tglitch_0': {'type': 'unif',
'lower': tstart+0.01*duration, 'lower': tstart+0.01*duration,
'upper': tstart+0.5*duration}, 'upper': tstart+0.5*duration},
...@@ -27,8 +30,8 @@ theta_prior = {'F0': {'type': 'norm', 'loc': F0, 'scale': abs(1e-6*F0)}, ...@@ -27,8 +30,8 @@ theta_prior = {'F0': {'type': 'norm', 'loc': F0, 'scale': abs(1e-6*F0)},
'upper': tstart+0.99*duration}, 'upper': tstart+0.99*duration},
} }
nwalkers = 50 nwalkers = 100
nsteps = [500, 500, 500] nsteps = [1000, 1000, 5000]
mcmc = pyfstat.MCMCGlitchSearch( mcmc = pyfstat.MCMCGlitchSearch(
'semi_coherent_twoglitch_search', 'data', sftfilepath='data/*twoglitch*sft', 'semi_coherent_twoglitch_search', 'data', sftfilepath='data/*twoglitch*sft',
......
Supports Markdown
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