{
  "cells": [
    {
      "attachments": {
        "generated/multilevel-models/complete-pooling.png": {
          "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAF2CAYAAAB6XrNlAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90\nbGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAO\nxAAADsQBlSsOGwAAJ05JREFUeJzt3QmcjWX/x/FLdtkKJUpIhCwJlUePSqEo7UqlhWwt6NFKG1kq\n9RRZSjSlelDSIioq0pOypKzRvuixF4kyOP/X9/q/7vs1w8yY5T7nus99Pu/X67zmnJkxc03xneu+\n7uv6/QrFYrGYAQBExiGuBwAACBbBDgARQ7ADQMQQ7AAQMQQ7AEQMwQ4AEUOwA0DEEOwAEDEEOwBE\nDMEOABFDsANAxBDsABAxBDsARAzBDgARQ7ADQMQQ7AAQMQQ7AEQMwQ4AEUOwA0DEEOwAEDEEOwBE\nDMEOABFDsANAxBDsABAxBDsARAzBDgARQ7ADQMQQ7AAQMQS7Q7169TKFChUymzdvPuBjZ5xxhqlY\nsaKTcQFIbgS7Q19++aWpWrVqlgG+bNky06hRIyfjApDcCHZHYrGYWb58eZbh/eOPP5rffvuNYEeo\nLV682F5xHuzxxRdfuB5qyiniegCp6ttvvzU7duwwjRs3znImLw0bNnQwMiB3ihYtau6//37/9fTp\n0+2V5p133mlKlChh36dgb9CggcNRpiaC3REvvLOalef0MSAs9Pcz49/R2bNnm/Lly5vhw4c7HRdY\ninEmp/DWpWuRIkVMvXr1HIwMyJ+VK1ea+vXrux4GCHa3wV6yZElz/PHHH/CxhQsXmjp16pjixYs7\nGRuQVz/99JPZtm2bOfHEE10PBQS722CvXLmyOeSQQw64IfXLL7+wDIOkoo0AQrCHA8HuwO+//253\nvqxbt84+PJrx9O7d2z5nGQbJZMWKFfYtwR4O3Dx1QDsHpFy5cqZZs2bmwgsvNDt37jTvvPOOKVWq\nlL/DoGbNmubKK690PFrg4NauXWvf1q1b1/VQwIzd7Y3T8ePHm+bNm5u0tDTz3nvvma5du9q3OrS0\nceNGu1QDJANdbUqZMmVcDwXM2N0Ge6tWrUzHjh0P+LjW2IFkUqtWLftWV5+nnnqq6d+/vylbtqzr\nYaWsQjEdgURCafllw4YNdicBEAWbNm0yXbp0MXPnzrUbAv74448DNgYgcQj2BNu7d68pXbq0ad26\ntZkxY4br4QCIIH6lOrjJ9Ndff1EuAEDcMGMHgIhhxg4AEUOwA0DEEOwAEDEEOwBEDMEOABFDsANA\nxBDsABAxBDsARAzBDgARQ7ADQMQQ7AAQMQQ7AEQMwR6B/qlTp051PQxExNKlS83ChQtdDwMFRLAn\nsUWLFtnmwZ06dTLvvvuu6+Egye3atcv+XfI6IO3bt8/1kJBPBHsS2759u1m3bp193r17d9u1Bsiv\nBx980Hz99ddGlbzVnpEOSMmLeuxJ7sYbbzTPPvusfX7zzTebUaNGuR4SktCSJUvMKaecYjt8VahQ\nwaxatcocccQRroeFfCLYI7DGXr9+ffPrr7+aQoUKmY8++si0bNnS9bCQRNLT020fXq/J+osvvmiu\nuuoq18NCAXCtleTKly9vxo4da5/rd3S3bt1s6z0gtx555BE/1Nu3b286d+7sekgoIGbsEXHFFVeY\nKVOm2Od33323GTp0qOshIQmsXr3aNG7c2OzevduUKVPGrFy50hxzzDGuh4UCItgjYuPGjaZevXpm\ny5YtpnDhwnbLWpMmTVwPCyGm9fTTTz/dLFiwwL7WlV/Pnj1dDwsBYCkmInSj68knn/T/wXbt2tWu\nnQLZGT16tB/qrVq1sjurEA3M2CNE/ys7dOhgZs6caV8PGTLE3HPPPa6HhRD64Ycf7E33nTt3mhIl\nSphly5aZ448/3vWwEBCCPWJ+/vln+w9We9qLFStmb4qdcMIJroeFENE/+TZt2pg5c+bY148++qg9\nkIToYCkmYnTjS7scRDfEtCTDCUJklJaW5od606ZNTd++fV0PCQFjxh5BCvIzzzzT7mmXkSNHmltu\nucX1sBAC//vf/+xNdp1/KFKkiPn8889NgwYNXA8LAWPGHkE6Cq7TqFo79bY/ak0VqU1zuN69e9tQ\nF91/IdSjiWCPKN0IGzRokH3+559/2h0PXJyltmnTppnXX3/dPtesnRvr0cVSTITt2bPHnHbaaWbx\n4sX29XPPPWeuu+4618OCAzrfoDDXeQeVntA2R9WGQTQxY48wraFOmDDBvpV+/fqZ9evXux4WHLjt\ntttsqItulhLq0UawR1zDhg3tGrtobfWmm25yPSQk2DvvvGNeeOEF+7xmzZpm8ODBroeEOGMpJgX8\n/ffftryASrHKq6++ai655BLXw0IC6DyDmrH89NNP9vX7779vzjrrLNfDQpwxY08BxYsXt0syWlsV\nzdq3bt3qelhIAF2teaGuyp+Eemog2FOE2p316dPHPt+wYYNdc0W0zZ8/39aDkSpVqtgTpkgNLMWk\nEG171L7l77//3l97bdu2rethIU79Sxs1amRb3ckbb7xhLrjgAtfDQoIwY08hhx56qBk/frz/mj6p\n0aUzDF6oq1Y/oZ5amLGnIPqkRpvKBDRv3pz+pSmMYE9B2vaowyqqG+KtxdInNRroXwphKSbF+6SK\nKkDSJzUa6F8KYcaewjL2Sb3rrrvMsGHDXA8JBUD/UngI9hRGn9To9i8dN26c6dGjh+thwRGWYlIY\nfVKj279UN8iRupixpzj6pCY/+pdifwQ76JOaxOhfiqywFIMs+6RqaQbhR/9SZIUZOyz6pCYf+pci\nO8zYYdEnNbloPqYqnfQvRVYIdvjok5pc/UunT59un+v+CDe8kRFLMTigT6pK/C5ZssS+pk9q+NC/\nFAfDjB2ZaK124sSJmfqkejVlEA70L8XBEOw4aJ9UVYBEONC/FLnBUgyyRJ/U8KF/KXKLGTuyRJ/U\ncPcvVckAQh3ZIdiRLfqkhrd/qXegDMgKSzHIU5/UWbNmmXbt2rkeVkqhfynyihk78tQnVaVg6ZOa\nWPQvRV4xY0ee+6Rqvf2pp55yPaSU7F+qZhqVKlVyPSyEHMGOfPVJVU0ZNXZA4vqXvvTSS7S6Q66w\nFIN89Unt1q0bfVIT3L/0yiuvdD0kJAlm7MgT+qQmBv1LURAEO/KEPqnxR/9SFBRLMcgT+qTGH/1L\nUVDM2FHgPqkPPfSQGTBggOthRbJ/6fLly02tWrVcDwtJhmBHvtAnNXj6p9i2bVsze/Zs+5r+pcgv\nlmKQL/RJjU//Ui/Utc2R/qXIL2bsyDf6pAaH/qUIEjN25Bt9UoNB/1IEjWBHgdAnteDoX4qgsRSD\nAqNPav6pxn3dunXt+QBdAX3yySe0ukOBMWNHgdEnNf/034r+pQgawY5A0Ce14P1LvSUtoKBYikFg\n6JOae/QvRTwxY0dg6JOae/QvRTwR7AgUfVIPjv6liDeWYhA4+qTm3L9U5XjXrl1rX9O/FPHAjB2B\no09q9nSD1At1+pciXpixI27UZUlr7qJdMqNGjTKpjP6lSBSCHQnrk6q15ZYtW5pUpJr1CvUvvvjC\nvqZ/KeKJpRgkrE+qKkCmap9UleD1Qp3+pYg3ZuyIu1Tvk0r/UiQawY6405F51UPRnvZU65NK/1K4\nwFIM4i6V+6SOGTOG/qVIOGbscNIndciQIZEvT6va9CoboH399C9FIhHsSJhU6pNK/1K4xFIMEiaV\n+qRm7F/atGlT+pcioZixI6FSoU8q/UvhGjN2JFQq9EnVKVv6l8Ilgh0JF+U+qapB/9prr9nn9C+F\nKyzFwIko9knN2L9UNem1zZFWd3CBGTuc0NqzCoRFqU8q/UsRFgQ7nGnUqFFk+qTu37908ODBroeE\nFMZSDELfJ1WFw7Zv325PsLqkLZpabqlcuXKm99O/FGHDjB2h7ZOqOYfW3o866ihTtWpVf4ukC9pv\nr7K7Gss111xjtmzZkmX/UtWgJ9ThnGbsgGt9+/bVlaN9XHvttbHvv/8+ds455/jv06N///7Oxrdm\nzZpMY6lUqVJsypQpsXnz5vnvq1KlSuy3335zNkbAw1IMQtknVfvc96/d3qlTJzN58mQn49Pyytln\nn51lG0CNXehfirBgKQahoIAcOHCg/zqrhhyqNeNKdt/bC/VmzZqZ888/P8GjArJGsMM5lfAdOnSo\n6dWrV46fF8Zg9yxatMicc8455rvvvkvYmIDsEOxwShUedVNywIABdtdJTn799VdnRcNy80tFyzVa\nTlLteVY44RLBDqdbHbVu7fUCPRiFuqtDTL/88kuuPm/nzp32cNIrr7wS9zEB2SHY4bQgWKlSpfL0\nZ1wtx+T1++qeAeAKwQ5nihYtaj7++GNz2WWXRSbYq1evbl5//XXTvn37uI8JyA7BDufNN6ZOnWo+\n+OADe3ozXsGuNW/Vgs8PnSzdtm1bjp9TsmRJW0ZAJ2g7duyYr+8DBOX/KzABjqn5xtKlS824cePM\nvffe69czz0uwK7xV210VIxcvXmzfqs+oyhFo+6Q+rqsEhbDquZx88sn+o2HDhn6N+Lx8T29/vVrf\n6ZcUEAr+USUgJDZu3Bjr3r17rFChQplOe+qh06j7++qrr2J9+vSJVaxY8YDPz+2jSJEisTZt2sSm\nT58eS09Pz/T1Z82aleWfadCgQWzu3LkJ/C8D5A4nTxFaaimntnmffPKJ/74KFSqYzZs323rub775\nphkzZozdZpgVNczW9kMVD9MsvXDhwnbmrqWVlStXmk2bNmX5544++mjTo0cPW/dFBb9uv/12M2LE\nCP/jhx12mF120ed4ZYeBMCHYEWr66/niiy+aG2+80W6PVPel0aNH21D1yg94qlSpYo/0e8sr6mCk\ncM/u62oLo5Zr9Jg3b56ZP39+ps/Rss1tt91mg97ry6obvfplUrFixTj+1EDBEOxICrp5mZaWZve8\n621GrVu3Nr1797ZH+hXG+aUbn2PHjjXPP/+8ndV71BWpa9eu5rTTTjMtWrQo0M8BJALBjqQwe/Zs\nuzTilccVtdK78847zQknnBDo99qxY4eZNGmSvYnrlefVnnstyTzwwAPZ3mQFwoJgR+gNHz7c77Qk\nWhp59tlnTdu2beP6fTds2GCvBLzm1F6xr5kzZ7IUg1BjHztCS3OOu+66K1Ooa9a+YsWKuIe6HHnk\nkbajk0oF66atV+yrVatWZt26dXH//kB+MWNHaCnQNVsXrZ3rJurll1/uZCyqUXPuuefaomVSp04d\n29HJdbs+ICsEO0Lp4YcftrN10VbF6dOnJ2SWnhMdmjrvvPPMggUL7OvGjRubuXPnmnLlyjkdF7A/\ngh2h8+GHH/p9QzVTf+utt5yHukenWHVKVnvs5aqrrrJXEkCYEOwIFe1I0aEilQaQ//znP+aKK64w\nYaKDTaeeeqrfVENFv6gPgzDh5ilCRcsvXqh36dIldKEulSpVsnvpCxUqZF/rsJS3LRIIA4IdoVqC\n0alSOeqoo8wTTzxhwur000/3T6NqW2SfPn1cDwnwsRSDUFBJXZ3wXLt2rX2tdfUOHTqYMFMj60aN\nGplvv/3WvlbpYa2/A64xY0covPfee36oa/kl7KHudUkaP368/3rUqFFOxwN4mLEjFFS8S7N0UV12\nbSVMFjqNqvrvKjvw448/2pOxgEvM2OGcbpbOmDHDPlehrWQKdVHZAW856ZlnnnE9HIBgh3sKQ+/C\n0QvJoH311Vf21KpKA5QtW9Yu9Wh2HQR1UFKNdtHSzO7duwP5ukB+Eexw7uWXX7ZvVVjr0ksvDfzr\nqxFHkyZNzGeffWYrQqrRtAp5tWnTJpAQLlWqlLn++uvt8/Xr19vTqIBLBDuc2rhxoz9z1unSoEvi\nql/pJZdcYpd3VDzsscces4ee+vXrZ2/WvvHGG4HdI/CoUBjgEsEOp9S9yKOuR0FTDfWdO3faY/9l\nypTJtHwiGdvuFcRJJ53kH1jK+DMBLhDscCpjCDZt2jTwm7KvvPKKufjii03NmjUzfezwww/3rxiC\noHX72rVrp2Swax//DTfcYKpVq2Zr++gX3P4P7XRC4tCJF055Iah//Jr1BmnKlCl2p0rnzp0P+Jj6\np0rx4sUD+3664lizZo3t8qSG26nQjEPF0HQoS4e1LrroItuTduHChX6D8erVq9uwr1evnuuhphSC\nHU6tXLnSvlUglC5dOvBDT6Im1V41Ro8aWXsNsIOiG7TejWCt559xxhkmyvTLUTuNdu3aZYNcDUg8\nukpSqeWnnnrK3qxGYhHscMprGq3CWkHS9knNHGXEiBHZfl6Q/VIzztBVpTLq1BdWyzADBw7MFOoZ\ng12/UAn2xGONHU799ddf9m3Qu2E0I1e4qkKkQn7/h7c8owNRQVFDEI9msVGn+xdaQsvq7IF39bV3\n714HIwPBDqe0Bi7ejpIgW9l5VSKz+p5aOtAN1eOOOy6w76mSAp5UCLRPP/3U1K9fP8v/xl5P2Bo1\najgYGQh2OOXN1L2bmUFJT0+3b4sVK3bAxxTqKrWb1U3Vgsg4S884e4+irVu32m5SlStXzvLj77zz\njn0b9fsMYUWwwykvAL219qB4a/ZZbWccOnSoPS3q1VMPSsafIerB7l1hKdz3t3r1anvjWrtljj32\nWAejA8EOp7ylEIWBN8sOQq1atcwRRxxh3nzzzUw3MocMGWKP/OutPh6k5cuXZ/r+UabaOKpiqe2q\n33zzjf/+bdu22fsaWu5SQ3K4QbDDKe+0qZZiVq1aFeh692233WbX2nWD9M477zTnnHOO3cHRvXv3\nuHQ88vbkK/RSYW1Z/011L+Ef//iH6du3r7n11lvtLqMvvvjCFkNTOWM4onrsgCuTJ09WWUf7mDBh\nQqBfe9++fbGHHnoodvTRR8dKlCgRa9KkSez555+PxcPu3btjxYsXtz9H69atY6li9OjRsbp168aK\nFSsWq1ChQuziiy+OLVq0yPWwUh6NNuCULuN1OEl69eplxowZY5KRZqneydk77riDZQg4xVIMnK+x\nly9f3m9mnazzDPU7jWcxMyAvCHY4313hnUxUM4yPPvrIJBvdKBw3bpxfe6Z169auh4QUR7DDuYwn\nF5NxKUaz9a+//tpvxK0uTYBLBDuc066VRo0a2eevvfaaf2o0WWT8ZXTTTTc5HQsgBDucy1hvZM+e\nPf6yRjJQ9yevC5PqybPFD2FAsCMUdLy/XLly9vkjjzxi29aFnW709uzZ0693E69G3EBeEewIBVUD\nHDRokF/xUc2hw15IKy0tza+Jop6qV199teshARb72BEamvmqaJQaY4gaT+v0aBipLLAqG6pWSpEi\nRczixYv9+wSAa8zYERoqAzBx4kS/gNaAAQMCLTMQFF1JdOvWzS+Ade+99xLqCBWCHaGi4lnDhw/3\nl2TatWtne4iGhS5wtfPl3Xff9Zdg7r77btfDAjJhKQahXJLp2LGjmTFjhn1du3Ztu1e8atWqTsel\nfyr9+/c3jz/+uH2tm73//e9/7ZIMECYEO0Jp586d5txzz/VPoqrb/Zw5cwLteJTX5RftennmmWfs\nay0XzZ4921Y2BMKGpRiEkhphvPXWW6ZFixb29Q8//GBOOeUUM3XqVCd71bUk5IW6xqa964Q6wopg\nR2iVLVvWduJp06aNfb1lyxbTqVMnc9lll2XZGSlouph9+umnzYknnmivFrzlF83UVds9iriAjwaC\nHaF26KGH2pn77bff7jeLfvXVV+269ssvvxy3ve46IKXw1gEkrwOTTpYuWLDAv4qIomHDhpkePXpk\n2fIOScRxPXgg1xYsWBCrU6eO35hDj+rVq8eGDx8e27hxY4G//p49e2JvvPFGrG3btpm+h5pIDB06\nNJaenh6LslWrVtmfVT9zvXr1Iv/zRhk3T5FUdu3aZe6//357eMk7yi/FihWzSzQXX3yxnVkfc8wx\nfsPlnGg2vnTpUjNv3jzbzm3/rZX6WjphGvWdL7ryOf300+0ViYwdO9ZerSA5EexISitWrDCjR482\nkyZNMn/++ecBH69YsaJteNGkSRPbtFq7WAoXLmz3xivM9efVo3TNmjVZriv/85//tB2dLr30Unuy\nNOpGjhzp94Ft1aqV3V7qLX0h+RDsSGrbtm2z4a7SuatXry5wvZouXbrYQNcN01ShHUe6ItEW0xIl\nSphly5b57QqRnAh2RIL+GmtJZdGiRbZui2bjmpWnp6dn+2e0XKNZvfdo2bKlKVOmjEm1/25t27a1\nO328ypq6UY3kRrAjsv7++2/bbk/LLdom6XU46tu3r6lRo4Zdokl1zz33nLnhhhsy7fpJhaWnqOP/\nICJL/UdVnCvjLLxatWr2oBOM7VTlVc9UmE+YMIFQjwjujgApyCtm9vvvv9vX99xzj2nYsKHrYSEg\nBDuQgqZNm2amT59un9erV88GO6KDYAdSzNatW/2m29rrryUYLVshOgh2IMX069fPr7WjG8mnnnqq\n6yEhYAQ7kELUo/WFF16wz2vWrGkGDx7sekiIA4IdSBF//PGHLfDlUQkFFVlD9BDsQIpQCz+vFo56\ntp511lmuh4Q4IdiBFDB//nxbW0eqVKliHn30UddDQhwR7EAKVMTUDN2jyo3ly5d3OibEF8EORNyg\nQYNs4xBRaYULLrjA9ZAQZwQ7EGGff/65v+xSoUIFW54X0UewAxGlypYq8OW1D3zyyScpfJYiCHYg\nojRT//LLL+3z8847z3Tu3Nn1kJAgBDsQQWo68uCDD9rnqm45bty4XLUKRDQQ7EDEaOmla9euZvfu\n3X7zDDUVQeog2IGI0X51rym1+pd2797d9ZCQYAQ7ELH+pTphKupfqrIBNKVOPfwfByLUPEOzczWl\n9vav05Q6NRHsQESkpaX5TanVv1TleZGaCHYgAuhfiowIdiDJ0b8U+yPYgSRH/1Lsj2AHkhj9S5EV\ngh1IYvQvRVYIdiBJ0b8U2SHYgSRE/1LkhGAHkhD9S5ETgh1IMvQvxcEQ7EASoX8pcoNgB5II/UuR\nGwQ7kCToX4rcItiBJED/UuQFwQ4kAfqXIi8IdiDk6F+KvCLYgRCjfynyg2AHQoz+pcgPgh0IKfqX\nIr/4WwKEEP1LURAEOxBC9C9FQRDsQMjQvxQFRbADIUL/UgSBYAdChP6lCALBDoQE/UsRFIIdCAn6\nlyIoBDsQAvQvRZAIdsAx+pciaAQ74Bj9SxE0gh1wiP6liAeCHXCE/qWIF4IdcIT+pYgXgh1wgP6l\niCeCHUgw+pci3gh2IMHoX4p4I9iBBKJ/KRKBYAcShP6lSBSCHUgQ+pciUQh2IAHoX4pE4m8WEGf0\nL0WiEexAnNG/FIlGsANxRP9SuECwA3FC/1K4QrADcUL/UrhCsANxQP9SuESwA3FA/1K4RLADAaN/\nKVwj2IEA0b8UYcC+K0SyJot3GEh27NjhP//7779t+HqKFi1qT4IGhf6lCINCMe3JAiLiq6++snVY\nvPXtgylcuLCttjhgwIBA+pf+85//9PuXrly5klZ3cIKlGETKZ599lutQ92b3b7/9doG/L/1LESYE\nOyKlQ4cOts55Xlx99dUF/r70L0WYsBSDyBk4cKAZMmRIrj63atWq5ttvv83VHvM5c+aYhQsX2rZ2\nlStXztS/tHnz5nb2r/6lq1atotUdnCLYETlbtmwxNWrUyHSTNKca6b179z7o523fvt2GuZZcDjvs\nMPPEE0+Ya665xuzZs8c0a9bMb3X34osvmquuuiqQnwPIL4IdKTtrz8tsfenSpaZJkyaZ3te2bVvT\noEEDM2LECL9/6YwZM2h1B+cIdqTsrD23s3V58803TceOHbP9eOnSpe0SDK3uEAbcPEUkaa371ltv\nzXG2rv6jufXzzz/n+PEjjzwy0955wCWCHZGu15LdDhlVWsxLUa6DBbuWdBo1amSGDx9u0tPT8zxW\nIEgEO1Ju1p7X2Xpugt071aqTpyr4paUgwBWCHZGftWv9uyCz9dwGe8btj9odA7hCsCPys/YuXbr4\nrxXyeZ2t5zXYq1WrxgElOEWwI/Luu+8+WxNGbrzxxjzP1vft22fWrVt30M/T19U2S+2O0Y4cwBW2\nOyKS9Ndas+zly5fbLY96/ssvv9jDRCVLlrTBe+KJJ5pixYrlqiG1inrl5MILLzSPPfaYrb8OuEbZ\nXkTCX3/9Zd577z175H/JkiVm8eLFZvPmzTn+GYW6DhidfPLJ9qEDR8cee+wBn6dfCNmpU6eOGTly\npGnTpk0gPwcQBIIdSe37778348aNsz1F87oTZffu3faXgB6iE6Pt27e3h5YU8occ8v8rlatXrz7g\nz2ob5QMPPGBuvvnmXM36gUQi2JGUZs+ebeu1zJo1yy67ZKT19Fq1apm6deuaE044wdZ20fq3glph\nroNE3333nQ1s1W/ftm2b/XP6OioJoIeWVHr27Gl69eplPv3000xf/7rrrjPDhg3LVAgMCBPW2JFU\ntLxyyy23mMmTJ2d6v8JbR/7PPPNMG+q5vUGqv/5aQ9esffr06WbZsmWZPl69enVz77332obUWpuf\nMmWKOeOMMwL9mYCgEexIGgpezaIzNtLQac/LLrvMtqALYklEM/hp06bZKwGt23u0PPPwww8fsCce\nCCOCHaGnkrlqEJ1xlq694ppJn3TSSXH5nlqe+fe//22XZTLO3idNmmRatmwZl+8JBIVgR6ht2rTJ\ntGvXzp7m9G5wXnnllXYGHWQT6ux8/PHHtvyvxiFa4tFyTE6VHgHXCHaE1vr16+2auZZHvBovajzd\nuHHjhF8xqOb6zJkz/ZuzmrnrFwwQRgQ7Qmnr1q2mVatWZsWKFfZ17dq1zVNPPWUOP/xwJ+PRP5Nn\nn33WPP300364v/baa5QOQCgR7AgdbUnUzpMFCxbY1/Xr17ehntcm1fGg4l7aZukty3zwwQemRYsW\nrocFZEKtGISO1rS9UD/uuOPsyc4whLpcffXV9kauV6ZXBcb+/PNP18MCMiHYESrqLTp06FD7XFsL\nFerlypUzYdKtWzd7Q9drsKEywECYEOwI1RKMTnXu2bPHvu7fv79tORc22plz++2325LAol8+H330\nkethAT6CHaGhmbp38lN7xVW3Jax0FaFuSZ7rr7+eJRmEBsGOUPjtt9/MI4884i/BaHlDM+Mw0w3e\nc8891z5X7ZmJEye6HhJgEewIhbS0NLNr1y5/DfuII44wydJ6zytlMGbMmAMKkgEuEOxwTh2Kxo4d\n628hTKa94dpXf/bZZ9vnOkg1d+5c10MCCHa49/7775uvv/7aPlcd9LJly5pkoiJkHs3aAdcIdjjn\nzdb3D8mC+Ne//mUPDnk7bOJJLfbUScmrQPnrr7/G/XsCOSHY4dTevXttSztROKo5RhDWrl1rm2UU\nKRL/XjK6yXvRRRf5P8+HH34Y9+8J5IRgh1MKYG+bYJMmTQL5mmpereYZxx9/vEmUjOWDvVZ7gCsE\nO5zKGIJqY1dQqtGuipDy1ltvmaZNm9qHjv7Hk2q1e2WECXa4Rs9TOLV48WL/eRDLMArx33//3daa\nOf/88/2+pKo5E0+q9qilpC+//NLWjtdOH68ZNpBoBDuc8ma36id67LHHFvjrqQHGmjVrbLDfdNNN\npmLFiiZR9ItJwb5jxw67xBTEFQiQH0wp4NSPP/7oL2Vo1huEb775xja3TmSoS40aNQ74uQAXCHY4\n5TWM1ow9KNoTn5cbp6pJoxZ4BZXxZ/BO0QIuEOwIRbAXLVo0kK+n3TDaFZPbYFfT6g0bNgSybOKV\nFsj4cwEuEOxwyttnrpuNQfBOsOY22PX5Kr8bxLJNxp8hEfvngezwtw9OeVsE1Y0oqPX1nIJd5QtG\njx5tNm/ebJo1a2ZveHqnRgsq488Q5NISkFfM2OGUVxdmy5YtgXw9La3IoYceesDHZs6caUaMGGHu\nu+8+M2/ePLvfXQ2qg9q9kvFnSLZ6N4gWgh1ONWzY0L5dt26d2b59e4G/Xu3ate3bgQMH2ho0Tz/9\ntNm4caOdTT/++OPmjjvuMI0bN7ZlADp06ODvPw+Cqjt6GjRoEMjXBPKDYIdTOlCUVTDmlxpfXHPN\nNTbMn3vuOTN+/Hi73KNeqgr3Vq1a+Z+rm6x6X1Az9tWrV9u3tWrVMuXLlw/kawL5wRo7nDr55JMz\nBWPz5s0L9PV02rNPnz72kdHWrVvt3vaMp0Hffvttu2RStWpVE8QSkK469v+ZABeYscOpjIW/vBlv\nPKikgLZCar+6Svmq+fS4ceP8pZuCyni1QbDDNWbscEqzaIWrjuCrDIAqPWZ147OgtI7es2dPc//9\n99utiFoDr1+/fmAVIOfMmeM/L+hVB1BQhWI0aYRjQ4YMsTc75a677jKXXnqpSSaqDdOuXTt7KEml\nEbTlMqjyCEB+sBQD57p27eqfPH311VeTriH0jBkz/JOmuiog1OEawQ7nVFr3kksusc8121WFxGSh\nX0L6ZeSVFLjhhhtcDwkg2BEOvXv39p9PnDgxaWbtc+fONT/88IN93qlTJ1OpUiXXQwIIdoRDy5Yt\n/d0kn3zyid2KGHba4jh8+HD/9a233up0PICHYEco6CToM888469PP/bYY2bTpk0mzFSewCsjcMst\nt2Q6bAW4RLAjVHva77nnHv9U6NChQ0O7JKNaM7NmzbLPa9asaYYNG+Z6SICPYEeoaNujV2dl/vz5\nZvLkySZsdMJUv3Qy3hOIx957IL8IdoSKdpakpaX59cy1JKPthGGhcr+60ZtxCSZj/RkgDAh2hHJJ\nRuV0PYMGDQpFuK9fv9706NHDrwnTunVr8+ijj7oeFnAATp4itEaOHJmpmFe/fv1M586d7Y3WRPvu\nu+/srheFu5xyyim2jEDp0qUTPhbgYAh2hNqoUaMybSNs0aKFvcGqQ02JsHfvXvPyyy/b2u67d+/2\nt2bqCqJcuXIJGQOQVwQ7Qu+FF14w3bp1M+np6fa1blRq9t6xY8e4zt518OjBBx80y5cv99/Xvn17\nM3XqVFOqVKm4fV+goAh2JIVly5aZ6667zjbMyFhF8dprr7Vvgwx4LbdMmzbNvPTSS/4sXUH+8MMP\n2xunGWu6A2FEsCNpaMauk56DBw/2Z+9SrVo1WxHy/PPPN2XKlMnX1963b59ZuHCheeWVV+w2S732\naNfLhAkTbE13IBkQ7EjK2bt2p3z66aeZ3l+8eHHTrFkzU7duXf+RXe0WVWNUDXg191CTjM8//9zf\n7eI5/PDD7VIMs3QkG4IdSWvRokVmzJgx9hCTVzZ3fxUqVLDNPBT6CmctrezatcuGuG6MZkVLOwrz\nyy+/3JQsWTLOPwUQPIIdSU+HhXSoadKkSWbFihXZBnZOtMvmvPPOM7169aLmC5IewY5I2blzp63n\nvmTJEvvQjhbVndEsXYFfokQJOwuvUaOGrSbpPapUqeJ66EBgCHYAiBjuCAFAxBDsABAxBDsARAzB\nDgARQ7ADQMQQ7AAQMQQ7AEQMwQ4AEUOwA0DEEOwAEDEEOwBEDMEOABFDsANAxBDsABAxBDsARAzB\nDgARQ7ADQMQQ7AAQMQQ7AEQMwQ4AEUOwA0DEEOwAEDEEOwBEDMEOABFDsANAxBDsABAxBDsARAzB\nDgARQ7ADQMQQ7AAQMQQ7AEQMwQ4AEUOwA0DEEOwAEDEEOwBEDMEOACZa/g8sYICwXspK/AAAAABJ\nRU5ErkJggg==\n"
        }
      },
      "cell_type": "markdown",
      "id": "a2575ec1-8882-49c8-b4ac-1c3bcadefdc4",
      "metadata": {},
      "source": [
        "# Multilevel models\n",
        "\n",
        "## Departmental tardiness\n",
        "\n",
        "Imagine that your university hires a new president. On her first week,\n",
        "she schedules a meeting with professors from three departments,\n",
        "${ D = \\{\\text{G}, \\text{E}, \\text{M} \\} }$. The professor from the\n",
        "Department of Government (${ d{=}\\text{G} }$) shows up 1 minute late,\n",
        "the professor from the Department of English (${ d{=}\\text{E} }$) shows\n",
        "up 10 minutes before the scheduled meeting time, and the professor from\n",
        "the Department of Math (${ d{=}\\text{M} }$) shows up 30 minutes late.\n",
        "From this single meeting the new president updates her beliefs about\n",
        "when professors will show up to scheduled meetings.\n",
        "\n",
        "But how should the president update her beliefs? It could be that\n",
        "departments have different relationships with punctuality and it would\n",
        "be useful to expect professors from the Math department to show up later\n",
        "than professors from the English department. It could also be the case\n",
        "that there is not meaningful between-department variance and she should\n",
        "just learn about the greater population of professors at the university.\n",
        "\n",
        "A model that infers information about the population without\n",
        "differentiating subpopulations is doing *“complete variance pooling”*,\n",
        "meaning that all of the variance in the data is treated as being\n",
        "generated by the same distribution. I.e. the data is pooled together and\n",
        "treated as iid observations drawn from a single distribution.\n",
        "\n",
        "If the president thinks of each department independently, what she\n",
        "learns about the tardiness of professors from one department will not\n",
        "change her expectation about professors from other departments. In this\n",
        "case, there would be *“no variance pooling”* between departments. The\n",
        "president’s observations of the three professors have no bearing on her\n",
        "belief about what to expect of a professor from a fourth department\n",
        "(e.g. Psychology).\n",
        "\n",
        "*“**Partial** variance pooling”* models represent uncertainty at\n",
        "multiple levels of abstraction. Based on her observation from this\n",
        "meeting, the president will update her beliefs about these departments\n",
        "and the greater population. In other words, what she observes about the\n",
        "Math department might strongly update her beliefs about that department,\n",
        "moderately update her beliefs about the broader population, and that in\n",
        "turn could weakly update her beliefs about the English department.\n",
        "\n",
        "## Complete pooling\n",
        "\n",
        "Let’s build a complete pooling model of the president’s belief about the\n",
        "tardiness of professors. The data observed by the president is\n",
        "${ t = \\langle t_G, t_E, t_M \\rangle = \\langle 1, -10, 30 \\rangle }$.\n",
        "We’ll model the president’s belief about the greater population of\n",
        "professors as a normal distribution centered at $\\theta$ with standard\n",
        "deviation $\\dot\\sigma$. Prior to meeting with anyone, the president\n",
        "thinks that professors will, on average, show up on time, but that some\n",
        "will show up early and some late. We can encode this prior belief as a\n",
        "normal distribution centered at zero, with some standard deviation,\n",
        "let’s say 20. I.e. the prior for $\\theta$ is:\n",
        "\n",
        "$$\n",
        "\\theta ~\\sim~ \\mathcal{N}\\left( \\dot{\\mu}{=}0, \\dot{\\tau}{=}20 \\right)\n",
        "$$\n",
        "\n",
        "If she observers that professors tend to be late, then her posterior\n",
        "estimate of $\\theta$ will have an expected value greater than zero\n",
        "(i.e. ${ \\operatorname{E}[\\theta \\mid t] > 0 }$), and if she ends up\n",
        "expecting professors to be early, then\n",
        "${ \\operatorname{E}[\\theta \\mid t] < 0 }$.\n",
        "\n",
        "In this tutorial I’ll use the dot above to indicate that a variable is\n",
        "fixed, not inferred (the dot is not standard notation, just what I’m\n",
        "using here). In the graphical schematic, fixed parameters are depicted\n",
        "as bare symbols rather than nodes. Note that “fixed” means the values\n",
        "are not updated during inference (but you can, of course, manually tune\n",
        "these values or learn suitable values by minimizing some loss function\n",
        "in conjunction with a model fitting algorithm).\n",
        "\n",
        "Since the model infers the distribution of $\\theta$, the graphical\n",
        "schematic depicts it as an open node (i.e. a latent random variable).\n",
        "The data in this model are $t$.\n",
        "\n",
        "We specify the likelihood as,\n",
        "\n",
        "$$\n",
        "t_d ~\\sim~ \\mathcal{N}(\\theta, \\dot\\sigma{=}15)\n",
        "$$\n",
        "\n",
        "I.e., for whatever the true value of $\\theta$ is for the population of\n",
        "professors, the president thinks that when professors actually show up\n",
        "will follow a normal distribution, centered on the true value of\n",
        "$\\theta$, with standard deviation 15.\n",
        "\n",
        "$$\n",
        "\\begin{align*}\n",
        "\\theta ~\\sim&~ \\mathcal{N}\\left( \\dot{\\mu}{=}0, \\dot{\\tau}{=}20 \\right) \\\\\n",
        "t_d ~\\sim&~ \\mathcal{N}(\\theta, \\dot\\sigma{=}15) \\\\\n",
        "d ~\\in&~ D,~~~ \\text{where}~~~ D = \\{\\text{G}, \\text{E}, \\text{M} \\}\n",
        "\\end{align*}\n",
        "$$\n",
        "\n",
        "![](attachment:generated/multilevel-models/complete-pooling.png)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "3fd04022",
      "metadata": {},
      "outputs": [],
      "source": [
        "%reset -f\n",
        "import jax\n",
        "import jax.numpy as jnp\n",
        "from memo import memo\n",
        "from enum import IntEnum\n",
        "from jax.scipy.stats.norm import pdf as normpdf\n",
        "from jax.scipy.stats.cauchy import pdf as cauchypdf\n",
        "from matplotlib import pyplot as plt\n",
        "\n",
        "normpdfjit = jax.jit(normpdf)\n",
        "\n",
        "class Department(IntEnum):\n",
        "    GOVERNMENT = 0\n",
        "    ENGLISH = 1\n",
        "    MATH = 2\n",
        "\n",
        "t = jnp.array([1, -10, 30])\n",
        "sigma = 15\n",
        "\n",
        "Theta = jnp.linspace(-40, 40, 200)\n",
        "\n",
        "@jax.jit\n",
        "def professor_arrival_likelihood(d, theta):\n",
        "    ### likelihood of a professor from department d \n",
        "    ### showing up t_d minutes early/late, \n",
        "    ### for a given theta and sigma.\n",
        "    return normpdf(t[d], theta, sigma)\n",
        "\n",
        "@memo\n",
        "def complete_pooling[\n",
        "    _theta: Theta,\n",
        "](mu=0, tau=1):\n",
        "    president: knows(_theta)\n",
        "    president: thinks[\n",
        "        department: chooses(d in Department, wpp=1),\n",
        "        department: chooses(theta in Theta, wpp=normpdfjit(theta, mu, tau))\n",
        "    ]\n",
        "    president: observes_event(wpp=professor_arrival_likelihood(department.d, department.theta))\n",
        "    return president[Pr[department.theta == _theta]]\n",
        "\n",
        "mu_ = 0\n",
        "tau_ = 20\n",
        "res = complete_pooling(mu=mu_, tau=tau_)\n",
        "\n",
        "### check the size and sum of the output\n",
        "# res.shape\n",
        "# res.sum()\n",
        "\n",
        "fig, ax = plt.subplots()\n",
        "\n",
        "ax.axvline(0, color=\"black\", linestyle=\"-\")\n",
        "theta_posterior = res\n",
        "theta_expectation = jnp.dot(Theta, theta_posterior)\n",
        "ax.plot(Theta, theta_posterior, label=r\"$P(\\theta \\mid t)$\")\n",
        "ax.axvline(\n",
        "    theta_expectation, \n",
        "    color='red', \n",
        "    linestyle='--', \n",
        "    label=(r\"$\\operatorname{E}\" + rf\"[\\theta \\mid t]={theta_expectation:0.2f}$\")\n",
        ")\n",
        "_ = ax.set_title(r\"Posterior of $\\theta$\")\n",
        "_ = ax.legend(bbox_to_anchor=(0.9, 0.5), loc='center left')\n",
        "\n",
        "_ = plt.suptitle(fr\"$\\dot\\tau$ = {tau_}\", y=1)\n",
        "plt.tight_layout()\n",
        "plt.show()"
      ]
    },
    {
      "attachments": {
        "generated/multilevel-models/no-pooling.png": {
          "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAF2CAYAAAB6XrNlAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90\nbGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAO\nxAAADsQBlSsOGwAALvFJREFUeJzt3Qd0VMX7PvCXKjX0IigQivSiFBWpIiBNUFAEEZBepKoIiAWQ\nIkWlCUiVoqIgRRBpioDSQYqUKEhVeu8g+z/PfP/3/jYkQMq9O7c8n3P2ZBNCdkLIs7NzZ943QSAQ\nCAgREXlGQt0DICIiazHYiYg8hsFOROQxDHYiIo9hsBMReQyDnYjIYxjsREQew2AnIvIYBjsRkccw\n2ImIPIbBTkTkMQx2IiKPYbATEXkMg52IyGMY7EREHsNgJyLyGAY7EZHHMNiJiDyGwU5E5DEMdiIi\nj2GwExF5DIOdiMhjGOxERB7DYCci8hgGOxGRxzDYiYg8hsFOROQxDHYiIo9hsGvUvn17SZAggZw6\ndSrKn1WqVEkyZsyoZVxE5G4Mdo22bdsm2bNnjzbAt2/fLsWLF9cyLiJyNwa7JoFAQHbs2BFteB88\neFDOnj3LYCdH27Rpk3rFeb/b77//rnuovpNY9wD8at++fXLp0iUpUaJEtDN5KFasmIaREcVMkiRJ\n5P333zffnzt3rnql+fbbb0uyZMnUxxDsRYsW1ThKf2Kwa2KEd3Sz8nv9GZFT4P9n8P/RZcuWSdq0\naWXw4MFax0VcitHmXuGNl66JEyeWQoUKaRgZUdz88ccfUrhwYd3DIAa73mBPnjy55MuXL8qfbdiw\nQfLnzy8PPPCAlrERxdahQ4fk/PnzUqRIEd1DIQa73mDPmjWrJEyYMMoFqSNHjnAZhlwFGwGAwe4M\nDHYNzp07p3a+HD16VN0MmPF06NBB3ecyDLnJzp071VsGuzPw4qkG2DkAadKkkdKlS0u9evXkypUr\n8uOPP0qKFCnMHQa5c+eWRo0aaR4t0f1FRESotwULFtQ9FOKMXe+F0wkTJkiZMmVk6tSpsnTpUmnZ\nsqV6i0NLJ06cUEs1RG6AV5uQOnVq3UMhztj1BnvFihWlbt26Uf4ca+xEbpI3b171Fq8+n3jiCXnz\nzTclLCxM97B8K0EARyAppLD8cvz4cbWTgMgLTp48KU2bNpWVK1eqDQEXL16MsjGAQofBHmL//fef\npEqVSqpUqSILFy7UPRwi8iA+pWq4yHTt2jWWCyAi23DGTkTkMZyxExF5DIOdiMhjGOxERB7DYCci\n8hgGOxGRxzDYiYg8hsFOROQxDHYiIo9hsBMReQyDnYjIYxjsREQew2AnIvIYBrsH+qd+8803uodB\nHrF161bZsGGD7mFQPDHYXWzjxo2qeXDDhg1lyZIluodDLnf16lX1f8nogHT79m3dQ6I4YrC72IUL\nF+To0aPqfps2bVTXGqK46tu3r/z555+CSt5oz8gOSO7Feuwu17p1a5k4caK6//rrr8uoUaN0D4lc\naPPmzfL444+rDl8ZMmSQXbt2SebMmXUPi+KIwe6BNfbChQvLP//8IwkSJJBVq1ZJuXLldA+LXOTm\nzZuqD6/RZH3GjBnyyiuv6B4WxQNfa7lc2rRpZezYseo+nqNbtWqlWu8RxdSQIUPMUK9Vq5Y0btxY\n95Aonjhj94iXX35ZZs2ape736tVLBg4cqHtI5AK7d++WEiVKyI0bNyR16tTyxx9/yMMPP6x7WBRP\nDHaPOHHihBQqVEhOnz4tiRIlUlvWHnvsMd3DIgfDenr58uVl7dq16n288mvXrp3uYZEFuBTjEbjQ\nNWLECPMXtmXLlmrtlOhuxowZY4Z6xYoV1c4q8gbO2D0EP8ratWvLDz/8oN4fMGCA9O7dW/ewyIEO\nHDigLrpfuXJFkiVLJtu3b5d8+fLpHhZZhMHuMYcPH1a/sNjTnjRpUnVRrECBArqHRQ6CX/lq1arJ\n8uXL1ftDhw5VB5LIO7gU4zG48IVdDoALYliS4QlCCjZ16lQz1EuVKiVdu3bVPSSyGGfsHoQgr1y5\nstrTDiNHjpROnTrpHhY5wL///qsusuP8Q+LEiWXLli1StGhR3cMii3HG7kE4Co7TqFg7NbY/Yk2V\n/A1zuA4dOqhQB1x/Yah7E4Pdo3AhrF+/fur+5cuX1Y4Hvjjztzlz5si8efPUfczaeWHdu7gU42G3\nbt2SJ598UjZt2qTenzJlijRv3lz3sEgDnG9AmOO8A0pPYJsjasOQN3HG7mFYQ500aZJ6C926dZNj\nx47pHhZp0L17dxXqgIulDHVvY7B7XLFixdQaO2BttWPHjrqHRCH2448/yrRp09T93LlzS//+/XUP\niWzGpRgfuH79uiovgFKsMHv2bKlfv77uYVEI4DwDmrEcOnRIvb9ixQp5+umndQ+LbMYZuw888MAD\nakkGa6uAWfuZM2d0D4tCAK/WjFBH5U+Guj8w2H0C7c66dOmi7h8/flytuZK3rV69WtWDgWzZsqkT\npuQPXIrxEWx7xL7lv//+21x7rV69uu5hkU39S4sXL65a3cH8+fPlueee0z0sChHO2H0kZcqUMmHC\nBPN99kn1LpxhMEIdtfoZ6v7CGbsPsU+qt6FMQJkyZdi/1McY7D6EbY84rIK6IcZaLPukegP7lxJw\nKcbnfVIBFSDZJ9Ub2L+UgDN2Hwvuk9qzZ08ZNGiQ7iFRPLB/KRkY7D7GPqne7V86btw4adu2re5h\nkSZcivEx9kn1bv9SXCAn/+KM3efYJ9X92L+U7sRgJ/ZJdTH2L6XocCmGou2TiqUZcj72L6XocMZO\nCvukug/7l9LdMNjJhCPoqN+OPe0oP7Bz507JlSuXGfwRERGqG9OpU6fk0qVLnNWHQIoUKSQsLEwK\nFiwoJUuWVD8XwK8tSi/PnTtXvf/ee+9J3759NY+WnILBTpFgjbZHjx7qPmbwzz77rPz000+yfv16\nNTNMkiSJZMyYUQUM7pN98GSKYl74d79w4YLakoonXrQ7xBOu8XPC9ZHNmzer8sxEwGCnKH1ScST9\n999/V+9jl8Xzzz+vyv4iUFAxEBdYKXTwK4qa6uvWrVM3VOXcs2eP+ef4GFvdUbD/NcMk+v9+/fVX\ntUvGgBAfPny4PPjgg1rH5WdokJIzZ051a9iwoTpQFhzsOD2MJ+OECbkXgv6H/xPItGDBAlWfHScY\njUYcWAJABUhyBszWp0+fru6Hh4erg0mjR4+WJk2aqB1NRAqWYoimTp0aSJQoUaBVq1aBW7duBa5d\nuxYoVKgQlunUbfbs2bqH6HsXLlwI5MiRw/yZrFixQn186dKlgZQpUwZq1KgRuHz5su5hkgMw2Cmw\ndu3aQIIECQI9evQI3L59O8rHESJZsmQJnD59Wus4/a5jx45mqLdu3TrSn61fvz6QLl26KB8nf2Kw\n+5wxM69SpUqkUDd07drVDJNmzZppGSMFAqtWrTJ/DtmyZQucPXs2yud89dVXkWby5F/cFeNz77//\nvgwbNkx27NghuXPnvm+f1MWLF6stkOS8/qX4Va5Xr546f4B6Mcaed/Ih3c8spE9EREQgceLEgU8+\n+eSen7d8+XJztog1Xqz1Uuj07NnT/Pd/+eWX7/m5R44cCYSFhQV69eoVsvGR8zDYfeyNN95QQY2L\npfeDi6pGuGCtl0Jj8+bN6qI2/t0zZMgQOHHixH3/zgcffKA+F8ts5E/c7uhT2Br3xRdfSIsWLdSJ\nxpicSDX2smOLHfqkkr1QGx8/H6N0A+r3ZMqU6b5/D3/nzJkzavsq+ROD3afwS4+DLq+99lqc+qS2\natWKfVJD3L+0UaNGMa7WiesgEydOtHmE5FS8eOpTNWrUMC+Gxgb7pLqjf+l3330nDRo0kP3795uF\n3Mg/OGP3Icy0V6xYoUI6trAckCFDBnN5BqViyVpGm0LjJCn+nWPblBpdsZInTy5Lly61aZTkZAx2\nH0KBL6zfoqhXbLFPqjv6l6LGDxpvoCon+Q+D3Yfwy54uXbo498Vs3Lix1KxZ03ySMLovkTX9S3v1\n6mVW1sQ6eVyLe6HiI4PdnxjsPmSUeUXVwLjA3xs3bpxa+4V+/fpFqjZIcYPLXW3atFFNqaF///6S\nN2/eOH89lFretWuXKuRG/sJg9yHM4uJbv5t9Uu3pX7ps2TJ1H2V449u/FD9jPFls3LjRohGSWzDY\nfQbH01EeADsu4guzywoVKqj7v/32m3z22WcWjNC//UuNUsnoXzpp0iT1Nj6yZ8+uul2hxAD5C4Pd\nZ44dO6beZsuWLd5fC2u/WAPGWjBgbRhrxBQ7mFV37NhRtcCD3r17W9aUGj/n48ePW/K1yD0Y7D4N\n9ixZsljy9XABFmvsRsEwzOJ5NCJ25syZYzalRv9SBLtVsmbNav7MyT8Y7D5jdbBDt27dpGTJkuo+\n1ohRqoBiBkf/MVs3XgFhCcbKptQMdn9isPsMXpajPICxfGIFrAVPnjzZXBNG0GPNmO4P/1YnTpxQ\n93Gx1Oqm1HgC51KM/zDYfQazNytn64ZixYqZ+6+xVsw+qTHrXzpt2jR1H7XwjSUtK3HG7k8Mdp85\ne/aspE+f3pav/c4770ihQoXMWiVYO6boXbx4Udq2bWu+P2HCBFsaY+BnjeUe8hcGu89gzzmOm9sB\na8NYIzYOPmHtmKESPby6OXTokLqPkgFPP/20LY+DnzVLPvgPg91n8EtuV7Abpx27dOmi7mNt19ib\nTf8HtexRD8bYjmhnSYYkSZKog2O3b9+27THIeRjsPnPr1q14H3y5nw8//FDCw8PVfeyQwVoy/d8B\nMdSyN6DGPS5m2xnswFm7vzDYfSiuNWJiCmvFWDM2YC0Za8r0v7o6ERER6j7KJkfXlNpNP2tyJgY7\n2aJKlSqqfgxgLdnKQzduhdr1qK0OqGmP2vZEdmCwk22GDRtm9kkdPXq0rFmzRvwKSyHBhdJi2r+U\nKC4Y7GSbO/ukItj82icVM3XUro9t/1KiuGCwk63q1q0rDRs2VPextty3b1/xY/9S4/tGDXs82XHt\nm+zEYCfbYdnBOBTltz6pVvQvJYotBjvZzs99UlGjPr79S4lii8FOIfHKK69E6pNq7A7xMiv7lxLF\nBv+XUUjc2ScVa85e7pNq9C9FjXor+pcSxQaDnULGT31Sg/uXlipVKt79S4lig8FOIeWHPql39i8N\nrlVPFAoMdgopP/RJRS16O/qXEsUUg51Czst9UmfPnq1q0dvRv5QophjspIUX+6QG9y/FxWKr+5cS\nxRSDnbTAmjOCz0t9Uu3uX0oUUwx20qZ48eKe6ZN6Z/9SbG8k0oXBTlrFpE8qCocZM2GdsEUzusbQ\noepfShRTDHbS6l59UnFBdcqUKar0b/bs2WXVqlXaxon99mXKlFFjefXVV+X06dPR9i9FdyS7+pcS\nxRSDnbSLrk8qtkBWr15dWrRooZZp0NLv+++/1zbGffv2ybZt29T9GTNmSMGCBeWbb75RTzbB/Uv9\nUCqBnI/BTo5wZ59UBKdxctNw+PBhTaOL+tgnT55U5YiN+jeh6F9KFFMMdnIErEn36dPHfD+6hhxO\nCnaDUQumdOnSUqdOnRCPiih6DHbSDiV8Bw4cKO3bt7/n5zkx2A0bN26UqlWryv79+0M2JqK7YbCT\nVli3xkVJ7I4xmlHczT///KOtaFhMnlRWrFihygeg9rxXTtKSOzHYSZvr16/LM888Y/YCvR+Euq5D\nTEeOHInR5125ckUdTvr2229tHxPR3TDYSWtBsBQpUsTq7+hajont43IfO+nEYCdtkiRJImvWrJEX\nX3zRM8GeK1cumTdvntSqVcv2MRHdDYOdtDffwH7wn376SYoUKWJbsGPN+/bt23H6uzhZev78+Xt+\nTvLkyVUZgV27dkndunXj9DhEVmH1f3KEypUry9atW1X7vHfffdesZx6bYEd442DT5s2bZdOmTert\njh075MKFC2r7JP4crxIQwqjnguqSxq1YsWJmjfjYPCZgPzsOJuFJisgRAuQrzZo1C9SsWTPgZCdO\nnAi0adMmkCBBAmwtiXSrWrVqlM/fs2dPoEuXLoGMGTNG+fyY3hInThyoVq1aYO7cuYGbN29G+vqL\nFy+O9u8ULVo0sHLlyoCTzZs3T4312rVruodCIcQZOzlOpkyZZPz48aqwVqdOnVQLPcOWLVvUW5QY\nWLBggWqth22G0UmaNKnafpg5c2Y1S0+UKJGauWNp5Y8//lCnRw34ekuXLlW3hx56SD026r5kzZo1\nytdPly6dWnbB57DlHTkR/1eSYz322GPq4ipqs7Ru3Vptj0yfPr0qNYBQ/fvvvyN9Pmq1PPfcc+by\nCjoYIdyjg2UZbGHEcg1uv/zyi6xevVr9GT6O5SB0eULdmpw5c5p/Dxd68WSSMWNGm797orhLgGl7\nPP4+uUzz5s3VTHXRokXiJrh4OXXqVLXnHW+DValSRTp06KCO9GMNPa5w4RP1XlCrBrN6A+rWtGzZ\nUp588kkpW7asuMn8+fOlXr166pUKuzn5B3fFkCts2LBBPv7440ihjiep3bt3y/Lly+WFF16IV6gD\n6sKPGjVKnXDFrDxDhgzq43iMHj16qKWf6GrYEDkNg50cb/DgwVKtWjWz5jnWwNGxCLXaCxQoYPnj\npUqVStWtwTo8njAAWyU/+ugjqVChgpw6dcryxySyEoOdHAurhD179jTb5wEuaO7cuVPVardblixZ\nZPbs2fL111+bs3cU+6pYsaIcPXrU9scniisGOzlW79691SwZsMwya9Ys1XYuTZo0IRsDOjthnzr2\nw6NHq7EWj3V9J7TrI4oOg50cCYGOJRjAVkV0T3rppZe0jQct8VauXKkuoMLevXvVq4b7nUgl0oHB\nTo7z888/qyUYY6Y+d+7ckCy93A+6I2FtH9swATt00KOVyGkY7OQoly5dUn1ODdOmTXNEqBvCwsJU\nuKMkAcycOVNtKSRyEgY7OQpm6qj3Ak2bNpWXX35ZnHgyFtsusf4OOCx1+vRp3cMiMjHYyVFLMGPG\njDHXtD/99FNxqvLly6tyB3D8+HHp0qWL7iERmRjs5AjYJ96uXTvz/c8//1zVZHEy9GnNkyePuSSD\nJyYiJ2CwkyOg+FZERIS6j+WX2rVri9OhSxK2XxpwapXICRjs5Ag4wm94++23xU115EuVKqXu4yJq\nTHujEtmJwU7a4WLpwoUL1X3sEy9RooS4CQqQGctJWEIi0o3BTtohDI0io0ZIWm3Pnj3qgBNKA2DL\nIpZ6Dh48aMnXxslU43oAlmZu3LhhydcliisGO2n35Zdfqreocd6gQQPLvz4aZeBQ0fr161VFSDSa\n/uGHH1RhMStCOEWKFPLaa6+p+8eOHVMnVIl0YrCTVqi3YsyccRDpbn1H4wr9SuvXr6+Wd1A8bPjw\n4fLVV19Jt27d1MVaqw4XocGHAYXCiHRisJNW6F5kQNcjq7311lty5coV1YUpderUkZZPILjtXnw8\n+uij5oGl4O+JSAcGO2kVHILG7hIrL8p+++23qqa6UQLAgBZ7YFWFRqzbP/LII44N9oQJ7ftV37dv\nnyoDkSNHDlXbB09wd962bt1q2+NTVOx5SloZIYhffsx6rYQyv9ip0rhx4yh/hv6pYGW7OLziQNVH\nNARBMw4n9UU1Xk1YDc3FseXz8uXL8vzzz0u+fPlUtyujAXiuXLlU2KM7FYUOg520QpciQCCgc5HV\nh54ATaoRQMGM/eZogG0VXKA1LgRjPb9SpUri5WDHkyN2Gl29elUFORqQGPAqCVU5R48erS5WU2gx\n2Ekro2k0CmtZCdsnMXOEYcOG3fXzrGytFzxDR5VKJ7Ej2KdPn66WYfr06RMp1IODHU+oDPbQ4xo7\naWU0h7Z6Nwxm5AhXVIhEyN95M5ZnjMYZVkBDEANmsV6H6xd4woju7IHx6uu///7TMDJisJNWWAO3\nY0b577//mlUio3tMLB3ggqpRxMvqC5R+CLR169ZJ4cKFo/03NnrChoeHaxgZMdhJK2OmblzMtMrN\nmzfV26RJk0b5M4Q6Su1Gd1E1PoJn6cGzdy86c+aMXLhwQbJmzRrtn6MZCTjpOoOfMNhJKyMAjbV2\nqxhr9tFtZ0S5XZwWNeqpWyX4e/B6sBuvsBDud9q9e7e6cI3dMjlz5tQwOmKwk1bGUgjCwJhlWyFv\n3rySOXNmWbBgQaQLmQMGDFBH/vEWf26lHTt2RHp8JzFq8VgFtXEeeughtV31r7/+Mj+O5t64roHl\nLjQkJz0Y7KSVcdoUSzG7du2ydL27e/fuaq0dF0hRCrhq1apqB0ebNm1s6Xhk7MlH6Dltbdm4lmEl\n/JviWsJTTz0lXbt2lc6dO6tdRmjyjWJopUuXtvwxKWa43ZG0Ci4jgGAsXry4ZV+7R48ecuvWLRk3\nbpyMHDlSHZL54osv1IzSani1sX37dnM/u10HgpwyY4fXX39dPYFir/rYsWNVyQZse+zVq5flp4gp\ndhjs5Jhg37RpkzqabhWE6zvvvKNuoThoZVwAtqPmjRODHbDV0a5SyxR3XIoh7WvsadOmVffRM9Su\nALLbTz/9ZN53YrCTvzDYSSvMqo2TiWiGsWrVKnEbrF9juceoPVOlShXdQyKfY7CTdsEv5YN7n7pp\ntv7nn3+ajbjRpYlIJwY7aYddK8ZF0++++848NeoWwU9GHTt21DoWImCwk3bB9UaMXSxuge5PRhcm\n7AThFj9yAgY7OQKO96dJk0bdHzJkiGpb53S40NuuXTtzjzh3h5BTMNjJEVANsF+/fmbFRzSHdnoh\nralTp5o1UdBTtUmTJrqHRKQw2MkxcOClfPnyZi/SESNGiFOhLDBOW0LixIlVyKNTEJETMNjJMXCK\ncfLkyWYBLRwssrLMgFXwSqJVq1ZmAax3333X0hOzRPHFYCdHQfGswYMHm0syzz77rOoh6qR1dex8\nWbJkibkEgyP0RE7CYCdHLsnUrl1b3T98+LAq3mU0btAd6m+++aaMHz9evY+LvTNmzOASDDkOg50c\nuSQza9YsqVChgnofO2TKlSun+mvqXH7BDpiPP/5YvY/lokWLFqkOQkROw2AnR0IjjO+//17Kli2r\n3j9w4IA8/vjj8s0332jZq44loc8//9wcG/auo1wtkRMx2MmxwsLCVCeeatWqqfdPnz4tDRs2lBdf\nfDHazkh2LL1g2aVIkSKyfPlyc/ll2bJlannIi9xahI0iY7CTo6VMmVLN3N966y2zWfTs2bPVEsiX\nX35p2153LP8gvLH8YnRgwsnStWvXmq8ivGjQoEHStm3baFvekXsw2Mnx0JAap1F//fVXyZ8/v/rY\nqVOn5JVXXlG7aNCC7eTJk/F+HDxJoJUell3wOGh6bTw++qQi1AsWLChehfaEffv2VUtOqN+D8g7k\nTgx2co0nnnhCtm7dGmn2jrX3nj17qv6bOPmJImLYHhnTJQXMxlevXi0ffvih5M6dW+rWrWtuZTRm\n6Vu2bFFbGnEQyavwpNayZUu5ceOGeh+Nvr38/Xodf3LkKtiNgtk72tuNGTNGpk+fLpcvX1aBNHPm\nTHWDjBkzqoYXaFOHptX4e4kSJVJ74xHmO3fuVK349u7dG+2TAHbktG/fXho0aOCLgMO/JV6RANrb\noS8suVeCAK+W+Erz5s3VsgW26nnB+fPnVbijdC6WEuJbrwZPGAh0XDD1AuzeqVevnnpCQxOQ6OBV\nD65ZXLlyRZIlS6Z6t+bLly/kYyXreH8qQp6GXSo40ITToFim2bhxo+qditk4ZuVoMn03Dz/8sJrV\nGzfslUdDZj/BvA6zc4Q6oBAbQ939GOzkmZruWHbBDbs6AM2l0W4Pyy3YJml0OELxrvDwcLVE43co\nXobtm8b1hG7duukeElmAwU6ehaUHFOcKnoXnyJFDHXQiUZ2qunfvru7jOsKkSZN8cT3BD7grhsiH\njGJm586dU+/37t1bihUrpntYZBEGO5EPzZkzR+bOnavuFypUSAU7eQeDnchnzpw5YzbdxrUJLMHc\nbccMuRODnchncIHUqLWDC8k4+EXewmAn8hH0aJ02bZq6j5O2/fv31z0ksgGDncgnLl68aG4FhQkT\nJqgia+Q9DHYin0C9G6PNIHq2Pv3007qHRDZhsBP5wJo1a1Q9GMiWLZsMHTpU95DIRgx2Ih9AXXnD\n2LFjJW3atFrHQ/ZisBP5wF9//aXeorTCc889p3s4ZDMGO5GHBTcAz5Ahg4wcOVLreCg0GOxEHoXK\nlqNGjTLfHzFiBAuf+QSDncijcIEUtdahevXq0rhxY91DohBhsBN5uH+pYfTo0ap8APkDg53I4/1L\njaYi5B8MdiIP9y9FyzvyHwY7kYdgTR0nTAH9S40qjuQvDHYiD/cvxSlT8h8GO5FHsH8pGRjsRB7A\n/qUUjMFO5HLsX0p3YrATuRz7l9KdGOxELsb+pRQdBjuRi7F/KUWHwU7kUuxfSnfDYCdyIfYvpXth\nsBO5EPuX0r0w2IlcZvXq1exfSvfEYCdykatXr6oZuoH9Syk6DHYiF0H9l4iICHWf/UvpbhjsRC6x\nZcsWc9mF/UvpXhjsRC7pX9qiRQvVRAPYv5TuhcFO5AKYqW/btk3dr1mzJvuX0j0x2Ilc1L80derU\nMm7cOPYvpXtisBO5qH/pkCFD2L+U7ovBTuSS/qUVK1ZUHZKI7ofBTuSS/qUoG5AwIX9l6f74v4TI\nJf1L8+XLp3tY5BIMdiIHYv9Sig8GO5HDsH8pxReDnchB2L+UrMBgJ3IQ9i8lKzDYiRyC/UvJKgx2\nIodg/1KyCoOdyAHYv5SsxGAn0oz9S8lqDHYizdi/lKzGYCfSiP1LyQ4MdiJN2L+U7MJgJ9KE/UvJ\nLgx2Ig3Yv5TsxGAnCjH2LyW7MdiJQoz9S8luDHaiEGL/UgoFBjtRiLB/KYUKg50oRNi/lEKFwU4U\nAuxfSqHE/1lENmP/Ugo1BjuRzdi/lEKNwU5kI/YvJR0Y7EQ2Yf9S0oXBTmQT9i8lXRjsRDZg/1LS\nicFOZAP2LyWdGOxEFmP/UtKNwU5kIfYvJSfgvivyZE0W4zAQXLp0ybx//fp1Fb6GJEmSqJOgVmH/\nUnICztjJU/bs2aN6h4aFhZm34sWLm3+O2ufBf5YqVSoZMGCAJY/N/qXkFAx28pT169ebFy1jOrtf\ntGhRvB+X/UvJSRjs5Cm1a9dWdc5jo0mTJvF+XPYvJSdhsJOnoH9o586dY/z52bNnVzXSY2L58uUy\ncOBAOXbsWKSPs38pOQ2DnTy5hzyms3acBo3JwaELFy6oWfg777yjTpFiOyNKBrB/KTkRg518O2uP\nzWx93759ah0dzp49K82aNZMaNWqoJwb2LyWnYbCTb2ftMZ2tw+HDh6N8bMmSJTJs2DB1H7tr2L+U\nnILBTr6ctcdmtn63YA+WJUuWSHvniXRisJMvZ+2xma3HJNixVIP98oMHD1br7kQ6MdjJd7P22M7W\nYxLsxqlWnDxFwa/Tp0/H6usTWYnBTp6ftWP9Oz6z9ZgGe/D2xxkzZsTq6xNZicFOnp+1N23a1Hwf\nIR/b2Xpsgz1Hjhw8oERaMdjJ89577z1JlCiRut+6detYz9Zv374tR48eve/n4ev26dNHdu3aJeHh\n4XEeL1F8sbojeRIOD2GWvWPHDlXNcdCgQXLkyBF57LHHVMs6BG+RIkUkadKk9/1ax48fv+8F0Xr1\n6snw4cNV/XUi3Rjs5AnXrl2TpUuXyoYNG2Tz5s2yadMmOXXq1D3/DkK9aNGiUrJkSXWrXr265MyZ\nM8rn4QnhbvLnz69KCFSrVs2S74PICgx2crW///5bHQxCT9HY7kS5ceOGehLADXC4qFatWtKhQwcV\n8gkT/m+lcvfu3VH+LrZRfvDBB/L666/HaNZPFEoMdnKlZcuWyaeffiqLFy9Wyy7BsJ6eN29eKViw\noBQoUEDSpUun1r8R1AhzHCTav3+/CmzUbz9//rz6e/g6CxcuVDcsqbRr107at28v69ati/T1mzdv\nrpZ2smbNGtLvmSimGOzkKlhe6dSpk3z99deRPo7wrlu3rlSuXFmFekwvkCLM//33XzVrnzt3rmzf\nvl19HMHfo0cP+eyzz+Tdd99VM/TkyZPLrFmzpFKlSrZ8b0RWYbCTayB4MYsObqSB054vvviiakEX\nlyURLL+g2xFuderUUTN4XFzFKwGs2x84cEBtj8TyzEcffRRlTzyRE3G7IzkeSuY2atRIXnjhBTPU\nsVccjaKxtv7ss89ats6NpRuU5kVXJTTtMGDmjguta9asseRxiOzEYCdHO3nypFpeMZZeMMNGadwv\nv/xSHn30UdseN02aNOriKNbxM2XKpD6G2fszzzwj8+fPt+1xiazAYCfHQqeiChUqqCP6Ro0XzNK7\nd+8uyZIlC8kYypUrp9bVUWvdqAdTv359+eqrr0Ly+ERxwWAnRzpz5oxUrVpVrXnDI488IlOmTJES\nJUqEfCxhYWHSt29fadu2rXof3ZJeffVVWbBgQcjHQhQTDHZyHGxJxPr2zp071fuFCxeW8ePHS/r0\n6bWNCUtAKEfQtWtXM9xfeukl+e2337SNiehuGOzkOAMGDJC1a9eq+3ny5FEnO2Paw9RuTZo0MWfu\nWJZBgbHLly/rHhZRJAx2cpStW7fKwIED1X1sLUSo40Kmk7Rq1UrtxDEabKAMMJGTMNjJUUswONV5\n69Yt9f6bb76pWs45DZZl3nrrLVUSGPDks2rVKt3DIjIx2MkxMFM3Tn5iNwrqtjgVXkWgW5Lhtdde\n45IMOQaDnRzh7NmzMmTIEHMJBssbmBk7GUoL1KhRwyxBMHnyZN1DIlIY7OQIU6dOlatXr5pr2Jkz\nZxa3tN4zTr3idOqdBcmIdGCwk3boUDR27Fh1H8W73NRWDlswcRoVsOd+5cqVuodExGAn/VasWCF/\n/vmnuo866DgQ5CYoQmbArJ1INwY7aWfM1u8MSbdAiz10UjIqUP7zzz+6h0Q+x2AnrXCCEy3tAOGI\n5hh2eOONN6Rs2bLmVkor4SLv888/b34/P//8s+WPQRQbDHbSKiIiwtwmiEbTdj4OuiIlTmxPC4Lg\nSpNGqz0iXRjspFVwCKIWuh0uXryouiTly5dP7JIrVy6z4iSDnXRjsJNWmzZtMu/bsQyDtnao5w7f\nf/+9lCpVSt1Q48VK6LNqrLOjzDB2+hDpwtZ4pJUxu0U/0Zw5c1r+9RHi586dU0XF0PrOaECN4mJW\nwxPTtm3b5NKlS2rpx65XIET3w2AnrQ4ePGguZWDWazU0uN67d68K9o4dO0rGjBnFLuHh4ZG+LwY7\n6cKlGNIKDaONGbtd/vrrL0mXLp2toX7n92CcoiXSgcFOjgj2JEmS2PYYOPxk54VTQ3BDbeP7ItKB\nwU5aGdsP7brYiN0w2BUTimAP/h7s2lZJFBMMdtLK2CKIbkR2MEoVhCLYg78HO5eWiO6HwU5aGXVh\nTp8+bdv6eqiCPfh7cFu9G/IWBjtpVaxYMfX26NGjcuHCBcu//vnz59XblClTit1Q3dFQtGhR2x+P\n6G4Y7KQV9plHF4xWeeSRR9TbPn36qGJj48ePlxMnTogddu/erd7mzZtX0qZNa8tjEMUEg520Klmy\nZJRgtBI6HL366qsqzKdMmSITJkww1/WtfmWAVx13fk9EOvDSPWkVXPjLjmBPmDChdOnSRd3sFPxq\ng8FOunHGTlrh4JCxXILToW5tCL18+XLzfpkyZbSOhYjBTtoZBbkQ6osXLxa3QW0YY9wojVCuXDnd\nQyKfY7CTdi1btjRPns6ePdt1DaEXLlxonjRt166dLTVviGKDwU7aoeJi/fr1zX3nqJDoFngSwpOR\nUVKgRYsWuodExGAnZ+jQoYN5f/Lkya6Zta9cuVIOHDig7jds2FAyZcqke0hEDHZyBqxLG7tJfvvt\nN1m0aJE4HbY4Dh482Hy/c+fOWsdDZGCwkyOgIfTnn39urk8PHz5cTp48KU42bNgws4xAp06dIh22\nItKJwU6O2tPeu3dvdR8VGQcOHOjYJZlffvnF3AmDJtmDBg3SPSQiE4OdHAVH/406K6tXr5avv/5a\nnAYnTPGkE3xNIBS1aIhiisFOjoKdJVOnTjXrmWNJBtsJneLUqVPqQm/wEkzFihV1D4soEgY7OXJJ\nZuLEieb7/fr1c0S4Hzt2TNq2bWvWhKlSpYoMHTpU97CIomCwkyM1a9ZMRowYYXYm+uCDD2TmzJna\n1tz3798vrVq1MptvP/744zJv3jx54IEHtIyH6F4Y7ORY2D44cuRI8/1PPvlEFfPCzDlU/vvvP5k+\nfbo0adLEfFxszVyyZImkSpUqZOMgig0GOzka1rC/+OILs+QA9rjjIBBmy3bP3nHwCLN0vHK4ceOG\n+litWrVUqKdJk8bWxyaKDwY7uaJI2KZNm+TRRx81i4V9+OGH0rFjR1m/fr3lAY+Z+ZgxY6Rx48ay\nY8cO9bEUKVLIqFGjZMGCBeo+kZOxHju5poUeQhwnPfv37y83b96UDRs2qFuOHDmkQYMGUqdOHUmd\nOnWcvj7W8fG1vv32W7XNEu8bsOtl0qRJkidPHgu/IyL7JAg49QQI2aJ58+bqRKcbjuzfzfbt29Xu\nlHXr1kX6OC5kli5dWgoWLGje7la7BdUYIyIiVHMPNMnYsmWLudvFkD59eunbt6/a3oiGHW40f/58\nqVevnvp+eaHXPzhjJ1fO3tGUY+PGjfLZZ5+pQ0wIruvXr8uaNWvUzZAhQwbVzAOhhnDGWvnVq1dV\niOPCaHTQKANh/tJLL0ny5MlD+J0RWYPBTq6F2Tn6mKJmCw41YffKzp07IwU2DhIZh4nuVzq4Zs2a\n0r59e9Z8IddjsJPrYVb+xhtvqNuVK1dUPffNmzerGy5+ou4MZukIfDSyxiw8PDxcVZM0btmyZdP9\nbRBZhsHuwyqKd1uC8ALsWHnyySfVjf63Dx/ceo2A4oY/bZ/BfnDsKCF/MH7WRu0d8gcGu88w2P0F\nF4vxM8crNfIPBrvPYH0Z69DkD7i2gOsK5C8Mdp/JnDmznDhxQvcwKETws86SJYvuYVCIMdh9Btv6\njh8/7tjORCSWl0fAz5z8hcHuM/glx7rruXPndA+FQoDB7k8Mdp8xfslDWfqWYqdSpUrqZgUGuz8x\n2H3GWG9lsPsDlt24xu4/DHafyZgxo6qbglrj5G3Y1nr48GF56KGHdA+FQozB7jOJEiVSPUVRopb0\nQwGzAgUKqCfbwoULy9y5cy372iingMJoqKlD/sLjaD6Efp0rV67UPQzfW758uWrmga5Mw4cPV+WU\n0foPM+38+fPH++ujrDHq0+OJg/yFwe7TYEc3IHQiSpkype7h+Nb777+vQhc1041aLngfdW6sCHY0\nJsFsHa/SyF+4FONDTzzxhCoOheqHpAf+/VFPHp2fggt04WeTK1cuSx4DM3Z8PfIfBrsP5cyZUx58\n8EFZsWKF7qH41qlTp9SSS3Q7VqzYxYKLpugQxSqX/sRg9yEUhGrUqJFqTuHlEr5O352E4lzYjnin\n6D4WW/jZok591apV4/21yH0Y7D7VsmVLOXToEGftmmDdG+vfs2fPjtQ4G+vi8d2Kiq+H5ttNmzZl\nn1OfYrD7VKFChaRs2bIyceJE3UPxLTTKRiNtNJtGc3HMstFnNb4nRfFkffDgQfXkTf7EYPexVq1a\nybx581QIUOg988wzMnPmTNm7d6+88MILMnToUPn000/jvSNm5MiR6qIp9sWTPyUIsMyfb127dk1K\nlCih+n/+8MMPbMbgAXPmzFE7bX788UepXr267uGQJgx2n/v111+lfPnyahkAa7LkXmfOnFFLbDVq\n1JApU6boHg5pxGAn6dSpk1oSwBH07Nmz6x4OxQF+jfHEjNOsu3btknTp0ukeEmnEYCe5ePGi2qGB\nfdVLly6VPHny6B4SxQK2rHbs2FEmTJigrpnUqVNH95BIM148JVVPZNWqVZI2bVp56qmnZNu2bbqH\nRDGEIl/GmQSsrzPUCRjsZPZC/fnnn1WtEqy5Dxo0SC5duqR7WHQXeKG9bNkyKVeunLpQunjxYrVt\nkggY7GQKCwtTIdG5c2cZOHCg5M6dWz7++GM2v3bYDB3r6OiwVK1aNXW6FBfAK1eurHto5CBcY6e7\n1jIZMmSIjB49Wq5evapCHnujy5Qpo+rMYPkGR+LJPjhBigqc58+fV8tjKOq1ZcsW1bMWwd6/f381\nYye6E4OdiMhjuBRDROQxDHYiIo9hsBMReQyDnYjIYxjsREQew2AnIvIYBjsRkccw2ImIPIbBTkQk\n3vL/APQ/aSRQT1CKAAAAAElFTkSuQmCC\n"
        }
      },
      "cell_type": "markdown",
      "id": "df8ef4c1-f4e9-4155-aec0-df2dee4c622e",
      "metadata": {},
      "source": [
        "Try changing the data to `t = jnp.array([-20, -10, 30])`. Notice how\n",
        "changing the data for one group updates the expectation of the\n",
        "population.\n",
        "\n",
        "## No pooling\n",
        "\n",
        "In contrast to complete pooling, which assumes all observations are\n",
        "generated by a shared population-level distribution paramaterized by\n",
        "$\\theta$, the no pooling model treats each department as having its own\n",
        "independent parameter $\\theta_d$. The key change is that instead of\n",
        "having a single $\\theta$ drawn from the prior distribution, we now have\n",
        "three independent parameters $\\theta_G$, $\\theta_E$, and $\\theta_M$. The\n",
        "practical meaning of this change is significant:\n",
        "\n",
        "In complete pooling, an observation of the Math professor being late\n",
        "would shift our expectations for all departments. In no pooling,\n",
        "learning that the Math professor is late only updates our beliefs about\n",
        "the Math department. Each department’s parameter is estimated using only\n",
        "data from that department.\n",
        "\n",
        "The key change in our model specification is that we have replaced the\n",
        "single $\\theta$ with three parameters, $\\theta_d$.\n",
        "\n",
        "$$\n",
        "\\begin{align*}\n",
        "\\theta_d ~\\sim&~ \\mathcal{N}\\left( \\dot{\\mu}{=}0, \\dot{\\tau}{=}20 \\right) \\\\\n",
        "t_d ~\\sim&~ \\mathcal{N}(\\theta_d, \\dot\\sigma{=}15) \\\\\n",
        "d ~\\in&~ D,~~~ \\text{where}~~~ D = \\{\\text{G}, \\text{E}, \\text{M} \\}\n",
        "\\end{align*}\n",
        "$$\n",
        "\n",
        "![](attachment:generated/multilevel-models/no-pooling.png)\n",
        "\n",
        "The plate (rectangular box) with label $d$ indicates that everything\n",
        "inside is repeated for each department $d \\in D$ Variables inside the\n",
        "plate have a unique value for each department. Variables outside the\n",
        "plate are shared across all departments.\n",
        "\n",
        "To change the `complete_pooling` model to the `no_pooling` model we make\n",
        "a simple change, adding the conditional statement\n",
        "\n",
        "``` python\n",
        "president: observes [department.d] is _d\n",
        "```\n",
        "\n",
        "and the query variable `_d` in the definition, `[_d: Department, ...]`."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "4bca2eca",
      "metadata": {},
      "outputs": [],
      "source": [
        "%reset -f\n",
        "import jax\n",
        "import jax.numpy as jnp\n",
        "from memo import memo\n",
        "from enum import IntEnum\n",
        "from jax.scipy.stats.norm import pdf as normpdf\n",
        "from jax.scipy.stats.cauchy import pdf as cauchypdf\n",
        "from matplotlib import pyplot as plt\n",
        "\n",
        "normpdfjit = jax.jit(normpdf)\n",
        "\n",
        "class Department(IntEnum):\n",
        "    GOVERNMENT = 0\n",
        "    ENGLISH = 1\n",
        "    MATH = 2\n",
        "\n",
        "t = jnp.array([1, -10, 30])\n",
        "sigma = 15\n",
        "\n",
        "Theta = jnp.linspace(-40, 40, 200)\n",
        "\n",
        "@jax.jit\n",
        "def professor_arrival_likelihood(d, theta):\n",
        "    ### likelihood of a professor from department d \n",
        "    ### showing up t_d minutes early/late, \n",
        "    ### for a given theta and sigma.\n",
        "    return normpdf(t[d], theta, sigma)\n",
        "\n",
        "@memo\n",
        "def no_pooling[\n",
        "    _d: Department,\n",
        "    _theta: Theta,\n",
        "](mu=0, tau=1):\n",
        "    president: knows(_theta)\n",
        "    president: thinks[\n",
        "        department: chooses(d in Department, wpp=1),\n",
        "        department: chooses(theta in Theta, wpp=normpdfjit(theta, mu, tau))\n",
        "    ]\n",
        "    president: observes [department.d] is _d ### new conditional statement\n",
        "    president: observes_event(wpp=professor_arrival_likelihood(department.d, department.theta))\n",
        "    return president[Pr[department.theta == _theta]]\n",
        "\n",
        "mu_ = 0\n",
        "tau_ = 20\n",
        "res = no_pooling(mu=mu_, tau=tau_)\n",
        "\n",
        "### check the size and sum of the output\n",
        "# res.shape\n",
        "# res.sum()\n",
        "# res[0].sum()\n",
        "# res[1].sum()\n",
        "# res[2].sum()\n",
        "\n",
        "fig, ax = plt.subplots()\n",
        "\n",
        "ax.axvline(0, color=\"black\", linestyle=\"-\")\n",
        "for d in Department:\n",
        "    department_name = Department._member_names_[d]\n",
        "    department_abbrev = department_name[0]\n",
        "    theta_posterior = res[d]\n",
        "    theta_expectation = jnp.dot(Theta, theta_posterior)\n",
        "    ax.plot(\n",
        "        Theta, \n",
        "        theta_posterior, \n",
        "        label=(\n",
        "            rf\"$p(\\theta_{department_abbrev} \\mid y),~ \" \n",
        "            + r\"\\operatorname{E}\" \n",
        "            + rf\"[\\theta_{department_abbrev} \\mid t]={theta_expectation:0.2f}$\")\n",
        "    )\n",
        "_ = ax.set_title(r\"Posterior of $\\theta_d$\")\n",
        "_ = ax.legend(bbox_to_anchor=(0.9, 0.5), loc='center left')\n",
        "\n",
        "_ = plt.suptitle(fr\"$\\dot\\tau$ = {tau_}\", y=1)\n",
        "plt.tight_layout()\n",
        "plt.show()"
      ]
    },
    {
      "attachments": {
        "generated/multilevel-models/partial-pooling.png": {
          "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAF2CAYAAAB6XrNlAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90\nbGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAO\nxAAADsQBlSsOGwAAOnRJREFUeJzt3Qd4FFX7NvADoUgLvQgaCEWkKx1EqoCAAkpTREB6ERFUBMRC\nR4qFjnQEldfwUhSQpryA9CJFSqT3Fnpv81338Zv5b0jbJDN7pty/61qzCcnOWTN59uwz5zxPEk3T\nNEFERK6RVPUAiIjIXAzsREQuw8BOROQyDOxERC7DwE5E5DIM7ERELsPATkTkMgzsREQuk0z1ALwG\n+8FOnz4ttm7dKv766y8REREhbt++Le7fvy9SpkwpUqVKJZ588klRqlQpUbJkSZEpUybVQyaKAufr\n33//LbZt2yb2798vbt68Ke7cuSP/7YknnhBp0qQRzz77rDyPixQpIpInT656yJ7CwB4AJ06cELNm\nzRLr16+Xfwjnzp3z+2dDQ0PlH0e1atXEW2+9JdKnT2/pWIlimpD8/vvvYt68eXJSsmvXLnH37l2/\nfhYTluLFi4vSpUuLRo0aierVq4skSZJYPmYvS8KSAtZ49OiR/EMYP368WLhwofw8sVKnTi1atGgh\nunTpIkqUKGHKOIlic+XKFTFz5kwxYcIEceDAAVMes2DBgvIcbtmypciQIYMpj0mRMbCbDAF82rRp\nYsSIESI8PDzSvwUFBcm3pZiBY/aCjyEhIfKta7JkycS9e/fErVu3xD///CNn9vrt4MGDUY5TsWJF\n8dlnn4natWsH8NmRV5w9e1Z8/vnnYvbs2fKc9IVgjHNXv2GSkTFjRnkeA1Iyly9fFjt37ox0HuNF\nIrqJSv/+/UWOHDkC+vxcD4GdzHH48GGtatWqeKGMdCtVqpQ2bdo07caNGwl63HPnzmlffvmllidP\nniiP3aZNG+3y5cumPxfypkePHmlz5szRMmXKFOk8S5kypdaqVStt48aN8nsS8rgbNmzQWrZsKR/L\n97FxrB9++CFBj0vRY2A3wcOHD7WxY8dqadKkMU7WFClSaK1bt9Y2b95s2nEePHigLV68WKtbt26k\nP4xcuXJpS5YsMe045E1nzpzRGjZsGOncwmRixIgR2sWLF007zoULF7Thw4dHmajg2BgDJR4DeyJF\nRERo1atXj3SCli1bVtu7d6+lx12xYoUWEhIS6bgdO3bU7t27Z+lxyZ0wMfCdpSdNmlTr3bu3dvv2\nbcuOicfGMXAs39n70qVLLTumVzCwJ8KpU6e0IkWKRHq7ipTJ/fv3A3L8q1evah06dIgU3OvXr2/p\nHyO5z48//qglS5bMOIcKFy5s6jvNuGzatEkrVKiQcXyM5aeffgrY8d2IgT2BTp48qeXLl884GYsW\nLWr5LD0my5cv17JkyWKM5aWXXmJwJ7/MmDFDS5IkiXHu9OjRQ8m5g2Pi2Po4MKaZM2cGfBxuwcCe\nwBwhZjX6SVi+fHmZklEJLyrItetjatCgQcDeOZAzzZs3L1IaZOTIkaqHJPP5vukgjJHij4E9ARcw\nX3jhBePkq1Chgnb9+nXNLqtyfIP7u+++q3pIZFNbtmyRF/j1c2XMmDGaXWAsvosQMFaKHwb2eMLV\nfP2kK168uHbp0iXNTjBzz5w5szHGVatWqR4S2QzSHr7vOAcNGqTZDcakjw/Xse7cuaN6SI7CwB4P\n+/btM9bgYmkjZsh2hCWR+h9F7ty5tWvXrqkeEtlInz59Ii0xtOP6cYwJ6UR9nH379lU9JEdhYI9H\nCga5dP1Ew7p1O8Maen2snTt3Vj0csgmsdtHz6lhaaOd146dPn9YyZswoxxoUFMSUTDwwsPtp1KhR\nRqDE7lJsSrIz7EbNmTOnMebVq1erHhIphovpvstzscPU7mbPnh0pJcMFAf5hYPfD3bt3tezZs9s+\nBRNbSqZ27dqqh0OKhYWF2T4FE1dKhqtk/MNGG36YP3++UWq3U6dOspSuE9StW1dUrlxZ3l+2bFm0\nxcTIO1BpVDds2DBHlM7FGDHW6J4DxYyB3Q++JxMCu5OgPKpu4sSJSsdC6uzbt0+WkYYaNWrI0rlO\ngYYdqOEOq1atko09KHYM7HHYs2ePWLNmjbyPErn58+cXTvLaa6+J7Nmzy/soJ4xuTeQ9vi/qvi/2\nTsEJSvwwsMdhypQpjv6DSJEihWjfvr28jxrZ6IBD3oJOR2iWATlz5hT169cXToMxY+wwY8YMv7s3\neRUDexz02XrmzJlFvXr1TH3szp07yxzixYsXo/xb1apVRZYsWUw5TuvWraM8H/IOtLG7evWqvI/G\nFmjq4jTomYqxA57L7t27VQ/J1hjYY4FZAVIxUKZMGdkByUzoMJMrV65oAzj+GM1qf5c3b175wgTo\nZEPe4vs7L1++vHCqcuXKGfd5HseOgT0WmBWgGzugBZiZsNQUjx9d8D527JhMm5gV2PGuQB8/jsm3\nsd7iGwTNPI/R1BrnVly3v/76y5Tj+Y6dgT12zntP5oI/CDh06JC4ceOGeO6556KdyQM6u5sF41++\nfLl8ocK7ELOfD9n/PMY7w6efftrU9Aj6ovouC8Y7zY8//tjof4rAXqxYMVOOh/7AeOcZERHBwB4H\nBnZFgV0P3tHNymP7t4RC82zf58XA7g1okK7no/E7N3PtOs5P33N0xYoVstG177pzM2HsOI+xJwPP\nCc8NiwMoKqZiYqFvSsLFJjNnOnEFb7x1xTELFy5s2vGQZ3/8eZH7IaX34MEDeT9fvnyWHuvvv/8W\nRYoUsfQY+nmMd55Xrlyx9FhOxsAeizt37siPqVKlMn2XHgI7HrdAgQJR/m3z5s1yA0nKlClNOx6O\n9fjzIvfz/V37ngNmO378uFytUrRoUWElnsf+YWCPBd7qgRVv9xDYc+TIIZImTRrlgtTJkydNTcOA\n74sEL5567xwGK9MWerrH6sDu+xx4HseMgd2Pk8j3j8MMeAuJlS+nTp2SNx1mPPomKDPTMI8/BzPf\nCZC9+QZCs89jX/qyYKsDO89j//DiaSz0K/vYho/liWalY7ByANKnTy/Xxzds2FDcunVL/PbbbyJ1\n6tTGCgPkE998801TjulbSkB/XuR+vr9rK8tJhIeHy4+FChUSVvJNv/A8jhln7LHIli2b/IiLT6dP\nnzb9wunkyZNF2bJl5RZpLEVs27at/IhNS+fPn5epGrMcOXLEuJ81a1bTHpfsDatU9I11R48etew4\n+s7WdOnSCSvp5zEWF2BiRNHjjD0Wzz//fKQlggi4Zgb2KlWqiAYNGkT5d+TYrVy6WbJkSdMfn+wJ\n6QqsVMG7RJwDZr7z9KUXx8O7T+xu/fDDD0VwcLCpx8DY9fMYz4mpmJhxxu7n2m9c1DQLAjuWT2I2\nFSj6HwRmOmZufCLnnMdY5up7TcdMH3zwgXj55ZfF2rVrxahRo0TatGlNPwbGjneywH0YsWNgjwUC\noF4wyaydbg8fPpQXmgIZXDHT0V+YcHGLuUlvCcRWfKT3li5dKvP4N2/ejLLaywy+kyvfSRdFxcAe\nCwRAfcPFli1bxKNHj0y5yIQLQIEM7FiBc+HCBXmfMx3v8f2db9q0STgV9nfoeB7HjoE9DpUqVZIf\nERixlTmxsGoAM+ghQ4aIQJk1a1aU50PegT0Rempkzpw58l2j02DMs2fPNi7QMp0YOwb2OLRr187R\n/Rax9XrSpEnyPi5mNWnSRPWQSME7z5YtWxo7RJcsWSKcZvHixeLEiRPyPp4L04mxY2CPA6ovVqhQ\nwTi5rFwyZoVffvnFWKqJhhtp0qRRPSRSAE1dnDxB8R2z73Oh6DGw+0HfDYoUij77dQr+QZB+0fzF\nF1+U97ERDmWjneLgwYNGGrRy5cqWFxpzAwZ2PzRu3NjocoRAqb8ltDt0dMcN0OUd3d7Ju3x79n7y\nySfCKfr16+fovsNKaOSXIUOGaPjfhVvt2rW1R48eaXZ27do1LXfu3MaYly1bpnpIpNjdu3e1AgUK\nGOfEvHnzNLsLCwszxoux4zlQ3BjY/XTv3j2tZMmSxkk2ZcoUzc46depkjPWdd95RPRyyiXXr1mlJ\nkiSR50W2bNm0CxcuaHZ1/vx5LWvWrHKsGDPGTv5hYI+HnTt3asmTJ5cnWnBwsHb8+HHNjlauXGkE\n9Vy5cmmXL19WPSSykR49ehjnR7NmzTS7atq0qTHOnj17qh6OozCwx9OAAQOMk61MmTIy5WEnBw8e\n1HLkyGGMccmSJaqHRDZz8+bNSCmZr7/+WrObr776KlIKBmMm/zGwJyAlU7p0aeOkq1atmnbr1i3N\nDk6ePKnlyZPHGFvbtm1VD4lsCmmNoKAg41yZOnWqZhcYiz4ujPHPP/9UPSTHYWBPgDNnzmj58+c3\nTr6qVasqn7kfOnRICw0NjXSB986dO0rHRPY2e/ZsI9+OjxMmTFA9JG38+PGRxjRnzhzVQ3IkBvYE\nOnLkiPb0008bgbRUqVIyuKqAGc2TTz5pjOWFF17Qbty4oWQs5CwIpPp5g9unn34q35UGGo7Zr1+/\nSGPB2ChhGNgT4ejRo5Fm7mnSpNHGjh2rPXz4MCDHR94RF5X0GQ5uNWvWZFCneKc+kiZNapxDWP21\na9eugB0fx/JdcYax2Ck15EQM7Il09uxZrXz58pFmGkjNHD582PIcqe8FMNxatGjB9AslCNa0p0uX\nzjiXsPpr0KBBls7e8dgDBw40VprhhjE4YX293TGwm+DBgwfa8OHDtZQpUxonaOrUqbWuXbtqe/bs\nMe042BS1evVqrUmTJpFm6bjhbSxRYhw4cECet77nVaFChWRKxMxrSFevXtXGjRsnH9v3WLVq1dKO\nHTtm2nG8jIHdRHv37tXKlSsX6WTFrUqVKtrcuXMTvGsOfwhI8RQuXDjKY+u3kJAQ5Rdwydk+/vhj\n43xKlixZpPMrbdq0WpcuXbTdu3cn+PHxs507d5aP5fvYmKVPnjzZ9ru5nSQJ/qOmmIE7oW702LFj\nZXuwx2vKpEqVSlaLRJMAdIDBx5CQEFmCNHny5OLu3buyA80///wju8Wg2w1uf//9t2yo/XiN7c8/\n/1yWYJ0yZYr82rvvvivGjBkT0OdL7oDzrFy5cvL8zZw5s2yqjnP4559/lqWffaH3L85d/YZzGm0e\n9VK6aCRz5coV8ddffxnnMG6Pt+XDOd+0aVPZmwB/B2Qi1a8sbnX//n1twYIF8u1lTLPs+N5SpEgh\n8+jr1683ZjfYVZozZ05jedjatWtVP3VyGOS6S5QoYZxnWAbpew1p8ODB8h2hWecxHgu1l/DYZA3O\n2AMAM/Dp06eL9evXi+3bt4vr16/7/bPZsmWTs6Jq1arJeuroLfm4RYsWiQYNGsj7BQsWlDMlNiIg\nfw0ePNiooFivXj1Zwz9JkiSRvgczefQjmDdvnnw3uX//fr9bRaL/KSqL4l1qo0aN5DGCgoIseS70\nLwb2AMMfAwI93pru2LFDXLp0Sb51vXfvngzGuOXIkcNI1eBt7+N/ZNF54403xNy5c+X9Pn36BLT1\nHjnXvn37ZCoF5x9aziHt9/TTT8f5c2hYradaEORv3bol04h6yjF16tQymOupGjZ4CSwGdpc4f/68\nKFy4sIiIiJCzITT+LVmypOphkY1hFo7mGxs2bJCfT5gwQXTq1En1sMgEbLThEkjZfPvtt8YfbNu2\nbaNc9CLyNW7cOCOoV6lSRXTo0EH1kMgknLG7CH6Vr7zyitGsGLnTvn37qh4W2RB696LFHFIoSP/t\n2rVLFChQQPWwyCQM7C6DJZb4g8UF2hQpUoidO3eyJR5Fgj/5WrVqiZUrV8rPR4wYIT788EPVwyIT\nMRXjMrjwNXz4cHkfF8SQkvF39QJ5w4wZM4ygjov077//vuohkck4Y3chBHIsj1yzZo38fPTo0aJb\nt26qh0U2cObMGXmRHRuIkiVLJpffFitWTPWwyGScsbsQ1g1jN6q+lh3LH5FTJW/DHK5Lly4yqAOu\nvzCouxMDu0vhQtiAAQOMNcdY8cA3Z96GzUULFiyQ9zFr54V192IqxsVQX6ZChQpypyBg9yt2r5L3\nYH8Dgjn2O2DDG5Y5ojYMuRNn7C6GHOrUqVPlR+jRo4c4e/as6mGRAj179pRBHXCxlEHd3RjYXa54\n8eIyxw7IrXbt2lX1kCjAfvvtNzFr1ix5P2/evGLgwIGqh0QWYyrGA1AOGOUF9u7dKz8PCwuTxZjI\n/bCfoWjRouL48ePy81WrVonq1aurHhZZjDN2D0iZMqVMyejFxDBrR/Excj+8W9ODert27RjUPYKB\n3SPKly8vunfvLu+fO3dO5lzJ3dauXSvrwUDOnDnlDlPyBqZiPATLHrFu+ciRI0butXbt2qqHRRZA\nCV102UKJaFi4cKGoX7++6mFRgHDG7iGoiT158mTjc6xtj0/TD3IO7GHQgzpq9TOoewtn7B7Uvn17\n9kl1MZQJKFu2rNG/FBfNUdaZvIOB3YOw7BGbVVA3RM/FVqpUSfWwyASowV+mTBlZ1RNmz54t3nrr\nLdXDogBjKsaD0FEe3XJ0qACJ9nzkfKjsqQd19BZt3ry56iGRApyxe5hvn9TevXuLoUOHqh4SKehf\nSu7DwO5h7JPq3v6lEydOFB07dlQ9LFKEqRgPY59U9/YvxQVy8i7O2D2OfVKdj/1L6XEM7MQ+qQ7G\n/qUUHaZiKNo+qUjNkP2xfylFhzN2ktgn1XnYv5RiwsBOBmxBR/12rGlH+YE9e/aIPHnyGIE/PDxc\ndmO6ePGiuHHjBmf1AZA6dWoRHBwsChUqJEqVKiV/L4A/W5Renj9/vvz8s88+E/3791c8WrILBnaK\nBDnaXr16yfuYwb/88svi999/F5s2bZIzw+TJk4ssWbLIAIP7ZB28mKKYF/6/X7t2TS5JxQsv2h3i\nBVf/PeH6yLZt22R5ZiJgYKcofVKxJf2vv/6Sn2OVxWuvvSbL/iKgoGIgLrBS4OBPFDXVN27cKG+o\nyrl//37j3/E1trojX/82wyT6//7880+5SkaHID5q1Cjx5JNPKh2Xl6FBSu7cueWtWbNmckOZb2DH\n7mG8GCdNyrUQ9C+eCWRYtGiRrM+OHYx6Iw6kAFABkuwBs/Xvv/9e3g8NDZUbk8aOHStatGghVzQR\nSUjFEM2YMUMLCgrS2rVrpz148EC7c+eOVrhwYaTp5C0sLEz1ED3v2rVrWkhIiPE7WbVqlfz68uXL\ntTRp0mh16tTRbt68qXqYZAMM7KRt2LBBS5IkidarVy/t0aNHUb6OIJI9e3YtIiJC6Ti9rmvXrkZQ\nb9++faR/27Rpk5YxY8YoXydvYmD3OH1mXqNGjUhBXff+++8bwaRVq1ZKxkiatmbNGuP3kDNnTu3y\n5ctRvufHH3+MNJMn7+KqGI/7/PPPxciRI8Xu3btF3rx54+yTunTpUrkEkuzXvxR/yg0bNpT7D1Av\nRl/zTh6k+pWF1AkPD9eSJUumff3117F+38qVK43ZInK8yPVS4PTu3dv4///GG2/E+r0nT57UgoOD\ntT59+gRsfGQ/DOwe9sEHH8hAjYulccFFVT24INdLgbFt2zZ5URv/3zNnzqydP38+zp/54osv5Pci\nzUbexOWOHoWlcTNnzhRt2rSROxr92ZGqr2XHEjv0SSVroTY+fj966QbU78maNWucP4efuXTpkly+\nSt7EwO5R+KPHRpd33nknQX1S27Vrxz6pAe5f+uabb/pdrRPXQaZMmWLxCMmuePHUo+rUqWNcDI0P\n9kl1Rv/S//73v6Jx48bi8OHDRiE38g7O2D0IM+1Vq1bJIB1fSAdkzpzZSM+gVCyZS29TqO8kxf/n\n+DalRlesVKlSieXLl1s0SrIzBnYPQoEv5G9R1Cu+2CfVGf1LUeMHjTdQlZO8h4Hdg/DHnjFjxgT3\nxWzevLmoW7eu8SKhd18ic/qX9unTx6isiTx5Qot7oeIjA7s3MbB7kF7mFVUDEwI/N3HiRJn7hQED\nBkSqNkgJg8tdHTp0kE2pYeDAgSJ//vwJfjyUWt67d68s5EbewsDuQZjFJbZ+N/ukWtO/dMWKFfI+\nyvAmtn8pfsd4sdiyZYtJIySnYGD3GGxPR3kArLhILMwuK1euLO+vX79ejB8/3oQRerd/qV4qGf1L\np06dKj8mRq5cuWS3K5QYIG9hYPeYs2fPyo85c+ZM9GMh94scMHLBgNwwcsQUP5hVd+3aVbbAg759\n+5rWlBq/53PnzpnyWOQcDOweDezZs2c35fFwARY5dr1gGGbx3BoRP/PmzTOaUqN/KQK7WXLkyGH8\nzsk7GNg9xuzADj169BClSpWS95EjRqkC8g+2/mO2rr8DQgrGzKbUDOzexMDuMXhbjvIAevrEDMgF\nT5s2zcgJI9AjZ0xxw/+r8+fPy/u4WGp2U2q8gDMV4z0M7B6D2ZuZs3Vd8eLFjfXXyBWzT6p//Utn\nzZol76MWvp7SMhNn7N7EwO4xly9fFpkyZbLksT/55BNRuHBho1YJcscUvevXr4uOHTsan0+ePNmS\nxhj4XSPdQ97CwO4xWHOO7eZWQG4YOWJ94xNyxwwq0cO7m+PHj8v7KBlQvXp1S46D3zVLPngPA7vH\n4I/cqsCu73bs3r27vI/crr42m/4PatmjHoy+HNHKkgzJkyeXG8cePXpk2THIfhjYPebBgweJ3vgS\nl0GDBonQ0FB5HytkkEum/9sghlr2OtS4x8VsKwM7cNbuLQzsHpTQGjH+Qq4YOWMdcsnIKdO/dXXC\nw8PlfZRNjq4ptZN+12RPDOxkiRo1asj6MYBcspmbbpwKtetRWx1Q0x617YmswMBOlhk5cqTRJ3Xs\n2LFi3bp1wquQCvEtlOZv/1KihGBgJ8s83icVgc2rfVIxU0ft+vj2LyVKCAZ2slSDBg1Es2bN5H3k\nlvv37y+82L9Uf96oYY8XO+a+yUoM7GQ5pB30TVFe65NqRv9SovhiYCfLeblPKmrUJ7Z/KVF8MbBT\nQLz11luR+qTqq0PczMz+pUTxwbOMAuLxPqnIObu5T6revxQ16s3oX0oUHwzsFDBe6pPq27+0dOnS\nie5fShQfDOwUUF7ok/p4/1LfWvVEgcDATgHlhT6pqEVvRf9SIn8xsFPAublPalhYmKxFb0X/UiJ/\nMbCTEm7sk+rbvxQXi83uX0rkLwZ2UgI5ZwQ+N/VJtbp/KZG/GNhJmRIlSrimT+rj/UuxvJFIFQZ2\nUsqfPqkoHKbPhFXCEs3oGkMHqn8pkb8Y2Emp2Pqk4oLq9OnTZenfXLlyiTVr1igbJ9bbly1bVo7l\n7bffFhEREdH2L0V3JKv6lxL5i4GdlIuuTyqWQNauXVu0adNGpmnQ0u+XX35RNsZDhw6JnTt3yvuz\nZ88WhQoVEv/5z3/ki41v/1IvlEog+2NgJ1t4vE8qAqe+c1N34sQJRaOLeuwLFy7IcsR6/ZtA9C8l\n8hcDO9kCctL9+vUzPo+uIYedArtOrwVTpkwZ8eqrrwZ4VETRY2An5VDCd8iQIaJz586xfp8dA7tu\ny5YtombNmuLw4cMBGxNRTBjYSSnkrXFREqtj9GYUMTl9+rSyomH+vKisWrVKlg9A7Xm37KQlZ2Jg\nJ2Xu3r0rXnrpJaMXaFwQ1FVtYjp58qRf33fr1i25Oennn3+2fExEMWFgJ6UFwVKnTh2vn1GVjonv\ncbmOnVRiYCdlkidPLtatWyeaNGnimsCeJ08esWDBAlGvXj3Lx0QUEwZ2Ut58A+vBf//9d1G0aFHL\nAjty3o8ePUrQz2Jn6dWrV2P9nlSpUskyAnv37hUNGjRI0HGIzMLq/2QL1apVEzt27JDt8z799FOj\nnnl8AjuCNzY2bdu2TWzdulV+3L17t7h27ZpcPol/x7sEBGHUc0F1Sf1WvHhxo0Z8fI4JWM+OjUl4\nkSKyBY08pVWrVlrdunU1Ozt//rzWoUMHLUmSJFhaEulWs2bNKN+/f/9+rXv37lqWLFmifL+/t2TJ\nkmm1atXS5s+fr92/fz/S4y9dujTanylWrJi2evVqzc4WLFggx3rnzh3VQ6EA4oydbCdr1qxi0qRJ\nsrBWt27dZAs93fbt2+VHlBhYtGiRbK2HZYbRSZEihVx+mC1bNjlLDwoKkjN3pFb+/vtvuXtUh8db\nvny5vD311FPy2Kj7kiNHjiiPnzFjRpl2wfew5R3ZEc9Ksq2SJUvKi6uozdK+fXu5PDJTpkyy1ACC\n6pEjRyJ9P2q11K9f30ivoIMRgnt0kJbBEkaka3D73//+J9auXSv/DV9HOghdnlC3Jnfu3MbP4UIv\nXkyyZMli8bMnSrgkmLYn4ufJYVq3bi1nqosXLxZOgouXM2bMkGve8dFXjRo1RJcuXeSWfuTQEwoX\nPlHvBbVqMKvXoW5N27ZtRYUKFUTFihWFkyxcuFA0bNhQvlNhNyfv4KoYcoTNmzeLr776KlJQx4vU\nvn37xMqVK8Xrr7+eqKAOqAs/ZswYucMVs/LMmTPLr+MYvXr1kqmf6GrYENkNAzvZ3rBhw0StWrWM\nmufIgaNjEWq1P/vss6YfL23atLJuDfLweMEALJX88ssvReXKlcXFixdNPyaRmRjYybaQJezdu7fR\nPg9wQXPPnj2yVrvVsmfPLsLCwsRPP/1kzN5R7KtKlSri1KlTlh+fKKEY2Mm2+vbtK2fJgDTL3Llz\nZdu59OnTB2wM6OyEdepYD48erXouHnl9O7TrI4oOAzvZEgI6UjCApYrontS0aVNl40FLvNWrV8sL\nqHDgwAH5riGuHalEKjCwk+388ccfMgWjz9Tnz58fkNRLXNAdCbl9LMMErNBBj1Yiu2FgJ1u5ceOG\n7HOqmzVrli2Cui44OFgGd5QkgDlz5sglhUR2wsBOtoKZOuq9QMuWLcUbb7wh7LgzFssukX8HbJaK\niIhQPSwiAwM72SoFM27cOCOn/c033wi7evHFF2W5Azh37pzo3r276iERGRjYyRawTrxTp07G5999\n952syWJn6NOaL18+IyWDFyYiO2BgJ1tA8a3w8HB5H+mXV155RdgduiRh+aUOu1aJ7ICBnWwBW/h1\nH3/8sXBSHfnSpUvL+7iI6m9vVCIrMbCTcrhY+uuvv8r7WCf+3HPPCSdBATI9nYQUEpFqDOykHIKh\nXmRUD5Jm279/v9zghNIAWLKIVM+xY8dMeWzsTNWvByA1c+/ePVMelyihGNhJuR9++EF+RI3zxo0b\nm/74aJSBTUWbNm2SFSHRaHrJkiWysJgZQTh16tTinXfekffPnj0rd6gSqcTATkqh3oo+c8ZGpJj6\njiYU+pU2atRIpndQPGzUqFHixx9/FD169JAXa83aXIQGHzoUCiNSiYGdlEL3Ih26Hpnto48+Erdu\n3ZJdmNKlSxcpfQK+bfcS4/nnnzc2LPk+JyIVGNhJKd8gqK8uMfOi7M8//yxrquslAHRosQdmVWhE\n3v6ZZ56xbWBPmtS6P/VDhw7JMhAhISGytg9e4B6/7dixw7LjU1TseUpK6UEQf/yY9ZoJZX6xUqV5\n8+ZR/g39U8HMdnF4x4Gqj2gIgmYcduqLqr+bMBuai2PJ582bN8Vrr70mChQoILtd6Q3A8+TJI4M9\nulNR4DCwk1LoUgQICOhcZPamJ0CTagQgX/p6czTANgsu0OoXgpHPr1q1qnBzYMeLI1Ya3b59WwZy\nNCDR4V0SqnKOHTtWXqymwGJgJ6X0ptEorGUmLJ/EzBFGjhwZ4/eZ2VrPd4aOKpV2YkVg//7772Ua\npl+/fpGCum9gxwsqA3vgMcdOSunNoc1eDYMZOYIrKkQiyD9+09MzeuMMM6AhiA6zWLfD9Qu8YES3\n90B/9/Xw4UMFIyMGdlIKOXArZpRnzpwxqkRGd0ykDnBBVS/iZfYFSi8EtI0bN4oiRYpE+/9Y7wkb\nGhqqYGTEwE5K6TN1/WKmWe7fvy8/pkiRIsq/Iaij1G50F1UTw3eW7jt7d6NLly6Ja9euiRw5ckT7\n72hGAna6zuAlDOyklB4A9Vy7WfScfXTLGVFuF7tF9XrqZvF9Dm4P7Po7LAT3x+3bt09euMZqmdy5\ncysYHTGwk1J6KgTBQJ9lmyF//vwiW7ZsYtGiRZEuZA4ePFhu+cdH/LuZdu/eHen4dqLX4jELauM8\n9dRTcrnqwYMHja+juTeuayDdhYbkpAYDOyml7zZFKmbv3r2m5rt79uwpc+24QIpSwDVr1pQrODp0\n6GBJxyN9TT6Cnt1yy/q1DDPh/ymuJbzwwgvi/fffF++9955cZYQm3yiGVqZMGdOPSf7hckdSyreM\nAAJjiRIlTHvsXr16iQcPHoiJEyeK0aNHy00yM2fOlDNKs+Hdxq5du4z17FZtCLLLjB3effdd+QKK\nteoTJkyQJRuw7LFPnz6m7yKm+GFgJ9sE9q1bt8qt6WZBcP3kk0/kLRAbrfQLwFbUvLFjYAcsdbSq\n1DIlHFMxpDzHniFDBnkfPUOtCkBW+/333437dgzs5C0M7KQUZtX6zkQ0w1izZo1wGuSvke7Ra8/U\nqFFD9ZDI4xjYSTnft/K+vU+dNFv/559/jEbc6NJEpBIDOymHVSv6RdP//ve/xq5Rp/B9MeratavS\nsRABAzsp51tvRF/F4hTo/qR3YcJKEC7xIztgYCdbwPb+9OnTy/vDhw+XbevsDhd6O3XqZKwR5+oQ\nsgsGdrIFVAMcMGCAUfERzaHtXkhrxowZRk0U9FRt0aKF6iERSQzsZBvY8PLiiy8avUi//fZbYVco\nC4zdlpAsWTIZ5NEpiMgOGNjJNrCLcdq0aUYBLWwsMrPMgFnwTqJdu3ZGAaxPP/3U1B2zRInFwE62\nguJZw4YNM1IyL7/8suwhaqe8Ola+LFu2zEjBYAs9kZ0wsJMtUzKvvPKKvH/ixAlZvEtv3KA6qH/4\n4Ydi0qRJ8nNc7J09ezZTMGQ7DOxky5TM3LlzReXKleXnWCFTqVIl2V9TZfoFK2C++uor+TnSRYsX\nL5YdhIjshoGdbAmNMH755RdRsWJF+fnRo0dFuXLlxH/+8x8la9WREvruu++MsWHtOsrVEtkRAzvZ\nVnBwsOzEU6tWLfl5RESEaNasmWjSpEm0nZGsSL0g7VK0aFGxcuVKI/2yYsUKmR5yI6cWYaPIGNjJ\n1tKkSSNn7h999JHRLDosLEymQH744QfL1roj/YPgjfSL3oEJO0s3bNhgvItwo6FDh4qOHTtG2/KO\nnIOBnWwPDamxG/XPP/8UBQsWlF+7ePGieOutt+QqGrRgu3DhQqKPgxcJtNJD2gXHQdNr/fjok4qg\nXqhQIeFWaE/Yv39/mXJC/R6UdyBnYmAnxyhfvrzYsWNHpNk7cu+9e/eW/Tex8xNFxLA80t+UAmbj\na9euFYMGDRJ58+YVDRo0MJYy6rP07du3yyWN2IjkVnhRa9u2rbh37578HI2+3fx83Y6/OXIUrEbB\n7B3t7caNGye+//57cfPmTRmQ5syZI2+QJUsW2fACberQtBo/FxQUJNfGI5jv2bNHtuI7cOBAtC8C\nWJHTuXNn0bhxY08EOPy/xDsSQHs79IUl50qi8WqJp7Ru3VqmLbBUzw2uXr0qgztK5yKVkNh6NXjB\nQEDHBVM3wOqdhg0byhc0NAGJDt714JrFrVu3xBNPPCF7txYoUCDgYyXzuH8qQq6GVSrY0ITdoEjT\nbNmyRfZOxWwcs3I0mY7J008/LWf1+g1r5dGQ2Uswr8PsHEEdUIiNQd35GNjJNTXdkXbBDas6AM2l\n0W4P6RYsk9Q7HKF4V2hoqEzReB2Kl2H5pn49oUePHqqHRCZgYCfXQuoBxbl8Z+EhISFyoxMJ2amq\nZ8+e8j6uI0ydOtUT1xO8gKtiiDxIL2Z25coV+Xnfvn1F8eLFVQ+LTMLATuRB8+bNE/Pnz5f3Cxcu\nLAM7uQcDO5HHXLp0yWi6jWsTSMHEtGKGnImBnchjcIFUr7WDC8nY+EXuwsBO5CHo0Tpr1ix5Hztt\nBw4cqHpIZAEGdiKPuH79urEUFCZPniyLrJH7MLATeQTq3ehtBtGztXr16qqHRBZhYCfygHXr1sl6\nMJAzZ04xYsQI1UMiCzGwE3kA6srrJkyYIDJkyKB0PGQtBnYiDzh48KD8iNIK9evXVz0cshgDO5GL\n+TYAz5w5sxg9erTS8VBgMLATuRQqW44ZM8b4/Ntvv2XhM49gYCdyKVwgRa11qF27tmjevLnqIVGA\nMLATubh/qW7s2LGyfAB5AwM7kcv7l+pNRcg7GNiJXNy/FC3vyHsY2IlcBDl17DAF9C/VqziStzCw\nE7m4fyl2mZL3MLATuQT7l5KOgZ3IBdi/lHwxsBM5HPuX0uMY2Ikcjv1L6XEM7EQOxv6lFB0GdiIH\nY/9Sig4DO5FDsX8pxYSBnciB2L+UYsPATuRA7F9KsWFgJ3KYtWvXsn8pxYqBnchBbt++LWfoOvYv\npegwsBM5COq/hIeHy/vsX0oxYWAncojt27cbaRf2L6XYMLATOaR/aZs2bWQTDWD/UooNAzuRA2Cm\nvnPnTnm/bt267F9KsWJgJ3JQ/9J06dKJiRMnsn8pxYqBnchB/UuHDx/O/qUUJwZ2Iof0L61SpYrs\nkEQUFwZ2Iof0L0XZgKRJ+SdLceNZQuSQ/qUFChRQPSxyCAZ2Ihti/1JKDAZ2Ipth/1JKLAZ2Ihth\n/1IyAwM7kY2wfymZgYGdyCbYv5TMwsBOZBPsX0pmYWAnsgH2LyUzMbATKcb+pWQ2BnYixdi/lMzG\nwE6kEPuXkhUY2IkUYf9SsgoDO5Ei7F9KVmFgJ1KA/UvJSgzsRAHG/qVkNQZ2ogBj/1KyGgM7UQCx\nfykFAgM7UYCwfykFCgM7UYCwfykFCgM7UQCwfykFEs8sIouxfykFGgM7kcXYv5QCjYGdyELsX0oq\nMLATWYT9S0kVBnYii7B/KanCwE5kAfYvJZUY2IkswP6lpBIDO5HJ2L+UVGNgJzIR+5eSHXDdFbmy\nJou+GQhu3Lhh3L97964MvrrkyZPLnaBmYf9SsgPO2MlV9u/fL3uHBgcHG7cSJUoY/47a577/ljZt\nWjF48GBTjs3+pWQXDOzkKps2bTIuWvo7u1+8eHGij8v+pWQnDOzkKq+88oqscx4fLVq0SPRx2b+U\n7ISBnVwF/UPfe+89v78/V65cska6P1auXCmGDBkizp49G+nr7F9KdsPATq5cQ+7vrB27Qf3ZOHTt\n2jU5C//kk0/kLlIsZ0TJAPYvJTtiYCfPztrjM1s/dOiQzKPD5cuXRatWrUSdOnXkCwP7l5LdMLCT\nZ2ft/s7W4cSJE1G+tmzZMjFy5Eh5H6tr2L+U7IKBnTw5a4/PbD2mwO4re/bskdbOE6nEwE6enLXH\nZ7buT2BHqgbr5YcNGybz7kQqMbCT52bt8Z2t+xPY9V2t2HmKgl8RERHxenwiMzGwk+tn7ch/J2a2\n7m9g913+OHv27Hg9PpGZGNjJ9bP2li1bGp8jyMd3th7fwB4SEsINSqQUAzu53meffSaCgoLk/fbt\n28d7tv7o0SNx6tSpOL8Pj9uvXz+xd+9eERoamuDxEiUWqzuSK2HzEGbZu3fvltUchw4dKk6ePClK\nliwpW9Yh8BYtWlSkSJEizsc6d+5cnBdEGzZsKEaNGiXrrxOpxsBOrnDnzh2xfPlysXnzZrFt2zax\ndetWcfHixVh/BkG9WLFiolSpUvJWu3ZtkTt37ijfhxeEmBQsWFCWEKhVq5Ypz4PIDAzs5GhHjhyR\nG4PQUzS+K1Hu3bsnXwRwA2wuqlevnujSpYsM8kmT/pup3LdvX5SfxTLKL774Qrz77rt+zfqJAomB\nnRxpxYoV4ptvvhFLly6VaRdfyKfnz59fFCpUSDz77LMiY8aMMv+NQI1gjo1Ehw8flgEb9duvXr0q\nfw6P8+uvv8obUiqdOnUSnTt3Fhs3boz0+K1bt5apnRw5cgT0ORP5i4GdHAXplW7duomffvop0tcR\nvBs0aCCqVasmg7q/F0gRzM+cOSNn7fPnzxe7du2SX0fg79Wrlxg/frz49NNP5Qw9VapUYu7cuaJq\n1aqWPDciszCwk2Mg8GIW7dtIA7s9mzRpIlvQJSQlgvQLuh3h9uqrr8oZPC6u4p0A8vZHjx6VyyOR\nnvnyyy+jrIknsiMudyTbQ8ncN998U7z++utGUMdacTSKRm795ZdfNi3PjdQNSvOiqxKadugwc8eF\n1nXr1plyHCIrMbCTrV24cEGmV/TUC2bYKI37ww8/iOeff96y46ZPn15eHEUeP2vWrPJrmL2/9NJL\nYuHChZYdl8gMDOxkW+hUVLlyZblFX6/xgll6z549xRNPPBGQMVSqVEnm1VFrXa8H06hRI/Hjjz8G\n5PhECcHATrZ06dIlUbNmTZnzhmeeeUZMnz5dPPfccwEfS3BwsOjfv7/o2LGj/Bzdkt5++22xaNGi\ngI+FyB8M7GQ7WJKI/PaePXvk50WKFBGTJk0SmTJlUjYmpIBQjuD99983gnvTpk3F+vXrlY2JKCYM\n7GQ7gwcPFhs2bJD38+XLJ3d2+tvD1GotWrQwZu5Iy6DA2M2bN1UPiygSBnaylR07doghQ4bI+1ha\niKCOC5l20q5dO7kSR2+wgTLARHbCwE62SsFgV+eDBw/k5x9++KFsOWc3SMt89NFHsiQw4MVnzZo1\nqodFZGBgJ9vATF3f+YnVKKjbYld4F4FuSbp33nmHKRmyDQZ2soXLly+L4cOHGykYpDcwM7YzlBao\nU6eOUYJg2rRpqodEJDGwky3MmDFD3L5928hhZ8uWTTil9Z6+6xW7Ux8vSEakAgM7KYcORRMmTJD3\nUbzLSW3lsAQTu1EBa+5Xr16tekhEDOyk3qpVq8Q///wj76MOOjYEOQmKkOkwaydSjYGdlNNn648H\nSadAiz10UtIrUJ4+fVr1kMjjGNhJKezgREs7QHBEcwwrfPDBB6JixYrGUkoz4SLva6+9ZjyfP/74\nw/RjEMUHAzspFR4ebiwTRKNpK4+DrkjJklnTgsC30qTeao9IFQZ2Uso3CKIWuhWuX78uuyQVKFBA\nWCVPnjxGxUkGdlKNgZ2U2rp1q3HfijQM2tqhnjv88ssvonTp0vKGGi9mQp9VPc+OMsNY6UOkClvj\nkVL67Bb9RHPnzm364yOIX7lyRRYVQ+s7vQE1iouZDS9MO3fuFDdu3JCpH6vegRDFhYGdlDp27JiR\nysCs12xocH3gwAEZ2Lt27SqyZMkirBIaGhrpeTGwkypMxZBSaBitz9itcvDgQZExY0ZLg/rjz0Hf\nRUukAgM72SKwJ0+e3LJjYPOTlRdOdb4NtfXnRaQCAzsppS8/tOpiI1bDYFVMIAK773OwalklkT8Y\n2EkpfYkguhFZQS9VEIjA7vscrEwtEcWFgZ2U0uvCREREWJZfD1Rg930OTqt3Q+7CwE5KFS9eXH48\ndeqUuHbtmumPf/XqVfkxTZo0wmqo7qgrVqyY5ccjigkDOymFdebRBUazPPPMM/Jjv379ZLGxSZMm\nifPnzwsr7Nu3T37Mnz+/yJAhgyXHIPIHAzspVapUqSiB0UzocPT222/LYD59+nQxefJkI69v9jsD\nvOt4/DkRqcBL96SUb+EvKwJ70qRJRffu3eXNSr7vNhjYSTXO2EkpbBzS0yXYHerUhtArV6407pct\nW1bpWIgY2Ek5vSAXgvrSpUuF06A2jD5ulEaoVKmS6iGRxzGwk3Jt27Y1dp6GhYU5riH0r7/+auw0\n7dSpkyU1b4jig4GdlEPFxUaNGhnrzlEh0SnwIoQXI72kQJs2bVQPiYiBneyhS5cuxv1p06Y5Zta+\nevVqcfToUXm/WbNmImvWrKqHRMTATvaAvLS+mmT9+vVi8eLFwu6wxHHYsGHG5++9957S8RDpGNjJ\nFtAQ+rvvvjPy06NGjRIXLlwQdjZy5EijjEC3bt0ibbYiUomBnWy1pr1v377yPioyDhkyxLYpmf/9\n73/GShg0yR46dKjqIREZGNjJVrD1X6+zsnbtWvHTTz8Ju8EOU7zo+F4TCEQtGiJ/MbCTrWBlyYwZ\nM4x65kjJYDmhXVy8eFFe6PVNwVSpUkX1sIgiYWAnW6ZkpkyZYnw+YMAAWwT3s2fPio4dOxo1YWrU\nqCFGjBihelhEUTCwky21atVKfPvtt0Znoi+++ELMmTNHWc798OHDol27dkbz7XLlyokFCxaIlClT\nKhkPUWwY2Mm2sHxw9OjRxudff/21LOaFmXOgPHz4UHz//feiRYsWxnGxNHPZsmUibdq0ARsHUXww\nsJOtIYc9c+ZMo+QA1rhjIxBmy1bP3rHxCLN0vHO4d++e/Fq9evVkUE+fPr2lxyZKDAZ2ckSRsK1b\nt4rnn3/eKBY2aNAg0bVrV7Fp0ybTAzxm5uPGjRPNmzcXu3fvll9LnTq1GDNmjFi0aJG8T2RnrMdO\njmmhhyCOnZ4DBw4U9+/fF5s3b5a3kJAQ0bhxY/Hqq6+KdOnSJejxkcfHY/38889ymSU+12HVy9Sp\nU0W+fPlMfEZE1kmi2XUHCFmidevWckenE7bsx2TXrl1ydcrGjRsjfR0XMsuUKSMKFSpk3GKq3YJq\njOHh4bK5B5pkbN++3VjtosuUKZPo37+/XN6Ihh1OtHDhQtGwYUP5fHmh1zs4YydHzt7RlGPLli1i\n/PjxchMTAtfdu3fFunXr5E2XOXNm2cwDQQ3BGbny27dvyyCOC6PRQaMMBPOmTZuKVKlSBfCZEZmD\ngZ0cC7Nz9DFFzRZsasLqlT179kQK2NhIpG8miqt0cN26dUXnzp1Z84Ucj4GdHA+z8g8++EDebt26\nJeu5b9u2Td5w8RN1ZzBLR8BHI2vMwkNDQ2U1Sf2WM2dO1U+DyDQM7B6sohhTCsINsGKlQoUK8kb/\nrsMHp14joIThb9tjsB4cK0rIG/TftV57h7yBgd1jGNi9BReL8TvHOzXyDgZ2j0F+GXlo8gZcW8B1\nBfIWBnaPyZYtmzh//rzqYVCA4HedPXt21cOgAGNg9xgs6zt37pxtOxORML08An7n5C0M7B6DP3Lk\nXa9cuaJ6KBQADOzexMDuMfofeSBL31L8VK1aVd7MwMDuTQzsHqPnWxnYvQFpN+bYvYeB3WOyZMki\n66ag1ji5G5a1njhxQjz11FOqh0IBxsDuMUFBQbKnKErUknooYPbss8/KF9siRYqI+fPnm/bYKKeA\nwmioqUPewu1oHoR+natXr1Y9DM9buXKlbOaBrkyjRo2S5ZTR+g8z7YIFCyb68VHWGPXp8cJB3sLA\n7tHAjm5A6ESUJk0a1cPxrM8//1wGXdRM12u54HPUuTEjsKMxCWbreJdG3sJUjAeVL19eFodC9UNS\nA///UU8enZ98C3Thd5MnTx5TjoEZOx6PvIeB3YNy584tnnzySbFq1SrVQ/GsixcvypRLdCtWzFjF\ngoum6BDFKpfexMDuQSgI9eabb8rmFG4u4Wv31UkozoXliI+L7mvxhd8t6tTXrFkz0Y9FzsPA7lFt\n27YVx48f56xdEeS9kf8OCwuL1DgbefHELkXF46H5dsuWLdnn1KMY2D2qcOHComLFimLKlCmqh+JZ\naJSNRtpoNo3m4phlo89qYneK4sX62LFj8sWbvImB3cPatWsnFixYIIMABd5LL70k5syZIw4cOCBe\nf/11MWLECPHNN98kekXM6NGj5UVTrIsnb0qiscyfZ925c0c899xzsv/nkiVL2IzBBebNmydX2vz2\n22+idu3aqodDijCwe9yff/4pXnzxRZkGQE6WnOvSpUsyxVanTh0xffp01cMhhRjYSXTr1k2mBLAF\nPVeuXKqHQwmAP2O8MGM36969e0XGjBlVD4kUYmAncf36dblCA+uqly9fLvLly6d6SBQPWLLatWtX\nMXnyZHnN5NVXX1U9JFKMF09J1hNZs2aNyJAhg3jhhRfEzp07VQ+J/IQiX/qeBOTXGdQJGNjJ6IX6\nxx9/yFolyLkPHTpU3LhxQ/WwKAZ4o71ixQpRqVIleaF06dKlctkkETCwkyE4OFgGiffee08MGTJE\n5M2bV3z11Vdsfm2zGTry6OiwVKtWLbm7FBfAq1WrpnpoZCPMsVOMtUyGDx8uxo4dK27fvi2DPNZG\nly1bVtaZQfoGW+LJOthBigqcV69elekxFPXavn277FmLwD5w4EA5Yyd6HAM7EZHLMBVDROQyDOxE\nRC7DwE5E5DIM7ERELsPATkTkMgzsREQuw8BOROQyDOxERC7DwE5EJNzl/wHgIM4lXK3bIAAAAABJ\nRU5ErkJggg==\n"
        }
      },
      "cell_type": "markdown",
      "id": "a178b540-dcec-4a13-a20a-042b0b79abd0",
      "metadata": {},
      "source": [
        "Try changing the data to `t = jnp.array([-20, -10, 30])`. Notice how\n",
        "changing the data for one group does not affect the expectation of the\n",
        "posterior of the other groups.\n",
        "\n",
        "Previously, the `complete_pooling` model returned an array of shape\n",
        "`(200,)` with 200 being the size of the sample space of `Theta`. With\n",
        "the addition of the conditional statement, what array shape does\n",
        "`no_pooling` return? Make sure you understand the shape change, what\n",
        "each dimension reflects, and why the values sum in the ways that they\n",
        "do.\n",
        "\n",
        "The key theoretical distinction is that in *no pooling*, we’re asserting\n",
        "that the timing behavior of professors in different departments is not\n",
        "informative about the behavior of professors from other departments.\n",
        "This is a unrealistically strong assumption. This highlights an\n",
        "important practical tradeoff: *no pooling* can better capture\n",
        "differences between departments but may suffer from high variance in its\n",
        "estimates due to small sample sizes within each department. Complete\n",
        "pooling has lower variance but may miss important systematic differences\n",
        "between departments.\n",
        "\n",
        "## Partial pooling\n",
        "\n",
        "In partial pooling, every observation updates beliefs across multiple\n",
        "levels of abstraction. The critical insight is that while departments\n",
        "may differ in their typical timing behavior, the departments are not\n",
        "statistically independent. They are influenced by common causes. To\n",
        "model the president’s reasoning about these multiple levels of\n",
        "uncertainty, we modify the `no_pooling` model by adding hyperpriors over\n",
        "the hyperparameters $\\mu$ and $\\tau$.\n",
        "\n",
        "Rather than being fixed values, $\\mu$ and $\\tau$ become latent\n",
        "parameters to be jointly inferred alongside with $\\theta_d$.\n",
        "\n",
        "The higher-level latent parameters $\\mu$ and $\\tau$ characterize a\n",
        "population-level distribution from which department-specific $\\theta_d$\n",
        "values are generated.",
        "The key change to the model is that rather than\n",
        "${ \\theta_d ~\\sim~ \\mathcal{N}\\left( \\dot{\\mu}{=}0, \\dot{\\tau}{=}20 \\right) }$,\n",
        "we have ${ \\theta_d ~\\sim~ \\mathcal{N}\\left(\\mu, \\tau \\right) }$ and\n",
        "specify hyperpriors on these new latent parameters:\n",
        "${ \\mu ~\\sim~ \\mathcal{N}(0, \\dot\\sigma_{\\mu}) }$ and\n",
        "${ \\tau ~\\sim~ \\text{HalfCauchy}(\\dot\\sigma_{\\tau}) }$.\n",
        "\n",
        "Note that there are new fixed parameters: `mu_scale`\n",
        "(${ \\dot\\sigma_\\mu}$) and `tau_scale` (${ \\dot\\sigma_\\tau }$), but I am\n",
        "omitting there from the graphical schematic for the sake of visual\n",
        "simplicity.\n",
        "\n",
        "$$\n",
        "\\begin{align*}\n",
        "\\mu ~\\sim&~ \\mathcal{N}(0, \\dot\\sigma_{\\mu}) \\\\\n",
        "\\tau ~\\sim&~ \\text{HalfCauchy}(\\dot\\sigma_{\\tau}) \\\\\n",
        "\\theta_d ~\\sim&~ \\mathcal{N}\\left(\\mu, \\tau \\right) \\\\\n",
        "t_d ~\\sim&~ \\mathcal{N}(\\theta_d, \\dot\\sigma{=}15) \\\\\n",
        "d ~\\in&~ D,~~~ \\text{where}~~~ D = \\{\\text{G}, \\text{E}, \\text{M} \\}\n",
        "\\end{align*}\n",
        "$$\n",
        "\n",
        "![](attachment:generated/multilevel-models/partial-pooling.png)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "04b64c86",
      "metadata": {},
      "outputs": [],
      "source": [
        "%reset -f\n",
        "import jax\n",
        "import jax.numpy as jnp\n",
        "from memo import memo\n",
        "from enum import IntEnum\n",
        "from jax.scipy.stats.norm import pdf as normpdf\n",
        "from jax.scipy.stats.cauchy import pdf as cauchypdf\n",
        "from matplotlib import pyplot as plt\n",
        "\n",
        "normpdfjit = jax.jit(normpdf)\n",
        "\n",
        "class Department(IntEnum):\n",
        "    GOVERNMENT = 0\n",
        "    ENGLISH = 1\n",
        "    MATH = 2\n",
        "\n",
        "t = jnp.array([1, -10, 30])\n",
        "sigma = 15\n",
        "\n",
        "Mu = jnp.linspace(-25, 25, 100)  ### sample space for new hyperprior\n",
        "Tau = jnp.linspace(1, 30, 100)  ### sample space for new hyperprior\n",
        "Theta = jnp.linspace(-40, 40, 200)\n",
        "\n",
        "### PDF for new hyperprior\n",
        "@jax.jit\n",
        "def half_cauchy(x, scale=1.0):\n",
        "    return 2 * cauchypdf(x, 0, scale)\n",
        "\n",
        "@jax.jit\n",
        "def professor_arrival_likelihood(d, theta):\n",
        "    ### likelihood of a professor from department d \n",
        "    ### showing up t_d minutes early/late, \n",
        "    ### for a given theta and sigma.\n",
        "    return normpdf(t[d], theta, sigma)\n",
        "\n",
        "@memo\n",
        "def department_model[_mu: Mu, _tau: Tau](d):\n",
        "    department: knows(_mu, _tau)\n",
        "    department: chooses(theta in Theta, wpp=normpdfjit(theta, _mu, _tau))\n",
        "    return E[ department[professor_arrival_likelihood(d, theta)] ]\n",
        "\n",
        "@memo\n",
        "def partial_pooling[_mu: Mu, _tau: Tau](mu_scale=5, tau_scale=5):\n",
        "    president: knows(_mu, _tau)\n",
        "    president: thinks[\n",
        "        population: chooses(mu in Mu, wpp=normpdfjit(mu, 0, mu_scale)),\n",
        "        population: chooses(tau in Tau, wpp=half_cauchy(tau, tau_scale)),\n",
        "    ]\n",
        "    president: observes_event(wpp=department_model[population.mu, population.tau]({Department.GOVERNMENT}))\n",
        "    president: observes_event(wpp=department_model[population.mu, population.tau]({Department.ENGLISH}))\n",
        "    president: observes_event(wpp=department_model[population.mu, population.tau]({Department.MATH}))\n",
        "    return president[Pr[population.mu == _mu, population.tau == _tau]]"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "705255c7-60af-4768-87da-0e17ee202d2e",
      "metadata": {},
      "source": [
        "Let’s think about what each variable represents.\n",
        "\n",
        "The likelihood model `professor_arrival_likelihood` gives the\n",
        "probability that a professor from department $d$ shows up $t_d$ minutes\n",
        "early/late, under a given $\\theta_d$ and $\\sigma$:\n",
        "\n",
        "$$\n",
        "t_d ~\\sim~ \\mathcal{N}(\\theta_d, \\dot\\sigma{=}15)\n",
        "$$\n",
        "\n",
        "We can write the likelihood as\n",
        "${ P(t_d \\mid \\theta_d, \\mu, \\tau ; \\sigma) }$. The marginal posterior\n",
        "${ P(\\theta_d \\mid t) }$ will thus represent the president’s belief\n",
        "about the true value of $\\theta_d$ for a given department $d$."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "6fb148d3",
      "metadata": {},
      "outputs": [],
      "source": [
        "@memo\n",
        "def department_model_theta[_theta: Theta](d, mu_scale, tau_scale):\n",
        "    obs: thinks[\n",
        "        department: chooses(mu in Mu, tau in Tau, wpp=partial_pooling[mu, tau](mu_scale, tau_scale)),\n",
        "        department: chooses(theta in Theta, wpp=normpdfjit(theta, mu, tau))\n",
        "    ]\n",
        "    obs: observes_event(wpp=professor_arrival_likelihood(d, department.theta))\n",
        "    obs: knows(_theta)\n",
        "    return obs[Pr[department.theta == _theta]]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "0b1f0d8b",
      "metadata": {},
      "outputs": [],
      "source": [
        "def plot_model(mu_scale=1, tau_scale=1, figsize=(10, 8)):\n",
        "    posterior = partial_pooling(mu_scale=mu_scale, tau_scale=tau_scale)\n",
        "\n",
        "    # Marginal over Tau (sum over Mu)\n",
        "    posterior_tau = posterior.sum(axis=0)\n",
        "    # Marginal over Mu (sum over Tau)\n",
        "    posterior_mu = posterior.sum(axis=1)\n",
        "\n",
        "    fig, axs = plt.subplots(3, 1, figsize=figsize)\n",
        "\n",
        "    ax = axs[0]\n",
        "    ax.axvline(0, color=\"black\", linestyle=\"-\")\n",
        "    ax.plot(Mu, posterior_mu, label=r\"$P(\\mu \\mid t)$\")\n",
        "    mu_expectation = jnp.dot(Mu, posterior_mu)\n",
        "    ax.axvline(\n",
        "        mu_expectation, \n",
        "        color='red', \n",
        "        linestyle='--', \n",
        "        label=r\"$\\operatorname{E}\" + rf\"[\\mu \\mid t]={mu_expectation:0.2f}$\")\n",
        "    _ = ax.set_title(r\"Posterior of $\\mu$\")\n",
        "    _ = ax.legend(bbox_to_anchor=(0.9, 0.5), loc='center left')\n",
        "\n",
        "    ax = axs[1]\n",
        "    ax.axvline(0, color=\"black\", linestyle=\"-\")\n",
        "    ax.plot(Tau, posterior_tau, label=r\"$P(\\tau \\mid t)$\")\n",
        "    tau_expectation = jnp.dot(Tau, posterior_tau)\n",
        "    ax.axvline(\n",
        "        tau_expectation, \n",
        "        color='red', \n",
        "        linestyle='--', \n",
        "        label=r\"$\\operatorname{E}\" + rf\"[\\tau \\mid t]={tau_expectation:0.2f}$\")\n",
        "    _ = ax.set_title(r\"Posterior of $\\tau$\")\n",
        "    _ = ax.legend(bbox_to_anchor=(0.9, 0.5), loc='center left')\n",
        "\n",
        "    ax = axs[2]\n",
        "    ax.axvline(0, color=\"black\", linestyle=\"-\")\n",
        "    for d in Department:\n",
        "        department_name = Department._member_names_[d]\n",
        "        department_abbrev = department_name[0]\n",
        "        theta_posterior = department_model_theta(d, mu_scale, tau_scale)\n",
        "        theta_expectation = jnp.dot(Theta, theta_posterior)\n",
        "        ax.plot(\n",
        "            Theta, \n",
        "            theta_posterior, \n",
        "            label=(\n",
        "                rf\"$P(\\theta_{department_abbrev} \\mid t),~ \" \n",
        "                + r\"\\operatorname{E}\" \n",
        "                + rf\"[\\theta_{department_abbrev} \\mid t]={theta_expectation:0.2f}$\"))\n",
        "    _ = ax.set_title(r\"Posterior of $\\theta_d$\")\n",
        "    _ = ax.legend(bbox_to_anchor=(0.9, 0.5), loc='center left')\n",
        "\n",
        "    _ = plt.suptitle(f\"mu_scale = {mu_scale}, tau_scale = {tau_scale}\", y=1)\n",
        "    plt.tight_layout()\n",
        "    plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "8bdd424a",
      "metadata": {},
      "outputs": [],
      "source": [
        "import ipywidgets as widgets\n",
        "from ipywidgets import interactive\n",
        "\n",
        "def plot_model_widget(mu_scale=1, tau_scale=1):\n",
        "    plot_model(mu_scale=mu_scale, tau_scale=tau_scale, figsize=(9, 7))\n",
        "\n",
        "interactive_plot = interactive(\n",
        "    plot_model_widget, \n",
        "    mu_scale=widgets.IntSlider(min=1, max=40, step=1, value=1),\n",
        "    tau_scale=widgets.IntSlider(min=1, max=40, step=1, value=1),\n",
        ")\n",
        "output = interactive_plot.children[-1]\n",
        "output.layout.height = '700px'\n",
        "interactive_plot"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "6e8b3d9d-80d3-45ab-8ccf-ea6223199749",
      "metadata": {},
      "source": [
        "### Multiple observations with missing data\n",
        "\n",
        "What if the president instead had a large meeting with 7 professors and\n",
        "the departments were not equally represented? How would she update her\n",
        "beliefs given unequal observations? Let’s say that she invited 3\n",
        "professors from the Department of Government, 3 from the Department of\n",
        "English, and 1 from the Department of Math. We can explore how\n",
        "modulating the higher-level hyperpriors changes her posterior inference\n",
        "about each department.\n",
        "\n",
        "$$\n",
        "\\begin{align*}\n",
        "\\mu ~\\sim&~ \\mathcal{N}(0, \\dot\\sigma_{\\mu}) \\\\\n",
        "\\tau ~\\sim&~ \\text{HalfCauchy}(\\dot\\sigma_{\\tau}) \\\\\n",
        "\\theta_d ~\\sim&~ \\mathcal{N}\\left(\\mu, \\tau \\right) \\\\\n",
        "t_{(d,i)} ~\\sim&~ \\mathcal{N}(\\theta_d, \\dot\\sigma{=}15) \\\\\n",
        "d ~\\in&~ D,~~~ \\text{where}~~~ D = \\{\\text{G}, \\text{E}, \\text{M} \\}\n",
        "\\end{align*}\n",
        "$$"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "a992789b",
      "metadata": {},
      "outputs": [],
      "source": [
        "%reset -f\n",
        "import jax\n",
        "import jax.numpy as jnp\n",
        "from memo import memo\n",
        "from enum import IntEnum\n",
        "from jax.scipy.stats.norm import pdf as normpdf\n",
        "from jax.scipy.stats.norm import logpdf as normlogpdf\n",
        "from jax.scipy.stats.cauchy import pdf as cauchypdf\n",
        "from matplotlib import pyplot as plt\n",
        "\n",
        "normpdfjit = jax.jit(normpdf)\n",
        "\n",
        "class Department(IntEnum):\n",
        "    GOVERNMENT = 0\n",
        "    ENGLISH = 1\n",
        "    MATH = 2\n",
        "\n",
        "###NEW\n",
        "t = jnp.array([\n",
        "    [10, 0, -10], \n",
        "    [-16, -15, -14], \n",
        "    [30, jnp.nan, jnp.nan],\n",
        "])\n",
        "\n",
        "sigma = 15\n",
        "\n",
        "Mu = jnp.linspace(-25, 25, 100)\n",
        "Tau = jnp.linspace(1, 30, 100)\n",
        "Theta = jnp.linspace(-40, 40, 200)\n",
        "\n",
        "@jax.jit\n",
        "def half_cauchy(x, scale=1.0):\n",
        "    return 2 * cauchypdf(x, 0, scale)\n",
        "\n",
        "@jax.jit\n",
        "def professor_arrival_likelihood(d, theta):\n",
        "    ### likelihood of a professor from department d \n",
        "    ### showing up t_d minutes early/late, \n",
        "    ### for a given theta and sigma.\n",
        "    return jnp.exp(jnp.nansum(normlogpdf(t[d], theta, sigma)))  ###NEW\n",
        "\n",
        "@memo\n",
        "def department_model[_mu: Mu, _tau: Tau](d):\n",
        "    department: knows(_mu, _tau)\n",
        "    department: chooses(theta in Theta, wpp=normpdfjit(theta, _mu, _tau))\n",
        "    return E[ department[professor_arrival_likelihood(d, theta)] ]\n",
        "\n",
        "@memo\n",
        "def partial_pooling[_mu: Mu, _tau: Tau](mu_scale=5, tau_scale=5):\n",
        "    president: knows(_mu, _tau)\n",
        "    president: thinks[\n",
        "        population: chooses(mu in Mu, wpp=normpdfjit(mu, 0, mu_scale)),\n",
        "        population: chooses(tau in Tau, wpp=half_cauchy(tau, tau_scale)),\n",
        "    ]\n",
        "    president: observes_event(wpp=department_model[population.mu, population.tau]({Department.GOVERNMENT}))\n",
        "    president: observes_event(wpp=department_model[population.mu, population.tau]({Department.ENGLISH}))\n",
        "    president: observes_event(wpp=department_model[population.mu, population.tau]({Department.MATH}))\n",
        "    return president[Pr[population.mu == _mu, population.tau == _tau]]\n",
        "\n",
        "@memo\n",
        "def department_model_theta[_theta: Theta](d, mu_scale, tau_scale):\n",
        "    obs: thinks[\n",
        "        department: chooses(mu in Mu, tau in Tau, wpp=partial_pooling[mu, tau](mu_scale, tau_scale)),\n",
        "        department: chooses(theta in Theta, wpp=normpdfjit(theta, mu, tau))\n",
        "    ]\n",
        "    obs: observes_event(wpp=professor_arrival_likelihood(d, department.theta))\n",
        "    obs: knows(_theta)\n",
        "    return obs[Pr[department.theta == _theta]]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "16da25fb",
      "metadata": {},
      "outputs": [],
      "source": [
        "def plot_model(mu_scale=1, tau_scale=1, figsize=(10, 8)):\n",
        "    posterior = partial_pooling(mu_scale=mu_scale, tau_scale=tau_scale)\n",
        "\n",
        "    # Marginal over Tau (sum over Mu)\n",
        "    posterior_tau = posterior.sum(axis=0)\n",
        "    # Marginal over Mu (sum over Tau)\n",
        "    posterior_mu = posterior.sum(axis=1)\n",
        "\n",
        "    fig, axs = plt.subplots(3, 1, figsize=figsize)\n",
        "\n",
        "    ax = axs[0]\n",
        "    ax.axvline(0, color=\"black\", linestyle=\"-\")\n",
        "    ax.plot(Mu, posterior_mu, label=r\"$P(\\mu \\mid t)$\")\n",
        "    mu_expectation = jnp.dot(Mu, posterior_mu)\n",
        "    ax.axvline(\n",
        "        mu_expectation, \n",
        "        color='red', \n",
        "        linestyle='--', \n",
        "        label=r\"$\\operatorname{E}\" + rf\"[\\mu \\mid t]={mu_expectation:0.2f}$\")\n",
        "    _ = ax.set_title(r\"Posterior of $\\mu$\")\n",
        "    _ = ax.legend(bbox_to_anchor=(0.9, 0.5), loc='center left')\n",
        "\n",
        "    ax = axs[1]\n",
        "    ax.axvline(0, color=\"black\", linestyle=\"-\")\n",
        "    ax.plot(Tau, posterior_tau, label=r\"$P(\\tau \\mid t)$\")\n",
        "    tau_expectation = jnp.dot(Tau, posterior_tau)\n",
        "    ax.axvline(\n",
        "        tau_expectation, \n",
        "        color='red', \n",
        "        linestyle='--', \n",
        "        label=r\"$\\operatorname{E}\" + rf\"[\\tau \\mid t]={tau_expectation:0.2f}$\")\n",
        "    _ = ax.set_title(r\"Posterior of $\\tau$\")\n",
        "    _ = ax.legend(bbox_to_anchor=(0.9, 0.5), loc='center left')\n",
        "\n",
        "    ax = axs[2]\n",
        "    ax.axvline(0, color=\"black\", linestyle=\"-\")\n",
        "    for d in Department:\n",
        "        department_name = Department._member_names_[d]\n",
        "        department_abbrev = department_name[0]\n",
        "        theta_posterior = department_model_theta(d, mu_scale, tau_scale)\n",
        "        theta_expectation = jnp.dot(Theta, theta_posterior)\n",
        "        ax.plot(\n",
        "            Theta, \n",
        "            theta_posterior, \n",
        "            label=(\n",
        "                rf\"$P(\\theta_{department_abbrev} \\mid t),~ \" \n",
        "                + r\"\\operatorname{E}\" \n",
        "                + rf\"[\\theta_{department_abbrev} \\mid t]={theta_expectation:0.2f}$\"))\n",
        "    _ = ax.set_title(r\"Posterior of $\\theta_d$\")\n",
        "    _ = ax.legend(bbox_to_anchor=(0.9, 0.5), loc='center left')\n",
        "\n",
        "    _ = plt.suptitle(f\"mu_scale = {mu_scale}, tau_scale = {tau_scale}\", y=1)\n",
        "    plt.tight_layout()\n",
        "    plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "74a1604f",
      "metadata": {},
      "outputs": [],
      "source": [
        "import ipywidgets as widgets\n",
        "from ipywidgets import interactive\n",
        "\n",
        "def plot_model_widget(mu_scale=1, tau_scale=1):\n",
        "    plot_model(mu_scale=mu_scale, tau_scale=tau_scale, figsize=(9, 7))\n",
        "\n",
        "interactive_plot = interactive(\n",
        "    plot_model_widget, \n",
        "    mu_scale=widgets.IntSlider(min=1, max=40, step=1, value=1),\n",
        "    tau_scale=widgets.IntSlider(min=1, max=40, step=1, value=1),\n",
        ")\n",
        "output = interactive_plot.children[-1]\n",
        "output.layout.height = '700px'\n",
        "interactive_plot"
      ]
    },
    {
      "attachments": {
        "generated/multilevel-models/partial-pooling-learnedsigma.png": {
          "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAF2CAYAAAB6XrNlAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90\nbGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAO\nxAAADsQBlSsOGwAAP69JREFUeJzt3QmcTeX/B/CHYWRs2SdqmBDGVnalkCxRiEiSFMlWlvSrZPiV\nklC/RHZCo/hlQrKFkkSSZF9S9t3YDTOW8399nl/n/O81Y8xyzn2ec+7n/Xrd3An3PPc693uf+z3P\n8/1mMgzDEERE5BmZVQ+AiIjsxcBOROQxDOxERB7DwE5E5DEM7EREHsPATkTkMQzsREQew8BOROQx\nWVQPINhgP9jhw4fFb7/9Jv744w8RFxcnLl26JK5cuSKyZcsmsmfPLu644w5RpUoVUblyZZEvXz7V\nQyZKAufr1q1bxfr168WOHTvExYsXxeXLl+Xv3XbbbSJHjhyiTJky8jwuV66cyJo1q+ohBxUG9gA4\ncOCAmD59uli9erV8Ixw7dizVfzcyMlK+OerVqyeeeeYZkSdPHkfHSnSzCcn3338vYmNj5aRk06ZN\nIiEhIVV/FxOWihUriqpVq4pWrVqJhx9+WGTKlMnxMQezTCwp4Izr16/LN8KYMWPEvHnz5M8ZFRYW\nJtq3by+6d+8uKlWqZMs4iVJy5swZMW3aNDF27Fixc+dOWx6zdOnS8hzu0KGDuP322215TPLHwG4z\nBPApU6aI4cOHi127dvn9XkhIiPxaihk4Zi/4NSIiQn51zZIli0hMTBTx8fHizz//lDN787Z79+4k\nx7n//vvFwIEDRaNGjQL47ChYHD16VAwaNEjExMTIc9IXgjHOXfOGSUbevHnleQxIyZw+fVps3LjR\n7zzGh0RyE5W3335bhIeHB/T5eR4CO9nj77//NurWrYsPSr9blSpVjClTphgXLlxI1+MeO3bM+OCD\nD4zixYsneewXXnjBOH36tO3PhYLT9evXjRkzZhj58uXzO8+yZctmPPfcc8Yvv/wi/0x6HnfNmjVG\nhw4d5GP5PjaO9cUXX6TrcSl5DOw2uHbtmjF69GgjR44c1skaGhpqdOzY0fj1119tO87Vq1eNBQsW\nGE2aNPF7YxQtWtRYuHChbceh4HTkyBGjRYsWfucWJhPDhw83Tp48adtxTpw4YQwbNizJRAXHxhgo\n4xjYMyguLs54+OGH/U7Q6tWrG9u2bXP0uEuXLjUiIiL8jvvSSy8ZiYmJjh6XvAkTA99ZeubMmY03\n3njDuHTpkmPHxGPjGDiW7+x90aJFjh0zWDCwZ8ChQ4eMcuXK+X1dRcrkypUrATn+2bNnjS5duvgF\n92bNmjn6ZiTv+fLLL40sWbJY51BUVJSt3zRvZe3atUbZsmWt42MsM2fODNjxvYiBPZ0OHjxolChR\nwjoZy5cv7/gs/Wa+++47o0CBAtZYHnnkEQZ3SpWpU6camTJlss6dPn36KDl3cEwc2xwHxjRt2rSA\nj8MrGNjTmSPErMY8CWvWrClTMirhQwW5dnNMzZs3D9g3B3Kn2NhYvzTIiBEjVA9J5vN900EYI6Ud\nA3s6LmA+8MAD1slXq1Yt4/z584Yuq3J8g3vPnj1VD4k0tW7dOnmB3zxXRo0aZegCY/FdhICxUtow\nsKcRruabJ13FihWNU6dOGTrBzD1//vzWGJcvX656SKQZpD18v3G+++67hm4wJnN8uI51+fJl1UNy\nFQb2NNi+fbu1BhdLGzFD1hGWRJpvimLFihnnzp1TPSTSyJtvvum3xFDH9eMYE9KJ5jj79++vekiu\nwsCehhQMcunmiYZ16zrDGnpzrN26dVM9HNIEVruYeXUsLdR53fjhw4eNvHnzyrGGhIQwJZMGDOyp\n9OGHH1qBErtLsSlJZ9iNWqRIEWvMK1asUD0kUgwX032X52KHqe5iYmL8UjJcEJA6DOypkJCQYBQu\nXFj7FExKKZlGjRqpHg4pNnv2bO1TMLdKyXCVTOqw0UYqzJkzxyq127VrV1lK1w2aNGkiHnroIXl/\nyZIlyRYTo+CBSqOmoUOHuqJ0LsaIsSb3HOjmGNhTwfdkQmB3E5RHNY0bN07pWEid7du3yzLSUL9+\nfVk61y3QsAM13GH58uWysQeljIH9FrZs2SJWrlwp76NEbsmSJYWbPPHEE6Jw4cLyPsoJo1sTBR/f\nD3XfD3u34AQlbRjYb2HSpEmufkOEhoaKF198Ud5HjWx0wKHggk5HaJYBRYoUEc2aNRNugzFj7DB1\n6tRUd28KVgzst2DO1vPnzy+aNm1q62N369ZN5hBPnjyZ5Pfq1q0rChQoYMtxOnbsmOT5UPBAG7uz\nZ8/K+2hsgaYuboOeqRg74Lls3rxZ9ZC0xsCeAswKkIqBatWqyQ5IdkKHmaJFiyYbwPFmtKv93d13\n3y0/mACdbCi4+P6b16xZU7hVjRo1rPs8j1PGwJ4CzArQjR3QAsxOWGqKx08ueO/bt0+mTewK7PhW\nYI4fx+TX2ODiGwTtPI/R1Brn1q1uf/zxhy3H8x07A3vK3PedzANvCPjrr7/EhQsXxL333pvsTB7Q\n2d0uGP93330nP6jwLcTu50P6n8f4ZnjXXXfZmh5BX1TfZcH4pvn6669b/U8R2CtUqGDL8dAfGN88\n4+LiGNhvgYFdUWA3g3dys/KUfi+90Dzb93kxsAcHNEg389H4N7dz7TrOT99zdOnSpbLRte+6czth\n7DiPsScDzwnPDYsDKCmmYlJgbkrCxSY7Zzq3Ct746opjRkVF2XY85NlvfF7kfUjpXb16Vd4vUaKE\no8faunWrKFeunKPHMM9jfPM8c+aMo8dyMwb2FFy+fFn+mj17dtt36SGw43FLlSqV5Pd+/fVXuYEk\nW7Zsth0Px7rxeZH3+f5b+54Ddtu/f79crVK+fHnhJJ7HqcPAngJ81QMnvu4hsIeHh4vMmTMnuSB1\n8OBBW9Mw4PshwYunwXcOg5NpCzPd43Rg930OPI9vjoE9FSeR75vDDvgKiZUvhw4dkjcTZjzmJig7\n0zA3Pgc7vwmQ3nwDod3nsS9zWbDTgZ3ncerw4mkKzCv72IaP5Yl2pWOwcgDy5Mkj18e3aNFCxMfH\ni8WLF4uwsDBrhQHyiU8//bQtx/QtJWA+L/I+339rJ8tJ7Nq1S/5atmxZ4STf9AvP45vjjD0FhQoV\nkr/i4tPhw4dtv3A6ceJEUb16dblFGksRO3XqJH/FpqXjx4/LVI1d9uzZY90vWLCgbY9LesMqFXNj\n3d69ex07jrmzNVeuXMJJ5nmMxQWYGFHyOGNPwX333ee3RBAB187AXqdOHdG8efMkv48cu5NLNytX\nrmz745OekK7AShV8S8Q5YOc3T19mcTx8+8Tu1n79+oncuXPbegyM3TyP8ZyYirk5zthTufYbFzXt\ngsCO5ZOYTQWK+YbATMfOjU/knvMYy1x9r+nY6dVXXxWNGzcWP/30k/jwww9Fzpw5bT8Gxo5vssB9\nGCljYE8BAqBZMMmunW7Xrl2TF5oCGVwx0zE/mHBxi7nJ4BKIrfhI7y1atEjm8S9evJhktZcdfCdX\nvpMuSoqBPQUIgOaGi3Xr1onr16/bcpEJF4ACGdixAufEiRPyPmc6wcf333zt2rXCrbC/w8TzOGUM\n7LdQu3Zt+SsCI7YyZxRWDWAGPWTIEBEo06dPT/J8KHhgT4SZGpkxY4b81ug2GHNMTIx1gZbpxJQx\nsN9C586dXd1vEVuvx48fL+/jYlbr1q1VD4kUfPPs0KGDtUN04cKFwm0WLFggDhw4IO/juTCdmDIG\n9ltA9cVatWpZJ5eTS8acMH/+fGupJhpu5MiRQ/WQSAE0dXHzBMV3zL7PhZLHwJ4K5m5QpFDM2a9b\n8A1B5kXzBx98UN7HRjiUjXaL3bt3W2nQhx56yPFCY17AwJ4KTz75pNXlCIHS/EqoO3R0xw3Q5R3d\n3il4+fbsfeutt4RbDBgwwNV9h5UwKFWGDBli4OXCrVGjRsb169cNnZ07d84oVqyYNeYlS5aoHhIp\nlpCQYJQqVco6J2JjYw3dzZ492xovxo7nQLfGwJ5KiYmJRuXKla2TbNKkSYbOunbtao31+eefVz0c\n0sSqVauMTJkyyfOiUKFCxokTJwxdHT9+3ChYsKAcK8aMsVPqMLCnwcaNG42sWbPKEy137tzG/v37\nDR0tW7bMCupFixY1Tp8+rXpIpJE+ffpY58dTTz1l6KpNmzbWOPv27at6OK7CwJ5G77zzjnWyVatW\nTaY8dLJ7924jPDzcGuPChQtVD4k0c/HiRb+UzH/+8x9DNx999JFfCgZjptRjYE9HSqZq1arWSVev\nXj0jPj7e0MHBgweN4sWLW2Pr1KmT6iGRppDWCAkJsc6VyZMnG7rAWMxxYYw///yz6iG5DgN7Ohw5\ncsQoWbKkdfLVrVtX+cz9r7/+MiIjI/0u8F6+fFnpmEhvMTExVr4dv44dO1b1kIwxY8b4jWnGjBmq\nh+RKDOzptGfPHuOuu+6yAmmVKlVkcFUBM5o77rjDGssDDzxgXLhwQclYyF0QSM3zBrfo6Gj5rTTQ\ncMwBAwb4jQVjo/RhYM+AvXv3+s3cc+TIYYwePdq4du1aQI6PvCMuKpkzHNwaNGjAoE5pTn1kzpzZ\nOoew+mvTpk0BOz6O5bviDGPRKTXkRgzsGXT06FGjZs2afjMNpGb+/vtvx3OkvhfAcGvfvj3TL5Qu\nWNOeK1cu61zC6q93333X0dk7Hnvw4MHWSjPcMAY3rK/XHQO7Da5evWoMGzbMyJYtm3WChoWFGT16\n9DC2bNli23GwKWrFihVG69at/WbpuOFrLFFG7Ny5U563vudV2bJlZUrEzmtIZ8+eNT799FP52L7H\natiwobFv3z7bjhPMGNhttG3bNqNGjRp+JytuderUMWbNmpXuXXN4IyDFExUVleSxzVtERITyC7jk\nbq+//rp1PmXJksXv/MqZM6fRvXt3Y/Pmzel+fPzdbt26ycfyfWzM0idOnKj9bm43yYT/qClm4E2o\nGz169GjZHuzGmjLZs2eX1SLRJAAdYPBrRESELEGaNWtWkZCQIDvQ/Pnnn7JbDLrd4LZ161bZUPvG\nGtuDBg2SJVgnTZok/1/Pnj3FqFGjAvp8yRtwntWoUUOev/nz55dN1XEOf/XVV7L0sy/0/sW5a95w\nTqPNo1lKF41kzpw5I/744w/rHMbtxrZ8OOfbtGkjexPgfUA2Uv3J4lVXrlwx5s6dK79e3myWndZb\naGiozKOvXr3amt1gV2mRIkWs5WE//fST6qdOLoNcd6VKlazzDMsgfa8hvffee/IboV3nMR4LtZfw\n2OQMztgDADPwzz77TKxevVr8/vvv4vz586n+u4UKFZKzonr16sl66ugteaNvvvlGNG/eXN4vXbq0\nnCmxEQGl1nvvvWdVUGzatKms4Z8pUya/P4OZPPoRxMbGym+TO3bsSHWrSPQ/RWVRfEtt1aqVPEZI\nSIgjz4X+h4E9wPBmQKDHV9MNGzaIU6dOya+uiYmJMhjjFh4ebqVq8LX3xjdZctq2bStmzZol77/5\n5psBbb1H7rV9+3aZSsH5h5ZzSPvdddddt/x7aFhtploQ5OPj42Ua0Uw5hoWFyWBupmrY4CWwGNg9\n4vjx4yIqKkrExcXJ2RAa/1auXFn1sEhjmIWj+caaNWvkz2PHjhVdu3ZVPSyyARtteARSNiNHjrTe\nsJ06dUpy0YvI16effmoF9Tp16oguXbqoHhLZhDN2D8E/5WOPPWY1K0butH///qqHRRpC7160mEMK\nBem/TZs2iVKlSqkeFtmEgd1jsMQSb1hcoA0NDRUbN25kSzzyg7d8w4YNxbJly+TPw4cPF/369VM9\nLLIRUzEegwtfw4YNk/dxQQwpmdSuXqDgMHXqVCuo4yJ97969VQ+JbMYZuwchkGN55MqVK+XPn3zy\niXj55ZdVD4s0cOTIEXmRHRuIsmTJIpffVqhQQfWwyGacsXsQ1g1jN6q5lh3LH5FTpeCGOVz37t1l\nUAdcf2FQ9yYGdo/ChbB33nnHWnOMFQ/8chbcsLlo7ty58j5m7byw7l1MxXgY6svUqlVL7hQE7H7F\n7lUKPtjfgGCO/Q7Y8IZljqgNQ97EGbuHIYc6efJk+Sv06dNHHD16VPWwSIG+ffvKoA64WMqg7m0M\n7B5XsWJFmWMH5FZ79OihekgUYIsXLxbTp0+X9++++24xePBg1UMihzEVEwRQDhjlBbZt2yZ/nj17\ntizGRN6H/Qzly5cX+/fvlz8vX75cPPzww6qHRQ7jjD0IZMuWTaZkzGJimLWj+Bh5H76tmUG9c+fO\nDOpBgoE9SNSsWVP06tVL3j927JjMuZK3/fTTT7IeDBQpUkTuMKXgwFRMEMGyR6xb3rNnj5V7bdSo\nkephkQNQQhddtlAiGubNmyeaNWumelgUIJyxBxHUxJ44caL1M9a2p6XpB7kH9jCYQR21+hnUgwtn\n7EHoxRdfZJ9UD0OZgOrVq1v9S3HRHGWdKXgwsAchLHvEZhXUDTFzsbVr11Y9LLIBavBXq1ZNVvWE\nmJgY8cwzz6geFgUYUzFBCB3l0S3HhAqQaM9H7ofKnmZQR2/Rdu3aqR4SKcAZexDz7ZP6xhtviPff\nf1/1kEhB/1LyHgb2IMY+qd7tXzpu3Djx0ksvqR4WKcJUTBBjn1Tv9i/FBXIKXpyxBzn2SXU/9i+l\nGzGwE/ukuhj7l1JymIqhZPukIjVD+mP/UkoOZ+wksU+q+7B/Kd0MA7vCui3obLRjxw5x7tw5mR9V\nDatjsL4dnZeyZs0qq0BizTvpB29bLFXF+WNeMMUHM6UN+gOj1EaePHnkUlH0L8C573YM7AH0xx9/\nyFotWL2AC1xId+TOnVsGz+zZs8uTTIfgbnbawQkfERGhekiUDEwGDh06ZJVljoyMtMoyU+rhPYhJ\n1unTp62Lz0hpYSc2JjZ33nmncCMG9gBArY5BgwbJBhe4SNmgQQNZRhc3BE6d3pCYrWNc69evlz+z\nT6p+2L/UmVTkjh075Gv5yy+/yFVieJ27du0qN++Fh4cLV0FgJ2dcu3bN6NOnj5EpUyajXLlyRmxs\nrHH9+nVDdxs3bjSyZMmCD3zj9ttvNw4fPqx6SOSjQ4cO8t8GN5xfZL9Lly4Zn3zyiREeHm5kz57d\nmDBhguEmDOwOSUhIMNq2bWuEhoYakyZNMq5evWq4SXR0tBU8WrZsqXo49I9FixZZ/y533323ceHC\nBdVD8rSLFy8ar7/+uny9hwwZ4oqJGTCwOwBvtsaNGxs5cuQwli5darjR5cuXjaioKCuIzJ49W/WQ\ngt65c+eMiIgI699k+fLlqocUNEaOHClf81dffdUVwZ2B3QGdO3c28ubNa6xdu9ZwszVr1sg0Ek7o\nwoULG3FxcaqHFNR69OhhBfUXX3xR9XCCTkxMjJE5c2Zj3Lhxhu4Y2G22bNky+cabOXOm4QW9e/e2\ngslzzz2nejhBa+XKlda/Q5EiRYzTp0+rHlJQ6tevn5ErVy7jwIEDhs64KsaBnqJYCztnzhytVrvY\n1Sd10aJFonHjxqqHFVTYv1Qf8fHx8t+idOnSYv78+dq+x9UvnPYQFNA6deqUGDNmjLb/4Bntk4pS\nsOyTGljsX6qPsLAw+X5YsGCBnLzpijN2myQkJIgiRYqIXr16iYEDBwov90nFxo3Ro0erHlJQ9i9F\nM42CBQuqHlbQa9asmbhw4YL4/vvvhY44Y7cJvh6jZsfzzz8vvAhVA++44w6r9jf6pJKzUBv/hRde\nsAqyoX4Pg7o+E50ffvhB7N69W+iIgd0mmM0i9+zVVmQ39knt3Lkz+6QGuH/p008/rXpI9I9HH31U\nTnSmTJkidMTAblOjg6VLl8pg52XNmzcXTz31lLy/a9cu8fbbb6sekmch5YLcOqB/KT5UvXLdxguy\nZMkiv52jbDLKcOiGgd0GS5YskRcZ0YnI65AOQK7XTM8gB0z2MtsUoja++Tp79Zugm7Vt21aWTja/\nVemEgd0Ga9eulRXhvFDu81bYJ9V57F/qDlFRUfLbFN7/umFgtwH+YYOpul67du1EkyZNrFLEZvcl\nsiet9+abb8r7KCGLazc6lHOmpEJCQkS1atUY2L3o7NmzMh+KUrfBArnecePGydkKIBdsNnyg9MPK\n4y5dulhNVwYPHixKliypeliUArzvUeZXNwzsGbRu3Tr5hgymGTuwT6r9cCEOF+EBM0H2L9VfjRo1\n5EICbEzUCQN7Bm3dulWuLcbmpGCD2eVDDz0k769evVruuKX0wUW4vn37WisuJk+eLH8lvd17771W\nMx2dMLBn0NGjR4MyqANyv8gBIxcMyA0jR0xpg2982M2LDW7Qv39/NqV2ifB/OisdO3ZM6ISB3YbA\n7rq2WTYqVaqUtd4aBcMwi2eVirSJjY216o6gdSICO7lDaGioXP6LOKATBvYMCvbADn369BFVqlSR\n95EjnjZtmuohuQZys5itm9+AkIJBc2pyj/DwcAZ2r8FXsMKFC4tghlwwtlabOWEEeuSM6dbwWqEp\nNeBiabBdhPeCwoULMxXjNZyx/w9q0Jvrr5Er7tmzp+ohaW/x4sVi+vTp8v7dd99tpbTIXcI5Y/ee\n06dPi3z58qkehhbeeustuRsPvv76a5k7puShpj1q25tQ4xtlKch98uXLx+WOXoM13LiAQkLmhpEj\nNotVIXes2wmvC3y72b9/v7yPkgEPP/yw6iFROuH9r1tZDQb2DMCGnOvXrwdFjZi07MRDsxFA3tFc\nm03/D7XsUQ8GsFSWJRncLWvWrFbBNl0wsGeAWa6TG0n8vfvuuyIyMlLexwoZ5JLp//uX+pZ3Rjle\n1Londwf2K5yxew/rZPtjn9SbwwVSbEEH9i/1hkwavv8Z2MkR9evXl/VjALlkbrr5X/9S1FYHbGpB\nbXsiJzCwk2NGjBhh9UlF8+tVq1aJYIWv6r6F0ti/lJzEwE4B65OKwBasfVIxU0ftemD/UnIaAzs5\nin1S/9e/1Hze7F9KgcDATo5D2sHcxBVsfVLZv5RUYGAnxwVzn1TUqGf/Ugo0BnYKiGeeecavT6q5\nOsTL2L+UVOFZRkr6pCLn7OU+qWb/UtSoB/YvpUBiYKeACaY+qb79S6tWrcr+pRRQDOwUUMHQJ/XG\n/qW+teqJAoGBnQIqGPqkohY9+5eSSgzsFHBe7pM6e/ZsWYse2L+UVGFgJyW82CfVt38pLhazfymp\nwsBOSiDnjMDnpT6p7F9KumBgJ2UqVarkmT6pN/YvxfJGIlUY2En7PqkoHGbOhFXCEs3kmhazfynp\nhoGdtO2Tiguqn332mSz9W7RoUbFy5Upl48R6++rVq8uxPPvssyIuLi7Z/qXojsT+paQaAztp2ScV\nSyAbNWokXnjhBZmmQRvC+fPnKxvjX3/9JTZu3Cjvx8TEiLJly4r//ve/8sPGt39pMJRKIP0xsJOW\nfVIROM2dm6YDBw4oGl3SY584cUKWIzbr3wD7l5IuGNhJC8hJDxgwwPo5uYYcOgV2k1kLplq1auLx\nxx8P8KiIksfATsqhhO+QIUNEt27dUvxzOgZ207p160SDBg3E33//HbAxEd0MAzsphbw1LkpidYzZ\njOJmDh8+rKxoWGo+VJYvXy7LB6D2vFd20pI7MbCTMgkJCeKRRx6xeoHeCoK6qk1MBw8eTNWfi4+P\nl5uTvvrqK8fHRHQzDOyktCBYWFhYmv6OqnRMWo/LdeykEgM7KZM1a1axatUq0bp1a88E9uLFi4u5\nc+eKpk2bOj4mopthYCflzTewHvz7778X5cuXdyywI+d9/fr1dP1d7Cw9e/Zsin8me/bssozAtm3b\nRPPmzUUwwmvMawt6YPV/0kK9evXEhg0bZPu86Ohoq555WgI7ggo2Nq1fv1789ttv8tfNmzeLc+fO\nyeWT+H18S0AQRj0XVJc0bxUrVrRqxKflmID17NiYhA8pr8NGMbQ0xGtr3vAzln3imom5mxipqDJl\nylivL7pI4eeQkBDVTyE4GJRuly9fxvTEmDt3ruqheMrx48eNLl26GJkyZZKvr++tQYMGSf78jh07\njF69ehkFChRI8udTe8uSJYvRsGFDY86cOcaVK1f8Hn/RokXJ/p0KFSoYK1asMLzu+vXrxqpVq4x2\n7doZYWFh6X6N8XfxGHgsPKZXREdHG+XKlTN0whk7aadgwYJi/PjxsrDWyy+/LFvomX7//Xdr5vjN\nN9/I1npYZpic0NBQufywUKFCcpaO2SJm7kitbN26Ve4eNeHxvvvuO3m788475bFR9yU8PDzJ4+fN\nm1emXfBnvNzy7sKFC2LGjBnyNd60aVOyfwbNyfEa4zXBawyXLl0Sp0+flt+W8Fr7rhj64osv5A2V\nPbt37y7atWsncubMGbDnFDRUf7K4GWfszsPMbvr06Ua2bNnka12qVCnju+++MyIjI5PMCIsUKWJ0\n7drVmDhxovH7778bCQkJKT7u/v375Qx9wIABxoMPPpjk8bJmzWq8/vrrxqhRo6z/17p1a+PEiROG\nl+G1GTdunJEnT54krwlmpn379jVmzJghvyldu3btpo+D38OfwZ/F34mKikryeDjG+PHjXT2Dj9Zw\nxs7AngEM7IFz5swZ4+OPPzY6duyYJDjUr1/fiI2NNRITEzN0jK1btxo9e/Y0cuXK5ff4ZcuWNUaM\nGGH8/PPPhtft3bvXeOSRR5Kkqdq2bWusXLkyQwEYf/fHH380nnrqKfmYvsfAMXFsN4pmYPcWBvbA\nwSw9IiLCLxggyG/fvt32Y50/f94YM2aMkT9/futYmTNnlrP3S5cuGV6epefMmdN6zqGhofLbzJEj\nR2w/Hh4Tj41jmMfDB6obZ+/RDOzewsAeGO+//75fQL/zzjuNxYsXO37co0ePGi1btvQ7drVq1TyX\nisE3HVzU9H2eVatWNbZs2eL4sTdv3iyP5XtsjCWj374CiYHdYxjYnYWZG2bJvm/6zp07y7RMIMcw\nc+ZMv9k7csUHDx40vCA+Pt5o2rSp3yx9yJAhSVYGOQnHwjF9Z++PPfaYHJsbRDOwewsDu7PeeOMN\nvwuZs2bNUjaWw4cPG5UqVbLGU7p0aePYsWOGm+Hi8qOPPmo9p4IFCxrr169XNh4cG2Mwx9OkSZMU\nL4DrIlrDwM6dp6SlDz74QAwdOlTexzI6dE9q06aNsvGgJd6KFStErVq15M87d+6UHZ5utSNVV9iF\n26FDB7Fo0SL5M5Z4/vTTT6Jy5crKxoRjYwwYCyxcuFA899xz6d4xHMwY2Ek7P/zwg3jjjTfkfewU\nnTNnjgyiqqE70uLFi63gh6qU6NHqRh9//LGYNWuWvI91/ijpULp0adXDkmPAvgGMCWbOnCnLIFPa\nMLCTdpti0OfUNH36dC2Cuil37twyuKMkAWADz7x584Sb4NsG6t+b34aWLFkiSpUqJXRxzz33yDGZ\nG5769+8vdu3apXpYrsLATlrBTB31XgCpgrZt2wodd8ZOnTpVZMqUSf6MHahxcXHCDVDT/vnnn7da\nDyLdde+99wrdYEzvv/++vI+xYsyqmqy4EQM7aZWC+fTTT62cNtIFunrwwQdluQM4duyY6NWrl3AD\npDXWrFljPYeePXsKXeH1xRgBZSU++eQT1UNyDQZ20gIukHXt2tX6ecKECbL+iM7Qp7VEiRJWSgYf\nTDrDB5BvCmbKlCmy2YmuMDaM0Tclg+dAt6bvvyoFFRTfMvOoSL889thjQncoTTtx4kTr51GjRgmd\nTZo0yUrBvPPOO6JkyZJCdxgjxgoY++TJk1UPyRUY2EkLqCBoev3114Wb6sij1jjgImpqe6MGGqpX\nomKmWZHR99uR7nANA2MG1Otnrv3WGNhJOVws/fbbb+V9rBPX8WJeSlB+1kwnIYWkowULFlgNQ7A2\n3E2lchHUcSEd8BzwXChlDOykHIKh2VLNDJJ2Q5cfbHDKnz+/XLKIVM++fftseWx0UDKvByA1k5iY\nKHT+RtStWzfhNr5j9n0ulDwGdlIOjRegQIEC4sknn7T98bHhBZuK1q5dKzp27CgbTWNXY8OGDW0J\nwmFhYXI5Hhw9elTuUNUJGorgGgbUrVtXREVFCbcpV66cqFOnjryPNe6+TVIoKQZ2Uur48ePWzBkb\nkW7WdzS98NW9VatWMr2zZcsW8eGHH4ovv/xS9OnTR16stWtzUbNmzaz769atEzpBX9Lkxuk2vmP3\nfU6UFAM7KeX7BkXTY7u99tprsiVbTEyMdQHOTJ+Ab9u9jLjvvvusDUu6BR009nbyNQ4U8yK1jq+x\nbhjYSSnfN6jvG9eui7JfffWVaNmypVUCwJQvXz7rG4MdkLfHVngdg445Hnzw4APIbn/99ZcsAxER\nESFr++A4N942bNjg6Q9P3Xi3Ey+5gpNBB0WusFIFDZNvlJCQIH/Nli2bbcfDbBh1WPbv3y9Onjwp\nrxno9BqjwJbvtxY7oLk4lnxevHhRPPHEE7LmzK+//mo1AC9evLgM9nbk9TF2fHjiNWZgTxkDOym1\ndetW+SsCgt1L8MwLhigFiwDky1xvXqRIEduOhwu05oVg5PNxoVK1c+fOWcsc7f7gxIcjVhpdunRJ\nBnLz4ibgWxKqco4ePVperLbzNTY/PM+fP2/7B5VXMLCTUnhzmoW17ITlk5g5wogRI27658qUKWPb\nMX1n6KhSqQPfcdj9Gn/++ecyDTNgwAC/oO4b2PGBamdgv/E1ZmBPHnPspJS5xd3u1TCYkeONj40t\n/3QK87uZ6RmzcYYdzJomgFmsTq+vE68xrl8ghZbc3gPz25fdu0R9n4PvcyN/DOyklNkdx7woZpcj\nR45YVSKTOyZSB7igahbxsoNvQS1dtr37dh+y+zX+5Zdf5Pry5F7jQ4cOyV8jIyNtPaaOr7GOGNhJ\nKXMGZl7MtMuVK1fkr6GhoUl+D0EdVQKTu6iaEb6zdN/Zu0q+M1w7X+NTp07J/H14eHiyv49mJGD3\ndQbfWbour7GOGNhJKfPNaeba7WLmk5Nbzohyu9gtatZTt4vvc9Al6PiOw87X2Jz9I7jfaPv27fLC\nNVbLFCtWTNjJ9znYnVryEgZ2UspMhSAYmLNsu8q9om/mN99843cB8b333pNb/vGr2VfTLps3b/Y7\nvg6wXh+9Wm8cX0ahNg6aTmPZ4e7du63/j+beuK6BFBAaktvNfA44vrkXgZJiYCelzJ2QSBNs27bN\n1lxs3759Za4dF0hRCrhBgwZyBUeXLl0c6Xhkrq1G0LE7t5yRmbXZfHvjxo22fnjiNUWe+4EHHhC9\ne/cWr7zyilxlhCbfKIZWrVo1YSeMfdOmTfI+npPd1wy8hIGdlPLd4m73ppN//etf4t133xVnzpyR\nbdWQF542bZqsS253UNA56Dj14Ym2emhliIqZY8eOlWv477//ftl6D8XWnNjzYF4ncHNphEDgOnZS\nyvcNipom2JpuFwRXtIIz28E5Seegc+NrXKlSJdseG0sdnSq1HOi6Ql7CGTspz7GbOWD0DDXrsrvN\n999/r23Q8a3B4ztOt9H5NdYNAzsphVm1uTMRzTBWrlwp3AYXCtGyzaw9U79+faETrNcvW7asvD97\n9mxZx8ZtUH8dYwc8lxuLupE/BnZSzvervBu742Am+eeff1qNuJFz1u3D0+xAhMYiU6ZMEW6DMZtN\nUXC+6HQNQ0cM7KQcVq2Yed+vv/7a2jXqFr4fRj169BA6whJErN13Y0NojNX8RpQjRw7x7LPPqh6S\n9hjYSTnfeiNXr1613sRugO5PZhcm5LLtXuJnlzx58oj27dvL+3v27HFVQ2g0OkdtfcBzwHOhlDGw\nkxawvd98ww4bNky2rdMdLvR27drVqscSqNUh6eU7Pqw716UCZUowRozVzY24VWBgJy2gGuA777xj\n1QNBc2jd0wVTp061aqKgp6o5I9YV0l1mGgOz9jfeeEPoDpugzNk60kl2LtX0MgZ20gY2vDz44INW\nL9KRI0cKXaEssDmTzJIliwzy6BSkO7ymZjVGbC7CElOdL0qb1y/QEOXjjz9WPSTXYGAnbaAMAFY/\nmIWrsLHIzp2SdsE3ic6dO1sFsKKjo10zk0S5A+y8NWFDGOq76AZj6tSpk/XzhAkT5NgpdRjYSSso\nnjV06FArJdO4cWPZBk2nvDpWvixZssRKwbz55pvCTR5//HErJYM0B36Oj48XusBYHnvsMb8UjJ1d\nmIIBAztpmZLBGxvQrxPFu8zGDaqDer9+/awZLy72xsTEuCIFk1xKBs2tzZ6waGWnQ9cnBHWMZdWq\nVfJnjJEpmLRjYCctUzKzZs0SDz30kPwZK2Rq164t+2uqTL9gBcxHH30kf0a6CEsG0UHIjZDWWLp0\nqYiIiJA/4xsIvh0lV189UHDsRx991Po2hFruy5YtYwomHRjYSUvYTDN//nxZLRDwtbxGjRriv//9\nr5K16gh6yPOaY8PadZSrdbO77rpLdpPCr4ByDniNzSbggbR27Vp5bLOkBMaEoI6a75R2DOykrdy5\nc8tOPA0bNpQ/x8XFiaeeekq0bt062c5ITqRekHYpX768DDJm+gUzXaSHvHJNA2mPe+65x6rXg53A\nWAoZiGbROAaWNOIDHMc20y8///yzNs1KXMmgdLt8+TJKERpz585VPRRPS0hIMF577TUjc+bM8vXG\nrUCBAsaMGTOMq1evOnLMnTt3GvXr17eOh1vVqlWNbdu2GV504sQJo2XLln7Pt0yZMsbKlSuN69ev\n2348POaPP/4oj+F7zFatWsmxuEl0dLRRrlw5QycM7BnAwB5Ya9asMUqXLu0XCIoXL24MHTrUOH78\neIYfHx8S8+bNMxo1auR3jNDQUGPIkCHGlStXDC/bunWrERIS4vfccatSpYoxefJk4+LFixk+Bh5j\n0qRJRuXKlf2OkT9/fmPmzJmOfIg4jYHdYxjYAy8+Pj7J7N0Mvs8884wRGxtr7Nu3L9UB4vz583JW\nOnjwYCMiIiJJUMMsfcuWLYbX4UOtVq1a1vO+9957k7wWt99+u9G7d29j8eLFaZpV48/i7+Dv4jFu\nfFzM0o8dO+bo8wu2wJ4J/1GdDnIrdMxBp/S5c+eK5s2bqx5OUNmyZYvcOfn555+LixcvJvn9AgUK\nyGYMaFOHptVYxRISEiJzuqg/gr+Pjjw7d+5MtrkHVuSgLsmTTz4pd5Z6HVoHmn1g69SpIy+q/vjj\nj/I1xoXi5Mo7YEUNXmOs5UezFHNjGZZNoh0hep/iNU5uHwL+LVq0aCHr19SrV8/VZXgHDhwoq5Li\nnNKG6k8WN+OMXb0zZ84Yo0aNMsqWLZtkJpjWW86cOY3u3bsbmzdvNoLJnj17jLCwMPka3Hbbbcau\nXbv8fv/AgQPGwIEDjfDw8Ay/xnfccYcxaNAg4+DBg4ZXRGs4Y/f+VIQ8DatUsKEJu0E3bNgg1q1b\nJ/t6YqaIGRSaTN8MltRhxmnesFY+V65cIpjg20qXLl2snacoxFaqVCm/P4Mlh2+//bYYMGCAbFRt\nvr64YY/Bzb70YxaO1Tbm64uyxlhx48YNXW7DwE6egCCCtAtuL730kpUqwxI6pFuwTNLscITiXZGR\nkTJFE+xQvAzLNwGBt0+fPjf9swjISFGZG8fMTUUI7vhgMHeuIiWDtf5YthhsH5S6YGAnz0L/URTn\n8g0uyAtjIwwJ2amqb9++8j6uI0yePDnN1xOw18C3WTbpgRuUiIKQWcwMFzmhf//+omLFiqqHRTZh\nYCcKQrGxsWLOnDnyflRUlAzs5B0M7ERB5tSpU1bTbVybQAoGaSvyDgZ2oiCDC6RmrR1cSK5Zs6bq\nIZHNGNiJggh6tE6fPl3ev/vuu8XgwYNVD4kcwMBOFCTOnz9vLQWFiRMnihw5cigdEzmDgZ0oSKCF\nn7m9Hz1bH374YdVDIocwsBMFAbS/Q90XKFKkiBg+fLjqIZGDGNiJPA47QjFDN40dO1YW7SLvYmAn\n8jjUf8G2f0BphWbNmqkeEjmMgZ3Iw37//Xcr7ZI/f35Znpe8j4GdyKNQ2fKFF16waqmPHDmShc+C\nBAM7kUdhpr5x40Z5v0mTJqJdu3aqh0QBwsBO5EHbt2+XNdQB1S3HjRvn6i5FlDYM7EQeg9RLp06d\nRGJiovx52LBhsqkIBQ8GdiKPwXp1dDoy+5eiQxIFFwZ2Ig/Zu3ev3GEKaLSOsgGZM/NtHmz4L04U\nRP1LKTgwsBMFYf9S8jYGdiIPsKN/KXkHAzuRy7F/Kd2IgZ3I5di/lG7EwE7kYuxfSslhYCdyMfYv\npeQwsBO5FPuX0s0wsBO5EPuXUkoY2IlciP1LKSUM7EQuw/6ldCsM7EQuwv6llBoM7EQuwv6llBoM\n7EQuwf6llFoM7EQuwP6llBYM7EQuwP6llBYM7ESaY/9SSisGdiKNsX8ppQcDO5HG2L+U0oOBnUhT\n7F9K6cWzhEhD7F9KGcHATqQh9i+ljGBgJ9IM+5dSRjGwE2mE/UvJDgzsRBph/1KyAwM7kSbYv5Ts\nwsBOpAn2LyW7MLATaYD9S8lODOxEirF/KdmNgZ1IMfYvJbsxsBMpxP6l5AQGdiJF2L+UnMLATqQI\n+5eSUxjYiRRg/1JyEgM7UYCxfyk5jYGdKMDYv5ScxsBOFEDsX0qBwMBOFCDsX0qBwsBOFCDsX0qB\nwsBOFADsX0qBxDOLyGHsX0qBxsBO5DD2L6VAY2AnchD7l5IKDOxEDmH/UlKFgZ3IIexfSqowsBM5\ngP1LSSUGdiIHsH8pqcTATmQz9i8l1RjYiWzE/qWkA667Ik/WZDE3A8GFCxes+wkJCTL4mrJmzSp3\ngtqF/UtJB5yxk6fs2LFD9g7NnTu3datUqZL1+6h97vt7OXPmFO+9954tx2b/UtIFAzt5ytq1a62L\nlqmd3S9YsCDDx2X/UtIJAzt5ymOPPSbrnKdF+/btM3xc9i8lnTCwk6egf+grr7yS6j9ftGhRWSM9\nNZYtWyaGDBkijh496vf/2b+UdMPATp5cQ57aWTt2g6Zm49C5c+fkLPytt96Su0ixnBElA9i/lHTE\nwE5BO2tPy2z9r7/+knl0OH36tHjuuefEo48+Kj8Y2L+UdMPATkE7a0/tbB0OHDiQ5P8tWbJEjBgx\nQt7H6hr2LyVdMLBTUM7a0zJbv1lg91W4cGG/tfNEKjGwU1DO2tMyW09NYEeqBuvlhw4dKvPuRCox\nsFPQzdrTOltPTWA3d7Vi5ykKfsXFxaXp8YnsxMBOnp+1I/+dkdl6agO77/LHmJiYND0+kZ0Y2Mnz\ns/YOHTpYPyPIp3W2ntbAHhERwQ1KpBQDO3newIEDRUhIiLz/4osvpnm2fv36dXHo0KFb/jk87oAB\nA8S2bdtEZGRkusdLlFGs7kiehM1DmGVv3rxZVnN8//33xcGDB0XlypVlyzoE3vLly4vQ0NBbPtax\nY8dueUG0RYsW4sMPP5T114lUY2AnT7h8+bL47rvvxK+//irWr18vfvvtN3Hy5MkU/w6CeoUKFUSV\nKlXkrVGjRqJYsWJJ/hw+EG6mdOnSsoRAw4YNbXkeRHZgYCdX27Nnj9wYhJ6iaV2JkpiYKD8EcANs\nLmratKno3r27DPKZM/8vU7l9+/YkfxfLKP/973+Lnj17pmrWTxRIDOzkSkuXLhUff/yxWLRokUy7\n+EI+vWTJkqJs2bKiTJkyIm/evDL/jUCNYI6NRH///bcM2KjffvbsWfn38DjffvutvCGl0rVrV9Gt\nWzfxyy+/+D1+x44dZWonPDw8oM+ZKLUY2MlVkF55+eWXxcyZM/3+P4J38+bNRb169WRQT+0FUgTz\nI0eOyFn7nDlzxKZNm+T/R+D/17/+JcaMGSOio6PlDD179uxi1qxZom7duo48NyK7MLCTayDwYhbt\n20gDuz1bt24tW9ClJyWC9Au6HeH2+OOPyxk8Lq7imwDy9nv37pXLI5Ge+eCDD5KsiSfSEZc7kvZQ\nMvfpp58WLVu2tII61oqjUTRy640bN7Ytz43UDUrzoqsSmnaYMHPHhdZVq1bZchwiJzGwk9ZOnDgh\n0ytm6gUzbJTG/eKLL8R9993n2HHz5MkjL44ij1+wYEH5/zB7f+SRR8S8efMcOy6RHRjYSVvoVPTQ\nQw/JLfpmjRfM0vv27Stuu+22gIyhdu3aMq+OWutmPZhWrVqJL7/8MiDHJ0oPBnbS0qlTp0SDBg1k\nzhvuuece8dlnn4l777034GPJnTu3ePvtt8VLL70kf0a3pGeffVZ88803AR8LUWowsJN2sCQR+e0t\nW7bIn8uVKyfGjx8v8uXLp2xMSAGhHEHv3r2t4N6mTRuxevVqZWMiuhkGdtLOe++9J9asWSPvlyhR\nQu7sTG0PU6e1b9/emrkjLYMCYxcvXlQ9LCI/DOyklQ0bNoghQ4bI+1haiKCOC5k66dy5s1yJYzbY\nQBlgIp0wsJNWKRjs6rx69ar8uV+/frLlnG6QlnnttddkSWDAh8/KlStVD4vIwsBO2sBM3dz5idUo\nqNuiK3yLQLck0/PPP8+UDGmDgZ20cPr0aTFs2DArBYP0BmbGOkNpgUcffdQqQTBlyhTVQyKSGNhJ\nC1OnThWXLl2yctiFChUSbmm9Z+56xe7UGwuSEanAwE7KoUPR2LFj5X0U73JTWzkswcRuVMCa+xUr\nVqgeEhEDO6m3fPly8eeff8r7qIOODUFugiJkJszaiVRjYCflzNn6jUHSLdBiD52UzAqUhw8fVj0k\nCnIM7KQUdnCipR0gOKI5hhNeffVVcf/991tLKe2Ei7xPPPGE9Xx++OEH249BlBYM7KTUrl27rGWC\naDTt5HHQFSlLFmdaEPhWmjRb7RGpwsBOSvkGQdRCd8L58+dll6RSpUoJpxQvXtyqOMnATqoxsJNS\nv/32m3XfiTQM2tqhnjvMnz9fVK1aVd5Q48VO6LNq5tlRZhgrfYhUYWs8Usqc3aKfaLFixWx/fATx\nM2fOyKJiaH1nNqBGcTG74YNp48aN4sKFCzL149Q3EKJbYWAnpfbt22elMjDrtRsaXO/cuVMG9h49\neogCBQoIp0RGRvo9LwZ2UoWpGFIKDaPNGbtTdu/eLfLmzetoUL/xOZi7aIlUYGAnLQJ71qxZHTsG\nNj85eeHU5NtQ23xeRCowsJNS5vJDpy42YjUMVsUEIrD7PgenllUSpQYDOyllLhFENyInmKUKAhHY\nfZ+Dk6klolthYCelzLowcXFxjuXXAxXYfZ+D2+rdkLcwsJNSFStWlL8eOnRInDt3zvbHP3v2rPw1\nR44cwmmo7miqUKGC48cjuhkGdlIK68yTC4x2ueeee+SvAwYMkMXGxo8fL44fPy6csH37dvlryZIl\nxe233+7IMYhSg4GdlKpSpUqSwGgndDh69tlnZTD/7LPPxMSJE628vt3fDPCt48bnRKQCL92TUr6F\nv5wI7JkzZxa9evWSNyf5fttgYCfVOGMnpbBxyEyXYHeoWxtCL1u2zLpfvXp1pWMhYmAn5cyCXAjq\nixYtEm6D2jDmuFEaoXbt2qqHREGOgZ2U69Spk7XzdPbs2a5rCP3tt99aO027du3qSM0borRgYCfl\nUHGxVatW1rpzVEh0C3wI4cPILCnwwgsvqB4SEQM76aF79+7W/SlTprhm1r5ixQqxd+9eef+pp54S\nBQsWVD0kIgZ20gPy0uZqktWrV4sFCxYI3WGJ49ChQ62fX3nlFaXjITIxsJMW0BB6woQJVn76ww8/\nFCdOnBA6GzFihFVG4OWXX/bbbEWkEgM7abWmvX///vI+KjIOGTJE25TMjz/+aK2EQZPs999/X/WQ\niCwM7KQVbP0366z89NNPYubMmUI32GGKDx3fawKBqEVDlFoM7KQVrCyZOnWqVc8cKRksJ9TFyZMn\n5YVe3xRMnTp1VA+LyA8DO2mZkpk0aZL18zvvvKNFcD969Kh46aWXrJow9evXF8OHD1c9LKIkGNhJ\nS88995wYOXKk1Zno3//+t5gxY4aynPvff/8tOnfubDXfrlGjhpg7d67Ili2bkvEQpYSBnbSF5YOf\nfPKJ9fN//vMfWcwLM+dAuXbtmvj8889F+/btreNiaeaSJUtEzpw5AzYOorRgYCetIYc9bdo0q+QA\n1rhjIxBmy07P3rHxCLN0fHNITEyU/69p06YyqOfJk8fRYxNlBAM7uaJI2G+//Sbuu+8+q1jYu+++\nK3r06CHWrl1re4DHzPzTTz8V7dq1E5s3b5b/LywsTIwaNUp888038j6RzliPnVzTQg9BHDs9Bw8e\nLK5cuSJ+/fVXeYuIiBBPPvmkePzxx0WuXLnS9fjI4+OxvvrqK7nMEj+bsOpl8uTJokSJEjY+IyLn\nMLCTayAdEx0dLZo3by5Xp/zyyy/y/+/fv1989NFHcpZdrVo1UbZsWet2s9otqMa4a9cu2dwDTTJ+\n//13a7WLKV++fOLtt9+WyxvRsIPILRjYyZWzdzTlWLdunRgzZozcxIRAnZCQIFatWiVvpvz588tm\nHli9guCMXPmlS5dkEMeF0eSgUQaCeZs2bUT27NkD+MyI7MHATq6F2Tn6mKJmCzY1YfXKli1b/AI2\nNhKZm4luVTq4SZMmolu3bqz5Qq7HwE6uh1n5q6++Km/x8fGynvv69evlDRc/UXcGs3QEfDSyxiw8\nMjJSVpM0b0WKFFH9NIhsw8CewYqE4HuhjdTCipVatWrJG1EgYMKg2zUYvUbjMubaaqzQIKLgdOXK\nFVnjSCcM7BmcsaNYlbl5hYiCT2JiojXJ0wUDewYhX4v8LREFp0uXLslrNzphYM+gQoUKiePHj6se\nBhEpcvz4cVG4cGGhEwb2DMIyuUAWpSIivRw9elTGAZ0wsGcQAztRcDvKwO49DOxE7lO3bl15yygU\noGNg9yDk1o4dO6Z6GESkwPnz52U5C+bYPebOO+8UBw4c4Fp2oiC0Z88eKw7ohIHdhnol+MRGjRIi\n0s/MmTNFmTJlZCG4cuXKiTlz5tj22Cj1jN3OqCSqEwb2DMI/KFqkmSVkiUgfy5Ytkw1TSpUqJb7+\n+mvx2muvyfaKO3futOXx8b5H0ThsVNSJXqNxoZCQEDlrRxMIVAYkIn0MGjRIztbnzZtn1XPBz6gl\nVLp06Qw/Pt73qAqqG87YbVCzZk3O2Ik0LM61bt062V3Lt0gX3q/FixfP8OOfO3dObNu2TT6ebhjY\nbXD//ffLr3YHDx5UPRQi+sfJkyfloobkVqzYsYplxYoVcrkjA7tHNWjQQLZRQ7MHItJDgQIFZHGu\n5JYj27FEecqUKXItvI61/BnYbYCr7R06dJANj1mbnUiv61+zZ8/2e18iL753794MPfaRI0fEt99+\nKzp37ix0xMBuk06dOsmT5fvvv1c9FCL6B5qRo1l5ixYtxIIFC+S3avSyzehO0WnTpolcuXKJli1b\nCh0xsNukfPnyokaNGuKTTz5RPRQi+scjjzwiZsyYIa+BIQgPHz5cfPzxxxlaEYN9KxMmTBDt27fX\nttl5JgPZf7LFkiVLROPGjUVsbKy2n+RElDFvvfWWnMBhU2KxYsWEjhjYbdaxY0cZ4LEMKm/evKqH\nQ0Q22rBhg8zbjxw5UvTo0UPoioHdZqdOnRJRUVGiYcOGMg9nNrwmIne7fPmyXNqMEgIrV67UroG1\nL+48tRmWPSL/9sQTT4gcOXKI0aNHy6vzROReZ8+eFc2bN5dFv7AZUeegDnqPzqWaNWsml1hhnSvq\nVLDZNZG7W9/Vq1dPXoD98ccfbSlF4DQGdodgxr548WKxcOFC8cADD8hiRMx6EbnH1atXZToVOXXM\n2FetWiUqVqwo3ICB3UH4lP/5559lega7U7FLDQE+ISFB9dCIKIXmGV988YUs8Yv9KfXr15fv4xIl\nSgi34MXTAMGnfXR0tKwvERoaKipXrixrTFSqVEnkyZNH5uN1z9sRebFQ2Pnz50VcXJxYv369zJ+b\nvRWefvppWR3ynnvuEW7DwE5E5DGcIhIReQwDOxGRxzCwExF5DAM7EZHHMLATEXkMAzsRkccwsBMR\neQwDOxGRxzCwExEJb/k/sNyry7XAGmMAAAAASUVORK5CYII=\n"
        }
      },
      "cell_type": "markdown",
      "id": "ecc57330-7185-42ab-ad4a-6fdbdaaef5b5",
      "metadata": {},
      "source": [
        "When `tau_scale` is small (e.g. ${ \\dot\\sigma_{\\tau}{=}1 }$), the\n",
        "distribution ${ \\theta_d ~\\sim&~ \\mathcal{N}(\\mu, \\tau) }$ will be\n",
        "narrow, which will pressure the individual ${ \\theta_d }$ to be similar\n",
        "to each other. When `tau_scale` is large\n",
        "(e.g. ${ \\dot\\sigma_{\\tau}{=}20 }$), the ${ \\theta_d }$ can diverge more\n",
        "freely resulting in values quite different from each other. Think about\n",
        "how this will interact with the number of observations and the values of\n",
        "those observations.\n",
        "\n",
        "When `mu_scale` is small (e.g. ${ \\dot\\sigma_{\\mu}{=}1 }$) the location\n",
        "of the distribution over $\\theta_d$ will be be centered near zero. When\n",
        "`mu_scale` is large (e.g. ${ \\dot\\sigma_{\\mu}{=}20 }$), the center of\n",
        "${ \\mathcal{N}(\\mu, \\tau) }$ can move to new locations.\n",
        "\n",
        "In this way, a large ${ \\dot\\sigma_{\\mu} }$ and small\n",
        "${ \\dot\\sigma_{\\tau} }$ would allow ${ \\mathcal{N}(\\mu, \\tau) }$ to move\n",
        "to easily move to a new location that accommodates the broader\n",
        "tendencies of departments and professors, but would constrain the\n",
        "dispersion between departments. Notice how increasing `mu_scale` allows\n",
        "the posterior ${ P(\\mu \\mid t) }$ to move away from zero towards the\n",
        "mean of the combined data, whereas decreasing `mu_scale` pulls the\n",
        "expectation ${ \\operatorname{E}[\\mu \\mid t] }$ towards zero. Also notice\n",
        "how this effect is stronger when `tau_scale` is small.\n",
        "\n",
        "Pay close attention to the interaction between `mu_scale`, `tau_scale`,\n",
        "and the observations for each department (including the expectation and\n",
        "variance of the observations of each department).\n",
        "\n",
        "### Multiple observations with missing data and learned variance\n",
        "\n",
        "$$\n",
        "\\begin{align*}\n",
        "\\mu ~\\sim&~ \\mathcal{N}(0, \\dot\\sigma_{\\mu}) \\\\\n",
        "\\tau ~\\sim&~ \\text{HalfCauchy}(\\dot\\sigma_{\\tau}) \\\\\n",
        "\\theta_d ~\\sim&~ \\mathcal{N}(\\mu, \\tau) \\\\\n",
        "\\sigma_d ~\\sim&~ \\text{HalfCauchy}(5) \\\\\n",
        "t_{(d,i)} ~\\sim&~ \\mathcal{N}(\\theta_d, \\sigma_d) \\\\\n",
        "d ~\\in&~ \\{ 0, 1, 2\\}\n",
        "\\end{align*}\n",
        "$$\n",
        "\n",
        "![](attachment:generated/multilevel-models/partial-pooling-learnedsigma.png)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "c13c06ca",
      "metadata": {},
      "outputs": [],
      "source": [
        "%reset -f\n",
        "import jax\n",
        "import jax.numpy as jnp\n",
        "from memo import memo\n",
        "from enum import IntEnum\n",
        "from jax.scipy.stats.norm import pdf as normpdf\n",
        "from jax.scipy.stats.norm import logpdf as normlogpdf\n",
        "from jax.scipy.stats.cauchy import pdf as cauchypdf\n",
        "from matplotlib import pyplot as plt\n",
        "\n",
        "normpdfjit = jax.jit(normpdf)\n",
        "\n",
        "class Department(IntEnum):\n",
        "    GOVERNMENT = 0\n",
        "    ENGLISH = 1\n",
        "    MATH = 2\n",
        "\n",
        "t = jnp.array([\n",
        "    [10, 0, -10], \n",
        "    [-16, -15, -14], \n",
        "    [30, jnp.nan, jnp.nan],\n",
        "])\n",
        "\n",
        "Mu = jnp.linspace(-25, 25, 100)\n",
        "Tau = jnp.linspace(1, 30, 100)\n",
        "Theta = jnp.linspace(-30, 40, 200)\n",
        "Sigma = jnp.linspace(1, 30, 100)  ###NEW\n",
        "\n",
        "@jax.jit\n",
        "def half_cauchy(x, scale=1.0):\n",
        "    return 2 * cauchypdf(x, 0, scale)\n",
        "\n",
        "@jax.jit\n",
        "def professor_arrival_likelihood(d, theta, sigma): ###NEW\n",
        "    ### likelihood of a professor from department d \n",
        "    ### showing up t_d minutes early/late, \n",
        "    ### for a given theta and sigma.\n",
        "    return jnp.exp(jnp.nansum(normlogpdf(t[d], theta, sigma)))\n",
        "\n",
        "@memo(cache=True)\n",
        "def department_model[_mu: Mu, _tau: Tau](d):\n",
        "    department: knows(_mu, _tau)\n",
        "    department: chooses(theta in Theta, wpp=normpdfjit(theta, _mu, _tau))\n",
        "    department: chooses(sigma in Sigma, wpp=half_cauchy(sigma, 5))  ###NEW\n",
        "    return E[ department[professor_arrival_likelihood(d, theta, sigma)] ]  ###NEW\n",
        "\n",
        "@memo(cache=True)\n",
        "def partial_pooling[_mu: Mu, _tau: Tau](mu_scale=5, tau_scale=5):\n",
        "    president: knows(_mu, _tau)\n",
        "    president: thinks[\n",
        "        population: chooses(mu in Mu, wpp=normpdfjit(mu, 0, mu_scale)),\n",
        "        population: chooses(tau in Tau, wpp=half_cauchy(tau, tau_scale)),\n",
        "    ]\n",
        "    president: observes_event(wpp=department_model[population.mu, population.tau]({Department.GOVERNMENT}))\n",
        "    president: observes_event(wpp=department_model[population.mu, population.tau]({Department.ENGLISH}))\n",
        "    president: observes_event(wpp=department_model[population.mu, population.tau]({Department.MATH}))\n",
        "    return president[Pr[population.mu == _mu, population.tau == _tau]]\n",
        "\n",
        "@memo(cache=True)\n",
        "def department_model_theta[_theta: Theta](d, mu_scale, tau_scale):\n",
        "    obs: thinks[\n",
        "        department: chooses(mu in Mu, tau in Tau, wpp=partial_pooling[mu, tau](mu_scale, tau_scale)),\n",
        "        department: chooses(theta in Theta, wpp=normpdfjit(theta, mu, tau)),\n",
        "        department: chooses(sigma in Sigma, wpp=half_cauchy(sigma, 5)),  ###NEW\n",
        "    ]\n",
        "    obs: observes_event(wpp=professor_arrival_likelihood(d, department.theta, department.sigma))  ###NEW\n",
        "    obs: knows(_theta)\n",
        "    return obs[Pr[department.theta == _theta]]\n",
        "\n",
        "@memo(cache=True)\n",
        "def department_model_sigma[_sigma: Sigma](d, mu_scale, tau_scale):  ###NEW\n",
        "    obs: thinks[\n",
        "        department: chooses(mu in Mu, tau in Tau, wpp=partial_pooling[mu, tau](mu_scale, tau_scale)),\n",
        "        department: chooses(theta in Theta, wpp=normpdfjit(theta, mu, tau)),\n",
        "        department: chooses(sigma in Sigma, wpp=half_cauchy(sigma, 5)),\n",
        "    ]\n",
        "    obs: observes_event(wpp=professor_arrival_likelihood(d, department.theta, department.sigma))\n",
        "    obs: knows(_sigma)\n",
        "    return obs[Pr[department.sigma == _sigma]]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "52eb8c0d",
      "metadata": {},
      "outputs": [],
      "source": [
        "# Cache for precomputed results\n",
        "cache = {}\n",
        "\n",
        "def compute_distributions(mu_scale, tau_scale, verbose=False):\n",
        "    \"\"\"Retrieve from cache or compute the JAX-based distributions.\"\"\"\n",
        "    key = (mu_scale, tau_scale)\n",
        "    if key in cache:\n",
        "        return cache[key]  # Use cached results\n",
        "\n",
        "    if verbose:\n",
        "        print(f\"{mu_scale=}, {tau_scale=}\")\n",
        "\n",
        "    # Cache results\n",
        "    cache[key] = (\n",
        "        partial_pooling(mu_scale=mu_scale, tau_scale=tau_scale).sum(axis=1), \n",
        "        partial_pooling(mu_scale=mu_scale, tau_scale=tau_scale).sum(axis=0), \n",
        "        {d: department_model_sigma(d, mu_scale, tau_scale) for d in Department}, \n",
        "        {d: department_model_theta(d, mu_scale, tau_scale) for d in Department},\n",
        "    )\n",
        "    return cache[key]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "0d9a6e0d",
      "metadata": {},
      "outputs": [],
      "source": [
        "def plot_model(mu_scale=1, tau_scale=1, figsize=(10, 8)):\n",
        "    posterior_mu, posterior_tau, sigma_posteriors, theta_posteriors = compute_distributions(mu_scale, tau_scale)\n",
        "\n",
        "    fig, axs = plt.subplots(4, 1, figsize=figsize)\n",
        "\n",
        "    ax = axs[0]\n",
        "    ax.axvline(0, color=\"black\", linestyle=\"-\")\n",
        "    ax.plot(Mu, posterior_mu, label=r\"$P(\\mu \\mid t)$\")\n",
        "    mu_expectation = jnp.dot(Mu, posterior_mu)\n",
        "    ax.axvline(\n",
        "        mu_expectation, \n",
        "        color='red', \n",
        "        linestyle='--', \n",
        "        label=r\"$\\operatorname{E}\" + rf\"[\\mu \\mid t]={mu_expectation:0.2f}$\")\n",
        "    _ = ax.set_title(r\"Posterior of $\\mu$\")\n",
        "    _ = ax.legend(bbox_to_anchor=(0.9, 0.5), loc='center left')\n",
        "\n",
        "    ax = axs[1]\n",
        "    ax.axvline(0, color=\"black\", linestyle=\"-\")\n",
        "    ax.plot(Tau, posterior_tau, label=r\"$P(\\tau \\mid t)$\")\n",
        "    tau_expectation = jnp.dot(Tau, posterior_tau)\n",
        "    ax.axvline(\n",
        "        tau_expectation, \n",
        "        color='red', \n",
        "        linestyle='--', \n",
        "        label=r\"$\\operatorname{E}\" + rf\"[\\tau \\mid t]={tau_expectation:0.2f}$\")\n",
        "    _ = ax.set_title(r\"Posterior of $\\tau$\")\n",
        "    _ = ax.legend(bbox_to_anchor=(0.9, 0.5), loc='center left')\n",
        "\n",
        "    ax = axs[2]\n",
        "    ax.axvline(0, color=\"black\", linestyle=\"-\")\n",
        "    for d in Department:\n",
        "        department_name = Department._member_names_[d]\n",
        "        department_abbrev = department_name[0]\n",
        "        sigma_posterior = sigma_posteriors[d]\n",
        "        sigma_expectation = jnp.dot(Sigma, sigma_posterior)\n",
        "        ax.plot(\n",
        "            Sigma, \n",
        "            sigma_posterior, \n",
        "            label=(\n",
        "                rf\"$P(\\sigma_{department_abbrev} \\mid t),~ \" \n",
        "                + r\"\\operatorname{E}\" \n",
        "                + rf\"[\\sigma_{department_abbrev} \\mid t]={sigma_expectation:0.2f}$\"))\n",
        "    _ = ax.set_title(r\"Posterior of $\\sigma_d$\")\n",
        "    _ = ax.legend(bbox_to_anchor=(0.9, 0.5), loc='center left')\n",
        "\n",
        "    ax = axs[3]\n",
        "    ax.axvline(0, color=\"black\", linestyle=\"-\")\n",
        "    for d in Department:\n",
        "        department_name = Department._member_names_[d]\n",
        "        department_abbrev = department_name[0]\n",
        "        theta_posterior = theta_posteriors[d]\n",
        "        theta_expectation = jnp.dot(Theta, theta_posterior)\n",
        "        ax.plot(\n",
        "            Theta, \n",
        "            theta_posterior, \n",
        "            label=(\n",
        "                rf\"$P(\\theta_{department_abbrev} \\mid t),~ \" \n",
        "                + r\"\\operatorname{E}\" \n",
        "                + rf\"[\\theta_{department_abbrev} \\mid t]={theta_expectation:0.2f}$\"))\n",
        "    _ = ax.set_title(r\"Posterior of $\\theta_d$\")\n",
        "    _ = ax.legend(bbox_to_anchor=(0.9, 0.5), loc='center left')\n",
        "\n",
        "    _ = plt.suptitle(f\"mu_scale = {mu_scale}, tau_scale = {tau_scale}\", y=1)\n",
        "    plt.tight_layout()\n",
        "    plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "934c8597",
      "metadata": {},
      "outputs": [],
      "source": [
        "import ipywidgets as widgets\n",
        "from ipywidgets import interactive\n",
        "\n",
        "for mu_scale in range(1, 20, 3):\n",
        "    for tau_scale in range(1, 20, 3):\n",
        "        compute_distributions(mu_scale, tau_scale, verbose=True)\n",
        "\n",
        "def plot_model_widget(mu_scale=1, tau_scale=1):\n",
        "    plot_model(mu_scale=mu_scale, tau_scale=tau_scale, figsize=(9, 7))\n",
        "\n",
        "interactive_plot = interactive(\n",
        "    plot_model_widget, \n",
        "    mu_scale=widgets.IntSlider(min=1, max=19, step=3, value=1),\n",
        "    tau_scale=widgets.IntSlider(min=1, max=19, step=3, value=1),\n",
        ")\n",
        "output = interactive_plot.children[-1]\n",
        "output.layout.height = '700px'\n",
        "interactive_plot"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": ".venv",
      "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.13.2"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 5
}