Alternative error analysis for BICEP2 result

Kyle Cranmer [email protected] Mar 17, 2014

The BICEP2 error quoted on r is 0.2+0.07-0.05 (green below lines), and it looks like it's really a Bayesian shortest interval using a uniform prior on r. Let's see what the error bar looks like by using -2lnL < 1 (intersection of blue and red lines). The likelihood curve is shown here and they made the data public -- thanks BICEP2!

In [2]:
%pylab inline
Populating the interactive namespace from numpy and matplotlib
In [3]:
import requests
In [4]:
import numpy as np
from matplotlib import pyplot as plt
In [5]:
r=requests.get('http://bicepkeck.org/B2_3yr_rlikelihood_20140314.txt')
In [6]:
header = r.text.split('\n')[:12]

Convert the text into an array of floats.

In [7]:
rawLikelihoodData = r.text.split('\n')[13:]
In [8]:
likelihoodScan = []
temp = [i.split('\t') for i in rawLikelihoodData]
temp
for x,y in temp[:-1]:
    likelihoodScan.append([float(x), float(y)])
likelihoodScan = np.array(likelihoodScan)

And the integral is ~1 (so it's like a normalized posterior including the bin width ie. prob/0.001)

In [9]:
np.sum(likelihoodScan[:,1])
Out[9]:
0.99999982967601886

Get shortest interval, it will be from [0.145,0.262] to get the 68% for 1$\sigma$

In [10]:
lmax = likelihoodScan[:,1].max()
rpeak = np.where(likelihoodScan[:,1]==lmax)[0][0]
In [11]:
rminBicep, rmaxBicep = rpeak,rpeak
integral,cutoff=0.,lmax-lmax/1000
while integral<0.68:
    rRange = np.where(likelihoodScan[:,1]>cutoff)
    rminBicep=np.min(rRange)
    rmaxBicep=np.max(rRange)
    integral = np.sum(likelihoodScan[rminBicep:rmaxBicep,1])
    cutoff-=lmax/1000
print '[%f, %f]'%( rminBicep, rmaxBicep)
[145.000000, 262.000000]

Remake the likelihood curve shown here

In [12]:
plt.plot(likelihoodScan[:,0],likelihoodScan[:,1])
plt.plot(likelihoodScan[rminBicep:rmaxBicep,0],likelihoodScan[rminBicep:rmaxBicep,1],'g')
plt.plot([rminBicep/1000.,rminBicep/1000.],[0.,likelihoodScan[rminBicep,1]],'g')
plt.plot([rmaxBicep/1000.,rmaxBicep/1000.],[0.,likelihoodScan[rmaxBicep,1]],'g')
plt.xlabel('r')
plt.ylabel('Likelihood')
Out[12]:
<matplotlib.text.Text at 0x110b42b10>

Fortunately, a constant factor on the likelihood cancels out in likelihood ratio (or difference of log likelihoods). So let's see what the typical $\Delta \log L$ or $\Delta \chi^2$ interval would look like.

In [23]:
minLogL = -2 * np.log(np.max(likelihoodScan[:,1]))
rmin, rmax = 0,0
for r, l in likelihoodScan:
    if -2*np.log(l)-minLogL < 1:
        rmin= r
        break
        
for r, l in likelihoodScan[200:]:
    if -2*np.log(l)-minLogL > 1:
        rmax=r
        break
print '[%f, %f]'%( rmin, rmax)
[0.147000, 0.260000]
In [24]:
plt.plot(likelihoodScan[:,0],-2* np.log(likelihoodScan[:,1])-minLogL)
plt.plot([0,0.8],[1.,1.],'r')
plt.plot([rmaxBicep/1000.,rmaxBicep/1000.],[0.,1.],'g')
plt.plot([rminBicep/1000.,rminBicep/1000.],[0.,1.],'g')
plt.plot([rmax,rmax],[0.,1.],'r')
plt.plot([rmin,rmin],[0.,1.],'r')
plt.xlabel('r')
plt.ylabel('-2 ln Likleihood')
plt.ylim([0,2])
plt.xlim([0.1,0.3])
Out[24]:
(0.1, 0.3)

So that's slightly different from the BICEP2 result. In terms of the Bayesian credibility, that interval corresponds to

In [25]:
np.sum(likelihoodScan[147:260,1])
Out[25]:
0.66418630000000012

Slightly different numbers and a different probabilistic interpretation, but luckily it doesn't change the profound result.

Assuming that the BICEP2 data are in the asymptotic regieme, then we can quickly estimate the statistical significance of the observation (testing the r=0 hypothesis).

In [26]:
significance = np.sqrt(-2*np.log(likelihoodScan[0,1]))
In [27]:
significance
Out[27]:
7.6406188623296876

This was done during a Software Carpentry bootcamp at NYU, it's more of a fun python project than a critique of the BICEP2 result. Congratulations to BICEP2!