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