From 0897b029624315cbcc1ff81f328ae292537e2e24 Mon Sep 17 00:00:00 2001
From: Hartmut Stadie <hartmut.stadie@cern.ch>
Date: Mon, 7 Oct 2024 13:51:33 +0200
Subject: [PATCH] first estimator parts
---
lecture_3.ipynb | 1098 ++++++++++++++++-------------------------------
1 file changed, 359 insertions(+), 739 deletions(-)
diff --git a/lecture_3.ipynb b/lecture_3.ipynb
index 9fbb8cf..eda3ad9 100644
--- a/lecture_3.ipynb
+++ b/lecture_3.ipynb
@@ -2,74 +2,77 @@
"cells": [
{
"cell_type": "markdown",
- "id": "977b88f5-3cb7-445b-adac-f44be4d69c90",
- "metadata": {
- "slideshow": {
- "slide_type": "slide"
- },
- "tags": []
- },
+ "id": "612072eb",
+ "metadata": {},
"source": [
- "## Wieder einmal Bundesligatore..."
+ "# Lecture 3\n",
+ "\n",
+ "---\n",
+ "\n",
+ "## Parameter estimation\n",
+ "<br>\n",
+ "<br>\n",
+ "\n",
+ " Hartmut Stadie\n",
+ "\n",
+ "hartmut.stadie@uni-hamburg.de"
]
},
{
- "cell_type": "code",
- "execution_count": 6,
- "id": "641ec6cc-6f1a-4399-b652-7e3da333782a",
+ "cell_type": "markdown",
+ "id": "a038c48f-7143-480d-b0ba-817f0e5211c1",
"metadata": {
"slideshow": {
- "slide_type": "subslide"
+ "slide_type": "slide"
},
"tags": []
},
- "outputs": [
- {
- "name": "stdout",
- "output_type": "stream",
- "text": [
- "(array([24., 45., 66., 67., 58., 25., 18., 1., 2., 0.]), array([-0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5]), <BarContainer object of 10 artists>)\n"
- ]
- },
- {
- "data": {
- "image/png": "iVBORw0KGgoAAAANSUhEUgAAAigAAAG0CAYAAAAYQdwgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAgPElEQVR4nO3df0yd5f3/8dextaegh+Na7TkQaXt0p/7CaleUlVZhU0iQNBoWnVK1xi1p11aLZKtFlohGD5VPJCxjYmiW2saw+oc/2qWzgtNRDWnEKrNDUzXFljmPRMfOwRYPsVyfP/z2fHuk+2wHDpyrN89Hcic717k5581J3Hnm4j49LmOMEQAAgEXOSvcAAAAA30WgAAAA6xAoAADAOgQKAACwDoECAACsQ6AAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOskFSgLFy6Uy+Uac6xfv16SZIxRXV2dcnJylJGRoeLiYvX29k7K4AAAwLlmJnNyd3e3Tpw4Eb/9t7/9TSUlJbr11lslSQ0NDWpsbNQzzzyjRYsW6bHHHlNJSYkOHTokj8fzXz3H6Oio/vGPf8jj8cjlciUzHgAASBNjjIaGhpSTk6OzzkrBH2jMBGzcuNFcfPHFZnR01IyOjhq/32+2bNkSv//rr782Xq/XPP300//1Y/b39xtJHBwcHBwcHGfg0d/fP5G0iEtqB+VUIyMjevbZZ1VdXS2Xy6XDhw8rHA6rtLQ0fo7b7VZRUZG6urq0Zs2a0z5OLBZTLBaL3zb/78uV+/v7lZWVNd7xAADAFIpGo8rNzf2v/2Lyn4w7UF566SX961//0j333CNJCofDkiSfz5dwns/n05EjR/7t49TX1+uRRx4Zs56VlUWgAABwhknV5Rnj/iPR73//e5WVlSknJydh/buDGWP+z2FramoUiUTiR39//3hHAgAADjGuHZQjR47o1Vdf1QsvvBBf8/v9kr7dScnOzo6vDwwMjNlVOZXb7Zbb7R7PGAAAwKHGtYOybds2zZs3T+Xl5fG1QCAgv9+vjo6O+NrIyIg6OztVWFg48UkBAMC0kfQOyujoqLZt26bVq1dr5sz//+Mul0tVVVUKhUIKBoMKBoMKhULKzMxUZWVlSocGAADOlnSgvPrqqzp69KjuvffeMfdt2rRJw8PDWrdunQYHB1VQUKD29vaUXdELAACmB5c5+bleS0SjUXm9XkUiET7FAwDAGSLV7998Fw8AALAOgQIAAKxDoAAAAOsQKAAAwDoECgAAsA6BAgAArEOgAAAA6xAoAADAOgQKAACwzri+zRhwooWb96R7hKR9sqX8P58EAGcgdlAAAIB1CBQAAGAdAgUAAFiHQAEAANbhIllMijPxglMAgD3YQQEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYJ2kA+XTTz/VnXfeqblz5yozM1NXX321Dhw4EL/fGKO6ujrl5OQoIyNDxcXF6u3tTenQAADA2ZIKlMHBQS1fvlxnn322Xn75Zb3//vt68skndd5558XPaWhoUGNjo5qbm9Xd3S2/36+SkhINDQ2lenYAAOBQM5M5+YknnlBubq62bdsWX1u4cGH8fxtj1NTUpNraWlVUVEiStm/fLp/Pp7a2Nq1Zs2bMY8ZiMcVisfjtaDSa7O8AAAAcJqkdlN27dys/P1+33nqr5s2bpyVLlmjr1q3x+/v6+hQOh1VaWhpfc7vdKioqUldX12kfs76+Xl6vN37k5uaO81cBAABOkVSgHD58WC0tLQoGg3rllVe0du1a3X///dqxY4ckKRwOS5J8Pl/Cz/l8vvh931VTU6NIJBI/+vv7x/N7AAAAB0nqTzyjo6PKz89XKBSSJC1ZskS9vb1qaWnR3XffHT/P5XIl/JwxZszaSW63W263O9m5AQCAgyW1g5Kdna3LL788Ye2yyy7T0aNHJUl+v1+SxuyWDAwMjNlVAQAA+HeSCpTly5fr0KFDCWsffvihFixYIEkKBALy+/3q6OiI3z8yMqLOzk4VFhamYFwAADAdJPUnngceeECFhYUKhUK67bbb9NZbb6m1tVWtra2Svv3TTlVVlUKhkILBoILBoEKhkDIzM1VZWTkpvwAAAHCepALlmmuu0Ysvvqiamho9+uijCgQCampq0qpVq+LnbNq0ScPDw1q3bp0GBwdVUFCg9vZ2eTyelA8PAACcyWWMMeke4lTRaFRer1eRSERZWVnpHgfjtHDznnSPMC18sqU83SMAgKTUv3/zXTwAAMA6BAoAALAOgQIAAKyT1EWyAOxyJl7rw3UzAP4b7KAAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOsQKAAAwDoECgAAsA6BAgAArEOgAAAA6xAoAADAOgQKAACwDoECAACsQ6AAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOsQKAAAwDoECgAAsA6BAgAArEOgAAAA6xAoAADAOgQKAACwDoECAACsQ6AAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOsQKAAAwDoECgAAsA6BAgAArEOgAAAA6xAoAADAOgQKAACwDoECAACsQ6AAAADrJBUodXV1crlcCYff74/fb4xRXV2dcnJylJGRoeLiYvX29qZ8aAAA4GxJ76BcccUV+uyzz+LHwYMH4/c1NDSosbFRzc3N6u7ult/vV0lJiYaGhlI6NAAAcLakA2XmzJny+/3x44ILLpD07e5JU1OTamtrVVFRoby8PG3fvl3Hjx9XW1tbygcHAADOlXSgfPTRR8rJyVEgENDtt9+uw4cPS5L6+voUDodVWloaP9ftdquoqEhdXV3/9vFisZii0WjCAQAAprekAqWgoEA7duzQK6+8oq1btyocDquwsFBffvmlwuGwJMnn8yX8jM/ni993OvX19fJ6vfEjNzd3HL8GAABwkqQCpaysTD/5yU905ZVX6sYbb9SePXskSdu3b4+f43K5En7GGDNm7VQ1NTWKRCLxo7+/P5mRAACAA03oY8bnnHOOrrzySn300UfxT/N8d7dkYGBgzK7Kqdxut7KyshIOAAAwvU0oUGKxmD744ANlZ2crEAjI7/ero6Mjfv/IyIg6OztVWFg44UEBAMD0MTOZk3/5y19q5cqVmj9/vgYGBvTYY48pGo1q9erVcrlcqqqqUigUUjAYVDAYVCgUUmZmpiorKydrfgAA4EBJBcrf//533XHHHfriiy90wQUX6Ic//KH279+vBQsWSJI2bdqk4eFhrVu3ToODgyooKFB7e7s8Hs+kDA8AAJzJZYwx6R7iVNFoVF6vV5FIhOtRzmALN+9J9wiw1CdbytM9AoBJkOr3b76LBwAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1ZqZ7APxnCzfvSfcIAABMKXZQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYJ0JBUp9fb1cLpeqqqria8YY1dXVKScnRxkZGSouLlZvb+9E5wQAANPIuAOlu7tbra2tWrx4ccJ6Q0ODGhsb1dzcrO7ubvn9fpWUlGhoaGjCwwIAgOlhXIHy1VdfadWqVdq6dau+973vxdeNMWpqalJtba0qKiqUl5en7du36/jx42pra0vZ0AAAwNnGFSjr169XeXm5brzxxoT1vr4+hcNhlZaWxtfcbreKiorU1dV12seKxWKKRqMJBwAAmN5mJvsDO3fu1DvvvKPu7u4x94XDYUmSz+dLWPf5fDpy5MhpH6++vl6PPPJIsmMAAAAHS2oHpb+/Xxs3btSzzz6r2bNn/9vzXC5Xwm1jzJi1k2pqahSJROJHf39/MiMBAAAHSmoH5cCBAxoYGNDSpUvjaydOnNC+ffvU3NysQ4cOSfp2JyU7Ozt+zsDAwJhdlZPcbrfcbvd4ZgcAAA6V1A7KDTfcoIMHD6qnpyd+5Ofna9WqVerp6dFFF10kv9+vjo6O+M+MjIyos7NThYWFKR8eAAA4U1I7KB6PR3l5eQlr55xzjubOnRtfr6qqUigUUjAYVDAYVCgUUmZmpiorK1M3NQAAcLSkL5L9TzZt2qTh4WGtW7dOg4ODKigoUHt7uzweT6qfCgAAOJTLGGPSPcSpotGovF6vIpGIsrKy0j2OFRZu3pPuEYCU+WRLebpHADAJUv3+zXfxAAAA6xAoAADAOgQKAACwDoECAACsQ6AAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOsQKAAAwDoECgAAsA6BAgAArEOgAAAA6xAoAADAOgQKAACwDoECAACsQ6AAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOsQKAAAwDoECgAAsA6BAgAArEOgAAAA6xAoAADAOgQKAACwDoECAACsQ6AAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOsQKAAAwDoECgAAsA6BAgAArEOgAAAA6xAoAADAOgQKAACwDoECAACsQ6AAAADrJBUoLS0tWrx4sbKyspSVlaVly5bp5Zdfjt9vjFFdXZ1ycnKUkZGh4uJi9fb2pnxoAADgbEkFyoUXXqgtW7bo7bff1ttvv60f//jHuvnmm+MR0tDQoMbGRjU3N6u7u1t+v18lJSUaGhqalOEBAIAzJRUoK1eu1E033aRFixZp0aJFevzxx3Xuuedq//79MsaoqalJtbW1qqioUF5enrZv367jx4+rra1tsuYHAAAONO5rUE6cOKGdO3fq2LFjWrZsmfr6+hQOh1VaWho/x+12q6ioSF1dXf/2cWKxmKLRaMIBAACmt6QD5eDBgzr33HPldru1du1avfjii7r88ssVDoclST6fL+F8n88Xv+906uvr5fV640dubm6yIwEAAIdJOlAuueQS9fT0aP/+/frFL36h1atX6/3334/f73K5Es43xoxZO1VNTY0ikUj86O/vT3YkAADgMDOT/YFZs2bp+9//viQpPz9f3d3d+s1vfqMHH3xQkhQOh5WdnR0/f2BgYMyuyqncbrfcbneyYwAAAAeb8L+DYoxRLBZTIBCQ3+9XR0dH/L6RkRF1dnaqsLBwok8DAACmkaR2UB566CGVlZUpNzdXQ0ND2rlzp/7yl79o7969crlcqqqqUigUUjAYVDAYVCgUUmZmpiorKydrfgAA4EBJBcrnn3+uu+66S5999pm8Xq8WL16svXv3qqSkRJK0adMmDQ8Pa926dRocHFRBQYHa29vl8XgmZXgAAOBMLmOMSfcQp4pGo/J6vYpEIsrKykr3OFZYuHlPukcAUuaTLeXpHgHAJEj1+zffxQMAAKxDoAAAAOsQKAAAwDoECgAAsA6BAgAArEOgAAAA6xAoAADAOgQKAACwDoECAACsQ6AAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOsQKAAAwDoECgAAsA6BAgAArEOgAAAA6xAoAADAOjPTPQCA6WXh5j3pHmFcPtlSnu4RgGmFHRQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYJ1p9ymeM/UTBAAATCfsoAAAAOsQKAAAwDoECgAAsA6BAgAArEOgAAAA6xAoAADAOgQKAACwDoECAACsQ6AAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOsQKAAAwDoECgAAsA6BAgAArJNUoNTX1+uaa66Rx+PRvHnzdMstt+jQoUMJ5xhjVFdXp5ycHGVkZKi4uFi9vb0pHRoAADhbUoHS2dmp9evXa//+/ero6NA333yj0tJSHTt2LH5OQ0ODGhsb1dzcrO7ubvn9fpWUlGhoaCjlwwMAAGeamczJe/fuTbi9bds2zZs3TwcOHND1118vY4yamppUW1uriooKSdL27dvl8/nU1tamNWvWpG5yAADgWBO6BiUSiUiS5syZI0nq6+tTOBxWaWlp/By3262ioiJ1dXWd9jFisZii0WjCAQAAprdxB4oxRtXV1VqxYoXy8vIkSeFwWJLk8/kSzvX5fPH7vqu+vl5erzd+5ObmjnckAADgEOMOlA0bNui9997TH/7whzH3uVyuhNvGmDFrJ9XU1CgSicSP/v7+8Y4EAAAcIqlrUE667777tHv3bu3bt08XXnhhfN3v90v6diclOzs7vj4wMDBmV+Ukt9stt9s9njEAAIBDJbWDYozRhg0b9MILL+i1115TIBBIuD8QCMjv96ujoyO+NjIyos7OThUWFqZmYgAA4HhJ7aCsX79ebW1t2rVrlzweT/y6Eq/Xq4yMDLlcLlVVVSkUCikYDCoYDCoUCikzM1OVlZWT8gsAAADnSSpQWlpaJEnFxcUJ69u2bdM999wjSdq0aZOGh4e1bt06DQ4OqqCgQO3t7fJ4PCkZGAAAOF9SgWKM+Y/nuFwu1dXVqa6ubrwzAQCAaY7v4gEAANYhUAAAgHUIFAAAYJ1x/TsoADDdLNy8J90jJO2TLeXpHgEYN3ZQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdZIOlH379mnlypXKycmRy+XSSy+9lHC/MUZ1dXXKyclRRkaGiouL1dvbm6p5AQDANJB0oBw7dkxXXXWVmpubT3t/Q0ODGhsb1dzcrO7ubvn9fpWUlGhoaGjCwwIAgOlhZrI/UFZWprKystPeZ4xRU1OTamtrVVFRIUnavn27fD6f2tratGbNmolNCwAApoWUXoPS19encDis0tLS+Jrb7VZRUZG6urpO+zOxWEzRaDThAAAA01tKAyUcDkuSfD5fwrrP54vf91319fXyer3xIzc3N5UjAQCAM9CkfIrH5XIl3DbGjFk7qaamRpFIJH709/dPxkgAAOAMkvQ1KP8Xv98v6dudlOzs7Pj6wMDAmF2Vk9xut9xudyrHAAAAZ7iU7qAEAgH5/X51dHTE10ZGRtTZ2anCwsJUPhUAAHCwpHdQvvrqK3388cfx2319ferp6dGcOXM0f/58VVVVKRQKKRgMKhgMKhQKKTMzU5WVlSkdHAAAOFfSgfL222/rRz/6Ufx2dXW1JGn16tV65plntGnTJg0PD2vdunUaHBxUQUGB2tvb5fF4Ujc1AABwNJcxxqR7iFNFo1F5vV5FIhFlZWWl/PEXbt6T8scEABt9sqU83SNgGkn1+zffxQMAAKxDoAAAAOsQKAAAwDoECgAAsA6BAgAArEOgAAAA6xAoAADAOgQKAACwDoECAACsQ6AAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOsQKAAAwDoECgAAsA6BAgAArEOgAAAA6xAoAADAOgQKAACwDoECAACsQ6AAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOsQKAAAwDoECgAAsA6BAgAArEOgAAAA6xAoAADAOgQKAACwDoECAACsQ6AAAADrECgAAMA6M9M9AAAAJy3cvCfdIyTtky3l6R7BkdhBAQAA1iFQAACAdQgUAABgHQIFAABYh4tkAcChzsQLToGT2EEBAADWmbRAeeqppxQIBDR79mwtXbpUb7zxxmQ9FQAAcJhJCZTnnntOVVVVqq2t1bvvvqvrrrtOZWVlOnr06GQ8HQAAcJhJCZTGxkb97Gc/089//nNddtllampqUm5urlpaWibj6QAAgMOk/CLZkZERHThwQJs3b05YLy0tVVdX15jzY7GYYrFY/HYkEpEkRaPRVI8mSRqNHZ+UxwUATE+T9X51pjn5OhhjUvJ4KQ+UL774QidOnJDP50tY9/l8CofDY86vr6/XI488MmY9Nzc31aMBAJBy3qZ0T2CXL7/8Ul6vd8KPM2kfM3a5XAm3jTFj1iSppqZG1dXV8dujo6P65z//qblz5572fFtFo1Hl5uaqv79fWVlZ6R7HsXidpwav89TgdZ4avM5TIxKJaP78+ZozZ05KHi/lgXL++edrxowZY3ZLBgYGxuyqSJLb7Zbb7U5YO++881I91pTJysriP4ApwOs8NXidpwav89TgdZ4aZ52VmstbU36R7KxZs7R06VJ1dHQkrHd0dKiwsDDVTwcAABxoUv7EU11drbvuukv5+flatmyZWltbdfToUa1du3Yyng4AADjMpATKT3/6U3355Zd69NFH9dlnnykvL09/+tOftGDBgsl4Oiu43W49/PDDY/5chdTidZ4avM5Tg9d5avA6T41Uv84uk6rPAwEAAKQI38UDAACsQ6AAAADrECgAAMA6BAoAALAOgZICTz31lAKBgGbPnq2lS5fqjTfeSPdIjlJfX69rrrlGHo9H8+bN0y233KJDhw6leyzHq6+vl8vlUlVVVbpHcaRPP/1Ud955p+bOnavMzExdffXVOnDgQLrHcpRvvvlGv/71rxUIBJSRkaGLLrpIjz76qEZHR9M92hlt3759WrlypXJycuRyufTSSy8l3G+MUV1dnXJycpSRkaHi4mL19vYm/TwEygQ999xzqqqqUm1trd59911dd911Kisr09GjR9M9mmN0dnZq/fr12r9/vzo6OvTNN9+otLRUx44dS/dojtXd3a3W1lYtXrw43aM40uDgoJYvX66zzz5bL7/8st5//309+eSTZ/S/om2jJ554Qk8//bSam5v1wQcfqKGhQf/zP/+j3/72t+ke7Yx27NgxXXXVVWpubj7t/Q0NDWpsbFRzc7O6u7vl9/tVUlKioaGh5J7IYEKuvfZas3bt2oS1Sy+91GzevDlNEznfwMCAkWQ6OzvTPYojDQ0NmWAwaDo6OkxRUZHZuHFjukdynAcffNCsWLEi3WM4Xnl5ubn33nsT1ioqKsydd96ZpomcR5J58cUX47dHR0eN3+83W7Zsia99/fXXxuv1mqeffjqpx2YHZQJGRkZ04MABlZaWJqyXlpaqq6srTVM5XyQSkaSUfSEVEq1fv17l5eW68cYb0z2KY+3evVv5+fm69dZbNW/ePC1ZskRbt25N91iOs2LFCv35z3/Whx9+KEn661//qjfffFM33XRTmidzrr6+PoXD4YT3RbfbraKioqTfFyft24yngy+++EInTpwY8yWIPp9vzJclIjWMMaqurtaKFSuUl5eX7nEcZ+fOnXrnnXfU3d2d7lEc7fDhw2ppaVF1dbUeeughvfXWW7r//vvldrt19913p3s8x3jwwQcViUR06aWXasaMGTpx4oQef/xx3XHHHekezbFOvved7n3xyJEjST0WgZICLpcr4bYxZswaUmPDhg1677339Oabb6Z7FMfp7+/Xxo0b1d7ertmzZ6d7HEcbHR1Vfn6+QqGQJGnJkiXq7e1VS0sLgZJCzz33nJ599lm1tbXpiiuuUE9Pj6qqqpSTk6PVq1enezxHS8X7IoEyAeeff75mzJgxZrdkYGBgTD1i4u677z7t3r1b+/bt04UXXpjucRznwIEDGhgY0NKlS+NrJ06c0L59+9Tc3KxYLKYZM2akcULnyM7O1uWXX56wdtlll+n5559P00TO9Ktf/UqbN2/W7bffLkm68sordeTIEdXX1xMok8Tv90v6diclOzs7vj6e90WuQZmAWbNmaenSpero6EhY7+joUGFhYZqmch5jjDZs2KAXXnhBr732mgKBQLpHcqQbbrhBBw8eVE9PT/zIz8/XqlWr1NPTQ5yk0PLly8d8VP7DDz909BeqpsPx48d11lmJb3MzZszgY8aTKBAIyO/3J7wvjoyMqLOzM+n3RXZQJqi6ulp33XWX8vPztWzZMrW2turo0aNau3ZtukdzjPXr16utrU27du2Sx+OJ71h5vV5lZGSkeTrn8Hg8Y67rOeecczR37lyu90mxBx54QIWFhQqFQrrtttv01ltvqbW1Va2trekezVFWrlypxx9/XPPnz9cVV1yhd999V42Njbr33nvTPdoZ7auvvtLHH38cv93X16eenh7NmTNH8+fPV1VVlUKhkILBoILBoEKhkDIzM1VZWZncE6XiY0bT3e9+9zuzYMECM2vWLPODH/yAj7+mmKTTHtu2bUv3aI7Hx4wnzx//+EeTl5dn3G63ufTSS01ra2u6R3KcaDRqNm7caObPn29mz55tLrroIlNbW2tisVi6Rzujvf7666f9/+TVq1cbY779qPHDDz9s/H6/cbvd5vrrrzcHDx5M+nlcxhiTiqICAABIFa5BAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEwJYqLi1VVVZXuMQCcIQgUAABgHQIFAABYh0ABkBZ79+6V1+vVjh070j0KAAsRKACm3M6dO3Xbbbdpx44duvvuu9M9DgALESgAptRTTz2ltWvXateuXbr55pvTPQ4AS81M9wAApo/nn39en3/+ud58801de+216R4HgMXYQQEwZa6++mpdcMEF2rZtm4wx6R4HgMUIFABT5uKLL9brr7+uXbt26b777kv3OAAsxp94AEypRYsW6fXXX1dxcbFmzpyppqamdI8EwEIECoApd8kll+i1115TcXGxZsyYoSeffDLdIwGwjMvwh2AAAGAZrkEBAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgnf8FQX9AbUh+Cg8AAAAASUVORK5CYII=",
- "text/plain": [
- "<Figure size 640x480 with 1 Axes>"
- ]
- },
- "metadata": {},
- "output_type": "display_data"
- }
- ],
"source": [
- "import numpy as np\n",
- "import matplotlib.pyplot as plt\n",
+ "## Estimator\n",
+ "\n",
+ "Estimator $\\hat a$ for parameter $a$ from sample $x_1,\\dots x_n$\n",
+ "\n",
+ "Criteria for good estimators:\n",
"\n",
- "summe = np.loadtxt(\"../summe.txt\")\n",
+ "- consistent: \n",
+ " $\\lim_{n\\to \\infty } \\hat a = a$\n",
"\n",
+ "- unbiased: \n",
+ " $E[\\hat a]= a$\n",
"\n",
- "freq = plt.hist(summe, bins=10, range=(-0.5,9.5))\n",
- "plt.xlabel(\"k\")\n",
- "print(freq)"
+ "- efficient: $V[\\hat a]$ is small\n",
+ "\n",
+ "Exercise: Estimate the mean $\\mu$ and variance $\\sigma^2$ of a p.d.f. $p(x)$ from a sample\n",
+ "$x_1,\\dots,x_n$"
]
},
{
"cell_type": "markdown",
- "id": "e6bcf92b-386d-492b-9f09-e29118e93505",
+ "id": "d32672c8-da0b-468b-bf73-df40aef7245a",
"metadata": {
"slideshow": {
- "slide_type": "subslide"
+ "slide_type": "slide"
},
"tags": []
},
"source": [
- "18-mal fielen sechs Tore. Wie groß ist der Fehler auf die 18? Wie groß ist das 68,27\\%-Konfidenzintervall für das 6-Tore-Bin?"
+ "### Estimate for the mean\n",
+ "\n",
+ "Estimator for the mean $\\mu$:\n",
+ "$$\\hat \\mu = \\bar x = \\frac{1}{n}\\sum_{i=1}^n x_i$$\n",
+ "\n",
+ "Tests:\n",
+ "\n",
+ "- consistent:\n",
+ " $$\\lim_{n\\to \\infty } \\hat \\mu = \\lim_{n\\to \\infty } \\frac{1}{n}\\sum_{i=1}^n x_i = \\int x p(x)\\,dx = \\mu$$\n",
+ "\n",
+ "- unbiased:\n",
+ " $$E[\\hat \\mu]= E[ \\frac{1}{n}\\sum_{i=1}^n x_i ] = \\frac{1}{n} \\sum_{i=1}^n E[x_i] = \\frac{1}{n} \\sum_{i=1}^n \\mu = \\mu$$"
]
},
{
"cell_type": "markdown",
- "id": "8c4a1016-7fc9-4eb6-b515-b51c0d16fed8",
+ "id": "4021d6cd-150e-4a6d-bdc4-ed049fd983d3",
"metadata": {
"slideshow": {
"slide_type": "subslide"
@@ -77,78 +80,83 @@
"tags": []
},
"source": [
- "naiv: $k=18$, Poisson mit $\\hat \\mu = 18$ und $\\sigma = \\sqrt{\\hat \\mu} = \\sqrt{18}$: $18\\pm 4.12$"
+ "### Estimator for the mean \n",
+ "\n",
+ "Estimator for the mean: $\\mu$:\n",
+ "$$\\hat \\mu = \\bar x = \\frac{1}{n}\\sum_{i=1}^n x_i$$\n",
+ "\n",
+ "Test:\n",
+ "\n",
+ "- efficient: $V[\\hat a]$ is small\n",
+ " $$V[\\hat \\mu] = V[ \\frac{1}{n}\\sum_{i=1}^n x_i ] = \\frac{1}{n^2} V[ \\sum_{i=1}^n x_i ] = \\frac{1}{n^2} \\sum_{i=1}^n\\sum_{j=1}^n \\text{cov}(x_i, x_j)$$\n",
+ " $$V[\\hat \\mu] = \\frac{1}{n^2} \\sum_{i=1}^n \\text{cov}(x_i, x_i)= \\frac{n\\sigma^2}{n^2} = \\frac{\\sigma^2}{n}$$\n",
+ " (or via the) central limit theorem)"
]
},
{
"cell_type": "markdown",
- "id": "b26cd06e-7700-4011-b5d8-6bf56e22489f",
+ "id": "cc6b6159",
"metadata": {
"slideshow": {
"slide_type": "slide"
- },
- "tags": []
- },
- "source": [
- "## Konfidenz: \n",
- "$P(x_- \\le x \\le x_+) = \\int_{x_-}^{x^+} P(x) dx$"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "1903f299-4d16-41b0-bb48-5f17a07a9006",
- "metadata": {
- "slideshow": {
- "slide_type": "subslide"
- },
- "tags": []
+ }
},
"source": [
- "## Aber was ist hier $x$? $\\mu$ oder $k$???"
+ "### Test via Monte Carlo\n",
+ "\n",
+ "check for bias and efficiency for two different estimators of the mean:\n",
+ "$$\\hat \\mu_1 = \\frac{1}{n}\\sum_{i=1}^n x_i\\text{ and } \\hat \\mu_2 = \\frac{x_\\text{max}+ x_\\text{min}}{2}$$\n",
+ "for samples with 10 elements drawn from a Gaussian p.d.f. with $\\mu=1$ and $\\sigma = 1$ or from a uniform p.d.f. with the range [0, 2]:\n",
+ " - draw 1000 samples\n",
+ " - compute the means\n",
+ " - plot $\\hat \\mu - \\mu_\\text{true}$ and compute the mean and variance of this variable\n"
]
},
{
- "cell_type": "markdown",
- "id": "b54b11f9-e897-4d2a-a5df-961c16cdbb3f",
+ "cell_type": "code",
+ "execution_count": 47,
+ "id": "3f4ad5a3",
"metadata": {
- "slideshow": {
- "slide_type": "slide"
- },
- "tags": []
+ "cell_style": "center"
},
+ "outputs": [],
"source": [
- "Konfidenz für $P(k=18) = 1$ (Messung!) \n",
- "\n",
- "Suchen Konfidenzintervall in $\\mu$.\n",
- "\n",
+ "import scipy\n",
+ "import numpy as np\n",
+ "import matplotlib.pyplot as plt\n",
"\n",
- "Also: $P(\\mu_- \\le \\mu \\le\\mu_+) = \\frac{\\int_{\\mu_-}^{\\mu_+} P(k, \\mu) d\\mu}{\\int_{0}^{\\infty} P(k, \\mu) d\\mu}$"
+ "xs = scipy.stats.norm.rvs(1,1,size=(1000,10))\n",
+ "ys = scipy.stats.uniform.rvs(0,2,size=(1000,10))"
]
},
{
"cell_type": "code",
- "execution_count": 7,
- "id": "5052bdb2-6e4b-4ea6-b08e-ba22f2a4a3d2",
+ "execution_count": 49,
+ "id": "89de9c74",
"metadata": {
+ "cell_style": "center",
"slideshow": {
- "slide_type": "subslide"
- },
- "tags": []
+ "slide_type": "skip"
+ }
},
"outputs": [
{
"data": {
"text/plain": [
- "(1.000000000000001, 4.6105475479254554e-11)"
+ "(array([ 5., 13., 36., 114., 301., 327., 128., 55., 15., 6.]),\n",
+ " array([-0.46779215, -0.37579571, -0.28379928, -0.19180284, -0.09980641,\n",
+ " -0.00780997, 0.08418646, 0.1761829 , 0.26817933, 0.36017577,\n",
+ " 0.4521722 ]),\n",
+ " <BarContainer object of 10 artists>)"
]
},
- "execution_count": 7,
+ "execution_count": 49,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
- "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiwAAAGxCAYAAABBZ+3pAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABE8ElEQVR4nO3deXyV5Z338e9ZkpM9IQtZIISA7KsEl6C4tBqKHa2tnUHbauep2KG1o8A4rbiMrY7iWMeHOgqMCjN1pqNMXVqfFltiK4gSq0BYRGQNSYCEkEAWsp2cc+7nj5NzICSB7PdZPu/X67wkd+5z8ju3JPlyXb/rui2GYRgCAAAIYFazCwAAALgYAgsAAAh4BBYAABDwCCwAACDgEVgAAEDAI7AAAICAR2ABAAABz252AQPF4/Ho+PHjio+Pl8ViMbscAADQA4ZhqKGhQVlZWbJaux9HCZnAcvz4cWVnZ5tdBgAA6IPy8nKNHDmy28+HTGCJj4+X5H3DCQkJJlcDAAB6or6+XtnZ2f7f490JmcDimwZKSEggsAAAEGQu1s5B0y0AAAh4BBYAABDwCCwAACDgEVgAAEDAI7AAAICAR2ABAAABj8ACAAACHoEFAAAEPAILAAAIeAQWAAAQ8AgsAAAg4BFYAABAwCOwAEGg1eXWifoWtbS5zS4FAEwRMndrBkLRzvJa/dufD+qDAyfldHlksUh5o4Zp4dwxmjcl/aJ3NwWAUEFgAQKQYRh6/k8H9Ys/7ZfH8B6zWCTDkLaWntbW0m36q+mZeuab0xUTybcxgNDHTzogwBiGoSd+t1drPyqRJH1tZpZ+eN0lGp8ep4q6Fv3PX8q0etMh/W5XhcpONel/7rlScQ6+lQGENnpYgACz5sMSf1h54mtT9IvbL9WEjHhZLBZlJUXrgXkT9Pr3r1RybKR2Ha3T91/dKqfLY3LVADC4CCxAAPn0yCktf/cLSdJjN0/Wnfmjuzxv9uhk/ef/uUyxkTZtOVSjZzfsG8IqAWDoEViAANHS5taP39glt8fQrTOz9LdzRl/w/Okjk/TcgpmSpJc+OKxN+08OfpEAYBICCxAgVm08pJLqRg2Pd+iJW6f2aAXQvCkZuis/R5L08Nu71exk2TOA0ERgAQJARV2zVm06JEl67OYpio+K6PFzf/KVicpKjNLR08164f0Dg1UiAJiKwAIEgBf+fFBOl0eXj07WTdMyevXcWIdd/3TzFEnSy5tLdLy2eTBKBABTEVgAk5WfatK6T8slSf9QML5Pm8HNm5KuK3KT5XR59Iv3GGUBEHoILIDJ/uOjI3J5DF19SaquGJPSp9ewWCz68VcmSpJ+va1ch0+eGcgSAcB0BBbARA0tbfrfrd7RlYVzc/v1Wnk5w/SlicPlMbxTQwAQSggsgIne2HZUZ1pdGpsWq2vGpfX79RZdO1aS9Ob2ozrZ0Nrv1wOAQEFgAUzi8Rj65ZYjkqS/vSpXVmv/b2R42ehhmpGdJKfLo1eLjvT79QAgUBBYAJN8XFKjIzVNio+y67ZZIwbkNS0Wi/7umjGSpP/6uJR9WQCEDAILYJK3th+TJN08I2tA77g8b0qGRg6LVm1Tm9bvrhiw1wUAMxFYABM0OV16tz1MfOPSgRld8bFZLbrj8lGSpNc+KRvQ1wYAsxBYABNs2HNCjU63RiXHKC9n2IC//l/njZTNatHW0tPaf6JhwF8fAIYagQUwwZvbj0qSvjFrRJ82iruY4QlRumHScEnS//yFURYAwY/AAgyxU41OfXSwWpL09QGeDjqXb1ro7eJjcro8g/Z1AGAoEFiAIfbe3hPyGNLkzATlpMQO2teZOy5Nw+Mdqmtu08Z9VYP2dQBgKBBYgCH2x88qJUlfmdq7mxz2ls1q0S0zsiRJv91xfFC/FgAMNgILMITOtLq0+YB3OmjelMENLJJ0a/uUU+HeE6pvaRv0rwcAg4XAAgyhjfuq5HR7lJsaq/HpcYP+9aZkJeiS4XFyujz6Q/vIDgAEIwILMIT+uOeEJO/oymCsDjqfxWLRrTO900L/byfTQgCCF4EFGCIut0cf7D8pSbpx8vAh+7pfne4NLEWHalTb5ByyrwsAA4nAAgyRnUfrVNfcpsToCM0YmTRkXzc3NVYTM+Ll8hh6by+rhQAEJwILMEQ2tS8tvnpcquy2of3W8zX4/uEz7i0EIDgRWIAhsql9Oui68WlD/rXnT/MGlg8OVOtMq2vIvz4A9BeBBRgCNWdatetYnSTpWhMCy4T0eI1OiZHT5WETOQBBicACDIHNB6plGNKkzAQNT4ga8q9vsVj0lamZkqR3Wd4MIAgRWIAh4NsszozRFR/fzrobv6hSq8ttWh0A0BcEFmAIfHy4RpI0Z2yKaTVMH5GotHiHGp1ufVJyyrQ6AKAvCCzAICs/1aRjtc2yWy3KyxlmWh1Wq0XXT/CO8Pz5C/pYAAQXAgswyIoOeUdXZmQnKdZhN7WWL030blj3PoEFQJAhsACDzDcdlD/GvOkgn6vHpSnCZtGRmiYdPnnG7HIAoMcILMAgMgxDRe2B5coACCxxDrsuz02WxLQQgOBCYAEGUWlNkyrqWhRhM7d/5VzXT2ifFmI/FgBBhMACDCLfdNDM7CRFR9pMrsbL18fySckpNbS0mVwNAPQMgQUYREUB1L/iMyYtTqNTYtTmNvwNwQAQ6AgswCD6tH2/kysCKLBI0txx3uXNHx6sNrkSAOgZAgswSI7XNut4XYtsVotmZieZXU4HV49LlSR9eIDAAiA4EFiAQbK97LQkaVJmvOn7r5wvf2yKbFaLDlc36lhts9nlAMBFEViAQbKt1BtY8kYFxuqgcyVERWjGyERJ0ocHTppcDQBcHIEFGCTb2wPLrABZzny+q9v7WDYzLQQgCBBYgEHQ7HRrz/F6SQqY/VfOd017H8tHB6vl8RgmVwMAF9anwLJy5Url5uYqKipKeXl52rx58wXP37Rpk/Ly8hQVFaUxY8Zo9erVnc5ZsWKFJkyYoOjoaGVnZ2vJkiVqaWnpS3mA6XYdrZXLYyg9waERSdFml9OlGdlJinPYdbqpzR+uACBQ9TqwrFu3TosXL9bDDz+s4uJizZ07V/Pnz1dZWVmX55eUlOimm27S3LlzVVxcrIceekj33Xef3nzzTf85v/rVr/Tggw/qscce0969e7VmzRqtW7dOy5Yt6/s7A0y0rb3hNi9nmCwWi8nVdC3CZvXfLmDzQfpYAAS2XgeW5557TnfffbcWLlyoSZMmacWKFcrOztaqVau6PH/16tUaNWqUVqxYoUmTJmnhwoX63ve+p2effdZ/TlFRka666ip961vf0ujRo1VQUKA77rhDW7du7fs7A0zk61/Jy0k2uZILm8vyZgBBoleBxel0atu2bSooKOhwvKCgQFu2bOnyOUVFRZ3OnzdvnrZu3aq2Nu+24FdffbW2bdumTz75RJJ0+PBhrV+/Xl/96le7raW1tVX19fUdHkAgMAzj7AqhAO1f8fHtx7K19LRa2twmVwMA3etVYKmurpbb7VZ6enqH4+np6aqsrOzyOZWVlV2e73K5VF3t/Vfd7bffrieeeEJXX321IiIiNHbsWF1//fV68MEHu61l+fLlSkxM9D+ys7N781aAQVNa06TTTW2KtFs1OTPB7HIuaExqrIbHO+R0ebSjvNbscgCgW31quj1/Tt4wjAvO03d1/rnHN27cqCeffFIrV67U9u3b9dZbb+l3v/udnnjiiW5fc9myZaqrq/M/ysvL+/JWgAG361idJGlyZoIi7YG9EM9isfj7WHw3agSAQNSr7TdTU1Nls9k6jaZUVVV1GkXxycjI6PJ8u92ulBTvD8pHH31Ud955pxYuXChJmjZtmhobG/X9739fDz/8sKzWzj/0HQ6HHA5Hb8oHhsSu9pGK6e0bswW6K8ek6J2dx1V0qEaLbzC7GgDoWq/++RcZGam8vDwVFhZ2OF5YWKg5c+Z0+Zz8/PxO52/YsEGzZ89WRESEJKmpqalTKLHZbDIMwz8aAwSLXUe9IyzTRyaZW0gPXTnG2xhcXF5LHwuAgNXr8eqlS5fqlVde0dq1a7V3714tWbJEZWVlWrRokSTvVM1dd93lP3/RokUqLS3V0qVLtXfvXq1du1Zr1qzRAw884D/n5ptv1qpVq/T666+rpKREhYWFevTRR3XLLbfIZrMNwNsEhobbY+iz497AMiNIRlhyz+ljKS6rNbscAOhSr+/ItmDBAtXU1Ojxxx9XRUWFpk6dqvXr1ysnJ0eSVFFR0WFPltzcXK1fv15LlizRiy++qKysLD3//PO67bbb/Oc88sgjslgseuSRR3Ts2DGlpaXp5ptv1pNPPjkAbxEYOodOnlGT062YSJvGpMWZXU6P+PpY3tl5XB8frlH+2BSzSwKATixGiMy51NfXKzExUXV1dUpICOyVGQhdv95arn98Y5cuz03W//5dvtnl9Nj//KVMD729W1fkJmtdENUNIPj19Pd3YC9hAILM7mPBNR3kQx8LgEBHYAEG0M72httpQdJw60MfC4BAR2ABBojT5dHe9psIBtsIC/uxAAh0BBZggOw/0SCn26PE6AiNSo4xu5xeI7AACGQEFmCA7DxaK8m7YVyg3qH5Qvx9LGX0sQAIPAQWYIDsKvdtGBdc00E+uamxSot3yOn2+De/A4BAQWABBohvw7hpI5LMLaSPLBaLLhvtvbv0p0dOmVwNAHREYAEGgNPl0YETZyRJU7KCdx+g2TneaaGtBBYAAYbAAgyAg1Vn5HR7FB9l18hh0WaX02eXjW4PLKWn5fGExJ6SAEIEgQUYAJ9XeJczT85MCMqGW59JmfGKibSpocWl/VUNZpcDAH4EFmAAfN6+/8rkIJ4OkiS7zapZo9r7WEqYFgIQOAgswAD4vMLbcDs5M7gDiyTN9jfenja5EgA4i8AC9JNhGCEzwiKd08dC4y2AAEJgAfrpWG2z6ltcirBZNG54vNnl9Nulo5Jks1p0vK5Fx2qbzS4HACQRWIB+842ujBser0h78H9LxUTaNbV9pIhRFgCBIvh/ugIm868QCoHpIJ/Z7dNCbCAHIFAQWIB+8vevhEDDrY9vx9utNN4CCBAEFqCf9oRQw61PXvuOt/tONKiuqc3kagCAwAL0S11Tm78xdVIIjbCkxTuUmxorw5C2lzHKAsB8BBagH3z9KyOHRSsxOsLkagbW7BzvtNAn9LEACAAEFqAfzt2SP9T49mPZRh8LgABAYAH6IZQ2jDvfrPYRll3HatXm9phcDYBwR2AB+iGUR1jGpMYqMTpCLW0e7W1/nwBgFgIL0Edtbo8Ott/ROJQabn2sVosuHZUkSdpeyrQQAHMRWIA+KqluVJvbUGykTSOHRZtdzqDw3bl5e1mtuYUACHsEFqCP9lV6R1fGZ8TLYrGYXM3gOBtYGGEBYC4CC9BHvsAyMSP4b3jYnRnZibJYpKOnm1VV32J2OQDCGIEF6KN9J9pHWNJDN7DER0VoQvv7Y5QFgJkILEAf+UZYJoTwCIt0dnkzfSwAzERgAfqgyelS2akmSfKPQIQqfx8LK4UAmIjAAvTBgRNnJEmpcQ6lxDlMrmZwzWpf2rzrWJ2cLjaQA2AOAgvQB2eng+JMrmTw5abGalhMhJwuj3+jPAAYagQWoA/CoeHWx2Kx+KeFtjEtBMAkBBagD8JhSfO5zjbeElgAmIPAAvRBOI2wSPJv0V/MCAsAkxBYgF461ejUyYZWSeETWGaMTJLVIh2va1FlHRvIARh6BBagl3zTQdnJ0Yp12E2uZmjEOuyamOG9wSPTQgDMQGABemlfpXelzIT00LtD84Xk5bAfCwDzEFiAXtrXvgdLOCxpPtesnCRJjLAAMAeBBeil/Sd8e7CE1wiLb2nzZ8fq1epym1wNgHBDYAF6wTAM7fdtGhcmDbc+o5JjlBIbKafbo8+OsYEcgKFFYAF64XhdixpaXYqwWZSbGmt2OUPKYrHo0vZRlmKmhQAMMQIL0Au+htsxqXGKtIfft08eG8gBMEn4/cQF+mFfpbfhdnyY7HB7Pv8GcmW1ptYBIPwQWIBeOFDVvsPt8PBaIeQzfWSibFaLKupaVFHXbHY5AMIIgQXohYNV3hGWS8I0sMRE2v33T2KUBcBQIrAAPWQYhg61B5Zx6eEZWKSzy5vZQA7AUCKwAD1UUdeiRqdbdqtFOSnhtULoXP4+lvJaU+sAEF4ILEAP+aaDclJiFGEL328d3wjL7mN1cro8JlcDIFyE709doJfCvX/FJyclRsNiIuR0efR5BRvIARgaBBaghw6eJLBIbCAHwBwEFqCHGGE5a1Z7H8t2VgoBGCIEFqCHfCuELkkLz03jzsUIC4ChRmABeuB0o1M1jU5J0tjh4btCyGdGdpIsFuno6WZVNbSYXQ6AMEBgAXrA178yIilaMZF2k6sxX5zD7r9bNRvIARgKBBagB3z9K2PpX/HzTQtxI0QAQ4HAAvSAv+E2jcDiw40QAQwlAgvQA6wQ6sy3UmjX0Vq53GwgB2BwEViAHiCwdDYmNU4JUXa1tHn0RWWD2eUACHEEFuAimpwuHattlkRgOZfVatFMljcDGCIEFuAiDp9slCQlx0YqOTbS5GoCCxvIARgqfQosK1euVG5urqKiopSXl6fNmzdf8PxNmzYpLy9PUVFRGjNmjFavXt3pnNraWt17773KzMxUVFSUJk2apPXr1/elPGBA0XDbPTaQAzBUeh1Y1q1bp8WLF+vhhx9WcXGx5s6dq/nz56usrKzL80tKSnTTTTdp7ty5Ki4u1kMPPaT77rtPb775pv8cp9OpG2+8UUeOHNEbb7yhffv26eWXX9aIESP6/s6AAcKS5u7NzE6SJB2paVLNmVZziwEQ0nq9A9Zzzz2nu+++WwsXLpQkrVixQn/84x+1atUqLV++vNP5q1ev1qhRo7RixQpJ0qRJk7R161Y9++yzuu222yRJa9eu1alTp7RlyxZFRERIknJycvr6noABdaDK21BK/0pnidERumR4nA5WndGO8lp9eVK62SUBCFG9GmFxOp3atm2bCgoKOhwvKCjQli1bunxOUVFRp/PnzZunrVu3qq2tTZL0zjvvKD8/X/fee6/S09M1depUPfXUU3K73d3W0traqvr6+g4PYDCwQujCLm0fZWEDOQCDqVeBpbq6Wm63W+npHf8VlZ6ersrKyi6fU1lZ2eX5LpdL1dXVkqTDhw/rjTfekNvt1vr16/XII4/oX//1X/Xkk092W8vy5cuVmJjof2RnZ/fmrQA90ub2qLSmSRKBpTuzcnx9LLXmFgIgpPWp6dZisXT42DCMTscudv65xz0ej4YPH66XXnpJeXl5uv322/Xwww9r1apV3b7msmXLVFdX53+Ul5f35a0AF1Ra0yiXx1BMpE1ZiVFmlxOQfDve7iyvldtjmFsMgJDVqx6W1NRU2Wy2TqMpVVVVnUZRfDIyMro83263KyUlRZKUmZmpiIgI2Ww2/zmTJk1SZWWlnE6nIiM7LyV1OBxyOBy9KR/oNX/DbVrcBUN5OBs3PF5xDrvOtLq0/0SDJmUmmF0SgBDUqxGWyMhI5eXlqbCwsMPxwsJCzZkzp8vn5Ofndzp/w4YNmj17tr/B9qqrrtLBgwfl8Zzd3nv//v3KzMzsMqwAQ8UXWMYxHdQtm9WiGdmJkpgWAjB4ej0ltHTpUr3yyitau3at9u7dqyVLlqisrEyLFi2S5J2queuuu/znL1q0SKWlpVq6dKn27t2rtWvXas2aNXrggQf85/zgBz9QTU2N7r//fu3fv1+///3v9dRTT+nee+8dgLcI9B1LmntmFnduBjDIer2secGCBaqpqdHjjz+uiooKTZ06VevXr/cvQ66oqOiwJ0tubq7Wr1+vJUuW6MUXX1RWVpaef/55/5JmScrOztaGDRu0ZMkSTZ8+XSNGjND999+vn/zkJwPwFoG+O3iSFUI9cfbOzQQWAIPDYvg6YINcfX29EhMTVVdXp4QE5tDRfx6PoSmP/VHNbW796R+u1Vh2uu3WqUanZj3hnfrd8U83KimGqVwAPdPT39/cSwjoxvG6ZjW3uRVhsygnOcbscgJacmykclNjJUk7ymvNLQZASCKwAN3w9a+MTomV3ca3ysX4NpCj8RbAYOCnMNANdrjtnUtzaLwFMHgILEA3DtFw2yu+EZYd5bXysIEcgAFGYAG6wQhL70zMiFd0hE0NLS5/2AOAgUJgAbpx7i63uDi7zarpI9lADsDgILAAXag506rTTW2yWAgsvXEpG8gBGCQEFqALvtGVEUnRio60XeRs+MzybyBXa2odAEIPgQXoAjvc9o1vhGV/VYPqW9pMrgZAKCGwAF3wN9wyHdQrafEOZSdHyzCkXeV1ZpcDIIQQWIAusEKo7y7N9o6ycF8hAAOJwAJ0gcDSd74+FhpvAQwkAgtwnjOtLlXUtUgisPSFr4+luLxWIXJvVQABgMACnOdQ++hKalwkdx3ug0mZCXLYraptalNJdaPZ5QAIEQQW4DxsGNc/kXarpo1gAzkAA4vAApyHJc39dyl9LAAGGIEFOA8Nt/03y9fHwggLgAFCYAHO4+thGTc83uRKgpev8faLyno1trpMrgZAKCCwAOdwujwqPdUkiRGW/shIjFJWYpQ8hrTrKBvIAeg/AgtwjiM1jXJ7DMU57EpPcJhdTlA7u7yZPhYA/UdgAc7hXyE0PE4Wi8XkaoKbv/G2tNbUOgCEBgILcA7uITRwfCMsO8pPs4EcgH4jsADnYIXQwJmSlaAIm0XVZ5wqP9VsdjkAghyBBTgHgWXgREXYNCWrfQM5+lgA9BOBBWjn8Rg6XE1gGUhn+1gILAD6h8ACtDtW26yWNo8ibVZlD4s2u5yQMOucGyECQH8QWIB2vumg3NRY2W18awwE3wjL58fr1dLmNrcYAEGNn8pAO/pXBt6IpGgNj3fI5TG0+xgbyAHoOwIL0O7cPVgwMCwWi3+UpZgbIQLoBwIL0O5AVYMkRlgGmm8/FjaQA9AfBBZAkmEYbBo3SHyNt9vL2EAOQN8RWABJJ8+0qr7FJYtFGpMWa3Y5IWXaiETZrRZVNbTqeF2L2eUACFIEFkBn+1eyh8UoKsJmcjWhJTrSpkmZCZLoYwHQdwQWQNIhVggNKm6ECKC/CCyAWNI82M5uIMcIC4C+IbAAkg6epOF2MPlGWPYcq1eriw3kAPQegQUQe7AMtlHJMUqOjZTT7dGe4/VmlwMgCBFYEPbqW9p0or5VElNCg8VisWiWfwO5WlNrARCcCCwIe76G27R4hxKjI0yuJnRdes5+LADQWwQWhD3fdNA4RlcGlX+L/lICC4DeI7Ag7Pkbbgksg2rGyCTZrBYdr2vR8dpms8sBEGQILAh77MEyNGIddk1u30Du0yOnTK4GQLAhsCDscQ+hoXPZ6GRJ0tYjTAsB6B0CC8JaS5tbZaeaJDHCMhQuG+1tvGWEBUBvEVgQ1kqqG+UxpIQou9LiHWaXE/Jmt4+w7DvRoLqmNpOrARBMCCwIawd8K4TS42WxWEyuJvSlxTuUmxorw5C2lTHKAqDnCCwIa/SvDL2z00L0sQDoOQILwtrBqgZJ0rh0AstQ8U0LfVrCCAuAniOwIKxxD6Ghd3l7YNl1tE4tbdwIEUDPEFgQtlxuj0qqGyWxy+1QykmJUWqcQ063R7uO1pldDoAgQWBB2Co91aQ2t6HoCJuyEqPNLidsWCwWXZ7L8mYAvUNgQdg6cOLsDrdWKyuEhpJvAzkCC4CeIrAgbB3iHkKm8QWWbUdOy+0xTK4GQDAgsCBsHTjhXSFEYBl6EzPiFeewq6HVpX2VDWaXAyAIEFgQtrhLs3nsNqsuHZUkiWkhAD1DYEFY8ngM/5JmVgiZ43L6WAD0AoEFYelYbbNa2jyKtFk1KjnG7HLC0uxzAoth0McC4MIILAhLvtGV3NRY2W18G5hhZnaSImwWnahv1dHTzWaXAyDA8ZMaYcl/DyGmg0wTHWnT1BGJkqRP2KYfwEUQWBCWDlSxQigQ+PpYCCwALobAgrDECEtguHJMiiTp45IakysBEOj6FFhWrlyp3NxcRUVFKS8vT5s3b77g+Zs2bVJeXp6ioqI0ZswYrV69uttzX3/9dVksFt166619KQ24KMMwdMC3Qoi7NJtq9uhhslktKq1p0vFa+lgAdK/XgWXdunVavHixHn74YRUXF2vu3LmaP3++ysrKujy/pKREN910k+bOnavi4mI99NBDuu+++/Tmm292Ore0tFQPPPCA5s6d2/t3AvTQyYZWNbS4ZLV4m25hnvioCH8fy18YZQFwAb0OLM8995zuvvtuLVy4UJMmTdKKFSuUnZ2tVatWdXn+6tWrNWrUKK1YsUKTJk3SwoUL9b3vfU/PPvtsh/Pcbre+/e1v62c/+5nGjBlz0TpaW1tVX1/f4QH0hG90JSclVg67zeRqcOUYbx9L0SECC4Du9SqwOJ1Obdu2TQUFBR2OFxQUaMuWLV0+p6ioqNP58+bN09atW9XW1uY/9vjjjystLU133313j2pZvny5EhMT/Y/s7OzevBWEMfpXAou/j+UwjbcAuterwFJdXS2326309PQOx9PT01VZWdnlcyorK7s83+Vyqbq6WpL00Ucfac2aNXr55Zd7XMuyZctUV1fnf5SXl/fmrSCMsUIosMzO8faxlJ1q0jH6WAB0o09NtxaLpcPHhmF0Onax833HGxoa9J3vfEcvv/yyUlNTe1yDw+FQQkJChwfQE2zJH1g69LEcZloIQNfsvTk5NTVVNput02hKVVVVp1EUn4yMjC7Pt9vtSklJ0Z49e3TkyBHdfPPN/s97PB5vcXa79u3bp7Fjx/amTOCCmBIKPPljUrSzvFYfH67RN2aNNLscAAGoVyMskZGRysvLU2FhYYfjhYWFmjNnTpfPyc/P73T+hg0bNHv2bEVERGjixInavXu3duzY4X/ccsstuv7667Vjxw56UzCgTjc6VX3GKUkam0ZgCRS+xlv6WAB0p1cjLJK0dOlS3XnnnZo9e7by8/P10ksvqaysTIsWLZLk7S05duyYXn31VUnSokWL9MILL2jp0qW65557VFRUpDVr1ui1116TJEVFRWnq1KkdvkZSUpIkdToO9NfBk97RlRFJ0Yp19PqvPwbJ7NHJHfpYRiRFm10SgADT65/YCxYsUE1NjR5//HFVVFRo6tSpWr9+vXJyciRJFRUVHfZkyc3N1fr167VkyRK9+OKLysrK0vPPP6/bbrtt4N4F0EP7T9BwG4jiHHZNG5GoHeW1+vhQjW7LY1oIQEcWI0Tu615fX6/ExETV1dXRgItu/fSdPfrPLUf0/WvG6KGbJpldDs7x9LtfaPWmQ/rrvJH6+V/PMLscAEOkp7+/uZcQwsq+Su8Iy/j0eJMrwfnyx3JfIQDdI7AgrPimhMZzD6GA49uPpfxUs46ebjK7HAABhsCCsFF9plU1jU5ZLPSwBKJYh13TR/r2Y2G1EICOCCwIG77RlVHJMYqJZIVQIPJt0//RoWqTKwEQaAgsCBv76V8JeFdf4t3t+qOD1QqR9QAABgiBBWFj3wnvHiz0rwSuvJxhctitOlHfqkPte+YAgERgQRg5cIIRlkAXFWHTZaO9u95+eIBpIQBnEVgQFgzD0L72wDIhg8ASyK5qnxb68CCBBcBZBBaEhcr6FjW0uGS3WjQmlSmhQDZ3nDewfHz4lNrcHpOrARAoCCwIC74N40anxirSzl/7QDY5M0FJMRE60+rSrqO1ZpcDIEDwkxth4UB7w+0E+lcCntVq0VVjvaMsm+ljAdCOwIKwsI+G26By9bizy5sBQCKwIEzs9zfc0r8SDHz7sRSX1epMq8vkagAEAgILQp7HY/inhMYxwhIUspNjNCo5Ri6Pob8c5maIAAgsCANHTzeruc2tSLtVOckxZpeDHmJ5M4BzEVgQ8nz9K5ekxclu4698sPAtb2YDOQASgQVhYD8bxgWl/DEpslikA1VnVFnXYnY5AExGYEHI8wWWcdxDKKgMi43UtBGJkqQP9p80uRoAZiOwIOT5No1jD5bgc934NEnSxv1VJlcCwGwEFoS0NrdHh082SmIPlmB07YThkrwbyLnYph8IawQWhLRDJ8/I6fYo3mHXyGHRZpeDXpqZnaSkmAg1tLi0vazW7HIAmIjAgpC2t6JekjQxM14Wi8XkatBbNqtF14xrnxbax7QQEM4ILAhpX1R4+1cmZSaYXAn66roJvsBC4y0QzggsCGmf+0ZYMggsweqa9sbbzyvqVVXP8mYgXBFYENL2+kdYaLgNVqlxDk0f6V3evJHlzUDYIrAgZJ1saFX1mVZZLGwaF+x8y5s3MS0EhC0CC0LWF5Xe6aDRKbGKibSbXA364+zy5pMsbwbCFIEFIcu3QojpoODnW95c3+JScXmt2eUAMAGBBSHLv0KIhtugd+7y5ve/YHkzEI4ILAhZ/hVCLGkOCb7lze/TxwKEJQILQpLT5dGhk2ckMSUUKq6bMFxWi3eq7+jpJrPLATDECCwISYdOnlGb21B8lF0jktiSPxQkx0Zqdk6yJOlPe5kWAsINgQUhyd9wm5HAlvwh5IbJ3tVChZ+fMLkSAEONwIKQ9EUlG8aFohsnZ0iSPj5co/qWNpOrATCUCCwISXtpuA1JuamxGpsWK5fHYBM5IMwQWBCSzu7BQmAJNTdMTpfEtBAQbggsCDlVDS2qPuP0bsmfzpRQqCloDyzv76tSG7veAmGDwIKQs+eYd3RlbFqcoiNtJleDgTYze5hSYiPV0OLSpyWnzC4HwBAhsCDkfHasTpI0NYvpoFBks1r05Une1UIbmBYCwgaBBSHns+PtgWVEosmVYLDcMOlsH4thGCZXA2AoEFgQcj5rnxIisISuuePSFB1h07HaZv//bwChjcCCkHK60aljtc2SpMlMCYWs6EibvjTROy20/rMKk6sBMBQILAgpvumg0SkxSoiKMLkaDKb507ybyL27u4JpISAMEFgQUnzTA1OYDgp5108YLofdqiM1Tdpb0WB2OQAGGYEFIcU3wjKNwBLyYh12XTchTZL0LtNCQMgjsCCknF3STGAJBzdNy5Qk/Z5pISDkEVgQMuqa21Ra0yRJmkLDbVj40sThirRZdfhko/afOGN2OQAGEYEFIePz497+lRFJ0RoWG2lyNRgK8VERumZ8qiRp/W6mhYBQRmBByNjj3zCO0ZVwMn+qd1qIPhYgtBFYEDJ2H6PhNhzdMDldETaL9p84o/0nWC0EhCoCC0KGL7BMoeE2rCRGR+ja8d5N5H6745jJ1QAYLAQWhIS65jYdPtkoSZo+ksASbr42M0uS9Nsdx1ktBIQoAgtCwq6jtZKkUckxSolzmFsMhtwNk9IVG2nT0dPN2lZ62uxyAAwCAgtCws7yWknSjOwkU+uAOaIjbZo31btV/2+YFgJCEoEFIWFHe2CZSWAJW7fOHCFJ+v2uCrW5PSZXA2CgEVgQ9AzD0I5yb8PtzGz6V8LVnLEpSo1z6HRTmzYfOGl2OQAGGIEFQe9YbbOqz7TKbrWwQiiM2W1W3TzDuyfLb4qPm1wNgIFGYEHQ29k+ujIxM15RETaTq4GZfNNChZ+fUGOry+RqAAwkAguC3s72FUL0r2D6yETlpsaquc3NVv1AiCGwIOjtKKuVJM0YmWRqHTCfxWLRN/NGSpJ+vfWoydUAGEgEFgQ1l9vj3+GWERZI0m2zRspqkT45ckqHT3IHZyBU9CmwrFy5Urm5uYqKilJeXp42b958wfM3bdqkvLw8RUVFacyYMVq9enWHz7/88suaO3euhg0bpmHDhumGG27QJ5980pfSEGb2nzij5ja34hx2jU2LM7scBICMxChdOz5NkvTrbYyyAKGi14Fl3bp1Wrx4sR5++GEVFxdr7ty5mj9/vsrKyro8v6SkRDfddJPmzp2r4uJiPfTQQ7rvvvv05ptv+s/ZuHGj7rjjDr3//vsqKirSqFGjVFBQoGPH2AAKF+brX5k+MlFWq8XcYhAwFlyWLUl6c9tRudiTBQgJFqOXN9644oorNGvWLK1atcp/bNKkSbr11lu1fPnyTuf/5Cc/0TvvvKO9e/f6jy1atEg7d+5UUVFRl1/D7XZr2LBheuGFF3TXXXd1eU5ra6taW1v9H9fX1ys7O1t1dXVKSEjozVtCEPvHX+/Ur7cd1Q+vG6sff2Wi2eUgQDhdHuUv/5NqGp165a7ZumFyutklAehGfX29EhMTL/r7u1cjLE6nU9u2bVNBQUGH4wUFBdqyZUuXzykqKup0/rx587R161a1tbV1+Zympia1tbUpOTm521qWL1+uxMRE/yM7O7s3bwUhYmv7fWNmjx5mciUIJJF2q75+qXeJ8/9uLTe5GgADoVeBpbq6Wm63W+npHf+1kp6ersrKyi6fU1lZ2eX5LpdL1dXVXT7nwQcf1IgRI3TDDTd0W8uyZctUV1fnf5SX80Mp3FSfaVVJtfcOzXmjug+3CE9/0z4t9OcvqnSyofUiZwMIdH1qurVYOvYKGIbR6djFzu/quCQ988wzeu211/TWW28pKiqq29d0OBxKSEjo8EB42XrEO7oyPj1OiTERJleDQDM+PV4zs5Pk8hh6g+ZbIOj1KrCkpqbKZrN1Gk2pqqrqNIrik5GR0eX5drtdKSkpHY4/++yzeuqpp7RhwwZNnz69N6UhDG0rPSVJysthdAVd+9YVoyRJ//1xqdyeXrXrAQgwvQoskZGRysvLU2FhYYfjhYWFmjNnTpfPyc/P73T+hg0bNHv2bEVEnP1X8c9//nM98cQT+sMf/qDZs2f3piyEKV//ymX0r6Abt8zIUlJMhI7VNuv9L6rMLgdAP/R6Smjp0qV65ZVXtHbtWu3du1dLlixRWVmZFi1aJMnbW3Luyp5FixaptLRUS5cu1d69e7V27VqtWbNGDzzwgP+cZ555Ro888ojWrl2r0aNHq7KyUpWVlTpzhk2f0LWWNrc+a98wbjYjLOhGVIRNC2Z7e1l+WXTE3GIA9EuvA8uCBQu0YsUKPf7445o5c6Y++OADrV+/Xjk5OZKkioqKDnuy5Obmav369dq4caNmzpypJ554Qs8//7xuu+02/zkrV66U0+nUN7/5TWVmZvofzz777AC8RYSineW1anMbSot3KDs52uxyEMC+c2WOLBZp84Fqdr4Fgliv92EJVD1dx43Q8OL7B/XzP+7TTdMytPLbeWaXgwB3939+qj99UaXvXZWrf7p5stnlADjHoOzDAgSKrUdouEXP3ZnvHQH+9bZyNTldJlcDoC8ILAg6LrdHn7Yvab58NIEFF3fNuDSNTolRQ4tLbxdzyw8gGBFYEHQ+O16vM60uJUTZNTmL6T9cnNVq0Z35oyVJaz4skYclzkDQIbAg6BQdqpEkXZ6bIhs3PEQPLbgsW/FRdh0+2aj39p4wuxwAvURgQdD5+LA3sOSPTbnImcBZcQ67vnOlt5flpQ8Om1wNgN4isCCotLk9+rS94TZ/DIEFvfN/5oxWpM2qraWn/TslAwgOBBYElV1H69TkdGtYTIQmZsSbXQ6CzPCEKP9dnP99E6MsQDAhsCCo+KaDrshNkZX+FfTBPdeMkSQV7j2hg1VsJAcECwILgoqv4Zb+FfTVJcPjdMOkdBmG9DK9LEDQILAgaLS63Nra3ndAYEF//OA67yjLm9uPqvxUk8nVAOgJAguCxrYjp9XS5lFavEPjhseZXQ6CWF5Osq6+JFUuj6GVGw+aXQ6AHiCwIGhs2n9SknfXUouF/hX0z+IbxkmSfr2VURYgGBBYEDT8gWV8qsmVIBTMHp2sueO8oywvvs8oCxDoCCwICifqW/RFZYMsFmnuuDSzy0GI8I2yvLGNURYg0BFYEBR8oyvTRyQqOTbS5GoQKvJyzo6y/N/C/WaXA+ACCCwICh+0B5ZrxzO6goH1j/MmSJLe3nFMe47XmVwNgO4QWBDw3B5Dmw9US5KuIbBggE0fmaRbZmTJMKSn3/3C7HIAdIPAgoC382it6prbFB9l18zsJLPLQQj6x3kTFGGzaPOBav9oHoDAQmBBwHvv8xOSvMuZ7Tb+ymLgZSfH6M4rR0uSlr/7hTwew9yCAHTCT38EvA3tgaVgSrrJlSCU/f2XLlF8lF17K+r1v1vLzS4HwHkILAhoh0+e0cGqM4qwWXT9xOFml4MQNiw2UotvGC9J+pc/fKHaJqfJFQE4F4EFAa2wfXTlyjEpSoiKMLkahLq78nM0Pj1Op5va9K8bWOYMBBICCwKafzpoMtNBGHwRNqt+dstUSdKv/lKqz46xzBkIFAQWBKyTDa3aXnZaknQDgQVDJH9sim6ekSWPIT3628/kpgEXCAgEFgSs9/aekGFIM0YmKjMx2uxyEEYevmmS4hx2FZfV6tWiI2aXA0AEFgSw3++qkCQVTMkwuRKEm4zEKD04f6Ik6Zk/7OM+Q0AAILAgIFU1tGjLIe/utrfMyDK5GoSjb10+SpfnJqu5za1lb+2WYTA1BJiJwIKA9PtdFfIY0qWjkpSdHGN2OQhDVqtF/3LbdDnsVn14sFqvf8reLICZCCwISL/dcVyS9DVGV2Ci3NRYPVDgvTniE7/7XCXVjSZXBIQvAgsCTmlNo3aU18pqkb46ncACc919da7yx6SoyenW/a8Xy+nymF0SEJYILAg4/2+nd3TlqktSlRbvMLkahDur1aLnFsxQYnSEdh2t0/99jw3lADMQWBBQDMPQW8XHJNFsi8CRmRitf7ltmiRp9aZD3NEZMAGBBQHlk5JTOnyyUTGRNs2flml2OYDfV6Zm6o7LR8kwpPteL2apMzDECCwIKL6VGLfMyFKcw25yNUBHj908WdNHJqq2qU0/+NU2tbS5zS4JCBsEFgSMuqY2rd/t3Szu9stHmVwN0FlUhE2rvpOn5NhIfXasXo/85jP2ZwGGCIEFAePt4qNqdXk0MSNeM0Ymml0O0KURSdF64Y5LZbVIb2w7qlWbDpldEhAWCCwICB6PoVeLSiVJ37pilCwWi8kVAd2bc0mqfnrLFEnerfvfaV/ZBmDwEFgQEDbur9Lh6kbFR9l126yRZpcDXNRd+aO18OpcSdID/7tTfzlcY3JFQGgjsCAgrPmwRJJ0+2XZiqXZFkHioZsm6StTMuR0e7Twl1u162it2SUBIYvAAtN9UVmvjw7WyGqRvjtntNnlAD1mtVq04vaZuiI3WQ2tLt219hN9UVlvdllASCKwwHQvvu9tWvzK1AyNHMaNDhFcoiJsWvO3l+nSUUmqbWrTd175iw6caDC7LCDkEFhgqoNVDfrdLm/D4r3XX2JyNUDfxDns+s//c7mmZCWo+oxTf/3vRdpZXmt2WUBIIbDAVP/254MyDOnGyemaksVSZgSvxOgI/ffdV2hGtnek5Vsvf6wth6rNLgsIGQQWmOZgVYP/Rof3f3mcydUA/TcsNlK/WniF5oxNUaPTrb/9j0/9myEC6B8CC0zz1Pov5DGkgsnpmjqC0RWEhjiHXWv/9jIVTE6X0+XRD3+1XSve2y+Phx1xgf4gsMAUmw+c1J+/qJLdatGD8yeaXQ4woKIibFr57Vn63lXefVpWvHdAP3ptu5qcLpMrA4IXgQVDrs3t0ZO/3yvJu/nWmLQ4kysCBp7dZtU/3TxZz9w2XRE2i9bvrtStL36kfZWsIAL6gsCCIffy5sP6orJBSTERuu/LrAxCaPuby7L12j1XKi3eof0nzuiWFz7Uf31cyk0TgV4isGBIlVQ3asV7ByRJj351spJiIk2uCBh8s0cn69375+q6CWlqdXn06G8+08JfblVFXbPZpQFBg8CCIeNye/TjN3bK6fJo7rhUfWPWCLNLAoZMapxDa797mR79q8mKsFn0py+qdONzH+i/io7QkAv0AIEFQ+YXfzqgT4+cVpzDrqe+Po07MiPsWK0W3X11rn7393N16agknWl16dHf7tE3V29Rcdlps8sDAhqBBUPiwwPVeuH9g5Kkp74xTdnJbMGP8DUhI15vLJqjn90yRbGRNm0vq9XXV27R379WrPJTTWaXBwQkAgsG3cGqBv3wV9tkGNIdl2frlhlZZpcEmM5mtei7c0brzw9cp7/OGymLRfp/O4/rS/+6Ucve2q2yGoILcC6LESKt6vX19UpMTFRdXZ0SEhLMLgftTja06hurPlL5qWbNGpWk/7nnSkVF2MwuCwg4e47X6an1e/XRwRpJ3kDztRlZ+rtrx2pCRrzJ1QGDp6e/vwksGDRVDS361st/0cGqMxqVHKO3fzhHKXEOs8sCAtonJaf0wvsH9cH+k/5jl+cm684rczRvSoYi7QyMI7QQWGCq8lNN+tv/+ESHTjYqMzFKr91zpUanxppdFhA0dh2t1aqNh7Th8xNyt68iSo2L1M0zsnTLjCzNzE6icR0hgcAC02w9ckp/91/bVNPoVFZilF77/pXKSSGsAH1RUdes1z4p1+uflKmqodV/fFRyjG6ekakbJqVr+sgk2ayEFwQnAguGnMvt0aqNh/SLPx2Qy2NoSlaCXvnubGUmRptdGhD02twebdp3Uu/sPK7Cz0+ouc3t/1xybKSuHZ+m6yakKX9MioYnRJlYKdA7BBYMqb8crtETv/9cnx2rlyT91fRMPfPN6YqJtJtcGRB6mpwuvbe3Su/urtCHB6rV0Nrxpoo5KTG6bHSyLhs9TLNGDdOYtDhGYBCwCCwYdIZhaFvpaf37B4dV+PkJSVK8w67Hb52iW2eOYH4dGAJtbo+2lZ7W+/uq9MH+an1RWa/zf6pHRVg1ISNBU7ISNDkzQRMy4pWbGquU2Ei+T2G6QQ0sK1eu1M9//nNVVFRoypQpWrFihebOndvt+Zs2bdLSpUu1Z88eZWVl6cc//rEWLVrU4Zw333xTjz76qA4dOqSxY8fqySef1Ne//vUe10RgGTqVdS1av7tCbxcf0+5jdZIkq0W64/JRWnLjeKWyEggwTX1Lm7aVntbWI6f0aclpfXa8Tk1Od5fnxjvsyk2L1eiUWI1OiVFWUrQyEqOUmej9b0KUnUCDQdfT39+9Hq9ft26dFi9erJUrV+qqq67Sv//7v2v+/Pn6/PPPNWrUqE7nl5SU6KabbtI999yj//7v/9ZHH32kH/7wh0pLS9Ntt90mSSoqKtKCBQv0xBNP6Otf/7refvtt/c3f/I0+/PBDXXHFFb0tEQOo1eXWkeom7TvRoK1HTumTklP6orLB//lIu1VfnzlCC+fmalw6e0UAZkuIitD1E4br+gnDJUluj6HSmkZ9XlGvPcfr9fnxeh2sOqPjdc1qaHVp19E67Tpa1+VrxUTalJEYpdQ4h4bFRCg5NlLDYiKVHBuppJhIJcdGKCEqQrEOu+IcdsU67Ip12OSws9cSBl6vR1iuuOIKzZo1S6tWrfIfmzRpkm699VYtX7680/k/+clP9M4772jv3r3+Y4sWLdLOnTtVVFQkSVqwYIHq6+v17rvv+s/5yle+omHDhum1117rUV2DNcJS1dAip8vTYYjV92dDxnkf+z5vnPex/5ndnN/59br7GgNVg9NlqLnNpWanR01Ol5rb3DrT6lLNGadONrTqZEOrKutbVHaqyb+k8lyXjR6mm6Zl6pYZWeytAgShlja3yk416fDJRh2paVRpTaMq6lpUWdeiyvoW1Ta19fm1I2wWb3iJ9AYZR4RVDrtVDrtNkXbfn63tfz73mE12m0U2q0V2q/e/vofdapHVYmn/vFU2S+fzrBaLrBZJFskiiywWySLJYrnAn6X2j9uPd/dn33ntz7mYnoxM9XTsqieDXD2paiAGy4YnOAY8kA7KCIvT6dS2bdv04IMPdjheUFCgLVu2dPmcoqIiFRQUdDg2b948rVmzRm1tbYqIiFBRUZGWLFnS6ZwVK1Z0W0tra6taW88u8auvr+/NW+mxv/uvbSouqx2U1w4W8Q67xg6P08zsJF2Rm6zLcpOZ9gGCXFSETePT4zW+m5HRZqdblfUtqqhr1qlGp043tel0o7P9z96PTzW2qqHFpcZWl860utTS5pEktbkN1Ta19Sv0IDC99cM5mjVqmClfu1eBpbq6Wm63W+np6R2Op6enq7KyssvnVFZWdnm+y+VSdXW1MjMzuz2nu9eUpOXLl+tnP/tZb8rvk0ibVVER3p0lfQnWl1J9YdWXpP3htbvPX+R5lvNeoPP5vo8v/Ho6//wunmezWhQTaVNMhF3RkTbFRNoUHWlTWpxDafFnH2NS45Se4GAeGwgz0ZE25abGKrcXGz663B41Ot1qbD0bYhpb3XK63Wpt86jV5ZHT5VGry61Wl8f/OPeY223I5THk9njkNuT9r8eQ2+M7fvbPnvOOeQxDhiHvf+UdaTYkyVCHj33n+Uem/ed7R609/s+1n3fuc7sYcT7XRacsLnLCxZ5/sUmRiz//Yl//widYTfxd0Kc1p+f/8jIM44K/0Lo6//zjvX3NZcuWaenSpf6P6+vrlZ2dffHie2nd3+UP+GsCQCiy26xKjLYqMTrC7FIQgnoVWFJTU2Wz2TqNfFRVVXUaIfHJyMjo8ny73a6UlJQLntPda0qSw+GQw8G0BAAA4aBXd9GKjIxUXl6eCgsLOxwvLCzUnDlzunxOfn5+p/M3bNig2bNnKyIi4oLndPeaAAAgvPR6Smjp0qW68847NXv2bOXn5+ull15SWVmZf1+VZcuW6dixY3r11VcleVcEvfDCC1q6dKnuueceFRUVac2aNR1W/9x///265ppr9C//8i/62te+pt/+9rd677339OGHHw7Q2wQAAMGs14FlwYIFqqmp0eOPP66KigpNnTpV69evV05OjiSpoqJCZWVl/vNzc3O1fv16LVmyRC+++KKysrL0/PPP+/dgkaQ5c+bo9ddf1yOPPKJHH31UY8eO1bp169iDBQAASGJrfgAAYKKe/v7uVQ8LAACAGQgsAAAg4BFYAABAwCOwAACAgEdgAQAAAY/AAgAAAh6BBQAABDwCCwAACHgEFgAAEPB6vTV/oPJt2FtfX29yJQAAoKd8v7cvtvF+yASWhoYGSVJ2drbJlQAAgN5qaGhQYmJit58PmXsJeTweHT9+XPHx8bJYLAP2uvX19crOzlZ5eTn3KBpEXOehw7UeGlznocF1HhqDeZ0Nw1BDQ4OysrJktXbfqRIyIyxWq1UjR44ctNdPSEjgm2EIcJ2HDtd6aHCdhwbXeWgM1nW+0MiKD023AAAg4BFYAABAwCOwXITD4dBjjz0mh8Nhdikhjes8dLjWQ4PrPDS4zkMjEK5zyDTdAgCA0MUICwAACHgEFgAAEPAILAAAIOARWAAAQMAjsFzEypUrlZubq6ioKOXl5Wnz5s1mlxTUPvjgA918883KysqSxWLRb37zmw6fNwxDP/3pT5WVlaXo6Ghdd9112rNnjznFBrHly5frsssuU3x8vIYPH65bb71V+/bt63AO17r/Vq1apenTp/s308rPz9e7777r/zzXeHAsX75cFotFixcv9h/jWvffT3/6U1kslg6PjIwM/+fNvsYElgtYt26dFi9erIcffljFxcWaO3eu5s+fr7KyMrNLC1qNjY2aMWOGXnjhhS4//8wzz+i5557TCy+8oE8//VQZGRm68cYb/feKQs9s2rRJ9957rz7++GMVFhbK5XKpoKBAjY2N/nO41v03cuRIPf3009q6dau2bt2qL33pS/ra177m/yHONR54n376qV566SVNnz69w3Gu9cCYMmWKKioq/I/du3f7P2f6NTbQrcsvv9xYtGhRh2MTJ040HnzwQZMqCi2SjLffftv/scfjMTIyMoynn37af6ylpcVITEw0Vq9ebUKFoaOqqsqQZGzatMkwDK71YBo2bJjxyiuvcI0HQUNDgzFu3DijsLDQuPbaa43777/fMAz+Pg+Uxx57zJgxY0aXnwuEa8wISzecTqe2bdumgoKCDscLCgq0ZcsWk6oKbSUlJaqsrOxwzR0Oh6699lqueT/V1dVJkpKTkyVxrQeD2+3W66+/rsbGRuXn53ONB8G9996rr371q7rhhhs6HOdaD5wDBw4oKytLubm5uv3223X48GFJgXGNQ+bmhwOturpabrdb6enpHY6np6ersrLSpKpCm++6dnXNS0tLzSgpJBiGoaVLl+rqq6/W1KlTJXGtB9Lu3buVn5+vlpYWxcXF6e2339bkyZP9P8S5xgPj9ddf1/bt2/Xpp592+hx/nwfGFVdcoVdffVXjx4/XiRMn9M///M+aM2eO9uzZExDXmMByERaLpcPHhmF0OoaBxTUfWD/60Y+0a9cuffjhh50+x7XuvwkTJmjHjh2qra3Vm2++qe9+97vatGmT//Nc4/4rLy/X/fffrw0bNigqKqrb87jW/TN//nz/n6dNm6b8/HyNHTtWv/zlL3XllVdKMvcaMyXUjdTUVNlstk6jKVVVVZ0SJgaGrxudaz5w/v7v/17vvPOO3n//fY0cOdJ/nGs9cCIjI3XJJZdo9uzZWr58uWbMmKFf/OIXXOMBtG3bNlVVVSkvL092u112u12bNm3S888/L7vd7r+eXOuBFRsbq2nTpunAgQMB8feZwNKNyMhI5eXlqbCwsMPxwsJCzZkzx6SqQltubq4yMjI6XHOn06lNmzZxzXvJMAz96Ec/0ltvvaU///nPys3N7fB5rvXgMQxDra2tXOMB9OUvf1m7d+/Wjh07/I/Zs2fr29/+tnbs2KExY8ZwrQdBa2ur9u7dq8zMzMD4+zwkrb1B6vXXXzciIiKMNWvWGJ9//rmxePFiIzY21jhy5IjZpQWthoYGo7i42CguLjYkGc8995xRXFxslJaWGoZhGE8//bSRmJhovPXWW8bu3buNO+64w8jMzDTq6+tNrjy4/OAHPzASExONjRs3GhUVFf5HU1OT/xyudf8tW7bM+OCDD4ySkhJj165dxkMPPWRYrVZjw4YNhmFwjQfTuauEDINrPRD+4R/+wdi4caNx+PBh4+OPPzb+6q/+yoiPj/f/zjP7GhNYLuLFF180cnJyjMjISGPWrFn+ZaHom/fff9+Q1Onx3e9+1zAM79K5xx57zMjIyDAcDodxzTXXGLt37za36CDU1TWWZPzHf/yH/xyudf9973vf8/98SEtLM7785S/7w4phcI0H0/mBhWvdfwsWLDAyMzONiIgIIysry/jGN75h7Nmzx/95s6+xxTAMY2jGcgAAAPqGHhYAABDwCCwAACDgEVgAAEDAI7AAAICAR2ABAAABj8ACAAACHoEFAAAEPAILAAAIeAQWAAAQ8AgsAAAg4BFYAASsDz/8UBEREWptbfUfKykpkcViUWlpqYmVARhqBBYAAWvHjh2aNGmSHA5Hh2NJSUnKyckxsTIAQ43AAiBg7dy5U5deemmHYzt27NCMGTNMqgiAWQgsAALWjh07NHPmzA7HiouLCSxAGCKwAAhIbrdbe/bs6TTCsn379k4hBkDoI7AACEj79u1Tc3OzsrKy/MeKiop07NgxRliAMERgARCQduzYIUn6t3/7Nx04cEDvvvuu7rrrLknqsGoIQHggsAAISDt27NCNN96okpISTZ06VQ899JCefvppJSQk6MUXXzS7PABDzGIYhmF2EQBwvnnz5mnWrFlavny52aUACACMsAAISDt37tT06dPNLgNAgCCwAAg4lZWVOnHiBIEFgB9TQgAAIOAxwgIAAAIegQUAAAQ8AgsAAAh4BBYAABDwCCwAACDgEVgAAEDAI7AAAICAR2ABAAABj8ACAAACHoEFAAAEvP8Pi5AZG0xgAUYAAAAASUVORK5CYII=",
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAigAAAGdCAYAAAA44ojeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAjMklEQVR4nO3de3CU1eH/8c82NyBNtoSQ3WyJEEfwlohtUCBj5R7IiCg4AxabQqUOlkuNgAjaGcN3HAJ0BJwidHQYUBRhdABtoZY4QIQGKkQZuXhBJQpDlggNuwmNmxjO7w9+bF0SLhuS7El4v2aeGfd5zu6e5wySN0/24jDGGAEAAFjkJ5GeAAAAwMUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWiY70BJri3LlzOnHihBISEuRwOCI9HQAAcBWMMaqqqpLH49FPfnL5ayRtMlBOnDihtLS0SE8DAAA0wbFjx9StW7fLjmmTgZKQkCDp/AkmJiZGeDYAAOBq+P1+paWlBX+OX06bDJQLv9ZJTEwkUAAAaGOu5uUZvEgWAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWiY70BADAGgXOSM/gvAJfpGcARBxXUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYJK1BWrFihO+64Q4mJiUpMTFT//v31j3/8I3jcGKOCggJ5PB517NhRAwcO1KFDh0IeIxAIaPr06UpOTlZ8fLxGjRql48ePN8/ZAACAdiGsQOnWrZsWLFigffv2ad++fRo8eLAeeOCBYIQsWrRIixcv1rJly7R371653W4NGzZMVVVVwcfIz8/Xxo0btW7dOu3atUvV1dUaOXKk6uvrm/fMAABAm+UwxphreYCkpCT9+c9/1qOPPiqPx6P8/Hw9/fTTks5fLXG5XFq4cKEmT54sn8+nrl27as2aNRo3bpwk6cSJE0pLS9OWLVs0fPjwq3pOv98vp9Mpn8+nxMTEa5k+APxPgTPSMzivwBfpGQAtIpyf39FNfZL6+nq99dZbOnv2rPr376+jR4/K6/UqJycnOCYuLk4DBgxQSUmJJk+erNLSUtXV1YWM8Xg8ysjIUElJySUDJRAIKBAIhJwggHbEljAAYI2wXyR74MAB/fSnP1VcXJwef/xxbdy4Ubfddpu8Xq8kyeVyhYx3uVzBY16vV7GxsercufMlxzSmsLBQTqczuKWlpYU7bQAA0IaEHSg333yz9u/frz179ugPf/iDJkyYoMOHDwePOxyOkPHGmAb7LnalMXPnzpXP5wtux44dC3faAACgDQk7UGJjY3XTTTepT58+KiwsVO/evfXiiy/K7XZLUoMrIRUVFcGrKm63W7W1taqsrLzkmMbExcUF3zl0YQMAAO3XNX8OijFGgUBA6enpcrvdKioqCh6rra1VcXGxsrOzJUlZWVmKiYkJGVNeXq6DBw8GxwAAAIT1ItlnnnlGubm5SktLU1VVldatW6cdO3bovffek8PhUH5+vubPn6+ePXuqZ8+emj9/vjp16qTx48dLkpxOpyZNmqSZM2eqS5cuSkpK0qxZs5SZmamhQ4e2yAkCAIC2J6xAOXnypPLy8lReXi6n06k77rhD7733noYNGyZJmj17tmpqajRlyhRVVlaqb9++2rp1qxISEoKPsWTJEkVHR2vs2LGqqanRkCFDtHr1akVFRTXvmQEAgDbrmj8HJRL4HBSgneFtxqH4HBS0U+H8/Oa7eAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYJ6xAKSws1F133aWEhASlpKTowQcf1Oeffx4yZuLEiXI4HCFbv379QsYEAgFNnz5dycnJio+P16hRo3T8+PFrPxsAANAuhBUoxcXFmjp1qvbs2aOioiL98MMPysnJ0dmzZ0PGjRgxQuXl5cFty5YtIcfz8/O1ceNGrVu3Trt27VJ1dbVGjhyp+vr6az8jAADQ5kWHM/i9994Lub1q1SqlpKSotLRU9957b3B/XFyc3G53o4/h8/m0cuVKrVmzRkOHDpUkvf7660pLS9P777+v4cOHh3sOAACgnbmm16D4fD5JUlJSUsj+HTt2KCUlRb169dJjjz2mioqK4LHS0lLV1dUpJycnuM/j8SgjI0MlJSWNPk8gEJDf7w/ZAABA+9XkQDHGaMaMGbrnnnuUkZER3J+bm6s33nhD27Zt0wsvvKC9e/dq8ODBCgQCkiSv16vY2Fh17tw55PFcLpe8Xm+jz1VYWCin0xnc0tLSmjptAADQBoT1K54fmzZtmj755BPt2rUrZP+4ceOC/52RkaE+ffqoe/fu2rx5s8aMGXPJxzPGyOFwNHps7ty5mjFjRvC23+8nUgAAaMeadAVl+vTpevfdd7V9+3Z169btsmNTU1PVvXt3HTlyRJLkdrtVW1urysrKkHEVFRVyuVyNPkZcXJwSExNDNgAA0H6FFSjGGE2bNk0bNmzQtm3blJ6efsX7nD59WseOHVNqaqokKSsrSzExMSoqKgqOKS8v18GDB5WdnR3m9AEAQHsU1q94pk6dqrVr1+qdd95RQkJC8DUjTqdTHTt2VHV1tQoKCvTQQw8pNTVVZWVleuaZZ5ScnKzRo0cHx06aNEkzZ85Uly5dlJSUpFmzZikzMzP4rh4AAHB9CytQVqxYIUkaOHBgyP5Vq1Zp4sSJioqK0oEDB/Taa6/pzJkzSk1N1aBBg7R+/XolJCQExy9ZskTR0dEaO3asampqNGTIEK1evVpRUVHXfkYAAKDNcxhjTKQnES6/3y+n0ymfz8frUYD2oMAZ6RnYpcAX6RkALSKcn998Fw8AALAOgQIAAKxDoAAAAOsQKAAAwDoECgAAsA6BAgAArEOgAAAA6xAoAADAOgQKAACwDoECAACsQ6AAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOsQKAAAwDoECgAAsA6BAgAArEOgAAAA6xAoAADAOgQKAACwDoECAACsQ6AAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOsQKAAAwDoECgAAsE50pCcAAAjVY87mVnuusgX3tdpzAeHgCgoAALAOgQIAAKxDoAAAAOsQKAAAwDoECgAAsA6BAgAArEOgAAAA6xAoAADAOgQKAACwDoECAACsQ6AAAADrECgAAMA6BAoAALBOWIFSWFiou+66SwkJCUpJSdGDDz6ozz//PGSMMUYFBQXyeDzq2LGjBg4cqEOHDoWMCQQCmj59upKTkxUfH69Ro0bp+PHj1342AACgXQgrUIqLizV16lTt2bNHRUVF+uGHH5STk6OzZ88GxyxatEiLFy/WsmXLtHfvXrndbg0bNkxVVVXBMfn5+dq4caPWrVunXbt2qbq6WiNHjlR9fX3znRkAAGizHMYY09Q7f/fdd0pJSVFxcbHuvfdeGWPk8XiUn5+vp59+WtL5qyUul0sLFy7U5MmT5fP51LVrV61Zs0bjxo2TJJ04cUJpaWnasmWLhg8ffsXn9fv9cjqd8vl8SkxMbOr0AdiiwBnpGVilx/drW+25yhbc12rPBYTz8/uaXoPi8/kkSUlJSZKko0ePyuv1KicnJzgmLi5OAwYMUElJiSSptLRUdXV1IWM8Ho8yMjKCYy4WCATk9/tDNgAA0H41OVCMMZoxY4buueceZWRkSJK8Xq8kyeVyhYx1uVzBY16vV7GxsercufMlx1yssLBQTqczuKWlpTV12gAAoA1ocqBMmzZNn3zyid58880GxxwOR8htY0yDfRe73Ji5c+fK5/MFt2PHjjV12gAAoA1oUqBMnz5d7777rrZv365u3boF97vdbklqcCWkoqIieFXF7XartrZWlZWVlxxzsbi4OCUmJoZsAACg/QorUIwxmjZtmjZs2KBt27YpPT095Hh6errcbreKioqC+2pra1VcXKzs7GxJUlZWlmJiYkLGlJeX6+DBg8ExAADg+hYdzuCpU6dq7dq1euedd5SQkBC8UuJ0OtWxY0c5HA7l5+dr/vz56tmzp3r27Kn58+erU6dOGj9+fHDspEmTNHPmTHXp0kVJSUmaNWuWMjMzNXTo0OY/QwAA0OaEFSgrVqyQJA0cODBk/6pVqzRx4kRJ0uzZs1VTU6MpU6aosrJSffv21datW5WQkBAcv2TJEkVHR2vs2LGqqanRkCFDtHr1akVFRV3b2QAAgHbhmj4HJVL4HBSgneFzUELwOShor1rtc1AAAABaAoECAACsQ6AAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOsQKAAAwDoECgAAsA6BAgAArEOgAAAA6xAoAADAOgQKAACwDoECAACsQ6AAAADrREd6AgBaVo85myM9hSsq6xDpGQCwDVdQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYJ+xA+eCDD3T//ffL4/HI4XBo06ZNIccnTpwoh8MRsvXr1y9kTCAQ0PTp05WcnKz4+HiNGjVKx48fv6YTAQAA7UfYgXL27Fn17t1by5Ytu+SYESNGqLy8PLht2bIl5Hh+fr42btyodevWadeuXaqurtbIkSNVX18f/hkAAIB2JzrcO+Tm5io3N/eyY+Li4uR2uxs95vP5tHLlSq1Zs0ZDhw6VJL3++utKS0vT+++/r+HDh4c7JQAA0M60yGtQduzYoZSUFPXq1UuPPfaYKioqgsdKS0tVV1ennJyc4D6Px6OMjAyVlJQ0+niBQEB+vz9kAwAA7VezB0pubq7eeOMNbdu2TS+88IL27t2rwYMHKxAISJK8Xq9iY2PVuXPnkPu5XC55vd5GH7OwsFBOpzO4paWlNfe0AQCARcL+Fc+VjBs3LvjfGRkZ6tOnj7p3767NmzdrzJgxl7yfMUYOh6PRY3PnztWMGTOCt/1+P5ECAEA71uJvM05NTVX37t115MgRSZLb7VZtba0qKytDxlVUVMjlcjX6GHFxcUpMTAzZAABA+9XigXL69GkdO3ZMqampkqSsrCzFxMSoqKgoOKa8vFwHDx5UdnZ2S08HAAC0AWH/iqe6ulpffvll8PbRo0e1f/9+JSUlKSkpSQUFBXrooYeUmpqqsrIyPfPMM0pOTtbo0aMlSU6nU5MmTdLMmTPVpUsXJSUladasWcrMzAy+qwcAAFzfwg6Uffv2adCgQcHbF14bMmHCBK1YsUIHDhzQa6+9pjNnzig1NVWDBg3S+vXrlZCQELzPkiVLFB0drbFjx6qmpkZDhgzR6tWrFRUV1QynBAAA2jqHMcZEehLh8vv9cjqd8vl8vB4FuIIeczZHegpXVNZhfKSnYJUe369ttecqW3Bfqz0XEM7Pb76LBwAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYJ2wv80YANB+tIUvk7wYX3B4feAKCgAAsA6BAgAArEOgAAAA6/AaFACwTFmH8ZGegiSpx/drIz0FXMe4ggIAAKxDoAAAAOsQKAAAwDoECgAAsA6BAgAArEOgAAAA6xAoAADAOgQKAACwDoECAACsQ6AAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOsQKAAAwDrRkZ4AgMgp6zA+0lMAgEZxBQUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYJ+xA+eCDD3T//ffL4/HI4XBo06ZNIceNMSooKJDH41HHjh01cOBAHTp0KGRMIBDQ9OnTlZycrPj4eI0aNUrHjx+/phMBAADtR9iBcvbsWfXu3VvLli1r9PiiRYu0ePFiLVu2THv37pXb7dawYcNUVVUVHJOfn6+NGzdq3bp12rVrl6qrqzVy5EjV19c3/UwAAEC7EfZ38eTm5io3N7fRY8YYLV26VM8++6zGjBkjSXr11Vflcrm0du1aTZ48WT6fTytXrtSaNWs0dOhQSdLrr7+utLQ0vf/++xo+fPg1nA4AAGgPmvU1KEePHpXX61VOTk5wX1xcnAYMGKCSkhJJUmlpqerq6kLGeDweZWRkBMcAAIDrW7N+m7HX65UkuVyukP0ul0vffPNNcExsbKw6d+7cYMyF+18sEAgoEAgEb/v9/uacNgAAsEyLvIvH4XCE3DbGNNh3scuNKSwslNPpDG5paWnNNlcAAGCfZg0Ut9stSQ2uhFRUVASvqrjdbtXW1qqysvKSYy42d+5c+Xy+4Hbs2LHmnDYAALBMswZKenq63G63ioqKgvtqa2tVXFys7OxsSVJWVpZiYmJCxpSXl+vgwYPBMReLi4tTYmJiyAYAANqvsF+DUl1drS+//DJ4++jRo9q/f7+SkpJ0ww03KD8/X/Pnz1fPnj3Vs2dPzZ8/X506ddL48eMlSU6nU5MmTdLMmTPVpUsXJSUladasWcrMzAy+qwcAAFzfwg6Uffv2adCgQcHbM2bMkCRNmDBBq1ev1uzZs1VTU6MpU6aosrJSffv21datW5WQkBC8z5IlSxQdHa2xY8eqpqZGQ4YM0erVqxUVFdUMpwQAANo6hzHGRHoS4fL7/XI6nfL5fPy6B7iCHnM2X/JYWYfxrTgTtDU9vl8b6Sk0qmzBfZGeApoonJ/ffBcPAACwDoECAACsQ6AAAADrECgAAMA6zfpR90B7d7kXnAIAmg9XUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHWiIz0BAICdyjqMj/QUJEk9vl8b6SkgApr9CkpBQYEcDkfI5na7g8eNMSooKJDH41HHjh01cOBAHTp0qLmnAQAA2rAW+RXP7bffrvLy8uB24MCB4LFFixZp8eLFWrZsmfbu3Su3261hw4apqqqqJaYCAADaoBYJlOjoaLnd7uDWtWtXSeevnixdulTPPvusxowZo4yMDL366qv673//q7VruYQHAADOa5FAOXLkiDwej9LT0/Xwww/r66+/liQdPXpUXq9XOTk5wbFxcXEaMGCASkpKLvl4gUBAfr8/ZAMAAO1Xs79Itm/fvnrttdfUq1cvnTx5Us8//7yys7N16NAheb1eSZLL5Qq5j8vl0jfffHPJxywsLNS8efOae6oAgDaox5zNkZ5C2MoW3BfpKbQ5zX4FJTc3Vw899JAyMzM1dOhQbd58/g/Sq6++GhzjcDhC7mOMabDvx+bOnSufzxfcjh071tzTBgAAFmnxz0GJj49XZmamjhw5Enw3z4UrKRdUVFQ0uKryY3FxcUpMTAzZAABA+9XigRIIBPTpp58qNTVV6enpcrvdKioqCh6vra1VcXGxsrOzW3oqAACgjWj216DMmjVL999/v2644QZVVFTo+eefl9/v14QJE+RwOJSfn6/58+erZ8+e6tmzp+bPn69OnTpp/Hg7PhAIAABEXrMHyvHjx/XrX/9ap06dUteuXdWvXz/t2bNH3bt3lyTNnj1bNTU1mjJliiorK9W3b19t3bpVCQkJzT0VAADQRjV7oKxbt+6yxx0OhwoKClRQUNDcTw0AANoJviwQAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWKfZP0kWwJWVdeC7pwDgcriCAgAArEOgAAAA6xAoAADAOgQKAACwDoECAACsQ6AAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOsQKAAAwDoECgAAsA6BAgAArMO3GSNieszZHOkpAAAsxRUUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANbhc1AAAFYr6zA+0lOQJPX4fm2kp3BdIVAAAGhhbfGDKcsW3BfR5ydQcF2x5V9iAIDL4zUoAADAOgQKAACwDoECAACsQ6AAAADrECgAAMA6BAoAALAObzNuJ9rie+wBALgUAgUAgKtgy+coXS+faEugoFXY8j82AKBtiOhrUJYvX6709HR16NBBWVlZ2rlzZySnAwAALBGxQFm/fr3y8/P17LPP6uOPP9avfvUr5ebm6ttvv43UlAAAgCUiFiiLFy/WpEmT9Pvf/1633nqrli5dqrS0NK1YsSJSUwIAAJaIyGtQamtrVVpaqjlz5oTsz8nJUUlJSYPxgUBAgUAgeNvn80mS/H5/i8wv47l/tsjjXs/8DhPpKQBAu3Au8N9WeZ6W+Bl74TGNufLPhIgEyqlTp1RfXy+XyxWy3+Vyyev1NhhfWFioefPmNdiflpbWYnNE83JGegIA0G6MbZVncS5tuceuqqqS03n5nwwRfRePw+EIuW2MabBPkubOnasZM2YEb587d07/+c9/1KVLl0bHtyd+v19paWk6duyYEhMTIz0da7FOV4d1ujqs09Vhna6MNQpljFFVVZU8Hs8Vx0YkUJKTkxUVFdXgaklFRUWDqyqSFBcXp7i4uJB9P/vZz1pyitZJTEzkD/dVYJ2uDut0dVinq8M6XRlr9D9XunJyQUReJBsbG6usrCwVFRWF7C8qKlJ2dnYkpgQAACwSsV/xzJgxQ3l5eerTp4/69++vl19+Wd9++60ef/zxSE0JAABYImKBMm7cOJ0+fVr/93//p/LycmVkZGjLli3q3r17pKZkpbi4OD333HMNfsWFUKzT1WGdrg7rdHVYpytjjZrOYa7mvT4AAACtKKIfdQ8AANAYAgUAAFiHQAEAANYhUAAAgHUIFAtVVlYqLy9PTqdTTqdTeXl5OnPmzBXv9+mnn2rUqFFyOp1KSEhQv3792vW3Qzd1nS6YPHmyHA6Hli5d2mJzjLRw16iurk5PP/20MjMzFR8fL4/Ho9/+9rc6ceJE6026FSxfvlzp6enq0KGDsrKytHPnzsuOLy4uVlZWljp06KAbb7xRf/3rX1tpppEVzjpt2LBBw4YNU9euXZWYmKj+/fvrn/+8Pr7XLNw/Txf861//UnR0tO68886WnWBbZWCdESNGmIyMDFNSUmJKSkpMRkaGGTly5GXv8+WXX5qkpCTz1FNPmY8++sh89dVX5u9//7s5efJkK8269TVlnS7YuHGj6d27t/F4PGbJkiUtO9EICneNzpw5Y4YOHWrWr19vPvvsM7N7927Tt29fk5WV1Yqzblnr1q0zMTEx5pVXXjGHDx82TzzxhImPjzfffPNNo+O//vpr06lTJ/PEE0+Yw4cPm1deecXExMSYt99+u5Vn3rrCXacnnnjCLFy40Hz44Yfmiy++MHPnzjUxMTHmo48+auWZt65w1+mCM2fOmBtvvNHk5OSY3r17t85k2xgCxTKHDx82ksyePXuC+3bv3m0kmc8+++yS9xs3bpz5zW9+0xpTtEJT18kYY44fP25+/vOfm4MHD5ru3bu320C5ljX6sQ8//NBIuuJfuG3F3XffbR5//PGQfbfccouZM2dOo+Nnz55tbrnllpB9kydPNv369WuxOdog3HVqzG233WbmzZvX3FOzSlPXady4ceZPf/qTee655wiUS+BXPJbZvXu3nE6n+vbtG9zXr18/OZ1OlZSUNHqfc+fOafPmzerVq5eGDx+ulJQU9e3bV5s2bWqlWbe+pqyTdH6t8vLy9NRTT+n2229vjalGTFPX6GI+n08Oh6NdfP9VbW2tSktLlZOTE7I/Jyfnkmuye/fuBuOHDx+uffv2qa6ursXmGklNWaeLnTt3TlVVVUpKSmqJKVqhqeu0atUqffXVV3ruuedaeoptGoFiGa/Xq5SUlAb7U1JSGny54gUVFRWqrq7WggULNGLECG3dulWjR4/WmDFjVFxc3NJTjoimrJMkLVy4UNHR0frjH//YktOzQlPX6Me+//57zZkzR+PHj28XX3R26tQp1dfXN/hSUpfLdck18Xq9jY7/4YcfdOrUqRabayQ1ZZ0u9sILL+js2bMaO3ZsS0zRCk1ZpyNHjmjOnDl64403FB0dsQ9zbxMIlFZSUFAgh8Nx2W3fvn2SJIfD0eD+xphG90vn/6UiSQ888ICefPJJ3XnnnZozZ45GjhzZ5l7M15LrVFpaqhdffFGrV6++5Ji2oCXX6Mfq6ur08MMP69y5c1q+fHmzn0ckXXz+V1qTxsY3tr+9CXedLnjzzTdVUFCg9evXNxrJ7c3VrlN9fb3Gjx+vefPmqVevXq01vTaLfGsl06ZN08MPP3zZMT169NAnn3yikydPNjj23XffNaj0C5KTkxUdHa3bbrstZP+tt96qXbt2NX3SEdCS67Rz505VVFTohhtuCO6rr6/XzJkztXTpUpWVlV3T3FtLS67RBXV1dRo7dqyOHj2qbdu2tYurJ9L5/1eioqIa/Ou2oqLikmvidrsbHR8dHa0uXbq02FwjqSnrdMH69es1adIkvfXWWxo6dGhLTjPiwl2nqqoq7du3Tx9//LGmTZsm6fw/MI0xio6O1tatWzV48OBWmXtbQKC0kuTkZCUnJ19xXP/+/eXz+fThhx/q7rvvliT9+9//ls/nU3Z2dqP3iY2N1V133aXPP/88ZP8XX3zR5r58sSXXKS8vr8FfmMOHD1deXp5+97vfXfvkW0lLrpH0vzg5cuSItm/f3q5+CMfGxiorK0tFRUUaPXp0cH9RUZEeeOCBRu/Tv39//e1vfwvZt3XrVvXp00cxMTEtOt9Iaco6SeevnDz66KN68803dd9997XGVCMq3HVKTEzUgQMHQvYtX75c27Zt09tvv6309PQWn3ObErnX5+JSRowYYe644w6ze/dus3v3bpOZmdngraE333yz2bBhQ/D2hg0bTExMjHn55ZfNkSNHzF/+8hcTFRVldu7c2drTbzVNWaeLted38RgT/hrV1dWZUaNGmW7dupn9+/eb8vLy4BYIBCJxCs3uwttCV65caQ4fPmzy8/NNfHy8KSsrM8YYM2fOHJOXlxccf+Ftxk8++aQ5fPiwWbly5XX1NuOrXae1a9ea6Oho89JLL4X8uTlz5kykTqFVhLtOF+NdPJdGoFjo9OnT5pFHHjEJCQkmISHBPPLII6aysjJkjCSzatWqkH0rV640N910k+nQoYPp3bu32bRpU+tNOgKauk4/1t4DJdw1Onr0qJHU6LZ9+/ZWn39Leemll0z37t1NbGys+eUvf2mKi4uDxyZMmGAGDBgQMn7Hjh3mF7/4hYmNjTU9evQwK1asaOUZR0Y46zRgwIBG/9xMmDCh9SfeysL98/RjBMqlOYz5/6/2AgAAsATv4gEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFjn/wEe95Ly+GbpewAAAABJRU5ErkJggg==",
"text/plain": [
"<Figure size 640x480 with 1 Axes>"
]
@@ -159,57 +167,45 @@
],
"source": [
"import scipy\n",
- "\n",
- "mus = np.linspace(0,50,1000)\n",
- "plt.plot(mus, scipy.stats.poisson.pmf(18,mus))\n",
- "plt.xlabel(r\"$\\mu$\")\n",
- "scipy.integrate.quad(lambda x: scipy.stats.poisson.pmf(18,x), 0, 100)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "71c0df91-c1e0-4fbf-a94b-fd88f839ef07",
- "metadata": {
- "slideshow": {
- "slide_type": "subslide"
- },
- "tags": []
- },
- "source": [
- "Also: $P(\\mu_- \\le \\mu \\le\\mu_+) = \\int_{\\mu_-}^{\\mu_+} P(k, \\mu) d\\mu$"
+ "import numpy as np\n",
+ "import matplotlib.pyplot as plt\n",
+ "#xs = ys\n",
+ "mu1 = np.mean(xs, axis=1)\n",
+ "mu2 = (np.max(xs, axis=1) + np.min(xs, axis=1))/2\n",
+ "plt.hist(mu1-1)\n",
+ "plt.hist(mu2-1)"
]
},
{
"cell_type": "markdown",
- "id": "d278b736-1f45-4b6c-8e1a-c851bb618aeb",
+ "id": "5b83fd0c",
"metadata": {
"slideshow": {
- "slide_type": ""
+ "slide_type": "slide"
},
"tags": []
},
"source": [
- "Test der naiven Antwort:"
+ "### Estimator for the variance\n",
+ "\n",
+ "Estimator for the variance $\\sigma^2$:\n",
+ "\n",
+ "Let's repeat the tests from before with the sample variance \n",
+ "$$\\hat \\sigma^2 = V[x] = \\frac{1}{n}\\sum_{i=1}^n (x_i - <x>)^2$$\n",
+ "as an estimator."
]
},
{
"cell_type": "code",
"execution_count": null,
- "id": "7c461ddf-0cb8-49b0-adc6-a7ae6f11212b",
- "metadata": {
- "slideshow": {
- "slide_type": ""
- },
- "tags": []
- },
+ "id": "b61361d1",
+ "metadata": {},
"outputs": [],
- "source": [
- "scipy.integrate.quad(lambda x: scipy.stats.poisson.pmf(18,x), 18-np.sqrt(18), 18+np.sqrt(18))"
- ]
+ "source": []
},
{
"cell_type": "markdown",
- "id": "88f176b2-deb2-4212-b4ba-d86b7964e545",
+ "id": "4dfd3266-a5a2-4db7-a96b-ef21118bc85e",
"metadata": {
"slideshow": {
"slide_type": "slide"
@@ -217,220 +213,145 @@
"tags": []
},
"source": [
- "## Example from Publication:\n",
- "\"Search for flavor changing neutral currents in decays of top quarks\" (D0)\n",
+ "### Estimator of the variance\n",
"\n",
- "Physics Letters B 701 (2011), pp. 313-320\n",
- "[https://arxiv.org/abs/1103.4574]\n",
+ "Test:\n",
"\n",
- "![wrong](./figures/12/ht_bq_1jets_comb_5pc.png \"Example\")"
+ "- consistent:\n",
+ " $\\lim_{n\\to \\infty } \\hat {\\sigma^2} = \\lim_{n\\to \\infty } \\frac{n-1}{n}\\sigma^2 =\\sigma^2$\n",
+ "\n",
+ "- unbiased:\n",
+ " $$E[\\hat {\\sigma^2}]= E[ \\frac{1}{n}\\sum_{i=1}^n (x_i^2 - \\bar x^2)] = \\frac{1}{n} \\sum_{i=1}^n E[x_i^2-\\bar x^2] = E[x^2] - E[\\bar x^2]$$\n",
+ " $$E[\\hat{\\sigma^2}]= E[x^2] - E[x]^2 + E[x]^2 - E[\\bar x^2] = E[x^2] - E[x]^2 - (E[\\bar x^2] - E[\\bar x]^2)$$\n",
+ " $$E[\\hat {\\sigma^2}] = V(x) - V(\\bar x) = \\sigma^2 - \\frac{\\sigma^2}{n} = \\frac{n-1}{n}\\sigma^2 \\ne \\sigma^2$$\n",
+ " \n",
+ "Bias corrected estimator:\n",
+ "$$\\hat \\sigma^2 = \\frac{n}{n-1}V[x] = \\frac{n}{n-1} \\frac{1}{n}\\sum_{i=1}^n (x_i - <x>)^2 = \\frac{1}{n-1}\\sum_{i=1}^n (x_i - <x>)^2$$\n"
]
},
{
"cell_type": "markdown",
- "id": "5135b8e6-711b-440e-8507-f049d017c137",
+ "id": "cb43d9dd",
"metadata": {
"slideshow": {
- "slide_type": ""
- },
- "tags": []
+ "slide_type": "slide"
+ }
},
"source": [
- "Hm, ok!\n",
- "\n",
- "Nächstes Bin: 1-mal fielen sieben Tore. \n",
- "\n",
- "Naives Konfidenzintzervall: $1\\pm 1$\n",
+ "### Test with Monte Carlo\n",
"\n",
- "Test:"
+ "Repeat the MC test with the new estimator:\n",
+ "$$\\hat \\sigma^2 = \\frac{1}{n-1}\\sum_{i=1}^n (x_i - <x>)^2$$\n",
+ "This is `np.var(xs, ddof=1)` in *numpy*"
]
},
{
"cell_type": "code",
"execution_count": null,
- "id": "e6102fec-1275-439a-b6b4-8cb1ee4960fc",
- "metadata": {
- "slideshow": {
- "slide_type": ""
- },
- "tags": []
- },
+ "id": "c7d6b2ac",
+ "metadata": {},
"outputs": [],
- "source": [
- "mus = np.linspace(0,10,1000)\n",
- "plt.plot(mus, scipy.stats.poisson.pmf(1,mus))\n",
- "plt.xlabel(r\"$\\mu$\")\n",
- "scipy.integrate.quad(lambda x: scipy.stats.poisson.pmf(1,x), 1-np.sqrt(1), 1+np.sqrt(1))"
- ]
+ "source": []
},
{
"cell_type": "markdown",
- "id": "3139be1a-8a4d-4ae9-9fa7-aa99d7efffe6",
+ "id": "7c7c6ecd-1611-4787-8985-da0881e90d9f",
"metadata": {
- "jupyter": {
- "source_hidden": true
- },
"slideshow": {
- "slide_type": ""
+ "slide_type": "slide"
},
"tags": []
},
"source": [
- "zu klein!\n",
- "\n",
- "Suche Intervall: $[\\mu_-, \\mu_+]$\n",
+ "# Methode der kleinsten Quadrate\n",
"\n",
- "hier:\n",
- "$\\mu_- = 0$"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "630c83e0-2be4-41dd-a124-2ac672a2c507",
- "metadata": {},
- "outputs": [],
- "source": [
- "def intervall(mu_plus):\n",
- " return scipy.integrate.quad(lambda x: scipy.stats.poisson.pmf(1,x), 0, mu_plus)[0]\n",
+ "## Herleitung\n",
"\n",
- "mus = np.linspace(0,10,100)\n",
- "plt.plot(mus, np.vectorize(intervall)(mus))\n",
- "plt.grid()\n",
- "plt.xlabel(r\"$\\mu$\")\n",
+ "### Methode der kleinsten Quadrate \n",
"\n",
- "scipy.optimize.brentq(lambda x: intervall(x)-0.6827, 0,10)"
+ "$y(x) = mx + a$: Finde $\\hat m$ und $\\hat a$!"
]
},
{
"cell_type": "markdown",
- "id": "1c4ad7b4-bbdb-4f9d-a795-0e2dae39a181",
+ "id": "73e6e46a-4e41-4e66-857b-2946946498d1",
"metadata": {
"slideshow": {
- "slide_type": "slide"
+ "slide_type": "skip"
},
"tags": []
},
"source": [
- "Korrektes Intervall: $[0; 2{,}36]$\n",
- "\n",
- "\n",
- "## Was ist denn nun das richtige Intervall für das 6-Tore-Bin?\n",
- "\n",
- "Welches?\n",
- "\n",
- "\n",
- "z.B. symmetrisch in der Konfidenz:\n"
+ "<img src=\"./figures/11/line.png\" style=\"width:90.0%\" alt=\"image\" />"
]
},
{
"cell_type": "code",
- "execution_count": null,
- "id": "5fffd76e-14da-4854-8af0-5262e0c3c9a3",
+ "execution_count": 1,
+ "id": "55bf0cad-8a50-4831-8a37-daf9bcb39b71",
"metadata": {
"slideshow": {
"slide_type": ""
},
"tags": []
},
- "outputs": [],
+ "outputs": [
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjMAAAG0CAYAAADHD6Y/AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA3VElEQVR4nO3deXgV5eH28e9JQhKCSQDZCwq4gIKoBVQURQVxQRGlra+11uXXVhBUsIugrYq1oG1FWhco1lqtdWlFEEVlUQF3AVHBBTcQFBBRTJBAQpJ5/5gaRNkSksyZk+/nus7V80xmknuuocntM8tJBEEQIEmSFFNpUQeQJEnaHZYZSZIUa5YZSZIUa5YZSZIUa5YZSZIUa5YZSZIUa5YZSZIUa5YZSZIUa5YZSZIUa5YZSZIUa5GWmblz53LaaafRqlUrEokEU6ZM2errQRBw7bXX0qpVK+rXr8+xxx7Lm2++GU1YSZKUlDKi/OEbNmzg4IMP5oILLmDgwIHf+fof//hHxo4dyz//+U/2339/rr/+ek444QSWLFlCbm7uLv2M8vJyVq5cSW5uLolEorp3QZIk1YAgCFi/fj2tWrUiLW0ncy9BkgCCyZMnV4zLy8uDFi1aBDfccEPFsk2bNgX5+fnBhAkTdvn7rlixIgB8+fLly5cvXzF8rVixYqd/6yOdmdmRpUuXsnr1avr27VuxLCsri169evHCCy9w0UUXbXO74uJiiouLK8bB/z4UfMWKFeTl5dVsaEmSVC0KCwtp06bNLp2JSdoys3r1agCaN2++1fLmzZvz0UcfbXe7MWPGMGrUqO8sz8vLs8xIkhQzu3KJSNLfzfTtnQiCYIc7NnLkSAoKCipeK1asqOmIkiQpQkk7M9OiRQsgnKFp2bJlxfI1a9Z8Z7bmm7KyssjKyqrxfJIkKTkk7cxMu3btaNGiBTNnzqxYVlJSwpw5czjyyCMjTCZJkpJJpDMzX331Fe+//37FeOnSpbz22ms0btyYvfbai2HDhjF69Gj2228/9ttvP0aPHk1OTg4//vGPI0wtSZKSSaRlZv78+Rx33HEV48svvxyA8847j3/+85/85je/YePGjVx88cWsW7eOww8/nBkzZuzyM2YkSVLqSwRf37ucogoLC8nPz6egoMC7mSRJionK/P1O2mtmJEmSdoVlRpIkxZplRpIkxZplRpIkxZplRpIkxZplRpIkxZplRpIkxZplRpIkxZplRpIkVVpRSSltR0yj7YhpFJWURprFMiNJkmLNMiNJkmLNMiNJkmLNMiNJkmLNMiNJkmLNMiNJkmLNMiNJkmLNMiNJkmLNMiNJkmLNMiNJkmLNMiNJkmLNMiNJknbL6oJNkf58y4wkSaq0SQs+rnjfZ+wcHpy3PLIslhlJklQpqwo2cs3UNyvG5QFc+fBiVhVsjCSPZUaSJFXK0rUbKA+2XlYWBCxbWxRJHsuMJEmqlHZNGpCW2HpZeiJB2yY5keSxzEiSpEppmV+fUf07VYzTEjD6zM60zK8fSR7LjCRJqrSBXVtXvJ91eS/O6r5XZFksM5Ikabe0yM+O9OdbZiRJUqxZZiRJUqxZZiRJUpV0SCwnk81Rx7DMSJKkSirbTMZzf+axzKsYljEp6jRkRB1AkiTFzFOjyHzhFkjAXolPISiPNI4zM5IkqXJ6DKW8UXsuLRnC0M2XQiLaOmGZkSRJO7b2fXh27JZxbgs2XfQiU8uPAhLb3ay2eJpJkiRtW3k5vPI3mDUKSjdC0w7QsV/4tbTkqRDJk0SSJCWPdcvgkaGw7Nlw3P5YaNElykTbZZmRJElbBAHM/wfM+B1s3gD1GkDf66Db/0Ei+lNK22KZkSRJW0wZDK/fH77f60gYcDs0bved1XIyM1h2Q79aDrdtXgAsSZK2OHAAZGTDiWPg/GnbLDLJxpkZSZLqsvWr4bMl0L5XOO5wElz2BuQ2jzZXJTgzI0lSXRQEsOghuO1w+M+5ULhyy9diVGTAmRlJkuqeDWvhseHw9tRw3PJg2Lwx2ky7wTIjSVJd8tYj8NjlULQ2fFbMMb+Boy+H9HpRJ6syy4wkSXVBeTlMvggW/SccN+sEZ4wPZ2VizjIjSVJdkJYGOY3Dz1HqORx6XQEZWVGnqhaWGUmSUtWmAijZAHmtwnHva6DLWfC970ebq5p5N5MkSanog6fh9iPhof8LTzEBZOakXJEBZ2YkSUotxV/BzKth/p3hOL0erF8J+a2jzVWDLDOSJKWKZc/BlIvhy4/C8WG/gD7XQmaDSGPVNMuMJElxt3kTzLoWXh4fjvPbwOm3bXmqb4qzzEiSFHeJNFj2bPj++z+Fvn+A7LxoM9Uiy4wkSXG0eROkpYfXxGRkwhkTws9Z2u+EqJPVOu9mkiQpblYuhIm94NmxW5a1OKhOFhmwzEiSFB+lJfDMaLijN3z2Drx6N5QURZ0qcp5mkiQpDlYvhimDYPWicNzpDDjlpvDZMXWcZUaSpGRWVgrP3wyzb4TyzVC/MfS7CTqfGXWypGGZkSSpFhWVlHLg1dMBeOu6E8nJ3Mmf4oIVMPemsMh06Aen3gy5zWshaXxYZiRJSjZBAIlE+L5xOzj5BsjIDj9X6evlquAFwJIkJZPPP4B/9oPlL21Z1vV8OPj/WWS2wzIjSVIyKC+HlyfChJ7w0fPw+K/DGRrtlKeZJEmK2pfL4ZEhsHRuOG57dPhxBM7E7BLLjCRJUQkCWPBPmH4VlHwF9XLghOug2/9BmidPdpVlRpKkiKR9+BQ8elk4aHMEDLgd9twn2lAxZJmRJCki5e17wwH9oc3hcMTg8LOWVGlJPYdVWlrKb3/7W9q1a0f9+vVp37491113HeXl5VFHkySp8tZ/Subjw8hjQzhOJOBH98CRQy0yuyGpZ2ZuvPFGJkyYwN13302nTp2YP38+F1xwAfn5+Vx22WVRx5Mkadctfhim/ZKMjV9wVcYKrij9Rbjci3x3W1KXmRdffJHTTz+dfv36AdC2bVvuv/9+5s+fH3EySZJ20YbP4fFfwpuTAShvfhB3LT8p4lCpJalPM/Xs2ZOnnnqKd999F4DXX3+d5557jlNOOSXiZJIk7YJ3psHth4dFJpEOva5g0/kzeCfYK+pkKSWpZ2auuOIKCgoK6NixI+np6ZSVlfGHP/yBs88+e7vbFBcXU1xcXDEuLCysjaiSJG1t3p0w7fLwfdMD4Izx0OpQKCmtWGV1wSbaN90jooCpI6lnZh588EHuvfde7rvvPl599VXuvvtu/vznP3P33Xdvd5sxY8aQn59f8WrTpk0tJpYk6X86nQG5reCoYXDRnLDIAJMWfFyxSp+xc3hw3vKIAqaORBAk77OS27Rpw4gRIxgyZEjFsuuvv557772Xd955Z5vbbGtmpk2bNhQUFJCXl1fjmSVJddSmQlj0n/CBd19f1Fv8FWRtmXlZVbCRo254mvJv/OVNTyR4bsRxtMyvX8uBk1thYSH5+fm79Pc7qU8zFRUVkfatJyCmp6fv8NbsrKwssrKyajqaJKmWFZWUcuDV0wF467oTyclMoj9hH84JP46gYAVk5UOXH4bLs7Y+hbR07YatigxAWRCwbG2RZWY3JNG/hO867bTT+MMf/sBee+1Fp06dWLhwIWPHjuXCCy+MOpokSVCyAWZdC69MDMcN94a8VttdvV2TBqQl+M7MTNsmOTWbM8UldZm55ZZb+N3vfsfFF1/MmjVraNWqFRdddBFXX3111NEkSXXdRy/ClMGwbmk47nYhnPD778zGfFPL/PqM6t+J3z3yJgBpCRh9ZmdnZXZTUpeZ3Nxcxo0bx7hx46KOIknSFs+NC2dkCCDve9D/Fti39y5tOrBr64oyM+vyXt7NVA2SusxIkpSUWh0KBHDIT+Ck0ZCdX6Vv0yI/u3pz1VGWGUmSdqa0GFYvhtZdw3H7XnDxS9DsgGhzCUjy58xIkhS5Va/DxOPg7tPgi6VblltkkoZlRpKkbSnbDLNvhDuOhzVvQr36UPDxzrdTrfM0kyRJ37bmbZh8UTgrA3BAfzj1ZmjQJNpc2ibLjCRJ3/TCLfDUdVBWAtkNod9N0Hnglqf6KulYZiRJ+qaiz8Mis/9JcNpfILdF1Im0E5YZSVLdVl4OG7/Ycgrp2JHQokv4QZHOxsSCFwBLkmJndcGm6vlGXyyFu0+Ff/8AykrDZRlZ0PlMi0yMWGYkSbEwacGWO4n6jJ3Dg/OWV/2bBQHM+zuMPwo+eh4+exc+XVQNKRUFy4wkKemtKtjINVPfrBiXB3Dlw4tZVbCx8t/syxXwrwEw7ZeweQPs3RMGP/+/p/oqjrxmRpKU9Jau3bDVJ00DlAUBy9YW7fqHNAYBvPZveHIkFBdCRjb0uRYOuwjS/G/7OLPMSJKSXrsmDUhLsFWhSU8kaNskZ9e/SVAO8+8Ki0zr7jBgAjTZt/rD7kROZgbLbuhX6z83lVlFJUlJr2V+fUb171QxTkvA6DM773xWJgigvOx/G6XDgPHQZxRcOD2SIqOaYZmRJMXCwK6tK97PurwXZ3Xfa8cbfPUZ/OdceGrUlmVN94eew8Jio5ThaSZJUuy0yM/e8QpvPQKPDQ8fgJeeBUdc7MPvUphlRpKUOoq+gMd/DYsfCsfNOsEZ4y0yKc4yI0lKDUuehEcvha8+hUQa9BwOva4IH4KnlGaZkaQUV1RSyoFXTwfgretOJCczBX/1F30Bk34GJeuhyf7hnUqtu0adSrUkBf9FS5LqnJzGcNIY+OwdOP63UG8Xnz2jlGCZkSTFT/F6mHEdHHAa7HN8uOz750abSZHx1mxJUqwcnnib7L8fA/P/AVMvhdLiqCMpYs7MSJLiYXMR12TczQUZ06EAyN8LTr/VC3xlmZEkxcDyl8meMpgLMj4AoPSQn5Jx0h8gOy/iYEoGnmaSJCW3T9+Eu04i7YsPWBU05rySKyg55WaLjCo4MyNJSm7NO8EB/SlNz+LEeb0pZI+oEynJODMjSUoupSUw98+wYe2WZQP/Tslpt1tktE3OzEiSksfqxTBlEKxeBKtehx/dA4kEpNcjJx2W3dAv6oRKQpYZSVL0ykrh+XEw+wYo3wz1G0OnM8IiI+2EZUaSFK3PlsDkQbDy1XDcoR+cNg72aBZpLMWHZUaS6pDVBZto3zSJrjt5fxbc/2MoK4bsfDj5j9DlLGdkVCleACxJKW7Sgo8r3vcZO4cH5y2PMM23tO4ODZrCvn3g4pfg4P9nkVGlWWYkKYWtKtjINVPfrBiXB3Dlw4tZVbAxmkDl5fD2oxAE4Tg7H342E855CPJaRZNJsWeZkaQUtnTtBsqDrZeVBQHL1hbVfph1H8E9/eHBn8Crd29ZntfK2RjtFq+ZkaQU1q5JA9ISbFVo0hMJ2jbJqb0QQRCWl+lXQclXUC8HsLyo+jgzI0kprGV+fUb171QxTkvA6DM70zK/fu0EKPgE7h0Ij14WFpm9esCg56DrebXz81UnWGYkKcUN7Nq64v2sy3txVve9aucHvzMNbu8BHzwF6VnQ93o4fxrsuU/t/HzVGZ5mkqQ6pEV+du39sAbNoGQ9fK8rDJgATfevvZ+tOsUyI0mqPp9/sGXmpU13+OkjsNeRkO6fG9UcTzNJknbfhs/hP+eFp5U+W7JlebtjLDKqcf4LkyTtnrcfg8eGwYbPIC0DVrwCTTtEnUp1iGVGklQ1G9fBE1fAGw+G46YHwBnjodWh0eZSnWOZkSRV3nuzYOpQWL8KEmlw5KVw3JWQkRV1MtVBlhlJUuWtei0sMnvuG96p1KZ71IlUh1lmJEm7ZvNGqPe/h+0dNSychen2f5BZi08TlrbBu5kkSTtWsgGm/QruOB5Ki8Nl6Rlw5CUWGSUFZ2YkSdv30YswZTCsWxqO35sJB5wabSbpWywzkqTv2rwRnr4eXrwNCCCvNZx+K+xzXNTJpO+wzEhSisvJzGDZDf12fYOP58PkQfD5e+H40J/AiaMhO79mAkq7yTIjSdra078Pi8weLaD/X2H/E6NOJO2QZUaStLXT/gpz/wQnXAc5jaNOI+2UZUaS6rKyzfDsTbDxSzj5hnBZo73D62OkmLDMSFJd9elbMGUQrHo9HB/yY2jZJdpMUhVYZiSprikrhRf+CrPHQFkJZDeEfjdBi4OiTiZViWVGkuqSte+Fdyp9Mj8c738SnPYXyG0RbS5pN1hmJKmuKC2Bu/vD+pWQlQcn3RCeWkokok4m7RbLjCTVFRmZ0Pf38Nq/of8tkN866kRStbDMSFKqKi+H+XeGpaXDyeGyzgPDl7MxSiGWGUlKRV+ugEeGwNI50KApDHklfGaMJUYpyDIjSakkCGDhvfDkSChZDxn14ehfhXcsSSnKMiNJqaJwFTx6Kbw3Ixy3PgwGjIcm+0abS6phlhlJSgXrP4Xbj4BNX0J6Jhz/W+gxFNLSo04m1TjLjCSlgtzm0PFUWPMmDJgAzTpGnUiqNZYZSYqrtx6B1t0hr1U4PuWP4axMer1oc0m1LC3qAJKkSir6Ah66EP7zU5h6aXjRL0BmA4uM6qSkLzOffPIJP/nJT9hzzz3JycnhkEMOYcGCBVHHkqRoLHkyvDZm8SRIpEOrQ6C8LOpUUqSS+jTTunXrOOqoozjuuON44oknaNasGR988AENGzaMOpok1a6NX4a3W79+Xzhu0gHOGA/f6xppLCkZJHWZufHGG2nTpg133XVXxbK2bdtGF0iSorDmbbh3IBR+AiTgyEvguKugXnbUyaSkkNSnmaZOnUq3bt344Q9/SLNmzTj00EO54447oo4lSbWr4d6QkQ2N28OF08PPV7LISBWSusx8+OGHjB8/nv3224/p06czaNAgLr30Uu65557tblNcXExhYeFWL0mKnU8WhJ+tBJCZAz/+Dwx6DvY6PNpcUhJKBMHXl8Enn8zMTLp168YLL7xQsezSSy9l3rx5vPjii9vc5tprr2XUqFHfWV5QUEBeXl6NZZWkalFSBLOuhVf+BieOhh5Dok4kRaKwsJD8/Pxd+vud1DMzLVu25MADD9xq2QEHHMDy5cu3u83IkSMpKCioeK1YsaKmY0pS9Vj+MkzoGRYZCD8sUtJOJfUFwEcddRRLlizZatm7777L3nvvvd1tsrKyyMrKquloklR9Nm+CZ66HF24FAshtBaffAvv2iTqZFAtJXWaGDx/OkUceyejRo/nRj37EK6+8wsSJE5k4cWLU0SSpeqx6Ayb9DNb+7z/cDj4bTroB6jeMNJYUJ0ldZrp3787kyZMZOXIk1113He3atWPcuHGcc845UUeTpOqRSMAXH0KDZnDaX6DjKVEnkmInqS8Arg6VuYBIkmpF0ReQ03jL+O3HYO8jt14m1XEpcwGwJKWUslKY8ye4uROsfG3L8gNOtchIuyGpTzNJUspY8w5MGQQrF4bjxZPCz1WStNssM5JUk8rL4MXb4OnroawYsvPh5D9Blx9FnUxKGZYZSdqOopJSDrx6OgBvXXciOZmV/JX5+QcwZTCseDkc73sC9P8r5LWq5qRS3WaZkaSa8t7MsMhk5sJJo+HQc8O7lyRVK8uMJFWn8nJI+9+9FYf9AtavhO4/g4Z7RZtLSmHezSRJ1SEIYP5dMPEYKNkQLktLgxOus8hINcwyI0m7q+ATuHcgPDYMVi+CBf+MOpFUp3iaSZKqKgjg9fvhiRFQXAAZ2dD7ajh8cNTJpDrFMiNJVbH+U3j0Mnj3iXD8va4wYAI03T/aXFIdZJmRpKqYcVVYZNIz4diRcOSlkO6vVCkK/j9PkqrihN/DV2vCT7hufmDUaaQ6zQuAJWkXpC95DKZftWVBXks4b6pFRkoCzsxI0g7k8RWj6t1N1qTnwwX7HA/79o42lKStWGYkaTvSPpjFjKwraJFYR5BII3HUMGjbM+pYkr7FMiNJ37apEKZfSfbCfxHQmIdKe9LtrCtpe9BRUSeTtA2WGUn6piCAfw2ATxbwQOmxjCz9GQFppN33JWPOXM5Z3X2ar5RsvABYkr4pkYCjf8WqvIO5suwXBP/7NVkewJUPL2ZVwcaIA0r6NsuMJH30Arwzbcu44yksPfVByoOtVysLApatLardbJJ2ytNMkuquzRvhqd/DS7dDdh5c/HJ4yzXQrnlD0hJsVWjSEwnaNsmJKKyk7XFmRlLd9PF8mHA0vHQbEMAB/SGzQcWXW+bXZ1T/ThXjtASMPrMzLfPrRxBW0o5Uusycf/75zJ07tyaySFLNKy2GWdfCnSfA5+/BHi3gx/+F028NZ2e+YWDX1hXvZ13ey4t/pSRV6TKzfv16+vbty3777cfo0aP55JNPaiKXJFW/zRvhjuPhuZshKIcuZ8GQl2D/vjvdtEV+di0ElFQVlS4zkyZN4pNPPmHo0KH897//pW3btpx88sk89NBDbN68uSYySlL1qFcf9j4ScprAWffCmROhfqOoU0naTVW6ZmbPPffksssuY+HChbzyyivsu+++nHvuubRq1Yrhw4fz3nvvVXdOSaqaT9+Cdcu2jPtcC0NehgNOiyqRpGq2WxcAr1q1ihkzZjBjxgzS09M55ZRTePPNNznwwAO5+eabqyujJFVeWSk8exNM7AVThkB5ebg8swE0aBJtNknVqtK3Zm/evJmpU6dy1113MWPGDLp06cLw4cM555xzyM3NBeCBBx5g8ODBDB8+vNoDS9JOffYuTBkMn8wPx9l5sHkDZOVGm0tSjah0mWnZsiXl5eWcffbZvPLKKxxyyCHfWefEE0+kYcOG1RBPkiqhvAxeGg9P/x5KN0FWPpx8Axx8dvhkX0kpqdJl5uabb+aHP/wh2dnbv7K/UaNGLF26dLeCSVKlfLUG/nMeLH8hHO9zPPS/BfJb73g7SbFX6TJz7rnn1kQOSdo92fmwqQDqNYATr4euFzgbI9URfpyBpPgq+AT2aAbp9SAjC37wD6iXDY3aRp1MUi3y4wwkxU8QwKv3wG2Hw3Pjtixv1tEiI9VBlhlJ8VK4Cu77EUy9BErWw4ezwwt/JdVZnmaSFA9BAG/8B574dXhtTHoWHH8V9BgKaelRp5MUIcuMpOT31Wfw2DB457Fw3OpQGDAhPK1Ug3IyM1h2Q78a/RmSdp9lRlLy2/QlvD8L0upBryug53BI99eXpJC/DSQlp9Li8A4lgCb7Qf9bw5mYFgdFm0tS0vECYEnVrqiklLYjptF2xDSKSkor/w3eeRz+cjAsf3nLsi4/tMhI2ibLjKTksfFLmDwYHjgb1q+C58dFnUhSDHiaSVJyeH8WPHIJrF8JJODIS+C4q6JOJSkGLDOSolW8HqZfBa/eHY4bt4cB42GvI6LNJSk2LDOSovXOtC1F5vBB0PsayMyJNpOkWLHMSIpWl7Ng+YvQeSC0OybqNJJiyAuAJdWu5S/D3aeFT/GF8JOtT/uLRUZSlVlmpCSy27c0J7PNm2DG7+AfJ8LSuTD7xqgTSUoRnmaSVPM+WRDecr12STg+5Bzo9ZtoM0lKGZYZSTWmHqXUmzMaXhgHQRns0Tw8pdTh5KijSUohlhlJNebyjP9S7/lHw0HnH8Apf4KcxtGGkpRyvGZGUo2ZWNqP8iYd4Yd3ww/utMhIqhGWGUnVZ8078PT1EAQArCOPZT+aCZ0GRJtLUkqzzEjafeVl8Pxf4G9Hw9w/8fBjj1Z8qc/Nz/LgvOURhpOU6iwzknbP2vfhHyfBzKuhrIRVbQdw9fysii+XB3Dlw4tZVbAxwpCSUpllRlLVlJfDSxNgQk/4+BXIzIX+t7K0502UB1uvWhYELFtbFE1OSSnPu5kkVc2UwfDGA+H79sdC/1uhYRvaFWwkLcFWhSY9kaBtEz9vSVLNcGZGUtUccjZk7gH9boJzp0DDNgC0zK/PqP6dKlZLS8DoMzvTMr9+REElpTrLjKRdU/AxvDtjy7j9sTBsEXT/Wfj5St8wsGvrivezLu/FWd33qqWQkuoiy4yUpFYXbIo6QigIYOG/4fYe8NAFsG7Zlq/twnNjWuRn11w2ScIyIyWVSQs+rnjfZ+yc6G9pXr8a7j8bHrkYiguhacfwNmxJSiKWGSlJrCrYyDVT36wYR3pLcxDAoofgtsPh3ScgPRN6XwMXToc996n9PJK0A97NJCWJpWs3bPeW5lq9eDYI4KEL4c2Hw3GLLnDG36D5gbWXQZIqwTIjJYl2TRokxy3NiQQ02hvSMuCYX8PRv4T0erWbQZIqwdNMUpKI9Jbmoi+2vrD32JHwi9lw7AiLjKSkZ5mRkkgktzS/Oz28U+m/F0BZabgsIwtaHFTzP1uSqoGnmaQkVeO3NG8qgOlXwsJ7w3FWLqxfVfHwO0mKC8uMVBd98DQ8cgkUfgwkoMcQOP63UM+n9EqKH8uMVJds3gjTr4L5d4bjRu1gwO2w95HR5pKk3RCra2bGjBlDIpFg2LBhUUeR4imtHqx6LXzf/ecw+HmLjKTYi83MzLx585g4cSJdunSJOooULyVFkJYeXtSbngEDJsD6leFnK0lSCojFzMxXX33FOeecwx133EGjRo2ijiPFx4p58Lej4ZnRW5Y13b/Gi0xOZgbLbujHshv6kZMZm/9mkhRTsSgzQ4YMoV+/fvTp02en6xYXF1NYWLjVS6pzSoth5jXwj77w+fvhRxOUbIg6lSTViKT/T6YHHniAV199lXnz5u3S+mPGjGHUqFE1nEpKYisXwuTB8Nnb4bjLWXDyjZDZINpcklRDknpmZsWKFVx22WXce++9ZGfv2jM3Ro4cSUFBQcVrxYoVNZxSShKlJeHppDt6h0WmQVM4699w5kSo7+lZSakrqWdmFixYwJo1a+jatWvFsrKyMubOncutt95KcXEx6enpW22TlZVFVlZWbUeVovfVp/DibRCUwYEDoN9YaLBn1KkkqcYldZnp3bs3ixYt2mrZBRdcQMeOHbniiiu+U2SkOicIwg+GhPDJvf3GhncsdR4YbS5JqkVJXWZyc3Pp3LnzVssaNGjAnnvu+Z3lUp3z2bvwyMVw3FWwz3HhsoPPijaTJEUgqcuMVNd8fUvzDpWXwUvj4enfQ+kmmPFbGPTclhkaSapjYldmZs+eHXUEKTqffwCPDIHlL4bjfXpD/1ssMpLqtNiVGalOKi8PP09p5tWwuQgy94C+10PX8y0ykuo8y4wUB0tnw+O/Ct+3PRpOvw0a7R1pJElKFpYZKQ7aHwcH/xhaHQrdfwZpSf2IKEmqVf5GlJJR4UqY9HPYsDYcJxJwxng4/BcWGUn6FmdmpGQSBPDGg/DEb2BTQbhs4B3RZpKkJGeZkZLFV2vg0WGwZFo4bvV9OPqXkUaSpDiwzEjJYPHDMO2XsPELSKsHx14BRw0Pn+YrSdohf1NKUZv/D3hsePi++UHhtTEtDoo2kyTFiFcSSlHrPBAatYVjfgM/f9oiI0mV5MyMVNs2roOF90KPoeFdStn5cPFLUK9+1MkkKZYsM1Jtem8WTB0K61dBVm74BF+wyEjSbrDMSLWheD1MvwpevTscN94Hmh0YbSZJShGWGammfTgHHhkKBcvD8eGDoffVkJkTbS5JShGWGakmPXczzLo2fN9wbxhwO7TtGWkkSUo1lhmlhKKSUg68ejoAb113IjmZSfJPe++jIJEOXc+DE34PWXtEnUiSUk6S/MaXUsTmjfDJgi2zL20Og0vmQ+P20eaSpBTmc2ak6vLxAvjbMfCvM+Gzd7cst8hIUo2yzEi7q7QYnvo93HkCrH0X6jeEDWuiTiVJdYanmaTdseoNmDIYPl0cjjv/AE75E+Q0jjaXJNUhlhmpqp4dC8/8AcpLIWdPOPVmOPD0qFNJUp1jmZGqqrwsLDIHnAb9boY9mkadSJLqJMuMtKvKy+CrNZDXMhz3HA4tOsP+J4WfsSRJioQXAEu7Yu378I8T4d4zwwt+AdIzoMPJFhlJiphlRtqR8nJ4aTxM6Akfz4OCj+HTN6NOJUn6Bk8zSduzbhlMGQIfPReO2x8H/W+Bhm0ijSVJ2pplRvq2IIAFd8H038LmDVCvAfT9PXS70FNKkpSELDPStwUBLJoUFpm9e8Lpt0LjdlGnkiRth2VGKWd1wSbaN63kBzoGQXibdXo9SEuDAbfBkifhsF+EY0lS0vK3tFLCpAUfV7zvM3YOD85bvusbr18N9/8/mH7llmWN2sIRgywykhQD/qZW7K0q2Mg1U7fcYVQewJUPL2ZVwcYdbxgE8MZ/4bbD4d0nYcHdULiyhtNKkqqbp5kUe0vXbqA82HpZWRCwbG0RLfPrb3ujrz6DacPh7UfDccuDYcAEyGtVs2ElSdXOMqPYa9ekAWkJtio06YkEbZvkbHuDt6bCY8OhaC2kZcAxv4GjLw+vl5EkxY6nmRR7LfPrM6p/p4pxWgJGn9l527MyG7+EqZeERaZZJ/j503DsFRYZSYoxZ2aUEgZ2bc3vHgmvm5l1ea/t381UvyH0uwnWvAW9roCMrNoLKUmqEZYZpZwW+dlbBpsKwruUOvSDjqeEyw76QTTBJEk1wjKj1PXB0/DIUCj8BN6bBfscD/Wyd76dJClWLDNKPSVfwYxRMP8f4bhROxgw3iIjSSnKMqOUcnjibbL/PhK+/Chc0P3ncMIoyGwQbTBJUo2xzChl7Jv4mPszryftywDy24SfqdT+2KhjSZJqmGVGKeP9oDWTyo7mjG57k3HSaMjOizqSJKkW+JwZxdfmTfD0H6BwVcWiK0p/Qckp4ywyklSHODOjeFq5ECYPgs/egVWvwQ/uA6Dcfi5JdY5lRvFSWgJz/wTP3gRBGTRoCt8/DxKJqJNJkiJimVF8rF4MUwbB6kXhuNMZcMpN0GBPKCmNNpskKTKWGcXD+0/BfWdB+Wao3zj8SILOZ0adSpKUBCwzioe9joCGbaBpRzh1HOQ23+rLOZkZLLuhXzTZJEmRsswoOZWXweKHofNASEsLH3p34Qxo0MTrYyRJW7HMKPl8/gFMuRhWvARFa+GIweHyPZpGm0uSlJQsM0oe5eUw/06YeTVsLoLMPSDL58VIknbMMqPksO4jmDoUls4Nx22PhtNvg0Z7R5tLkpT0LDOK3ltTYcrg8NOu6+VAn1HQ/WfhtTKSJO2EZUbRa9QWSjdBmyNgwO2w5z5RJ5IkxYhlRrUvCGDN29D8wHDcsgtc8AR8ryukpUebTZIUO87jq3at/xQeOAf+djSsemPL8jaHWWQkSVVimVHtWfww3H4ELJkGJGD1GzvdRJKknfE0k2rehs/h8V/Cm5PDcYuDYMAEaNE52lySpJRgmVHNWvIETL0ENnwGiXQ45ldw9K8gIzPqZJKkFGGZUc36YmlYZJoeAGeMh1aHRp1IkpRiLDOqfsXrISs3fH/4IEivB4eeC/Wyo80lSUpJXgCs6rOpMDyl9LdeUFIULktLg8N+bpGRJNUYy0wdV1RSStsR02g7YhpFJaVV/0YfzoHxR8Gr98AXH8IHT1VfSEmSdsDTTNo9JRtg5jUw745w3HBvGDAe2h4VbS5JUp1hmVHVffRi+JlK65aG427/BydcB1l7RJtLklSnWGZUdc+PC4tM3vfg9Fthn+OjTiRJqoMsM6qcIIBEInx/6jh49s/Q+2rIzo80liSp7krqC4DHjBlD9+7dyc3NpVmzZgwYMIAlS5ZEHatuKi2Gp66DRy/dsiyvJfS7ySIjSYpUUpeZOXPmMGTIEF566SVmzpxJaWkpffv2ZcOGDVFHq1tWvQ4Tj4NnbwrvVlq5MOpEkiRVSOrTTE8++eRW47vuuotmzZqxYMECjjnmmIhS1SFlm8MCM/dPUF4KOU3g1LE+xVeSlFSSusx8W0FBAQCNGzfe7jrFxcUUFxdXjAsLC2s8V0pa8zZMviiclQE4oD/0Gwt7NI02lyRJ35LUp5m+KQgCLr/8cnr27Ennztv/tOUxY8aQn59f8WrTpk0tpkwRZZvh3z8Mi0x2Qxh4J/zoHouMJCkpxabMDB06lDfeeIP7779/h+uNHDmSgoKCiteKFStqKWEKSa8HJ98I+58EQ16Gg36w5Q4mSZKSTCxOM11yySVMnTqVuXPn0rp16x2um5WVRVZWVi0lSw0JyrkgfTrpb22EQ34QLuzYDzqcYomRJCW9pC4zQRBwySWXMHnyZGbPnk27du2ijpRyEuuWcVu9v9CQDXz++LO02u9YaNDkf1+0yEiSkl9Sl5khQ4Zw33338cgjj5Cbm8vq1asByM/Pp379+hGni7kggPl38shjM7hy82WUk0ba5oAxbxZx1mFRh5MkadclgiAIog6xPYntzAzcddddnH/++bv0PQoLC8nPz6egoIC8vLxqTBdjX66AqUNZ9cEbHFX8V8q/celUeiLBcyOOo2W+ZVGSFJ3K/P1O6pmZJO5Z8bVhLYw/CooLWJo4ZKsiA1AWBCxbW2SZkSTFRmzuZlI1adAEuvwIWh9Gu/MmkPatya/0RIK2TXKiySZJUhVYZlJdEMAb/4V1y7Ys63s9XPgkLdsdwKj+nSoWpyVg9JmdnZWRJMWKZSaVffUZPPgTePhnMGUIlJeHy+tlQ1o6AAO7brnVfdblvTir+15RJJUkqcqS+poZ7Ya3HoHHhkPR55CWAe17QVDOjvpri/zs2ssnSVI1scykmqIv4PFfw+KHwnHzzjBgPLTsEm0uSZJqiGUmlax5G+45Hb76FBLp0HM49LoCMjKjTiZJUo2xzKSSxu0hZ8/wwyHPGA/f6xp1IkmSapxlJu4+ehFad4f0DMjIgrMfgD2ahxf5SpJUB3g3UxUVlZTSdsQ02o6YRlFJae0HKF4Pjw6Du06C52/esrzR3hYZSVKd4sxMHC19Fh65GL5cHo43fhlpHEmSomSZiZOSInhqFLw8IRzn7wWn3xredi1JUh1lmYmLla/BQxfCFx+E467nh0/yzcqNMpUkSZGzzMRFvRwo+BhyW0H/W2C/PlEnkiQpKVhmktn61ZDbInzfdH/4f/dB625Qv2GksSRJSibezZSMSkvg6eth3EGw4pUty/frY5GRJOlbnJlJNqsXweTB8OmicPzOY9DmsBr7cTmZGSy7oV+NfX9JkmqaZSZZlJWGz4uZfSOUb4b6jaHfTdD5zKiTSZKU1CwzyeCzJTB5EKx8NRx3PBVOvRn2aBZtLkmSYsAykwyWvxgWmex8OPlP0OVHkEhEnUqSpFiwzESlrDT8PCWA758X3rn0/Z9CXqtoc0mSFDPezVTbysvh5b/B+B6wqSBclkjAsSMsMpIkVYFlphqsLti0ayuu+wju6Q9P/AbWvgsL7q7ZYJIk1QGWmSqatODjivd9xs7hwXnLt79yEMD8u2D8kbDs2fBpvqf8GXoMrYWkkiSlNq+ZqYJVBRu5ZuqbFePyAK58eDHH7N+Ulvn1t1654BOYegl88FQ43qsHnH4b7LlPLSaWJCl1WWaqYOnaDZQHWy8rCwKWrS36bpl55g9hkUnPgt5XwxGDIS299sJKkpTiLDNV0K5JA9ISbFVo0hMJ2jbJ+e7KJ/weij4P/7fp/rUXUpKkOsJrZqqgZX59RvXvVDFOS8DoMzuHszKLJ8Gjl4XXyQA02BN+/KBFRpKkGmKZqaKBXVtXvJ91eS/OOrAB/Oc8eOhCWPBPePfJ6MJJklSHeJqpGnxv9VPw5C9hw2eQlgFH/wr27RN1LEmS6gTLzG7I4yuurXcPWZOeCxc0PQDOmACtDok0lyRJdYllpqqCgH9m/pHvp71PkEgjcdRlcOxIyMiKOpkkSXWK18xUVSLBTaU/5P3yVhT/9HHoc61FRpKkCFhmdsPz5QdxYsmNlH+ve9RRJEmqsywzu6kMH4AnSVKULDOSJCnWLDOSJCnWLDOSJCnWLDOSJCnWLDOSJCnWfGheFeVkZrDshn5Rx5Akqc5zZkaSJMWaZUaSJMWaZUaSJMWaZUaSJMWaZUaSJMWaZUaSJMWaZUaSJMWaZUaSJMWaZUaSJMWaZUaSJMWaZUaSJMWaZUaSJMWaZUaSJMWaZUaSJMWaZUaSJMVaRtQBaloQBAAUFhZGnESSJO2qr/9uf/13fEdSvsysX78egDZt2kScRJIkVdb69evJz8/f4TqJYFcqT4yVl5ezcuVKcnNzSSQS1fq9CwsLadOmDStWrCAvL69av3cycP/iL9X30f2Lv1TfR/ev6oIgYP369bRq1Yq0tB1fFZPyMzNpaWm0bt26Rn9GXl5eSv4j/Zr7F3+pvo/uX/yl+j66f1WzsxmZr3kBsCRJijXLjCRJijXLzG7IysrimmuuISsrK+ooNcL9i79U30f3L/5SfR/dv9qR8hcAS5Kk1ObMjCRJijXLjCRJijXLjCRJijXLjCRJijXLzA7cfvvttGvXjuzsbLp27cqzzz67w/XnzJlD165dyc7Opn379kyYMKGWklZdZfZx9uzZJBKJ77zeeeedWky86+bOnctpp51Gq1atSCQSTJkyZafbxOkYVnb/4nb8xowZQ/fu3cnNzaVZs2YMGDCAJUuW7HS7uBzDquxf3I7h+PHj6dKlS8UD1Xr06METTzyxw23icvyg8vsXt+P3bWPGjCGRSDBs2LAdrhfFMbTMbMeDDz7IsGHDuOqqq1i4cCFHH300J598MsuXL9/m+kuXLuWUU07h6KOPZuHChVx55ZVceumlTJo0qZaT77rK7uPXlixZwqpVqype++23Xy0lrpwNGzZw8MEHc+utt+7S+nE7hpXdv6/F5fjNmTOHIUOG8NJLLzFz5kxKS0vp27cvGzZs2O42cTqGVdm/r8XlGLZu3ZobbriB+fPnM3/+fI4//nhOP/103nzzzW2uH6fjB5Xfv6/F5fh907x585g4cSJdunTZ4XqRHcNA23TYYYcFgwYN2mpZx44dgxEjRmxz/d/85jdBx44dt1p20UUXBUcccUSNZdxdld3HZ555JgCCdevW1UK66gUEkydP3uE6cTyGX9uV/Yvz8QuCIFizZk0ABHPmzNnuOnE+hruyf3E/hkEQBI0aNQr+/ve/b/NrcT5+X9vR/sX1+K1fvz7Yb7/9gpkzZwa9evUKLrvssu2uG9UxdGZmG0pKSliwYAF9+/bdannfvn154YUXtrnNiy+++J31TzzxRObPn8/mzZtrLGtVVWUfv3booYfSsmVLevfuzTPPPFOTMWtV3I5hVcX1+BUUFADQuHHj7a4T52O4K/v3tTgew7KyMh544AE2bNhAjx49trlOnI/fruzf1+J2/IYMGUK/fv3o06fPTteN6hhaZrZh7dq1lJWV0bx5862WN2/enNWrV29zm9WrV29z/dLSUtauXVtjWauqKvvYsmVLJk6cyKRJk3j44Yfp0KEDvXv3Zu7cubURucbF7RhWVpyPXxAEXH755fTs2ZPOnTtvd724HsNd3b84HsNFixaxxx57kJWVxaBBg5g8eTIHHnjgNteN4/GrzP7F8fg98MADvPrqq4wZM2aX1o/qGKb8p2bvjkQisdU4CILvLNvZ+ttankwqs48dOnSgQ4cOFeMePXqwYsUK/vznP3PMMcfUaM7aEsdjuKvifPyGDh3KG2+8wXPPPbfTdeN4DHd1/+J4DDt06MBrr73Gl19+yaRJkzjvvPOYM2fOdv/gx+34VWb/4nb8VqxYwWWXXcaMGTPIzs7e5e2iOIbOzGxDkyZNSE9P/84MxZo1a77TOL/WokWLba6fkZHBnnvuWWNZq6oq+7gtRxxxBO+99151x4tE3I5hdYjD8bvkkkuYOnUqzzzzDK1bt97hunE8hpXZv21J9mOYmZnJvvvuS7du3RgzZgwHH3wwf/nLX7a5bhyPX2X2b1uS+fgtWLCANWvW0LVrVzIyMsjIyGDOnDn89a9/JSMjg7Kysu9sE9UxtMxsQ2ZmJl27dmXmzJlbLZ85cyZHHnnkNrfp0aPHd9afMWMG3bp1o169ejWWtaqqso/bsnDhQlq2bFnd8SIRt2NYHZL5+AVBwNChQ3n44Yd5+umnadeu3U63idMxrMr+bUsyH8NtCYKA4uLibX4tTsdve3a0f9uSzMevd+/eLFq0iNdee63i1a1bN8455xxee+010tPTv7NNZMewRi8vjrEHHnggqFevXnDnnXcGb731VjBs2LCgQYMGwbJly4IgCIIRI0YE5557bsX6H374YZCTkxMMHz48eOutt4I777wzqFevXvDQQw9FtQs7Vdl9vPnmm4PJkycH7777brB48eJgxIgRARBMmjQpql3YofXr1wcLFy4MFi5cGADB2LFjg4ULFwYfffRREATxP4aV3b+4Hb/BgwcH+fn5wezZs4NVq1ZVvIqKiirWifMxrMr+xe0Yjhw5Mpg7d26wdOnS4I033giuvPLKIC0tLZgxY0YQBPE+fkFQ+f2L2/Hblm/fzZQsx9AyswO33XZbsPfeeweZmZnB97///a1umTzvvPOCXr16bbX+7Nmzg0MPPTTIzMwM2rZtG4wfP76WE1deZfbxxhtvDPbZZ58gOzs7aNSoUdCzZ89g2rRpEaTeNV/fBvnt13nnnRcEQfyPYWX3L27Hb1v7BgR33XVXxTpxPoZV2b+4HcMLL7yw4vdL06ZNg969e1f8oQ+CeB+/IKj8/sXt+G3Lt8tMshzDRBD878ocSZKkGPKaGUmSFGuWGUmSFGuWGUmSFGuWGUmSFGuWGUmSFGuWGUmSFGuWGUmSFGuWGUmSFGuWGUmSFGuWGUmSFGuWGUmx8tlnn9GiRQtGjx5dsezll18mMzOTGTNmRJhMUlT8bCZJsfP4448zYMAAXnjhBTp27Mihhx5Kv379GDduXNTRJEXAMiMploYMGcKsWbPo3r07r7/+OvPmzSM7OzvqWJIiYJmRFEsbN26kc+fOrFixgvnz59OlS5eoI0mKiNfMSIqlDz/8kJUrV1JeXs5HH30UdRxJEXJmRlLslJSUcNhhh3HIIYfQsWNHxo4dy6JFi2jevHnU0SRFwDIjKXZ+/etf89BDD/H666+zxx57cNxxx5Gbm8tjjz0WdTRJEfA0k6RYmT17NuPGjeNf//oXeXl5pKWl8a9//YvnnnuO8ePHRx1PUgScmZEkSbHmzIwkSYo1y4wkSYo1y4wkSYo1y4wkSYo1y4wkSYo1y4wkSYo1y4wkSYo1y4wkSYo1y4wkSYo1y4wkSYo1y4wkSYo1y4wkSYq1/w+B+I5fUZ/qLQAAAABJRU5ErkJggg==",
+ "text/plain": [
+ "<Figure size 640x480 with 1 Axes>"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
"source": [
- "def intervall_minus(mu):\n",
- " return scipy.integrate.quad(lambda x: scipy.stats.poisson.pmf(18,x), mu, 18)[0]\n",
+ "#hideme\n",
+ "import numpy as np\n",
+ "import scipy.stats as stats\n",
+ "import matplotlib.pyplot as plt \n",
"\n",
- "def intervall_plus(mu):\n",
- " return scipy.integrate.quad(lambda x: scipy.stats.poisson.pmf(18,x), 18, mu)[0]\n",
+ "def f(x):\n",
+ " return 2*x + 1\n",
"\n",
- "print(\"naiv:\", 18-np.sqrt(18), 18+np.sqrt(18))\n",
- "mu_minus = scipy.optimize.brentq(lambda x: intervall_minus(x)-0.6827/2, 0,18)\n",
- "mu_plus = scipy.optimize.brentq(lambda x: intervall_plus(x)-0.6827/2, 18,40)\n",
- "print(\"symmetrisch in P:\", mu_minus, mu_plus) \n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "52f8f412-cb94-4c32-95b6-03f178cbcdb6",
- "metadata": {},
- "outputs": [],
- "source": [
- "scipy.integrate.quad(lambda x: scipy.stats.poisson.pmf(18,x), mu_minus, mu_plus)"
+ "n = 10\n",
+ "xs = np.linspace(0,4,n)\n",
+ "sigma_y=0.4\n",
+ "ys = stats.multivariate_normal.rvs(f(xs), np.eye(n)*sigma_y**2, 1, random_state=42)\n",
+ " \n",
+ "x_axis = np.linspace(0,4,100)\n",
+ "plt.errorbar(xs,ys,yerr=sigma_y,fmt=\".\")\n",
+ "plt.plot(x_axis, f(x_axis),'--')\n",
+ "plt.xlabel(\"x\")\n",
+ "plt.ylabel(\"y\")\n",
+ "plt.savefig(\"line.png\")\n",
+ "plt.show()"
]
},
{
"cell_type": "markdown",
- "id": "f7ef04e8-1dc9-4e8b-b597-39d277bde591",
- "metadata": {
- "slideshow": {
- "slide_type": "slide"
- },
- "tags": []
- },
- "source": [
- "## Konfidenz-Intervalle: Poisson"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "e0a6084d-dd5d-4ffc-8849-a61a8327fb4e",
+ "id": "64e67452-e6bd-442c-a174-e12bdb18dba0",
"metadata": {
"slideshow": {
"slide_type": ""
},
"tags": []
},
- "outputs": [],
- "source": [
- "def intervall_minus(k, mu):\n",
- " return scipy.integrate.quad(lambda x: scipy.stats.poisson.pmf(k,x), mu, k)[0]\n",
- "\n",
- "def intervall_plus(k, mu):\n",
- " return scipy.integrate.quad(lambda x: scipy.stats.poisson.pmf(k,x), k, mu)[0]\n",
- "\n",
- "def conv_limits(k, c):\n",
- " c_low = intervall_minus(k, 0)\n",
- " c_up = c / 2\n",
- " if c_low < c/2:\n",
- " mu_minus = 0\n",
- " else:\n",
- " c_low = c / 2\n",
- " mu_minus = scipy.optimize.brentq(lambda x: intervall_minus(k, x)-c_low, 0,k)\n",
- " c_up = c - c_low\n",
- " mu_plus = scipy.optimize.brentq(lambda x: intervall_plus(k, x)-c_up, k, 2*k + 20)\n",
- " return mu_minus, mu_plus\n",
- "\n",
- "conv_limits(3, 0.6827)\n",
- "\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "382e5a13-5f5f-44e5-90dd-7272fbdb171a",
- "metadata": {},
- "outputs": [],
"source": [
- "ks = np.arange(0, 20)\n",
- "\n",
- "limits = np.zeros((len(ks), 2))\n",
- "for i,k in enumerate(ks):\n",
- " limits[i] = conv_limits(k, 0.6827)\n",
- "#print(limits)\n",
- "plt.plot(ks, limits[:,0], 'b_')\n",
- "plt.plot(ks, limits[:,1], 'r_')\n",
- "plt.plot(ks, ks-np.sqrt(ks), 'b.')\n",
- "plt.plot(ks, ks+np.sqrt(ks), 'r.')\n",
- "plt.xlabel(\"$k$\")\n",
- "plt.ylabel(r\"$\\mu$\")\n",
- "plt.grid()\n",
- "\n",
- "print(limits[2])\n"
+ "### Methode der kleinsten Quadrate \n",
+ "\n",
+ "$$\\chi^2 = \\sum_i \\left(\\frac{y_i - \\hat y(x)}{\\sigma_i}\\right)^2$$\n",
+ "quantifiziert die Übereinstimmung von Modell zu Daten \n",
+ "$\\rightarrow$ $\\hat m$ und $\\hat a$ sollten $\\chi^2$ minimieren.\n",
+ "<img src=\"./figures/11/line.png\" style=\"width:80.0%\"\n",
+ "alt=\"image\" />"
]
},
{
"cell_type": "markdown",
- "id": "61d94d51-fc1c-4ddb-9cbc-c416a28158f8",
+ "id": "bdc0949d-791e-4e29-9adc-169578e62821",
"metadata": {
"slideshow": {
"slide_type": "slide"
@@ -438,82 +359,49 @@
"tags": []
},
"source": [
- "## Konfidenzintervall: Gauß"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "cb5826ae-29df-41c3-ad8f-dde61b5d90ba",
- "metadata": {},
- "outputs": [],
- "source": [
- "def intervall_minus(x, mu, sigma):\n",
- " return scipy.integrate.quad(lambda y: scipy.stats.norm.pdf(x,y, sigma), mu, x)[0]\n",
+ "### Methode der kleinsten Quadrate II \n",
"\n",
- "def intervall_plus(x, mu, sigma):\n",
- " return scipy.integrate.quad(lambda y: scipy.stats.norm.pdf(x, y, sigma), x, mu)[0]\n",
+ "Minimiere\n",
+ "$\\chi^2 = \\sum_i \\left(\\frac{y_i - \\hat y(x)}{\\sigma_i}\\right)^2 = \\sum_i \\frac{(y_i - m x_i - a)^2}{\\sigma_i^2}$:\n",
"\n",
- "def conv_limits(x, c, sigma):\n",
- " c_low = c / 2\n",
- " mu_minus = scipy.optimize.brentq(lambda mu: intervall_minus(x, mu, sigma)-c_low, x - 10*sigma,x)\n",
- " c_up = c - c_low\n",
- " mu_plus = scipy.optimize.brentq(lambda mu: intervall_plus(x, mu, sigma)-c_up, x, x + 10*sigma)\n",
- " return mu_minus, mu_plus\n",
+ "Erste Ableitung ist Null: \n",
"\n",
- "conv_limits(3, 0.6827, 1)\n"
+ "$$\\begin{aligned}\n",
+ " \\frac{d\\chi^2}{dm} &=& -2\\sum_i x_i\\frac {y_i -\\hat m x_i - \\hat a}{\\sigma_i^2} = 0\\\\\n",
+ " \\frac{d\\chi^2}{da} &=& -2\\sum_i \\frac{y_i - \\hat m x_i - \\hat a}{\\sigma_i^2} = 0 \\\\\n",
+ " \\sum_i\\frac{x_iy_i}{\\sigma_i^2} - \\hat m \\sum_i\\frac{x_i^2}{\\sigma_i^2}- \\hat a \\sum_i \\frac{x_i}{\\sigma_i^2} &=& 0 \\\\\n",
+ " \\sum_i\\frac{y_i}{\\sigma_i^2} - \\hat m \\sum_i\\frac{x_i}{\\sigma_i^2}- \\hat a \\sum_i \\frac{1}{\\sigma_i^2} &=& 0 \n",
+ "\\end{aligned}$$"
]
},
{
- "cell_type": "code",
- "execution_count": null,
- "id": "e87da118-0c11-4345-a656-096de6f980e9",
+ "cell_type": "markdown",
+ "id": "322c67fc-ddd1-4830-a5f5-6f118ef54c4c",
"metadata": {
"slideshow": {
- "slide_type": ""
+ "slide_type": "slide"
},
"tags": []
},
- "outputs": [],
- "source": [
- "xs = np.linspace(0,20,100)\n",
- "sigma = 1\n",
- "limits = np.zeros((len(xs), 2))\n",
- "for i,k in enumerate(xs):\n",
- " limits[i] = conv_limits(k, 0.6827, sigma)\n",
- "#print(limits)\n",
- "plt.plot(xs, limits[:,0], 'b')\n",
- "plt.plot(xs, limits[:,1], 'r')\n",
- "plt.plot(xs, xs-sigma, 'b.')\n",
- "plt.plot(xs, xs+sigma, 'r.')\n",
- "plt.xlabel(\"$x$\")\n",
- "plt.ylabel(r\"$\\mu$\")\n",
- "plt.grid()\n",
- "\n",
- "print(limits[2])"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "202144ae-1057-41fc-8757-7f6f40c75868",
- "metadata": {},
- "outputs": [],
"source": [
- "print(scipy.integrate.quad(lambda y: scipy.stats.norm.pdf(3, y, 1), 2, 3)[0], scipy.integrate.quad(lambda x: scipy.stats.norm.pdf(x, 3, 1), 2,3)[0])"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "24359dd4-7289-4235-872e-1138d5fe0b0b",
- "metadata": {},
- "source": [
- "Für Gauß $G(x, \\mu, \\sigma) = G(\\mu, x, \\sigma)$: $\\int_{\\mu_-}^{\\mu_+} G(x, \\mu, \\sigma) d\\mu = \\int_{\\mu_-}^{\\mu_+} G(\\mu, x, \\sigma) d\\mu$"
+ "### Methode der kleinsten Quadrate III \n",
+ "\n",
+ "Minimiere\n",
+ "$\\chi^2 = \\sum_i \\left(\\frac{y_i - \\hat y(x)}{\\sigma_i}\\right)^2 = \\sum_i \\frac{(y_i - m x_i - a)^2}{\\sigma_i^2}$: \n",
+ "$$\\begin{aligned}\n",
+ " \\sum_i\\frac{x_iy_i}{\\sigma_i^2} - \\hat m \\sum_i\\frac{x_i^2}{\\sigma_i^2}- \\hat a \\sum_i \\frac{x_i}{\\sigma_i^2} &=& 0 \\\\\n",
+ " \\sum_i\\frac{y_i}{\\sigma_i^2} - \\hat m \\sum_i\\frac{x_i}{\\sigma_i^2}- \\hat a \\sum_i \\frac{1}{\\sigma_i^2} &=& 0 \n",
+ "\\end{aligned}$$ mit\n",
+ "$\\frac{1}{\\sum_i 1/\\sigma_i^2} \\sum_i \\frac{f}{\\sigma_i^2} = \\langle f \\rangle$: \n",
+ "$$\\begin{aligned}\n",
+ " \\langle xy \\rangle -\\langle x^2 \\rangle \\hat m& - \\langle x \\rangle \\hat a&= 0\\\\\n",
+ " \\langle y \\rangle - \\langle x \\rangle \\hat m& - \\hat a& = 0 \n",
+ "\\end{aligned}$$"
]
},
{
"cell_type": "markdown",
- "id": "10f482d3-6667-4250-aefa-5df896e46936",
+ "id": "fb25687c-4410-4281-b540-39369732fb26",
"metadata": {
"slideshow": {
"slide_type": "slide"
@@ -521,315 +409,153 @@
"tags": []
},
"source": [
- "# Konfidenzregionen\n",
+ "### Methode der kleinsten Quadrate IV \n",
"\n",
- "Erst einmal Konfidenzen für Gauß-Verteilung zum Intervall $[\\mu - z\\sigma, \\mu - z\\sigma]$:"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "580b83dc-3cb0-40df-bc77-8ad48d90fce5",
- "metadata": {},
- "outputs": [],
- "source": [
- "def conv_gaus(z):\n",
- " return scipy.stats.norm.cdf(z) - scipy.stats.norm.cdf(-z)\n",
- "\n",
- "for z in [1,2,3,4,5]:\n",
- " print(z, conv_gaus(z))"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "3904ad68-2b4d-4207-a7e8-35def0372238",
- "metadata": {},
- "source": [
- "$z$ für 68\\%, 90\\%, 95\\% und 99\\%: "
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "5a71ee06-db32-4349-a3da-ae0e394fbd74",
- "metadata": {},
- "outputs": [],
- "source": [
- "for c in [0.68, 0.90, 0.95, 0.99]:\n",
- " print(c, scipy.optimize.brentq(lambda z: conv_gaus(z)-c,0, 10))"
+ "$$\\begin{aligned}\n",
+ " \\hat m&=&\\frac{\\langle xy \\rangle - \\langle y \\rangle\\langle x \\rangle}{\\langle x^2 \\rangle - \\langle x \\rangle^2} = \\frac{1}{\\sum_i 1/\\sigma_i^2} \\sum_i \\frac{x_i - \\langle x \\rangle}{\\sigma_i^2(\\langle x^2 \\rangle - \\langle x \\rangle^2)}y_i\\\\\n",
+ " \\hat a &=& \\frac{ \\langle y \\rangle \\langle x^2 \\rangle- \\langle y \\rangle \\langle x \\rangle^2- \\langle x \\rangle \\langle xy \\rangle+ \\langle y \\rangle \\langle x \\rangle^2}{ \\langle x^2 \\rangle- \\langle x \\rangle^2}\\\\\n",
+ " &=& \\frac{ \\langle y \\rangle \\langle x^2 \\rangle - \\langle x \\rangle \\langle xy \\rangle}{ \\langle x^2 \\rangle - \\langle x \\rangle^2} = \\frac{1}{\\sum_i 1/\\sigma_i^2} \\sum_i \\frac{\\langle x^2 \\rangle - \\langle x \\rangle x_i}{\\sigma_i^2(\\langle x^2 \\rangle - \\langle x \\rangle^2)}y_i\n",
+ "\\end{aligned}$$"
]
},
{
"cell_type": "markdown",
- "id": "2d5f1dab-cebd-45fe-934e-36cd06ba2782",
+ "id": "b08ade05-e5eb-4a7e-9315-31a45c51aadb",
"metadata": {
"slideshow": {
- "slide_type": "slide"
+ "slide_type": ""
},
"tags": []
},
"source": [
- "## Regionen in zwei oder drei Dimensionen:\n",
+ "## Fehler\n",
"\n",
- "$\\vec x^T = (x_1, x_2, \\dots, x_n)$ \n",
- "\n",
- "$x_i$ seien unabhängige Zufallsvariablen und Gauß verteilt:"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "73c11f8d-c8e4-4dc7-ad37-4ca4a33c55ea",
- "metadata": {},
- "outputs": [],
- "source": [
- "def conv_gaus_nd(z, n):\n",
- " return np.power(scipy.stats.norm.cdf(z) - scipy.stats.norm.cdf(-z),n)\n",
+ "### Fehler \n",
"\n",
- "for z in [1,2,3,4,5]:\n",
- " print(z, conv_gaus_nd(z, 1), conv_gaus_nd(z, 2), conv_gaus_nd(z, 3))"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "8dd1036c-edf9-40c4-be16-56025f49134b",
- "metadata": {},
- "source": [
- "$z$ für 68\\%, 90\\%, 95\\% und 99\\%: "
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "7f14d38e-1e41-460e-9876-2afb62bed08b",
- "metadata": {},
- "outputs": [],
- "source": [
- "for c in [0.68, 0.90, 0.95, 0.99]:\n",
- " print(c, scipy.optimize.brentq(lambda z: conv_gaus_nd(z, 1)-c,0, 10), scipy.optimize.brentq(lambda z: conv_gaus_nd(z, 2)-c,0, 10), scipy.optimize.brentq(lambda z: conv_gaus_nd(z,3)-c,0, 10))"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "9aca662e-bdd6-408e-8f28-852614e51577",
- "metadata": {},
- "source": [
- "1,2,3-$\\sigma$-Äquivalente in zwei und drei Dimensionen:"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "702daac8-81e8-406c-b70f-c3146c58a493",
- "metadata": {},
- "outputs": [],
- "source": [
- "for k in [1, 2, 3]:\n",
- " c = conv_gaus_nd(k, 1)\n",
- " print(c, scipy.optimize.brentq(lambda z: conv_gaus_nd(z, 2)-c,0, 10), scipy.optimize.brentq(lambda z: conv_gaus_nd(z,3)-c,0, 10))"
+ "$$\\begin{aligned}\n",
+ "V(\\hat m) = \\sum_i \\left(\\frac{d\\hat m}{y_i}\\sigma_i\\right)^2\\text{; }\\frac{d\\hat m}{y_i} & = & \\frac{1}{\\sum_i 1/\\sigma_i^2} \\frac{x_i - \\langle x \\rangle}{\\sigma_i^2(\\langle x^2 \\rangle - \\langle x \\rangle^2)} \\\\\n",
+ "V(\\hat a) = \\sum_i \\left(\\frac{d\\hat a}{y_i}\\sigma_i\\right)^2\\text{; }\\frac{d\\hat a}{y_i} & = & \\frac{1}{\\sum_i 1/\\sigma_i^2} \\frac{\\langle x^2 \\rangle - \\langle x \\rangle x_i}{\\sigma_i^2(\\langle x^2 \\rangle - \\langle x \\rangle^2)}\n",
+ "\\end{aligned}$$ $$\\begin{aligned}\n",
+ "V(\\hat m) &=& \\left(\\frac{1}{\\sum_i 1/\\sigma_i^2}\\right)^2 \\sum_i \\left(\\frac{x_i - \\langle x \\rangle}{\\sigma_i^2(\\langle x^2 \\rangle - \\langle x \\rangle^2)}\\right)^2 \\sigma_i^2 \\\\\n",
+ "&=& \\frac{1}{\\sum_i 1/\\sigma_i^2} \\frac{\\langle x^2 \\rangle - 2\\langle x \\rangle \\langle x \\rangle + \\langle x \\rangle^2}{(\\langle x^2 \\rangle - \\langle x \\rangle^2)^2} \n",
+ "= \\frac{1}{\\sum_i 1/\\sigma_i^2} \\frac{1}{\\langle x^2 \\rangle - \\langle x \\rangle^2} \\\\\n",
+ "V(\\hat a) &=& \\frac{1}{\\sum_i 1/\\sigma_i^2} \\frac{\\langle x^2 \\rangle^2 - 2\\langle x^2 \\rangle\\langle x \\rangle^2 + \\langle x^2 \\rangle\\langle x \\rangle^2}{(\\langle x^2 \\rangle - \\langle x \\rangle^2)^2}\n",
+ "= \\frac{1}{\\sum_i 1/\\sigma_i^2} \\frac{\\langle x^2 \\rangle}{\\langle x^2 \\rangle - \\langle x \\rangle^2}\n",
+ "\\end{aligned}$$"
]
},
{
"cell_type": "markdown",
- "id": "bba9c7fc-9e11-4441-90e0-f1d186a1726e",
+ "id": "d8828d04-0af8-4dc3-a846-24acbc9a0f8e",
"metadata": {
"slideshow": {
- "slide_type": "slide"
+ "slide_type": ""
},
"tags": []
},
"source": [
- "# Hypothesentest\n",
- "\n",
- "\n",
- "Hypothese: \"Die $k_i$ Tore pro Bundesligaspiel $i$ sind Poisson-verteilt mit einem gemeinsamen $\\mu = <k>$.\"\n",
- "\n",
- "Benötigt für den Test eine alternative Hypothese: \"Die Tore pro Bundesligaspiel $k_i$ sind Poisson verteilt mit $\\mu_i = ki$.\"\n",
- "\n",
- "Fehler 1. und 2. Art\n",
- "* Fehler 1. Art: Die Hypothese stimmt, wird aber verworfen.\n",
- " \n",
- " Irrtumswahrscheinlichkeit: $\\alpha$ (Signifikanzniveau, significance)\n",
- " \n",
- " Spezifität: $1-\\alpha$ (Effizienz) \n",
- "\n",
+ "### Korrelation \n",
"\n",
- "* Fehler 2. Art: Die Hypothese wird angenommen, ist aber falsch (falsch positiv).\n",
+ "$$\\begin{aligned}\n",
+ "V(\\hat m) &=& \\frac{1}{\\sum_i 1/\\sigma_i^2} \\frac{1}{\\langle x^2 \\rangle - \\langle x \\rangle^2} \\\\\n",
+ "V(\\hat a) &=& \\frac{1}{\\sum_i 1/\\sigma_i^2} \\frac{\\langle x^2 \\rangle}{\\langle x^2 \\rangle - \\langle x \\rangle^2}\\\\\n",
+ "\\text{cov}(\\hat m, \\hat a) &=&= \\frac{1}{\\sum_i 1/\\sigma_i^2} \\frac{\\langle (x-\\langle x \\rangle)(\\langle x^2 \\rangle - \\langle x \\rangle x)\\rangle}{(\\langle x^2 \\rangle - \\langle x \\rangle^2)^2}\\\\\n",
+ "&=& \\frac{1}{\\sum_i 1/\\sigma_i^2} \\frac{\\langle x^2 \\rangle \\langle x \\rangle - \\langle x \\rangle \\langle x^2 \\rangle - \\langle x \\rangle \\langle x^2 \\rangle + \\langle x \\rangle^2\\langle x \\rangle}{(\\langle x^2 \\rangle - \\langle x \\rangle^2)^2}\\\\\n",
+ "&=& - \\frac{1}{\\sum_i 1/\\sigma_i^2} \\frac{\\langle x \\rangle}{\\langle x^2 \\rangle - \\langle x \\rangle^2}\n",
+ "\\end{aligned}$$\n",
"\n",
- " Wahrscheinlichkeit für Fehler: $\\beta$\n",
- " \n",
- " Trennschärfe, Sensitivität: $1-\\beta$ (power)"
+ "### Beispiel in Jupyter "
]
},
{
"cell_type": "markdown",
- "id": "4e67fff2-211b-4eef-95e2-44e3b363098c",
+ "id": "424bdd1f-53bf-422b-bc4a-0702231b976d",
"metadata": {
"slideshow": {
- "slide_type": "slide"
+ "slide_type": ""
},
"tags": []
},
"source": [
- "## Beispiel:\n",
- "Gauß verteilte Zufallsvariable $x$ ($\\sigma = 1$)\n",
+ "### Minimales $\\chi^2$ \n",
"\n",
- "Hypothese: $\\mu = 0$\n",
- "\n",
- "Für welchen Bereich $x_- < x < x_+$ sollte man die Hypothese annehmen für Fehler 1. Art $\\alpha = 5\\%$?"
+ "$$\\begin{aligned}\n",
+ " \\chi^2 &=& \\sum_i \\frac{(y_i - \\hat m x_i - \\hat a)^2}{\\sigma_i^2} = \\sum_i \\frac{\\left[y_i - \\frac{\\langle xy \\rangle - \\langle y \\rangle\\langle x \\rangle}{\\langle x^2 \\rangle - \\langle x \\rangle^2} x_i - \\frac{ \\langle y \\rangle \\langle x^2 \\rangle - \\langle x \\rangle \\langle xy \\rangle}{ \\langle x^2 \\rangle - \\langle x \\rangle^2} \\right]^2}{\\sigma_i^2}\\\\\n",
+ " & = & \\sum_i \\frac{\\left[(\\langle x^2 \\rangle - \\langle x \\rangle^2)y_i - (\\langle xy \\rangle - \\langle y \\rangle\\langle x \\rangle)x_i - \\langle y \\rangle \\langle x^2 \\rangle + \\langle x \\rangle \\langle xy \\rangle\\right]^2}{\\sigma_i^2 ( \\langle x^2 \\rangle - \\langle x \\rangle^2)^2} \\\\\n",
+ " &=& \\dots\\\\\n",
+ "& =& (\\sum_i \\frac{1}{\\sigma_i^2}) V(y) ( 1- \\rho^2_{xy})\n",
+ "\\end{aligned}$$"
]
},
{
- "cell_type": "code",
- "execution_count": null,
- "id": "fb7d561b-51c8-4857-b0ea-81dec8facbda",
- "metadata": {},
- "outputs": [],
- "source": []
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "85e25147-c902-4434-b3dc-908486408c5b",
+ "cell_type": "markdown",
+ "id": "e89f415e-8836-4b97-893c-a7335c3a21e5",
"metadata": {
"slideshow": {
"slide_type": ""
},
"tags": []
},
- "outputs": [],
"source": [
- "x = np.linspace(-5,5,100)\n",
- "xout1 = np.linspace(-5, -1.9599639845401602)\n",
- "xout2 = np.linspace(1.9599639845401602, 5)\n",
- "plt.plot(x, scipy.stats.norm.pdf(x))\n",
- "plt.fill_between(xout1, scipy.stats.norm.pdf(xout1), np.zeros_like(xout1),color=\"c\")\n",
- "plt.fill_between(xout2, scipy.stats.norm.pdf(xout2), np.zeros_like(xout2),color=\"c\")\n",
+ "### Beispiel in Jupyter \n",
"\n",
- "plt.xlabel(\"$x$\")\n",
- "\n",
- "plt.grid()"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "6cbb0d56-981d-4dfa-9a7a-0d6d4a51a357",
- "metadata": {},
- "source": [
- "Fehler 2. Art:?"
+ "## In Python"
]
},
{
"cell_type": "markdown",
- "id": "cb79afd6-ba74-448b-9096-05cedeb842a6",
- "metadata": {},
- "source": [
- "Beispiel: wahres $\\mu = 3$:"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "603a31e1-1425-4eb6-9b12-c44ba31d06a6",
- "metadata": {},
- "outputs": [],
- "source": [
- "x = np.linspace(-5,8,200)\n",
- "xin = np.linspace(-1.9599639845401602, 1.9599639845401602)\n",
- "xout1 = np.linspace(-5, -1.9599639845401602)\n",
- "xout2 = np.linspace(1.9599639845401602, 5)\n",
- "plt.plot(x, scipy.stats.norm.pdf(x))\n",
- "plt.fill_between(xout1, scipy.stats.norm.pdf(xout1), np.zeros_like(xout1),color=\"c\")\n",
- "plt.fill_between(xout2, scipy.stats.norm.pdf(xout2), np.zeros_like(xout2),color=\"c\")\n",
- "plt.plot(x, scipy.stats.norm.pdf(x, 3))\n",
- "plt.fill_between(xin, scipy.stats.norm.pdf(xin, 3), np.zeros_like(xin),color=\"orange\")\n",
- "\n",
- "plt.xlabel(\"$x$\")\n",
- "\n",
- "plt.grid()\n",
- "print(\"Fehler 2. Art: beta:\", scipy.stats.norm.cdf(1.9599639845401602, 3) - scipy.stats.norm.cdf(-1.9599639845401602, 3))"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "f6a454b7-fce4-4dca-ad3c-9e220956706b",
+ "id": "7311c0ff-0ce0-4427-a50e-b6698100e454",
"metadata": {
"slideshow": {
- "slide_type": "slide"
+ "slide_type": ""
},
"tags": []
},
"source": [
- "# Neyman-Pearson Test\n",
- "\n",
- "Hypothese: $H$\n",
- "\n",
- "Alternative Hypothese: $A$\n",
+ "### Mit Python I\n",
"\n",
- "Suche Kriterium, das $\\alpha$ und $\\beta$ minimiert. $H$ wird für $x$ im Bereich $V$ verworfen:\n",
+ "Mit scipy.optimize:\n",
"\n",
- "$\\int_{V} P_{H}(x) dx = \\alpha$ (klein)\n",
+ "```\n",
+ "import scipy.optimize as opti\n",
+ "def fitf(x, m , a):\n",
+ " return m*x + a\n",
+ "pfit, Vfit = opti.curve_fit(fitf , xs, ys, \n",
+ " sigma=[sigma_y]*len(ys),absolue_sigma=True)\n",
+ "print(pfit, Vfit)\n",
+ "```\n",
"\n",
- "$\\int_{V} P_{A}(x) dx = 1 - \\beta$ (groß)\n",
- "\n",
- "In $V$ sind die $x$-Werte, für die $\\frac{P_{A}(x)}{P_{H}(x)}$ groß ist.\n",
- "\n",
- "Neyman-Pearson-Kriterium: $\\frac{P_{A}(x)}{P_{H}(x)} > c$"
+ "Vorsicht! Falsche Unsicherheit ohne `absolute_sigma=True`"
]
},
{
"cell_type": "markdown",
- "id": "39da5505-d3e1-4be6-a1b5-ecef10c58bce",
- "metadata": {},
- "source": []
- },
- {
- "cell_type": "markdown",
- "id": "b636ff7e-5ba1-4fc5-bb40-fb29b63ca88c",
+ "id": "a5fec52e-bd3a-4437-bb43-620de44939b2",
"metadata": {
"slideshow": {
- "slide_type": "slide"
+ "slide_type": ""
},
"tags": []
},
"source": [
- "# Endlich: Passt das Modell?\n",
+ "### Mit Python II\n",
"\n",
- "Hypothese: \"Die $k_i$ Tore pro Bundesligaspiel $i$ sind Poisson-verteilt mit einem gemeinsamen $\\mu = <k>$.\"\n",
+ "Mit scipy.optimize:\n",
"\n",
- "Benötigt für den Test eine alternative Hypothese: \"Die Tore pro Bundesligaspiel $k_i$ sind Poisson verteilt mit $\\mu_i = ki$.\"\n",
- "\n",
- "Wie gut unser Modell einer Poissonverteilung für alle Spiele zu den Daten passt, lässt sich durch den Vergleich der log-Likelihood unseres Modells $-lnL(\\mu; \\vec k)$ zur log-Likelihood eines saturierten Modells (je Spiel ein eigener $\\mu$-Parameter mit $\\mu_i = k_i$), also $-ln\\hat L(\\vec k; \\vec k)$, abschätzen.\n",
- "\n",
- "$P_{H} (\\vec k) = \\prod_{i} P(k_i,\\mu)$; $P_{A} (\\vec k) = \\prod_{i} P(k_i, k_i)$\n",
- "\n",
- "Neyman-Pearson: Likelihoodquotient: $\\frac{P_{A}(x)}{P_{H}(x)} = \\frac{\\hat L(\\vec k; \\vec k)}{L(\\mu; \\vec k)} > c$\n",
- "\n",
- "$d = \\ln \\frac{P_{A}(\\vec k)}{P_{H}(\\vec k)} = -\\ln \\frac{P_{H}(\\vec k)}{P_{A}(\\vec k)} = -\\ln L(\\mu; \\vec k)-(-\\ln\\hat L(\\vec k; \\vec k))$\n",
- "\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "75234801-341f-4474-a2bf-93aafe9add77",
- "metadata": {},
- "outputs": [],
- "source": [
- "mu = np.mean(summe)\n",
- "\n",
- "print(\"mu:\", mu)\n",
- "\n",
- "print(\"P(H):\", np.prod(scipy.stats.poisson.pmf(summe,mu)))\n",
- "print(\"P(A):\", np.prod(scipy.stats.poisson.pmf(summe,summe)))\n",
- "print(\"-ln P(H):\", -np.sum(scipy.stats.poisson.logpmf(summe,mu)))\n",
- "print(\"-ln P(A):\", -np.sum(scipy.stats.poisson.logpmf(summe,summe)))\n",
- "print(\"d:\", -np.sum(scipy.stats.poisson.logpmf(summe,mu)) + np.sum(scipy.stats.poisson.logpmf(summe,summe)))"
+ "```\n",
+ "def chi2(x, y, sy, a, m):\n",
+ " my = m * x + a\n",
+ " r = (y - my)/sy\n",
+ " return np.sum(r**2)\n",
+ " \n",
+ "res = opti.minimize( lambda p: chi2(xs, ys, sigma_y, p[1], p[0]),x0=np.zeros(2))\n",
+ "print(res.x, res.hess_inv*2)\n",
+ "```"
]
},
{
"cell_type": "markdown",
- "id": "e840a509-4f5b-48ef-9594-eb1dc0bbdc4a",
+ "id": "84b00301-c0d3-4595-858e-4f784d69c876",
"metadata": {
"slideshow": {
"slide_type": "slide"
@@ -837,76 +563,48 @@
"tags": []
},
"source": [
- "Ist das gut?\n",
+ "### Inverse Hesse-Matrix und $\\chi^2$ \n",
+ "\n",
+ "$\\Delta \\chi2$ und Kovarianz Ellipse um Minimum gemäß Kovarianzmatrix\n",
+ "genau bei $\\Delta \\chi^2 = 1$. \n",
+ "$$1 = \\delta \\chi^2 = (\\vec a -\\hat \\vec a)^T V^{-1} (\\vec a-\\hat \\vec a)$$\n",
+ "Mit\n",
+ "$\\chi^2(\\vec a) = \\chi^2(\\hat \\vec a) + (\\vec a -\\hat \\vec a)^T V^{-1} (\\vec a-\\hat \\vec a)$\n",
+ "und\n",
+ "$H_{ij} = \\frac{\\partial^2 \\chi^2(\\vec a)}{\\partial a_i \\partial a_j}$ \n",
+ "$$H_{ij} = \\frac{\\partial^2 (a_k -\\hat a_k) V^{-1}_{kl} (a_l -\\hat a_l)}{\\partial a_i \\partial a_j} = \\frac{\\partial( \\delta_{ik}V^{-1}_{kl} (a_l -\\hat a_l) + (a_k -\\hat a_k) V^{-1}_{kl} \\delta_{il})}{\\partial a_j}$$\n",
+ "$$H_{ij} = \\delta_{ik}V^{-1}_{kl}\\delta_{lj} + \\delta_{jk}V^{-1}_{kl}\\delta_{il} = 2V^{-1}_{ij} \\text{ und } V_{ij} = 2 * H^{-1}_{ij}$$\n",
"\n",
- "Wie sieht die Verteilung von d aus, wenn das Modell stimmt und $\\mu=<k>$?"
+ "Vorsicht! Manche Algorithmen in `minimize` berechnen keine inverse\n",
+ "Hesse-Matrix."
]
},
{
- "cell_type": "code",
- "execution_count": null,
- "id": "90e27f87-36fd-4b11-add1-c82a4eeb79dd",
+ "cell_type": "markdown",
+ "id": "632583d4-1e97-498b-b8f3-91e259baa24d",
"metadata": {
- "jupyter": {
- "source_hidden": true
- },
"slideshow": {
- "slide_type": ""
+ "slide_type": "slide"
},
"tags": []
},
- "outputs": [],
"source": [
- "#simuliere 1000 Spielzeiten\n",
- "muobs = np.mean(summe)\n",
+ "# Maximum-Likelihood \n",
"\n",
- "tore = scipy.stats.poisson.rvs(muobs,size=(1000,306))\n",
+ "Maximum-Likelihood (ML) Daten: $x_1,...,x_N$ \n",
+ "Wahrscheinlichkeit der Daten für Modell mit Parametern $a$:\n",
+ "$$P(x_1,...,x_N; a) = \\prod_i P(x_i ; a)$$\n",
"\n",
- "d = np.zeros(len(tore))\n",
- "mu = np.zeros(len(tore))\n",
- "for i in range(len(tore)):\n",
- " mu[i] = np.mean(tore[i,:])\n",
- " d[i] = np.sum(-scipy.stats.poisson.logpmf(tore[i],mu[i]) + scipy.stats.poisson.logpmf(tore[i],tore[i]))\n",
- " \n",
- "plt.hist(mu, bins=50)\n",
- "muobs = np.mean(summe)\n",
- "plt.plot([muobs, muobs], [0, 100], linestyle = 'dotted')\n",
- "plt.grid()\n",
- "plt.xlabel(\"$\\hat \\mu$\")\n",
- "plt.show()\n",
- "dobs = -np.sum(scipy.stats.poisson.logpmf(summe,muobs)) + np.sum(scipy.stats.poisson.logpmf(summe,summe))\n",
- "plt.hist(d, bins=50)\n",
- "plt.plot([dobs, dobs], [0, 100], linestyle = 'dotted')\n",
- "\n",
- "plt.grid()\n",
- "plt.xlabel(\"$-\\ln(P(H)/P(A))$\")\n",
- "plt.show()\n",
- "\n"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "f6db4bd4-de67-48e1-add9-4fdc05dbbcf7",
- "metadata": {},
- "source": [
- "$p$-Wert \n",
+ "Likelihoodfunktion: $$L(a) = \\prod_i P(x_i ; a)$$\n",
"\n",
- "Anteil der erwarteter Werte, die höher als der Daten-Wert sind."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "6647f14e-2282-4080-9139-c332db4af23a",
- "metadata": {},
- "outputs": [],
- "source": [
- "print(\"p:\", np.sum(d > dobs)/len(d))"
+ "ML-Schätzer $\\hat a$: Maximum von $L(a)$:\n",
+ "$$\\left.\\frac{dL}{da}\\right|_{a = \\hat a} = 0$$ (praktischer:\n",
+ "Log-Likelihood: $-\\ln L = \\sum_i -\\ln P(x_i; a)$)"
]
},
{
"cell_type": "markdown",
- "id": "c835358c-1291-431c-b0ab-792a4b940d68",
+ "id": "dd0afbf7-8504-46f8-b3da-4fb33487a635",
"metadata": {
"slideshow": {
"slide_type": "slide"
@@ -914,177 +612,91 @@
"tags": []
},
"source": [
- "## Geht es auch ohne MC?\n",
+ "### Beispiel\n",
"\n",
- "Einfacherer Fall: Gauß-Verteilungen:\n",
"\n",
- "$d = -\\ln L(\\mu; \\vec x)-(-\\ln\\hat L(\\vec x; \\vec x) = -\\ln \\prod_{i} \\frac{G(x_i,\\mu, \\sigma)}{G(x_i, x_i, \\sigma)}=-\\ln \\prod_{i} \\frac{\\exp(-\\frac{1}{2}(\\frac{x_i-\\mu}{\\sigma})^2)}{\\exp(-\\frac{1}{2}(\\frac{x_i-x_i}{\\sigma})^2)} = -\\ln \\exp(-\\frac{1}{2}(\\frac{x_i-\\mu}{\\sigma})^2) = \\frac{1}{2} \\sum_i\\left(\\frac{x_i - \\mu}{\\sigma}\\right)^2$\n",
+ "$y(x) = mx + a$: Finde $\\hat m$ und $\\hat a$ Daten: $y_1,...,y_N$ und\n",
+ "Modell: $$P(y_i; m, a) = G(y_i; \\mu = m x_i + a, \\sigma=\\sigma_i)$$\n",
+ "$$L(m, a) = \\prod_i G(y_i; \\mu = m x_i + a, \\sigma=\\sigma_i)$$\n",
"\n",
- "$\\chi^2 = -2 \\ln\\frac{P(H}{P(A)}$ (Wilks Theorem)\n",
+ "0.5 <img src=\"./figures/11/line.png\" alt=\"image\" />\n",
"\n",
- "$-2 \\ln\\frac{P(H}{P(A)}$ sollte gemäß $\\chi^2$-Verteilung verteilt sein, mit n Freiheitsgraden n = Zahl der Messpunkte - Zahl der Modellparameter"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "9d667a3b-6a89-419a-8be4-79cb7de96151",
- "metadata": {},
- "outputs": [],
- "source": [
- "plt.hist(2*d, bins=50, density=True)\n",
- "plt.plot([2*dobs, 2*dobs], [0, 0.02], linestyle = 'dotted')\n",
- "ds = np.linspace(200, 425, 100)\n",
- "plt.plot(ds,scipy.stats.chi2.pdf(ds, 305))\n",
+ "<img src=\"./figures/11/like_a.png\" style=\"width:49.0%\"\n",
+ "alt=\"image\" />\n",
+ "<img src=\"./figures/11/loglike_a.png\" style=\"width:49.0%\"\n",
+ "alt=\"image\" />\n",
"\n",
- "plt.grid()\n",
- "plt.xlabel(\"$-2\\ln(P(H)/P(A))$\")\n",
- "plt.show()\n",
- "\n",
- "\n",
- "print(\"p-Wert über Chi2:\", scipy.stats.chi2.sf(2*dobs, 305))"
+ "ML-Schätzer für Poisson $\\mu$ $$\\begin{aligned}\n",
+ " L(\\mu) & = & \\prod_i^N P(k_i; \\mu) = \\prod_i^N \\frac{\\mu^{k_i}e^{-\\mu}}{k_i!}\\\\\n",
+ " \\ln L(\\mu) & = & \\sum_{i=1}^N \\left( \\ln \\mu^{k_i} + \\ln e^{-\\mu} - \\ln k_i!\\right)\\\\\n",
+ " & = & \\sum_{i=1}^N \\left( k_i \\ln \\mu -\\mu - \\ln k_i!\\right)\\\\\n",
+ " 0 \\stackrel{!}{=} \\frac{d \\ln L(\\mu)}{d\\mu} \\Big|_{\\hat \\mu}& = & \\sum_{i=1}^N \\left( \\frac{k_i}{\\hat \\mu} - 1\\right) = \\sum_{i=1}^N \\frac{k_i}{\\hat\\mu} - N\\\\\n",
+ " N & = & \\frac{1}{\\hat\\mu} \\sum_{i=1}^N k_i \\rightarrow \\hat\\mu = \\frac{1} {N} \\sum_{i=1}^N k_i\n",
+ " \n",
+ "\\end{aligned}$$"
]
},
{
"cell_type": "markdown",
- "id": "eb7b32c7-a9a7-44de-aa70-ce9744218a35",
- "metadata": {},
- "source": [
- "Stimmt die Gaußsche Näherung in unserem Fall?"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "cda0957a-414f-4bea-ab33-e61218a4a4f5",
- "metadata": {},
- "outputs": [],
- "source": [
- "plt.hist( scipy.stats.chi2.sf(2*d, 305), bins=50)\n",
- "plt.plot([ scipy.stats.chi2.sf(2*dobs, 305), scipy.stats.chi2.sf(2*dobs, 305)], [0, 200], linestyle = 'dotted')\n",
- "\n",
- "plt.grid()\n",
- "plt.xlabel(\"p\")\n",
- "plt.show()"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "557559e1-8e4a-412a-9486-a2fb9b97de90",
- "metadata": {},
- "source": [
- "Für große $\\mu$? Handball??"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "472aac93-e2dd-424a-ba0e-51c10d409f8e",
+ "id": "d9571970-772e-4e95-8474-82e0afb1dd68",
"metadata": {
"slideshow": {
"slide_type": ""
},
"tags": []
},
- "outputs": [],
"source": [
- "#simuliere 5000 Spielzeiten\n",
- "muobs = 35\n",
- "\n",
- "tore = scipy.stats.poisson.rvs(muobs,size=(5000,306))\n",
+ "Varianz des ML-Schätzers\n",
"\n",
- "d = np.zeros(len(tore))\n",
- "mu = np.zeros(len(tore))\n",
- "for i in range(len(tore)):\n",
- " mu[i] = np.mean(tore[i,:])\n",
- " d[i] = np.sum(-scipy.stats.poisson.logpmf(tore[i],mu[i]) + scipy.stats.poisson.logpmf(tore[i],tore[i]))\n",
- " \n",
- "plt.hist(mu, bins=50)\n",
- "plt.grid()\n",
- "plt.xlabel(\"$\\hat \\mu$\")\n",
- "plt.show()\n",
- "\n",
- "plt.hist(2 * d, bins=50, density=True)\n",
- "ds = np.linspace(200, 425, 100)\n",
- "plt.plot(ds,scipy.stats.chi2.pdf(ds, 305))\n",
- "plt.grid()\n",
- "plt.xlabel(\"$-2\\ln(P(H)/P(A))$\")\n",
- "plt.show()\n",
- "\n",
- "plt.hist( scipy.stats.chi2.sf(2*d, 305), bins=50)\n",
- "plt.grid()\n",
- "plt.xlabel(\"p\")\n",
- "plt.show()"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "450edccb-2c13-4400-9613-4a1f3e069ce6",
- "metadata": {
- "slideshow": {
- "slide_type": "slide"
- },
- "tags": []
- },
- "source": [
- "# $\\chi^2$-Test\n",
+ "Rao-Cramér-Frechet-Ungleichung: Schätzer $\\hat a$ mit Bias (Verzerrung)\n",
+ "$b$\n",
+ "$$V(\\hat a) \\geq \\frac{\\left(1+ \\frac{\\partial b}{\\partial a} \\right)^2}{E\\left[-\\frac{\\partial^2 \\ln L}{\\partial a^2}\\right]}$$\n",
+ "Fisher-Information:\n",
+ "$$I(\\hat a) = E\\left [-\\frac{\\partial^2 \\ln L}{\\partial a^2}\\right]$$\n",
"\n",
+ "ML-Schätzer für Poisson $V(\\hat \\mu)$ $$\\begin{aligned}\n",
+ "V(\\hat \\mu) & \\geq &\\frac{\\left(1+ \\frac{\\partial b}{\\partial \\mu} \\right)^2}{E\\left[-\\frac{\\partial^2 \\ln L}{\\mu^2}\\right]} \\\\\n",
+ " & = & \\frac{1}{E\\left[-\\frac{\\partial(\\sum_{i=1}^N \\frac{k_i}{\\mu} - N)}{\\partial \\mu^2}\\right]} \\\\\n",
+ " & = & \\frac{1}{E\\left[-\\sum_{i=1}^N \\frac{-k_i}{\\hat \\mu^2}\\right]} = \\frac{1}{E\\left[\\sum_{i=1}^N \\frac{k_i}{\\hat \\mu^2}\\right]} \\\\\n",
+ " & = & \\frac{1}{\\frac{1}{\\hat \\mu^2}E\\left[\\sum_{i=1}^N k_i \\right]} = \\frac{1}{\\frac{1}{\\hat \\mu^2}E\\left[N \\hat \\mu \\right]}\\\\\n",
+ " & = & \\frac{\\hat \\mu}{N}\n",
+ "\\end{aligned}$$\n",
"\n",
- "Haben wir schon die ganze Zeit gemacht...\n",
+ "Varianz für mehrere Parameter $\\vec \\theta$\n",
"\n",
- "Wie sind die Konfidenzintervalle $n$ Freiheitsgrade/Dimensionen?\n",
+ "Für effizienten und erwartungstreuen Schätzer:\n",
+ "$$\\left(V^{-1}\\right)_{ij} = E\\left[ -\\frac{\\partial^2 \\ln L(\\theta)}{\\partial \\theta_i \\partial \\theta_j}\\right]$$\n",
"\n",
- "$$\\chi^2 = (\\overrightarrow{y\\ } - \\vec{f})^{T} V(\\vec{y}- \\vec{f})$$\n",
- "\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "2fd5c302-732b-4c6b-b224-5b9f79df57cc",
- "metadata": {
- "slideshow": {
- "slide_type": ""
- },
- "tags": []
- },
- "outputs": [],
- "source": [
- "def conv_chi2_nd(z, n):\n",
- " return scipy.stats.chi2.cdf(z,n)\n",
+ "Näherung für große Datensätze:\n",
+ "$$\\left(\\hat V^{-1}\\right)_{ij} = -\\frac{\\partial^2 \\ln L(\\theta)}{\\partial \\theta_i \\partial \\theta_j}\\Big|_{\\theta=\\hat \\theta} =$$\n",
"\n",
+ "Graphisch:\n",
+ "$$\\ln L(\\theta) \\approx \\ln L(\\hat \\theta) + \\frac{\\partial \\ln L}{\\partial \\theta}\\Big|_{\\hat \\theta}(\\theta - \\hat \\theta) + \\frac{1}{2} \\frac{\\partial^2 \\ln L}{\\partial \\theta^2}(\\theta - \\hat \\theta)^2$$\n",
+ "$$\\ln L(\\hat \\theta + \\sigma_\\theta) \\approx \\ln L(\\hat \\theta) + \\frac{1}{2} \\frac{\\partial^2 \\ln L}{\\partial \\theta^2}(\\sigma_\\theta)^2 = \\ln L(\\hat \\theta) - \\frac{1}{2}$$\n",
"\n",
- "for z in [1,2,3,4,5]:\n",
- " print(z, conv_chi2_nd(z, 1), conv_chi2_nd(z, 2), conv_chi2_nd(z, 3))\n",
+ "Zusammenhang ML und $\\chi^2$\n",
"\n",
+ "Likelihood-Quotient:\n",
+ "$$\\lambda = -2 \\ln \\frac{L(\\hat \\theta)}{L(\\hat \\theta^\\prime_\\text{saturiert})}$$\n",
"\n",
- "print(\"Kritische chi2-Werte\")\n",
- "for k in [1, 2, 3]:\n",
- " c = conv_gaus_nd(k, 1)\n",
- " print(c, scipy.optimize.brentq(lambda z: conv_chi2_nd(z, 1)-c,0, 30), scipy.optimize.brentq(lambda z: conv_chi2_nd(z, 2)-c,0, 30), scipy.optimize.brentq(lambda z: conv_chi2_nd(z,3)-c,0, 30))"
+ "Mit Normalverteilung: $$\\begin{aligned}\n",
+ "\\lambda &=& -2 \\ln \\frac{L(\\hat \\theta)}{L(\\hat \\theta^\\prime_\\text{saturiert})} = -2 \\ln \\frac{\\prod_i G(x_i; \\hat \\mu, \\sigma_i)}{\\prod_i G(x_i; x_i, \\sigma_i)}\\\\\n",
+ "& = & -2 \\ln \\frac{\\frac{1}{\\sqrt{2\\pi}\\sigma_i}exp\\left(\\frac{(x_i-\\hat \\mu)^2}{2\\sigma_i^2}\\right)}{\\frac{1}{\\sqrt{2\\pi}\\sigma_i}exp\\left(\\frac{(x_i-x_i)^2}{2\\sigma_i^2}\\right)} = -2\\ln exp\\left(\\frac{(x_i-\\hat \\mu)^2}{2\\sigma_i^2}\\right) \\\\\n",
+ "& = & -2 \\frac{(x_i-\\hat \\mu)^2}{2\\sigma_i^2} = \\chi^2 \\text{; also } \\ln L(\\theta) = - \\chi^2(\\theta) / 2 \n",
+ "\\end{aligned}$$"
]
},
{
"cell_type": "code",
"execution_count": null,
- "id": "23626cec-4198-44f4-a1c8-81e36992268f",
- "metadata": {},
- "outputs": [],
- "source": [
- "for c in [0.68, 0.90, 0.95, 0.99]:\n",
- " print(c, scipy.optimize.brentq(lambda z: conv_chi2_nd(z, 1)-c,0, 30), scipy.optimize.brentq(lambda z: conv_chi2_nd(z, 2)-c,0, 30),scipy.optimize.brentq(lambda z: conv_chi2_nd(z,3)-c,0, 30))"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "d79940bc-7c64-42c8-a564-898995f51ddc",
+ "id": "a73415de",
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
+ "celltoolbar": "Slideshow",
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
@@ -1101,6 +713,14 @@
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.9.20"
+ },
+ "livereveal": {
+ "autolaunch": true,
+ "overlay": "<div class='myfooter'><h2 style='text-align:center;'>Hartmut Stadie - University of Hamburg</h2></div>",
+ "scroll": true
+ },
+ "toc": {
+ "base_numbering": 43
}
},
"nbformat": 4,
--
GitLab