Note
Go to the end to download the full example code.
Long transient MCMC search
MCMC search for a long transient CW signal.
By default, the standard persistent-CW 2F-statistic and the transient max2F statistic are compared.
You can turn on either BSGL = True or BtSG = True (not both!) to test alternative statistics.
14 import os
15
16 import numpy as np
17 import PyFstat_example_make_data_for_long_transient_search as data
18
19 import pyfstat
20 from pyfstat.utils import get_predict_fstat_parameters_from_dict
21
22 if not os.path.isdir(data.outdir) or not np.any(
23 [f.endswith(".sft") for f in os.listdir(data.outdir)]
24 ):
25 raise RuntimeError(
26 "Please first run PyFstat_example_make_data_for_long_transient_search.py !"
27 )
28
29 label = "PyFstatExampleLongTransientMCMCSearch"
30 logger = pyfstat.set_up_logger(label=label, outdir=data.outdir)
31
32 tstart = data.tstart
33 duration = data.duration
34
35 inj = {
36 "tref": data.tstart,
37 "F0": data.F0,
38 "F1": data.F1,
39 "F2": data.F2,
40 "Alpha": data.Alpha,
41 "Delta": data.Delta,
42 "transient_tstart": data.transient_tstart,
43 "transient_duration": data.transient_duration,
44 }
45
46 DeltaF0 = 6e-7
47 DeltaF1 = 1e-13
48
49 # to make the search cheaper, we exactly target the transientStartTime
50 # to the injected value and only search over TransientTau
51 theta_prior = {
52 "F0": {
53 "type": "unif",
54 "lower": inj["F0"] - DeltaF0 / 2.0,
55 "upper": inj["F0"] + DeltaF0 / 2.0,
56 },
57 "F1": {
58 "type": "unif",
59 "lower": inj["F1"] - DeltaF1 / 2.0,
60 "upper": inj["F1"] + DeltaF1 / 2.0,
61 },
62 "F2": inj["F2"],
63 "Alpha": inj["Alpha"],
64 "Delta": inj["Delta"],
65 "transient_tstart": tstart + 0.25 * duration,
66 "transient_duration": {
67 "type": "halfnorm",
68 "loc": 0.001 * duration,
69 "scale": 0.5 * duration,
70 },
71 }
72
73 ntemps = 2
74 log10beta_min = -1
75 nwalkers = 100
76 nsteps = [100, 100]
77
78 transientWindowType = "rect"
79 BSGL = False
80 BtSG = False
81
82 mcmc = pyfstat.MCMCTransientSearch(
83 label=label + ("BSGL" if BSGL else "") + ("BtSG" if BtSG else ""),
84 outdir=data.outdir,
85 sftfilepattern=os.path.join(data.outdir, f"*{data.label}*sft"),
86 theta_prior=theta_prior,
87 tref=inj["tref"],
88 nsteps=nsteps,
89 nwalkers=nwalkers,
90 ntemps=ntemps,
91 log10beta_min=log10beta_min,
92 transientWindowType=transientWindowType,
93 BSGL=BSGL,
94 BtSG=BtSG,
95 )
96 mcmc.run(walker_plot_args={"plot_det_stat": True, "injection_parameters": inj})
97 mcmc.print_summary()
98 mcmc.plot_corner(add_prior=True, truths=inj)
99 mcmc.plot_prior_posterior(injection_parameters=inj)
100
101 # plot cumulative 2F, first building a dict as required for PredictFStat
102 d, maxtwoF = mcmc.get_max_twoF()
103 for key, val in mcmc.theta_prior.items():
104 if key not in d:
105 d[key] = val
106 d["h0"] = data.h0
107 d["cosi"] = data.cosi
108 d["psi"] = data.psi
109 PFS_input = get_predict_fstat_parameters_from_dict(
110 d, transientWindowType=transientWindowType
111 )
112 mcmc.plot_cumulative_max(PFS_input=PFS_input)