import pymc3 as pm
import numpy as np
import theano

from theano.compile.ops import as_op
def background_rate_f(b, T, n):
    """
    
    """
    out = 0
    #n = int(n)
    for i in range(n+1):
        out += ((b*T)**i * np.exp(- b*T)) / np.math.factorial(i)
    return out

def log_background_rate(b, T, n):
    return np.log(background_rate_f(b, T, n))

def signal_rate_part(s, n, b, T):
    top_a = T * ((s + b) * T)**n 
    top_b = np.exp(-(s + b)*T)
    p = (top_a * top_b) / np.math.factorial(n)
    return theano.tensor.switch(theano.tensor.le(s, 0), 0, p)

#@as_op(itypes=[T.dscalar, T.dscalar, T.dscalar, T.dscalar], otypes=[T.dscalar])
def log_signal_rate(s,n,b,T):
    #if theano.tensor.lt(0, s): return np.array([[0.0]])
    p = -log_background_rate(b,T,n) + np.log(signal_rate_part(s,n,b,T))
    
    return p
def number_mweg(horizon):
    return 4./3 * np.pi * horizon**3 *(2.26)**-3  #* (0.0116) #* horizon**3
import theano.tensor as T
from pymc3 import DensityDist, Uniform, Normal
from pymc3 import Model
from pymc3 import distributions


with Model() as model:
    number_events = 0
    background_rate = 0.01 #
    observation_time = 46.1/365. # 
    duty_cycle = 0.2
    horizon = 73.2    
    grb_rate = 1.0
    
    #trunc = pm.Bound(pm.DensityDist, lower=0)
    signal_rate = pm.DensityDist('signal_rate', 
                        logp=lambda value: log_signal_rate(value, number_events, background_rate, observation_time),
                       testval=10)
    
    n_galaxy = number_mweg(horizon)
    
    cbc_rate = pm.Deterministic('cbc_rate', signal_rate * n_galaxy * horizon)
    
    efficiency = pm.Uniform('efficiency', 0,1)
   
    costheta = pm.Deterministic('cos_angle', 1.0 - ((grb_rate/cbc_rate*efficiency)))
    
    angle = pm.Deterministic("angle", theano.tensor.arccos(costheta))
import theano.tensor as T
from pymc3 import DensityDist, Uniform, Normal
from pymc3 import Model
from pymc3 import distributions

def grb_model(number_events, background_rate, 
              observation_time, horizon, grb_rate,
             efficiency_prior = "uniform"):
    with Model() as model:
        signal_rate = pm.DensityDist('signal_rate', 
                            logp=lambda value: log_signal_rate(value, number_events, background_rate, observation_time),
                           testval=50)

        n_galaxy = number_mweg(horizon)

        cbc_rate = pm.Deterministic('cbc_rate', signal_rate * n_galaxy)

        # Allow the efficiency prior to be switched-out
        if efficiency_prior == "uniform":
            efficiency = pm.Uniform('efficiency', 0,1)
        elif efficiency_prior == "jeffreys":
            efficiency = pm.Beta('efficiency', 0.5, 0.5)
        elif isinstance(efficiency_prior, float):
            efficiency = efficiency_prior

        costheta = pm.Deterministic('cos_angle', 1.0 - ((grb_rate/cbc_rate*efficiency)))

        angle = pm.Deterministic("angle", theano.tensor.arccos(costheta))
        
        return model
priors = ["uniform", "jeffreys", 1.0, 0.5]
# O1 Scenarios
number_events = 0 # There were no BNS detections in O1
background_rate = 0.01 # We take the FAR to be 1/100 yr
observation_time = 46.1/365. # The number of days of analysis conducted by gstLAL
horizon = 73.2  # The O1 BNS horizon distance in O1 BNS paper
grb_rate = 10.0
o1_models = []
for prior in priors:
    o1_models.append( grb_model(number_events, background_rate, observation_time, horizon, grb_rate, prior))
# O2 Scenarios
number_events = 1 # Assume O2 will see one BNS observation
background_rate = 0.01 # We take the FAR to be 1/100 yr
observation_time = 0.5 # The number of days of analysis conducted by gstLAL
duty_cycle = 0.4
observation_time *= duty_cycle
horizon = 80.0  # The O1 BNS horizon distance in O1 BNS paper
grb_rate = 10.0
o2_models = []
for prior in priors:
    o2_models.append( grb_model(number_events, background_rate, observation_time, horizon, grb_rate, prior))
samples = 1000000
o1_traces = []
for model in o1_models:
    with model:
        step = pm.Metropolis()
        o1_traces.append(pm.sample(samples, step))
    
o2_traces = []    
for model in o2_models:
    with model:
        step = pm.Metropolis()
        o2_traces.append(pm.sample(samples, step))
 100%|██████████| 1000000/1000000 [02:34<00:00, 6473.32it/s]
100%|██████████| 1000000/1000000 [02:39<00:00, 6270.46it/s]
100%|██████████| 1000000/1000000 [01:17<00:00, 12977.07it/s]
100%|██████████| 1000000/1000000 [01:18<00:00, 12661.71it/s]
100%|██████████| 1000000/1000000 [02:27<00:00, 6763.74it/s]
100%|██████████| 1000000/1000000 [02:36<00:00, 6376.05it/s]
100%|██████████| 1000000/1000000 [01:16<00:00, 13125.98it/s]
100%|██████████| 1000000/1000000 [01:17<00:00, 12848.71it/s]

import matplotlib.pyplot as plt
%matplotlib inline
plt.style.use("/home/daniel/papers/thesis/thesis-style.mpl")
priors = ["U(0,1)", "jeffreys", "$\delta(1)$", "$\delta(0.5)$"]
f, ax = plt.subplots(4,len(priors), figsize = (15, 15))
for i in range(len(priors)):

    o1_trace = o1_traces[i]
    o2_trace = o2_traces[i]

    ax[0,i].hist(o1_trace[2000:]['signal_rate'], histtype="step", label="O1 Signal Rate",
                 lw=2, normed=True, bins=100);
    ax[0,i].hist(o2_trace[2000:]['signal_rate'], histtype="step", label="O2 Signal Rate",
                 lw=2, normed=True, bins=100);
    #ax[0,0].vlines(np.percentile(trace['signal_rate'][2000:], 90), 0, 0.085,
    #              label = "90% lower limit");
    ax[0,i].set_title("$\epsilon \sim$ {}".format(priors[i]));
    ax[0,i].set_xlabel("$s$ /yr")
    if i==3:ax[0,i].legend()

    ax[1,i].hist(o1_trace[2000:]['cbc_rate'], histtype="step", lw=2, bins=100,
                 normed=True, label="O1 CBC Rate");
    ax[1,i].hist(o2_trace[2000:]['cbc_rate'], histtype="step", lw=2, bins=100,
                 normed=True, label="O2 CBC Rate");
    ax[1,i].set_xlabel(u"$s \mathrm{/Gpc}^3\mathrm{/yr}$");
    #ax[1,0].vlines(np.percentile(o1_trace['cbc_rate'][2000:], 90), 0, 5e-5,
    #              label = "90% lower limit, O1");
    #ax[1,0].vlines(np.percentile(o2_trace['cbc_rate'][2000:], 90), 0, 5e-5,
    #              label = "90% lower limit, O2", color="orange");
    if i==3:ax[1,i].legend()
    try:
        ax[2,i].hist(o1_trace[2000:]['efficiency'], histtype="step", lw=2, label="O1 Efficiency");
        ax[2,i].hist(o2_trace[2000:]['efficiency'], histtype="step", lw=2, label="O2 Efficiency");
        #ax[1,0].set_xlabel(u"$s$ /Gpc³/yr")
        if i==0: ax[2,i].legend()
        ax[2,i].set_ylim(0,2.5e4)
    except KeyError:
            ax[2,i].spines['top'].set_visible(False)
            ax[2,i].spines['bottom'].set_visible(False)
            ax[2,i].spines['left'].set_visible(False)
            ax[2,i].spines['right'].set_visible(False)
            ax[2,i].grid(False)
            ax[2,i].set_yticklabels([])
            ax[2,i].set_xticklabels([])
            ax[2,i].yaxis.set_ticks_position('none')
            ax[2,i].xaxis.set_ticks_position('none')

    ax[3,i].hist(np.rad2deg(o1_trace[2000:]['angle'][np.isfinite(o1_trace[2000:]['angle'])]),
               histtype="step", lw=2, label="O1 Angle", log=True);
    ax[3,i].hist(np.rad2deg(o2_trace[2000:]['angle'][np.isfinite(o2_trace[2000:]['angle'])]),
               histtype="step", lw=2, label="O2 Angle", log=True);
    #ax[3,i].vlines(np.percentile(np.rad2deg(o1_trace[2000:]['angle'][np.isfinite(o1_trace[2000:]['angle'])]), 90), 0, 3e5,
    #              label = "90\% lower limit, O1");
    #ax[3,i].vlines(np.percentile(np.rad2deg(o2_trace[2000:]['angle'][np.isfinite(o2_trace[2000:]['angle'])]), 90), 0, 3e5,
    #              label = "90\% lower limit, O2");
    #ax[1,0].set_xlabel(u"$s$ /Gpc³/yr")
    if i==3: ax[3,i].legend()
    if i > 0: 
        ax[3,i].set_yticklabels([]), ax[2,i].set_yticklabels([])
        ax[0,i].set_yticklabels([]), ax[1,i].set_yticklabels([])
    ax[3,i].set_xlabel("Beaming angle [deg]")    
f.subplots_adjust(0.07, 0.05, .99, .95, wspace=0.05)
#f.savefig("/home/daniel/papers/thesis/figures/grbbeamingpymc.pdf")
png

O1 Trace Summaries

from pymc3 import summary

for i, trace in enumerate(o1_traces):
    print "# {}".format(priors[i])
    print summary(trace)
 # U(0,1)

signal_rate:

  Mean             SD               MC Error         95% HPD interval
  -------------------------------------------------------------------
  
  7.950            7.919            0.024            [0.000, 23.735]

  Posterior quantiles:
  2.5            25             50             75             97.5
  |--------------|==============|==============|--------------|
  
  0.203          2.287          5.513          11.033         29.261


cbc_rate:

  Mean             SD               MC Error         95% HPD interval
  -------------------------------------------------------------------
  
  1131505.753      1127094.477      3466.742         [3.684, 3378192.102]

  Posterior quantiles:
  2.5            25             50             75             97.5
  |--------------|==============|==============|--------------|
  
  28821.999      325449.128     784717.103     1570299.722    4164778.112


efficiency:

  Mean             SD               MC Error         95% HPD interval
  -------------------------------------------------------------------
  
  0.501            0.289            0.001            [0.008, 0.958]

  Posterior quantiles:
  2.5            25             50             75             97.5
  |--------------|==============|==============|--------------|
  
  0.025          0.250          0.502          0.752          0.975


cos_angle:

  Mean             SD               MC Error         95% HPD interval
  -------------------------------------------------------------------
  
  1.000            0.002            0.000            [1.000, 1.000]

  Posterior quantiles:
  2.5            25             50             75             97.5
  |--------------|==============|==============|--------------|
  
  1.000          1.000          1.000          1.000          1.000


angle:

  Mean             SD               MC Error         95% HPD interval
  -------------------------------------------------------------------
  
  0.005            0.008            0.000            [0.000, 0.013]

  Posterior quantiles:
  2.5            25             50             75             97.5
  |--------------|==============|==============|--------------|
  
  0.001          0.002          0.003          0.005          0.019

None
# jeffreys

signal_rate:

  Mean             SD               MC Error         95% HPD interval
  -------------------------------------------------------------------
  
  7.916            7.928            0.028            [0.000, 23.676]

  Posterior quantiles:
  2.5            25             50             75             97.5
  |--------------|==============|==============|--------------|
  
  0.201          2.284          5.489          10.973         29.248


cbc_rate:

  Mean             SD               MC Error         95% HPD interval
  -------------------------------------------------------------------
  
  1126724.550      1128323.743      3986.330         [2.938, 3369861.014]

  Posterior quantiles:
  2.5            25             50             75             97.5
  |--------------|==============|==============|--------------|
  
  28624.764      325012.929     781314.584     1561783.717    4162849.244


efficiency:

  Mean             SD               MC Error         95% HPD interval
  -------------------------------------------------------------------
  
  0.500            0.353            0.001            [0.006, 1.000]

  Posterior quantiles:
  2.5            25             50             75             97.5
  |--------------|==============|==============|--------------|
  
  0.002          0.147          0.500          0.854          0.998


cos_angle:

  Mean             SD               MC Error         95% HPD interval
  -------------------------------------------------------------------
  
  1.000            0.005            0.000            [1.000, 1.000]

  Posterior quantiles:
  2.5            25             50             75             97.5
  |--------------|==============|==============|--------------|
  
  1.000          1.000          1.000          1.000          1.000


angle:

  Mean             SD               MC Error         95% HPD interval
  -------------------------------------------------------------------
  
  nan              nan              nan              [0.000, nan]

  Posterior quantiles:
  2.5            25             50             75             97.5
  |--------------|==============|==============|--------------|
  
  0.000          0.002          0.003          0.005          0.019

None
# $\delta(1)$

signal_rate:

  Mean             SD               MC Error         95% HPD interval
  -------------------------------------------------------------------
  
  7.893            7.920            0.026            [0.000, 23.623]

  Posterior quantiles:
  2.5            25             50             75             97.5
  |--------------|==============|==============|--------------|
  
  0.195          2.266          5.461          10.919         29.107


cbc_rate:

  Mean             SD               MC Error         95% HPD interval
  -------------------------------------------------------------------
  
  1123373.830      1127308.384      3663.674         [37.517, 3362282.283]

  Posterior quantiles:
  2.5            25             50             75             97.5
  |--------------|==============|==============|--------------|
  
  27813.904      322500.459     777304.749     1554127.145    4142861.795


cos_angle:

  Mean             SD               MC Error         95% HPD interval
  -------------------------------------------------------------------
  
  1.000            0.003            0.000            [1.000, 1.000]

  Posterior quantiles:
  2.5            25             50             75             97.5
  |--------------|==============|==============|--------------|
  
  1.000          1.000          1.000          1.000          1.000


angle:

  Mean             SD               MC Error         95% HPD interval
  -------------------------------------------------------------------
  
  0.007            0.012            0.000            [0.002, 0.019]

  Posterior quantiles:
  2.5            25             50             75             97.5
  |--------------|==============|==============|--------------|
  
  0.002          0.004          0.005          0.008          0.027

None
# $\delta(0.5)$

signal_rate:

  Mean             SD               MC Error         95% HPD interval
  -------------------------------------------------------------------
  
  7.910            7.930            0.028            [0.000, 23.702]

  Posterior quantiles:
  2.5            25             50             75             97.5
  |--------------|==============|==============|--------------|
  
  0.199          2.275          5.484          10.946         29.235


cbc_rate:

  Mean             SD               MC Error         95% HPD interval
  -------------------------------------------------------------------
  
  1125793.043      1128630.389      4047.984         [11.412, 3373502.460]

  Posterior quantiles:
  2.5            25             50             75             97.5
  |--------------|==============|==============|--------------|
  
  28322.959      323735.745     780507.886     1557943.353    4161086.755


cos_angle:

  Mean             SD               MC Error         95% HPD interval
  -------------------------------------------------------------------
  
  1.000            0.002            0.000            [1.000, 1.000]

  Posterior quantiles:
  2.5            25             50             75             97.5
  |--------------|==============|==============|--------------|
  
  1.000          1.000          1.000          1.000          1.000


angle:

  Mean             SD               MC Error         95% HPD interval
  -------------------------------------------------------------------
  
  0.005            0.009            0.000            [0.001, 0.013]

  Posterior quantiles:
  2.5            25             50             75             97.5
  |--------------|==============|==============|--------------|
  
  0.002          0.003          0.004          0.006          0.019

None

O2 Trace Summaries

from pymc3 import summary

for i, trace in enumerate(o2_traces):
    print "# {}".format(priors[i])
    summary(trace)
 # U(0,1)

signal_rate:

  Mean             SD               MC Error         95% HPD interval
  -------------------------------------------------------------------
  
  9.980            7.064            0.020            [0.205, 23.766]

  Posterior quantiles:
  2.5            25             50             75             97.5
  |--------------|==============|==============|--------------|
  
  1.198          4.782          8.389          13.435         27.788


cbc_rate:

  Mean             SD               MC Error         95% HPD interval
  -------------------------------------------------------------------
  
  1854225.165      1312543.607      3629.764         [38137.171, 4415671.052]

  Posterior quantiles:
  2.5            25             50             75             97.5
  |--------------|==============|==============|--------------|
  
  222492.665     888475.131     1558557.412    2496206.608    5162875.206


efficiency:

  Mean             SD               MC Error         95% HPD interval
  -------------------------------------------------------------------
  
  0.500            0.289            0.001            [0.013, 0.962]

  Posterior quantiles:
  2.5            25             50             75             97.5
  |--------------|==============|==============|--------------|
  
  0.025          0.250          0.501          0.750          0.975


cos_angle:

  Mean             SD               MC Error         95% HPD interval
  -------------------------------------------------------------------
  
  1.000            0.000            0.000            [1.000, 1.000]

  Posterior quantiles:
  2.5            25             50             75             97.5
  |--------------|==============|==============|--------------|
  
  1.000          1.000          1.000          1.000          1.000


angle:

  Mean             SD               MC Error         95% HPD interval
  -------------------------------------------------------------------
  
  0.003            0.002            0.000            [0.000, 0.006]

  Posterior quantiles:
  2.5            25             50             75             97.5
  |--------------|==============|==============|--------------|
  
  0.001          0.002          0.002          0.003          0.007

# jeffreys

signal_rate:

  Mean             SD               MC Error         95% HPD interval
  -------------------------------------------------------------------
  
  9.990            7.066            0.021            [0.221, 23.834]

  Posterior quantiles:
  2.5            25             50             75             97.5
  |--------------|==============|==============|--------------|
  
  1.211          4.806          8.381          13.442         27.819


cbc_rate:

  Mean             SD               MC Error         95% HPD interval
  -------------------------------------------------------------------
  
  1856020.465      1312804.485      3915.780         [41031.631, 4428279.026]

  Posterior quantiles:
  2.5            25             50             75             97.5
  |--------------|==============|==============|--------------|
  
  224933.794     892954.709     1557054.878    2497383.629    5168646.346


efficiency:

  Mean             SD               MC Error         95% HPD interval
  -------------------------------------------------------------------
  
  0.501            0.354            0.001            [0.000, 0.994]

  Posterior quantiles:
  2.5            25             50             75             97.5
  |--------------|==============|==============|--------------|
  
  0.002          0.148          0.502          0.854          0.999


cos_angle:

  Mean             SD               MC Error         95% HPD interval
  -------------------------------------------------------------------
  
  1.000            0.000            0.000            [1.000, 1.000]

  Posterior quantiles:
  2.5            25             50             75             97.5
  |--------------|==============|==============|--------------|
  
  1.000          1.000          1.000          1.000          1.000


angle:

  Mean             SD               MC Error         95% HPD interval
  -------------------------------------------------------------------
  
  0.003            0.002            0.000            [0.000, 0.006]

  Posterior quantiles:
  2.5            25             50             75             97.5
  |--------------|==============|==============|--------------|
  
  0.000          0.001          0.002          0.003          0.007

# $\delta(1)$

signal_rate:

  Mean             SD               MC Error         95% HPD interval
  -------------------------------------------------------------------
  
  9.957            7.045            0.020            [0.217, 23.716]

  Posterior quantiles:
  2.5            25             50             75             97.5
  |--------------|==============|==============|--------------|
  
  1.216          4.786          8.363          13.389         27.734


cbc_rate:

  Mean             SD               MC Error         95% HPD interval
  -------------------------------------------------------------------
  
  1849875.187      1308929.760      3778.905         [40407.083, 4406231.756]

  Posterior quantiles:
  2.5            25             50             75             97.5
  |--------------|==============|==============|--------------|
  
  225915.338     889216.982     1553870.163    2487656.769    5152789.812


cos_angle:

  Mean             SD               MC Error         95% HPD interval
  -------------------------------------------------------------------
  
  1.000            0.000            0.000            [1.000, 1.000]

  Posterior quantiles:
  2.5            25             50             75             97.5
  |--------------|==============|==============|--------------|
  
  1.000          1.000          1.000          1.000          1.000


angle:

  Mean             SD               MC Error         95% HPD interval
  -------------------------------------------------------------------
  
  0.004            0.002            0.000            [0.002, 0.008]

  Posterior quantiles:
  2.5            25             50             75             97.5
  |--------------|==============|==============|--------------|
  
  0.002          0.003          0.004          0.005          0.009

# $\delta(0.5)$

signal_rate:

  Mean             SD               MC Error         95% HPD interval
  -------------------------------------------------------------------
  
  10.019           7.100            0.018            [0.192, 23.869]

  Posterior quantiles:
  2.5            25             50             75             97.5
  |--------------|==============|==============|--------------|
  
  1.202          4.803          8.394          13.507         27.934


cbc_rate:

  Mean             SD               MC Error         95% HPD interval
  -------------------------------------------------------------------
  
  1861419.908      1319102.530      3375.802         [35760.213, 4434726.009]

  Posterior quantiles:
  2.5            25             50             75             97.5
  |--------------|==============|==============|--------------|
  
  223386.568     892339.046     1559636.760    2509531.379    5190025.393


cos_angle:

  Mean             SD               MC Error         95% HPD interval
  -------------------------------------------------------------------
  
  1.000            0.000            0.000            [1.000, 1.000]

  Posterior quantiles:
  2.5            25             50             75             97.5
  |--------------|==============|==============|--------------|
  
  1.000          1.000          1.000          1.000          1.000


angle:

  Mean             SD               MC Error         95% HPD interval
  -------------------------------------------------------------------
  
  0.003            0.002            0.000            [0.001, 0.006]

  Posterior quantiles:
  2.5            25             50             75             97.5
  |--------------|==============|==============|--------------|
  
  0.001          0.002          0.003          0.003          0.007