arviz-psis-discrepancy.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0f988172-961e-4345-9ea3-c2903527254c",
   "metadata": {},
   "source": [
    "# Discrepancy of Generalized Pareto parameter estimate $\\hat k$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "9b17b87e-165a-4847-a4ea-bbc3ea708efa",
   "metadata": {},
   "outputs": [],
   "source": [
    "# uncomment to install arviz in the active environment\n",
    "#!pip install arviz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "2a0d52fc-c273-47f4-afcf-5abdba31a4bc",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import arviz as az\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import py.psis as psis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "881a5b6f-cfca-411b-a126-ca375fa3d7a3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.75904637, -1.42879456, -1.4685748 , -1.6957333 , -1.59647378,\n",
       "       -1.53114609, -1.55555131, -1.61090202])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data = az.load_arviz_data(\"non_centered_eight\")\n",
    "log_weights = data.log_likelihood[\"obs\"].stack(\n",
    "    __sample__=[\"chain\", \"draw\"]\n",
    ")\n",
    "\n",
    "lw_new, khat = az.psislw(log_weights)\n",
    "khats_az = khat.values\n",
    "\n",
    "khats_az"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "792d4c6b-a462-4be9-aa35-69ca3fc05608",
   "metadata": {},
   "source": [
    "We reload the dataset to make sure the computation gets the same inputs. This is necessary since while I write this `az.psislw` is mutating its input. See the [issue about it here](https://github.com/arviz-devs/arviz/issues/2376)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "b9a2e41b-58a1-4e6a-a01e-02f2263464b7",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/ho/code/PSIS/py/psis.py:297: RuntimeWarning: overflow encountered in exp\n",
      "  np.exp(temp, out=temp)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([-0.51473359, -2.94129762, -2.62328168, -3.22422551, -2.29807491,\n",
       "       -3.10027427, -0.95715744, -2.76054336])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data = az.load_arviz_data(\"non_centered_eight\")\n",
    "log_weights = data.log_likelihood[\"obs\"].stack(\n",
    "    __sample__=[\"chain\", \"draw\"]\n",
    ")\n",
    "khats = np.zeros(len(log_weights.values))\n",
    "for i,lw in enumerate(log_weights.values):\n",
    "    khat, sigmahat = psis.gpdfitnew(np.exp(lw))\n",
    "    khats[i] = khat\n",
    "\n",
    "khats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "e81eeb8b-46d2-4382-be5b-8a2197f234fd",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkwAAAG6CAYAAAD6RxlaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA0QUlEQVR4nO3de1zUdb7H8feAXARhCEVnKDTUPSFLZVoqHk/R0ZK9mB7Nbl42M9vYzfLWap3MJdvjVmZ2OavZWubJrlt2Ms968tbdoJXISKUwlUKQEgGVQGR+5w8fzHHi8gP8MTMMr+fjMY9Hv+/v+/t+PzPjPnjv7/Idm2EYhgAAANCkIF8XAAAA4O8ITAAAACYITAAAACYITAAAACYITAAAACYITAAAACYITAAAACYITAAAACa6+LqAQOByuXTo0CFFRUXJZrP5uhwAANAChmHo2LFjio+PV1BQ8+eQCEwWOHTokBISEnxdBgAAaINvv/1W5513XrN9CEwWiIqKknT6A4+OjvZxNQAAoCUqKyuVkJDg/jveHAKTBeovw0VHRxOYAADoYFpyOw03fQMAAJggMAEAAJggMAEAAJggMAEAAJggMAEAAJggMAEAAJggMAEAAJgIqMBUVlamSZMmKTo6WjExMZo+fbqOHz/e7DFpaWmy2Wwer9tvv91LFQMAgI4goBaunDRpkoqLi7V582bV1tZq2rRpuu222/Tiiy82e9yMGTP0wAMPuLcjIiLau1QAANCBBExg2rNnjzZt2qRPP/1Ul156qSTpySef1C9/+UstXbpU8fHxTR4bEREhh8PhrVIBAEAL1bkMZe8vU+mxavWMCteQxFgFB3n/h+4D5pLcjh07FBMT4w5LkjRq1CgFBQUpKyur2WPXrVunHj16KCUlRffcc4+qqqrau1wAAGBiU16xRjy0TTc+84nuejlXNz7ziUY8tE2b8oq9XkvAnGEqKSlRz549Pdq6dOmi2NhYlZSUNHncTTfdpD59+ig+Pl67du3S/PnzlZ+frzfeeKPJY2pqalRTU+PerqysPPs3AAAA3DblFSvjhRwZP2kvqahWxgs5WjF5kNJTnF6rx+/PMC1YsKDBTdk/fe3du7fN4992220aPXq0LrzwQk2aNElr167V+vXrtW/fviaPWbJkiex2u/uVkJDQ5vkBAICnOpehzA27G4QlSe62zA27VedqrEf78PszTHPnztXNN9/cbJ++ffvK4XCotLTUo/3UqVMqKytr1f1JQ4cOlSQVFBSoX79+jfa55557NGfOHPd2ZWUloQkAAItk7y9TcUV1k/sNScUV1creX6bUft29UpPfB6a4uDjFxcWZ9ktNTVV5ebl27typwYMHS5K2bdsml8vlDkEtkZubK0lyOps+zRcWFqawsLAWjwkAAFqu9FjTYakt/azg95fkWmrAgAFKT0/XjBkzlJ2drY8++kh33HGHbrjhBvcTckVFRUpKSlJ2drYkad++fVq8eLF27typAwcO6K233tLUqVN1+eWX66KLLvLl2wEAoNPqGRVuaT8rBExgkk4/7ZaUlKSRI0fql7/8pUaMGKFVq1a599fW1io/P9/9FFxoaKi2bNmiq6++WklJSZo7d64mTJigDRs2+OotAADQ6Q1JjJXTHq6mFg+wSXLaTy8x4C02wzC8d8dUgKqsrJTdbldFRYWio6N9XQ4AAB1e/VNykjxu/q4PUVY8Jdeav98BdYYJAAAEhvQUp1ZMHiSH3fOym8Me7vUlBaQOcNM3AADonNJTnLoq2eEXK30TmAAAgN8KDrJ5bemA5nBJDgAAwASBCQAAwASBCQAAwASBCQAAwASBCQAAwASBCQAAwASBCQAAwASBCQAAwASBCQAAwASBCQAAwASBCQAAwASBCQAAwASBCQAAwASBCQAAwASBCQAAwASBCQAAwASBCQAAwASBCQAAwASBCQAAwASBCQAAwASBCQAAwASBCQAAwASBCQAAwASBCQAAwASBCQAAwASBCQAAwASBCQAAwASBCQAAwASBCQAAwASBCQAAwASBCQAAwASBCQAAwASBCQAAwASBCQAAwASBCQAAwASBCQAAwASBCQAAwASBCQAAwASBCQAAwASBCQAAwASBCQAAwASBCQAAwASBCQAAwEQXXxcAtESdy1D2/jKVHqtWz6hwDUmMVXCQzddlAQA6CQIT/N6mvGJlbtit4opqd5vTHq5FY5KVnuL0YWUAgM6CS3Lwa5vyipXxQo5HWJKkkopqZbyQo015xT6qDADQmRCY4LfqXIYyN+yW0ci++rbMDbtV52qsBwAA1iEwwW9l7y9rcGbpTIak4opqZe8v815RAIBOKaAC05/+9CcNHz5cERERiomJadExhmHo/vvvl9PpVNeuXTVq1Ch9/fXX7VsoWqT0WNNhqS39AABoq4AKTCdPntTEiROVkZHR4mMefvhhPfHEE1q5cqWysrIUGRmp0aNHq7qaP8K+1jMq3NJ+AAC0VUA9JZeZmSlJWrNmTYv6G4ah5cuX67777tPYsWMlSWvXrlWvXr305ptv6oYbbmivUtECQxJj5bSHq6SiutH7mGySHPbTSwwAANCeAuoMU2vt379fJSUlGjVqlLvNbrdr6NCh2rFjhw8rgyQFB9m0aEyypNPh6Ez124vGJLMeEwCg3XXqwFRSUiJJ6tWrl0d7r1693PsaU1NTo8rKSo8X2kd6ilMrJg+Sw+552c1hD9eKyYNYhwkA4BV+f0luwYIFeuihh5rts2fPHiUlJXmpImnJkiXuy39of+kpTl2V7GClbwCAz/h9YJo7d65uvvnmZvv07du3TWM7HA5J0uHDh+V0/v+ZisOHD2vgwIFNHnfPPfdozpw57u3KykolJCS0qQa0THCQTan9uvu6DABAJ+X3gSkuLk5xcXHtMnZiYqIcDoe2bt3qDkiVlZXKyspq9km7sLAwhYWFtUtNAADA/wTUPUyFhYXKzc1VYWGh6urqlJubq9zcXB0/ftzdJykpSevXr5ck2Ww2zZo1Sw8++KDeeustffHFF5o6dari4+M1btw4H70LAADgb/z+DFNr3H///Xr++efd25dccokkafv27UpLS5Mk5efnq6Kiwt3nD3/4g06cOKHbbrtN5eXlGjFihDZt2qTwcNb2AQAAp9kMw+CHuM5SZWWl7Ha7KioqFB0d7etyAABAC7Tm73dAXZIDAABoDwQmAAAAEwQmAAAAEwQmAAAAEwQmAAAAEwQmAAAAEwQmAAAAEwQmAAAAEwQmAAAAEwQmAAAAEwQmAAAAEwQmAAAAEwQmAAAAEwQmAAAAEwQmAAAAEwQmAAAAEwQmAAAAEwQmAAAAEwQmAAAAEwQmAAAAEwQmAAAAEwQmAAAAEwQmAAAAEwQmAAAAEwQmAAAAEwQmAAAAEwQmAAAAEwQmAAAAEwQmAAAAEwQmAAAAEwQmAAAAEwQmAAAAEwQmAAAAEwQmAAAAEwQmAAAAEwQmAAAAEwQmAAAAEwQmAAAAEwQmAAAAEwQmAAAAEwQmAAAAEwQmAAAAEwQmAAAAEwQmAAAAEwQmAAAAEwQmAAAAEwQmAAAAEwQmAAAAEwQmAAAAEwQmAAAAEwQmAAAAEwEVmP70pz9p+PDhioiIUExMTIuOufnmm2Wz2Txe6enp7VsoAADoULr4ugArnTx5UhMnTlRqaqpWr17d4uPS09P13HPPubfDwsLaozwAANBBBVRgyszMlCStWbOmVceFhYXJ4XC0Q0UAACAQBNQlubZ699131bNnT11wwQXKyMjQkSNHfF0SAADwIwF1hqkt0tPTNX78eCUmJmrfvn2699579Ytf/EI7duxQcHBwo8fU1NSopqbGvV1ZWemtcgEAgA/4/RmmBQsWNLgp+6evvXv3tnn8G264Qddcc40uvPBCjRs3Tm+//bY+/fRTvfvuu00es2TJEtntdvcrISGhzfMDAAD/5/dnmObOnaubb7652T59+/a1bL6+ffuqR48eKigo0MiRIxvtc88992jOnDnu7crKSkITAAABzNLAVF5ertWrV2vPnj2SpJ///Oe65ZZbZLfb2zxmXFyc4uLirCrR1HfffacjR47I6XQ22ScsLIwn6QAA6EQsuyT3j3/8Q/369dNjjz2msrIylZWVadmyZerXr59ycnKsmqZZhYWFys3NVWFhoerq6pSbm6vc3FwdP37c3ScpKUnr16+XJB0/flx33323PvnkEx04cEBbt27V2LFj1b9/f40ePdorNQMAAP9n2Rmm2bNn65prrtEzzzyjLl1OD3vq1CndeuutmjVrlt5//32rpmrS/fffr+eff969fckll0iStm/frrS0NElSfn6+KioqJEnBwcHatWuXnn/+eZWXlys+Pl5XX321Fi9ezBkkAADgZjMMw7BioK5du+qzzz5TUlKSR/vu3bt16aWXqqqqyopp/FJlZaXsdrsqKioUHR3t63IAAEALtObvt2WX5KKjo1VYWNig/dtvv1VUVJRV0wAAAHidZYHp+uuv1/Tp0/XKK6/o22+/1bfffquXX35Zt956q2688UarpgEAAPA6y+5hWrp0qWw2m6ZOnapTp05JkkJCQpSRkaE///nPVk0DAADgdZbdw1SvqqpK+/btkyT169dPERERVg7vl7iHCQCAjqc1f78tX7gyIiJCKSkpkiSbzWb18AAAAF5n6U+jrF69WikpKQoPD1d4eLhSUlL017/+1copAAAAvM6yM0z333+/li1bppkzZyo1NVWStGPHDs2ePVuFhYV64IEHrJoKAADAqyy7hykuLk5PPPFEgyfiXnrpJc2cOVM//PCDFdP4Je5hAgCg4/HJOky1tbW69NJLG7QPHjzY/dQcAABAR2RZYJoyZYpWrFjRoH3VqlWaNGmSVdMAAAB4naVPya1evVrvvPOOhg0bJknKyspSYWGhpk6dqjlz5rj7LVu2zMppAQAA2pVlgSkvL0+DBg2SJPc6TD169FCPHj2Ul5fn7sdSAwAAoKOxLDBt377dqqEAAAD8iqXrMH3wwQeaPHmyhg8frqKiIknSf/3Xf+nDDz+0choAAACvsiwwvf766xo9erS6du2qnJwc1dTUSJIqKir0H//xH1ZNAwAA4HWWBaYHH3xQK1eu1DPPPKOQkBB3+z//8z8rJyfHqmkAAAC8zrLAlJ+fr8svv7xBu91uV3l5uVXTAAAAeJ1lgcnhcKigoKBB+4cffqi+fftaNQ0AAIDXWRaYZsyYobvuuktZWVmy2Ww6dOiQ1q1bp3nz5ikjI8OqaQAAALzOsmUFFixYIJfLpZEjR6qqqkqXX365wsLCNG/ePM2cOdOqaQAAALzOsh/frXfy5EkVFBTo+PHjSk5OVrdu3awc3i/x47sAAHQ8Xv/x3draWo0cOVJff/21QkNDlZycrCFDhnSKsAQAAAKfJYEpJCREu3btsmIoAAAAv2PZTd+TJ0/W6tWrrRoOAADAb1h20/epU6f07LPPasuWLRo8eLAiIyM99i9btsyqqQAAALzKssCUl5enQYMGSZK++uorj302m82qaQAAALzOssC0fft2q4YCAADwK5bdwwQAABCoCEwAAAAmCEwAAAAmCEwAAAAmCEwAAAAmLAtMP/74o6qqqtzbBw8e1PLly/XOO+9YNQUAAIBPWBaYxo4dq7Vr10qSysvLNXToUD366KMaO3asVqxYYdU0AAAAXmdZYMrJydG//Mu/SJL+9re/qVevXjp48KDWrl2rJ554wqppAAAAvM6ywFRVVaWoqChJ0jvvvKPx48crKChIw4YN08GDB62aBgAAwOssC0z9+/fXm2++qW+//Vb/+7//q6uvvlqSVFpaqujoaKumAQAA8DrLAtP999+vefPm6fzzz9fQoUOVmpoq6fTZpksuucSqaQAAALzOZhiGYdVgJSUlKi4u1sUXX6ygoNNZLDs7W9HR0UpKSrJqGr9TWVkpu92uiooKzqYBANBBtObvt2U/vitJDodDDofDo23IkCFWTgEAAOB1ZxWY5syZo8WLFysyMlJz5sxptu+yZcvOZioAAACfOavA9Nlnn6m2ttb93wAAAIHI0nuYOivuYQIAoOPx2T1MW7du1datW1VaWiqXy+Vut9lsWr16tZVTAQAAeI1lgSkzM1MPPPCALr30UjmdTtlsNquGBgAA8CnLAtPKlSu1Zs0aTZkyxaohAQAA/IJlC1eePHlSw4cPt2o4AAAAv2FZYLr11lv14osvWjUcAACA3zjrdZjquVwurVq1Slu2bNFFF12kkJAQj76swwQAADqqs16H6UwDBw6UJOXl5Xm0cwM4AADoyM4qMG3fvt2qOgAAAPyWZfcwFRYWqqk1MAsLC62apkkHDhzQ9OnTlZiYqK5du6pfv35atGiRTp482exx1dXV+v3vf6/u3burW7dumjBhgg4fPtzu9QIAgI7DssCUmJio77//vkH7kSNHlJiYaNU0Tdq7d69cLpeefvppffnll3rssce0cuVK3Xvvvc0eN3v2bG3YsEGvvfaa3nvvPR06dEjjx49v93oBAEDHYdlPowQFBenw4cOKi4vzaD948KCSk5N14sQJK6ZplUceeUQrVqzQN9980+j+iooKxcXF6cUXX9S1114r6XTwGjBggHbs2KFhw4a1aB5+GgUAgI7Hqz+NUv+knM1m08KFCxUREeHeV1dXp6ysLPfN4N5WUVGh2NjYJvfv3LlTtbW1GjVqlLstKSlJvXv3blVgAgAAge2sA1P9k3KGYeiLL75QaGioe19oaKguvvhizZs372ynabWCggI9+eSTWrp0aZN9SkpKFBoaqpiYGI/2Xr16qaSkpMnjampqVFNT496urKw863oBAID/OuvAVP+k3LRp0/T4449bfklqwYIFeuihh5rts2fPHiUlJbm3i4qKlJ6erokTJ2rGjBmW1iNJS5YsUWZmpuXjAgAA/2TZPUzt5fvvv9eRI0ea7dO3b1/3ma1Dhw4pLS1Nw4YN05o1axQU1PR97du2bdPIkSN19OhRj7NMffr00axZszR79uxGj2vsDFNCQgL3MAEA0IF49R6m9hYXF9fgRvKmFBUV6corr9TgwYP13HPPNRuWJGnw4MEKCQnR1q1bNWHCBElSfn6+CgsLlZqa2uRxYWFhCgsLa/mbAAAAHZplywr4WlFRkdLS0tS7d28tXbpU33//vUpKSjzuRSoqKlJSUpKys7MlSXa7XdOnT9ecOXO0fft27dy5U9OmTVNqaio3fAMAADe/P8PUUps3b1ZBQYEKCgp03nnneeyrv+pYW1ur/Px8VVVVufc99thjCgoK0oQJE1RTU6PRo0frL3/5i1drBwAA/s3v72HqCFiHCQCAjsdn9zBt3bpVW7duVWlpqVwul8e+Z5991sqpAAAAvMaywJSZmakHHnhAl156qZxOp2w2m1VDAwAA+JRlgWnlypVas2aNpkyZYtWQAAAAfsGyp+ROnjyp4cOHWzUcAACA37AsMN1666168cUXrRoOAADAb1h2Sa66ulqrVq3Sli1bdNFFFykkJMRj/7Jly6yaCgAAwKssC0y7du3SwIEDJUl5eXke+7gBHAAAdGSWBab6H+EFAAAINAHz0ygAAADt5azOMM2ZM0eLFy9WZGSk5syZ02xf7mECAAAd1VkFps8++0y1tbXu/24K9zABAICOjN+SswC/JQcAQMfTmr/f3MMEAABggsAEAABggsAEAABgwrLAVFpaatVQAAAAfsWywHTttdeqrq6u0X2nTp2yahoAAACva1Vgeumll5rcFxMTozvvvLNB+5EjRzRq1KjWVwYAAOAnWhSYSkpKNH78eG3evLnJPmvXrtXmzZv17LPPutv27NmjIUOGKDIy8uwrBQAA8JEWLVy5atUq1dbWeoShn4qJidHrr7+utLQ0paSk6OjRo7r++us1ffp0PfLII5YVDAAA4G0tWriyvLxcd911l44fP67XX3/d3T5+/HgNHDjQ/erdu7deeuklzZw5U9XV1XryySc1bdq0dn0D/oCFKwEA6Hha8/e7RWeYYmJi9Pzzz+t//ud/PNr79eunDz74QE899ZR++OEHnXPOObr44otlGIZuuukmDRo0SLW1tQoJCWn7uwEAAPAxy34apaioSLm5uR6vb775Rl26dFFSUpI+//xzK6bxS5xhAgCg47H8DFNLnHvuuTr33HP1q1/9yt12/Phx5ebmBnRYAgAAgY8f37UAZ5gAAOh4+PFdAAAAC1l2SQ44U53LUPb+MpUeq1bPqHANSYxVcJDN12UBANAmBCZYblNesTI37FZxRbW7zWkP16IxyUpPcfqwMgAA2oZLcrDUprxiZbyQ4xGWJKmkoloZL+RoU16xjyoDAKDtCEywTJ3LUOaG3WrsKYL6tswNu1Xn4jkDAEDHQmCCZbL3lzU4s3QmQ1JxRbWy95d5rygAACxAYIJlSo81HZba0g8AAH9BYIJlekaFW9oPAAB/QWCCZYYkxsppD1dTiwfYdPppuSGJsd4sCwCAs0ZggmWCg2xaNCZZkhqEpvrtRWOSWY8JANDhEJhgqfQUp1ZMHiSH3fOym8MerhWTB7EOEwCgQ2LhSlguPcWpq5IdrPQNAAgYBCa0i+Agm1L7dfd1GQAAWIJLcgAAACYITAAAACYITAAAACYITAAAACYITAAAACYITAAAACYITAAAACYITAAAACYITAAAACYITAAAACYITAAAACYITAAAACYITAAAACYCJjAdOHBA06dPV2Jiorp27ap+/fpp0aJFOnnyZLPHpaWlyWazebxuv/12L1UNAAA6gi6+LsAqe/fulcvl0tNPP63+/fsrLy9PM2bM0IkTJ7R06dJmj50xY4YeeOAB93ZERER7lwsAADqQgAlM6enpSk9Pd2/37dtX+fn5WrFihWlgioiIkMPhaO8SAQBABxUwl+QaU1FRodjYWNN+69atU48ePZSSkqJ77rlHVVVVXqgOAAB0FAFzhumnCgoK9OSTT5qeXbrpppvUp08fxcfHa9euXZo/f77y8/P1xhtvNHlMTU2Nampq3NuVlZWW1Q0AAPyPzTAMw9dFNGfBggV66KGHmu2zZ88eJSUlubeLiop0xRVXKC0tTX/9619bNd+2bds0cuRIFRQUqF+/fo32+eMf/6jMzMwG7RUVFYqOjm7VfAAAwDcqKytlt9tb9Pfb7wPT999/ryNHjjTbp2/fvgoNDZUkHTp0SGlpaRo2bJjWrFmjoKDWXXU8ceKEunXrpk2bNmn06NGN9mnsDFNCQgKBCQCADqQ1gcnvL8nFxcUpLi6uRX2Liop05ZVXavDgwXruuedaHZYkKTc3V5LkdDqb7BMWFqawsLBWjw0AADqmgLnpu6ioSGlpaerdu7eWLl2q77//XiUlJSopKfHok5SUpOzsbEnSvn37tHjxYu3cuVMHDhzQW2+9palTp+ryyy/XRRdd5Ku3AgAA/Izfn2Fqqc2bN6ugoEAFBQU677zzPPbVX3Wsra1Vfn6++ym40NBQbdmyRcuXL9eJEyeUkJCgCRMm6L777vN6/QAAwH/5/T1MHUFrroECAAD/0Jq/3wFzSQ4AAKC9EJgAAABMEJgAAABMEJgAAABMBMxTcoA/qHMZyt5fptJj1eoZFa4hibEKDrL5uiwAwFkiMAEW2ZRXrMwNu1VcUe1uc9rDtWhMstJTml4IFQDg/7gkB1hgU16xMl7I8QhLklRSUa2MF3K0Ka/YR5UBAKxAYALOUp3LUOaG3WpsQbP6tswNu1XnYskzAOioCEzAWcreX9bgzNKZDEnFFdXK3l/mvaIAAJYiMAFnqfRY02GpLf0AAP6HwAScpZ5R4Zb2AwD4HwITcJaGJMbKaQ9XU4sH2HT6abkhibHeLAsAYCECE3CWgoNsWjQmWZIahKb67UVjklmPCQA6MAITYIH0FKdWTB4kh93zspvDHq4VkwexDhMAdHAsXAlYJD3FqauSHaz0DQABiMAEWCg4yKbUft19XQYAwGJckgMAADBBYAIAADBBYAIAADBBYAIAADBBYAIAADDBU3I4K3Uug8foAQABj8CENtuUV6zMDbtVXPH/PyrrtIdr0ZhkFmoEAAQULsmhTTblFSvjhRyPsCRJJRXVynghR5vyin1UGQAA1iMwodXqXIYyN+yW0ci++rbMDbtV52qsBwAAHQ+BCa2Wvb+swZmlMxmSiiuqlb2/zHtFAQDQjghMaLXSY02Hpbb0AwDA3xGY0Go9o8It7QcAgL8jMKHVhiTGymkPV1OLB9h0+mm5IYmx3iwLAIB2Q2BCqwUH2bRoTLIkNQhN9duLxiSzHhMAIGAQmNAm6SlOrZg8SA6752U3hz1cKyYPYh0mAEBAYeFKtFl6ilNXJTtY6RsAEPAITDgrwUE2pfbr7usyAABoV1ySAwAAMEFgAgAAMEFgAgAAMEFgAgAAMEFgAgAAMEFgAgAAMEFgAgAAMEFgAgAAMEFgAgAAMEFgAgAAMEFgAgAAMMFvyQF+qs5l8MPGAOAnCEyAH9qUV6zMDbtVXFHtbnPaw7VoTLLSU5w+rAwAOicuyQF+ZlNesTJeyPEIS5JUUlGtjBdytCmv2EeVAUDnRWAC/Eidy1Dmht0yGtlX35a5YbfqXI31AAC0FwIT4Eey95c1OLN0JkNScUW1sveXea8oAACBCfAnpceaDktt6QcAsAaBCfAjPaPCLe0HALBGQAWma665Rr1791Z4eLicTqemTJmiQ4cONXtMdXW1fv/736t79+7q1q2bJkyYoMOHD3upYsDTkMRYOe3hamrxAJtOPy03JDHWm2UBQKcXUIHpyiuv1Kuvvqr8/Hy9/vrr2rdvn6699tpmj5k9e7Y2bNig1157Te+9954OHTqk8ePHe6liwFNwkE2LxiRLUoPQVL+9aEwy6zEBgJfZDMMI2Mdt3nrrLY0bN041NTUKCQlpsL+iokJxcXF68cUX3cFq7969GjBggHbs2KFhw4a1aJ7KykrZ7XZVVFQoOjra0veAzol1mACg/bXm73fALlxZVlamdevWafjw4Y2GJUnauXOnamtrNWrUKHdbUlKSevfu3arABFgtPcWpq5IdrPQNAH4i4ALT/Pnz9dRTT6mqqkrDhg3T22+/3WTfkpIShYaGKiYmxqO9V69eKikpafK4mpoa1dTUuLcrKyvPum7gp4KDbErt193XZQAA1AHuYVqwYIFsNluzr71797r733333frss8/0zjvvKDg4WFOnTpXVVx2XLFkiu93ufiUkJFg6PgAA8C9+fw/T999/ryNHjjTbp2/fvgoNDW3Q/t133ykhIUEff/yxUlNTG+zftm2bRo4cqaNHj3qcZerTp49mzZql2bNnNzpfY2eYEhISuIcJAIAOJKDuYYqLi1NcXFybjnW5XJLkEW7ONHjwYIWEhGjr1q2aMGGCJCk/P1+FhYWNBqx6YWFhCgsLa1NNAACg4/H7wNRSWVlZ+vTTTzVixAidc8452rdvnxYuXKh+/fq5w09RUZFGjhyptWvXasiQIbLb7Zo+fbrmzJmj2NhYRUdHa+bMmUpNTeWGbwAA4BYwgSkiIkJvvPGGFi1apBMnTsjpdCo9PV333Xef+2xQbW2t8vPzVVVV5T7uscceU1BQkCZMmKCamhqNHj1af/nLX3z1NgAAgB/y+3uYOgLWYQIAoONpzd9vv39KDgAAwNcC5pIccKY6l8GijwAAyxCYEHD4WREAgNW4JIeAsimvWBkv5HiEJUkqqahWxgs52pRX7KPKAAAdGYEJAaPOZShzw2419hRDfVvmht2qc/GcAwCgdQhMCBjZ+8sanFk6kyGpuKJa2fvLvFcUACAgEJgQMEqPNR2W2tIPAIB6BCYEjJ5R4Zb2AwCgHoEJAWNIYqyc9nA1tXiATaeflhuSGOvNsgAAAYDAhIARHGTTojHJktQgNNVvLxqTzHpMAIBWIzAhoKSnOLVi8iA57J6X3Rz2cK2YPIh1mAAAbcLClQg46SlOXZXsYKVvAIBlCEwISMFBNqX26+7rMgAAAYJLcgAAACYITAAAACYITAAAACYITAAAACa46duP1bkMnvQCAMAPEJj81Ka8YmVu2O3xY7JOe7gWjUlmLSEAALyMS3J+aFNesTJeyPEIS5JUUlGtjBdytCmv2EeVAQDQORGY/Eydy1Dmht0yGtlX35a5YbfqXI31AAAA7YHA5Gey95c1OLN0JkNScUW1sveXea8oAAA6OQKTnyk91nRYaks/AABw9ghMfqZnVLh5p1b0AwAAZ4/A5GeGJMbKaQ9XU4sH2HT6abkhibHeLAsAgE6NwORngoNsWjQmWZIahKb67UVjklmPCQAALyIw+aH0FKdWTB4kh93zspvDHq4VkwexDhMAAF7GwpV+Kj3FqauSHaz0DQCAHyAw+bHgIJtS+3X3dRkAAHR6XJIDAAAwQWACAAAwQWACAAAwQWACAAAwQWACAAAwQWACAAAwQWACAAAwQWACAAAwwcKVFjAMQ5JUWVnp40oAAEBL1f/drv873hwCkwWOHTsmSUpISPBxJQAAoLWOHTsmu93ebB+b0ZJYhWa5XC4dOnRIUVFRstn4rTfpdGpPSEjQt99+q+joaF+X0+nw+fsWn7/v8R34Vkf5/A3D0LFjxxQfH6+goObvUuIMkwWCgoJ03nnn+boMvxQdHe3X/2MJdHz+vsXn73t8B77VET5/szNL9bjpGwAAwASBCQAAwASBCe0iLCxMixYtUlhYmK9L6ZT4/H2Lz9/3+A58KxA/f276BgAAMMEZJgAAABMEJgAAABMEJgAAABMEJgAAABMEJljiT3/6k4YPH66IiAjFxMS0+Lg9e/bommuukd1uV2RkpC677DIVFha2X6EBrK3fQb3bb79dNptNy5cvt7y2zqC1n39tba3mz5+vCy+8UJGRkYqPj9fUqVN16NCh9i82ALXl379hGLr//vvldDrVtWtXjRo1Sl9//XX7FhrAysrKNGnSJEVHRysmJkbTp0/X8ePHmz2mpKREU6ZMkcPhUGRkpAYNGqTXX3/dSxW3DoEJljh58qQmTpyojIyMFh+zb98+jRgxQklJSXr33Xe1a9cuLVy4UOHh4e1YaeBqy3dQb/369frkk08UHx/fDpV1Dq39/KuqqpSTk6OFCxcqJydHb7zxhvLz83XNNde0c6WBqS3//h9++GE98cQTWrlypbKyshQZGanRo0erurq6HSsNXJMmTdKXX36pzZs36+2339b777+v2267rdljpk6dqvz8fL311lv64osvNH78eF133XX67LPPvFR1KxiAhZ577jnDbre3qO/1119vTJ48uX0L6oRa8x0YhmF89913xrnnnmvk5eUZffr0MR577LF2q60zaO3nf6bs7GxDknHw4EFri+pEWvr5u1wuw+FwGI888oi7rby83AgLCzNeeumldqwwMO3evduQZHz66afutr///e+GzWYzioqKmjwuMjLSWLt2rUdbbGys8cwzz7RbrW3FGSb4hMvl0saNG/VP//RPGj16tHr27KmhQ4fqzTff9HVpnYrL5dKUKVN099136+c//7mvy+n0KioqZLPZ2nRJFa2zf/9+lZSUaNSoUe42u92uoUOHaseOHT6srGPasWOHYmJidOmll7rbRo0apaCgIGVlZTV53PDhw/XKK6+orKxMLpdLL7/8sqqrq5WWluaFqluHwASfKC0t1fHjx/XnP/9Z6enpeuedd/Rv//ZvGj9+vN577z1fl9dpPPTQQ+rSpYvuvPNOX5fS6VVXV2v+/Pm68cYb/f7HSgNBSUmJJKlXr14e7b169XLvQ8uVlJSoZ8+eHm1dunRRbGxss5/nq6++qtraWnXv3l1hYWH67W9/q/Xr16t///7tXXKrEZjQpAULFshmszX72rt3b5vGdrlckqSxY8dq9uzZGjhwoBYsWKBf//rXWrlypZVvo0Nrz+9g586devzxx7VmzRrZbDaLKw8M7fn5n6m2tlbXXXedDMPQihUrLKg8MHjr80fT2vs7WLhwocrLy7Vlyxb94x//0Jw5c3Tdddfpiy++sPBdWKOLrwuA/5o7d65uvvnmZvv07du3TWP36NFDXbp0UXJyskf7gAED9OGHH7ZpzEDUnt/BBx98oNLSUvXu3dvdVldXp7lz52r58uU6cOBAm8YNJO35+derD0sHDx7Utm3bOLt0hvb8/B0OhyTp8OHDcjqd7vbDhw9r4MCBbRozELX0O3A4HCotLfVoP3XqlMrKytyf9U/t27dPTz31lPLy8ty3BFx88cX64IMP9J//+Z9+93+eCUxoUlxcnOLi4tpl7NDQUF122WXKz8/3aP/qq6/Up0+fdpmzI2rP72DKlCke929I0ujRozVlyhRNmzatXebsaNrz85f+Pyx9/fXX2r59u7p3795uc3VE7fn5JyYmyuFwaOvWre6AVFlZqaysrDY9aRqoWvodpKamqry8XDt37tTgwYMlSdu2bZPL5dLQoUMbPaaqqkqSFBTkebErODjYfRXCn3BJDpYoLCxUbm6uCgsLVVdXp9zcXOXm5nqswZGUlKT169e7t++++2698soreuaZZ1RQUKCnnnpKGzZs0O9+9ztfvIUOr7XfQffu3ZWSkuLxCgkJkcPh0AUXXOCrt9Fhtfbzr62t1bXXXqt//OMfWrdunerq6lRSUqKSkhKdPHnSV2+jw2rt52+z2TRr1iw9+OCD7kfap06dqvj4eI0bN85H76LjGjBggNLT0zVjxgxlZ2fro48+0h133KEbbrjBvVxJUVGRkpKSlJ2dLen099G/f3/99re/VXZ2tvbt26dHH31Umzdv9s/vwNeP6SEw/OY3vzEkNXht377d3UeS8dxzz3kct3r1aqN///5GeHi4cfHFFxtvvvmmdwsPIG39Ds7EsgJt19rPf//+/Y32/+kxaJm2/Pt3uVzGwoULjV69ehlhYWHGyJEjjfz8fO8XHyCOHDli3HjjjUa3bt2M6OhoY9q0acaxY8fc++v/zZ/5nXz11VfG+PHjjZ49exoRERHGRRdd1GCZAX9hMwzD8EIuAwAA6LC4JAcAAGCCwAQAAGCCwAQAAGCCwAQAAGCCwAQAAGCCwAQAAGCCwAQAAGCCwAQAAGCCwAR0IkePHlVmZqaKi4t9XQoAdCis9A10IlOnTtWRI0cUEhKiN99809flAECHwRkmoJPYuHGjjh07po0bNyomJkbr1q3zdUkA0GFwhgkAAMAEZ5gAAABMEJgAAABMEJiATmLevHkaN25ci/qmpaVp1qxZ7VqPN+fxNy153/7+2fh7fYDVuvi6AADekZubqxEjRrSo7xtvvKGQkJB2rqhjSUtL08CBA7V8+XKvzGfld7Bjxw6NGDFC6enp2rhxoyVjnlnfu+++qyuvvLLJvmlpadq+fbsl8wK+whkmoJP4/PPPNXDgwBb1jY2NVVRUVPsWhGZZ+R2sXr1aM2fO1Pvvv69Dhw412/fkyZMtGvPM+oYPH67i4uIGr6efflo2m02/+93vzvo9AL5GYAI6ge+++04//PCDOzCVl5drzJgxGjFihEpKShr0/+nllrS0NM2cOVOzZs3SOeeco169eumZZ57RiRMnNG3aNEVFRal///76+9//7nHMHXfcoTvuuEN2u109evTQwoUL1dyDuS6XS0uWLFFiYqK6du2qiy++WH/7298a1NbaWlo67p133qk//OEPio2NlcPh0B//+EdJ0s0336z33ntPjz/+uGw2m2w2mw4cONCisU+cOKGpU6eqW7ducjqdevTRR5v9rpr7DpqqrznHjx/XK6+8ooyMDP3qV7/SmjVrGsxzxx13aNasWerRo4dGjx6tVatWKT4+Xi6Xy6Pv2LFjdcsttzSoLzQ0VA6Hw+N19OhRzZs3T/fee68mTpzYovcM+DMCE9AJ5ObmKiYmRueff76++OILXXbZZTr33HO1fft2ORyOFo3x/PPPq0ePHsrOztbMmTOVkZGhiRMnavjw4crJydHVV1+tKVOmqKqqyuOYLl26KDs7W48//riWLVumv/71r03OsWTJEq1du1YrV67Ul19+qdmzZ2vy5Ml67733zqqW1owbGRmprKwsPfzww3rggQe0efNmPf7440pNTdWMGTPcZ08SEhJaNPbdd9+t9957T//93/+td955R++++65ycnJa9Jk39h00Vl9zXn31VSUlJemCCy7Q5MmT9eyzzzYIrc8//7xCQ0P10UcfaeXKlZo4caKOHDnicRmtrKxMmzZt0qRJk0zrLC8v19ixY5WWlqbFixe36b0CfscAEPAWL15sXHHFFca6deuMc845x1i1alWz/a+44grjrrvu8tgeMWKEe/vUqVNGZGSkMWXKFHdbcXGxIcnYsWOH+5gBAwYYLpfL3Wf+/PnGgAEDGp2nurraiIiIMD7++GOPWqZPn27ceOONba6lreMahmFcdtllxvz58xv9TFpS87Fjx4zQ0FDj1Vdfde87cuSI0bVr1wZj/ZTZd/DT+poyfPhwY/ny5YZhGEZtba3Ro0cPY/v27R7jXnLJJQ2OGzt2rHHLLbe4t59++mkjPj7eqKura7S+enV1dcYvfvELY8CAAUZlZWWztQEdCTd9A51Abm6udu3apTvuuEMbN25Uampqq8e46KKL3P8dHBys7t2768ILL3S39erVS5JUWlrqbhs2bJhsNpt7OzU1VY8++qjq6uoUHBzsMX5BQYGqqqp01VVXebSfPHlSl1xySZtraeu4kuR0Oj3ez0+Zjb1v3z6dPHlSQ4cOde+LjY3VBRdc0OSYzWltffn5+crOztb69eslSV26dNH111+v1atXKy0tzd1v8ODBDY6dNGmSZsyYob/85S8KCwvTunXrdMMNNygoqPkLE/fee6927Nih7Oxs7oNDQCEwAZ1Abm6uxo8frxdffFHl5eVtGuOnT2zZbDaPtvpg9NP7Xlrq+PHjkk7/hMu5557rsS8sLKzNtZztuM29H7Oxy8rKmjy2LVpb3+rVq3Xq1CnFx8e72wzDUFhYmJ566inZ7XZJUmRkZINjx4wZI8MwtHHjRl122WX64IMP9NhjjzVb38svv6ylS5dq48aN+tnPftaatwb4PQITEOCOHTumb775Rq+99pqGDx+uG264QR9//LF+/vOft/vcWVlZHtuffPKJfvaznzU4uyRJycnJCgsLU2Fhoa644grLarBq3NDQUNXV1bVq7HPOOUchISHKyspS7969JUlHjx7VV199Zel7bMypU6e0du1aPfroo7r66qs99o0bN04vvfSSbr/99iaPDw8P1/jx47Vu3ToVFBToggsu0KBBg5rsn5ubq+nTp+vPf/6zRo8ebdn7APwFgQkIcJ9//rmCg4OVnJysSy65RHl5eRozZoyys7PVo0ePdp27sLBQc+bM0W9/+1vl5OToySefbPIpsaioKM2bN0+zZ8+Wy+XSiBEjVFFRoY8++kjR0dH6zW9+06YarBr3/PPPV1ZWlg4cOKBu3bq5H6s3G3v69Om6++671b17d/Xs2VP//u//bnpZywpvv/22jh49qunTp7vPJNWbMGGCVq9e3Wxgkk5flvv1r3+tL7/8UpMnT26y3w8//KBx48YpLS1NkydPbvDkZXBwsOLi4tr+ZgA/QGACAlxubq6SkpLcl58eeeQR7dmzR+PHj9eWLVsUGhrabnNPnTpVP/74o4YMGaLg4GDddddduu2225rsv3jxYsXFxWnJkiX65ptvFBMTo0GDBunee+89qzqsGHfevHn6zW9+o+TkZP3444/av3+/zj//fNOxH3nkER0/flxjxoxRVFSU5s6dq4qKirN6Py2xevVqjRo1qkFYkk4Hpocffli7du1qdox//dd/VWxsrPLz83XTTTc12W/jxo06ePCgDh48KKfT2WB/nz593MswAB2VzTCaWRQFANrI2ytjA0B7Yh0mAAAAEwQmAAAAE1ySAwAAMMEZJgAAABMEJgAAABMEJgAAABMEJgAAABMEJgAAABMEJgAAABMEJgAAABMEJgAAABMEJgAAABMEJgAAABMEJgAAABMEJgAAABP/BwRSOHPPEilhAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(khats_az, khats)\n",
    "plt.xlabel(\"$\\hat k$ implemented in ArviZ\")\n",
    "plt.ylabel(\"$\\hat k$ in this repo\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b977452c-262b-47b1-a37b-e0aa39dbb8ae",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.14"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}