{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "#
Particle Physics II - Physics Beyond the Standard Model
\n", "\n", "##
Exercise sheet 1
\n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# import material\n", "from material import neyman_construction, tram_time, \\\n", " poisson_cdf, find_limit, find_cls_limit\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Confidence intervals and limits" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Confidence intervals quantify our degree of belief in the value of a given parameter $\\mu$ to lie within certain boundaries in the possible value space of $\\mu$. They play a crucial role in all physics measurements. We distinguish _double sided_ confidence intervals (for point estimates) from _single sided_ confidence intervals (equivalent to _upper_ or _lower boundaries_ on $\\mu$).\n", "\n", "
\"Drawing\"
\n", "\n", "In particle physics, we often call upper or lower boundaries _limits_. With the \"upper limit on $\\mu$\" we have, e.g., the following scenario in mind: \n", "\n", " * We consider a counting experiment, for which we expect a new signal in addition to an already known background.\n", " * Such a scenario corresponds to a \"two hypothesis test\" in statistics; we compare the (established) $H_{0}$ \"background-only\" hypothesis with the (alternative) $H_{1}$ \"signal plus background\" hypothesis. For this purpose we define $\\mu\\geq0$ to be the strength of the new signal;\n", " * Imagine the outcome of the experiment with an observation of $N_\\mathrm{obs}$.\n", " * We usually set a limit when opting for the $H_{0}$ hypothesis, i.e., $N_\\mathrm{obs}$ can be explained by the known background only, no extra signal is needed to explain such an outcome of the experiment.\n", " * With the limit we answer the following question: **What is the largest value of $\\mu$, for which the probability to observe $N_\\mathrm{obs}$ is smaller than a value $\\beta$?** \n", "\n", "We denote this limit as $\\mu_{1-\\beta}$, you might know $\\beta$ as the error of second kind, for the classical two hypothesis test. The value $1-\\beta$ we call _confidence level_ (CL). The idea is that a signal with strength $\\mu\\geq\\mu_{1-\\beta}$ we would **not** have missed with our experiment with a \"degree of belief\" (i.e. \"a confidence\") of $1-\\beta$. In particle physics we very often choose $\\beta=0.05$ (in astroparticle physics also $\\beta=0.10$ is used). We will thus refer to the (upper) limit as $\\mu_{0.95}$, which we call \"an (upper) limit at 95% CL\".\n", "\n", "When searching for new physics we usually set upper limits, e.g. on signal strengths, masses of hypothesized new particles, or energy scales, while this does not have to be the case. Check what kind of limit we are setting in **Exercise 1**." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Exercise 1: Limits in daily life**\n", "\n", "Imagine you are new in Karlsruhe planning your way to KIT by tram. You take a huge number of test rides and measure the time $t$ required to get from your home to KIT. You observe, that it is very unlikely that the tram is ahead of time, the ride rather takes a little longer than scheduled, e.g., due to traffic disruptions. NB: For this reason, we have modeled the probability distribution of $t$ using an asymmetric LogNormal density distribution. To be in time for your lecture and exams, you want to determine the time needed to get from home to KIT in order to find the latest point in time to leave home. In the following we are going to measure $t$ with negative values, i.e. $t=-10\\,\\mathrm{min}$ refers to \"10 mins before lecture (exam) starts\". You want to be in time for your lecture in 95% of all rides, i.e. at 95% confidence ($t_{0.95}$). For your exams you are adding some contingency, to be in time with 99% confidence ($t_{0.99}$).\n", " \n", "a) Comparing $t_{0.95}$ with $t_{0.99}$, which one do you expect to be larger (note the sign of $t$)? \n", "\n", "b) Use the function `tram_time(n_rides=number_of_test_rides, cl=CL, type='single')`, provided by the module _material.py_, to determine the exact values of $t_{0.95}$ and $t_{0.99}$. \n", "\n", "c) When you leave home at the times you calculated for lecture and exam in (b), what is the probability of being late for a lecture (exam)? \n", "\n", "In a second scenario, a fellow student wants to move to the same city district that you are living in. He wants to know how long it usually takes to get to university by tram. To give him a good idea of the distribution $t$ you provide him with a double sided range covering most of your rides.\n", "\n", "d) If you compare the width of a confidence interval $\\Delta t_{0.95}$ that covers 95% of your rides with the confidence interval $\\Delta t_{0.90}$, which one do you expect to be larger? \n", "\n", "e) Use the function `tram_time(n_rides=number_of_test_rides, cl=CL, type='double')` to determine $\\Delta t_{0.95}$ and $\\Delta t_{0.90}$. \n", "\n", "f) Comparing the upper bound of $\\Delta t_{0.95}$ with $t_{0.95}$, which one is smaller and why (note the sign of $t$)?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "tram_time(n_rides=??, cl=??, type='single')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "g) Determine $t_{0.95}$ for your test rides, but with different number_of_test_rides (see table below). Determine the standard deviation $\\sigma$ of this estimate from an ensembe test of length 10000. The argument `type='ensemble'` in the `upper_limit=tram_time(n_rides=number_of_test_rides, cl=0.95, type='ensemble')` function will give you $t_{0.95}$. What happens if you increase number_of_test_rides? \n", "\n", "| number_of_test_rides | $\\sigma$ |\n", "| --- | --- |\n", "| 10000 | |\n", "| 100 | |\n", "| 10 | |" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "\n", "for i in range(10000):\n", " upper_limit = tram_time(n_rides=??, cl=??, type='ensemble')\n", "\n", "print('Standard Deviation: ', ??)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Neyman construction for model parameters\n", "\n", "We will now discuss the more particle-physics related example of an upper limit on $\\mu$, that has been given above, in more detail. For a counting experiment, we can assume that the number of counts $k$ is distributed according to the Poisson distribution \n", "\n", "$$P(\\mu, k)=\\frac{\\mu^{k}}{k!}e^{-\\mu},$$ \n", "\n", "with the parameter $\\mu$ coinciding with the expectation value of $P(\\mu, k)$. Note that $k$ is a running parameter that we usually \"integrate over\", e.g., to obtain $\\beta$, equivalent to the variable $x$ of a continuous probability density function, like the normal distribution. \n", "\n", "The physics model that we are having in mind would be implemented in our expectation for $\\mu$ which, e.g., could be a sum of signal ($s$) and background ($b$) contributions ($\\mu=s+b$), or just $b$. In the course of the following exercises our goal will be to better understand the following points: \n", "\n", " * The meaning of $\\mu_{0.95}$ and how to construct it.\n", " * Identifying a few difficulties that you usually encounter when following such a procedure close to physical boundaries." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Exercise 2: Neyman construction**\n", "\n", "In a first step we do not care whether $\\mu=s+b$ or just $b$, and take it as given with some random value, which we freely choose to be $\\mu=8.5$. You can use the function `neyman_construction(columns='one', cl=0.95)`, provided by the module _material.py_, to display $P(\\mu=8.5, k)$. On the x-axis of this display you can read off the probabilities to observe $k=1,\\,2,\\,3,\\,\\ldots$ counts in the experiment for a given $\\mu$ (shown on the y-axis). \n", "\n", "a) What is the probability to observe $k=6$?\n", "\n", "b) What is the probability to observe $k\\leq6$? \n", " \n", "c) The red line shown in the display represents the boundary, at which $\\sum_k P(\\mu, k)\\leq\\beta$ holds, where each sum of this kind is assumed to start from $k=1$ (in this example, in general $k$=0) up to the value $k=k'$ for which the inequality is still true (assume $\\beta=0.05$). What is the exact probability to observe a value of $k$ below the red line for $\\mu=8.5$? Can you explain what the red line means when interpreting this setup as a limit on $\\mu$? \n", "\n", "In practice we have a statistical model of our measurement, here represented by $P(\\mu, k)$ and an outcome of the experiment $N_\\mathrm{obs}$ observations. To obtain an upper limit on $\\mu$ at 95% CL run the following procedure: \n", "\n", " * Scan $\\mu$, as fine granular as possible, from small to large values; \n", " * For each value of $\\mu$ calculate $\\sum_k P(\\mu, k)\\leq\\beta$;\n", " * The upper limit at 95% CL, $\\mu_{0.95}$, is the smallest value of $\\mu$ with $\\sum_k P(\\mu, k)\\leq\\beta$.\n", " \n", "This procedure is called the **Neyman construction**. You can use the function `neyman_construction(columns='all', cl=0.95)` to complete the display, you produced, before to turn it into a scatter plot for varying values of $\\mu$ (on the y-axis). In the scatter plot, the red line corresponds to the 95% CL that we have chosen for given observations $N_\\mathrm{obs}$ on the x-axis, as described above.\n", "\n", "d) The procedure above describes how to obtain an upper limit $\\mu_{0.95}$. What would you have to do to determine a lower limit?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "neyman_construction(columns='one', cl=??)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Exercise 3: Probability to observe less than $N_{\\mathrm{obs}}$ events for a given $\\mu$**\n", "\n", "Calculate the probability $\\mathcal{P}(k\\leq N_\\mathrm{obs} | P(\\mu=20, k))$, for a Poisson distributed random variable $k$ and $N_{\\mathrm{obs}}=10,\\,20,\\,25,\\,30$. Note: to determine the cumulative distribution function (CDF) of $P(\\mu, k)$, you can use the function `poisson_cdf(k=N_obs, mu)`, provided by the module _material.py_.\n", "\n", "| No. of counts | Probability |\n", "| --- | --- |\n", "| 10 | |\n", "| 20 | |\n", "| 25 | |\n", "| 30 | |" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "poisson_cdf(k=??, mu=??)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Exercise 4: Upper limit for a given $N_{\\mathrm{obs}}$ and varying confidence $1-\\beta$**\n", "\n", "Calculate the following upper limits for a simple counting experiment with the outcome $N_\\mathrm{obs}=20$: $\\mu_{0.90}$, $\\mu_{0.95}$, $\\mu_{0.99}$. For this purpose you can use the function `find_limit(obs=N_obs, cl=CL)`, provided by the module _material.py_, where $\\mathrm{CL}=1-\\beta$.\n", "\n", "| CL | $\\mu_{\\mathrm{x}}$ |\n", "| --- | --- |\n", "| 0.90 | |\n", "| 0.95 | |\n", "| 0.99 | |" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "find_limit(obs=??, cl=??)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Exercise 5: Upper limit for a signal $s$ for a given background expectation $b$**\n", "\n", "We now imagine a counting experiment with an outcome of $N_\\mathrm{obs}=2$. We want to calculate a 95% CL upper limit $s_{0.95}$ on the signal contribution $s$ for a _perfectly known_ background expectation $b$ and $\\mu=s+b$. Compute $s_{0.95}$ for $b=0,\\,2,\\,5,\\,10$.\n", "\n", "| b | Limit from **Exercise 5 a)** | Limit from **Exercise 5 b) & c)** |\n", "| :---: | :---: | :---: |\n", "| 0 | | |\n", "| 2 | | |\n", "| 5 | | |\n", "| 10 | | |\n", "\n", "a) Naive approach: compute the limit as $s_{0.95}=\\mu_{0.95}-b$. Are all results meaningful? Use the `find_limit(obs=N_obs, cl=CL)` function again. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "find_limit(obs=??, cl=??)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For a counting experiment 0 is a (natural) boundary. To prevent uninterpretable results near such boundaries, in particle physics we are using the _modified frequentist_ $\\mathrm{CL}_{\\mathrm{s}}$ confidence interval, calculated under the assumption of the \"signal plus background\" ($\\mathrm{CL}_{\\mathrm{s+b}}$) and \"background-only\" ($\\mathrm{CL}_{\\mathrm{b}}$) hypotheses (see figure below), according to:\n", "\n", "$$\\mathrm{CL}_{\\mathrm{s}}=\\frac{\\mathrm{CL}_\\mathrm{s+b}}{\\mathrm{CL}_\\mathrm{b}}=\\frac{\\mathcal{P}(k\\leq N_\\mathrm{obs}|{P(\\mu=\\mathrm{s+b}, k))}}{\\mathcal{P}(k\\leq N_\\mathrm{obs}|{P(\\mu=\\mathrm{b}, k))}}.$$\n", "\n", "In contrast to **Exercise 5 a)**, where $b$ has been assumed to be perfectly known, with this approach $b$ is introduced as the expectation value of another Poisson distribution, i.e. that $b$ is also subject to fluctuations with a variance of $\\sqrt{b}$. \n", "\n", "Effectively $\\mathcal{P}(k\\leq N_\\mathrm{obs}|{P(\\mu=\\mathrm{b}, k))}\\leq 1$, which means that $\\mathcal{P}(k\\leq N_\\mathrm{obs}|{P(\\mu=\\mathrm{s+b}, k))}$ will always be gauged to a larger value $\\mathrm{CL}_{\\mathrm{s}}\\geq \\mathrm{CL}_{\\mathrm{s+b}}$. At the same time $\\mathrm{CL}_{\\mathrm{s}}\\leq1$, since $\\mathcal{P}(k\\leq N_\\mathrm{obs}|{P(\\mu=\\mathrm{s+b}, k))}\\leq \\mathcal{P}(k\\leq N_\\mathrm{obs}|{P(\\mu=\\mathrm{b}, k))}$ by construction. A $\\mathrm{CL}_{\\mathrm{s}}$ limit is always more conservative than a limit based on $\\mathrm{CL}_{\\mathrm{s+b}}$. With this approach a limit on $s$ can never turn negative near boundaries. Also a tiny signal on a huge expected background can never be excluded.\n", "\n", "
\"Drawing\"
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "b) Use the function `poisson_cdf(k=N_obs, mu)` to calculate $\\mathrm{CL}_{\\mathrm{s}}$ and determine $s_{0.95}$ using the modified frequentist approach decribed above. For this purpose use the third row in table above of exercise 5a) and assume $N_\\mathrm{obs}=2$ and $b=5$. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "s = ??\n", "b = ??\n", "obs = ??\n", "\n", "poisson_cdf(k=??, mu=??)\n", "\n", "print('CLs: ', ??,'< 0.05')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "c) The procedure of **Exercise 5 b)** is automated by the function `find_cls_limits(obs=N_obs, b=background, cl=CL)`, provided by the module _material.py_. Use this function to determine $s_{0.95}$ for the observed events given in the table of **Exercise 5 a)**, using the modified frequentist $\\mathrm{CL}_{\\mathrm{s}}$ approach." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "find_cls_limit(obs=??, b=??, cl=??)" ] } ], "metadata": { "interpreter": { "hash": "002896ba1c76203336bd9e9fc74aac49db99ae768a34cad7088e5540b4dfcef2" }, "kernelspec": { "display_name": "Python 3", "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.8.2" } }, "nbformat": 4, "nbformat_minor": 2 }