diff --git a/README.md b/README.md index 3cfa083381ff560f165f1a74d12445f5ce1c1f25..3e5b51a6841eb8e248159f9d6c4ce88c52a560d8 100644 --- a/README.md +++ b/README.md @@ -4,6 +4,14 @@ This is a python package containing basic wrappers of the `lalpulsar` module with capabilities to perform a variety of searches, primarily focussing on semi-coherent glitch searches. +## Examples + +All examples can be run from their source scripts in [examples](examples), or +for each example there is descriptive documentation: + +* [Making fake data with and without glitches](docs/make_fake_data.md) +* [Fully coherent searches MCMC](docs/fully_coherent_search.md) + ## Installation The script can be installed system wide via diff --git a/docs/fully_coherent_search.md b/docs/fully_coherent_search.md new file mode 100644 index 0000000000000000000000000000000000000000..6f274d60077c5f63fa3080a308e6d14ddec88eec --- /dev/null +++ b/docs/fully_coherent_search.md @@ -0,0 +1,116 @@ +# Fully coherent search using MCMC + +In this example, we will show the basics of setting up and running an MCMC +search for a fully-coherent search. This is based on the example +[fully_coherent_search.py](../example/fully_coherent_search.py). We will run +the search on the `basic` data generated in the +[make_fake_data](make_fake_data.md) example. + +We first need to import the search tool, in this example we will use the +`MCMCSearch`, but one could equally use `MCMCGlitchSearch` with `nglitch=0`. +We first import this, + +``` +from pyfstat import MCMCSearch +``` + +Next, we define some variables defining the center of our search, namely the +Crab parameters (with which the `basic` data was produced): + +``` +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 +``` + +Next, we specify our prior, this is done using a dictionary: + +``` +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} +``` +Each key and value of the `theta_prior` contains an instruction to the MCMC +search. If the value is a scalar, the MCMC search holds these fixed (as is the +case for `F2`, `Alpha`, and `Delta` here). If instead the value is a dictionary +describing a distribution, this is taken as the prior and the variable is +simulated in the MCMC search (as is the case for `F0` and `F1`). Note that +for `MCMCSearch`, `theta_prior` must contain at least all of the variables +given here (even if they are zero), and if `binary=True`, it must also contain +the binary parameters. + +Next, we define the parameters of the MCMC search: + +``` +ntemps = 1 +nwalkers = 100 +nsteps = [100, 500, 1000] +``` + +These can be considered the *tuning parameters* of the search. In this instance +we use a single temperature with 100 walkers. For the number of steps we use +a 2-step process. First we run the simulation starting from values picked +about the prior for `100` steps, afterwhich we take the maximum likelihood +chain and pick a new set of walkers scattered about this. This step removes +walkers which get lost is islands of low probability (it can be removed by +simply called `nsteps = [500, 1000]`). Finally, the simulation run for `500` +steps of burn-in then `1000` steps of production to estimate the posterior. + +Passing all this to the MCMC search, we also need to give it a label and +directory to save the data and provide `sftlabel` and `sftdir` which defines +which data to use in the search + +``` +mcmc = MCMCSearch('fully_coherent', 'data', sftlabel='basic', sftdir='data', + theta_prior=theta_prior, tref=tref, tstart=tstart, tend=tend, + nsteps=nsteps, nwalkers=nwalkers, ntemps=ntemps, + scatter_val=1e-10) +``` + +To run code, we call + +``` +mcmc.run() +``` + +This produces two `png` images. The first is the chains after the +initialisation step + +Here one can see there are multiple evenly spaced solutions (at least in +frequency) corresponding to different peaks +in the likelihood. Taking the largest peak, the walkers are reinitialised about +this point (scattered in a relative way by the `scatter_val`) and run again +producing + +Now once can see that the chains have converged onto the single central peak. + +To get posteriors, we call + +``` +mcmc.plot_corner() +``` +which produces a corner plot + +illustrating the tightly constrained posteriors on `F0` and `F1` and their +covariance. Furthermore, one may wish to get a summary which can be printed +to the terminal via + +``` +mcmc.print_summary() +``` +which gives the maximum twoF value, median and standard-deviation, in this case +this is +``` +Max twoF: 1756.44177246 +F0 = 2.999999974e+01 +/- 6.363964377e-08 +F1 = -9.999960025e-11 +/- 9.920455272e-17 +``` diff --git a/docs/img/fully_coherent_corner.png b/docs/img/fully_coherent_corner.png new file mode 100644 index 0000000000000000000000000000000000000000..c9a97c32f7515e7eccc3cc8d428eee548b7c68aa Binary files /dev/null and b/docs/img/fully_coherent_corner.png differ diff --git a/docs/img/fully_coherent_init_0_walkers.png b/docs/img/fully_coherent_init_0_walkers.png new file mode 100644 index 0000000000000000000000000000000000000000..b2780c0dcff011fa323b6d16d2ef9e81aa213773 Binary files /dev/null and b/docs/img/fully_coherent_init_0_walkers.png differ diff --git a/docs/img/fully_coherent_walkers.png b/docs/img/fully_coherent_walkers.png new file mode 100644 index 0000000000000000000000000000000000000000..128f1fe57414b62d3c764391e6342766c12d4c93 Binary files /dev/null and b/docs/img/fully_coherent_walkers.png differ diff --git a/examples/Makefile b/examples/Makefile index eccd363185bc8cb4695b8caab89b1f25057ff8ff..888f4dc85e2b88bef574605d288fd4b926a09506 100644 --- a/examples/Makefile +++ b/examples/Makefile @@ -1,5 +1,8 @@ basic_sft = data/H-4800_H1_1800SFT_basic-1000000000-8640000.sft glitch_sft = data/H-4800_H1_1800SFT_glitch-1000000000-8640000.sft +data/fully_coherent_corner.png : $(basic_sft) fully_coherent_search.py + python fully_coherent_search.py + $(basic_sft) $(glitch_sft): make_fake_data.py python make_fake_data.py diff --git a/examples/fully_coherent_search.py b/examples/fully_coherent_search.py new file mode 100644 index 0000000000000000000000000000000000000000..8e6b4227d262f8ec67d3cb9a8089a2381662f345 --- /dev/null +++ b/examples/fully_coherent_search.py @@ -0,0 +1,31 @@ +from pyfstat import MCMCSearch + +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 + } + +ntemps = 1 +nwalkers = 100 +nsteps = [100, 500, 1000] + +mcmc = MCMCSearch('fully_coherent', 'data', sftlabel='basic', sftdir='data', + theta_prior=theta_prior, tref=tref, tstart=tstart, tend=tend, + nsteps=nsteps, nwalkers=nwalkers, ntemps=ntemps, + scatter_val=1e-10) +mcmc.run() +mcmc.plot_corner() +mcmc.print_summary() diff --git a/pyfstat.py b/pyfstat.py index 22c919f6efcae2fcb8378051dfbfc59144ad4328..4a594f295d0db996c9a447b75cad5d7dec028e74 100755 --- a/pyfstat.py +++ b/pyfstat.py @@ -711,13 +711,13 @@ class MCMCSearch(BaseSearchClass): if ndim > 1: for i in range(ndim): + axes[i].ticklabel_format(useOffset=False, axis='y') cs = chain[:, start:stop, i].T axes[i].plot(cs, color="k", alpha=alpha) if symbols: axes[i].set_ylabel(symbols[i]) if draw_vline is not None: axes[i].axvline(draw_vline, lw=2, ls="--") - axes[i].ticklabel_format(useOffset=False, axis='y') else: cs = chain[:, start:stop, 0].T