Note

Click here to download the full example code

MCMC search: Semicoherent F-statisticΒΆ

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

  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
 import pyfstat
 import numpy as np
 import os

 label = "PyFstat_example_semi_coherent_MCMC_search"
 outdir = os.path.join("PyFstat_example_data", label)

 # Properties of the GW data
 data_parameters = {
     "sqrtSX": 1e-23,
     "tstart": 1000000000,
     "duration": 100 * 86400,
     "detectors": "H1",
 }
 tend = data_parameters["tstart"] + data_parameters["duration"]
 mid_time = 0.5 * (data_parameters["tstart"] + tend)

 # Properties of the signal
 depth = 10
 signal_parameters = {
     "F0": 30.0,
     "F1": -1e-10,
     "F2": 0,
     "Alpha": np.radians(83.6292),
     "Delta": np.radians(22.0144),
     "tref": mid_time,
     "h0": data_parameters["sqrtSX"] / depth,
     "cosi": 1.0,
 }

 data = pyfstat.Writer(
     label=label, outdir=outdir, **data_parameters, **signal_parameters
 )
 data.make_data()

 # The predicted twoF, given by lalapps_predictFstat can be accessed by
 twoF = data.predict_fstat()
 print("Predicted twoF value: {}\n".format(twoF))

 DeltaF0 = 1e-7
 DeltaF1 = 1e-13
 VF0 = (np.pi * data_parameters["duration"] * DeltaF0) ** 2 / 3.0
 VF1 = (np.pi * data_parameters["duration"] ** 2 * DeltaF1) ** 2 * 4 / 45.0
 print("\nV={:1.2e}, VF0={:1.2e}, VF1={:1.2e}\n".format(VF0 * VF1, VF0, VF1))

 theta_prior = {
     "F0": {
         "type": "unif",
         "lower": signal_parameters["F0"] - DeltaF0 / 2.0,
         "upper": signal_parameters["F0"] + DeltaF0 / 2.0,
     },
     "F1": {
         "type": "unif",
         "lower": signal_parameters["F1"] - DeltaF1 / 2.0,
         "upper": signal_parameters["F1"] + DeltaF1 / 2.0,
     },
 }
 for key in "F2", "Alpha", "Delta":
     theta_prior[key] = signal_parameters[key]

 ntemps = 1
 log10beta_min = -1
 nwalkers = 100
 nsteps = [300, 300]

 mcmc = pyfstat.MCMCSemiCoherentSearch(
     label=label,
     outdir=outdir,
     nsegs=10,
     sftfilepattern=os.path.join(outdir, "*{}*sft".format(label)),
     theta_prior=theta_prior,
     tref=mid_time,
     minStartTime=data_parameters["tstart"],
     maxStartTime=tend,
     nsteps=nsteps,
     nwalkers=nwalkers,
     ntemps=ntemps,
     log10beta_min=log10beta_min,
 )
 mcmc.transform_dictionary = dict(
     F0=dict(subtractor=signal_parameters["F0"], symbol="$f-f^\\mathrm{s}$"),
     F1=dict(
         subtractor=signal_parameters["F1"], symbol="$\\dot{f}-\\dot{f}^\\mathrm{s}$"
     ),
 )
 mcmc.run(
     walker_plot_args={"plot_det_stat": True, "injection_parameters": signal_parameters}
 )
 mcmc.print_summary()
 mcmc.plot_corner(add_prior=True, truths=signal_parameters)
 mcmc.plot_prior_posterior(injection_parameters=signal_parameters)
 mcmc.plot_chainconsumer(truth=signal_parameters)
 mcmc.plot_cumulative_max(
     savefig=True,
     custom_ax_kwargs={"title": "Cumulative 2F for the best MCMC candidate"},
 )

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

Gallery generated by Sphinx-Gallery