{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# SMIPP 21/22 - Exercise Sheet 6\n", "\n", "## Prof. Dr. K. Reygers, Dr. R. Stamen, Dr. M. Völkl\n", "\n", "## Hand in by: Thursday, December 2nd: 12:00\n", "### Submit the file(s) through the Übungsgruppenverwaltung\n", "\n", "\n", "### Names (up to two):\n", "### Points: " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 6.1 The lighthouse problem (10 points)\n", "\n", "A lighthouse is somewhere off a piece of straight coastline at a position $\\alpha$ along the shore and a distance $\\beta$ out at sea (see figure). Its lamp rotates with constant velocity but it is not on at all the times: it emits a series of short highly collimated flashes at random intervals and hence at random angles $\\theta_k$. These pulses are intercepted on the coast by photo-detectors that record only the fact that a flash has occurred, but not the angle from which it came. $N = 1000$ flashes have so far been recorded with positions $x_k$.\n", "\"Drawing\"\n", "\n", "Our goal is to estimate the lighthouse position (i.e. the parameters $\\alpha$ and $\\beta$) using the maximum likelihood method.\n", "\n", "This problem is discussed, e.g., in Sivia, Skilling, \"Data Analysis: A Bayesian Tutorial\".\n", "\n", "a) Plot the distribution of the recorded positions as a histogram.\n", "\n", "b) Show that the normalzed probability density function for measuring the light flash at position $x$ for given values $\\alpha$ and $\\beta$ is given by\n", "\n", "$$ p(x | \\alpha, \\beta) = \\frac{1}{\\pi}\\frac{\\beta}{(x - \\alpha)^2 + \\beta^2} $$\n", "\n", "Hint: Start from $x_k(\\theta_k) = \\alpha + \\beta\\tan\\theta_k$. The angle $\\theta$ is uniformly distributed between $-\\pi/2$ and $\\pi/2$. With the formula for a change of variables one arrives at $p(x | \\alpha, \\beta)$.\n", "\n", "c) Define a python function that takes the values of $x_k$ (as a numpy array) and the parameters $\\alpha$ and $\\beta$ as input and returns the corresponding values of the probability density function.\n", "\n", "d) Write a python function that returns the negative log-likelihood for the given dataset as a function of $\\alpha$ and $\\beta$.\n", "\n", "e) Plot the negative log-logarithm of the likelihood as a function of $\\alpha$ (in the range [−10,40]) for the following values of $\\beta = 20, 30, 40, 50$ (4 curves on the same plot). As an optional task, you can also plot a 2d contour or a 3d visualization of $- \\mathrm{ln} L(\\alpha, \\beta)$, see e.g. the matplotlib [contour demo](https://matplotlib.org/3.1.1/gallery/images_contours_and_fields/contour_demo.html).\n", "\n", "f) Determine the maximum likelihood estimate for $\\alpha$ and $\\beta$ with the aid of iminuit. Compare the maximum likelihood estimate for $\\alpha$ with the mean and median of the $x_k$ as obtained by `['xk'].mean()` and`d['xk'].median()`, respectively." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Load the data as follows: " ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0 31.045519\n", "1 39.151935\n", "2 1.813172\n", "3 -10.046478\n", "4 -2.235928\n", " ... \n", "995 -11.792170\n", "996 508689.961890\n", "997 -15.009988\n", "998 121.086618\n", "999 106.612406\n", "Name: xk, Length: 1000, dtype: float64\n" ] } ], "source": [ "import pandas as pd\n", "d = pd.read_csv('https://www.physi.uni-heidelberg.de/~reygers/lectures/2021/smipp/lighthouse.csv')\n", "print(d['xk'])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Solution:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 6.2 Binned likelihood fit (15 points)\n", "\n", "In order to reduce the large data size of the measurement of many particles, histograms are often used to store data in particle physics. The given data contain a signal peak above a constant background. You will try to extract an estimate for the number of signal particles in the sample. The signal is known to be a normal distribution centered around $0.5$ with a width of $0.05$\n", "\n", "a) Read in the bin edges and bin contents from the file with the given code. Plot the data. How many numbers are needed i) to store each individual event in the histogram and ii) to store the bin counts and histogram edges.\n", "\n", "b) For a simple counting estimate, assume that the signal is fully in the mass region $(0.35,0.65)$, corresponding to $3\\sigma$ of the signal width. Estimate the background from the histogram entries outside this region from a histogram. Sum up the bins in the signal region and subtract the estimated background. What is the uncertainty assuming Poissonian statistical fluctuations? \n", "\n", "Hint: `np.sum(n_i[(binCenters > 0.35) & (binCenters < 0.65)]` gives you the number of counts in $(0.35,0.65)$.\n", "\n", "c) Define a function which returns the binned negative log-likelihood for a particular contribution of signal and background. Find the signal and background counts using iminuit. Compare to the simple estimate from b).\n", "\n", "Here is the code snippet to read the data:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "\n", "n_i = np.loadtxt(\"https://www.physi.uni-heidelberg.de/~reygers/lectures/2021/smipp/HistoEntries.txt\")\n", "binEdges = np.loadtxt(\"https://www.physi.uni-heidelberg.de/~reygers/lectures/2021/smipp/BinEdges.txt\")\n", "\n", "binCenters = (binEdges[:-1]+binEdges[1:])/2\n", "binwidth = binEdges[1]-binEdges[0]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Solution:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 6.3 Linear least squares and error propagation (15 points)\n", "Perform a least squares fit of a parabola\n", "\n", "$$ y(x) = \\theta_0 + \\theta_1 x + \\theta_2 x^2$$\n", "\n", "for the four independent measurments $(x_i, y_i)$ given by $(-0.6, 5 \\pm 2)$, $(-0.2, 3 \\pm 1)$, $(0.2, 5 \\pm 1)$, $(0.6, 8 \\pm 2)$. \n", "\n", "a) Determine the best fit parameters $\\hat \\theta_i$ and their covariances using the formula for linear least squares fits.\n", "\n", "b) Plot the fitted parabola and the $1\\sigma$ error band around it as obtaind from the error propagation formula $\\sigma_y^2 = A^\\mathrm{T} U A$. What is the predicted value $y$ at $x=1$ and its uncertainty?\n", "\n", "c) Compare the fit result and the covariance matrix with the results obtained with iminuit. Do the results agree?\n", "\n", "This problem is taken from [from W. Metzger, Statistical Methods in Data Analysis](www.hef.ru.nl/~wes/stat_course/statist.pdf).\n", "\n", "Hints: \n", "\n", "1) The following numpy functions might come in handy\n", "* diagonal matrix from vector v in numpy: *A = np.diagflat(v)*\n", "* matrix from column vectors v0, v1, v2: *A = np.column_stack((v0, v1, v2))*\n", "* multiplication of matrices A and B in numpy: *C = A.dot(B)*\n", "* transposed matrix: *A_T = np.transpose(A)*\n", "* inverse matrix: *A_inv = inv(A)*, this requires *from numpy.linalg import inv*\n", "\n", "2) It might be useful to write a function which returns the uncertainty $\\sigma_y$ for a numpy array of $x$ values. A function $f$ that only works for a scalar value $x$ can be vectorized (i.e., it works also for numpy arrays) with *numpy.vectorize*.\n", "\n", "3) Data points with error bars can be drawn with *plt.errorbar(x, y, yerr=sigma_y)*\n", "\n", "4) A band between $y$ values can be drawn with the aid of *matplotlib.pyplot.fill_between*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Solution" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "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.3" } }, "nbformat": 4, "nbformat_minor": 4 }