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