<jupyter_start><jupyter_text>Look Elsewhere Effect in 2-d<NAME>, Nov 19, 2015Based on*Estimating the significance of a signal in a multi-dimensional search* by  <NAME> and <NAME> http://arxiv.org/pdf/1105.4355v1.pdfThis is for the special case of a likelihood function of the form $L(\mu, \nu_1, \nu_2)$ where $\mu$ is a single parameter of interest and$\nu_1,\nu_2$ are two nuisance parameters that are not identified under the null.For example, $\mu$ is the signal strength of a new particle and $\nu_1,\nu_2$ are theunknown mass and width of the new particle. Under the null hypothesis, those parameters don't mean anything... aka they "are not identified under the null" in the statistics jargon.This introduces a 2-d look elsewhere effect.The LEE correction in this case is based on \begin{equation}E[ \phi(A_u) ] = P(\chi^2_1 > u) + e^{-u/2} (N_1 + \sqrt{u} N_2) \,\end{equation}where    * $A_u$ is the 'excursion set above level $u$ (eg. the set of parameter points in $(\nu_1,\nu_2)$ that have a -2 log-likelihood ratio greater than $u$ )   * $\phi(A_u)$ is the Euler characteristic of the excursion set   * $E[ \phi(A_u) ]$ is the expectation of the Euler characteristic of those excursion sets under the null   * $P(\chi^2_1 > u)$ is the standard chi-square probability    * and $N_1$ and $N_2$ are two coefficients that characterize the chi-square random field.    structure of the notebookThe notebook is broken into two parts.  * create a ROOT file with some toy histograms for this example  * read in those toy histograms used to estimate $N_1$ and $N_2$  * calculate the global p-value (eg. correct for the look-elsewhere effect)<jupyter_code>%pylab inline --no-import-all
from lee2d import *<jupyter_output><empty_output><jupyter_text>Create some test histogramsThese are toy histograms, but they are supposed to represent scans of $ q(\nu_1, \nu_2)$ where\begin{equation}q(\nu_1, \nu_2) = -2 \log \frac{ \max_\theta L(\mu=0, \nu_1, \nu_2, \theta)}{ \max_{\mu, \theta} L(\mu, \nu_1, \nu_2, \theta)}\end{equation}and $\theta$ are nuisance parameters for the background model, $\mu$ is parameter of interest, and $\nu_1, \nu_2$ are nuisance parameters for the signal that are not meaningful for the null (eg. the mass and width of a hypothetical particle).In this setting, the maximum local significance is given by \begin{equation}Z_{local} = \max_{\nu_1, \nu_2} \sqrt{q(\nu_1, \nu_2)}\end{equation}<jupyter_code>from create_test_histograms import *
create_test_histograms()
#check to make sure test_hists.root exists
!ls *root<jupyter_output>HiggsLWA.root    HiggsLWANew.root test_hists.root<jupyter_text>Calculate global p-value based on maximum local significance and several toy likelihood scans<jupyter_code># Specify the necessary info
max_local_significance = 4. 
# choose u1, u2 thresholds for doing the scan.
# these are arbitrary 
# if there are enough toys the choice shouldn't matter, but 
# we may want to do some tests with other choices
this_u1, this_u2 = 0.1, 0.9
# Specify the root file with the histograms and their names
root_file_name = 'test_hists.root'
num_toy_scans = 25
names_of_toy_likleihood_scans = [('scan_toy_%d' %(i)) for i in range(11,num_toy_scans)]
def convert_hist_to_numpy(hist):
    """a little helper script"""
    temp = np.zeros((hist.GetNbinsX(), hist.GetNbinsY()))
    for i in range(temp.shape[0]):
        for j in range(temp.shape[1]):
            temp[i,j] = hist.GetBinContent(i+1, j+1)
    return temp
# Read in histograms, convert them to numpy arrays
inFile = ROOT.TFile(root_file_name, 'READ')
likelihoodScans = []
for histName in names_of_toy_likleihood_scans:
    inHist = inFile.Get(histName)
    temp = convert_hist_to_numpy(inHist)
    likelihoodScans.append(temp)   
from scipy.ndimage import grey_closing, binary_closing
def fill_holes(array):
    zero_array = array==0.
    temp = grey_closing(array, size=2)*zero_array
    return temp+array
def get_euler_characteristics(listOfScans, u1=0.1, u2=0.9):
    """
    loop through the likleihood scans and calculate expectation
    of Euler characteristic for excursion sets above levels u1, u2
    """
    
    n_plots = 3
    plt.figure(figsize=(9,n_plots*3))
    phis = np.zeros((len(listOfScans),2))
    for scan_no, scan in enumerate(listOfScans):

        # fill holes from failures in original likelihood
        scan = fill_holes(scan)
        
        #get excursion sets above those two levels
        exc1 = (scan>u1) + 0. #add 0. to convert from bool to double
        exc2 = (scan>u2) + 0.
        #print '\nu1,u2 = ', u1, u2

        if scan_no < n_plots:
            plt.subplot(n_plots,3,3*scan_no+1)
            aspect = 1.*scan.shape[0]/scan.shape[1]
            plt.imshow(scan.T, cmap='gray', aspect=aspect)
            plt.subplot(n_plots,3,3*scan_no+2)
            plt.imshow(exc1.T, cmap='gray', aspect=aspect)
            plt.subplot(n_plots,3,3*scan_no+3)
            plt.imshow(exc2.T, cmap='gray', aspect=aspect)
        
        phi1 = calculate_euler_characteristic(exc1)
        phi2 = calculate_euler_characteristic(exc2)
        #print 'phi1, phi2 = ', phi1, phi2
        phis[scan_no] = [phi1, phi2]

    plt.savefig('islands.png')
    print 'Exp phi_0=%f, phi_2=%f' %(mean(phis[:,0]), mean(phis[:,1]))
    return mean(phis[:,0]), mean(phis[:,1])<jupyter_output><empty_output><jupyter_text>The columns of the below are for $\left| \hspace{2em} q(\nu_1, \nu_2) \hspace{2em}\right| \hspace{2em}q(\nu_1, \nu_2)>u_1 \hspace{2em} \left| \hspace{2em} q(\nu_1, \nu_2) > u_2 \hspace{2em} \right|$<jupyter_code>expphi1, expphi2 = get_euler_characteristics(likelihoodScans, u1=this_u1, u2=this_u2)
global_p_value = do_LEE_correction(max_local_significance, this_u1, this_u2, expphi1, expphi2)<jupyter_output>n1, n2 = 3.15777649996 3.66198277791
 local p_value = 0.000032,  local significance = 4.000000
global p_value = 0.006036, global significance = 2.510004
{'fraction_non_alphanumeric': 0.1145374449339207, 'fraction_numerical': 0.03456455438834293, 'mean_word_length': 4.0281090289608175, 'pattern_counts': {'":': 0, '<': 19, '<?xml version=': 0, '>': 24, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '434437', 'license': 'MIT', 'n_tokens_mistral': 2159, 'n_tokens_neox': 2026, 'n_words': 620, 'path': 'root_TH2F_example.ipynb', 'repo_name': 'cranmer/look-elsewhere-2d', 'size': '5919'}