Note
Go to the end to download the full example code.
MCMC search: Semicoherent F-statistic with initialisation
Directed MCMC search for an isolated CW signal using the fully-coherent F-statistic. Prior to the burn-in stage, walkers are initialized with a certain scattering factor.
10 import os
11
12 import numpy as np
13
14 import pyfstat
15
16 label = "PyFstatExampleMCMCSearchUsingInitialisation"
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 "Alpha": np.radians(83.6292),
37 "Delta": np.radians(22.0144),
38 "tref": mid_time,
39 "h0": data_parameters["sqrtSX"] / depth,
40 "cosi": 1.0,
41 }
42
43 data = pyfstat.Writer(
44 label=label, outdir=outdir, **data_parameters, **signal_parameters
45 )
46 data.make_data()
47
48 # The predicted twoF (expectation over noise realizations) can be accessed by
49 twoF = data.predict_fstat()
50 logger.info("Predicted twoF value: {}\n".format(twoF))
51
52 DeltaF0 = 1e-7
53 DeltaF1 = 1e-13
54 VF0 = (np.pi * data_parameters["duration"] * DeltaF0) ** 2 / 3.0
55 VF1 = (np.pi * data_parameters["duration"] ** 2 * DeltaF1) ** 2 * 4 / 45.0
56 logger.info("\nV={:1.2e}, VF0={:1.2e}, VF1={:1.2e}\n".format(VF0 * VF1, VF0, VF1))
57
58 theta_prior = {
59 "F0": {
60 "type": "unif",
61 "lower": signal_parameters["F0"] - DeltaF0 / 2.0,
62 "upper": signal_parameters["F0"] + DeltaF0 / 2.0,
63 },
64 "F1": {
65 "type": "unif",
66 "lower": signal_parameters["F1"] - DeltaF1 / 2.0,
67 "upper": signal_parameters["F1"] + DeltaF1 / 2.0,
68 },
69 }
70 for key in "F2", "Alpha", "Delta":
71 theta_prior[key] = signal_parameters[key]
72
73 ntemps = 1
74 log10beta_min = -1
75 nwalkers = 100
76 nsteps = [100, 100]
77
78 mcmc = pyfstat.MCMCSearch(
79 label=label,
80 outdir=outdir,
81 sftfilepattern=data.sftfilepath,
82 theta_prior=theta_prior,
83 tref=mid_time,
84 minStartTime=data_parameters["tstart"],
85 maxStartTime=tend,
86 nsteps=nsteps,
87 nwalkers=nwalkers,
88 ntemps=ntemps,
89 log10beta_min=log10beta_min,
90 )
91 mcmc.setup_initialisation(100, scatter_val=1e-10)
92 mcmc.run(
93 walker_plot_args={"plot_det_stat": True, "injection_parameters": signal_parameters}
94 )
95 mcmc.print_summary()
96 mcmc.plot_corner(add_prior=True, truths=signal_parameters)
97 mcmc.plot_prior_posterior(injection_parameters=signal_parameters)