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

Total running time of the script: ( 0 minutes 0.000 seconds)

Gallery generated by Sphinx-Gallery