{ "cells": [ { "cell_type": "code", "execution_count": 1, "id": "2b23f073", "metadata": {}, "outputs": [], "source": [ "# Program to find measure of expected significance as a function\n", "# of a cut value x_cut applied to measured variable x\n", "# by Monte Carlo simulation of the likelihood ratio statistic.\n", "# G. Cowan / RHUL Physics / December 2022\n", "\n", "import numpy as np\n", "import scipy.stats as stats\n", "import matplotlib\n", "import matplotlib.pyplot as plt\n", "plt.rcParams[\"font.size\"] = 14" ] }, { "cell_type": "code", "execution_count": 2, "id": "3b565a2b", "metadata": {}, "outputs": [], "source": [ "# Define pdfs and likelihood-ratio statisic\n", "s_tot = 10.\n", "b_tot = 100.\n", "ps = s_tot/(s_tot+b_tot)\n", "def f_s(x):\n", " return 3.*(1.-x)**2\n", "def f_b(x):\n", " return 3.*x**2\n", "def q(x):\n", " return -2.*np.log(1. + (s_tot/b_tot)*f_s(x)/f_b(x))" ] }, { "cell_type": "code", "execution_count": 3, "id": "5880ac4d", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "....................................................................................................\n", "\n" ] } ], "source": [ "# Generate data under b and s+b hypotheses\n", "qb = []\n", "qsb = []\n", "numExp = 1000000\n", "np.random.seed(seed=1234567) # fix random seed\n", "for i in range(numExp):\n", " n = np.random.poisson(b_tot) # first b only\n", " r = np.random.uniform(0., 1., n)\n", " xb = r**(1./3.)\n", " qb.append(np.sum(q(xb)))\n", " n = np.random.poisson(s_tot+b_tot) # then s+b\n", " r1 = np.random.uniform(0., 1., n)\n", " r2 = np.random.uniform(0., 1., n)\n", " xsb = [1. - r1[j]**(1./3.) if r2[j]" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Plot histograms of q\n", "binLo, binHi = bin_edges[:-1], bin_edges[1:]\n", "xPlot = np.array([binLo, binHi]).T.flatten()\n", "ybPlot = np.array([qbHist, qbHist]).T.flatten()\n", "ysbPlot = np.array([qsbHist, qsbHist]).T.flatten()\n", "fig, ax = plt.subplots(1,1)\n", "plt.yscale(\"log\")\n", "plt.gcf().subplots_adjust(bottom=0.15)\n", "plt.gcf().subplots_adjust(left=0.15)\n", "plt.xlabel(r'$q$', labelpad=5)\n", "plt.ylabel(r'$f(q)$', labelpad=10)\n", "plt.plot(xPlot, ybPlot, label=r'$f(q|b)$', color='dodgerblue')\n", "plt.plot(xPlot, ysbPlot, label=r'$f(q|s+b)$', color='orange')\n", "ax.axvline(med_q_sb, color=\"black\", linestyle=\"dashed\", label = r'median$[q|s+b]$')\n", "ax.legend(loc='upper left', frameon=False)\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": 6, "id": "fb2a1953", "metadata": {}, "outputs": [], "source": [ "# Add code to calculate the p-value of the b-only hypothesis\n", "# for the median q from data generated according to s+b.\n", "# Find the corresponding significance Z." ] }, { "cell_type": "code", "execution_count": null, "id": "edb93068", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.8" } }, "nbformat": 4, "nbformat_minor": 5 }