pycbc icon indicating copy to clipboard operation
pycbc copied to clipboard

Question about inference with dynesty.

Open yanyuechuixue opened this issue 4 years ago • 0 comments

This is just a question about how to use pycbc inference.

I am use pycbc_inference to estimate the parameters of GW150914. I copy the run files from https://github.com/gwastro/pycbc-inference-paper/tree/master/run_files/GW150914 (with small changes, because the required config file format has changed). To accelerate the sampling, I use dynesty instead of pt_emcee. The full config file is pasted below.

How ever, when I use different random seed, it will give different posterior probalities. For example, with one seed, I got $d_L = 510 +120/ -170 Mpc$, with another, I got $d_L = 470 +140/-180$. However, the probality of a physical parameter should be uniqe. Is it a normal phenomenon for dynesty?

The pycbc version is 1.16.1 My config is :

data.ini

[data]
instruments = H1 L1
trigger-time = 1126259462.43
; See the documentation at
; http://pycbc.org/pycbc/latest/html/inference.html#simulated-bbh-example
; for details on the following settings:
analysis-start-time = -6
analysis-end-time = 2
psd-estimation = median-mean
psd-start-time = -256
psd-end-time = 256
psd-inverse-length = 8
psd-segment-length = 8
psd-segment-stride = 4
; The frame files must be downloaded from GWOSC before running. Here, we
; assume that the files have been downloaded to the same directory. Adjust
; the file path as necessary if not.
frame-files = H1:/media/DATA/gwf_files/H-H1_GWOSC_16KHZ_R1-1126257415-4096.gwf L1:/media/DATA/gwf_files/L-L1_GWOSC_16KHZ_R1-1126257415-4096.gwf
channel-name = H1:GWOSC-16KHZ_R1_STRAIN L1:GWOSC-16KHZ_R1_STRAIN
; this will cause the data to be resampled to 2048 Hz:
sample-rate = 4096
; We'll use a high-pass filter so as not to get numerical errors from the large
; amplitude low frequency noise. Here we use 15 Hz, which is safely below the
; low frequency cutoff of our likelihood integral (20 Hz)
strain-high-pass = 15
; The pad-data argument is for the high-pass filter: 8s are added to the
; beginning/end of the analysis/psd times when the data is loaded. After the
; high pass filter is applied, the additional time is discarded. This pad is
; *in addition to* the time added to the analysis start/end time for the PSD
; inverse length. Since it is discarded before the data is transformed for the
; likelihood integral, it has little affect on the run time.
pad-data = 8

sampler.ini

[sampler]
name = dynesty
nlive = 3000
dlogz = 0.001
checkpoint-interval = 2000

model.ini

[model]
name = marginalized_phase
low-frequency-cutoff = 20.0

[variable_params]
; waveform parameters that will vary in MCMC
delta_tc =
mass1 =
mass2 =
spin1_a =
spin1_azimuthal =
spin1_polar =
spin2_a =
spin2_azimuthal =
spin2_polar =
distance =
inclination =
polarization =
ra =
dec =

[static_params]
; waveform parameters that will not change in MCMC
approximant = IMRPhenomPv2
f_lower = 20
f_ref = 20
; we'll set the tc by using the trigger time in the data
; section of the config file + delta_tc
trigger_time = ${data|trigger-time}

[prior-delta_tc]
; coalescence time prior
name = uniform
min-delta_tc = -0.1
max-delta_tc = 0.1

[waveform_transforms-tc]
; we need to provide tc to the waveform generator
name = custom
inputs = delta_tc
tc = ${data|trigger-time} + delta_tc

[prior-mass1]
name = uniform
min-mass1 = 10.
max-mass1 = 80.

[prior-mass2]
name = uniform
min-mass2 = 10.
max-mass2 = 80.

[prior-spin1_a]
name = uniform
min-spin1_a = 0.0
max-spin1_a = 0.99

[prior-spin1_polar+spin1_azimuthal]
name = uniform_solidangle
polar-angle = spin1_polar
azimuthal-angle = spin1_azimuthal

[prior-spin2_a]
name = uniform
min-spin2_a = 0.0
max-spin2_a = 0.99

[prior-spin2_polar+spin2_azimuthal]
name = uniform_solidangle
polar-angle = spin2_polar
azimuthal-angle = spin2_azimuthal

[prior-distance]
; following gives a uniform volume prior
name = uniform_radius
min-distance = 10
max-distance = 1000


[prior-inclination]
; inclination prior
name = sin_angle

[prior-ra+dec]
; sky position prior
name = uniform_sky

[prior-polarization]
; polarization prior
name = uniform_angle

run.sh

#!/bin/sh

# configuration files
PRIOR_CONFIG=model.ini
DATA_CONFIG=data.ini
SAMPLER_CONFIG=sampler.ini

OUTPUT_PATH=inference.hdf

# the following sets the number of cores to use; adjust as needed to
# your computer's capabilities
NPROCS=56

# run sampler
# Running with OMP_NUM_THREADS=1 stops lalsimulation
# from spawning multiple jobs that would otherwise be used
# by pycbc_inference and cause a reduced runtime.
OMP_NUM_THREADS=1 \
pycbc_inference --verbose \
    --seed 1897234 \
    --config-file ${PRIOR_CONFIG} ${DATA_CONFIG} ${SAMPLER_CONFIG} \
    --output-file ${OUTPUT_PATH} \
    --nprocesses ${NPROCS} \
    --force

Thank you very much!

yanyuechuixue avatar Jun 03 '20 15:06 yanyuechuixue