{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Spiegelhalter's Z-test\n",
    "\n",
    "## Theoretical background\n",
    "\n",
    "Spiegelhalter's Z-test is a statistical test that tests whether a probabilistic model is calibrated. It is named after the statistician David Spiegelhalter, who proposed it in 1986. It is a non-parametric test that does not require any binning.\n",
    "\n",
    "The Spiegelhalter's Z-test was inspired by the fact that the Brier score (mean squared error) can be decomposed into reliability and resolution. In fact, any proper scoring rule can be decomposed into reliability and resolution, as shown by Brocker (2008). For example, the cross-entropy can be decomposed into KL-divergence (reliability) and entropy (resolution).\n",
    "\n",
    "The Brier score can be decomposed into reliability and resolution as follows:\n",
    "$$\n",
    "B = \\frac{1}{N} \\sum_{i=1}^N (x_i - p_i)^2 = \\frac{1}{N} \\sum_{i=1}^N (x_i - p_i)(1-2p_i) + \\frac{1}{N} \\sum_{i=1}^N p_i(1-p_i)\n",
    "$$\n",
    "Where the first term measure the reliability/calibration and the second term measure the resolution/discrimination.\n",
    "\n",
    "The Variance of the Brier score is:\n",
    "$$\n",
    "\\text{Var}(B) = \\frac{1}{N^2} \\sum_{i=1}^N (1-2p_i)^2 p_i (1-p_i)\n",
    "$$\n",
    "\n",
    "and the Speigelhalter's Z-test is defined as:\n",
    "$$\n",
    "Z = \\frac{B - E(B)}{\\sqrt{\\text{Var}(B)}} = \\frac{ \\sum_{i=1}^N (x_i - p_i)(1-2p_i)}{\\sum_{i=1}^N (1-2p_i)^2 p_i (1-p_i)}\n",
    "$$\n",
    "\n",
    "and $Z$ is approximately standard normal distributed under the null hypothesis of calibration. Spiegelhalter's Z-test has the right size in many situations and it is powerful in many situations. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pros of Spiegelhalter's Z test\n",
    "\n",
    "Spiegelhalter's Z test is a statistical test which can provide statistical evidence that the null hypothesis (well-calibrated) is true or false. It is a non-parametric test and doesn't require any hyperparameter tuning. It also doesn't require any binning of data, which is extremely useful compared to the Hosmer-Lemeshow test.\n",
    "\n",
    "## Cons of Spiegelhalter's Z test\n",
    "\n",
    "The power of Spiegelhalter's Z test is limited for some cases of miscalibration, such as prevalence shift. However, it is a very powerful test for many other cases of miscalibration."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculating the Spieegelhalter Z score and p-value using calzone\n",
    "\n",
    "We can call functions from the calzone package to calculate the Spiegelhalter Z score and p-value directly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Z-score: 0.3763269161877356, p-value: 0.7066738713391099\n"
     ]
    }
   ],
   "source": [
    "from calzone.utils import reliability_diagram,data_loader\n",
    "from calzone.metrics import spiegelhalter_z_test\n",
    "import numpy as np\n",
    "\n",
    "wellcal_dataloader = data_loader(data_path=\"../../../example_data/simulated_welldata.csv\")\n",
    "z,p_value = spiegelhalter_z_test(wellcal_dataloader.labels,wellcal_dataloader.probs,class_to_calculate=1)\n",
    "print(f\"Z-score: {z}, p-value: {p_value}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also use the CalibrationMetrics class"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'SpiegelhalterZ score': 0.3763269161877356,\n",
       " 'SpiegelhalterZ p-value': 0.7066738713391099}"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from calzone.metrics import CalibrationMetrics\n",
    "calmetrics = CalibrationMetrics(class_to_calculate=1)\n",
    "calmetrics.calculate_metrics(wellcal_dataloader.labels, wellcal_dataloader.probs, metrics=['SpiegelhalterZ'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## P-value of Spiegelhalter's z test\n",
    "\n",
    "Like to HL test, we can check whether the Spiegelhalter's z test has the correct size."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "### The size of HL Test\n",
    "from calzone.utils import fake_binary_data_generator\n",
    "np.random.seed(123)\n",
    "fakedata_generator = fake_binary_data_generator(alpha_val=0.5, beta_val=0.5)\n",
    "cal_metrics = CalibrationMetrics()\n",
    "sample_size = 1000\n",
    "simulation_size = 10000\n",
    "results = []\n",
    "# generate data\n",
    "for i in range(simulation_size):\n",
    "    X, y = fakedata_generator.generate_data(sample_size)\n",
    "    if i == 0:\n",
    "        tempresult = cal_metrics.calculate_metrics(y, X, ['SpiegelhalterZ'],return_numpy=False)\n",
    "        keys = list(tempresult.keys())\n",
    "        results.append(np.array(list(tempresult.values())))\n",
    "    else:\n",
    "        tempresult = cal_metrics.calculate_metrics(y, X, ['SpiegelhalterZ'],return_numpy=True)\n",
    "        results.append(tempresult)\n",
    "results = np.array(results)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The size of Spiegelhalter's z test is : 0.049\n"
     ]
    }
   ],
   "source": [
    "### Showing the size of the model\n",
    "import matplotlib.pyplot as plt\n",
    "z_scores = results[:,0]\n",
    "p_values = results[:,1]\n",
    "size = np.mean(p_values < 0.05)\n",
    "print(\"The size of Spiegelhalter's z test is :\", round(size,3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, \"P-value distribution of the Spiegelhalter's z test\")"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAHFCAYAAAAOmtghAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABWJElEQVR4nO3deVwU9f8H8NeysNyggJwi4I23Yiqa4YkHmlYmlomomX7NSpEss7yV8ko70CwFLQ/yrJRUMs888qpM7PDEA0Q0uZRr9/P7gx+r6y7XCnuMr+fjMY/H7Gc/M/OeYXfnzefzmRmZEEKAiIiISCIsjB0AERERUVVickNERESSwuSGiIiIJIXJDREREUkKkxsiIiKSFCY3REREJClMboiIiEhSmNwQERGRpDC5ISIiIklhcmPm4uPjIZPJ1JOlpSVq166NESNG4Pr16waPx9/fH5GRkQbfbmn27dsHmUyGffv2qcsiIyPh7+9fqfXcuHEDM2bMwG+//Vap5XRtSyaTYfz48ZVaT3liY2MRHx+vVX758mXIZDKd75mShIQENG3aFLa2tpDJZKUe5+TkZMyYMQOXL1/Weq9Lly5o1qxZ9QYKQAiBDRs2oHPnznB3d4eNjQ1q166NXr164auvvtJrnSXfY137ZQyG+Nzq+mwePnwYM2bMwN27d6tsOw+LjIxEly5dqmXdFTVv3jxs27atWrdR1vfkScHkRiLi4uJw5MgRJCUlYfTo0Vi/fj06d+6M3NxcY4dmcj744ANs3bq1UsvcuHEDM2fOrHRyo8+29FFacuPl5YUjR44gLCys2mPQ161btzBs2DDUq1cPO3fuxJEjR9CwYUOddZOTkzFz5kyj/mhPmTIFL730EgIDA/HVV1/hxx9/xJw5c+Dh4YHvvvtOr3WGhYXhyJEj8PLyquJozcvhw4cxc+bMaktuTIGhkhtjf0+MzdLYAVDVaNasGdq2bQsA6Nq1K5RKJWbPno1t27Zh6NChRo7OtNSrV6/at3Hv3j3Y2dkZZFtlsba2RocOHYwaQ3n++ecfFBYW4pVXXkFISIixwynT/fv3sWTJEkRERGDFihUa70VGRkKlUum13lq1aqFWrVpVESLpUPJ9pCcHW24kquSEduXKFZ3vFxYWwt3dHcOGDdN67+7du7C1tUVUVBQAIC8vD5MmTUKrVq3g7OwMFxcXBAcHV+i/1NKa23V1FwHATz/9hO7du8PJyQl2dnbo1KkT9uzZU4E9Bv766y/07t0bdnZ2cHNzw9ixY5Gdna1VT1eT+8aNG9G+fXs4OzvDzs4OdevWxciRI9WxPvXUUwCAESNGqLsAZ8yYoV6fg4MDzpw5g9DQUDg6OqJ79+6lbqvEF198gYYNG8La2hpNmjTBhg0bNN6fMWMGZDKZ1nKPHlN/f3+cPXsW+/fvV8dWss3SuqUOHTqE7t27w9HREXZ2dujYsSN27Nihczt79+7F//73P7i5ucHV1RXPP/88bty4oXOfHvX9998jODgYdnZ2cHR0RM+ePXHkyBH1+5GRkXj66acBAOHh4ZDJZKV2G8THx+PFF18EUJzAl+zro/t2/PhxdO7cWf13/PDDD7WSjqysLERHRyMgIAAKhQI+Pj6YMGFCuS2dubm5yM/PL7WFxcLiwU9qybGfP38+5s6dizp16sDGxgZt27bV+kyX9j2p6Pfhu+++Q4sWLWBtbY26deti6dKlOj8/QgjExsaiVatWsLW1Rc2aNTFo0CBcvHixzP1+2Ndff43AwEDY2dmhZcuW2L59u8b758+fx4gRI9CgQQPY2dnBx8cH/fv3x5kzZ8pc74wZM/D2228DAAICAtR/34d/IxISEhAcHAx7e3s4ODigV69eOH36tMZ6yvo+6rJs2TK0bNkSDg4OcHR0ROPGjfHee++VGWtkZKTGcICHp5LfBV1kMhlyc3OxevVqdf2HP+9paWkYM2YMateuDYVCgYCAAMycORNFRUUVjrmi3xPJE2TW4uLiBABx/PhxjfKlS5cKAGLFihWlLjtx4kRha2srMjMzNcpjY2MFAPHHH38IIYS4e/euiIyMFF9//bX4+eefxc6dO0V0dLSwsLAQq1ev1ljWz89PDB8+XCu+S5cuadTbu3evACD27t2rLvv666+FTCYTAwcOFFu2bBE//PCD6Nevn5DL5eKnn34q8zikpaUJd3d34ePjI+Li4kRiYqIYOnSoqFOnjtZ2hg8fLvz8/NSvDx8+LGQymRgyZIhITEwUP//8s4iLixPDhg0TQgiRmZmp3o/3339fHDlyRBw5ckRcvXpVvT4rKyvh7+8vYmJixJ49e8SuXbt0bksIIQAIX19f0aRJE7F+/Xrx/fffi969ewsAYuPGjep606dPF7q+oo8e01OnTom6deuK1q1bq2M7deqUEEKIS5cuCQAiLi5Ovfy+ffuElZWVCAoKEgkJCWLbtm0iNDRUyGQysWHDBq3t1K1bV7zxxhti165d4quvvhI1a9YUXbt2LfPvIYQQa9euFQBEaGio2LZtm0hISBBBQUFCoVCIgwcPCiGEOH/+vPj8888FADFv3jxx5MgRcfbsWZ3rS09PF/PmzRMAxOeff67e1/T0dCGEECEhIcLV1VU0aNBALF++XCQlJYlx48YJABqf09zcXNGqVSvh5uYmFi9eLH766SexdOlS4ezsLLp16yZUKlWZ+1W/fn3h6OgoFi1aJM6dO1dq/ZJj7+vrK55++mmxefNmsXHjRvHUU08JKysrcfjwYa1j/fD3pKLfhx9//FFYWFiILl26iK1bt4qNGzeK9u3bC39/f63Pz+jRo4WVlZWYNGmS2Llzp1i3bp1o3Lix8PDwEGlpaep6pX1u/f39Rbt27cS3334rEhMTRZcuXYSlpaW4cOGCut7+/fvFpEmTxKZNm8T+/fvF1q1bxcCBA4Wtra3466+/tI5PyWfz6tWr4o033hAAxJYtW9R/35Lfp7lz5wqZTCZGjhwptm/fLrZs2SKCg4OFvb29xmemrO/jo9avXy8AiDfeeEPs3r1b/PTTT2L58uXizTff1Fm/xPnz59XxlUyvvPKKACASEhJKXe7IkSPC1tZW9O3bV71cSeypqanC19dX+Pn5iS+++EL89NNPYvbs2cLa2lpERkZWOObyvidPCiY3Zq7kR/Ho0aOisLBQZGdni+3bt4tatWoJR0dHjR+sR/3xxx86E6B27dqJoKCgUpcrKioShYWFYtSoUaJ169Ya7+mb3OTm5goXFxfRv39/jXpKpVK0bNlStGvXroyjIMQ777wjZDKZ+O233zTKe/bsWW5ys3DhQgFA3L17t9T1Hz9+XCtJeHh9AMSqVat0vqfrJGFra6vxtykqKhKNGzcW9evXV5dVNLkRQoimTZuKkJAQrbq6kpsOHToId3d3kZ2drbH9Zs2aidq1a6tP1iXbGTdunMY658+fLwCI1NRUre2VUCqVwtvbWzRv3lwolUp1eXZ2tnB3dxcdO3ZUl5V8Fh5O7EqzceNGrb9niZCQEAFAHDt2TKO8SZMmolevXurXMTExwsLCQusfgk2bNgkAIjExscwYfv31V3XSDEA4OjqKfv36iTVr1mgkOiXH3tvbW9y/f19dnpWVJVxcXESPHj3UZY/+TSvzfXjqqaeEr6+vyM/PV5dlZ2cLV1dXjc/PkSNHBACxaNEijXVevXpV2NraismTJ6vLSvvcenh4iKysLHVZWlqasLCwEDExMaUer6KiIlFQUCAaNGggJk6cqHV8Hv5sLliwQOfvRUpKirC0tBRvvPGGRnl2drbw9PQUgwcP1oi9tO/jo8aPHy9q1KhRbr3yfPvtt0Imk4n33nuv3Lr29vYav5ElxowZIxwcHMSVK1c0ykt+n0qSoIrEXNb35EnBbimJ6NChA6ysrODo6Ih+/frB09MTP/74Izw8PCCEQFFRkcYEAM2bN0dQUBDi4uLU6zl37hx+/fVXdZdMiY0bN6JTp05wcHCApaUlrKyssHLlSpw7d65K4j98+DDu3LmD4cOHa8SpUqnQu3dvHD9+vMwug71796Jp06Zo2bKlRvnLL79c7rZLupwGDx6Mb7/9Vu+rzF544YUK1+3evTs8PDzUr+VyOcLDw3H+/Hlcu3ZNr+1XRG5uLo4dO4ZBgwbBwcFBY/vDhg3DtWvX8Pfff2ss8+yzz2q8btGiBYDSuzwB4O+//8aNGzcwbNgwja4aBwcHvPDCCzh69Cju3btXFbukwdPTE+3atdOK9+FYt2/fjmbNmqFVq1Yan7VevXrp7Cp91FNPPYXz589j586deO+99xAcHIw9e/YgIiICzz77LIQQGvWff/552NjYqF87Ojqif//+OHDgAJRKpc5tVPT7kJubixMnTmDgwIFQKBTq5R0cHNC/f3+NdW7fvh0ymQyvvPKKxjo9PT3RsmXLcvcbKO7mcHR0VL/28PCAu7u7xvEtKirCvHnz0KRJEygUClhaWkKhUODff//V+/di165dKCoqQkREhEbsNjY2CAkJ0Rl7Rb6P7dq1w927d/HSSy/hu+++Q0ZGRqVj279/P4YNG4ZXXnkFc+fOrfTyJbZv346uXbvC29tbYx/79Omj3k5Vxfwk4IBiiVizZg0CAwNhaWkJDw8PjTEBq1evxogRIzTql/wAjxw5Eq+//jr++usvNG7cGHFxcbC2tsZLL72krrtlyxYMHjwYL774It5++214enrC0tISy5Ytw6pVq6ok/ps3bwIABg0aVGqdO3fuwN7eXud7t2/fRkBAgFa5p6dnudt+5plnsG3bNnzyySeIiIhAfn4+mjZtiqlTp2och7LY2dnBycmpQnVLi6uk7Pbt26hdu3aF11UZ//33H4QQOseMeHt7q7f/MFdXV43X1tbWAIoH15amZB2lbUelUuG///6r8kGej8YKFMf7cKw3b97E+fPnYWVlpXMdFTlZWFlZoVevXujVqxeA4v0dNGgQtm/fjh9//BF9+/ZV1y3tb11QUICcnBw4OztrvV/R74NMJoMQQiNRLvFo2c2bN0utCwB169YtdVslKnJ8o6Ki8Pnnn+Odd95BSEgIatasCQsLC7z66qtlfmbKUnI8Sv4RedTDCTRQ8e/jsGHDUFRUhC+//BIvvPACVCoVnnrqKcyZMwc9e/Ysd/mzZ89i4MCB6Ny5M1auXFmBPSndzZs38cMPP5T7uXzcmJ8UTG4kIjAwUH211KP69++P48eP63zvpZdeQlRUFOLj4zF37lx8/fXXGDhwIGrWrKmu88033yAgIAAJCQkaAxTz8/PLjavkP9ZH6z56AnFzcwMAfPrpp6Ve3VPajzJQ/KOblpamVa6rTJcBAwZgwIAByM/Px9GjRxETE4OXX34Z/v7+CA4OLnd5XQN/y1JWrCUnkIePXUlCAVTs5FuakhNNamqq1nslg4RL/haPo2QfStuOhYWFxmfMkNzc3GBra1tqYq7P/ru6umLChAnYt28f/vzzT43kprS/tUKh0Gg90xVDed+HwsJCyGQy9cn/0W08uk6ZTIaDBw9qfJ5K6CrTxzfffIOIiAjMmzdPozwjIwM1atTQa50lx2PTpk3w8/Mrt35lvo8jRozAiBEjkJubiwMHDmD69Ono168f/vnnnzK3de3aNfTu3Rt16tTB5s2bS01KKsrNzQ0tWrQotfWn5J+Px4n5ScLk5gng6uqq8z8uoPhkN3DgQKxZswbBwcFIS0vT6pKSyWRQKBQaPxhpaWkVulqq5KqdP/74A40aNVKXf//99xr1OnXqhBo1aiA5OVmvG4V17doV8+fPx++//67RNbVu3bpKrcfa2hohISGoUaMGdu3ahdOnTyM4OLhCrRWVsWfPHty8eVOdsCmVSiQkJKBevXrqVpuHj93D/7H+8MMPOuOuSGz29vZo3749tmzZgoULF8LW1hYAoFKp8M0336B27dql3mOmMho1agQfHx+sW7cO0dHR6s9Obm4uNm/erL6CqrKq4u/Qr18/zJs3D66urjpb+8pSWFiIrKwsnd+nki6Xh09CQHHL54IFC9TJanZ2Nn744Qd07twZcrlc53Yq+n1QKBRo27Yttm3bhoULF6q7pnJycrSuYurXrx8+/PBDXL9+HYMHD674TleSTCbTSpR27NiB69evo379+mUuW9rft1evXrC0tMSFCxcq1f1bGfb29ujTpw8KCgowcOBAnD17ttREITMzE3369IFMJkNiYmKlWm1L+67269cPiYmJqFevXoUT/9JirurfK3PE5IYwcuRIJCQkYPz48ahduzZ69Oih8X6/fv2wZcsWjBs3DoMGDcLVq1cxe/ZseHl54d9//y1z3U899RQaNWqE6OhoFBUVoWbNmti6dSsOHTqkUc/BwQGffvophg8fjjt37mDQoEFwd3fHrVu38Pvvv+PWrVtYtmxZqduZMGECVq1ahbCwMPUN1dauXYu//vqr3P2fNm0arl27hu7du6N27dq4e/culi5dCisrK/V9V+rVqwdbW1usXbsWgYGBcHBwgLe3t9aJrKLc3NzQrVs3fPDBB7C3t0dsbCz++usvjcvB+/btCxcXF4waNQqzZs2CpaUl4uPjcfXqVa31NW/eHBs2bEBCQgLq1q0LGxsbNG/eXOe2Y2Ji0LNnT3Tt2hXR0dFQKBSIjY3Fn3/+ifXr11e6FUoXCwsLzJ8/H0OHDkW/fv0wZswY5OfnY8GCBbh79y4+/PBDvdZbcgfiFStWwNHRETY2NggICCg1eddlwoQJ2Lx5M5555hlMnDgRLVq0gEqlQkpKCnbv3o1Jkyahffv2OpfNzMyEv78/XnzxRfTo0QO+vr7IycnBvn37sHTpUgQGBuL555/XWEYul6Nnz56IioqCSqXCRx99hKysLMycObPUGCvzfZg1axbCwsLQq1cvvPXWW1AqlViwYAEcHBxw584d9To7deqE1157DSNGjMCJEyfwzDPPwN7eHqmpqTh06BCaN2+O//3vfxU+jqXp168f4uPj0bhxY7Ro0QInT57EggULKtTVWvKZXbp0KYYPHw4rKys0atQI/v7+mDVrFqZOnYqLFy+id+/eqFmzJm7evIlff/0V9vb2ZR7P0owePRq2trbo1KkTvLy8kJaWhpiYGDg7O5faBQYUj+VLTk7GihUrcPXqVY3vZO3atcvc1+bNm2Pfvn344Ycf4OXlBUdHRzRq1AizZs1CUlISOnbsiDfffBONGjVCXl4eLl++jMTERCxfvhy1a9euUMxV8T0xe8YczUyPr7RLwStDqVQKX19fAUBMnTpVZ50PP/xQ+Pv7C2traxEYGCi+/PJLnVfzPHq1lBBC/PPPPyI0NFQ4OTmJWrVqiTfeeEPs2LFD52j+/fv3i7CwMOHi4iKsrKyEj4+PCAsLq9CVNMnJyaJnz57CxsZGuLi4iFGjRonvvvuu3Kultm/fLvr06SN8fHyEQqEQ7u7uom/fvurLlUusX79eNG7cWFhZWQkAYvr06er12dvb64yptKtOXn/9dREbGyvq1asnrKysROPGjcXatWu1lv/1119Fx44dhb29vfDx8RHTp08XX331ldYVJZcvXxahoaHC0dFRAFBvU9cVKUIIcfDgQdGtWzdhb28vbG1tRYcOHcQPP/ygUae0z5auy/hLs23bNtG+fXthY2Mj7O3tRffu3cUvv/yic30V+RsLIcSSJUtEQECAkMvlGvsWEhIimjZtqlVf198gJydHvP/++6JRo0ZCoVAIZ2dn0bx5czFx4sQyrzDMz88XCxcuFH369BF16tQR1tbWwsbGRgQGBorJkyeL27dvq+uWHPuPPvpIzJw5U9SuXVsoFArRunVrrUuTS7uqsKLfh61bt4rmzZsLhUIh6tSpIz788EPx5ptvipo1a2rtw6pVq0T79u3Vf/t69eqJiIgIceLEiTKPWcnn9lGPfuf/++8/MWrUKOHu7i7s7OzE008/LQ4ePChCQkI0rugr7bM5ZcoU4e3tLSwsLLQ+Z9u2bRNdu3YVTk5OwtraWvj5+YlBgwZpXBpf1vfxUatXrxZdu3YVHh4eQqFQCG9vbzF48GD1bTBK4+fnp75a7tGp5HehNL/99pvo1KmTsLOzEwA0jsmtW7fEm2++KQICAoSVlZVwcXERQUFBYurUqSInJ6dSMZf2PXlSyIR4ZGg/ERE9tsuXLyMgIAALFixAdHS0QbddWFiIVq1awcfHB7t37zbotolMAbuliIjM3KhRo9CzZ091N8Xy5ctx7tw5LF261NihERkFkxsiIjOXnZ2N6Oho3Lp1C1ZWVmjTpg0SExO1xs8RPSnYLUVERESSwjsUExERkaQwuSEiIiJJYXJDREREkvLEDShWqVS4ceMGHB0dq+RmZURERFT9hBDIzs6Gt7e31vPEHvXEJTc3btyAr6+vscMgIiIiPVy9erXcO14/ccmNo6MjgOKDU5nngRAREZHxZGVlwdfXV30eL8sTl9yUdEU5OTkxuSEiIjIzFRlSwgHFREREJClMboiIiEhSmNwQERGRpDxxY26IiAxFpVKhoKDA2GEQmQ2FQlHuZd4VweSGiKgaFBQU4NKlS1CpVMYOhchsWFhYICAgAAqF4rHWw+SGiKiKCSGQmpoKuVwOX1/fKvlPlEjqSm6ym5qaijp16jzWjXaZ3BARVbGioiLcu3cP3t7esLOzM3Y4RGajVq1auHHjBoqKimBlZaX3evjvBBFRFVMqlQDw2E3rRE+aku9MyXdIX0xuiIiqCZ9fR1Q5VfWdYXJDREREksIxN0REBpKVV4i8gsdrbq8MG4UcTjb6j1swhBkzZmDbtm347bffjB0KSQiTGyIiA8jKK8Saw5dRqBQG26aVXIaIjv4mn+AQVTUmN0REBpBXoEShUqB3M0+42lf/QOPbuQXY+Wca8gqUTG7oicMxN0REBuRqr4C7k021T/okUF26dMH48eMxfvx41KhRA66urnj//fchhHZrU2ZmJmxtbbFz506N8i1btsDe3h45OTkAgHfeeQcNGzaEnZ0d6tatiw8++ACFhYVlxjBhwgSNsoEDByIyMlL9uqCgAJMnT4aPjw/s7e3Rvn177Nu3r9L7S9LFlhuiR1TXuAhzGP9AtHr1aowaNQrHjh3DiRMn8Nprr8HPzw+jR4/WqOfs7IywsDCsXbsWvXv3VpevW7cOAwYMgIODAwDA0dER8fHx8Pb2xpkzZzB69Gg4Ojpi8uTJesc4YsQIXL58GRs2bIC3tze2bt2K3r1748yZM2jQoIHe6yXpYHJD9JDqHBfB8Q+GUZ2Ddp+EBNXX1xcff/wxZDIZGjVqhDNnzuDjjz/WSm4AYOjQoYiIiMC9e/dgZ2eHrKws7NixA5s3b1bXef/999Xz/v7+mDRpEhISEvRObi5cuID169fj2rVr8Pb2BgBER0dj586diIuLw7x58/RaL0kLkxuih1TXuAiOfzCM6h60+yQkqB06dNC410hwcDAWLVqEuXPnIiYmRl2enJyMsLAwWFpa4vvvv8eQIUOwefNmODo6IjQ0VF1v06ZNWLJkCc6fP4+cnBwUFRXByclJ7/hOnToFIQQaNmyoUZ6fnw9XV1e910vSwuTG1OXmAv/fvIucHMDe3rjxPCFKxkWQeanOQbtPeoI6duxYhIeHq197e3vD0tISgwYNwrp16zBkyBCsW7cO4eHhsLQsPrUcPXoUQ4YMwcyZM9GrVy84Oztjw4YNWLRoUanbsbCw0Brj8/AYHZVKBblcjpMnT0Iul2vUK+kKI2JyQ0SSw+RUf0ePHtV63aBBA7i6uupsGRk6dChCQ0Nx9uxZ7N27F7Nnz1a/98svv8DPzw9Tp05Vl125cqXM7deqVQupqanq10qlEn/++Se6du0KAGjdujWUSiXS09PRuXNnvfaRpI9XSxERkdrVq1cRFRWFv//+G+vXr8enn36Kt956q9T6ISEh8PDwwNChQ+Hv748OHTqo36tfvz5SUlKwYcMGXLhwAZ988gm2bt1a5va7deuGHTt2YMeOHfjrr78wbtw43L17V/1+w4YN1WN9tmzZgkuXLuH48eP46KOPkJiY+Nj7T9LAlhsiIgO6nVtg0tuJiIjA/fv30a5dO8jlcrzxxht47bXXSq0vk8nw0ksvYcGCBZg2bZrGewMGDMDEiRMxfvx45OfnIywsDB988AFmzJhR6vpGjhyJ33//HREREbC0tMTEiRPVrTYl4uLiMGfOHEyaNAnXr1+Hq6srgoOD0bdv33L3T6lSQVUNQ7IsZIDcgu0FpkImdN3AQMKysrLg7OyMzMzMxxrUZjAcc2NQ6Vl5WHssBUPb16nSbo3qWi9pqs7jXJl15+Xl4dKlSwgICICNTXFdc7hDcZcuXdCqVSssWbKkegMzEqVKhds5BaiOv4AMgKuDggnOY9L13SlRmfM3W26IiAzAycYKER39+WwpI1IJQKD4b2Epr7onthcpBbLyCqESgLz86mQATG6IiAzEycaKyYYJsJTLYCWvyhYWVRWuy3Cqq4sOMH43HZMbUyeXAyX9yHL+T0CGxRviPVlM5REG1XXSLTJgl6Cpq84uOsD43XRMbkydjQ2wY4exo6AqUl2DSasjUeAN8QxHyv9BV5YhTroWVdcjZbaqq4sOMI1uOiY3RAZgo5DDSi7Dzj/TqmX9VnIZ+rXwhp2i6n5KbucW8IZ4BiD1/6ArqzpPuoD5JXvVreq76ABT6KYzanJz4MABLFiwACdPnkRqaiq2bt2KgQMHlrnM/v37ERUVhbNnz8Lb2xuTJ0/G2LFjDRMwkZ6qczDpvQIltv9xA1tPX6/ydVvJZfCpafvEJyDVSer/Qeurek669KQwanKTm5uLli1bYsSIEXjhhRfKrX/p0iX07dsXo0ePxjfffINffvkF48aNQ61atSq0vFnKzQXc3Yvn09N5Kfj/q66xINV5D5LqHExaXYlTdY+Lqerjbah7yFQHqf4HTWQMRk1u+vTpgz59+lS4/vLly1GnTh31PRgCAwNx4sQJLFy4ULrJDQDcu2fsCEyKIcaC2FRh944hmNtVONXZTWeOfz8iqlpmNebmyJEjGk+bBYBevXph5cqVKCwshJWV+fy4k/6q8+GIAK/iMYTq7Kbj34+IzCq5SUtLg4eHh0aZh4cHioqKkJGRAS8vL61l8vPzkZ+fr36dlZVV7XGSYfDhiObN3FqbiMh8mN1oLZlMc8BdydMjHi0vERMTA2dnZ/Xk6+tb7TESERGR8ZhVcuPp6Ym0NM0++vT0dFhaWsLV1VXnMlOmTEFmZqZ6unr1qiFCJSKSpC5dumDChAlVus7bt2/D3d0dly9frtL1mrpBgwZh8eLFxg5DkswquQkODkZSUpJG2e7du9G2bdtSx9tYW1vDyclJYyIiIt3S09MxZswY1KlTB9bW1vD09ESvXr1w5MgRAMCWLVswe/bsKt1mTEwM+vfvD39//ypdb0XFxsaqH9QYFBSEgwcPVmqZ9u2ewtHDhzTenzFjBmQymcbk6empUWfatGmYO3dutQ+X0Gf/KrtcTEwMZDKZRuKbnZ2NCRMmwM/PD7a2tujYsSOOHz/+uLtTIUYdc5OTk4Pz58+rX1+6dAm//fYbXFxcUKdOHUyZMgXXr1/HmjVrAABjx47FZ599hqioKIwePRpHjhzBypUrsX79emPtQvWzsABCQlCkEriTUwAoq/4qEA7AJKISL7zwAgoLC7F69WrUrVsXN2/exJ49e3Dnzh0AgIuLS5Vu7/79+1i5ciUSExOrdL0VlZCQgAkTJiA2NhadOnXCF198gT59+iA5ORl16tSp0DLLli/Hy4MG4vRvZxAQ4AcAUKoEmjRtip27dquXk8vlKFQ+uDw/sGkz+Pn7Y83XX2PM2P+VGuPj3HiwtP37/cyfcHD1rPRyuo7L8ePHsWLFCrRo0UKjfMxro5F89iy+/vpreHt745tvvkGPHj2QnJwMHx8fvfanomSiZNCKEezbtw9du3bVKh8+fDji4+MRGRmJy5cvazzvZP/+/Zg4caL6Jn7vvPNOpW7iV5lHppsK3gZfU3pWHtYeS8HQ9nU4oJgMpjKfu7y8PFy6dEn9X69abq5W3UKlCndyC+DiZAsre7sy66pZWAC2tuXWLbSxLV63vaJC99C5e/cuatasiX379iEkJERnnS5duqBVq1ZYsmQJLl++jICAAK06ISEh2LdvH4QQWLBgAZYvX47U1FQ0bNgQH3zwAQYNGqSuu2XLFowZMwa3bt0qjvn/j8eFs7/hvSnv4tixY/Dz88PXX3+N06dPY/v27fj+++/L3ZeKat++Pdq0aYNly5apywIDAzFw4EDExMRUaBmlSoXAwCboHdYfU2cUt2otiJmDnTt+wJ5Dx8rc/sIP5+LQ/r3Y9uNPpdZ5nLtMl7Z//Z8dgElTZ5T62ajoccnJyUGbNm0QGxuLOXPmoFWrVliwaDGuZ2Sivk8tfPfddwgLC1PXb9WqFfr164c5c+bojLfU7w4qd/42astNly5dUFZuFR8fr1UWEhKCU6dOVWNUpqc6L33mbfCJDMjBQavICoAHAFWfPsDDrRfu7qXf4yokBHj4IZf+/kBGhna9ospdau/g4AAHBwds27YNHTp0gLW1dZn1fX19kZqaqn6dlpaGHj164JlnngEAvP/++9iyZQuWLVuGBg0a4MCBA3jllVdQq1YtdfJ04MABtG3bVmO9J48fwwv9emP69On48ssv8c4772DGjBn4559/8O2332rFMW/ePMybN6/MWH/88Ud07txZo6ygoAAnT57Eu+++q1EeGhqKw4cP61yPrmXkFhbo3SsUv534FS7//xttayXHpQvn0bpxXVhbW+Opdu0we85c1K1bV2N9z3TsgE8XL4C9pdB5vOfNm4f5H32Isu5d/cP2HXj6kX17ONboyZM1Wox69OyJI0d0719p+wjoPi6vv/46wsLC0KNHD42ERVlUBKVSqZWg2Nra4tAhzS686mBWl4I/6XjpMxFVJ0tLS8THx2P06NFYvnw52rRpg5CQEAwZMkSrywEo7mYpGUeSl5eHgQMHIjg4GDNmzEBubi4WL16Mn3/+GcHBwQCAunXr4tChQ/jiiy/Uyc3ly5fh7e2tsd7p772D555/Xn1yHTJkCF566SUMGDAArVu31opj7NixGDx4cJn7pqsbJCMjA0qlUuctRh69eKW8ZTw9PXHz5i51K0jH4A5otWYNGjZsiJs3b2LOnDkI6fw0zp49q3EBjF8dX+Tn5+P2rXT4+flpbW/c/8aid7/nynz2mKe3N+7ouDt3WuoNKJVK2Dq5aLzvWMMNqalppT5EtKLHZcOGDTh16pTOcTQOjo7o0CEYs2fPRmBgIDw8PLB+/XocO3YMDRo0KGNvqgaTG1OXmwu3uv4YU6REzrl/ASY3ROYrJ0er6OFuKY3OgfT00tfzaPdEOVcZFSkFKvoohmcHPofQ3n1w6OBBHD16FLt37cL8+fPxxYoViBgeCSEAlRAoVKo0xoKMGjUK2dnZSEpKgoWFBZKTk5GXl4eePXtqrL+goEAjQbl//77Gf/fXrl3DiV+PYeGCBeoyhUIBIQRmzpypM2YXF5fHGguk6xYjpd1epKLLPHz3/ebNmyM4OBj16tXD6tWrERUVpX7P9v+7F++V0kpXy80NLi4uej01Ps+uuBXJ2VahblECABsrC8jlFuV2dZW1j1evXsVbb72F3bt3a7XOlIhbvRpjRr8KHx8fyOVytGnTBi+//LJBel+Y3FSxKn/mUW4e3G9nwA6A9s8iEZkVXc+GU6oAWAE2ivLrVma9ACxUKshQ/LtUORZo0zEEbTqGYFzUO4ga/z9MnzED/Qa9jEKVCnmFStzJLVCPBYmZNw87d+7Er7/+CkdHRwCASlWcTO3YsUOr1eTh7hc3Nzf8999/6td/nTsHAAh6qKvq77//Rrt27dC8eXOd0erbLeXm5ga5XK7zFiOPtlo8zjIAYG9vj+bNm+Pff//VKC8ZqF2rVi2dy+m7bwDg5eEOuVyOjFvpGuNqbmfcgqeHR6mJTUX28eTJk0hPT0dQUJD6faVSiQMHDuCzzz5Dyq1M1KtXD/v370dubi6ysrLg5eWF8PBwneO0qhqTmypUHQN/Le/fwxv/P1+dz8uprgcO8kosIuORWxT/d67Pf/0Pa9WiGXYl/lA8+NTCAjZWxd/rrLxCbNq8GbNmzcKPP/6IevXqqZdp0qQJrK2tkZKSUurgZABo3bo1vvnmG/XrzKxMyOUPfuvu3LmD+fPno1mzZqWuQ99uKYVCgaCgICQlJeG5555TlyclJWHAgAE616PPMkDx3fLPnTunlYT8+eefqF27Ntzc3Kp03x4n1oos1717d5w5c0ZjuREjRqBx48aIin5b429ob28Pe3t7/Pfff9j1/y2B1Y3JTRWqloG/D10FUR1JQnU+wBAwvyuxiKRGbmGBiv5bdPv2bbz44osYOXIkWrRoAUdHR5w4cQKLFi7AgAEDYCW3gEwGWMhksJTLcC75LEZGRuKdd95B06ZN1f/pKxQKuLi4IDo6GhMnToRKpcLTTz+NrKwsHD58GA4ODhg+fDiA4ucDTpkyBf/99x9q1qyJli1bQalUYuGCBRgSPhhvvfUW/Pz8cO7cOVy5ckXnuJTH6ZaKiorCsGHD0LZtWwQHB2PFihVISUlRX4X72WefYevWrdizZ0+FlwGA6Oho9O/fH3Xq1EF6ejrmzJmDrKws9X6XOHjwoNYzE6tq3yoaqz776OjoqJVw2tvbw9XVFc2aNcOd3ALs3rULcgsZGjVqhPPnz+Ptt99Go0aNMGLECL33p6KY3FSDKh34K6/6Bws+rDofYFhyJdb1/+4jz77q1l9drUxEFVGRz19hQT6UquJxKXJl2WNdiqrpFg/6cHBwQPv27fHxxx/jwoULKCwshK+vL0aPHo333ntPq/7vp0/h3r17mDNnjsaVMiWXgs+ePRvu7u6IiYnBxYsXUaNGDbRp00ZjXc2bN0fbtm3x7bffYsyYMahfvz4mT52Gzz79BB99GIPw8HCsXbsWvXr1Qo8ePbS6dR5XeHg4bt++jVmzZiE1NRXNmjVDYmKiOonKyMjAhQsXKrUMUDx26KWXXkJGRgZq1aqFDh064OjRoxp18vLysHXrVuzatatK96myseq7j+XJzMrEB1On4tq1a3BxccELL7yAuXPnGuQh10a9z40xVOd9bqrl/iu5uQ8uH83JqVw/vJFV5/152CJEhlaZz7OthRKtaxbCp44frBRlX04NPN59TIxFPRC6gvfQKUtiYiKio6Px559/QilQZes1dZ9//jm+++477N69u/zKZuRxPhuSuM8NSVt1tgpxLA8ZWmU+z4UF+fjv5nU421rB5tGBwjo8zh1oja0yV2KVpmev3nj1739wOeUqvLxrV01gZsDKygqffvqpscOQJCY3ps7CAii5asAMf/ycbKyYhJBkVPTznJcHZN2SwUpuIdnWBwsZ9LwSS7eho4rHcmTlFZZ6/xWpee2114wdgmQxuTF1traAgR40RkRUUVV1JZYu5tySRaaByQ0REemlMldiERkSU2MiomryhF2vQfTYquo7w+TG1N27V/xQPH//0h+iR0QmpeQGZgUFvG0BUWWUfGcevgmgPtgtZeqEAK5ceTBPRCbP0tISdnZ2uHXrFqysrGDB8SP0BClUqlBYUIg8uQrKSgyoV6lUuHXrFuzs7GBp+XjpCZMbIqIqJpPJ4OXlhUuXLuFKyT8nRE8IpUogN78Id60tIa/kZW8WFhaoU6dOuQ8uLQ+TGyKiaqBQKNCgQQN2TdET53ZOPg79fgP9W3rD1aH8m1g+TKFQVElLJ5MbIqJqYmFhoXWXVSKpsyoA7qvksFJYG+3zz45gIiIikhQmN0RERCQp7JYydTIZ0KTJg3kiIiIqE5MbU2dnB5w9a+woiIiIzAa7pYiIiEhSmNwQERGRpDC5MXX37gFNmxZPfPwCERFRuTjmxtQJASQnP5gnIiKiMrHlhoiIiCSFyQ0RERFJCpMbIiIikhQmN0RERCQpTG6IiIhIUni1lKmTyQA/vwfzREREVCYmN6bOzg64fNnYURAREZkNdksRERGRpDC5ISIiIklhcmPq7t8HnnqqeLp/39jREBERmTyOuTF1KhVw4sSDeSIiIioTW26IiIhIUpjcEBERkaQwuSEiIiJJYXJDREREksLkhoiIiCSFV0uZAzc3Y0dARERkNpjcmDp7e+DWLWNHQUREZDbYLUVERESSwuSGiIiIJIXJjam7fx/o0qV44uMXiIiIysUxN6ZOpQL2738wT0RERGViyw0RERFJCpMbIiIikhQmN0RERCQpTG6IiIhIUpjcEBERkaTwailzYGdn7AiIiIjMBpMbU2dvD+TmGjsKIiIis8FuKSIiIpIUJjdEREQkKUZPbmJjYxEQEAAbGxsEBQXh4MGDZdZfu3YtWrZsCTs7O3h5eWHEiBG4ffu2gaI1grw8ICyseMrLM3Y0REREJs+oyU1CQgImTJiAqVOn4vTp0+jcuTP69OmDlJQUnfUPHTqEiIgIjBo1CmfPnsXGjRtx/PhxvPrqqwaO3ICUSiAxsXhSKo0dDRERkckzanKzePFijBo1Cq+++ioCAwOxZMkS+Pr6YtmyZTrrHz16FP7+/njzzTcREBCAp59+GmPGjMGJEycMHDkRERGZKqMlNwUFBTh58iRCQ0M1ykNDQ3H48GGdy3Ts2BHXrl1DYmIihBC4efMmNm3ahLCwMEOETERERGbAaMlNRkYGlEolPDw8NMo9PDyQlpamc5mOHTti7dq1CA8Ph0KhgKenJ2rUqIFPP/201O3k5+cjKytLYyIiIiLpMvqAYplMpvFaCKFVViI5ORlvvvkmpk2bhpMnT2Lnzp24dOkSxo4dW+r6Y2Ji4OzsrJ58fX2rNH4iIiIyLUZLbtzc3CCXy7VaadLT07Vac0rExMSgU6dOePvtt9GiRQv06tULsbGxWLVqFVJTU3UuM2XKFGRmZqqnq1evVvm+EBERkekwWnKjUCgQFBSEpKQkjfKkpCR07NhR5zL37t2DhYVmyHK5HEBxi48u1tbWcHJy0piIiIhIuoz6+IWoqCgMGzYMbdu2RXBwMFasWIGUlBR1N9OUKVNw/fp1rFmzBgDQv39/jB49GsuWLUOvXr2QmpqKCRMmoF27dvD29jbmrlQfe3uglMSNiIiItBk1uQkPD8ft27cxa9YspKamolmzZkhMTISfnx8AIDU1VeOeN5GRkcjOzsZnn32GSZMmoUaNGujWrRs++ugjY+0CERERmRiZKK0/R6KysrLg7OyMzMzMKu+iSs/Kw9pjKRjavg7cnWyqdN1ERETmoLrOhZU5fxv9aikqR14e8OKLxRMfv0BERFQuJjemTqkENm0qnvj4BSIionIxuSEiIiJJYXJDREREksLkhoiIiCSFyQ0RERFJCpMbIiIikhQmN0RERCQpRr1DMVWAnR2Qk/NgnoiIiMrE5MbUyWTFz5ciIiKiCmG3FBEREUkKkxtTl58PREYWT/n5xo6GiIjI5DG5MXVFRcDq1cVTUZGxoyEiIjJ5TG6IiIhIUpjcEBERkaQwuSEiIiJJYXJDREREksLkhoiIiCSFyQ0RERFJCu9QbOrs7ID09AfzREREVCYmN6ZOJgNq1TJ2FERERGaD3VJEREQkKUxuTF1+PvD668UTH79ARERULiY3pq6oCIiNLZ74+AUiIqJyMbkhIiIiSWFyQ0RERJLC5IaIiIgkhckNERERSQqTGyIiIpIUJjdEREQkKbxDsamztQUuXXowT0RERGVicmPqLCwAf39jR0FERGQ22C1FREREksLkxtQVFABvv108FRQYOxoiIiKTx+TG1BUWAgsXFk+FhcaOhoiIyOQxuSEiIiJJYXJDREREksLkhoiIiCSFyQ0RERFJCpMbIiIikhQmN0RERCQpvEOxqbO1Bf7888E8ERERlYnJjamzsACaNjV2FERERGaD3VJEREQkKWy5MXUFBcC8ecXz770HKBTGjYeIiMjEMbkxdYWFwMyZxfNvv83khoiIqBzsliIiIiJJYXJDREREksLkhoiIiCSFyQ0RERFJCpMbIiIikhQmN0RERCQpvBTc1NnYAL/++mCeiIiIysTkxtTJ5cBTTxk7CiIiIrOhV7fUpUuXqjoOIiIioiqhV3JTv359dO3aFd988w3y8vKqOiZ6WEEBsGBB8VRQYOxoiIiITJ5eyc3vv/+O1q1bY9KkSfD09MSYMWPwa8m4kEqKjY1FQEAAbGxsEBQUhIMHD5ZZPz8/H1OnToWfnx+sra1Rr149rFq1Sq9tm4XCQmDy5OKpsNDY0RAREZk8vZKbZs2aYfHixbh+/Tri4uKQlpaGp59+Gk2bNsXixYtx69atCq0nISEBEyZMwNSpU3H69Gl07twZffr0QUpKSqnLDB48GHv27MHKlSvx999/Y/369WjcuLE+u0FEREQSJBNCiMddSX5+PmJjYzFlyhQUFBTAysoK4eHh+Oijj+Dl5VXqcu3bt0ebNm2wbNkydVlgYCAGDhyImJgYrfo7d+7EkCFDcPHiRbi4uOgVa1ZWFpydnZGZmQknJye91lGa9Kw8rD2WgqHt68DdqYqubMrNBRwciudzcgB7+6pZLxERUTWolnMhKnf+fqz73Jw4cQLjxo2Dl5cXFi9ejOjoaFy4cAE///wzrl+/jgEDBpS6bEFBAU6ePInQ0FCN8tDQUBw+fFjnMt9//z3atm2L+fPnw8fHBw0bNkR0dDTu37//OLtBREREEqLXpeCLFy9GXFwc/v77b/Tt2xdr1qxB3759YWFRnCsFBATgiy++KLO7KCMjA0qlEh4eHhrlHh4eSEtL07nMxYsXcejQIdjY2GDr1q3IyMjAuHHjcOfOnVLH3eTn5yM/P1/9Oisrq7K7S0RERGZEr+Rm2bJlGDlyJEaMGAFPT0+dderUqYOVK1eWuy6ZTKbxWgihVVZCpVJBJpNh7dq1cHZ2BlCcaA0aNAiff/45bG1ttZaJiYnBzJkzy42DiIiIpEGvbqmkpCS88847WomNEEI9GFihUGD48OGlrsPNzQ1yuVyrlSY9PV2rNaeEl5cXfHx81IkNUDxGRwiBa9eu6VxmypQpyMzMVE9Xr16t0D4SERGRedIrualXrx4yMjK0yu/cuYOAgIAKrUOhUCAoKAhJSUka5UlJSejYsaPOZTp16oQbN24gJydHXfbPP//AwsICtWvX1rmMtbU1nJycNCazYmMD7N1bPPHxC0REROXSK7kp7QKrnJwc2FTiBBwVFYWvvvoKq1atwrlz5zBx4kSkpKRg7NixAIpbXSIiItT1X375Zbi6umLEiBFITk7GgQMH8Pbbb2PkyJE6u6QkQS4HunQpnuRyY0dDRERk8io15iYqKgpA8TiZadOmwc7OTv2eUqnEsWPH0KpVqwqvLzw8HLdv38asWbOQmpqKZs2aITExEX5+fgCA1NRUjXveODg4ICkpCW+88Qbatm0LV1dXDB48GHPmzKnMbhAREZGEVSq5OX36NIDilpszZ85AoVCo31MoFGjZsiWio6MrFcC4ceMwbtw4ne/Fx8drlTVu3FirK0vSCguBFSuK5197DbCyMm48REREJq5Syc3evXsBACNGjMDSpUvNb/yKOSooAMaPL56PjGRyQ0REVA69LgWPi4ur6jiIiIiIqkSFk5vnn38e8fHxcHJywvPPP19m3S1btjx2YERERET6qHBy4+zsrL653sP3mSEiIiIyJRVObh7uimK3FBEREZkqve5zc//+fdy7d0/9+sqVK1iyZAl2795dZYERERER6UOv5GbAgAFYs2YNAODu3bto164dFi1ahAEDBmDZsmVVGiARERFRZeiV3Jw6dQqdO3cGAGzatAmenp64cuUK1qxZg08++aRKA3ziWVsD27cXT9bWxo6GiIjI5Ol1Kfi9e/fg6OgIANi9ezeef/55WFhYoEOHDrhy5UqVBvjEs7QEwsKMHQUREZHZ0Kvlpn79+ti2bRuuXr2KXbt2ITQ0FEDxE715Yz8iIiIyJr2Sm2nTpiE6Ohr+/v5o3749goODARS34rRu3bpKA3ziFRYC8fHFU2GhsaMhIiIyeXp1Sw0aNAhPP/00UlNT0bJlS3V59+7d8dxzz1VZcITixy+MGFE8/+KLfPwCERFROfRKbgDA09MTnp6eGmXt2rV77ICIiIiIHodeyU1ubi4+/PBD7NmzB+np6VCpVBrvX7x4sUqCIyIiIqosvZKbV199Ffv378ewYcPg5eWlfiwDERERkbHpldz8+OOP2LFjBzp16lTV8RARERE9Fr2ulqpZsyZcXFyqOhYiIiKix6ZXcjN79mxMmzZN4/lSRERERKZAr26pRYsW4cKFC/Dw8IC/vz+sHrk8+dSpU1USHKH4kQvffvtgnoiIiMqkV3IzcODAKg6DSmVpWXx/GyIiIqoQvZKb6dOnV3UcRERERFVCrzE3AHD37l189dVXmDJlCu7cuQOguDvq+vXrVRYcASgqAjZuLJ6KiowdDRERkcnTq+Xmjz/+QI8ePeDs7IzLly9j9OjRcHFxwdatW3HlyhWsWbOmquN8cuXnA4MHF8/n5BR3UxEREVGp9Gq5iYqKQmRkJP7991/Y2Nioy/v06YMDBw5UWXBERERElaVXcnP8+HGMGTNGq9zHxwdpaWmPHRQRERGRvvRKbmxsbJCVlaVV/vfff6NWrVqPHRQRERGRvvRKbgYMGIBZs2ahsLAQACCTyZCSkoJ3330XL7zwQpUGSERERFQZeiU3CxcuxK1bt+Du7o779+8jJCQE9evXh6OjI+bOnVvVMRIRERFVmF6X3jg5OeHQoUPYu3cvTp48CZVKhTZt2qBHjx5VHR8RERFRpVQ6uVGpVIiPj8eWLVtw+fJlyGQyBAQEwNPTE0IIyGSy6ojzyaVQAHFxD+aJiIioTJVKboQQePbZZ5GYmIiWLVuiefPmEELg3LlziIyMxJYtW7Bt27ZqCvUJZWUFREYaOwoiIiKzUankJj4+HgcOHMCePXvQtWtXjfd+/vlnDBw4EGvWrEFERESVBklERERUUZUaULx+/Xq89957WokNAHTr1g3vvvsu1q5dW2XBEYofubBjR/HExy8QERGVq1LJzR9//IHevXuX+n6fPn3w+++/P3ZQ9JD8fKBfv+IpP9/Y0RAREZm8SiU3d+7cgYeHR6nve3h44L///nvsoIiIiIj0VankRqlUwrKMBzfK5XIUseuEiIiIjKjSV0tFRkbC2tpa5/v57DYhIiIiI6tUcjN8+PBy6/BKKSIiIjKmSiU3cSU3kyMiIiIyUXo9W4qIiIjIVOn1bCkyIIUC+OyzB/NERERUJiY3ps7KCnj9dWNHQUREZDbYLUVERESSwpYbU6dUAgcPFs937gzI5caNh4iIyMQxuTF1eXlAybO8cnIAe3vjxkNERGTi2C1FREREksLkhoiIiCSFyQ0RERFJCpMbIiIikhQmN0RERCQpTG6IiIhIUngpuKmzsgLmz38wT0RERGVicmPqFArg7beNHQUREZHZYLcUERERSQpbbkydUgmcOlU836YNH79ARERUDiY3pi4vD2jXrniej18gIiIql9G7pWJjYxEQEAAbGxsEBQXhYMlDIsvxyy+/wNLSEq1atareAImIiMisGDW5SUhIwIQJEzB16lScPn0anTt3Rp8+fZCSklLmcpmZmYiIiED37t0NFCkRERGZC6MmN4sXL8aoUaPw6quvIjAwEEuWLIGvry+WLVtW5nJjxozByy+/jODgYANFSkRERObCaMlNQUEBTp48idDQUI3y0NBQHD58uNTl4uLicOHCBUyfPr26QyQiIiIzZLQBxRkZGVAqlfDw8NAo9/DwQFpams5l/v33X7z77rs4ePAgLC0rFnp+fj7y8/PVr7OysvQPmoiIiEye0QcUy2QyjddCCK0yAFAqlXj55Zcxc+ZMNGzYsMLrj4mJgbOzs3ry9fV97JiJiIjIdBmt5cbNzQ1yuVyrlSY9PV2rNQcAsrOzceLECZw+fRrjx48HAKhUKgghYGlpid27d6Nbt25ay02ZMgVRUVHq11lZWeaV4FhZASVdcHz8AhERUbmMltwoFAoEBQUhKSkJzz33nLo8KSkJAwYM0Krv5OSEM2fOaJTFxsbi559/xqZNmxAQEKBzO9bW1rC2tq7a4A1JoQBmzDB2FERERGbDqDfxi4qKwrBhw9C2bVsEBwdjxYoVSElJwdixYwEUt7pcv34da9asgYWFBZo1a6axvLu7O2xsbLTKiYiI6Mll1OQmPDwct2/fxqxZs5CamopmzZohMTERfn5+AIDU1NRy73kjeSoVcO5c8XxgIGBh9GFSREREJk0mhBDGDsKQsrKy4OzsjMzMTDg5OVXputOz8rD2WAqGtq8Ddyebqllpbi7g4FA8z8cvEBGRiauWcyEqd/5mMwARERFJCpMbIiIikhQmN0RERCQpTG6IiIhIUpjcEBERkaQwuSEiIiJJMep9bqgCrKyA6OgH80RERFQmJjemTqEAFiwwdhRERERmg91SREREJClsuTF1KhVQ8giKOnX4+AUiIqJyMLkxdffvAyVPPOfjF4iIiMrFZgAiIiKSFCY3REREJClMboiIiEhSmNwQERGRpDC5ISIiIklhckNERESSwkvBTZ2lJTBu3IN5IiIiKhPPlqbO2hr4/HNjR0FERGQ22C1FREREksKWG1MnBJCRUTzv5gbIZMaNh4iIyMQxuTF19+4B7u7F83z8AhERUbnYLUVERESSwuSGiIiIJIXJDREREUkKkxsiIiKSFCY3REREJClMboiIiEhSeCm4qbO0BIYPfzBPREREZeLZ0tRZWwPx8caOgoiIyGywW4qIiIgkhS03pk6I4rsUA4CdHR+/QEREVA623Ji6e/cAB4fiqSTJISIiolIxuSEiIiJJYXJDREREksLkhoiIiCSFyQ0RERFJCpMbIiIikhQmN0RERCQpvM+NqZPLgUGDHswTERFRmZjcmDobG2DjRmNHQUREZDbYLUVERESSwuSGiIiIJIXJjanLzS1+npRMVjxPREREZWJyQ0RERJLC5IaIiIgkhckNERERSQqTGyIiIpIUJjdEREQkKUxuiIiISFJ4h2JTJ5cDffs+mCciIqIyMbkxdTY2wI4dxo6CiIjIbLBbioiIiCSFyQ0RERFJCpMbU5ebC9jbF098/AIREVG5OObGHNy7Z+wIiIiIzIbRW25iY2MREBAAGxsbBAUF4eDBg6XW3bJlC3r27IlatWrByckJwcHB2LVrlwGjJSIiIlNn1OQmISEBEyZMwNSpU3H69Gl07twZffr0QUpKis76Bw4cQM+ePZGYmIiTJ0+ia9eu6N+/P06fPm3gyImIiMhUyYQQwlgbb9++Pdq0aYNly5apywIDAzFw4EDExMRUaB1NmzZFeHg4pk2bVqH6WVlZcHZ2RmZmJpycnPSKuzTpWXlYeywFQ9vXgbuTTdWsNDcXcHAons/JKR57Q0REZKKq5VyIyp2/jdZyU1BQgJMnTyI0NFSjPDQ0FIcPH67QOlQqFbKzs+Hi4lIdIRIREZEZMtqA4oyMDCiVSnh4eGiUe3h4IC0trULrWLRoEXJzczF48OBS6+Tn5yM/P1/9OisrS7+AiYiIyCwYfUCxTCbTeC2E0CrTZf369ZgxYwYSEhLg7u5ear2YmBg4OzurJ19f38eO2aAsLICQkOLJwuh/LiIiIpNntLOlm5sb5HK5VitNenq6VmvOoxISEjBq1Ch8++236NGjR5l1p0yZgszMTPV09erVx47doGxtgX37iidbW2NHQ0REZPKMltwoFAoEBQUhKSlJozwpKQkdO3Ysdbn169cjMjIS69atQ1hYWLnbsba2hpOTk8ZERERE0mXUm/hFRUVh2LBhaNu2LYKDg7FixQqkpKRg7NixAIpbXa5fv441a9YAKE5sIiIisHTpUnTo0EHd6mNrawtnZ2ej7QcRERGZDqMmN+Hh4bh9+zZmzZqF1NRUNGvWDImJifDz8wMApKamatzz5osvvkBRURFef/11vP766+ry4cOHIz4+3tDhG0ZuLuDvXzx/+TIvBSciIiqH0R+/MG7cOIwbN07ne48mLPv27av+gExRRoaxIyAiIjIbvPyGiIiIJIXJDREREUkKkxsiIiKSFCY3REREJClMboiIiEhSjH61FJXDwgJo2/bBPBEREZWJyY2ps7UFjh83dhRERERmg00BREREJClMboiIiEhSmNyYunv3ih+/4O9fPE9ERERl4pgbUycEcOXKg3kiIiIqE1tuiIiISFKY3BAREZGkMLkhIiIiSWFyQ0RERJLC5IaIiIgkhVdLmTqZDGjS5ME8ERERlYnJjamzswPOnjV2FERERGaD3VJEREQkKUxuiIiISFKY3Ji6e/eApk2LJz5+gYiIqFwcc2PqhACSkx/MExERUZnYckNERESSwuSGiIiIJIXJDREREUkKkxsiIiKSFCY3REREJCm8WsrUyWSAn9+DeSIiIioTkxtTZ2cHXL5s7CiIiIjMBruliIiISFKY3BAREZGkMLkxdffvA089VTzdv2/saIiIiEwex9yYOpUKOHHiwTwRERGViS03REREJClMboiIiEhSmNwQERGRpDC5ISIiIklhckNERESSwqulzIGbm7EjICIiMhtMbkydvT1w65axoyAiIjIb7JYiIiIiSWFyQ0RERJLC5MbU3b8PdOlSPPHxC0REROXimBtTp1IB+/c/mCciIqIyseWGiIiIJIXJDREREUkKkxsiIiKSFCY3REREJClMboiIiEhSeLWUObCzM3YEREREZoPJjamztwdyc40dBRERkdlgtxQRERFJCpMbIiIikhQmN6YuLw8ICyue8vKMHQ0REZHJ45gbU6dUAomJD+aJiIioTGy5ISIiIkkxenITGxuLgIAA2NjYICgoCAcPHiyz/v79+xEUFAQbGxvUrVsXy5cvN1CkREREZA6MmtwkJCRgwoQJmDp1Kk6fPo3OnTujT58+SElJ0Vn/0qVL6Nu3Lzp37ozTp0/jvffew5tvvonNmzcbOHIiIiIyVUZNbhYvXoxRo0bh1VdfRWBgIJYsWQJfX18sW7ZMZ/3ly5ejTp06WLJkCQIDA/Hqq69i5MiRWLhwoYEjJyIiIlNltOSmoKAAJ0+eRGhoqEZ5aGgoDh8+rHOZI0eOaNXv1asXTpw4gcLCwmqLlYiIiMyH0a6WysjIgFKphIeHh0a5h4cH0tLSdC6Tlpams35RUREyMjLg5eWltUx+fj7y8/PVrzMzMwEAWVlZj7sLWrKz8pCXm4PsrCzYoKBqVvrw3YmzsnjFFBERmbRqORfiwXlbCFFuXaNfCi6TyTReCyG0ysqrr6u8RExMDGbOnKlV7uvrW9lQK+y96lqxt3d1rZmIiKhKVde5MDs7G87OzmXWMVpy4+bmBrlcrtVKk56ertU6U8LT01NnfUtLS7i6uupcZsqUKYiKilK/VqlUuHPnDlxdXctMovSRlZUFX19fXL16FU5OTlW6bnqAx9kweJwNh8faMHicDaO6jrMQAtnZ2fCuwD/6RktuFAoFgoKCkJSUhOeee05dnpSUhAEDBuhcJjg4GD/88ING2e7du9G2bVtYWVnpXMba2hrW1tYaZTVq1Hi84Mvh5OTEL44B8DgbBo+z4fBYGwaPs2FUx3Eur8WmhFGvloqKisJXX32FVatW4dy5c5g4cSJSUlIwduxYAMWtLhEREer6Y8eOxZUrVxAVFYVz585h1apVWLlyJaKjo421C0RERGRijDrmJjw8HLdv38asWbOQmpqKZs2aITExEX5+fgCA1NRUjXveBAQEIDExERMnTsTnn38Ob29vfPLJJ3jhhReMtQtERERkYow+oHjcuHEYN26czvfi4+O1ykJCQnDq1Klqjko/1tbWmD59ulY3GFUtHmfD4HE2HB5rw+BxNgxTOM4yUZFrqoiIiIjMhNGfLUVERERUlZjcEBERkaQwuSEiIiJJYXJDREREksLkppJiY2MREBAAGxsbBAUF4eDBg2XW379/P4KCgmBjY4O6deti+fLlBorUvFXmOG/ZsgU9e/ZErVq14OTkhODgYOzatcuA0Zqvyn6eS/zyyy+wtLREq1atqjdAiajscc7Pz8fUqVPh5+cHa2tr1KtXD6tWrTJQtOatssd67dq1aNmyJezs7ODl5YURI0bg9u3bBorW/Bw4cAD9+/eHt7c3ZDIZtm3bVu4yRjkPCqqwDRs2CCsrK/Hll1+K5ORk8dZbbwl7e3tx5coVnfUvXrwo7OzsxFtvvSWSk5PFl19+KaysrMSmTZsMHLl5qexxfuutt8RHH30kfv31V/HPP/+IKVOmCCsrK3Hq1CkDR25eKnucS9y9e1fUrVtXhIaGipYtWxomWDOmz3F+9tlnRfv27UVSUpK4dOmSOHbsmPjll18MGLV5quyxPnjwoLCwsBBLly4VFy9eFAcPHhRNmzYVAwcONHDk5iMxMVFMnTpVbN68WQAQW7duLbO+sc6DTG4qoV27dmLs2LEaZY0bNxbvvvuuzvqTJ08WjRs31igbM2aM6NChQ7XFKAWVPc66NGnSRMycObOqQ5MUfY9zeHi4eP/998X06dOZ3FRAZY/zjz/+KJydncXt27cNEZ6kVPZYL1iwQNStW1ej7JNPPhG1a9euthilpCLJjbHOg+yWqqCCggKcPHkSoaGhGuWhoaE4fPiwzmWOHDmiVb9Xr144ceIECgsLqy1Wc6bPcX6USqVCdnY2XFxcqiNESdD3OMfFxeHChQuYPn16dYcoCfoc5++//x5t27bF/Pnz4ePjg4YNGyI6Ohr37983RMhmS59j3bFjR1y7dg2JiYkQQuDmzZvYtGkTwsLCDBHyE8FY50Gj36HYXGRkZECpVGo9sdzDw0PrSeUl0tLSdNYvKipCRkYGvLy8qi1ec6XPcX7UokWLkJubi8GDB1dHiJKgz3H+999/8e677+LgwYOwtORPR0Xoc5wvXryIQ4cOwcbGBlu3bkVGRgbGjRuHO3fucNxNGfQ51h07dsTatWsRHh6OvLw8FBUV4dlnn8Wnn35qiJCfCMY6D7LlppJkMpnGayGEVll59XWVk6bKHucS69evx4wZM5CQkAB3d/fqCk8yKnqclUolXn75ZcycORMNGzY0VHiSUZnPs0qlgkwmw9q1a9GuXTv07dsXixcvRnx8PFtvKqAyxzo5ORlvvvkmpk2bhpMnT2Lnzp24dOmS+uHNVDWMcR7kv18V5ObmBrlcrvUfQHp6ulZWWsLT01NnfUtLS7i6ulZbrOZMn+NcIiEhAaNGjcLGjRvRo0eP6gzT7FX2OGdnZ+PEiRM4ffo0xo8fD6D4JCyEgKWlJXbv3o1u3boZJHZzos/n2cvLCz4+PnB2dlaXBQYGQgiBa9euoUGDBtUas7nS51jHxMSgU6dOePvttwEALVq0gL29PTp37ow5c+awdb0KGOs8yJabClIoFAgKCkJSUpJGeVJSEjp27KhzmeDgYK36u3fvRtu2bWFlZVVtsZozfY4zUNxiExkZiXXr1rG/vAIqe5ydnJxw5swZ/Pbbb+pp7NixaNSoEX777Te0b9/eUKGbFX0+z506dcKNGzeQk5OjLvvnn39gYWGB2rVrV2u85kyfY33v3j1YWGieBuVyOYAHrQv0eIx2HqzW4coSU3KZ4cqVK0VycrKYMGGCsLe3F5cvXxZCCPHuu++KYcOGqeuXXAI3ceJEkZycLFauXMlLwSugssd53bp1wtLSUnz++eciNTVVPd29e9dYu2AWKnucH8WrpSqmssc5Oztb1K5dWwwaNEicPXtW7N+/XzRo0EC8+uqrxtoFs1HZYx0XFycsLS1FbGysuHDhgjh06JBo27ataNeunbF2weRlZ2eL06dPi9OnTwsAYvHixeL06dPqy+1N5TzI5KaSPv/8c+Hn5ycUCoVo06aN2L9/v/q94cOHi5CQEI36+/btE61btxYKhUL4+/uLZcuWGThi81SZ4xwSEiIAaE3Dhw83fOBmprKf54cxuam4yh7nc+fOiR49eghbW1tRu3ZtERUVJe7du2fgqM1TZY/1J598Ipo0aSJsbW2Fl5eXGDp0qLh27ZqBozYfe/fuLfP31lTOgzIh2PZGRERE0sExN0RERCQpTG6IiIhIUpjcEBERkaQwuSEiIiJJYXJDREREksLkhoiIiCSFyQ0RERFJCpMbIpK8GTNmoFWrVsYOg4gMhMkNERERSQqTGyIiIpIUJjdEZFBdunTB+PHjMX78eNSoUQOurq54//33dT6FOTMzE7a2tti5c6dG+ZYtW2Bvb69+cvY777yDhg0bws7ODnXr1sUHH3yAwsLCMmOYMGGCRtnAgQMRGRmpfl1QUIDJkyfDx8cH9vb2aN++Pfbt26f3fhOR4TC5ISKDW716NSwtLXHs2DF88skn+Pjjj/HVV19p1XN2dkZYWBjWrl2rUb5u3ToMGDAADg4OAABHR0fEx8cjOTkZS5cuxZdffomPP/74sWIcMWIEfvnlF2zYsAF//PEHXnzxRfTu3Rv//vvvY62XiKqfpbEDIKInj6+vLz7++GPIZDI0atQIZ86cwccff4zRo0dr1R06dCgiIiJw79492NnZISsrCzt27MDmzZvVdd5//331vL+/PyZNmoSEhARMnjxZr/guXLiA9evX49q1a/D29gYAREdHY+fOnYiLi8O8efP0Wi8RGQZbbojI4Dp06ACZTKZ+HRwcjH///Rdz586Fg4ODekpJSUFYWBgsLS3x/fffAwA2b94MR0dHhIaGqpfftGkTnn76aXh6esLBwQEffPABUlJS9I7v1KlTEEKgYcOGGvHs378fFy5c0H/Hicgg2HJDRCZj7NixCA8PV7/29vaGpaUlBg0ahHXr1mHIkCFYt24dwsPDYWlZ/PN19OhRDBkyBDNnzkSvXr3g7OyMDRs2YNGiRaVux8LCQmuMz8NjdFQqFeRyOU6ePAm5XK5Rr6QrjIhMF5MbIjK4o0ePar1u0KABXF1d4erqqlV/6NChCA0NxdmzZ7F3717Mnj1b/d4vv/wCPz8/TJ06VV125cqVMrdfq1YtpKamql8rlUr8+eef6Nq1KwCgdevWUCqVSE9PR+fOnfXaRyIyHnZLEZHBXb16FVFRUfj777+xfv16fPrpp3jrrbdKrR8SEgIPDw8MHToU/v7+6NChg/q9+vXrIyUlBRs2bMCFCxfwySefYOvWrWVuv1u3btixYwd27NiBv/76C+PGjcPdu3fV7zds2FA91mfLli24dOkSjh8/jo8++giJiYmPvf9EVL2Y3BCRwUVEROD+/fto164dXn/9dbzxxht47bXXSq0vk8nw0ksv4ffff8fQoUM13hswYAAmTpyI8ePHo1WrVjh8+DA++OCDMrc/cuRIDB8+HBEREQgJCUFAQIC61aZEXFwcIiIiMGnSJDRq1AjPPvssjh07Bl9fX/13nIgMQiZ03VyCiKiadOnSBa1atcKSJUuMHQoRSRRbboiIiEhSmNwQERGRpLBbioiIiCSFLTdEREQkKUxuiIiISFKY3BAREZGkMLkhIiIiSWFyQ0RERJLC5IaIiIgkhckNERERSQqTGyIiIpIUJjdEREQkKf8HXGXrz5Zf56AAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "value,_,_=plt.hist(p_values, bins=20, density=True, alpha=0.5, label='p-value',histtype='step')\n",
    "plt.vlines(size,ymin=0,ymax=np.max(value),linestyle='--',label=r'Size($\\alpha=0.05$) = '+str(round(size,3)),colors='r')\n",
    "plt.legend(loc='upper right')\n",
    "plt.xlabel('p-value')\n",
    "plt.ylabel('Density')\n",
    "plt.title(\"P-value distribution of the Spiegelhalter's z test\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that Spiegelhalter's Z test has a accurate size."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reference\n",
    "\n",
    "Spiegelhalter, D. J. (1986). Probabilistic prediction in patient management and clinical trials.\n",
    "\n",
    "Bröcker, J. (2009). Reliability, Sufficiency, and the Decomposition of Proper Scores. Quarterly Journal of the Royal Meteorological Society, 135(643), 1512–1519. https://doi.org/10.1002/qj.456\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "uq",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}