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