Using the hierarchical model

The HierarchicalModel is for performing a joint Bayesian analysis over two or more events that share one or more common parameters. For example, you may wish to analyze two different binary neutron star mergers assuming the neutron stars share the same equation of state, or you want to analyze two signals that arrived at different times under the assumption that they are lensed copies of the same event.

The events are assumed to not share the same data, and may not even be observations from the same type of detector. What type of data is read for each event is determined by the model used for that event. Upon initialization, the hierarchical model passes the relevant parameters and sections in the config file to each sub-model’s from_config method. During the analysis, the sub-models’ loglikelihood functions are called and summed over. In that regard, the hierarchical model treats the sub-models as black boxes.

To set up the hierarchical model, you provide a list of sub-model labels in the [model] section of your config file that represent each event to analyze. For each sub-model label provided, you provide a [{label}__model] section that in turn specifies how to initialize that event’s model.

In the [variable_params] and [static_params] sections you specify which parameters belong to which event by prepending the parameter name with the event’s label, followed by a __; i.e., {label}__{param}. To specify that a parameter is common to a sub-set of events, you prepend each event’s label (separated by a _). Parameters that have no event labels prepended to them are treated as common to all events.

For example, say we have three events a, b, and c, who’s models take two parameters, foo and bar. Say we want foo to be common between a and b, but unique to c, and bar to common all events. Our analysis would therefore have three parameters, a_b__foo, c__foo, and bar.

Additional sections that are required by each model should have the model label prepended to them in the same manner. In the above example, if the models for events a, b, and c require a data section, then the config file should have sections [a__data], [b__data], and [c__data].

When sub-models are initialized, the event labels are stripped from the parameter names (and section headers), then passed to the sub-model. Sub-models do not need any special features to be run in a hierarchical analysis; any model that inherits from BaseModel can be used as a sub-model.

Lensing example

To illustrate how to use the HierarchicalModel here we perform a simple lensing analysis of two binary black hole simulations. The two simulations are injected into fake noise ~one week apart, with slightly different apparent sky locations (the difference may not be physically possible; this is only for illustration purposes). We will analyze the two events allowing them to have different sky locations and coalescence times, but sharing the same masses.

First, we need to create the simulations using pycbc_create_injections. To create the two injections we’ll use the configuration files:

# A GW150914-like injection, with approximately the same time and location.

[variable_params]

[static_params]
approximant = IMRPhenomD
tc = 1126259462.430
srcmass1 = 35
srcmass2 = 35
distance = 500
ra = 2.2
dec = -1.2
inclination = 0
polarization = 0
f_ref = 18
f_lower = 18
taper = start

[waveform_transforms-mass1]
name = custom
inputs = srcmass1, distance
mass1 = srcmass1 * (1+redshift(distance))

[waveform_transforms-mass2]
name = custom
inputs = srcmass2, distance
mass2 = srcmass2 * (1+redshift(distance))

Download

# A "lensed" version of event1. This has the same parameters except for the sky
# location and time: the event happens ~one week later, with a slightly
# different sky location. The sky location and time were picked arbitrarily;
# no attempt was made at making sure these are actually physically possible.

[variable_params]

[static_params]
approximant = IMRPhenomD
tc = 1126859462.0
srcmass1 = 35
srcmass2 = 35
distance = 500
ra = 2.0
dec = -1.
inclination = 0
polarization = 0
f_ref = 18
f_lower = 18
taper = start

[waveform_transforms-mass1]
name = custom
inputs = srcmass1, distance
mass1 = srcmass1 * (1+redshift(distance))

[waveform_transforms-mass2]
name = custom
inputs = srcmass2, distance
mass2 = srcmass2 * (1+redshift(distance))

Download

Create the injection hdf files by running:

pycbc_create_injections --verbose --force \
    --ninjections 1 \
    --config-file event1_inj.ini \
    --output-file event1_inj.hdf

pycbc_create_injections --verbose --force \
    --ninjections 1 \
    --config-file event2_inj.ini \
    --output-file event2_inj.hdf

Download

Now we’ll setup the configuration files to run the hierarchical analysis on these two events. First we’ll set the [model] section to use the HierarchicalModel, and tell it to analyze two events called event1 and event2:

[model]
name = hierarchical
submodels = event1 event2

Download

We now need to specify the model to use for each event. For this analysis we’ll use the Relative model for both events:

[event1__model]
name = relative
low-frequency-cutoff = 18.0
high-frequency-cutoff = 1024.0
epsilon = 0.005
mass1_ref = 35
mass2_ref = 35
tc_ref = FROM_INJECTION:tc
ra_ref = FROM_INJECTION:ra
dec_ref = FROM_INJECTION:dec

Download

[event2__model]
name = relative
low-frequency-cutoff = ${event1__model|low-frequency-cutoff}
high-frequency-cutoff = ${event1__model|high-frequency-cutoff}
epsilon = ${event1__model|epsilon}
mass1_ref = ${event1__model|mass1_ref}
mass2_ref = ${event1__model|mass2_ref}
tc_ref = FROM_INJECTION:tc
ra_ref = FROM_INJECTION:ra
dec_ref = FROM_INJECTION:dec

Download

We also need to provide data sections for each event:

# We'll inject the simulated signals into fake noise generated with the aLIGO
# design sensitivity. We need two separate data sections, one for each event.
[event1__data]
instruments = H1 L1
trigger-time = 1126259462.430
analysis-start-time = -6
analysis-end-time = 2
sample-rate = 2048
fake-strain = H1:aLIGOaLIGODesignSensitivityT1800044 L1:aLIGOaLIGODesignSensitivityT1800044
fake-strain-seed = H1:237 L1:82
psd-estimation = median-mean
psd-inverse-length = 8
psd-segment-length = 8
psd-segment-stride = 4
psd-start-time = -256
psd-end-time = 256
channel-name = H1:STRAIN L1:STRAIN
injection-file = event1_inj.hdf
strain-high-pass = 15
pad-data = 8

# For event2, we'll reuse many of the same settings, just changing the seed,
# injection file, and the trigger time.
[event2__data]
instruments = ${event1__data|instruments}
trigger-time = 1126859462.0
analysis-start-time = ${event1__data|analysis-start-time}
analysis-end-time = ${event1__data|analysis-end-time}
sample-rate = ${event1__data|sample-rate}
fake-strain = ${event1__data|fake-strain}
fake-strain-seed = H1:918 L1:6610
psd-estimation = ${event1__data|psd-estimation}
psd-inverse-length = ${event1__data|psd-inverse-length}
psd-segment-length = ${event1__data|psd-segment-length}
psd-segment-stride = ${event1__data|psd-segment-stride}
psd-start-time = ${event1__data|psd-start-time}
psd-end-time = ${event1__data|psd-end-time}
channel-name = ${event1__data|channel-name}
injection-file = event2_inj.hdf
strain-high-pass = ${event1__data|strain-high-pass}
pad-data = ${event1__data|pad-data}

Download

Now the prior:

[variable_params]
# common parameters
srcmchirp =
q =
# parameters unique to event1
event1__delta_tc =
event1__ra =
event1__dec =
# parameters unique to event2
event2__delta_tc =
# We'll make event2's sky location be dependent on
# event1. We'll do that by making the variable
# parameter be a deviation from event1's sky location.
# The ra and dec for event2 will be set with a
# transform, below.
event2__dra =
event2__ddec =

[static_params]
approximant = IMRPhenomD
f_lower = 20
f_ref = 20
distance = 500
inclination = 0
polarization = 0

# Prior bounds taken from 4-OGC
[prior-srcmchirp]
name = mchirp_from_uniform_mass1_mass2
min-srcmchirp = 23
max-srcmchirp = 42

[prior-q]
name = q_from_uniform_mass1_mass2
min-q = 1.
max-q = 4.

[waveform_transforms-mass1+mass2]
name = custom
inputs = srcmchirp, q
mass1 = mass1_from_mchirp_q(srcmchirp,q) * (1 + 0.105)
mass2 = mass2_from_mchirp_q(srcmchirp,q) * (1 + 0.105)

[prior-event1__delta_tc]
name = uniform
min-event1__delta_tc = -0.1
max-event1__delta_tc = 0.1

[prior-event2__delta_tc]
name = uniform
min-event2__delta_tc = -0.1
max-event2__delta_tc = 0.1

# Note that the output of waveform transforms that output parameters specific
# to a particular sub-model must include that model's label in the parameter
# name, just as the variable and static params do.
[waveform_transforms-event1__tc]
name = custom
inputs = event1__delta_tc
event1__tc = ${event1__data|trigger-time} + event1__delta_tc

[waveform_transforms-event2__tc]
name = custom
inputs = event2__delta_tc
event2__tc = ${event2__data|trigger-time} + event2__delta_tc

# We'll use a uniform prior for the sky location of event1...
[prior-event1__ra+event1__dec]
name = uniform_sky
azimuthal-angle = event1__ra
polar-angle = event1__dec

# ...and a Gaussian centered on event1 for event2.
[prior-event2__dra+event2__ddec]
name = gaussian
event2__dra_mean = 0
event2__dra_var = 1
event2__ddec_mean = 0
event2__ddec_var = 1

[waveform_transforms-event2__ra+event2__dec]
name = custom
inputs = event1__ra, event1__dec, event2__dra, event2__ddec
event2__ra = event1__ra + event2__dra
event2__dec = event1__dec + event2__ddec

[constraint-event2dec]
name = custom
constraint_arg = (event1__dec + event2__ddec >= -pi/2) & (event1__dec + event2__ddec <= pi/2)

Download

And finally the sampler configuration:

[sampler]
name = dynesty
dlogz = 0.1
nlive = 2000
checkpoint_time_interval = 1800
maxcall = 10000

Download

Now run:

pycbc_inference \
    --config-files model.ini model-event1_relbin.ini model-event2_relbin.ini prior.ini data.ini dynesty.ini \
    --nprocesses 8 \
    --output-file hierarchical.hdf \
    --seed 10 \
    --force \
    --verbose

Download

When it is done, you will have a file called hierarchical.hdf which contains the results.