Note
Click here to download the full example code
MCMC search with fully coherent BSGL statistic
Targeted MCMC search for an isolated CW signal using the fully coherent line-robust BSGL-statistic.
9 import pyfstat
10 import numpy as np
11 import os
12
13 label = os.path.splitext(os.path.basename(__file__))[0]
14 outdir = os.path.join("PyFstat_example_data", label)
15
16 # Properties of the GW data - first we make data for two detectors,
17 # both including Gaussian noise and a coherent 'astrophysical' signal.
18 data_parameters = {
19 "sqrtSX": 1e-23,
20 "tstart": 1000000000,
21 "duration": 100 * 86400,
22 "detectors": "H1,L1",
23 "SFTWindowType": "tukey",
24 "SFTWindowBeta": 0.001,
25 }
26 tend = data_parameters["tstart"] + data_parameters["duration"]
27 mid_time = 0.5 * (data_parameters["tstart"] + tend)
28
29 # Properties of the signal
30 depth = 10
31 signal_parameters = {
32 "F0": 30.0,
33 "F1": -1e-10,
34 "F2": 0,
35 "Alpha": np.radians(83.6292),
36 "Delta": np.radians(22.0144),
37 "tref": mid_time,
38 "h0": data_parameters["sqrtSX"] / depth,
39 "cosi": 1.0,
40 }
41
42 data = pyfstat.Writer(
43 label=label, outdir=outdir, **data_parameters, **signal_parameters
44 )
45 data.make_data()
46
47 # Now we add an additional single-detector artifact to H1 only.
48 # For simplicity, this is modelled here as a fully modulated CW-like signal,
49 # just restricted to the single detector.
50 SFTs_H1 = data.sftfilepath.split(";")[0]
51 data_parameters_line = data_parameters.copy()
52 signal_parameters_line = signal_parameters.copy()
53 data_parameters_line["detectors"] = "H1"
54 data_parameters_line["sqrtSX"] = 0 # don't add yet another set of Gaussian noise
55 signal_parameters_line["F0"] += 1e-6
56 signal_parameters_line["h0"] *= 10.0
57 extra_writer = pyfstat.Writer(
58 label=label,
59 outdir=outdir,
60 **data_parameters_line,
61 **signal_parameters_line,
62 noiseSFTs=SFTs_H1,
63 )
64 extra_writer.make_data()
65
66 # The predicted twoF, given by lalapps_predictFstat can be accessed by
67 twoF = data.predict_fstat()
68 print("Predicted twoF value: {}\n".format(twoF))
69
70 # MCMC prior ranges
71 DeltaF0 = 1e-5
72 DeltaF1 = 1e-13
73 theta_prior = {
74 "F0": {
75 "type": "unif",
76 "lower": signal_parameters["F0"] - DeltaF0 / 2.0,
77 "upper": signal_parameters["F0"] + DeltaF0 / 2.0,
78 },
79 "F1": {
80 "type": "unif",
81 "lower": signal_parameters["F1"] - DeltaF1 / 2.0,
82 "upper": signal_parameters["F1"] + DeltaF1 / 2.0,
83 },
84 }
85 for key in "F2", "Alpha", "Delta":
86 theta_prior[key] = signal_parameters[key]
87
88 # MCMC sampler settings - relatively cheap setup, may not converge perfectly
89 ntemps = 2
90 log10beta_min = -0.5
91 nwalkers = 50
92 nsteps = [100, 100]
93
94 # we'll want to plot results relative to the injection parameters
95 transform_dict = dict(
96 F0=dict(subtractor=signal_parameters["F0"], symbol="$f-f^\\mathrm{s}$"),
97 F1=dict(
98 subtractor=signal_parameters["F1"], symbol="$\\dot{f}-\\dot{f}^\\mathrm{s}$"
99 ),
100 )
101
102 # first search: standard F-statistic
103 # This should show a weak peak from the coherent signal
104 # and a larger one from the "line artifact" at higher frequency.
105 mcmc_F = pyfstat.MCMCSearch(
106 label=label + "_twoF",
107 outdir=outdir,
108 sftfilepattern=os.path.join(outdir, "*{}*sft".format(label)),
109 theta_prior=theta_prior,
110 tref=mid_time,
111 minStartTime=data_parameters["tstart"],
112 maxStartTime=tend,
113 nsteps=nsteps,
114 nwalkers=nwalkers,
115 ntemps=ntemps,
116 log10beta_min=log10beta_min,
117 BSGL=False,
118 )
119 mcmc_F.transform_dictionary = transform_dict
120 mcmc_F.run(
121 walker_plot_args={"plot_det_stat": True, "injection_parameters": signal_parameters}
122 )
123 mcmc_F.print_summary()
124 mcmc_F.plot_corner(add_prior=True, truths=signal_parameters)
125 mcmc_F.plot_prior_posterior(injection_parameters=signal_parameters)
126
127 # second search: line-robust statistic BSGL activated
128 mcmc_F = pyfstat.MCMCSearch(
129 label=label + "_BSGL",
130 outdir=outdir,
131 sftfilepattern=os.path.join(outdir, "*{}*sft".format(label)),
132 theta_prior=theta_prior,
133 tref=mid_time,
134 minStartTime=data_parameters["tstart"],
135 maxStartTime=tend,
136 nsteps=nsteps,
137 nwalkers=nwalkers,
138 ntemps=ntemps,
139 log10beta_min=log10beta_min,
140 BSGL=True,
141 )
142 mcmc_F.transform_dictionary = transform_dict
143 mcmc_F.run(
144 walker_plot_args={"plot_det_stat": True, "injection_parameters": signal_parameters}
145 )
146 mcmc_F.print_summary()
147 mcmc_F.plot_corner(add_prior=True, truths=signal_parameters)
148 mcmc_F.plot_prior_posterior(injection_parameters=signal_parameters)
Total running time of the script: ( 0 minutes 0.000 seconds)