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