EX3: GW150914#

In this example we load an LVK posterior from zenodo and compute the FI evidence for the CBC model.

Download posteriors from zenodo

%load_ext autoreload
%autoreload 2
%matplotlib inline

import h5py
import numpy as np
import pandas as pd
from collections import namedtuple
from funnel.plotting import plot_fi_evidence_results
from funnel.fi_core import get_fi_lnz_list
import shutil, os

np.random.seed(42)

Load GW150914 posterior + Nested Sampling LnZ#

# load the LVK posterior
FPATH = 'IGWN-GWTC2p1-v2-GW150914_095045_PEDataRelease_mixed_cosmo.h5'

LVK_data = namedtuple("LVK_data", "posterior, lnz, lnz_err, lnBF")


def load_lvk_data(fpath):
    with h5py.File(fpath, 'r') as f:
        sampling_params = list(f['C01:IMRPhenomXPHM/priors/analytic'].keys())
        # remove mass_1, mass_2 (as we already have chirp_mass, mass_ratio)
        sampling_params = [p for p in sampling_params if p not in ['mass_1', 'mass_2']]
        sampler_data = f['C01:IMRPhenomXPHM/meta_data/sampler']

        lnz, lnz_err =  sampler_data['ln_evidence'][0], sampler_data['ln_evidence_error'][0]
        lnBF = sampler_data['ln_bayes_factor'][0]
        post = f['C01:IMRPhenomXPHM']['posterior_samples'][()]
        post = pd.DataFrame({name: post[name][:] for name in post.dtype.names})
        post = post[sampling_params + ['log_likelihood', 'log_prior']]
        # only keep parameters with more than one unique value (ie drop deltas)
        post = post.loc[:, post.nunique() > 1]

    return LVK_data(post, lnz, lnz_err, lnBF)


GW150914_data = load_lvk_data(FPATH)

def get_configs(fpath):
    configs_dict = {}
    with h5py.File(fpath, 'r') as f:
        configs = f['C01:IMRPhenomXPHM/meta_data/meta_data']
        for k, v in configs.items():
            v = str(v[()][0])
            if 'marginalization' in k:
                configs_dict[k] = v
    configs_dict = {k:configs_dict[k] for k in sorted(configs_dict.keys())}
    df=  pd.DataFrame(configs_dict, index=['vals']).T
    return df


configs = get_configs(FPATH)
configs
vals
distance_marginalization b'True'
phase_marginalization b'False'
time_marginalization b'True' 
# GW
# 150914_data.posterior.columns WITHOUT log_likelihood and log_prior
sampling_params = set(list(GW150914_data.posterior.columns))-set(['log_likelihood', 'log_prior'])

cbc_params = []
print("CBC PARAMS:")
count = 0
for i in sorted(sampling_params):
    if 'calib' not in i:
        count += 1
        print(count, ")", i)
        cbc_params.append(i)
print("CALIB PARAMS:")
calib_params = []
count = 0
for i in sorted(sampling_params):
    if 'calib' in i:
        count += 1
        print(count, ")", i)
        calib_params.append(i)

CBC PARAMS:
1 ) a_1
2 ) a_2
3 ) azimuth
4 ) chirp_mass
5 ) geocent_time
6 ) luminosity_distance
7 ) mass_ratio
8 ) phase
9 ) phi_12
10 ) phi_jl
11 ) psi
12 ) theta_jn
13 ) tilt_1
14 ) tilt_2
15 ) time_jitter
16 ) zenith
CALIB PARAMS:
1 ) recalib_H1_amplitude_0
2 ) recalib_H1_amplitude_1
3 ) recalib_H1_amplitude_2
4 ) recalib_H1_amplitude_3
5 ) recalib_H1_amplitude_4
6 ) recalib_H1_amplitude_5
7 ) recalib_H1_amplitude_6
8 ) recalib_H1_amplitude_7
9 ) recalib_H1_amplitude_8
10 ) recalib_H1_amplitude_9
11 ) recalib_H1_phase_0
12 ) recalib_H1_phase_1
13 ) recalib_H1_phase_2
14 ) recalib_H1_phase_3
15 ) recalib_H1_phase_4
16 ) recalib_H1_phase_5
17 ) recalib_H1_phase_6
18 ) recalib_H1_phase_7
19 ) recalib_H1_phase_8
20 ) recalib_H1_phase_9
21 ) recalib_L1_amplitude_0
22 ) recalib_L1_amplitude_1
23 ) recalib_L1_amplitude_2
24 ) recalib_L1_amplitude_3
25 ) recalib_L1_amplitude_4
26 ) recalib_L1_amplitude_5
27 ) recalib_L1_amplitude_6
28 ) recalib_L1_amplitude_7
29 ) recalib_L1_amplitude_8
30 ) recalib_L1_amplitude_9
31 ) recalib_L1_phase_0
32 ) recalib_L1_phase_1
33 ) recalib_L1_phase_2
34 ) recalib_L1_phase_3
35 ) recalib_L1_phase_4
36 ) recalib_L1_phase_5
37 ) recalib_L1_phase_6
38 ) recalib_L1_phase_7
39 ) recalib_L1_phase_8
40 ) recalib_L1_phase_9

data = h5py.File(FPATH, 'r')
count = 0
skipped = []
for k,v in data['C01:IMRPhenomXPHM/priors/analytic'].items():
    if "DeltaFunction" not in str(v[()]):
        count += 1
        print(f"{count:2d}) {k}: {v[()]}")
    else:
        skipped.append(f"SKIP) {k}:  {v[()]}")
data.close()

for i in skipped:
    print(i)

 1) a_1: [b"Uniform(minimum=0, maximum=0.99, name='a_1', latex_label='$a_1$', unit=None, boundary=None)"]
 2) a_2: [b"Uniform(minimum=0, maximum=0.99, name='a_2', latex_label='$a_2$', unit=None, boundary=None)"]
 3) azimuth: [b"Uniform(minimum=0, maximum=6.283185307179586, name=None, latex_label='$\\\\epsilon$', unit=None, boundary='periodic')"]
 4) chirp_mass: [b"UniformInComponentsChirpMass(minimum=21.418182160215295, maximum=41.97447913941358, name='chirp_mass', latex_label='$\\\\mathcal{M}$', unit='$M_{\\\\odot}$', boundary=None)"]
 5) geocent_time: [b"Uniform(minimum=1126259462.2910001, maximum=1126259462.491, name='geocent_time', latex_label='$t_c$', unit='$s$', boundary=None)"]
 6) luminosity_distance: [b"PowerLaw(alpha=2, minimum=10, maximum=10000, name='luminosity_distance', latex_label='$d_L$', unit='Mpc', boundary=None)"]
 7) mass_1: [b"Constraint(minimum=1, maximum=1000, name='mass_1', latex_label='$m_1$', unit=None)"]
 8) mass_2: [b"Constraint(minimum=1, maximum=1000, name='mass_2', latex_label='$m_2$', unit=None)"]
 9) mass_ratio: [b"UniformInComponentsMassRatio(minimum=0.05, maximum=1.0, name='mass_ratio', latex_label='$q$', unit=None, boundary=None)"]
10) phase: [b"Uniform(minimum=0, maximum=6.283185307179586, name='phase', latex_label='$\\\\phi$', unit=None, boundary='periodic')"]
11) phi_12: [b"Uniform(minimum=0, maximum=6.283185307179586, name='phi_12', latex_label='$\\\\Delta\\\\phi$', unit=None, boundary='periodic')"]
12) phi_jl: [b"Uniform(minimum=0, maximum=6.283185307179586, name='phi_jl', latex_label='$\\\\phi_{JL}$', unit=None, boundary='periodic')"]
13) psi: [b"Uniform(minimum=0, maximum=3.141592653589793, name='psi', latex_label='$\\\\psi$', unit=None, boundary='periodic')"]
14) recalib_H1_amplitude_0: [b"Gaussian(mu=0.0022497139999999582, sigma=0.032482292999999995, name='recalib_H1_amplitude_0', latex_label='$A^H1_0$', unit=None, boundary='reflective')"]
15) recalib_H1_amplitude_1: [b"Gaussian(mu=-0.0022525476996070083, sigma=0.03909717651863342, name='recalib_H1_amplitude_1', latex_label='$A^H1_1$', unit=None, boundary='reflective')"]
16) recalib_H1_amplitude_2: [b"Gaussian(mu=0.004632777733005889, sigma=0.04078174073734512, name='recalib_H1_amplitude_2', latex_label='$A^H1_2$', unit=None, boundary='reflective')"]
17) recalib_H1_amplitude_3: [b"Gaussian(mu=0.003987700439142069, sigma=0.03287682593561539, name='recalib_H1_amplitude_3', latex_label='$A^H1_3$', unit=None, boundary='reflective')"]
18) recalib_H1_amplitude_4: [b"Gaussian(mu=0.0018107999614999065, sigma=0.01936965624646965, name='recalib_H1_amplitude_4', latex_label='$A^H1_4$', unit=None, boundary='reflective')"]
19) recalib_H1_amplitude_5: [b"Gaussian(mu=0.003528123018086278, sigma=0.012591404744135145, name='recalib_H1_amplitude_5', latex_label='$A^H1_5$', unit=None, boundary='reflective')"]
20) recalib_H1_amplitude_6: [b"Gaussian(mu=-0.0018296964612414133, sigma=0.01075481852283716, name='recalib_H1_amplitude_6', latex_label='$A^H1_6$', unit=None, boundary='reflective')"]
21) recalib_H1_amplitude_7: [b"Gaussian(mu=-0.005059244326670863, sigma=0.008407555954093059, name='recalib_H1_amplitude_7', latex_label='$A^H1_7$', unit=None, boundary='reflective')"]
22) recalib_H1_amplitude_8: [b"Gaussian(mu=-0.0010763014001435082, sigma=0.01391150436119593, name='recalib_H1_amplitude_8', latex_label='$A^H1_8$', unit=None, boundary='reflective')"]
23) recalib_H1_amplitude_9: [b"Gaussian(mu=-0.0017064540000000237, sigma=0.016027046000000045, name='recalib_H1_amplitude_9', latex_label='$A^H1_9$', unit=None, boundary='reflective')"]
24) recalib_H1_phase_0: [b"Gaussian(mu=-0.009355752999999998, sigma=0.030123145, name='recalib_H1_phase_0', latex_label='$\\\\phi^H1_0$', unit=None, boundary='reflective')"]
25) recalib_H1_phase_1: [b"Gaussian(mu=0.0050309026057799855, sigma=0.04177354702981138, name='recalib_H1_phase_1', latex_label='$\\\\phi^H1_1$', unit=None, boundary='reflective')"]
26) recalib_H1_phase_2: [b"Gaussian(mu=0.0043774230879091115, sigma=0.045027737962814354, name='recalib_H1_phase_2', latex_label='$\\\\phi^H1_2$', unit=None, boundary='reflective')"]
27) recalib_H1_phase_3: [b"Gaussian(mu=-0.0023545243645386326, sigma=0.036339899119904644, name='recalib_H1_phase_3', latex_label='$\\\\phi^H1_3$', unit=None, boundary='reflective')"]
28) recalib_H1_phase_4: [b"Gaussian(mu=-0.00444209518471385, sigma=0.032837260835226936, name='recalib_H1_phase_4', latex_label='$\\\\phi^H1_4$', unit=None, boundary='reflective')"]
29) recalib_H1_phase_5: [b"Gaussian(mu=-0.006872683631006688, sigma=0.02806419946997496, name='recalib_H1_phase_5', latex_label='$\\\\phi^H1_5$', unit=None, boundary='reflective')"]
30) recalib_H1_phase_6: [b"Gaussian(mu=-0.01508554971905464, sigma=0.022470973360749408, name='recalib_H1_phase_6', latex_label='$\\\\phi^H1_6$', unit=None, boundary='reflective')"]
31) recalib_H1_phase_7: [b"Gaussian(mu=-0.009105656795490394, sigma=0.0186470280190447, name='recalib_H1_phase_7', latex_label='$\\\\phi^H1_7$', unit=None, boundary='reflective')"]
32) recalib_H1_phase_8: [b"Gaussian(mu=-0.01327344249855343, sigma=0.019250972190353967, name='recalib_H1_phase_8', latex_label='$\\\\phi^H1_8$', unit=None, boundary='reflective')"]
33) recalib_H1_phase_9: [b"Gaussian(mu=-0.006712540999999999, sigma=0.020072625, name='recalib_H1_phase_9', latex_label='$\\\\phi^H1_9$', unit=None, boundary='reflective')"]
34) recalib_L1_amplitude_0: [b"Gaussian(mu=0.005476730999999901, sigma=0.08406088199999993, name='recalib_L1_amplitude_0', latex_label='$A^L1_0$', unit=None, boundary='reflective')"]
35) recalib_L1_amplitude_1: [b"Gaussian(mu=-0.012104387120680028, sigma=0.07651583147867605, name='recalib_L1_amplitude_1', latex_label='$A^L1_1$', unit=None, boundary='reflective')"]
36) recalib_L1_amplitude_2: [b"Gaussian(mu=-0.013219120050478253, sigma=0.06508029899775643, name='recalib_L1_amplitude_2', latex_label='$A^L1_2$', unit=None, boundary='reflective')"]
37) recalib_L1_amplitude_3: [b"Gaussian(mu=-0.006657329180423924, sigma=0.05665180922845485, name='recalib_L1_amplitude_3', latex_label='$A^L1_3$', unit=None, boundary='reflective')"]
38) recalib_L1_amplitude_4: [b"Gaussian(mu=-0.0026286630232402308, sigma=0.05041154190241824, name='recalib_L1_amplitude_4', latex_label='$A^L1_4$', unit=None, boundary='reflective')"]
39) recalib_L1_amplitude_5: [b"Gaussian(mu=-0.008989482013408759, sigma=0.0391433695986004, name='recalib_L1_amplitude_5', latex_label='$A^L1_5$', unit=None, boundary='reflective')"]
40) recalib_L1_amplitude_6: [b"Gaussian(mu=-0.0112041341409286, sigma=0.02176964617289734, name='recalib_L1_amplitude_6', latex_label='$A^L1_6$', unit=None, boundary='reflective')"]
41) recalib_L1_amplitude_7: [b"Gaussian(mu=-0.007323119754866439, sigma=0.01662038622968457, name='recalib_L1_amplitude_7', latex_label='$A^L1_7$', unit=None, boundary='reflective')"]
42) recalib_L1_amplitude_8: [b"Gaussian(mu=-0.007118273569718694, sigma=0.017534092484121667, name='recalib_L1_amplitude_8', latex_label='$A^L1_8$', unit=None, boundary='reflective')"]
43) recalib_L1_amplitude_9: [b"Gaussian(mu=-0.00654931000000003, sigma=0.018560549000000037, name='recalib_L1_amplitude_9', latex_label='$A^L1_9$', unit=None, boundary='reflective')"]
44) recalib_L1_phase_0: [b"Gaussian(mu=-0.027121368999999992, sigma=0.065822439, name='recalib_L1_phase_0', latex_label='$\\\\phi^L1_0$', unit=None, boundary='reflective')"]
45) recalib_L1_phase_1: [b"Gaussian(mu=-0.015666055283929895, sigma=0.06971337917415898, name='recalib_L1_phase_1', latex_label='$\\\\phi^L1_1$', unit=None, boundary='reflective')"]
46) recalib_L1_phase_2: [b"Gaussian(mu=-0.0021521761445731018, sigma=0.06846083085274324, name='recalib_L1_phase_2', latex_label='$\\\\phi^L1_2$', unit=None, boundary='reflective')"]
47) recalib_L1_phase_3: [b"Gaussian(mu=0.0022536976421022887, sigma=0.06194252192559962, name='recalib_L1_phase_3', latex_label='$\\\\phi^L1_3$', unit=None, boundary='reflective')"]
48) recalib_L1_phase_4: [b"Gaussian(mu=-0.002461783883218465, sigma=0.048859883075223744, name='recalib_L1_phase_4', latex_label='$\\\\phi^L1_4$', unit=None, boundary='reflective')"]
49) recalib_L1_phase_5: [b"Gaussian(mu=-0.0067988571593666286, sigma=0.03330895148538794, name='recalib_L1_phase_5', latex_label='$\\\\phi^L1_5$', unit=None, boundary='reflective')"]
50) recalib_L1_phase_6: [b"Gaussian(mu=0.0015671495057928866, sigma=0.0194164413063868, name='recalib_L1_phase_6', latex_label='$\\\\phi^L1_6$', unit=None, boundary='reflective')"]
51) recalib_L1_phase_7: [b"Gaussian(mu=0.0016729416638600334, sigma=0.013796249113051186, name='recalib_L1_phase_7', latex_label='$\\\\phi^L1_7$', unit=None, boundary='reflective')"]
52) recalib_L1_phase_8: [b"Gaussian(mu=0.003950310144720625, sigma=0.014338015652101466, name='recalib_L1_phase_8', latex_label='$\\\\phi^L1_8$', unit=None, boundary='reflective')"]
53) recalib_L1_phase_9: [b"Gaussian(mu=0.008061858000000002, sigma=0.015019266, name='recalib_L1_phase_9', latex_label='$\\\\phi^L1_9$', unit=None, boundary='reflective')"]
54) theta_jn: [b"Sine(minimum=0, maximum=3.141592653589793, name='theta_jn', latex_label='$\\\\theta_{JN}$', unit=None, boundary=None)"]
55) tilt_1: [b"Sine(minimum=0, maximum=3.141592653589793, name='tilt_1', latex_label='$\\\\theta_1$', unit=None, boundary=None)"]
56) tilt_2: [b"Sine(minimum=0, maximum=3.141592653589793, name='tilt_2', latex_label='$\\\\theta_2$', unit=None, boundary=None)"]
57) time_jitter: [b"Uniform(minimum=-0.00048828125, maximum=0.00048828125, name=None, latex_label=None, unit=None, boundary='periodic')"]
58) zenith: [b"Sine(minimum=0, maximum=3.141592653589793, name=None, latex_label='$\\\\kappa$', unit=None, boundary=None)"]
SKIP) recalib_H1_frequency_0:  [b"DeltaFunction(peak=19.999999999999996, name='recalib_H1_frequency_0', latex_label='$f^H1_0$', unit=None)"]
SKIP) recalib_H1_frequency_1:  [b"DeltaFunction(peak=30.514434824844255, name='recalib_H1_frequency_1', latex_label='$f^H1_1$', unit=None)"]
SKIP) recalib_H1_frequency_2:  [b"DeltaFunction(peak=46.55653663398341, name='recalib_H1_frequency_2', latex_label='$f^H1_2$', unit=None)"]
SKIP) recalib_H1_frequency_3:  [b"DeltaFunction(peak=71.03232013940804, name='recalib_H1_frequency_3', latex_label='$f^H1_3$', unit=None)"]
SKIP) recalib_H1_frequency_4:  [b"DeltaFunction(peak=108.37555516757195, name='recalib_H1_frequency_4', latex_label='$f^H1_4$', unit=None)"]
SKIP) recalib_H1_frequency_5:  [b"DeltaFunction(peak=165.35094073835938, name='recalib_H1_frequency_5', latex_label='$f^H1_5$', unit=None)"]
SKIP) recalib_H1_frequency_6:  [b"DeltaFunction(peak=252.27952521936786, name='recalib_H1_frequency_6', latex_label='$f^H1_6$', unit=None)"]
SKIP) recalib_H1_frequency_7:  [b"DeltaFunction(peak=384.9083564974527, name='recalib_H1_frequency_7', latex_label='$f^H1_7$', unit=None)"]
SKIP) recalib_H1_frequency_8:  [b"DeltaFunction(peak=587.2630478939719, name='recalib_H1_frequency_8', latex_label='$f^H1_8$', unit=None)"]
SKIP) recalib_H1_frequency_9:  [b"DeltaFunction(peak=895.9999999999999, name='recalib_H1_frequency_9', latex_label='$f^H1_9$', unit=None)"]
SKIP) recalib_L1_frequency_0:  [b"DeltaFunction(peak=19.999999999999996, name='recalib_L1_frequency_0', latex_label='$f^L1_0$', unit=None)"]
SKIP) recalib_L1_frequency_1:  [b"DeltaFunction(peak=30.514434824844255, name='recalib_L1_frequency_1', latex_label='$f^L1_1$', unit=None)"]
SKIP) recalib_L1_frequency_2:  [b"DeltaFunction(peak=46.55653663398341, name='recalib_L1_frequency_2', latex_label='$f^L1_2$', unit=None)"]
SKIP) recalib_L1_frequency_3:  [b"DeltaFunction(peak=71.03232013940804, name='recalib_L1_frequency_3', latex_label='$f^L1_3$', unit=None)"]
SKIP) recalib_L1_frequency_4:  [b"DeltaFunction(peak=108.37555516757195, name='recalib_L1_frequency_4', latex_label='$f^L1_4$', unit=None)"]
SKIP) recalib_L1_frequency_5:  [b"DeltaFunction(peak=165.35094073835938, name='recalib_L1_frequency_5', latex_label='$f^L1_5$', unit=None)"]
SKIP) recalib_L1_frequency_6:  [b"DeltaFunction(peak=252.27952521936786, name='recalib_L1_frequency_6', latex_label='$f^L1_6$', unit=None)"]
SKIP) recalib_L1_frequency_7:  [b"DeltaFunction(peak=384.9083564974527, name='recalib_L1_frequency_7', latex_label='$f^L1_7$', unit=None)"]
SKIP) recalib_L1_frequency_8:  [b"DeltaFunction(peak=587.2630478939719, name='recalib_L1_frequency_8', latex_label='$f^L1_8$', unit=None)"]
SKIP) recalib_L1_frequency_9:  [b"DeltaFunction(peak=895.9999999999999, name='recalib_L1_frequency_9', latex_label='$f^L1_9$', unit=None)"]

print(f"Posterior shape: {GW150914_data.posterior.shape}")
print(f"LnZ: {GW150914_data.lnz:.2f} +/- {GW150914_data.lnz_err:.2f}")
print(f"LnBF: {GW150914_data.lnBF}")
GW150914_data.posterior.head().T

Posterior shape: (147634, 58)
LnZ: -6984.67 +/- 0.14
LnBF: 303.45
0 1 2 3 4
a_1 9.242101e-01 6.473688e-01 2.056780e-01 7.112512e-01 2.506404e-01
a_2 3.310924e-01 3.133050e-01 8.750084e-01 4.800071e-03 2.271518e-01
azimuth 3.278707e+00 2.602855e+00 2.715020e+00 2.380721e+00 2.783846e+00
chirp_mass 2.918072e+01 2.995305e+01 3.143389e+01 3.074103e+01 3.127060e+01
geocent_time 1.126259e+09 1.126259e+09 1.126259e+09 1.126259e+09 1.126259e+09
luminosity_distance 3.237332e+02 5.109554e+02 5.009286e+02 5.780114e+02 5.871190e+02
mass_ratio 7.879585e-01 8.641748e-01 8.520294e-01 9.803410e-01 9.301303e-01
phase 3.840744e+00 4.337254e+00 1.757881e+00 4.723087e+00 5.542601e+00
phi_12 5.394340e+00 3.868170e+00 2.916934e+00 4.591901e+00 6.435515e-01
phi_jl 1.210051e+00 2.486261e-02 5.576937e+00 9.338661e-01 6.019726e+00
psi 1.800080e+00 7.571376e-01 1.272124e+00 8.861790e-01 2.418739e+00
recalib_H1_amplitude_0 2.999532e-02 -5.472740e-03 6.900000e-03 6.081133e-03 4.977923e-03
recalib_H1_amplitude_1 -3.383615e-02 -5.474793e-02 1.200648e-02 -4.799967e-02 6.687969e-03
recalib_H1_amplitude_2 -1.688680e-02 2.209578e-02 6.758393e-02 1.357331e-02 2.604320e-02
recalib_H1_amplitude_3 1.698306e-02 3.626471e-02 -1.131780e-02 -1.152776e-02 -4.927244e-02
recalib_H1_amplitude_4 -3.305150e-02 1.121430e-02 5.399314e-03 2.998931e-03 -1.089422e-02
recalib_H1_amplitude_5 7.041142e-03 -1.582923e-02 -1.287890e-02 -1.776959e-03 -1.428993e-02
recalib_H1_amplitude_6 6.394128e-04 8.274103e-03 -1.444985e-02 -1.033461e-02 -5.928833e-03
recalib_H1_amplitude_7 -1.628103e-02 1.094559e-02 -1.589106e-02 -7.269532e-04 -1.055759e-02
recalib_H1_amplitude_8 -1.309187e-02 -2.280925e-03 4.110290e-03 -9.604852e-03 3.810975e-03
recalib_H1_amplitude_9 3.221316e-02 -1.607438e-02 2.675564e-02 1.992964e-03 7.210042e-03
recalib_H1_phase_0 -2.719320e-02 -4.060107e-02 -7.166049e-03 6.323885e-03 -1.256244e-01
recalib_H1_phase_1 -1.345580e-02 5.485112e-02 -6.924302e-02 8.663208e-02 -2.985714e-02
recalib_H1_phase_2 7.148825e-02 9.037495e-03 -2.528264e-02 3.690566e-02 -1.084741e-03
recalib_H1_phase_3 5.833606e-02 -1.235240e-04 4.403403e-02 1.949062e-03 -4.622115e-02
recalib_H1_phase_4 -1.974057e-02 4.740506e-02 2.063412e-02 -3.390864e-02 2.773227e-02
recalib_H1_phase_5 -4.427238e-03 -4.936114e-03 -4.498003e-04 -2.420854e-02 -2.673151e-02
recalib_H1_phase_6 2.205619e-02 -1.312319e-02 -3.284825e-02 3.316393e-03 -8.395785e-04
recalib_H1_phase_7 1.259316e-03 -1.696564e-02 -1.007685e-02 -3.100947e-02 2.684205e-03
recalib_H1_phase_8 1.409629e-02 9.934225e-03 -6.075666e-03 2.741046e-02 1.409607e-03
recalib_H1_phase_9 7.438222e-03 1.741232e-02 -1.067432e-02 -1.260970e-02 4.605238e-03
recalib_L1_amplitude_0 -1.149459e-01 -1.533687e-01 5.256302e-02 -4.926279e-03 2.115813e-01
recalib_L1_amplitude_1 -4.511059e-02 -4.943587e-02 9.095780e-02 3.852026e-03 -1.819107e-01
recalib_L1_amplitude_2 -7.530716e-02 8.390149e-02 -4.031972e-02 -1.972729e-02 6.843934e-02
recalib_L1_amplitude_3 -4.488562e-02 4.857405e-02 1.939737e-02 -3.916733e-05 1.209214e-01
recalib_L1_amplitude_4 -4.366417e-02 -3.900651e-02 -5.958842e-02 -7.124533e-02 5.328298e-02
recalib_L1_amplitude_5 3.566428e-02 -4.372492e-02 5.438995e-02 2.581318e-02 -7.917951e-03
recalib_L1_amplitude_6 -3.862231e-03 -2.591929e-02 -4.370625e-02 2.947175e-03 -3.395698e-02
recalib_L1_amplitude_7 -5.773761e-03 -4.311472e-03 4.900587e-03 -1.826212e-02 -1.684315e-02
recalib_L1_amplitude_8 9.266035e-04 -1.677954e-02 -3.511803e-02 -6.485830e-03 -2.241593e-02
recalib_L1_amplitude_9 1.614393e-02 -2.712452e-02 1.511530e-02 3.721168e-03 -3.419746e-02
recalib_L1_phase_0 2.737868e-02 -4.073929e-03 6.091229e-02 -2.545232e-02 -1.465972e-01
recalib_L1_phase_1 5.689561e-02 -1.122678e-02 5.452316e-02 -1.197389e-02 -1.613542e-01
recalib_L1_phase_2 -3.519969e-03 -8.114218e-02 -3.653476e-02 -5.637367e-02 -6.982944e-02
recalib_L1_phase_3 -6.189141e-02 2.536081e-02 -1.243102e-01 7.364353e-02 -2.765010e-02
recalib_L1_phase_4 4.976252e-02 -4.251552e-02 -1.996706e-02 -3.091622e-02 6.611963e-02
recalib_L1_phase_5 2.161828e-02 -3.238076e-02 9.706364e-03 -5.364514e-02 -5.986921e-03
recalib_L1_phase_6 -3.086198e-04 6.332725e-03 -2.956626e-03 -1.210257e-02 -5.876330e-03
recalib_L1_phase_7 7.216118e-03 -3.571173e-03 2.202440e-03 -9.887829e-03 9.418501e-04
recalib_L1_phase_8 -8.043966e-03 1.301929e-03 5.764206e-03 5.924480e-03 -1.411924e-02
recalib_L1_phase_9 1.215333e-02 -7.934553e-03 -2.742098e-03 1.089411e-02 1.026025e-02
theta_jn 2.775318e+00 2.660558e+00 2.493798e+00 2.967835e+00 3.074509e+00
tilt_1 1.913304e+00 1.841044e+00 2.365895e+00 1.672429e+00 1.373358e+00
tilt_2 1.676054e+00 1.839976e+00 1.369656e+00 6.271606e-01 1.525285e+00
time_jitter 2.283661e-04 3.247132e-04 1.200043e-04 -2.701632e-04 4.351720e-04
zenith 2.326161e+00 2.334247e+00 2.345314e+00 2.369513e+00 2.339495e+00
log_likelihood 3.223043e+02 3.203015e+02 3.206540e+02 3.242657e+02 3.227294e+02
log_prior 7.693796e+01 7.862083e+01 7.603158e+01 8.132604e+01 6.226660e+01

NOTE: the log-likelihood column is actually the lnl-noise Lnl (the log-likelihood-ratio).

Compute FI LnZ#

def run_fi_computer(CBC_data, params, frac_samp=0.1, n_ref_points=1000, outdir='out_inj_downsampled', clean=False, ):
    if os.path.exists(outdir) and clean:
        shutil.rmtree(outdir)
    os.makedirs(outdir, exist_ok=True)

    n_total = len(CBC_data.posterior)
    n_samp = int(n_total * frac_samp)

    post = CBC_data.posterior.sample(
        n_samp,
        weights=np.exp(CBC_data.posterior.log_likelihood)
    )
    print(f"Using {100 * frac_samp:.2f}% of posterior, and {n_ref_points} reference points.")
    lnzs, r_vals, _ = get_fi_lnz_list(
        post[params + ['log_likelihood', 'log_prior']],
        num_ref_params=n_ref_points,
        r_vals=np.geomspace(50, 1e4, 100),
        cache_fn=f'{outdir}/lnzs.npz',
        weight_samples_by_lnl=True,
    )
    return lnzs, r_vals

Using a downsampled posterior sample for speed#

parms = sampling_params.copy() # remove geocent_time
parms.remove('geocent_time')
parms.remove('luminosity_distance')
parms.remove('time_jitter')
parms = list(parms)

OUTDIR = 'out_GW150914_downsampled'
CLEAN = True
lnzs, r_vals = run_fi_computer(GW150914_data, parms, outdir=OUTDIR, clean=CLEAN, n_ref_points=10, frac_samp=0.1)
plt_kwgs = dict(lnzs=lnzs, r_vals=r_vals,sampling_lnz=[GW150914_data.lnBF])
fig = plot_fi_evidence_results(**plt_kwgs, plot_all_lnzs=True)
fig.gca().set_ylabel("lnBF");

|funnel|INFO| Calculating FI LnZ with 10 reference points and a posterior of size:(14763, 53) 
|funnel|INFO| Posterior columns:
a_1
a_2
azimuth
chirp_mass
log_likelihood
log_prior
mass_ratio
phase
phi_12
phi_jl
psi
recalib_H1_amplitude_0
recalib_H1_amplitude_1
recalib_H1_amplitude_2
recalib_H1_amplitude_3
recalib_H1_amplitude_4
recalib_H1_amplitude_5
recalib_H1_amplitude_6
recalib_H1_amplitude_7
recalib_H1_amplitude_8
recalib_H1_amplitude_9
recalib_H1_phase_0
recalib_H1_phase_1
recalib_H1_phase_2
recalib_H1_phase_3
recalib_H1_phase_4
recalib_H1_phase_5
recalib_H1_phase_6
recalib_H1_phase_7
recalib_H1_phase_8
recalib_H1_phase_9
recalib_L1_amplitude_0
recalib_L1_amplitude_1
recalib_L1_amplitude_2
recalib_L1_amplitude_3
recalib_L1_amplitude_4
recalib_L1_amplitude_5
recalib_L1_amplitude_6
recalib_L1_amplitude_7
recalib_L1_amplitude_8
recalib_L1_amplitude_9
recalib_L1_phase_0
recalib_L1_phase_1
recalib_L1_phase_2
recalib_L1_phase_3
recalib_L1_phase_4
recalib_L1_phase_5
recalib_L1_phase_6
recalib_L1_phase_7
recalib_L1_phase_8
recalib_L1_phase_9
theta_jn
tilt_1
tilt_2
zenith 

Using 10.00% of posterior, and 10 reference points.
../_images/229345c673d92428d60dc201ca2f484a3a935d6de1ebd4be3a0fb419f9ec8ad2.png 
plt_kwgs = dict(lnzs=lnzs, r_vals=r_vals,sampling_lnz=[GW150914_data.lnBF])
fig = plot_fi_evidence_results(**plt_kwgs, plot_all_lnzs=True)
fig.gca().set_ylabel("lnBF");
../_images/d46c21c68e5234b770f922179af61cd83464a15bf3cb0a0fa0b4f932e0da0efe.png 
plt_kwgs = dict(lnzs=lnzs, r_vals=r_vals,sampling_lnz=[GW150914_data.lnBF])
fig = plot_fi_evidence_results(**plt_kwgs)
fig.gca().set_ylabel("lnBF");
../_images/7431bc76980f9a0916105e299d62dab0c52eebc18d7bc5ee68a44ac307412c96.png 
plt_kwgs = dict(lnzs=lnzs, r_vals=r_vals,sampling_lnz=[GW150914_data.lnBF])
fig = plot_fi_evidence_results(**plt_kwgs, plot_all_lnzs=True)
fig.gca().set_ylabel("lnBF");
../_images/38ced40abfcd2deb0c2b26ee0e07c7026ed1ff8410a5855fa7a7d0f02da2d6b7.png 
# OLD
plt_kwgs = dict(lnzs=lnzs, r_vals=r_vals,sampling_lnz=[GW150914_data.lnBF])
fig = plot_fi_evidence_results(**plt_kwgs)
fig.gca().set_ylabel("lnBF");
../_images/c99ba539ea659fb68417a9ed8fe90375cf8fd38c64eec6e6357e67215134f68c.png 
OUTDIR = 'out_GW150914_downsampled_no_calib'
CLEAN = True
if os.path.exists(OUTDIR) and CLEAN:
    shutil.rmtree(OUTDIR)
os.makedirs(OUTDIR, exist_ok=True)





# USING DOWNSAMPLED POSTERIOR SAMPLES
N_SAMP = 10000
N_REF_POINTS = 10
weights = GW150914_data.posterior.log_likelihood - GW150914_data.posterior.log_likelihood.max()
weights = np.exp(weights)
post = GW150914_data.posterior.sample(N_SAMP, weights=weights)
print(f"Using {100*(len(post)/len(GW150914_data.posterior)):.2f}% of posterior samples, and trying out {N_REF_POINTS} reference points.")

# post = post.drop(columns=calib_params)

lnzs, r_vals, fi_samps = get_fi_lnz_list(
    post,
    num_ref_params=N_REF_POINTS,
    r_vals=np.geomspace(1e2, 1e5, 10),
    cache_fn=f'{OUTDIR}/lnzs.npz',
    weight_samples_by_lnl=True,
)

plt_kwgs = dict(lnzs=lnzs, r_vals=r_vals,sampling_lnz=[GW150914_data.lnBF])
fig = plot_fi_evidence_results(**plt_kwgs, plot_all_lnzs=False, plt_kwgs=dict(alpha=.1))
fig.gca().set_ylabel("lnBF");

|funnel|INFO| Calculating FI LnZ with 10 reference points and a posterior of size: (10000, 55) 

Using 6.77% of posterior samples, and trying out 10 reference points.
../_images/4a3eadc50e6dc23cc3ddb2a11f3aedbd8897f259f48a49032521fa2aa309d495.png 
from funnel.plotting import plot_corner_and_mark_samples

p =  GW150914_data.posterior
weights = p.log_likelihood - p.log_likelihood.max()
weights = np.exp(weights)
samp = p.sample(5, weights=weights)
p = p[cbc_params]
fi_samps = samp[cbc_params].values
fig = plot_corner_and_mark_samples(p, fi_samps)
../_images/b3e4f3a902b24581fa9e03a4d09bb129a8c34dc44c0b365015a5e6f2424b1ed6.png 
sampling_params

{'a_1',
 'a_2',
 'azimuth',
 'chirp_mass',
 'geocent_time',
 'luminosity_distance',
 'mass_ratio',
 'phase',
 'phi_12',
 'phi_jl',
 'psi',
 'recalib_H1_amplitude_0',
 'recalib_H1_amplitude_1',
 'recalib_H1_amplitude_2',
 'recalib_H1_amplitude_3',
 'recalib_H1_amplitude_4',
 'recalib_H1_amplitude_5',
 'recalib_H1_amplitude_6',
 'recalib_H1_amplitude_7',
 'recalib_H1_amplitude_8',
 'recalib_H1_amplitude_9',
 'recalib_H1_phase_0',
 'recalib_H1_phase_1',
 'recalib_H1_phase_2',
 'recalib_H1_phase_3',
 'recalib_H1_phase_4',
 'recalib_H1_phase_5',
 'recalib_H1_phase_6',
 'recalib_H1_phase_7',
 'recalib_H1_phase_8',
 'recalib_H1_phase_9',
 'recalib_L1_amplitude_0',
 'recalib_L1_amplitude_1',
 'recalib_L1_amplitude_2',
 'recalib_L1_amplitude_3',
 'recalib_L1_amplitude_4',
 'recalib_L1_amplitude_5',
 'recalib_L1_amplitude_6',
 'recalib_L1_amplitude_7',
 'recalib_L1_amplitude_8',
 'recalib_L1_amplitude_9',
 'recalib_L1_phase_0',
 'recalib_L1_phase_1',
 'recalib_L1_phase_2',
 'recalib_L1_phase_3',
 'recalib_L1_phase_4',
 'recalib_L1_phase_5',
 'recalib_L1_phase_6',
 'recalib_L1_phase_7',
 'recalib_L1_phase_8',
 'recalib_L1_phase_9',
 'theta_jn',
 'tilt_1',
 'tilt_2',
 'zenith'}

Using the full posterior sample but a few refernce points#

OUTDIR = 'out_GW150914'
CLEAN = True
if os.path.exists(OUTDIR) and CLEAN:
    shutil.rmtree(OUTDIR)
os.makedirs(OUTDIR, exist_ok=True)

N_REF_POINTS = 10
post = GW150914_data.posterior
lnzs, r_vals, _ = get_fi_lnz_list(
    post,
    num_ref_params=N_REF_POINTS,
    r_vals=np.geomspace(50, 1e5, 100),
    cache_fn=f'{OUTDIR}/lnzs.npz',
)

|funnel|INFO| Calculating FI LnZ with 10 reference points and a posterior of size: (147634, 55) 

---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
Cell In[6], line 9
      7 N_REF_POINTS = 10
      8 post = GW150914_data.posterior
----> 9 lnzs, r_vals, _ = get_fi_lnz_list(
     10     post,
     11     num_ref_params=N_REF_POINTS,
     12     r_vals=np.geomspace(50, 1e5, 100),
     13     cache_fn=f'{OUTDIR}/lnzs.npz',
     14 )

File ~/Documents/projects/funnel/src/funnel/fi_core.py:142, in get_fi_lnz_list(posterior_samples, r_vals, num_ref_params, weight_samples_by_lnl, cache_fn)
    135 post_mask = get_post_mask(post, refi)
    136 fi_kwargs = dict(
    137     posterior_samples=post[post_mask],
    138     ref_samp=post[refi],
    139     ref_lnpri=ln_pri[refi],
    140     ref_lnl=ln_lnl[refi],
    141 )
--> 142 lnzs[i] = np.array([fi_ln_evidence(**fi_kwargs, r=ri) for ri in r_vals])
    143 median_lnzs[i] = np.nanmedian(lnzs[i])
    146 pbar.set_postfix_str(f"FI LnZ: {med_:.2f}")

File ~/Documents/projects/funnel/src/funnel/fi_core.py:142, in <listcomp>(.0)
    135 post_mask = get_post_mask(post, refi)
    136 fi_kwargs = dict(
    137     posterior_samples=post[post_mask],
    138     ref_samp=post[refi],
    139     ref_lnpri=ln_pri[refi],
    140     ref_lnl=ln_lnl[refi],
    141 )
--> 142 lnzs[i] = np.array([fi_ln_evidence(**fi_kwargs, r=ri) for ri in r_vals])
    143 median_lnzs[i] = np.nanmedian(lnzs[i])
    146 pbar.set_postfix_str(f"FI LnZ: {med_:.2f}")

File ~/Documents/projects/funnel/src/funnel/fi_core.py:62, in fi_ln_evidence(posterior_samples, ref_samp, r, ref_lnpri, ref_lnl)
     43 def fi_ln_evidence(
     44     posterior_samples: np.ndarray,
     45     ref_samp: np.array,
   (...)
     48     ref_lnl: float,
     49 ):
     50     """
     51     Returns the approx log-evidence of some posterior samples (using a reference parameter value).
     52     The approximation is based on the 'density estimation' method described in
   (...)
     60     :return: The log of the approximated log-evidence
     61     """
---> 62     approx_ln_post = fi_ln_posterior(posterior_samples, ref_samp, r)
     63     return ref_lnpri + ref_lnl - approx_ln_post

File ~/Documents/projects/funnel/src/funnel/fi_core.py:35, in fi_ln_posterior(posterior_samples, reference_sample, r)
     29 # patricio's suggestion
     30 # different r for each parameter
     31 # max_diff = np.nanmax(diff_from_ref, axis=0)
     32 # r = np.pi/max_diff
     33 # r is now a vector
     34 sin_diff = np.sin(r * diff_from_ref)
---> 35 sum_res = np.nansum(np.nanprod(sin_diff / diff_from_ref, axis=1))
     36 n_samp, n_dim = posterior_samples.shape
     37 const = 1 / (n_samp * np.power(np.pi, n_dim))

KeyboardInterrupt: 

plt_kwgs = dict(lnzs=lnzs, r_vals=r_vals,sampling_lnz=[GW150914_data.lnBF])
fig = plot_fi_evidence_results(**plt_kwgs)
fig.gca().set_ylabel("lnBF");
../_images/b5be6669a7be93f5af56b5e7f6e046a396e6f50dca268d645879a248c197769c.png 
fig = plot_fi_evidence_results(**plt_kwgs, plot_all_lnzs=True, plt_kwgs=dict(alpha=0.5, lw=2))
fig.gca().set_ylabel("lnBF");
../_images/e8f7a40f39912b1460fff107b9040632d03d364df27291df2f67e4ce50ab2559.png