Note
Go to the end to download the full example code.
MCMC search: Fully coherent F-statistic
Directed MCMC search for an isolated CW signal using the fully coherent F-statistic.
9 import os
10
11 import numpy as np
12
13 import pyfstat
14 from pyfstat.utils import get_predict_fstat_parameters_from_dict
15
16 label = "PyFstatExampleFullyCoherentMCMCSearch"
17 outdir = os.path.join("PyFstat_example_data", label)
18 logger = pyfstat.set_up_logger(label=label, outdir=outdir)
19
20 # Properties of the GW data
21 data_parameters = {
22 "sqrtSX": 1e-23,
23 "tstart": 1000000000,
24 "duration": 100 * 86400,
25 "detectors": "H1",
26 }
27 tend = data_parameters["tstart"] + data_parameters["duration"]
28 mid_time = 0.5 * (data_parameters["tstart"] + tend)
29
30 # Properties of the signal
31 depth = 10
32 signal_parameters = {
33 "F0": 30.0,
34 "F1": -1e-10,
35 "F2": 0,
36 "F3": 0,
37 "Alpha": np.radians(83.6292),
38 "Delta": np.radians(22.0144),
39 "tref": mid_time,
40 "h0": data_parameters["sqrtSX"] / depth,
41 "cosi": 1.0,
42 }
43
44 data = pyfstat.Writer(
45 label=label,
46 outdir=outdir,
47 **data_parameters,
48 signal_parameters=signal_parameters,
49 )
50 data.make_data()
51
52 # The predicted twoF (expectation over noise realizations) can be accessed by
53 twoF = data.predict_fstat()
54 logger.info("Predicted twoF value: {}\n".format(twoF))
55
56 DeltaF0 = 1e-7
57 DeltaF1 = 1e-13
58 VF0 = (np.pi * data_parameters["duration"] * DeltaF0) ** 2 / 3.0
59 VF1 = (np.pi * data_parameters["duration"] ** 2 * DeltaF1) ** 2 * 4 / 45.0
60 logger.info("\nV={:1.2e}, VF0={:1.2e}, VF1={:1.2e}\n".format(VF0 * VF1, VF0, VF1))
61
62 theta_prior = {
63 "F0": {
64 "type": "unif",
65 "lower": signal_parameters["F0"] - DeltaF0 / 2.0,
66 "upper": signal_parameters["F0"] + DeltaF0 / 2.0,
67 },
68 "F1": {
69 "type": "unif",
70 "lower": signal_parameters["F1"] - DeltaF1 / 2.0,
71 "upper": signal_parameters["F1"] + DeltaF1 / 2.0,
72 },
73 }
74 for key in "F2", "F3", "Alpha", "Delta":
75 theta_prior[key] = signal_parameters[key]
76
77 ntemps = 2
78 log10beta_min = -0.5
79 nwalkers = 100
80 nsteps = [100, 100]
81
82 mcmc = pyfstat.MCMCSearch(
83 label=label,
84 outdir=outdir,
85 sftfilepattern=data.sftfilepath,
86 theta_prior=theta_prior,
87 tref=mid_time,
88 minStartTime=data_parameters["tstart"],
89 maxStartTime=tend,
90 nsteps=nsteps,
91 nwalkers=nwalkers,
92 ntemps=ntemps,
93 log10beta_min=log10beta_min,
94 )
95 mcmc.transform_dictionary = dict(
96 F0=dict(subtractor=signal_parameters["F0"], symbol="$f-f^\\mathrm{s}$"),
97 F1=dict(
98 subtractor=signal_parameters["F1"], symbol="$\\dot{f}-\\dot{f}^\\mathrm{s}$"
99 ),
100 )
101 mcmc.run(
102 walker_plot_args={"plot_det_stat": True, "injection_parameters": signal_parameters}
103 )
104 mcmc.print_summary()
105 mcmc.plot_corner(add_prior=True, truths=signal_parameters)
106 mcmc.plot_prior_posterior(injection_parameters=signal_parameters)
107
108 mcmc.generate_loudest()
109
110 # plot cumulative 2F, first building a dict as required for PredictFStat
111 d, maxtwoF = mcmc.get_max_twoF()
112 for key, val in mcmc.theta_prior.items():
113 if key not in d:
114 d[key] = val
115 d["h0"] = data.h0
116 d["cosi"] = data.cosi
117 d["psi"] = data.psi
118 PFS_input = get_predict_fstat_parameters_from_dict(d)
119 mcmc.plot_cumulative_max(PFS_input=PFS_input)