Note
Click here to download the full example code
Glitch robust MCMC search
MCMC search employing a signal hypothesis allowing for a glitch to be present in the data. The setup corresponds to a targeted search, and the simulated signal contains a single glitch.
10 import numpy as np
11 import pyfstat
12 import time
13 from PyFstat_example_make_data_for_search_on_1_glitch import (
14 tstart,
15 duration,
16 tref,
17 F0,
18 F1,
19 F2,
20 Alpha,
21 Delta,
22 delta_F0,
23 dtglitch,
24 outdir,
25 )
26 import os
27
28 label = "PyFstat_example_glitch_robust_directed_MCMC_search_on_1_glitch"
29
30 Nstar = 1000
31 F0_width = np.sqrt(Nstar) * np.sqrt(12) / (np.pi * duration)
32 F1_width = np.sqrt(Nstar) * np.sqrt(180) / (np.pi * duration**2)
33
34 theta_prior = {
35 "F0": {"type": "unif", "lower": F0 - F0_width / 2.0, "upper": F0 + F0_width / 2.0},
36 "F1": {"type": "unif", "lower": F1 - F1_width / 2.0, "upper": F1 + F1_width / 2.0},
37 "F2": F2,
38 "delta_F0": {"type": "unif", "lower": 0, "upper": 1e-5},
39 "delta_F1": 0,
40 "tglitch": {
41 "type": "unif",
42 "lower": tstart + 0.1 * duration,
43 "upper": tstart + 0.9 * duration,
44 },
45 "Alpha": Alpha,
46 "Delta": Delta,
47 }
48
49 ntemps = 3
50 log10beta_min = -0.5
51 nwalkers = 100
52 nsteps = [250, 250]
53
54 mcmc = pyfstat.MCMCGlitchSearch(
55 label=label,
56 outdir=outdir,
57 sftfilepattern=os.path.join(outdir, "*1_glitch*sft"),
58 theta_prior=theta_prior,
59 tref=tref,
60 minStartTime=tstart,
61 maxStartTime=tstart + duration,
62 nsteps=nsteps,
63 nwalkers=nwalkers,
64 ntemps=ntemps,
65 log10beta_min=log10beta_min,
66 nglitch=1,
67 )
68 mcmc.transform_dictionary["F0"] = dict(
69 subtractor=F0, multiplier=1e6, symbol="$f-f_\\mathrm{s}$"
70 )
71 mcmc.unit_dictionary["F0"] = "$\\mu$Hz"
72 mcmc.transform_dictionary["F1"] = dict(
73 subtractor=F1, multiplier=1e12, symbol="$\\dot{f}-\\dot{f}_\\mathrm{s}$"
74 )
75 mcmc.unit_dictionary["F1"] = "$p$Hz/s"
76 mcmc.transform_dictionary["delta_F0"] = dict(
77 multiplier=1e6, subtractor=delta_F0, symbol="$\\delta f-\\delta f_\\mathrm{s}$"
78 )
79 mcmc.unit_dictionary["delta_F0"] = "$\\mu$Hz/s"
80 mcmc.transform_dictionary["tglitch"]["subtractor"] = tstart + dtglitch
81 mcmc.transform_dictionary["tglitch"][
82 "label"
83 ] = "$t^\\mathrm{g}-t^\\mathrm{g}_\\mathrm{s}$\n[d]"
84
85 t1 = time.time()
86 mcmc.run(save_loudest=False) # uses CFSv2 which doesn't support glitch parameters
87 dT = time.time() - t1
88 mcmc.print_summary()
89
90 print("Making corner plot...")
91 mcmc.plot_corner(
92 label_offset=0.25,
93 truths={"F0": F0, "F1": F1, "delta_F0": delta_F0, "tglitch": tstart + dtglitch},
94 quantiles=(0.16, 0.84),
95 hist_kwargs=dict(lw=1.5, zorder=-1),
96 truth_color="C3",
97 )
98
99 mcmc.plot_cumulative_max(savefig=True)
100
101 print(("Prior widths =", F0_width, F1_width))
102 print(("Actual run time = {}".format(dT)))
Total running time of the script: ( 0 minutes 0.000 seconds)