Note
Click here to download the full example code
MCMC search: Semicoherent F-statistic
Directed MCMC search for an isolated CW signal using the semicoherent F-statistic.
8 import pyfstat
9 import numpy as np
10 import os
11
12 label = "PyFstat_example_semi_coherent_MCMC_search"
13 outdir = os.path.join("PyFstat_example_data", label)
14
15 # Properties of the GW data
16 data_parameters = {
17 "sqrtSX": 1e-23,
18 "tstart": 1000000000,
19 "duration": 100 * 86400,
20 "detectors": "H1",
21 }
22 tend = data_parameters["tstart"] + data_parameters["duration"]
23 mid_time = 0.5 * (data_parameters["tstart"] + tend)
24
25 # Properties of the signal
26 depth = 10
27 signal_parameters = {
28 "F0": 30.0,
29 "F1": -1e-10,
30 "F2": 0,
31 "Alpha": np.radians(83.6292),
32 "Delta": np.radians(22.0144),
33 "tref": mid_time,
34 "h0": data_parameters["sqrtSX"] / depth,
35 "cosi": 1.0,
36 }
37
38 data = pyfstat.Writer(
39 label=label, outdir=outdir, **data_parameters, **signal_parameters
40 )
41 data.make_data()
42
43 # The predicted twoF, given by lalapps_predictFstat can be accessed by
44 twoF = data.predict_fstat()
45 print("Predicted twoF value: {}\n".format(twoF))
46
47 DeltaF0 = 1e-7
48 DeltaF1 = 1e-13
49 VF0 = (np.pi * data_parameters["duration"] * DeltaF0) ** 2 / 3.0
50 VF1 = (np.pi * data_parameters["duration"] ** 2 * DeltaF1) ** 2 * 4 / 45.0
51 print("\nV={:1.2e}, VF0={:1.2e}, VF1={:1.2e}\n".format(VF0 * VF1, VF0, VF1))
52
53 theta_prior = {
54 "F0": {
55 "type": "unif",
56 "lower": signal_parameters["F0"] - DeltaF0 / 2.0,
57 "upper": signal_parameters["F0"] + DeltaF0 / 2.0,
58 },
59 "F1": {
60 "type": "unif",
61 "lower": signal_parameters["F1"] - DeltaF1 / 2.0,
62 "upper": signal_parameters["F1"] + DeltaF1 / 2.0,
63 },
64 }
65 for key in "F2", "Alpha", "Delta":
66 theta_prior[key] = signal_parameters[key]
67
68 ntemps = 1
69 log10beta_min = -1
70 nwalkers = 100
71 nsteps = [300, 300]
72
73 mcmc = pyfstat.MCMCSemiCoherentSearch(
74 label=label,
75 outdir=outdir,
76 nsegs=10,
77 sftfilepattern=os.path.join(outdir, "*{}*sft".format(label)),
78 theta_prior=theta_prior,
79 tref=mid_time,
80 minStartTime=data_parameters["tstart"],
81 maxStartTime=tend,
82 nsteps=nsteps,
83 nwalkers=nwalkers,
84 ntemps=ntemps,
85 log10beta_min=log10beta_min,
86 )
87 mcmc.transform_dictionary = dict(
88 F0=dict(subtractor=signal_parameters["F0"], symbol="$f-f^\\mathrm{s}$"),
89 F1=dict(
90 subtractor=signal_parameters["F1"], symbol="$\\dot{f}-\\dot{f}^\\mathrm{s}$"
91 ),
92 )
93 mcmc.run(
94 walker_plot_args={"plot_det_stat": True, "injection_parameters": signal_parameters}
95 )
96 mcmc.print_summary()
97 mcmc.plot_corner(add_prior=True, truths=signal_parameters)
98 mcmc.plot_prior_posterior(injection_parameters=signal_parameters)
99 mcmc.plot_chainconsumer(truth=signal_parameters)
100 mcmc.plot_cumulative_max(
101 savefig=True,
102 custom_ax_kwargs={"title": "Cumulative 2F for the best MCMC candidate"},
103 )
Total running time of the script: ( 0 minutes 0.000 seconds)