Note
Click here to download the full example code
Follow up example
Multi-stage MCMC follow up of a CW signal produced by an isolated source using a ladder of coherent times.
8 import pyfstat
9 import numpy as np
10 import matplotlib.pyplot as plt
11 import os
12
13 label = "PyFstat_example_semi_coherent_directed_follow_up"
14 outdir = os.path.join("PyFstat_example_data", label)
15
16 # Properties of the GW data
17 data_parameters = {
18 "sqrtSX": 1e-23,
19 "tstart": 1000000000,
20 "duration": 100 * 86400,
21 "detectors": "H1",
22 }
23 tend = data_parameters["tstart"] + data_parameters["duration"]
24 mid_time = 0.5 * (data_parameters["tstart"] + tend)
25
26 # Properties of the signal
27 depth = 40
28 signal_parameters = {
29 "F0": 30.0,
30 "F1": -1e-10,
31 "F2": 0,
32 "Alpha": np.radians(83.6292),
33 "Delta": np.radians(22.0144),
34 "tref": mid_time,
35 "h0": data_parameters["sqrtSX"] / depth,
36 "cosi": 1.0,
37 }
38
39 data = pyfstat.Writer(
40 label=label, outdir=outdir, **data_parameters, **signal_parameters
41 )
42 data.make_data()
43
44 # The predicted twoF, given by lalapps_predictFstat can be accessed by
45 twoF = data.predict_fstat()
46 print("Predicted twoF value: {}\n".format(twoF))
47
48 # Search
49 VF0 = VF1 = 1e5
50 DeltaF0 = np.sqrt(VF0) * np.sqrt(3) / (np.pi * data_parameters["duration"])
51 DeltaF1 = np.sqrt(VF1) * np.sqrt(180) / (np.pi * data_parameters["duration"] ** 2)
52 theta_prior = {
53 "F0": {
54 "type": "unif",
55 "lower": signal_parameters["F0"] - DeltaF0 / 2.0,
56 "upper": signal_parameters["F0"] + DeltaF0 / 2,
57 },
58 "F1": {
59 "type": "unif",
60 "lower": signal_parameters["F1"] - DeltaF1 / 2.0,
61 "upper": signal_parameters["F1"] + DeltaF1 / 2,
62 },
63 }
64 for key in "F2", "Alpha", "Delta":
65 theta_prior[key] = signal_parameters[key]
66
67
68 ntemps = 3
69 log10beta_min = -0.5
70 nwalkers = 100
71 nsteps = [100, 100]
72
73 mcmc = pyfstat.MCMCFollowUpSearch(
74 label=label,
75 outdir=outdir,
76 sftfilepattern=os.path.join(outdir, "*{}*sft".format(label)),
77 theta_prior=theta_prior,
78 tref=mid_time,
79 minStartTime=data_parameters["tstart"],
80 maxStartTime=tend,
81 nwalkers=nwalkers,
82 nsteps=nsteps,
83 ntemps=ntemps,
84 log10beta_min=log10beta_min,
85 )
86
87 NstarMax = 1000
88 Nsegs0 = 100
89 walkers_fig, walkers_axes = plt.subplots(nrows=2, figsize=(3.4, 3.5))
90 mcmc.run(
91 NstarMax=NstarMax,
92 Nsegs0=Nsegs0,
93 plot_walkers=True,
94 walker_plot_args={
95 "labelpad": 0.01,
96 "plot_det_stat": False,
97 "fig": walkers_fig,
98 "axes": walkers_axes,
99 "injection_parameters": signal_parameters,
100 },
101 )
102 walkers_fig.savefig(os.path.join(outdir, label + "_walkers.png"))
103 plt.close(walkers_fig)
104
105 mcmc.print_summary()
106 mcmc.plot_corner(add_prior=True, truths=signal_parameters)
107 mcmc.plot_prior_posterior(injection_parameters=signal_parameters)
Total running time of the script: ( 0 minutes 0.000 seconds)