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