Upload New File

%% Cell type:code id: tags:

``` python

import numpy as np

import matplotlib.pyplot as plt

```

%% Cell type:code id: tags:

``` python

# load files from 3 PE runs.

# The first m1 means mass ratio q=1.

# NT1 means the first noise test.

# m1,m2,d are three parameters.

m1_1 = np.loadtxt('m1_NT1_m1_posterior_samples.txt')

m2_1 = np.loadtxt('m1_NT1_m2_posterior_samples.txt')

d_1 = np.loadtxt('m1_NT1_d_posterior_samples.txt')

m1_2 = np.loadtxt('m1_NT2_m1_posterior_samples.txt')

m2_2 = np.loadtxt('m1_NT2_m2_posterior_samples.txt')

d_2 = np.loadtxt('m1_NT2_d_posterior_samples.txt')

m1_3 = np.loadtxt('m1_NT3_m1_posterior_samples.txt')

m2_3 = np.loadtxt('m1_NT3_m2_posterior_samples.txt')

d_3 = np.loadtxt('m1_NT3_d_posterior_samples.txt')

```

%% Cell type:code id: tags:

``` python

nsample_1 = len(m1_1)

nsample_1

```

%% Output

    33762

%% Cell type:code id: tags:

``` python

nsample_2 = len(m1_2)

nsample_2

```

%% Output

    28767

%% Cell type:code id: tags:

``` python

nsample_3 = len(m1_3)

nsample_3

```

%% Output

    32259

%% Cell type:code id: tags:

``` python

def jsd(p1_sample,p2_sample,nbin1,nbin2,nbin):

    # First we make histogram.

    hist_1, edge_1 = np.histogram(p1_sample, bins=nbin1)

    x_1 = np.linspace(0.5*(edge_1[0]+edge_1[1]),0.5*(edge_1[nbin1-1]+edge_1[nbin1]),nbin1)

    hist_2, edge_2 = np.histogram(p2_sample, bins=nbin2)

    x_2 = np.linspace(0.5*(edge_2[0]+edge_2[1]),0.5*(edge_2[nbin2-1]+edge_2[nbin2]),nbin2)

    # Then we make a common value axis.

    x_min = min(x_1.min(),x_2.min())

    x_max = min(x_1.max(),x_2.max())

    x = np.linspace(x_min,x_max,nbin)

    y_1 = np.linspace(0,0,nbin)

    y_2 = np.linspace(0,0,nbin)

    # Then we interpolate the curve.

    for i in range(len(x)):

        t=0

        for j in range(len(x_1)):

            if x[i]< x_1[j]:

                t=j

                break

        y_1[i] = 0.5*(hist_1[t-1]+hist_1[t])

    for i in range(len(x)):

        t=0

        for j in range(len(x_2)):

            if x[i]< x_2[j]:

                t=j

                break

        y_2[i] = 0.5*(hist_2[t-1]+hist_2[t])

    # Finally we calculate JS divergence.

    from scipy.spatial import distance

    js_divergence=distance.jensenshannon(y_1,y_2)

    return js_divergence

```

%% Cell type:code id: tags:

``` python

nbin1=100

nbin2=100

nbin=100

```

%% Cell type:code id: tags:

``` python

jsd_m1_12 = jsd(m1_1,m1_2,nbin1,nbin2,nbin)

jsd_m1_13 = jsd(m1_1,m1_3,nbin1,nbin2,nbin)

jsd_m2_12 = jsd(m2_1,m2_2,nbin1,nbin2,nbin)

jsd_m2_13 = jsd(m2_1,m2_3,nbin1,nbin2,nbin)

jsd_d_12 = jsd(d_1,d_2,nbin1,nbin2,nbin)

jsd_d_13 = jsd(d_1,d_3,nbin1,nbin2,nbin)

```

%% Output

    ---------------------------------------------------------------------------

    NameError                                 Traceback (most recent call last)

    <ipython-input-7-8af5791cf77c> in <module>

    ----> 1 jsd_m1_12 = jsd(m1_1,m1_2,nbin1,nbin2,nbin)

          2 jsd_m1_13 = jsd(m1_1,m1_3,nbin1,nbin2,nbin)

          3 jsd_m2_12 = jsd(m2_1,m2_2,nbin1,nbin2,nbin)

          4 jsd_m2_13 = jsd(m2_1,m2_3,nbin1,nbin2,nbin)

          5 jsd_d_12 = jsd(d_1,d_2,nbin1,nbin2,nbin)

    NameError: name 'nbin1' is not defined

%% Cell type:code id: tags:

``` python

jsd_m1_12

```

%% Output

    0.49853152880729623

%% Cell type:code id: tags:

``` python

jsd_m1_13

```

%% Output

    0.29887917498036515

%% Cell type:code id: tags:

``` python

jsd_m2_12

```

%% Output

    0.3507920042290908

%% Cell type:code id: tags:

``` python

jsd_m2_13

```

%% Output

    0.2921521879260992

%% Cell type:code id: tags:

``` python

jsd_d_12

```

%% Output

    0.6309011050754164

%% Cell type:code id: tags:

``` python

jsd_d_13

```

%% Output

    0.5018860973907884

%% Cell type:code id: tags:

``` python

```