Note
Click here to download the full example code
Binary CW example: Comparison between MCMC and grid search
Comparison of the semicoherent F-statistic MCMC search algorithm to a simple grid search accross the parameter space corresponding to a CW source in a binary system.
10 import pyfstat
11 import os
12 import numpy as np
13 import matplotlib.pyplot as plt
14
15 # Set to false to include eccentricity
16 circular_orbit = False
17
18 label = "PyFstat_example_binary_mcmc_vs_grid_comparison" + (
19 "_circular_orbit" if circular_orbit else ""
20 )
21 outdir = os.path.join("PyFstat_example_data", label)
22
23 # Parameters to generate a data set
24 data_parameters = {
25 "sqrtSX": 1e-22,
26 "tstart": 1000000000,
27 "duration": 90 * 86400,
28 "detectors": "H1,L1",
29 "Tsft": 3600,
30 "Band": 4,
31 }
32
33 # Injected signal parameters
34 tend = data_parameters["tstart"] + data_parameters["duration"]
35 mid_time = 0.5 * (data_parameters["tstart"] + tend)
36 depth = 10.0
37 signal_parameters = {
38 "tref": data_parameters["tstart"],
39 "F0": 40.0,
40 "F1": 0,
41 "F2": 0,
42 "Alpha": 0.5,
43 "Delta": 0.5,
44 "period": 85 * 24 * 3600.0,
45 "asini": 4.0,
46 "tp": mid_time * 1.05,
47 "argp": 0.0 if circular_orbit else 0.54,
48 "ecc": 0.0 if circular_orbit else 0.7,
49 "h0": data_parameters["sqrtSX"] / depth,
50 "cosi": 1.0,
51 }
52
53
54 print("Generating SFTs with injected signal...")
55 writer = pyfstat.BinaryModulatedWriter(
56 label="simulated_signal",
57 outdir=outdir,
58 **data_parameters,
59 **signal_parameters,
60 )
61 writer.make_data()
62 print("")
63
64 print("Performing Grid Search...")
65
66 # Create ad-hoc grid and compute Fstatistic around injection point
67 # There's no class supporting a binary search in the same way as
68 # grid_based_searches.GridSearch, so we do it by hand constructing
69 # a grid and using core.ComputeFstat.
70 half_points_per_dimension = 2
71 search_keys = ["period", "asini", "tp", "argp", "ecc"]
72 search_keys_label = [
73 r"$P$ [s]",
74 r"$a_p$ [s]",
75 r"$t_{p}$ [s]",
76 r"$\omega$ [rad]",
77 r"$e$",
78 ]
79
80 grid_arrays = np.meshgrid(
81 *[
82 signal_parameters[key]
83 * (
84 1
85 + 0.01
86 * np.arange(-half_points_per_dimension, half_points_per_dimension + 1, 1)
87 )
88 for key in search_keys
89 ]
90 )
91 grid_points = np.hstack(
92 [grid_arrays[i].reshape(-1, 1) for i in range(len(grid_arrays))]
93 )
94
95 compute_f_stat = pyfstat.ComputeFstat(
96 sftfilepattern=os.path.join(outdir, "*simulated_signal*sft"),
97 tref=signal_parameters["tref"],
98 binary=True,
99 minCoverFreq=-0.5,
100 maxCoverFreq=-0.5,
101 )
102 twoF_values = np.zeros(grid_points.shape[0])
103 for ind in range(grid_points.shape[0]):
104 point = grid_points[ind]
105 twoF_values[ind] = compute_f_stat.get_fullycoherent_twoF(
106 F0=signal_parameters["F0"],
107 F1=signal_parameters["F1"],
108 F2=signal_parameters["F2"],
109 Alpha=signal_parameters["Alpha"],
110 Delta=signal_parameters["Delta"],
111 period=point[0],
112 asini=point[1],
113 tp=point[2],
114 argp=point[3],
115 ecc=point[4],
116 )
117 print(f"2Fstat computed on {grid_points.shape[0]} points")
118 print("")
119 print("Plotting results...")
120 dim = len(search_keys)
121 fig, ax = plt.subplots(dim, 1, figsize=(10, 10))
122 for ind in range(dim):
123 a = ax.ravel()[ind]
124 a.grid()
125 a.set(xlabel=search_keys_label[ind], ylabel=r"$2 \mathcal{F}$", yscale="log")
126 a.plot(grid_points[:, ind], twoF_values, "o")
127 a.axvline(signal_parameters[search_keys[ind]], label="Injection", color="orange")
128 plt.tight_layout()
129 fig.savefig(os.path.join(outdir, "grid_twoF_per_dimension.png"))
130
131
132 print("Performing MCMCSearch...")
133 # Fixed points in frequency and sky parameters
134 theta_prior = {
135 "F0": signal_parameters["F0"],
136 "F1": signal_parameters["F1"],
137 "F2": signal_parameters["F2"],
138 "Alpha": signal_parameters["Alpha"],
139 "Delta": signal_parameters["Delta"],
140 }
141
142 # Set up priors for the binary parameters
143 for key in search_keys:
144 theta_prior.update(
145 {
146 key: {
147 "type": "unif",
148 "lower": 0.999 * signal_parameters[key],
149 "upper": 1.001 * signal_parameters[key],
150 }
151 }
152 )
153 if circular_orbit:
154 for key in ["ecc", "argp"]:
155 theta_prior[key] = 0
156 search_keys.remove(key)
157
158 # ptemcee sampler settings - in a real application we might want higher values
159 ntemps = 2
160 log10beta_min = -1
161 nwalkers = 100
162 nsteps = [100, 100] # [burnin,production]
163
164 mcmcsearch = pyfstat.MCMCSearch(
165 label="mcmc_search",
166 outdir=outdir,
167 sftfilepattern=os.path.join(outdir, "*simulated_signal*sft"),
168 theta_prior=theta_prior,
169 tref=signal_parameters["tref"],
170 nsteps=nsteps,
171 nwalkers=nwalkers,
172 ntemps=ntemps,
173 log10beta_min=log10beta_min,
174 binary=True,
175 )
176 # walker plot is generated during main run of the search class
177 mcmcsearch.run(
178 plot_walkers=True,
179 walker_plot_args={"plot_det_stat": True, "injection_parameters": signal_parameters},
180 )
181 mcmcsearch.print_summary()
182
183 # call some built-in plotting methods
184 # these can all highlight the injection parameters, too
185 print("Making MCMCSearch {:s} corner plot...".format("-".join(search_keys)))
186 mcmcsearch.plot_corner(truths=signal_parameters)
187 print("Making MCMCSearch prior-posterior comparison plot...")
188 mcmcsearch.plot_prior_posterior(injection_parameters=signal_parameters)
189 print("")
190
191 print("*" * 20)
192 print("Quantitative comparisons:")
193 print("*" * 20)
194
195 # some informative command-line output comparing search results and injection
196 # get max twoF and binary Doppler parameters
197 max_grid_index = np.argmax(twoF_values)
198 max_grid_2F = twoF_values[max_grid_index]
199 max_grid_parameters = grid_points[max_grid_index]
200
201 # same for MCMCSearch, here twoF is separate, and non-sampled parameters are not included either
202 max_dict_mcmc, max_2F_mcmc = mcmcsearch.get_max_twoF()
203 print(
204 "Grid Search:\n\tmax2F={:.4f}\n\tOffsets from injection parameters (relative error): {:s}.".format(
205 max_grid_2F,
206 ", ".join(
207 [
208 "\n\t\t{1:s}: {0:.4e} ({2:.4f}%)".format(
209 max_grid_parameters[search_keys.index(key)]
210 - signal_parameters[key],
211 key,
212 100
213 * (
214 max_grid_parameters[search_keys.index(key)]
215 - signal_parameters[key]
216 )
217 / signal_parameters[key],
218 )
219 for key in search_keys
220 ]
221 ),
222 )
223 )
224 print(
225 "Max 2F candidate from MCMC Search:\n\tmax2F={:.4f}"
226 "\n\tOffsets from injection parameters (relative error): {:s}.".format(
227 max_2F_mcmc,
228 ", ".join(
229 [
230 "\n\t\t{1:s}: {0:.4e} ({2:.4f}%)".format(
231 max_dict_mcmc[key] - signal_parameters[key],
232 key,
233 100
234 * (max_dict_mcmc[key] - signal_parameters[key])
235 / signal_parameters[key],
236 )
237 for key in search_keys
238 ]
239 ),
240 )
241 )
242 # get additional point and interval estimators
243 stats_dict_mcmc = mcmcsearch.get_summary_stats()
244 print(
245 "Mean from MCMCSearch:\n\tOffset from injection parameters (relative error): {:s}"
246 "\n\tExpressed as fractions of 2sigma intervals: {:s}.".format(
247 ", ".join(
248 [
249 "\n\t\t{1:s}: {0:.4e} ({2:.4f}%)".format(
250 stats_dict_mcmc[key]["mean"] - signal_parameters[key],
251 key,
252 100
253 * (stats_dict_mcmc[key]["mean"] - signal_parameters[key])
254 / signal_parameters[key],
255 )
256 for key in search_keys
257 ]
258 ),
259 ", ".join(
260 [
261 "\n\t\t{1:s}: {0:.4f}%".format(
262 100
263 * np.abs(stats_dict_mcmc[key]["mean"] - signal_parameters[key])
264 / (2 * stats_dict_mcmc[key]["std"]),
265 key,
266 )
267 for key in search_keys
268 ]
269 ),
270 )
271 )
272 print(
273 "Median from MCMCSearch:\n\tOffset from injection parameters (relative error): {:s},"
274 "\n\tExpressed as fractions of 90% confidence intervals: {:s}.".format(
275 ", ".join(
276 [
277 "\n\t\t{1:s}: {0:.4e} ({2:.4f}%)".format(
278 stats_dict_mcmc[key]["median"] - signal_parameters[key],
279 key,
280 100
281 * (stats_dict_mcmc[key]["median"] - signal_parameters[key])
282 / signal_parameters[key],
283 )
284 for key in search_keys
285 ]
286 ),
287 ", ".join(
288 [
289 "\n\t\t{1:s}: {0:.4f}%".format(
290 100
291 * np.abs(stats_dict_mcmc[key]["median"] - signal_parameters[key])
292 / (
293 stats_dict_mcmc[key]["upper90"]
294 - stats_dict_mcmc[key]["lower90"]
295 ),
296 key,
297 )
298 for key in search_keys
299 ]
300 ),
301 )
302 )
Total running time of the script: ( 0 minutes 0.000 seconds)