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