Note
Go to the end to download the full example code.
Binary CW example: Semicoherent MCMC search
MCMC search of a CW signal produced by a source in a binary system using the semicoherent F-statistic.
9 import os
10
11 import numpy as np
12
13 import pyfstat
14
15 # If False, sky priors are used
16 directed_search = True
17 # If False, ecc and argp priors are used
18 known_eccentricity = True
19
20 label = "PyFstatExampleSemiCoherentBinarySearchUsingMCMC"
21 outdir = os.path.join("PyFstat_example_data", label)
22 logger = pyfstat.set_up_logger(label=label, outdir=outdir)
23
24 # Properties of the GW data
25 data_parameters = {
26 "sqrtSX": 1e-23,
27 "tstart": 1000000000,
28 "duration": 10 * 86400,
29 "detectors": "H1",
30 }
31 tend = data_parameters["tstart"] + data_parameters["duration"]
32 mid_time = 0.5 * (data_parameters["tstart"] + tend)
33
34 # Properties of the signal
35 depth = 0.1
36 signal_parameters = {
37 "F0": 30.0,
38 "F1": 0,
39 "F2": 0,
40 "Alpha": 0.15,
41 "Delta": 0.45,
42 "tp": mid_time,
43 "argp": 0.3,
44 "asini": 10.0,
45 "ecc": 0.1,
46 "period": 45 * 24 * 3600.0,
47 "tref": mid_time,
48 "h0": data_parameters["sqrtSX"] / depth,
49 "cosi": 1.0,
50 }
51
52 data = pyfstat.BinaryModulatedWriter(
53 label=label, outdir=outdir, **data_parameters, **signal_parameters
54 )
55 data.make_data()
56
57 theta_prior = {
58 "F0": signal_parameters["F0"],
59 "F1": signal_parameters["F1"],
60 "F2": signal_parameters["F2"],
61 "asini": {
62 "type": "unif",
63 "lower": 0.9 * signal_parameters["asini"],
64 "upper": 1.1 * signal_parameters["asini"],
65 },
66 "period": {
67 "type": "unif",
68 "lower": 0.9 * signal_parameters["period"],
69 "upper": 1.1 * signal_parameters["period"],
70 },
71 "tp": {
72 "type": "unif",
73 "lower": mid_time - signal_parameters["period"] / 2.0,
74 "upper": mid_time + signal_parameters["period"] / 2.0,
75 },
76 }
77
78 if directed_search:
79 for key in "Alpha", "Delta":
80 theta_prior[key] = signal_parameters[key]
81 else:
82 theta_prior.update(
83 {
84 "Alpha": {
85 "type": "unif",
86 "lower": signal_parameters["Alpha"] - 0.01,
87 "upper": signal_parameters["Alpha"] + 0.01,
88 },
89 "Delta": {
90 "type": "unif",
91 "lower": signal_parameters["Delta"] - 0.01,
92 "upper": signal_parameters["Delta"] + 0.01,
93 },
94 }
95 )
96
97
98 if known_eccentricity:
99 for key in "ecc", "argp":
100 theta_prior[key] = signal_parameters[key]
101 else:
102 theta_prior.update(
103 {
104 "ecc": {
105 "type": "unif",
106 "lower": signal_parameters["ecc"] - 5e-2,
107 "upper": signal_parameters["ecc"] + 5e-2,
108 },
109 "argp": {
110 "type": "unif",
111 "lower": signal_parameters["argp"] - np.pi / 2,
112 "upper": signal_parameters["argp"] + np.pi / 2,
113 },
114 }
115 )
116
117 ntemps = 3
118 log10beta_min = -1
119 nwalkers = 150
120 nsteps = [100, 200]
121
122 mcmc = pyfstat.MCMCSemiCoherentSearch(
123 label=label,
124 outdir=outdir,
125 nsegs=10,
126 sftfilepattern=data.sftfilepath,
127 theta_prior=theta_prior,
128 tref=signal_parameters["tref"],
129 minStartTime=data_parameters["tstart"],
130 maxStartTime=tend,
131 nsteps=nsteps,
132 nwalkers=nwalkers,
133 ntemps=ntemps,
134 log10beta_min=log10beta_min,
135 binary=True,
136 )
137
138 mcmc.run(
139 plot_walkers=True,
140 walker_plot_args={"plot_det_stat": True, "injection_parameters": signal_parameters},
141 )
142 mcmc.plot_corner(add_prior=True, truths=signal_parameters)
143 mcmc.plot_prior_posterior(injection_parameters=signal_parameters)
144 mcmc.print_summary()