Note

Go to the end 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)

Gallery generated by Sphinx-Gallery