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