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