MCMC search: Semicoherent F-statisticΒΆ

Directed MCMC search for an isolated CW signal using the semicoherent F-statistic.

  8 import os
  9
 10 import numpy as np
 11
 12 import pyfstat
 13
 14 label = "PyFstatExampleSemiCoherentMCMCSearch"
 15 outdir = os.path.join("PyFstat_example_data", label)
 16 logger = pyfstat.set_up_logger(label=label, outdir=outdir)
 17
 18 # Properties of the GW data
 19 data_parameters = {
 20     "sqrtSX": 1e-23,
 21     "tstart": 1000000000,
 22     "duration": 100 * 86400,
 23     "detectors": "H1",
 24 }
 25 tend = data_parameters["tstart"] + data_parameters["duration"]
 26 mid_time = 0.5 * (data_parameters["tstart"] + tend)
 27
 28 # Properties of the signal
 29 depth = 10
 30 signal_parameters = {
 31     "F0": 30.0,
 32     "F1": -1e-10,
 33     "F2": 0,
 34     "Alpha": np.radians(83.6292),
 35     "Delta": np.radians(22.0144),
 36     "tref": mid_time,
 37     "h0": data_parameters["sqrtSX"] / depth,
 38     "cosi": 1.0,
 39 }
 40
 41 data = pyfstat.Writer(
 42     label=label, outdir=outdir, **data_parameters, **signal_parameters
 43 )
 44 data.make_data()
 45
 46 # The predicted twoF (expectation over noise realizations) can be accessed by
 47 twoF = data.predict_fstat()
 48 logger.info("Predicted twoF value: {}\n".format(twoF))
 49
 50 DeltaF0 = 1e-7
 51 DeltaF1 = 1e-13
 52 VF0 = (np.pi * data_parameters["duration"] * DeltaF0) ** 2 / 3.0
 53 VF1 = (np.pi * data_parameters["duration"] ** 2 * DeltaF1) ** 2 * 4 / 45.0
 54 logger.info("\nV={:1.2e}, VF0={:1.2e}, VF1={:1.2e}\n".format(VF0 * VF1, VF0, VF1))
 55
 56 theta_prior = {
 57     "F0": {
 58         "type": "unif",
 59         "lower": signal_parameters["F0"] - DeltaF0 / 2.0,
 60         "upper": signal_parameters["F0"] + DeltaF0 / 2.0,
 61     },
 62     "F1": {
 63         "type": "unif",
 64         "lower": signal_parameters["F1"] - DeltaF1 / 2.0,
 65         "upper": signal_parameters["F1"] + DeltaF1 / 2.0,
 66     },
 67 }
 68 for key in "F2", "Alpha", "Delta":
 69     theta_prior[key] = signal_parameters[key]
 70
 71 ntemps = 1
 72 log10beta_min = -1
 73 nwalkers = 100
 74 nsteps = [300, 300]
 75
 76 mcmc = pyfstat.MCMCSemiCoherentSearch(
 77     label=label,
 78     outdir=outdir,
 79     nsegs=10,
 80     sftfilepattern=data.sftfilepath,
 81     theta_prior=theta_prior,
 82     tref=mid_time,
 83     minStartTime=data_parameters["tstart"],
 84     maxStartTime=tend,
 85     nsteps=nsteps,
 86     nwalkers=nwalkers,
 87     ntemps=ntemps,
 88     log10beta_min=log10beta_min,
 89 )
 90 mcmc.transform_dictionary = dict(
 91     F0=dict(subtractor=signal_parameters["F0"], symbol="$f-f^\\mathrm{s}$"),
 92     F1=dict(
 93         subtractor=signal_parameters["F1"], symbol="$\\dot{f}-\\dot{f}^\\mathrm{s}$"
 94     ),
 95 )
 96 mcmc.run(
 97     walker_plot_args={"plot_det_stat": True, "injection_parameters": signal_parameters}
 98 )
 99 mcmc.print_summary()
100 mcmc.plot_corner(add_prior=True, truths=signal_parameters)
101 mcmc.plot_prior_posterior(injection_parameters=signal_parameters)
102 mcmc.plot_chainconsumer(truth=signal_parameters)
103 mcmc.plot_cumulative_max(
104     savefig=True,
105     custom_ax_kwargs={"title": "Cumulative 2F for the best MCMC candidate"},
106 )

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

Gallery generated by Sphinx-Gallery