{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Optimisation: Maximising a log-likelihood\n", "\n", "As well as minimising error functions, PINTS optimisation can be used to find the maximum of a loglikelihood (or of any [pints.LogPDF object](https://pints.readthedocs.io/en/latest/log_pdfs.html#pints.LogPDF)).\n", "\n", "Following on from the [first example](optimisation-first-example.ipynb), we can define an inference problem using the [logistic model](https://pints.readthedocs.io/en/latest/toy/logistic_model.html#module-pints.toy):" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzs3Xd8XMW1wPHfbFNb9WZJlmRblruxjSvFgMEYCD2U0CEhEAJJID0hyXskvCSkEeAFAiQQmmmP0INpBmyqjXuXbMmSrN610kraOu+Pu7vWWsVyWcuWzvfzyUfau7v3zsrhnp0zM2eU1hohhBAjl2moGyCEEGJoSSAQQogRTgKBEEKMcBIIhBBihJNAIIQQI5wEAiGEGOEkEAghxAgngUAIIUY4CQRCCDHCWYa6AYORlpamx4wZM9TNEEKIY8ratWsbtdbp+3vdMREIxowZw5o1a4a6GUIIcUxRSpUP5nWSGhJCiBFOAoEQQoxwEgiEEGKEk0AghBAjnAQCIYQY4SQQCCHECCeBQAghRjgJBEKIEcPR7eHldZXIFr3hJBAIIUaMNzfW8IMXN1LS4BzqphxVJBAIIUaMpg4XAKUNHUPckqOLBAIhxIjR3OkGoLRRegQ9SSAQQowYrZ0eQHoE+5JAIIQYMVoCPYLdQ9wj0FofVQPWEgiEECNGizOQGhriweJHVpZy1n0rh7QNPUkgEEKMGC2B1FCT001b4PehsK68heK6Djrd3iFrQ08SCIQQI0ZLp5vRyTEAlDYO3ThBRXMnAFUtXUPWhp4kEAghRgSPz097t5fZ+cnA4MYJtNb8Y2Up9Y7ug7pmdWsXFU2dvc65JxAIKiUQCCFGkuK69v0OkDpdXnz+yAyiBmcMHTc6CbNJhY0TuLw+nvh0Ny6vL+w9Va1d/Pat7by+sfqAr+fx+bnqH1/wnefWhR1v6fTgdBvXqWyVQCCEGCF21DpY8teVvLO1tt/XeHx+Tv3TR/z9o10RaUNrYMZQenwUeSmxYT2CtzbXcNcb2/hsV1PYe5o6jPcEZxsdiOdWV1DW1ElxXTv+HsEtmBYCSQ0JIUaQssBNd9mW/gPB9hoHjR0uXttw4N++B6M5MGMoJdbGuLQ4SnqsJfhgRwMAjYGVx/u+J/hzsDpcXu5/fyc2s4luj5/qtr03/GAgMJsUVdIjEEIcaT94ccNBpTkOVU2bkWP/cEc9Hp+/z9esLW8BYGd9RyhwHE7BGUNJsVbGpsVR1uTE79d4fX5WFNUDxmyinpoOMhA8uqKEJqebH501ASCstlFwfGB6TiKVLZ19vv9Ik0AgxDDj9fn5zRvbet1MPT4/L6+r4g/LduDt52Y8kOXb6/jdW9sPqk21gUDg6PbyZVlzn69ZU95CfLQFgPe31x3UdQYSTA0lx9kYl26n2+OnxtHNuopWHN3GNM59b/jNTqOH0OIc/FTTekc3//h4N+cel8VXjx8NQEn93t7HnuZO0uxRFGbYJTUkhIiMrdUOHv90N+9uC0/DBNMeVa1dvLftwG+0j32ym0dXlvZKnwxGTVs3GfFR2Cymfq+9rryF0yZmMGlUPO8eRPv2J1hnKCXWxti0OMAoNbF8Rx0WkyI51trrswV7BE3OwX/mv7xbjMfn58dLJpIaZyMxxhqWhqpo7iQvJYac5Bjq2129BqiHggQCIYaZTZWtANQ5wm9e9T0eP/7p7gM6p8vrC6VuVu/u+xv9QGrbuhmTFsfJ49N4f3tdr9lD1a1d1LR1MzsviTOnZLKmrDm0Cvhwae30EGUxEWMzU5BuBILdjU4+3FHP3DEp5KbEhgaHg5pDg8WD6xG8+OUeXlizhxsXjmVMWhxKKcZn2NlVHx4IclNiGZ0cC0B168FNTT2cJBAIMcxsrGwDoG6fue/17UYguHBmNl+WtYQCxmCsr2jF5TXSSV+UNu3n1b3VOLrISoxm8eRM9jR3UVwXvpgrGGRm56eweHImfg0f7Kg/4OsMpMXpJiXOBhgzh+JsZlYWN1Jc18HpkzJIjbP1kRoyHrd2uvc7rXVNWTO/eHUzCwvT+PGSiaHjBelxoTECj89PdWsXeSmx5CQZC9uOhvRQRAOBUqpMKbVZKbVBKbUmcCxFKfWeUmpn4GdyJNsgxEgTvMHX79MjaAgEgu8sGk+czcy/Pi0b9Dk/L2nCpGBWXhKflxxYINBaU9fmYlRiNIsnZwC9xwDWlrcQYzUzKSue6TmJZCZE9ZtC+su7Rdz//s4Br7muooW/fRD+mpZON0mxRiBQSjEu3c7yHcY1Tp+cQUpcVGi/gqBgasivwdEV3ivw+vx0e3z4/Zqq1i5ueWYtOUkx/O3K47GY995aC9LtNHa4aOv0UNPajV8T6BEEAkHr0A8YH4kewSKt9Uyt9ZzA458By7XWhcDywGMhxGHQ4fKyM5CGqG/ft0fQjVIwJi2Oy+bk8uam6l69hv58XtrEtJxElkwZxc76jgMaJ2h2unH7/IxKiCYjIZoZuUm9bvLrKlqYmZuE1WzCZFIsnpzJyp0NdHvC8+d+v+bpL8p5dUPVgNd8+KMS/vxuMV3uve9v6fSQHGsNPR6bFofWkJ8ay7i0ONLsNhqd7rC0VbPTjcWkgN4zis594BMm/eptxt35Fif/4QNcHj//vH4uiT2uAUYgAChp7AhNHc1LiWVUYjQmNQJ6BP24EHgy8PuTwEVD0AYhhqUtVW1obaQj6hyusJtafbuLlFgbVrOJG04cg9evWbqqYr/n7HL7WF/RwgnjUlkwLgWAVaWDHycITh3NSowG4MzJGWzY0xoKVJ1uL1urHaHSDwCLp2TS6fb16n2UNnbQ2umhorkTt7fvmU8enz/0vj09pme2dLpJDqSGAMYFxgkWTcxAKUWq3Ybb6w+t+gUjEAQHlnsuKuv2+Ciqa+fUCencsbiQb59awNKb5jM+w96rPQWBYyX1ewNBbkosVrOJUQnRR0WZCUuEz6+Bd5VSGnhEa/0okKm1rgHQWtcopTIi3AYhRoxgWujMKaN4eEUJ7S4vCdHGN9R6h4v0+CjA6BXMzE0aVL5/bXkLHp9mQUEq03ISibOZ+aK0iXOPyxpUm4JTR0clGqmQJVNH8ed3i3lg+U7+56LpbKpsw+fXYYHgxIJU4mxm3ttex6JJe28RwbEEn19T0exkfEZ8r+tt3NNKu8uYDlrR1MmETOM1LU53WI+gMPDeMwLpqpRYGxa8NLe0YE+y4XK7iHE1MDc5la6GerpriyE2Cfxemls7OE6VcF2unzMKveD3gasGdvnA7wftA+0Hv488n5eLLOuI2r6LaO3jcksdWaXNgOb6qGJMVcDqL0FrQBvvC/1Pw6xrIDZlUH/rgxXpQHCS1ro6cLN/Tym1Y7BvVErdDNwMkJeXF6n2CTGsbKxsIycphslZxk2u3tEdCgQN7d1kJESHXjs9J5GX11Xh92tMgfRHXz4vbcRiUswdk4LVbGLOmJR+A4jfr9nV0BG6+QLUOMJ7BBMy47n5lHE8urKUgnQ7nYFv4LPykkLvibKYOXV8Iuu370KfFo9yO8HtpG3LZpaY9hCDC9eqEkizgKcTPF3g7QZPF9EV9TxgbSQKDxOXR8FqE9rr4ilvPVlFJngA8Hn4itfFzgQXlhe94HNzqd/DpdHAw4E2AF9GA+WBB2/v/ZzZwOtRwKeB/w3ADNxnAUqMx+dbgNeN378VfNFbA5xgwtnHdiDQWlcHftYrpV4B5gF1SqmsQG8gC+hzakCg9/AowJw5c46erXyEOIptqmxlRm4iGfHGTbfO4Qp9a65vd1HY4wY9LTuRpz4vp6zJybj03imNoM9LmjhudCL2KON2cUJBKvcs20Fjh4s0e1TYa1/dUMUP/28jy39wqnFOn5f2hiommStJa1wNlS3Q2czP7C3MztxOy9v/YIy1m5fjOkl6+o/gckC3A1ztPOQLjEPcv/f8NwM3B7M7a3teWYE1FqzRZHabSLZZcXotmLriIC4ZrzmGOp2MPT6VtIwksEShzFasZhuYbWC2Uuf08/Tqai48Pp/CrGSq2708uLKcS+eO4ZnV1Zx9XA5nTssBk4XV5Q4e/riMuy+eSU6KHUxmMFlAmY3flWnvT2Xmrjd3UNrcRYzNSly0jXu/djyg+PvK3TzxRSWf/vR0LBaL8Xow3osyHltjD+z/BAchYoFAKRUHmLTW7YHflwC/wYiF1wP3BH6+Fqk2CDGSNDvd7Gnu4ur5+WQmGDfoYB7e79c0tLvIiN97456akwDAlmpHv4Ggw+VlY2Ubt5w6LnRswbhUzPjYuGUrZ+R4wVEFjhpor2H05u08Y6km7an/Ap9x078Vza1W4Km95zUBS8xRNFviaPLGYLMngz0bUsdDdAJExdOuo/nLihpOnT6WRdPH4vBHc+Oz27j8xEm8sLGJ48Zl818XzzZulGYbKEVbp4f5d7/LdxaN573t9YxKiOJfX5/HnoYOvvmXFfz1hBmMnTW6z8/qbe3ib59/wOjc6RTOy2NXcQNLP1zNRTNP4K21q0m253HmtCkAbGws5QN/Cn+dtgRirH2eryf7aDOflpQQazNzwYxsSDKyHEmjNHX+dmr9CYyOi/wNvz+R7BFkAq8opYLXeVZr/bZS6kvgRaXUjUAFcFkE2yDEMe/LsmbufnMbL9x8AjE2c7+vC44PHDc6MZQCCi4qa+l04/Xr0BgBGDlym9nE1qo24+a0r85mitZ/yQWs5PKOL+CVemgpZ0ZbBUVR1VjeDh+s1WYbOb4k6lQCNeZsEiYsBHsGT2xsp9Efz48uOhFiUyEmBWKSUbZYvI5u7ly6jtsXFzKmMD3sfPHA6u0fs63NwqKpJ/DF1lq+1H5+Mv0ErDXFrG/1QUz47PPPSxvxa1g4IZ2iuvbQQq69dYZs9Cc1MJC8b32hlDgbKXG20MpkgFpHNzFWMwnRg7uFFmTE4fNr2ru95KXsveH3XEswOjmWL0qbeHhFCQ9fM5toa///1odbxAKB1roUmNHH8SbgjEhdV4ijRZfbh8lk5LsPxb/XVrKpso3djU6mZCf0+7pNlW0oZeT+7VEW7FGW0PTQ4GKyYMoIwGYxMWmUneqKnVC8Bxq2Q0MxNBZB407obmU2MNsGerOChGxIykfln8Rb2kyJO4nvf3WRcTwhm63NJs77m5EwP8meytLzFwDw1LqPmDQqHsbN7tXmzIRoXvr2if1+pkWT0nl4RSltXR7WlrdgM5uYnpNIQUYcr2+oRmtN4MsmACt3NmKPsjAzN4n81Dg+LGrA79ehVcopAwSCaKsZe5QlNDU2GBBSg4HAGR4IRiVGh117IAU9elxhgSC0lqALv19z1+tb2VHbztryFk4anzaocx8OkR4sFmLEuvqfXzAu3c6fL+v1fWjQtNZ8vLMRgFpH134CQSvj0uKIDwwOZyREhRaVNbS7MOEnz1cOGz6Cmo1Qs5EXWjYR4++AZwMnsWdC2gSY9lV0SgF3ruzGk5jPn2+6ACx7exP1H5dy/3+2c17iCaFxh0+/NEZDF01MZ8Oe1tDU1dq2bk6bcHCTAxdNzODBD0v4ZGcja8pbmJaTQLTVzLg0O45uL01Od2icQmvNyuIGTihIxWo2kZsSi9vrp669OzT1M3mAQACE3fCbOlyYTYqEaCvJcbawkhd1bd2h9Ntg9Ey95fbRI6hs6eLNzTXsqG0H4LOSxiMaCKTEhBAR4PNrtlQ5eGdrbb9llwejrKkzVLN+oJo0Wms2VrYxY/TemTfj4jzkNKyE93/NpHeuYmPUTUx/7Sx49RZY+wT4vezJOYdfer5O3aWvwk/L4EfFcMObcN5f2TbmWp5rnczs2fPDggDARbNysJoVz3+5J3Tsk12NTMi0s2hSBi2dHmod3bS7vHS6faEZQwdqZm4SiTFW3t5ay+bKNuaMMWbP9JybH1Te1EllSxcLC40baH7ghlvR1BnanSwpbuB8fqrdFqo31Ox0kxxrw2RSpMRae6WGRiUM/jPZoyyh1/cMBNFWM2n2KMqbOrnvvWImjYpnVl4Sn+468DIeh0J6BEJEQE1bF26fH7fPz7ryFuaPSz2o83yysyH0e3A+fl921nfgam/mXFsFLHsOdn/MP+u3Gk9+ZoHYQl7xncwVF1+MLW+OMShrMtO1p5VnHvyUk/Qkztkn3/7mphrMJsVZU0f1ul6aPYolU0fx73WV/Pgso67O6t3NXDU/j6mBXsvWKgd5qcZNb9RBBgKL2cQpE9L5z6Zq/JrQWoNxweqhjc7Q3/bjwN9qYWCsIZiCKW/upLnTWCEcHzXwLS81Liq0R0CT0x0aN0iJiwoVoNNaU+9wkXmAn6kgI44uj4/EfQaXc5JjeGNTNW6vn0evnc3mqjYe/HAXjm5PaOpvpEmPQIgIKO+xYfmHRQ1hz7m8vl6lE8DIE3+wI7z0wie7GslJiiErMTq0QjfE74fKtbiX/x7vo4tZH30zZ2z8Pqx9EuwZrMi5mWt9v0L/rIKHJjzGn8w3YZt9NaRPDExPhImj4rGYFJur2sJOrbXmP5tqOGl8WqhQ276umpdHa6eHd7bWsra8BZfXz8LCNCaOSkAp2Fbj6LWq+GCcPimdYL23YCDISYohymIK6xG8u62OvJRYxgSCT05yDCZl1P9vDdQZ2l9Ov2fhueYeRepS4qw43ca/W0unJ1Qy40BcMz+fW04t6HV8dFIMbq+fGblG5dUTC9Lwa1h9AKu3D5X0CISIgLImo9rk2LQ4Piqq52fnTAo9d9vS9TR2uHj1tpPC3vO7t7bz1uYalt2+kEmjEvD6/HxW0sS507Mormunpq0LvC4o/Qh2/AeKloGzHgsKt38clTNuI3/OuZAzGyxR7Py4lI9LtuPw2mhod5HeR0472mqmMDOeLdWOsONbqhxUNHfynUXj+/2MJ4xLJT81lmdXVXB8fjIWk2Le2FTsURbGpMaxrdoRyqNnHuBNs6dTCtONGkmpcaHxAJPJKBoXrPNf5+jm012N3LZofOhmbzWbyE6KobzJKEeRsp+0EBipoWanG79f0+x0h8ZkgqUpWjs9oUBxoIHgnOl9r8QOFp/78ZKJKKWYlZdElMXEpyWNLJ6SeUDXOFgSCISIgPKmTmwWE1+bm8s9y3ZQ09ZFVmIMW6vbQpU3i2rbmTjKGGh1urws316H1vDX94p55No5bKpqo73by8KCRPKaPqGg7l3401pj0ZXNDuMXs8I0lzvWpPLNs+Zy2z437eDNt669m/r2btLtfQ9uTstO4IMd9WEzcN7cVI3V3HdaKMhkUnxtbi5/fLuI8qZOjs9LDi06m5KVwKaqViYFVjgfSiBItUdxwYxsxqWFr3UYlx7HlkBP5rUNVfg1XDwrJ+w1+amxVDQb/xYDTR3teS2vX+Po9tDU4SIt2CMIvLfZ6Q7NxMo4hM/U09Xz88lPjeOk8UaKK9pqZu6YlAOu8nooJDUkRASUNTrJT4nljECdnI8C6aFHV5YSazNjNqmwCprLd9TT7TFSK+9srWPznlaK1n3Mf1ue5Jx3FnFr9c9Z4FmNnnw+XP0S/KSU3ac/xI3rxjF7slH0bF+hQODopr7d1e+Na/roRJqcbmoDNzitNW9uquHk8Wm9Kmnu69LZo7GYFLWO7rBZLlOyE9jT3MXOug7S7MbOZIfi/itmcfviwrBjBel29jR34vL6eHldFTNzk3otjMtLiQ2lhpL381lg71qCOocLR7eXlDgjeAZTRM09/k4HO+6xr7zUWK6anxeWtjpxfCo7atsPaje4gyGBQIgIKG/qJD81jvEZdnKSYvhwRz17mjt5c1MNV8/P4+Txaby+oRp/IPn9xsZqMhOi+Nsl47k55gMSn1rElRuu5RrLckxjF/Lucfcxx/V3HEvuh8IzwRLF+ooWvH7NT8+e2GetoOAq4jqHi3pH+KrinqZmJwKwObChzYY9rVS1dnHucX0sMut1jWjODKQvTi4MDwQAK3c2HNL4wEAK0uPwa3h7Sy07atu55PicXq/JS4mjyemmurW737GOnlLtxmt21hvTOFPswTGCQCDodFPbZpTz7u/veTicWGD8LY9Ur0ACgRCHmd+vKWtyMiY1FqUUp01M59Ndjfx9RQkmBd84eSwXzcqmqrWLNeUttHV5qCxaz0MJT5H44DTu1P/E4fLzS883+NvsZXD5k3jGn40HC9Vte0sWlzY4MZsU+alxfbYjI5CfL23ooMvj6/fGNTkrHpOCZVtqeeaLcv70ThE2syl0g9+f2xcXcu2CfGaMTgwdm5plBIL2bu9h++a8r+AirXvfK8ZqVpzXR+AKzhzqcHkHlxoK9AB2BnZQC/YQgmMELYHUUGpcFFZz5G6f03MSiY+28FlJY8Su0ZOMEQhxmNW1d+Py+hmTtrfe/dJVFTy7qoJLZ48mKzGGJVNGEWPdzJaPXyWv/f9YZv0Mf4sNjruc7hnXc8PTrTS6PSydbNT4yUoybqa1bd1MDtxkSxo6yE+J7TftEmuzEB9tCc0IyuhnAVSszcLEUQm8sr6KV9ZXEW01cfWCvF7THPszaVQCd180LexYenyUsdFLhztiPYLgPgHlTZ0smZIZttdAUM9VvINKDe3bIwicMynGilKEUmijEiPXGwAwmxQLxqUesfUEEgiEOMzKGo2po2MC39RPHJ+KzWzC7fPzrVPGgd9PXMl/eCf2f8grLabFlMKjlqu46Y7fgD2daOAnZ+/h7ytKQtMlgzfTnlNISxo6BqwaCsY4wdbAjKCe5SX29cg1s6ls6WRsehyZ8dEDlqUeDKUUU7ITWVnccEgDxQOJi7KEptV+tY+0EBBaxwD7X1Xc8zXF+/QILGYTiTFWWpxGaig40yeSTixI5b1tdewJbHYfSRIIhDjMygNTR/MDN6FYm4XzZ2SjtJ/Cpg/gpd9Bww7S7Pn8pOMmXvWdzDdOnYSy7y26dvncXC6fmxt6nG6PwqSMhWpgrFwua+wM27SlL5kJUaHCa+kD5LTzUmPDbpqHw5SsBFYWR26MAGB8hrGfQX9/h8QYK4kxVtq6PIMKBDaLccMvazT+DXuOK6TEGoXn6ttdYZvoRMrpkzLodPuOSPE5CQRCHCCny8v3X9jA2PQ4rp6X3+sGWtbUiS0whx0ArfnLrHpY/ht4cZNRy+eSx7BOupD3fv8h7k4P588YeLcvi9lERvzeRWWVLZ24fX4K0vbTI+jRC4jk4GZfgiuMsxIj9+35zq9MxtHlGbCwX35qLJsq20gexDoCMHoBbV0elAqvVpocZ6OurZtmp/uA1xAcjPzUuF5TgiNFAoEQB+i+94t5d1sdZpPi0ZWlnDohnbvOnxoaEyhrdJKbEoPZpKB2C7z7C2MRWPIYuPgRmH4ZmMxYgWsW5PNFaRNTsvovJheUlRQdKjMRXEhVkNH3QHFQcBFZ8JvukXTW1FH84ZLpzBsbud21Jg/i75abEggEg+gRgDFOUNroJDnWZvwbBqTE2UKzeCKV7hoqEgiEOABbq9t4/NMyrpyXy+1nTOC51RU89slu7lm2g4evNcoslzU5mZbkhde/B+ueguhEOPsemHMjWMJvRj9cMnHQ185KjKYoUJ2ypN5IXey7yGpfwR5Buj1q0CWTDxdjQd3QbzMbLD436ECwz9qBoJRYGx2BvZAPtM7Q0U4CgRCD5PNr7nxlC8mxVn569iSSYm18/8wJdLi8PPV5mVGtMsbMvObX+Vn7C1DVCQu+Daf8+LDsOTsqIYaPihrQWlPS0EFKnK3PmTI9Bb+59jdjaCS4dPZoYm1mkgYxawh6rx0I6vm3PhKpoSNJ1hEIMUhLV5WzcU8rvzx3Slju+LI5o/H4NB99vALvo4v5jekfOBImwC2fwNm/P2wbj2cnRdPp9uHo9lLa4KQgfeC0EBCq9XOkxweOJuPS7Xzn9MJB94iCZSVS9+0R9BhjkEAgxAjU4nTzp7eLOHl8GhfODF+4NCktmt8m/4fzv7gCWsq4w30rRWc/BxmTD2sbgguzatu6KWnoCNv1qj+hHsEAU0dFuFR7P6mhQMoo2moiIWZ4JVMkEAgxCK9uqKLd5eXnX5kU/s2ybiv8YxFXdy3lLd88/jrhaV71n8zY/eTuD0ZwGub2GgdNTjfjBtEjSI+PIsZqDk1lFfuXsp8ewaiEwW9ReawYXmFNiAh5aW0l03ISQnV58Pth1cPw/l0QnUDHxU/x4xdtsK4di0mRnXT4v4EHp2F+ussoOzCYHkG01cyy2xdGrMzDcJTa3xhBIB043GYMgfQIhNivbdUOtlY7uPT40caB9jp45mJ45+cw/gz49ufYZ1zImVMzcfv85KbEYolAHZr0eGNR2YEEAoAxaXFHZFHScJGbHItS9FofEgwMwzGoSo9AiP14aW0lVrPiwpk5sHslvHQjuNrhvPtg9g0QSBNcNns0/9lUE7E0jNVsIj0+iuq2bmxm0xEpczAS5abE8sEPTwvtdBYUDATDsUcggUCIAbi9fl7dUMWZk9JJXnM/fPQ7Y7/f616DzClhr11YmE5Behyz8yJXfiArMYY6h4v81Mj0OoQhWNCuJ3uUhe8sGs/Z0/rfrOdYJYFAjEjvb6sztoa8feGAaZMPi+pxO1v5ddff4cPlxqrg8+6DqN5pGbNJ8d73Tz3kgm0DyUqMZsOewaeFxOGjlOJHZw1+AeCxRAKBGJHe3lpLaaOTXfUdTMvZW0e/rdPDX98vZlpOIosmprPyiy94Pfou0mpq4Jw/wrybQ6mgvkQyCMDe/PT+SksIcSAkEIgRaV1FC2DU7OkZCN7fXscTn5UBcKJpK3+3/hWr1Yq6+lUYe8pQNDVMdmDm0P5KSwhxICQQiBGntdNNaYNRqydYojmouL4dm9nE+6dXkvPJPdRYcjFf8yKx+UdHSiA3xQgEwU3vhTgcIh4IlFJmYA1QpbU+Tyk1FngeSAHWAddqrd2RbocQQev3tAJGhidYxTNoZ207/21/hbyPX4BxpzH68kDRuKPE4smZPH3jvLBejBCH6khMO7gd2N7j8R+Av2qtC4EW4MYj0AYhQtaXt2BSMH9sSqiKJwA+L1+t/D1Xu16AWdfC1S8dVUEAjH0JFham7/+FQhyAiAYCpdRo4Fzgn4HHCjgdeCnwkieBiyLZBiH2tX5PKxNHJTBjdBK7G514fX7wdON9/hrO833A6ryb4IIIXhrsAAAgAElEQVT/BfORrd8vxFCJdI/gPuAngD/wOBVo1Vp7A48rgb43GxUiAvx+zYaKVo7PS6Ig3Y7b56eqrh6euQTLzmX8l+d6Wuf/aMCZQUIMNxELBEqp84B6rfXanof7eKnu5/03K6XWKKXWNDQ0RKSNYuTZ1dBBu8vLrLxkCjLsJOAk6aXLYM8XfDHzDzzlO4sJmTIQK0aWSPYITgIuUEqVYQwOn47RQ0hSSgUHqUcD1X29WWv9qNZ6jtZ6Tnq65ETF4bGu3Jg2enxeEuPtbpbafou9ZRtc/hTvW04h2moiN0UqdYqRJWKBQGv9c631aK31GOAK4AOt9dXAh8ClgZddD7wWqTYIsa/1Fa0kxVoZG9tN4ouXMMFUxZN5v4VJ51JU1874DHvYPrVCjARDUazkp8APlFK7MMYMHhuCNogRal1FCydnm1BPXwRNO/lzyl282TUdgJ11HZIWEiPSEVlQprX+CPgo8HspMO9IXFeIntq6PNTW1/O0/17o2gVXPodzcyYlm2po6/RQ6+iWQCBGJClfKEaMzaVVPG77IxnOYrj8KRi/mPHpdtq6PHxe2gTARAkEYgSSEhNiZPC6yHvvZnLUTrov/CexE88BoCDDqNmzbEsNAIWZUsNHjDzSIxDDn98Hr3yLvNZV/CnqO8TOvDT01PhAIFi+vZ44m5mcJNnsRYw8EgjE8KY1vPVj2PoKD9tuYHdu+EL2rIRoYqxmOlxeCjPjh92m5EIMhgQCMbyt+COseQzPgu/xh/YlTBqVEPa0yaQYl27U9p8gaSExQkkgEMPXhmeNrSVnXMX2qT9Aa5jUR/nmYHpIZgyJkUoCgRieSj6E178L406D8+9nR51RbrqvOv7BbR8lEIiRSgKBGH7qtsKL10HaRGOaqMVGUW070VYT+am9t3g8uTCNcWlxHDf66Co5LcSRItNHxfDSUQ/Pfg2ssXD1i6H9BHbUOijMiO+zfMTxecl88KPTjnBDhTh6SI9AHHW01vzPm9vYWt12YG/0dMPzV4OzEa56HhJHh54qqm3vc3xACCE9AnEUqmzp4p+f7MZsUkzNHmS6Rmt443aoXA2XPQnZs0JPNXa4aOxwyz6/QvRDegTiqFNU2w5ARXPn4N/06X2w6XlY9AuYGr5WIHi+faeOCiEMEgjEUaeo7gADwc734P1fw9Svwik/7vX09hoH0PeMISGEBAJxFCoOBoKmTrTucwO7vZpK4N83QuY0uPDBPreYLKptJ81uIz0+KhLNFeKYJ4FAHHWCqZx2l5e2Lk//L3S1G4PDygRXPAO2vncWK6prl96AEAOQQCCOKh6fn5KGjlDZh37TQ1rDa7dBYxFc+i9IHtPny3x+TXFdOxMzZXxAiP5IIBBHlbJGJx6f5swpmcAAgeDzB2Hba+w67kd84JnCBzvqqHN093pZeZOTbo9fpo4KMQCZPiqOKjsCaaHFkzN5ZEVp34Gg/DN477/4xHoi16yaAavWAJASZ+PFb50Qqh0EPWYMZUkgEKI/0iMQR5XiunZMCqbnJJJmt7Fn30DQXgf/dwP+pHxu7fgG1yzI57XbTmLpN+djUnD946upaesCjDTT8h31KAWFGRIIhOiPBAJxVCmqbWdMWhzRVjO5KbHhPQKf15gh1O1gy8n/i0PHcsakTGbkJnHS+DSe+Po82ro8XPfYat7aXMM593/MS2sruWhmDjE289B9KCGOchIIxFGluG5vKYi8fQPBij9A2cdw3r2scmYDML1HobhpOYn847o5lDd3cuvSdfj8mn9eN4d7L59xRD+DEMcaGSMQR41Ot5fy5k4umpUDGIHgzU01eHx+rGUrYOWfYObVMPMqNj23npykGNLs4WsDTihI5Ykb5lLa6OTyObnYLPJdR4j9OaBAoJQyAXattSNC7REj2K76DrSGiYF9AXJTYvH5NXVV5Yx++SZInwhf+RMAmypb+y0bfeL4NE4cn3bE2i3EsW6/X5eUUs8qpRKUUnHANqBIKdV7Hb8Qhyg4w2dij9SQCT/2/9wCbqdRTM4WR2unm/KmzrC0kBDi4A2m3zwl0AO4CHgLyAOujWirxIhUVNuOzbJ385i8lFi+bX6dpLovjJ5AxiQANlcZ5alnjE4asrYKMZwMJhBYlVJWjEDwmtbaA+ynAIwQB66orp3CDHto85jMto183/ISW1OXGGMDAZsqjUAwLUd6BEIcDoMJBI8AZUAcsFIplQ/IGIE47Ip71gTqasX88k3Um9J5PPF7YcXkNu5pZWxaHIkx1iFqqRDDy34Dgdb6Aa11jtb6K9pQDiza3/uUUtFKqdVKqY1Kqa1KqV8Hjo9VSq1SSu1USr2glLIdhs8hjnGNHS7qHC5joFhrePMOaK/mHxm/oLgtvKLo5qo2pktvQIjDZjCDxZlKqceUUssCj6cA1w/i3C7gdK31DGAmcLZSagHwB+CvWutCoAW48aBbL4aNt7fUAsZG8mxYCltfgUW/wJs1O2wtQX17NzVt3bLRvBCH0WBSQ08A7wDZgcfFwB37e1Og99AReGgN/E8DpwMvBY4/iTH2IEa41zZUUZhhZ0pUI7z1ExizEE66g7yUWNq6PLR1GuWoNwfGB2bkykCxEIfLYAJBmtb6RcAPoLX2Ar7BnFwpZVZKbQDqgfeAEqA1cA6ASiDngFsthpXKlk6+LGvhqzMyUC/fDGYrXPwwmEzkphh7DOxpMXoFGyvbMCmYmi1lpYU4XAYTCJxKqVQCM4UC6Z22wZxca+3TWs8ERgPzgMl9vayv9yqlblZKrVFKrWloaBjM5cQx6vWN1QBc5XoBqtbA+fdB4mjAmEIKRjnqTreX1bubKMyIJ9Ymi+KFOFwG81/TD4DXgQKl1KdAOnDpgVxEa92qlPoIWAAkKaUsgV7BaKC6n/c8CjwKMGfOHJmuOoy9tr6aq7KqSfzyfph5DUy9OPRcbkoMAHe9vpVmpxuvX3P1/LyhaqoQw9J+A4HWep1S6lRgIqCAosBaggEppdIBTyAIxACLMQaKP8QIJM9jDDq/dgjtF8e47TUOKuvqeSn5PkjKg3PuCXs+PtrK6ZMycHR5uHT2aOaNTeGEgtQhaq0Qw9N+A4FS6rp9Dh2vlEJr/dR+3poFPKmUMmOkoF7UWr+plNoGPK+U+h9gPfDYwTRcDA+vbqjiLuvT2Lur4cq3Iar3vgGP3zB3CFomxMgxmNRQz/8Ko4EzgHXAgIFAa70JmNXH8VKM8QIxwvn9mra1r3CZ+SM4+UeQN3+omyTEiDSY1NB3ez5WSiUCT0esRWJE0FrzwOuf8GPPQ7QmTyHp1J8OdZOEGLEOplh7J1B4uBsiRg6tNb9+fSvT1v6KBJOLxKv+BRZZYC7EUBnMGMEb7J3iaQKmAC9GslHi2Nfe7eHTXU2cNTUT1aNOkN+v+eVrW/CteZLF1vXoJb9HBaqKCiGGxmDGCP7c43cvUK61roxQe8Qw8fzqPfz2re1846Sx/Oq8ySil8Ps1d76ymU/XrGF5zFJ0/imo+bcMdVOFGPEGM0aw4kg0RAwvGytbUQoe/3Q3sTYzPzhzAne+spkXvyxnZfoTWN0W1IUPgUm2khRiqPUbCJRS7fS96ldhlBKSNf6iX1uq2jhzciap9ij+9uEuVhQ3sLmqjScnrmJ0+Ua46GFIyh3qZgohGCAQaK17T+gWYhAc3R7Kmjq5dPZobj1tPN0eH6+sr+Ku+YpTtjwCk86DGVcMdTOFEAGDLtiilMrAWEcAgNa6IiItEse8LVV7dxAzmRR/vmwGt5ycy4Q3LkRFJcD594dtNCOEGFqD2Y/gAqXUTmA3sAJjt7JlEW6XOIYFA0Fw8xizSTGx6O+o2s1GEIhLG8rmCSH2MZiRursxisUVa63HYqws/jSirRLHtM1VDrITo0m1RxkHqtbCx/fCjCth8nlD2zghRC+DCQQerXUTYFJKmbTWH2LsOCZEn7ZUte3dWN7TBa/cAvGj4Ox7Bn6jEGJIDGaMoFUpZQc+BpYqpeox1hMI0Yuj28PuRieXHB/Yb2j5b6CxGK59FWJkVzEhjkb99giUUn9TSp0EXIhRVuIO4G2MXcbOPzLNE8earVUOwBgoZvfH8MVDMPcmKFg0xC0TQvRnoB7BToxVxVnAC8BzWusnj0irxDFrc1UrANPTzfDUrZAyDs789RC3SggxkH57BFrr+7XWJwCnAs3Av5RS25VSv1JKTThiLRTHlNBA8Sd3QVulsXDMFjfUzRJCDGC/g8Va63Kt9R+01rOAq4CvAtsj3jJxTNpS1caVydth3VNw4vdkjwEhjgGDWUdgVUqdr5RairF+oBi4JOItE8ccR7eHlsZavtF0L2RMhUV3DnWThBCDMFCtoTOBK4FzgdUYewzfrLV2HqG2iWPM1ioHd1v/RYzPARe/CpaooW6SEGIQBhosvhN4FviR1rr5CLVHHMMca57nfPMXOE/8OXFZxw11c4QQgzRQ0TmZ7ycGzd9axQnbf8dO60QKF/1oqJsjhDgAUgxeHDqtaX3hW1i1m/JT/wrmQdcyFEIcBSQQiEO35jFSaj7mftO1LFywYKhbI4Q4QBIIxKFpKsH/zi/52D8d5t5ElMU81C0SQhwgCQRiUJ5fXcHDK0rCD/q88PLNuLWFn3hu5uoF+UPTOCHEIZFkrtgvr8/Pn98twtHl5Yq5uSTF2ownPrkXqtZwt/n7TJowidyU2KFtqBDioEiPQOzXqt3NNHa4cfv8vLGx2jhYtRY+uofq3PNY6pzLNdIbEOKYJYFA7Nebm2qItZkZn2HnpXVV4O6El29Gx4/ih85ryE2J4bSJGUPdTCHEQYpYIFBK5SqlPgwUqtuqlLo9cDxFKfWeUmpn4GdypNogDp3H5+ftLTUsnpzJFXNz2binldbXfwZNu1g363d8Xu3nu6cXYjbJHsRCHKsi2SPwAj/UWk/G2OryNqXUFOBnwHKtdSGwPPBYHKU+K2mipdPDecdlceHMHBab15O05Un0gtv41aZU8lNj+eqsnKFuphDiEEQsEGita7TW6wK/t2NULM3B2OgmuK/Bk8BFkWqDOHT/2VRNfJSFUyakk67auDfqHxSrMbyd+S221Ti4/YxCLGbJMApxLDsis4aUUmOAWcAqIFNrXQNGsFBK9ZlcVkrdDNwMkJeXdySaKfbh9vp5e0stZ07JJNpighduJY4ubuu+k7rXdjAuPY4LZmQPdTOFEIco4l/lAvsd/xu4Q2vtGOz7tNaPaq3naK3npKenR66Bol+f7GrA0e3lvBlZsPofsOs9/GfeTV3UGBzdXukNCDFMRPS/YqWUFSMILNVavxw4XKeUygo8nwXUR7IN4uC9ubGGhGgLC+Pr4N1fQuESrAtu5roTxjA7P5nzjpPegBDDQcRSQ0opBTwGbNda39vjqdeB64F7Aj9fi1QbxMHrcvt4Z2stF09LwfrKNyEmCS58CJTiR2dNHOrmCSEOo0iOEZwEXAtsVkptCBy7EyMAvKiUuhGoAC6LYBvEQXp3Wy1Ot4/veh6HxmK49hWwS4pOiOEoYoFAa/0J0N/k8jMidV1xePx7XRVXx28gc+dzcNIdUCDbUwgxXEmtIdFLvaObsp1beTT2Icg+Hk7/5VA3SQgRQTLlQ/TyxrpyHrD+LzazgksfB7N1qJskhIgg6RGIXpI+/x0zTSVw0VOQMnaomyOEiDDpEYgwFZ//m0tcr7Ij92sw5cKhbo4Q4giQQHAM2d3o5N73itFaR+YCLeWkv387W/UYMi79c2SuIYQ46kggOIa8tHYPDyzfSVVr1+E/uacb/wvX4fX5eDbvblISEw7/NYQQRyUJBMeQknonABVNnYf/5G//DFPtBr7v/jYXLDrp8J9fCHHUkkBwDClp6ACg7HAHgg3Pwdp/8YzlYppHn8G8sSmH9/xCiKOaBIJjhNfnp6zJ6BGUNzsP34lrNsGb36cxdS7/3fFVvn3aeIzqIEKIkUICwTFiT0sXHp8xSHzYUkPOJnj+anRsCt/xfJdxGYmcMUm2nBRipJFAcIwoqTfSQilxNsoPRyDweeGlr0NHHWsXPMAX9RZuObUAk2w5KcSII4HgGBEcHzh1QjoVzZ39TiGtbu3i7PtWUt60n/TR8rtg9wo4717+tDmO7MRoLpgpZaWFGIkkEBwjSho6SLNHMS0nkQ6Xl2anu8/XLd9Rz47adtaWt/R/sg3PwWf/C3O/SdXYS1i1u5mrF+RjlU1mhBiR5L/8Y0RJg5OC9DjyU2IBKG/uOz20enczADVt3X2fqGIVvPE9GLMQzr6Hd7fWAnDOtFGHv9FCiGOCBIJjgNaaXfUdFGTYyU81AkFfA8Zaa1bvbgKMFFEvrRXwwtWQOBoufwrMVt7eUsuETDvj0u0R/QxCiKOXBIJjQLPTTVuXh4J0O7kpsShFnwPGFc2d1DlcQB89Alc7PHcVeN1w5QsQm0JTh4svy5o5a6r0BoQYyaT66DGgpMEY+C1IjyPaamZUQnSfawlWlRppofzU2PAegc8L/3cD1G+Dq1+E9AkALN9ej18jgUCIEU56BMeA4IyhgkD6Ji8lts/U0KrdzaTE2VhYmLY3EGgNb/0Qdr0P590L4xeHXv/21lpykmKYmi11hYQYySQQHGH3vV/M0lXlB/SekvoOoiwmcpJiAOMbf19lJlaXNTFvTArZSTE4ur04XV749D5Y+wSc/AOYfUPotR0uL5/sbOSsqaNkJbEQI5wEgiPs/9ZU8ubGmgN6T0lDB+PS7aHFXvmpcTR2uIwbfUB1axd7mruYNzaF7EQjYLSvfgbevwumXQqn/yrsnB/uqMft83PW1MxD+0BCiGOeBIIjyO/X1Ld309jhOqD3BaeOBuUFppBW9JhC+mWZMT4wb2wKWYnRnGZaT8YHP4Sxp8BFD4Ep/J/6na21pMbZmDNGCswJMdJJIDiCWjrdeHz6gAJBt8dHZUtnaHwAYEyqERR6zhxatbuZ+CgLk7MSGNO1lb9b76c1vhCueBYsUWHn9Pr8fFTUwJlTMjFLSQkhRjwJBEdQcGpnS6cHj88/qPeUN3Xi11CQsTcQ5AXXEvSYObSqtIk5Y5IxN2wj441rqdXJvDjxPoiK73XOHbXtdLi8nFCQeigfRwgxTEggOILq2vfO7e+vRMS+9s4Y2psaSoyxkhRrDfUIGjtclDQ4OSvTAU9diLLGcoftvyntiunznOsrjPITx+clH9TnEEIMLxIIjqB6x95A0NA+uPRQsOro2LS4sOP5KbFUNHfi82t+95/t5Kk6Ltlyq/Hk9a9DUn6/ZSbWVbSSHh/F6OS+A4UQYmSRQHAEBVNDwKDHCXY3OslKjCbWFr72Ly81jtIGJz94cQOr1m/gzYQ/YtVuuO41SCskKzGm7zITGD2CWblJMm1UCAFEMBAopR5XStUrpbb0OJailHpPKbUz8HNE5Sbqe6SGGjsGlxra3eQMDQ73lJ8SS1VrF+s2rmdZ4j0k0AnXvgqZUwHISoqmpq27V7nqpg4XZU2dHJ8/ov70QogBRLJH8ARw9j7HfgYs11oXAssDj0eMOoeL3BQjHXMgPYIxab0DwYRR8eSrWpYl3EOC6jLSQdkzQ89nJ8bQ6fbh6PKGvW/DnlZAxgeEEHtFLBBorVcCzfscvhB4MvD7k8BFkbr+4aS1Zl1FS7+bwQxWvaObsWl2oq0mGgcxRtDa6aa108O4PgLBuVntLE/5I3aTB254MywIAGQHViFXt4Wnh9ZVtGAxKabnJB7CJxFCDCdHeowgU2tdAxD4eUxskPtFaTNffegz1lUMsNnLINQ5XGTGR5Fmj6JpELOGdjca00N79Qiq12N+4hws+IwgMGp6r/dmJUUDULNvIChvZXJWAjE280F+CiHEcHPUDhYrpW5WSq1RSq1paGgY0raUNhozd3Y3HvxewT6/pqHDRWZCNGn2qEGlhsoC202OTYvde3D3x/DE+WCNg2+8ExoT2FewzER1695xCZ9fs7GylePzkg76cwghhp8jHQjqlFJZAIGf9f29UGv9qNZ6jtZ6Tnp6+hFrYF+qWrrCfh6MJqcLn1+TmWD0CPadPvr052Xc/ea2sGO7GzsxKcgNlJRg2+vwzCWQmAM3vgOpBf1eLz0+CotJhc0cKqptp9PtY5aMDwghejjSgeB14PrA79cDrx3h6x+UykAA6G865mDUB6aOZiREkx5v6zVr6LUN1TzzRTlu794Vx2WNTrKTYogym+DzB+HF6yDrOPj6MkgYeKN5s0mRmRAdtpZgnSwkE0L0IZLTR58DPgcmKqUqlVI3AvcAZyqldgJnBh4f9aoCAWDfgdcDURdYTBZMDTUHeghBZU1OXF4/O2odYcfGpUbDsp/AO3fC5PPh+jcgdnCF4rISo8OC1/qKVtLsttDMJSGEgAjuUKa1vrKfp86I1DUjpbLFGBuoOoQeQXAxWTA15NdGEbo0exTt3Z5QD2F9RSvHjU5Ca01jYz1/jn8EKr+AE78Hi3/dq4roQLKSYtgYmC5qnLuFmbnJspBMCBHmqB0sPlq4vX7q212YlJEaOtgppHWObpSCNLsRCGDvWoKyHoPQwTpArRVbedp/J+M71sB598GSuw8oCABkJ0ZT29aN369ZUdxAaaOTOWMkLSSECCeBYD9q2rrQGiZnJdDt8Q+6WNy+6tu7SY2Lwmo2kWq3AdDYbpxrd2h2UBzr97RC0TISnjmbROVk/aInYc7XD+qaWYnRuH1+3t1Wy7eeXsOUrASump93UOcSQgxfEgj2IzhQPG+skZfvOR3zQNQ5XGQmGD2B3j0CIxBceFwGV7Y9Bs9dgSM2lwtc/0PKlEUH3faswKKy255dT3ZiDE/dOI+EaOtBn08IMTxJINiP4JTR+YFAUNV6cGsJ6hzdZCYYi7zS+wgExyU4+cau73GL5Q0qC67kX5Meps6UfkgVQoNrCUYlRPPMN+eHApAQQvQUscHi4aKyxZjLHyzSVnUIPYJgWYeEGAs2s4mGQCDIrH6PX3sfwN7i4w7vd8jNvI7SBie5yTFYzQcfqydnxfPt0wr42pzcUMkJIYTYlwSC/ahs7WJUQjTp9ihirOaDWkvg8flpcrrICPQIlFKk2m042lrhtdv4adszVMZMJP7Gpex8tpbGilaane4+i80dCIvZxE/PnnRI5xBCDH+SGtqPypYucpJjUEqRk9x/jf+BNHa40JrQGAHA6VE7+P7OG9Drl/Kg9wLenv8UpBUyKy+JDXtaKe+n/LQQQhxuEgj2o6qli5xAWiU7KabftQT/2VTDaX/6kLe31PZ6LrSGID4auh3wxh381nEnHq0oOfdF/uS9gtx0o/7PrNxkOlxenG5fr13JhBAiEiQQDMDr81Pr6GZ0slHrJycpulePoLHDxa1L13Lbs+soa+rk3+sqe53H2KJSU9j4Pjw4H9Y9yUepV3Cl+S9stRpF44I3/Zk9CsIdampICCEGQ8YIBlDr6Mbn1+QEZu5kJ8bQ2OGm2+Mj2mqmrcvDOfd/TFunh5+cPZGyRifLNtfi8fnDBnk7a4t5ynoP+R9sNkpGf+1pVm2xU72ylNIGJ0pBXqCw3NjUOBJjrLR1eRgrqSEhxBEgPYIBBKeOBqdwBgNCsFfwUVE9De0uHr9hLreeNp7TJmbQ7vLuLevQ1QLv/ILzP7mYWaZd+M/+A9z0EYyeQ5o9Cm+gLHR2YgzRVmN/AJNJMTM3CatZkR3YU0AIISJJegQDCC4m6zlGAMaisnHpdj4qaiAlzsYJBakAnFiQiknBZ0XVzKl9EVbcA12trE3+Cv/dcTHLFlwWOndaYHXx2vIWjhsdvlvYNxeOZd7YFCyHMHVUCCEGSwLBAIIDw8EAkJO0t0fgC9TvOXVCOmaTUcQtKcrE91NXccWqO8BfD2NPhbN+y4NvdWE1h5emCC4qa+/29podtLAwnYWFQ7sHgxBi5JBAMICqli7S46NCaZvMhGiUMgLEpkpjrv9pE9PB64ZNL8An9/LdjlI2+AuIvfIh4iYtBqWoc6wMDTgHpcXvnUoqs4OEEENJcg8DqGztDPUCAGwWE5nx0VS1dvFhUQN21c2Z7a/AA7Pg9e+AzU7xaY9wkfs3fOKfDkqxpaqNnfUdFGbaw87ds9xDvgwKCyGG0IgMBMV17by6vmq/r6tq6epV6yc7KRpXYzlj1v6eVdHfJXb5LyApD67+N3xrJWNOvpxYm4WPdzbg9fn5+cubSY61ccsp4dtKJsVYQymlsD2JhRDiCBv2qSGtda+NWO5fvpNlm2s4fXJGv9U4/X5NdWs3Z00bFTjgg5IPuMt5L1Odn6O1YnfGGRRe8BPInRd6n82iWDAulU92NvLEZ2Vsrmrjb1fNIjE2/DomkyIlzkZTh2vvnsRCCDEEhnWP4BevbOY7z60PO6a1ZvXuZvwaVpc29/vehg4Xbp+fSVFN8OHv4P6ZsPRSClzbecR7Pqe47sN98WNhQSBoYWEaZU2d/PGdIk6flMG507P6vEaaPcrYk9hiPrQPKoQQh2BY9wjMJsUH2+txeX2hm21ZUycN7UbJh89Lm1g8JbP3Gzsa6P78OV60PcO8lUWAgnGnwpLf8ErrNP74RjEZ8VFMyUro87oLC9MAsJgUd180rd+tIU+dkI7X5+/zOSGEOFKGdSA4dUI6T31ezpqyFk4ab9ycv9xt9AJykmL4vKRp74sd1VD0Fmx7HV32Mfnaj5scWk/4OUkLroHE0QCM2lYHwKKJGf3e4AvS7Zw1NZMlU0aFDTbv62fnSGVQIcTQG9aBYMG4VKxmxcrihlAgWLW7mZQ4G1+bnc0HH75D17sriSl7H6rXAdAel8+z/gt5V53EjZeeS+Fx2WHnnJAZj9mk+Mpxfad7wCgz/ci1cyL3wYQQ4jAa1oEgLsrC3DEprChu4OfnTILmUjJ3Pcej0duZsWYj37O1oD9TkDMbzvgvlnlm8+13O5g3NpX//drMPjdzyUuNZf1/nSlbPgohho1hHQio28atsctpqViJ73iAJlUAAAccSURBVC/lmDtq+AngVOmYJizhhxszyZhxFj+95GS8Pj+//8sKZuQm89xNC0JTO/siQUAIMZwM70Dw7i84ueQDakwpVCXMp6XgJr6/OoEHvnkZ00YnUd+6is0V3fwUWLallormTu78yuQBg4AQQgw3wzsQLPkt2hbHhQ8WM9eeSiJW6m3VTM42irydUJDKH98uoqHdxd8/KmFcehxL+ppFJIQQw9jwDgSZU1DAKRNbeW9bHal2G7Pzk0Pf+E8YZ1QN/cu7RWyrcfDHS47DJL0BIcQIM6wXlAWdOiGdti4PpQ1O5o1NCR2fnpOIPcrC81/uITMhigtnZQ9wFiGEGJ6GJBAopc5WShUppXYppX4W6eudPD6N4Bf9+T0CgcVsYu6YZAC+efI4WeErhBiRjnggUEqZgQeBc4Ap8P/t3VuIVXUUx/HvTyctrfCSXZwxtTArumtlJWElURYZ2IMi5UMQQXeKKHqptwLpBiFE9wi7WJSIlGFRT5l2QU1Ds+uU5XS/Upmrh/0fOo1z0sbZs8f//n3gcM5/zx7OWqwzZ8357332nzmSjizzOYcPHcQxbcMY3DKAo7ssAnP+MaMZO3IIc04+uMwQzMz6rSqOEZwEfBARHwJIehKYCawr80lvOHsiH3/zy3b/9c+a1MasSW1lPrWZWb9WRSNoBT5rGLcDJ5f9pFMn7MfUdA0gMzP7RxXHCLo7LSe220m6TNIqSas6Ojr6ICwzs3qqohG0A2Maxm3AF113ioj7I2JyREweNcrr95qZlaWKRrASmCBpvKRBwGxgcQVxmJkZFRwjiIitkq4EXgIGAg9FxHt9HYeZmRUq+WZxRCwFllbx3GZm9m+1+GaxmZk150ZgZlZzbgRmZjWniO1O4e93JHUAn/Tw1/cDvu7FcHYHzrkenHP+djXfsRGxw/Pvd4tGsCskrYqIWi0g7JzrwTnnr6/y9dSQmVnNuRGYmdVcHRrB/VUHUAHnXA/OOX99km/2xwjMzOy/1eETgZmZ/YesG0FfL4nZ1ySNkfSqpPWS3pN0Tdo+QtLLkjam++FVx9rbJA2U9I6kJWk8XtKKlPNT6YKG2ZA0TNIiSe+nep+Se50lXZde12slLZS0Z251lvSQpC2S1jZs67auKtyb3s9WSzqht+LIthFUsSRmBbYC10fEEcAU4IqU403A8oiYACxP49xcA6xvGN8B3JVy/g64tJKoynMP8GJEHA4cS5F7tnWW1ApcDUyOiKMoLlA5m/zq/AhwTpdtzep6LjAh3S4DFvRWENk2AhqWxIyIP4DOJTGzERGbI+Lt9PgnijeHVoo8H027PQpcWE2E5ZDUBpwHPJDGAs4EFqVdsspZ0r7A6cCDABHxR0R8T+Z1prgo5l6SWoAhwGYyq3NEvA5822Vzs7rOBB6LwhvAMEkH9UYcOTeC7pbEbK0oltJJGgccD6wADoiIzVA0C2D/6iIrxd3AjcC2NB4JfB8RW9M4t1ofAnQAD6fpsAckDSXjOkfE58B84FOKBvAD8BZ517lTs7qW9p6WcyPYqSUxcyBpb+BZ4NqI+LHqeMok6XxgS0S81bi5m11zqnULcAKwICKOB34ho2mg7qR58ZnAeGA0MJRiaqSrnOq8I6W9znNuBDu1JObuTtIeFE3giYh4Lm3+qvMjY7rfUlV8JTgNuEDSxxTTfWdSfEIYlqYQIL9atwPtEbEijRdRNIac6zwd+CgiOiLiT+A54FTyrnOnZnUt7T0t50aQ/ZKYaW78QWB9RNzZ8KPFwLz0eB7wQl/HVpaIuDki2iJiHEVNX4mIucCrwEVpt9xy/hL4TNLEtOksYB0Z15liSmiKpCHpdd6Zc7Z1btCsrouBS9LZQ1OAHzqnkHZZRGR7A2YAG4BNwC1Vx1NCflMpPhquBt5NtxkUc+bLgY3pfkTVsZaU/zRgSXp8CPAm8AHwDDC46vh6OdfjgFWp1s8Dw3OvM3Ab8D6wFngcGJxbnYGFFMdA/qT4j//SZnWlmBq6L72fraE4o6pX4vA3i83Mai7nqSEzM9sJbgRmZjXnRmBmVnNuBGZmNedGYGZWcy073sWsPiR1nroHcCDwF8XlHQB+jYhTKwnMrEQ+fdSsCUm3Aj9HxPyqYzErk6eGzHaSpJ/T/TRJr0l6WtIGSbdLmivpTUlrJB2a9hsl6VlJK9PttGozMOueG4FZzxxLsSbC0cDFwGERcRLFpbGvSvvcQ3Ht/BOBWelnZv2OjxGY9czKSNd5kbQJWJa2rwHOSI+nA0cWl8oBYF9J+0SxdoRZv+FGYNYzvzc83tYw3sY/f1cDgFMi4re+DMzs//LUkFl5lgFXdg4kHVdhLGZNuRGYledqYHJaaHwdcHnVAZl1x6ePmpnVnD8RmJnVnBuBmVnNuRGYmdWcG4GZWc25EZiZ1ZwbgZlZzbkRmJnVnBuBmVnN/Q1Jgs3h1FmKxAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import pints\n", "import pints.toy as toy\n", "\n", "# Create a model\n", "model = toy.LogisticModel()\n", "\n", "# Set some parameters\n", "real_parameters = [0.1, 50]\n", "\n", "# Create fake data\n", "times = model.suggested_times()\n", "values = model.simulate(real_parameters, times)\n", "\n", "sigma = 3\n", "noisy_values = values + np.random.normal(0, sigma, times.shape)\n", "\n", "# Create an inference problem\n", "problem = pints.SingleOutputProblem(model, times, noisy_values)\n", "\n", "# Show the generated data\n", "plt.figure()\n", "plt.xlabel('Time')\n", "plt.ylabel('Values')\n", "plt.plot(times, noisy_values)\n", "plt.plot(times, values)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As before, we can define an error function and minimise it:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Estimated parameters:\n", "[ 0.10028851 49.80588037]\n" ] } ], "source": [ "score = pints.SumOfSquaresError(problem)\n", "\n", "boundaries = pints.RectangularBoundaries([0, 5], [1, 500])\n", "\n", "x0 = np.array([0.5, 200])\n", "opt = pints.OptimisationController(score, x0, method=pints.XNES)\n", "opt.set_log_to_screen(False)\n", "x1, f1 = opt.run()\n", "print('Estimated parameters:')\n", "print(x1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We could also however define the inference problem statistically, by specifying a likelihood." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "log_likelihood = pints.GaussianLogLikelihood(problem)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then use an optimisation routine to find the maximum likelihood estimates. (Note that the below is supposed to fail.)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "ename": "ValueError", "evalue": "Starting point must have same dimension as function to optimise.", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mValueError\u001b[0m Traceback (most recent call last)", "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mopt\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpints\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mOptimisationController\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlog_likelihood\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mx0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmethod\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mpints\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mXNES\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2\u001b[0m \u001b[0mx2\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mf2\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mopt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrun\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;32m~/dev/pints/pints/_optimisers/__init__.py\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, function, x0, sigma0, boundaries, method)\u001b[0m\n\u001b[1;32m 303\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mfunction\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mn_parameters\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m!=\u001b[0m \u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 304\u001b[0m raise ValueError(\n\u001b[0;32m--> 305\u001b[0;31m \u001b[0;34m'Starting point must have same dimension as function to'\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 306\u001b[0m ' optimise.')\n\u001b[1;32m 307\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;31mValueError\u001b[0m: Starting point must have same dimension as function to optimise." ] } ], "source": [ "opt = pints.OptimisationController(log_likelihood, x0, method=pints.XNES)\n", "x2, f2 = opt.run()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Uh oh! What's happened here?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The likelihood function we used requires a `sigma` parameter, an extra parameter it adds to the inference problem that describes the estimated noise level in the data.\n", "\n", "This means the number of parameters in our log-likelihood has gone up by one from the problem's number of parameters:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "2\n" ] } ], "source": [ "print(model.n_parameters())" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "2\n" ] } ], "source": [ "print(problem.n_parameters())" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "3\n" ] } ], "source": [ "print(log_likelihood.n_parameters())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As a result, we need to update our initial point (and boundaries) with a guess for what sigma may be.\n", "\n", "In a realistic situtation, we could try to find a flat bit of signal to obtain a first estimate. In this example, we'll just start off by guessing `sigma=1`:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Estimated parameters:\n", "[ 0.10028851 49.8058805 2.9848238 ]\n" ] } ], "source": [ "y0 = np.array([0.5, 200, 1])\n", "\n", "boundaries_3d = pints.RectangularBoundaries([0, 5, 1e-3], [1, 500, 10])\n", "\n", "opt = pints.OptimisationController(log_likelihood, y0, boundaries=boundaries_3d, method=pints.XNES)\n", "opt.set_log_to_screen(False)\n", "y1, g1 = opt.run()\n", "print('Estimated parameters:')\n", "print(y1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that the noise has introduced a slight bias into the outcome, and the estimated sigma is different to the true value of 3.\n", "\n", "As before, we can now plot a simulation with the obtained parameters, and see how it matches the data:" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzs3Xd4XMd18OHfbK9YLHoH2CSxiaREdVu9ucpFTuzEjpwocazEjp3EdhzXOHby2VbiHjkuiq1IsiU7LpJsySpUL5REUqQoNgAkelkA23u78/1xdxdYAiBBigAJYt7nwQPs3bu7A5Z77sycOSOklCiKoihLl+FkN0BRFEU5uVQgUBRFWeJUIFAURVniVCBQFEVZ4lQgUBRFWeJUIFAURVniVCBQFEVZ4lQgUBRFWeJUIFAURVniTCe7AXNRU1MjOzo6TnYzFEVRFpXt27dPSClrj3beoggEHR0dbNu27WQ3Q1EUZVERQvTN5Tw1NKQoirLEqUCgKIqyxKlAoCiKssSpQKAoirLEqUCgKIqyxKlAoCiKssSpQKAoirLEqUCgKMqSEUll+fWOQdQWveVUIFAUZcl4YNcw//CLXRwcj53sppxSVCBQFGXJ8McyAHSPxU9yS04tKhAoirJkBOJ6IFA9gnIqECiKsmSEEnogODSuegRTqUCgKMqSEUhkgZPfI0hl80zE0ie1DVOpQKAoypIRjBd7BLGTmjn0nS1dvO27z560zz+cCgSKoiwZxTmCSCrHRGHi+GTYNxJhJJwinMyetDZMpQKBoihLRjCR4Yx6F6D3Ck6WvkACgMFg4qS1YSoVCBRFWRJS2TyJTJ5z26sAODiHCeNMTuPDd27ntaHwcX3mM13jPLF/rOyYpkkGA0kABoPJ43rfE00FAkVR5l1ek/zkuR4SmdwRz3vxkB//PE2iBgsZQ2ubKrCZDWU9gvFomr/6322MR8s/eyiU5A97Rnm2e+KYPy8Qz/A3d+/gy7/bW3Z8NJIik9cAFQgURVlCtvcF+dIDe/nltsFZz4mmsvzpj1/k1ocPzEsbgnF9PL7GZWFZjassc+j+XcM8utfHjv5g2WsC8XTh+7HPJ3xnSxfRVI6+QIJMTisd7/NPDgedKkNDi2LPYkVRToxwIovTasRkXNh7wOGQfuf78J5Rbrq4Y8ZzXukPkdMkj+718W/vlBgN4qjvm9ckUko0CVohC0grPC4ep/BcX0AfChJC0OK1sXc4ii+SQpOSB3ePANA1FmVNYwUAUsL+kSgAff44PRNxpJQUc430j5Oln6ceHwwmuHNrH3VuK2PRNI/v99FW5QRg6yG9d+GwGNk3HCkbdjo8kamp0ka1y3rUP4fXSwUCRVlCrv/201y7pp4v3bDumF+b1yTZvIbNbDzm1w6H9UDwYk8AfyyNy2ZC0yAvJXlN/3qqcxwAfzzDH14bYX1zJZqU5GX5xV6/+E+/aB5NT2FOIJeX1DitjITHGAwkyeQ1dg6EABgJpQglJjN5hkMpvU2xDLHUkYe1pvrh04ewGA387eUr+eIDe+idSNBc6Si9p9EgOKuhAl80fcy/x3xQQ0OKcpoJJ7Kc++VHp01SxtM5RsIp7nl54LjG4b/x6AGu/ebTaFr5lUvTJKlsnlg6RziRZSKWZjScYjCYoM8f5+B4rHRnndckd23tp3M0RvdYjJ7xOP3+BEPBJC/3BGiutGMyCB7Z4yOczBJN5Uik8yQzGumsRjYn0bRjDwIAkUKqZoXNRLPXgSZhJJxkR1+QfOF3CifK0znDKf3xsaR57h4K82JPgBvPbWFNUwUCGJgyBDQaSVHnttLkseGLpE6JSqgqECjKaWZHfxB/PMOrg+WZLmOFidB0TuPOrX1zei8pJZmcRiyd48Hdo/QHEjzTPU7PRJxOX5Q9w2H2DEfo8hUu6oEEI6EU49E0wXiWSFK/kPsiKdqqHFQ7LWw95J/2Obm8xgFflHPaKtnYWskLh/wn/AIZKdzRu21mWrx2QJ+sfak3QIXNRKPHRuiwC364MME810CQyOT4/lMHqXFZuWFjEzazkVq3lYHAlEAQTlFfYaOuwkoikyeezp+IX+91UYFAUU4zxWGO0Uiq7Liv8LjaaeHOF/pIZScvQNm8frH3x9IMh5L0TsQ5MBplz3CEA6NRtvUE6JnQh1ae7ZoglsqRzmpoGnMyEUtT67Zy4fJqtvcHyz4b4NBEnHROY3VjBRcur2Ysmi593okSSerzI0aDoLlSDwT9gQTb+gKc11FFldNSqkVUFE7qwSM6h2GhvCb5+sMHGAom+PhVq7Ca9CG01ioHA1Oyg0bDKRo9NuorbAD4oqkZ328hqUCgKKeZYiAYOywQFHsEt1y+An88w0+e6+HgeIw9w2H2j0TpGY8zHErhj2WIpnJkclppCObVwoSmxWSY1tOYi4lYhhqXHggyOY1dg6Gy5/eORABY01jBBcuqMAhm7DkUf6+xo1w885okmSkPNpFUlgqbGQCb2UiNy8pj+3zE03nOX1aFx26educfTuqBIZnNl2X+zOT2Zw+xvS/ILZetZENrZel4q9fOUDBJXpPE0jmi6RwNFTbq3LbS73OyzWsgEEL0CiF2CyF2CiG2FY5VCSEeFUJ0Fb5757MNirKUSClLF9nRSIpkJo8/lmYgkGBP4WK+ttHD8lonP3uxn1gqN6e7+l0DIdxWE5edUcvuoXBpTH0uMjmNcDJLjcvCuqYKnBbjtIv8vpEI9RVWql1WKh0WVjdW8MIsgeBrD+/nW491HfEzf7l9gA/fvb1seCmSypUCAUCL185YNI3JINjU6sVjN08fGpry+PAg8VTnOHdt7ePel/v576cO8sCrI7xjYxPXr2soO6+1ykEmrzEe1edOABo8Nuor9GwgX/TkF59biB7BFVLKjVLKzYXHnwa2SClXAVsKjxVFeZ00TbJvJEIokcVqMjAUStI9FmO4kAkzFk1jMRpwWo28c2MzA0F9ovRo9OASZn2Lh40tlSQy+WOq3ukv5OLXuKyYjAbO66jixZ5AKZhIKdk7EmF1IW0T4MLl1fT6E6ULZ1Eik9MnmQupnLN5rnuCQDxDcMrkbySVpcI+mSjZUhgeOrulErvFiNdhIZrKlQW5cDJLhc1U+nnqn8m3t3Ry77YB7nqxn9/vHuGSlTV88OJl09rS6tWzhQaCidJwXaPHhstqwm42nv49glncANxR+PkO4B0noQ2KclpIZvKMRVMcGo+xdyTCln16ptCmtkrCiSy5/OTtfiCeocppQQjBG1bW4HWYefzA2GxvXTISTjERS7OhpZL1zR5Az4yZq2Jxt5pCPvyFy6uJpnK8MhAsvX8okS3l7xfPAXjhUPmK3v2jUTSJnqE0ywSuP5amt7BoayQ8OTYfSZb3CJoLE8bnL9NLTnjs+nPF99WkJJzM0lrlKLx+alDJkc1L/uqNy/nt31zCL//6Ij59/Vkzrn0oBYJAotSe+gobQgjqK6ylIbuTab7XEUjgESGEBH4gpfwhUC+lHAGQUo4IIermuQ2KctrQNEk0nSOa0lMrc/nyu+JOXxSrycCmVi9bDwUIJrLUuvULcDCeweu0AGAyGljTWEGX7+h39sWhpg0tlXidFlqrHLw6GOLd57TMqc3Fuvs1Lv2zz2330lBh41uPdXHrjWezb8r8QFFDhY2Oagcv9QR456bJz9k7HCn9PBhMUumwTPu8VwYm5x9GwinWNunBS+8RTAaCTa1ezmpwc/EKPehMBgI9YMbTOTQJ7V4bncMBEvEohowRoWlE/GHqCNJuCmCL57FreZB5hMyD1BBSKzzWsEuNy+yHsI74MCG52hamdhyE1LjG0ks4kMbdPwJSAyRCytLPFocJ1l0PVvec/qyP13wHgkuklMOFi/2jQoj9c32hEOJDwIcA2tra5qt9inLK0zRJJKWnYkZS2SPm0Hf6Yqysc5Uu/v54esrPGTpqnKVzV9S5eO6gn1gqh8s2+6Vg10CIGpeFpkp9cnNDs4dH9/nI5jXMh61Q7p2I8+NnD/GZN6/GYdHfcyI6OTQE+kTtv7xtLZ/8v1186YG9tFU5cFqNpTvvovM7qrh/Rw+p8DhuQxpDLk62fzfX2iNk03Ec3b1UJmwYcklELo0hn8KQS3HGwWH+ny2GyKdZtcdC67AZmUvzQzFG24CR+gcMCC3LynyWy8kgfpdFaDlWZtP8kTWF63dgIg/5HIesOQwHJdiA5wtfwBrgzTbghcLXUdwBMDLlwEP6t88WHz9yhBe3vQy1izgQSCmHC9/HhBC/Ac4HfEKIxkJvoBGYsW9a6D38EGDz5s0nf8WFoiwgKfU7/3AiSzh55It/UTavcWgixlvPbqKqcOc/tUZOMJHhHMdkNsuqOv3i0j0eY+OULJepNCl5dSjMeR1VCKEPe5zd4uF3u0fo9EVLd9tFW/aPsWswzL6RKOe2e0HLkg6PssE6QvXEyxjTIYzpMDXpEPeuHGNHZy+ueJxP2LOs/J3EkIlhzMYwZGN8NRPj69Yc/HLy/b9R/MECdBW+pv65IbhSmtGMVhKYIGrBJl1kMVEpsliEG2m0oZndaEYzGExIgwXNYCaeFTzRHeTs5hqaq92Mx/M82hng0jMbeOxAgNXNVZzbUYMURl4bjfNkV4CbLl6Ow24DYUAKo/7dYARhRAoDFL4e2O1j52AEq9lMS7WLPz6/HRA81R3g/l2jfPmG9ditJv01iML7GahzW6n0th/9L/91mrdAIIRwAgYpZbTw87XAvwL3AzcBXy18v2++2qAoi00qmyeYyBBKZEvDPuFklqc6x3nb2Y2li/FMeifiZPOSM+rdVBcCgb8wPp/M6CWYq5yTdWtW1up1+bvHZg8EPRNxoqkcG1omn1/X7MFNgtHuV7hAM2BK+DAnfJgSY9zQ1cn7LAGWP53ESxhTOsTXAATw+/L3bhBGllnd+DJWLMZKpKGaTEUbmtmFZnGTMzm599UgVV4vl6/rYCAm+NGLPt55/ir+0BnGbHPx4avXI002NKMNabKxbyzFJ3+1m09ddyaP7PURT+f4xns20j0W4+9/sZPPnLeaiwrzD4eLpXN8Yf9Wbm5cxjs2NfNc9wTf3LufZes3cnfPHi5wVtG+fhUAj8X6+D9tgD9afwmZOdREigWH2XLoEGTgj9e0kmjQL+65+AQ7d+7noPVMVtS6CCezvNIf5PIz9RHzfKUNTIu71lA98JvCP1wT8DMp5R+EEC8DvxBC3Az0A++ZxzYoyilPSkkkmcMfT8+4yvSBXcPcu22Aze1emgqZLjPp9OllHM6od1FhN2MyiFKPoFiCuco5OUbusploqLDRPRad9l4im8Aa6SO9ewd/aXyNd/seoqJvCEtsCHNsmN22KHSifxXkzG7asi4m8HDQ0MaqZcvJ2ar4+Z4EeVsV77x4PXmbl7zFQ97qQTO7QAgGggkaKmz4ZyiEtzPSyUs9AdatuoBHdg7xuNbLB9ecTyzQw96RCDdXlN8tvzLgQ6DPZ+weCvNslz7ZHElNlpeYjdNixGQQpRTS4qSxx26hwmYulZsAfULa67DMqTAeTE4Ygz7/UVRXMbmWYEWtix89c4inOsdZVuOkvdo57X3my7wFAinlIWDDDMf9wFXz9bmKcqr46kP7aa2y86cXzNy1z2sSfzyNP5aZNuk7VbE08ng0fZRAEKPSYabWZUUIgddpKaVuFgPC1B4BUrK5Oo0c3Un1a89jDR3EGj6ENdyNOaGP2K4CMEO+r4KMu5V0RTuxxot4fsLG48Mm/vJNl2D2NpN1NPBYV5hvbemiudJOPi/50cV6xvj/7NjKxS01xJtXztjuqRfJw53fUcXj+8fYPxphz3CE5ko7XoeFVq+dpzrHSWby2C2TRfC29wU5o95Nhd1Mo8dGNJ0jlsqVVgZPnSw+nBCisKisvKxEhc1UOD65utgfz5TmPOaiWNIC9NTRonp3cXVxmkPjsVLhvZ0DodMjECjKUial5O4X+6h0mPmT89vKhnQyOY2JWJpAPHPUsf9IMkv3mJ7ZM36UQnEHfFHOrHeXPqvKYSkFgEAsxQoxxLqJARqGurEH9mLz7+V76cI6gq2Qs3rIeFYQa76UtGc5fkszn34iyuZN5/Cuiw+rVjoR47f37KQ62ME72psBeGVA/32vOKuOu7b2EUvnsBgNRFK5UsbQsdrUVonJIHixJ8C+kUgprbSlEDyGQklW1rlKf1ZdY1H+eHMrAA0e/eI7Ek5OKTg3eyAAqHSYS9VHw8ksLqsJk9FAhd1M75SSF/5Yetrk9pFUOS04LEYSmXxZj8BpNeKwGPFFUtzRH8JlNeGwGNk5EOKGjc1zfv/XSwUCRZkHgbhepiGaynFwPMbKOjeZnMZ4LE1wDgGgaNdgqFTnfuIIgSCSzDIUSnLFWfrYsjHp5xrjdhpCr7Ls9yOc4dvF31oTsB00o5WU90wi7ddx0NDGN181ccM1V7Fm5UqYErDu2znEbtnDLauXT/u8ZTUuzmpw84fXRrlhQxMSeKU/yDnt3tLcQ894rFRL/1junqdyWEysa/bwyN5R4uk8a5r0FNPJonGJUiDYNRhCk3BOu16soLFwwR2NpIiksgjAZT3yJc9jt5SGhkLJbCml1GM3l60jmIhlZp1XmYkQglavg56JeCmFt3i8zm1l6yE/E7EMf35xB77C/gXZ/BwLOZ0AKhAoyjzo9U/ePT6y14fdYjpqAEjn8kSSuVK6J+jDQk6rEaMQpTTMw2lS8r+PbePtxue5OXQfjf+3HVuomzVAVhrJ5tbxkuc6fjdez83veSdp7yow6P/1c+kcz+/cysqwjTWHTUQ/0zXB8hrnrEM3b1rXwDcf62L3UBiHxUQklWNTq5fltfqQxsEpd9DHGwgAzuvwluonFdcaNFXaMQjKirm91BPAZTWVsqEaCkMww+EUkUKK7NHG9Cvt5lLJ6MjUQGAzEU3rq47TuTzJbP6Yf6dz271UOswYDvtzrq+w0dsToMZl4S1nN7KjL8iDhays9uq59zpeDxUIFGUe9EzoFxOP3cyDr47wxpW1pef+94VeQoksf3fVqrLX/PS5Xh7b7+MH799MldOClJJX+kNsbKnEF0kzHptSGVPL4hzdhnvgcbJdj/ODVJc+lj/kIlF/HqFV7+b3oQ6+vtvOnW+6jNuf7GZfPML7q1eXfabLqpdf7horX1g2Gk5xwBflg7PsJgZwycoafvxMDw++NsqKwsV/U2slXoeFKqeFg+Ox0lDM6wsEVfzomR68DnNpfN1sNNDosZe2ekxkcjx/yM+VZ9aVLvY2s5Eqh4XRcJJUVjvqsBCAx2EmnMgipSSUzJbKUBTnFiKpbGmDmirnsQ13ve/8mddDFauQvu/8NqwmI+tbKjEIfZ7gmjX1x/QZx0sFAkWZB/uGIxgEXL26nt+8Mkg0lcVtMzMWTfHrV4bQNMn7zm8r3f3n8hpPdY2Tymrcu22AWy5bQX8ggT+eYVObl219AULBAJ6D91HR+zDuoacxZiLkhYmXcmfwfNWfs+Gyd5CqWV+6248fGCO9u5NAPFMoLzHzxXhVnYv9o+WZQ0936ZOWb1xVM+vvaDUZuWp1HQ+8OsJQMMGyGmdp2GNFrZND4/FSb6L6OOcIABo9dlbUOumodpbNtbR47aUewXPdE2RyGletLi9U0OCxMRLWdwQ7UsZQUaXdTCavkczmiSSzVBaGooo9g0gyW5p3eT3BbaqLV1STzOa56iz9ol/s1ewcCB3llSeOKkOtKCdQsSjafl+U+gobFy6vQpOwo1//T/2bV4YAvfbKk52Tayl3DoSIpnK0VTl4eM8oo+EUr/SHcJHgTfkn+Fz4i9yf+DPanvgoztEXCXe8mdfe+H0uyN/O5yr/Hx3v/AKpuk2lIABMWVSWJpDIUOWY+Y54ZZ2LsWi6rHbP053jrG6sKJVKns31axvJa5Jef4JNU8bMV9S6GAwmGAolcVtNx7W95VRffdfZ/O0V5VlHLV4HIyG9vPOW/WM0V9o5s758BW4xEESS2SNmDBVVFv6MAvFM2WsqptQh8h9DIBCi/Mtg0L+MBoHRIDAZBRtaK/mHa87AZjFgNgnMJsE57ZV0+qIkMnPfHvP1UD0CRTkB8prEF0mVFnCNhJI0euysqnPjsZvZ1htgQ4uHR/b4uOLMWoZCKZ44MM6N57QghOCpznFcVhNffOsaPvKzl3ntyV9wZeQPfMr2EtaXMoQsDfw0fx2XvPWDaM3ngTDwdOc4E5kDfObylaVNUKYqLSor9Ag2zTK5OXVh2bntXvr8cfoCCT586fRJ4sM1e+1saPGwazDMprbJivLLa116AOwLUuN+/XfOxUBSuqAKwbIaBzlNsndE3yXtr964DE9hDF4/T7Cyzsnj+8fIaRrrWzw0eGwYCs+J4vshQIBBwJkNeiDJ5jUkem9pVb0LrTC547KZGC/M1Vy8ogqr2Vh6r8n2zW1twZG8fUMzP39pgP2jMVqr5j+NVAUCRXmdwsksw6FkaS2AlJLhcJK1TRUYDYJz27283BOg2mUhm9d49zkt7B4Kc9uTBzk4HqfFa2drj593L8uy/sC32Gq7F8+4H7+sYKv3LbS88c/4Q6iVWx/t5Hv2dbQLvSM/GExgEMyab17sEQyHkiQy+bJslalWFLJuXu4NkMnleeLAOAahzwHMxXvPa8NkHGRt02TRuOKcQSiZLWX1FO+ETQaB0WDAKETZ3bFBiMKx4mP9gm8Qk4+nXmRT2Tz/8Ugnd23txyDg5jcsL00QF5UKziVzNFXayybiZ1JMS+0vbC1ZX2HDZjbSWKizFElmGYulqXJacM5hzuF4ndNeid1s5Nmu8QWZJ1CBQFGOU16TDIeSpbzzomAiSyqrlRZ/FRdF/eaVIS5eUU2L10Gl3cIPnz7Ek/tHuMa8m+/zEy7vfxUGDASbLuejg5t4KH02nz3vbKrqq6jR9KqbE7FM6cI/EExSX2HDYpp5hNdh0evdFyeCq2cJBA6LiRavnd/vHuH3u/XKaJeuqp2xsufhjAbB5g4vF62oxmgQmI0GTEZBi9dWyslfVe9iXXPFCblTnqoYwPaORLj0jNppQQDKg2TVHH6f4lzGgUJV1uJjr2Oyd+Ur7Dk8n6wmI+cvq+LZ7omjn3wCqECgKMchmsoyGEzOuCJ4OKRPYDYVFjRtbK3EaBDkNcmN5+qLndymHP9c9zxXd/4f7YwwbvTi2/Qxgmf9CTlnAw27hnBtG2Rd4Y62uCBr6lqCwWCitPfubKqcFroLm8h4j3Ah/NR1ZzIY1IezGgqbpgCYTQKL0YDFZMBiNGA2GjCbDJiN+vEjXdzXNXl4tnuC5kr7CQ8CoC8Oq6+w4oukufHcmUtid0xJv5ytRzRVdWFCvbMweV6cBzAbDVTYTATiGUYKew7PtzesrOHfHtzHSFj/e5lPKhAoyjGQUjIaSTERzcx6znBh85Fij8BpNekrYqXkDE+emle+TfWen7AuFWCntpxbcx/FevY7+OC5k+mkb9/QzFvWN5VSIaudVgxicnVxXpMMhZJl4/IzqXZZGBrU2zNbuqPBAGubPZzT7sVqMmI16xd9q+nIF/qjWdtcwbPdE/N6EVtZ5yKRyXPtLMMnlQ4LFTZ9jcNcegQWk37BPzShB8+pf2bVLqveI4ik2Ng298Vkx+vK1XX445lp6w7mgwoEijJH6VyegUCCQCzLfz56gKZKO9evbZhW/2c4lMRkEGXj0Z+5ooGa1/6Hmnv/B2MmQqT1SkbX/TUfeAiiWp5vndk47fOmLn4yGgReh6W0qGwsmiKbl7R6j94jmPqz1WzAbjaWvtvMxml7Cpwoxd7MkeojvV6fvn414WT2iFlJ7dVOdg+F8TrnNqZf47JyaCKOEOW9qCqnhdFwCn88U1YmYr6sqHXx6TedNe+fAyoQKMqchJNZBoMJNA3uerGPF3sCGISeDnp2s4cPX76ilDM/HErR4LFhNAhELkn1np9Qt+s2jJkI4fbrGdv0d6Rq9No9b13fx+6hMMtrjp4ZUuOyloaGBgL6Xf6RCraZjKI0PGUxGbhgWRXGebroz+S6tQ189V3rS1tBzof1LZ6jntNe7dADwRx6BKD3og5NxKk6rLpoldPCCwf9ADPORyxmKhAoyhEcPhR0YDTKA7uGecv6Rv5ocyuP7fPx6x2D3PF8L597yxpA7xE0V1jwHriH+u3/iTnhI9J6Jb7NnyJVvabs/f9klsqkM6lxW0uFz4oralumBAKLSd+Y3mkx4bAasZqMrKrXJ1Tr3NYFDQLF9rx3ltW0C6lYpmGuK4GL8wKHn1/lsBBL63n9C9EjWEgqECjKLHJ5jf5AorRHQDav8d3Hu6h2Wfizi9pxWEz80eZWEpk8v3llkEA8Q6XDTFNkJ1/P/4yWZzqJ151L/5XfI9FwwetuT63Lwsu9AaSUDAQTVDrMNHvtuGwmnIUL/+GKF6z5znI5lb33vDZqXdY5ZUHBZKbQ4auhq6Y8XojJ4oWkAoGizCCVzdPnT5DJTVaA/PWOQfoCCT7/lsn9eAGuWV3Pr3YM8uLuvfxZ5Ef83PQAEVlH/xXfI7z8bWUVPV+PGpeVTE7DYjIwHstwVoObtqMUJStufFJ3AhZ1LVatVQ4+eMmyOZ9fzByqPqwkx9T0WzU0pCinuWgqS39Anw8omoilueflAd6wsobzl5VvddhcaeUfq57jg6/9BIfI8u3cO6m7+p9Y19F0QtpjMRmodJhZ36yPh+elpGc8xts2HP39ixespdwjOFY1s/UICoHAaTHinsfFZCeDCgSKMkUwnmEolJxWLnrLPh85TXLTRR1lxy3hQ7Q8/QnWJ7bxfH4Nv23+R37Ra+d/qmfeF3eujAZBpcNMpcNc6n10FCaUdw+GiaRypRW7R1LntlLttHBWg/uo5yq64hzB4T2CYiA43XoDoAKBopSMRVL4ItNr/mtS8ug+H2cXatUAIDWq9/yEhpe/hmaycuiSW/nLZ1pJ9GpYjIbjrrbptpnwOixU2E3TcviLZQ6eKezDO5dAYDYaeO7TV2JZ4Inixay4mU7VLD0CFQgU5TQ1FEoSiM28SGz3UBhfJM37Cxk+5tgwLU9+HNfoViKtVzL0hq+SczZwma+bh14bpdFjO6ZFQEaDoMqp1/CfrVwEQI3TitkoSmUH5hJYUxxrAAAgAElEQVQIgNdd+XOpWVbjxG0zlTbBKSoFgor5XeV7MqhAoCxpUkoGAsmyEsyHe3SvD6fVyEUrqqnoeYjmZz+F0HIMXvofBFe9pzQZfO2aBh56bXTOC6isZgM1LiuVdjOGo+ycBWAwCBo8NgYCSZwW42mXwniqqHVb2f0v1007Xu20IgQ0V55+f+4qEChLUjav0TcRx2Q0EE3NXvM9lsrx/MEJ3nxWFcu2fp7q/XeRqN3AwBXfJVPRUXbuilon162pZ+NRyj7YLUZq3dbSZifHotFjZyCQZEWda17q9yizs1uM3H7TZs5umf/yEgtNBQJlSbrj+V6++tB+7vjz84+4YclTnWPUaz6+NPFlKsN7GV//1/g2fxJpnD4HIITgI1eumuFddA6rkTq39XVlnDQVxqeLewgoC+vKsxZm68iFpgKBsuRomixsViLpDyRY1zxZpqA/kOCL97/G8hoX53VUEdj1ex60fhNn0kDvNbcTbb/mmD/PbjFQX2E7ISmHjYVhpxVznB9QlLlQgUBZUjRNcmgixp5hvb7/QLA8ELw6GGIilgEZ49yBn/JV072MO1cy8tbbpw0FHY3ZJGiosM15RetcFHsEK1SPQDmBVCBQlgxNk/T44xwci5cmhwcKO1EV9foTVFs1Hu24G+/B33Cg9joS138Lk3Xu2wUKoefv17isc5oEPhYbW71UOy1snGXbSUU5HvMeCIQQRmAbMCSlfKsQYhlwD1AF7AA+IKWcvbi7opwAUkr6AgkS6Tz7R/XegNNqZCCYLDsvND7E3aav4D14gNFzP0lm40cwHcOkrMdupsEz+65hr9f6Fg/bP3/sw1OKciQLscrkY8C+KY+/BnxTSrkKCAI3L0AblCVMSn0uIFbIDto3EsVpMXJeR1VZj8AcOsit4U+wXOul76r/ZnzTR+dcJ8hsErTXOGirdsxbEFCU+TKv/2KFEC3AW4AfFx4L4Erg/wqn3AG8Yz7boCiDwSSR5GSK6P7RCGc2VNBW5cAfz5DI5LCP7WD5/e/ESYJ71n6fyLI3z/n9q10WzqhzU3Ga1Z9Rlo75vnX5FvApoFi+qxoISSmL/ysHgeaZXiiE+JAQYpsQYtv4+Pg8N1M5XR2+uXw8naPPn2B1o5u2Kr1yZ3bfH1j++/eSMrp5V+ZLWNvPn9N7W0wGltc6aaq0n/C5AEVZSPMWCIQQbwXGpJTbpx6e4dTpu38DUsofSik3Syk319bWzksbldPbWDSF/7CyEQd8USSwuqGCVq+Dtxpe4KJtHyXlXcX3V3yfPtlAe9WRSzuDXodmVZ0Lp1XlWyiL33z+K74EeLsQ4s2ADahA7yFUCiFMhV5BCzA8j21QlqhgPIMvPL2A3P6RCAYBq+pdNB78BdeYv0ePbT2pN/+cfU8OU+OyHvHibjQImr3241oVrCinqnnrEUgp/1lK2SKl7ADeCzwupfxT4AngxsJpNwH3zVcblKUpmsoyFErO+Ny+0Sgd1U5au+6k7blP87JpE1+o+BKaxU2vP07HETZ6cViNrKxzqSCgnHZORnrDPwH/IIToRp8zuP0ktEE5TaWyefoDiWn7CQDkNcmB0Si32B6h6YUvEm6/ntsavszBoCSX1xgMJmmvnnm9QI3bwvIap8oIUk5LCzLAKaV8Eniy8PMhYG6zcYpyDLJ5jV5/vGxnsan6Awnep/2O90zcSbjjTfRf+T2at43w9MF+ev0JcposbXReZDBAS6UDj0P1ApTTl5rpUk4Lmibp88fJ5mbMPQDA8cqP+IL5TnzN1zJ25ffAYKbFa0cCzxVq/E8dGrKYDLRXO1Q9f+W0pwKBcloYDCZJZmbpCgDeA/eyvu8bbOE8aq+9DWHQ7/CLKaTPdk9gENDi1R87rUbaqhyY1M5eyhKg/pUri54vkjrixjKeg/fT/MyneNm4iduqP4uYUkK6qdKOQcBoJEVzpR2z0YDXaWZZjVMFAWXJUP/SlUUtlMgwNsM+w0Xu/sdoffLjxOvP4+b0x2mpKS/WZjYaaPTopZ3bqp3UV1hp8TrUpi/KkqICgbJoJTN5BoMzp4kCOHzbaNtyC8nq1bxwwW1EcmY6aqZnBbV49UBwdrOHOrX9o7IEqUCgLEq5vEZfID5jmiiANdhJ+yN/TtbZRO91d9AV1u/wl80QCFoL8wQb21RpZ2VpUoFAWXSK1URnyxAyx4bp+MMHkAYLPW+6k7y9hl5/AoOAVu/09NDLzqjFbBSsn7JBjaIsJSprSFl0RiMp4un8jM8ZMhE6Hr4JYybKobf+kqy7DYCeiRgt3vIS0UaDYFmNkzWNFVy8oppql3VB2q8opxrVI1BOSXdt7Zu2exjok8MT0Vn2MdKytG25BWvoIH1X/4BU9drSU73+BB1TVg2bjILltU7sFiNCCBUElCVNBQLllOOLpPjcb1/j5y/1lx1PZY8wOSwlzc9+FvfQMwy98avEm99YeiqWyjEeTZfmB0xGvSegFoopik4FAuWUs29E30qyb0qPQNPkrDWEAGpfvY2qznsY2/hRgmf8Udlzvf44AB01DhUEFGUGKhAop5wDo1EA+v2TgWAwmCSdnXnlsLvvERpe/hqh5W/Hd+4npj3fM6EHgpW1LhUEFGUGKhAop5z9hUDQV7iTn4ilZ105bPPvo/WJvyNRu4HBS/9jxj2Ge/1xKmwmNnd4VRBQlBmorCHllFMMBJFUjpFQEn985slhY9JP+6M3o1nc9F39I6Rp5sVgPRNxVjdWYLeof+6KMhPVI1BOKdm8RvdYlOW1+sTuy73BGecFRD5D+5a/xpQcp++aH5NzNsz4fpqUDAQTrG1SawQUZTYqECinlJ6JONm85No1+oV9phRSgIYXv4Jz9CUeW/V5Hg428XTnOCPhmTOKUlmNsxrd89ZmRVnsVF9ZOaUUM4YuWOblv5/SF48drrLzl9Ts/Sk/N93AP+9cDhwAwGEx8m/vWM/KOhcATZU2th4KALCmsWJhfgFFWYRUj0A5pRwYjWIyCCodFrwOM6Ph8kBgH99F83OfIdxwEZ+L3cg7NjZz25+ew3++ZwMuq4kv3v8aA8EEdRVWLCYDD+waxmgQpeCgKMp0KhAop5R9IxFavHZMBgMNHnvZcI8xFaDtsb8mZ6/l4dVfJY+Rc9u9tHodnFHv5ss3rMNgEHzxvj3ct3OIy299kj/sGeUvLulQ2UKKcgQqECinlL0jEdqq9Inixgrb5NCQlqf1ib/DlPLTd/V/szuoj2pOvdNvqrRz641nk87l+fcH97OizsX9H7mEz75lzYL/HoqymKg5AuWU0R+I44ukuW6tXiG0wWPjiQMZMjmNll3fwj30NINv+BqpmrPpfHEvzZV2XNbJf8JWs4Erz6rnV7dczEg4xRtX1agNZhRlDlQgUE4J6VyeFw76AUrF4Ro9NiSgdT5C/SvfJrDqPQTPfC9SSjp9UTa2Tu4fYDQIOqqdGA2CVfVuVtWrLCFFmSs1NKScdFJKBgJJeib0VNFiIGiosNHEBOtf+iTJqjUMX/IVEAJ/PEMwkeWMwsVeCGivLi8xrSjK3B3T/xwhhEEIofLwlBNqLJommcnTOxHHaTFS49I3l290G/mu5bsILUv/Vd9HmvQtJTt9+srjYiBorrTjtKrOraIcr6MGAiHEz4QQFUIIJ7AXOCCE+OT8N01ZCuLpXGnz+T5/nI4aZ2lc/8w93+RcQxd31/8jGc+y0ms6fTFMhU1latwWvE7LSWm7opwu5tIjWCOljADvAB4E2oAPHO1FQgibEOIlIcQuIcQeIcSXCseXCSFeFEJ0CSHuFUKo/8VLVF7Tyz+APjzU60/QXhgWcvc9Su3uH3Kf6U08qF1c9rouX5RlNU68TguNHvuCt1tRTjdzCQRmIYQZPRDcJ6XMArNUhS+TBq6UUm4ANgLXCyEuBL4GfFNKuQoIAjcfX9OVxW4omCztO9znT5DM5lle48QcG6bl6X8gWb2WX9X+Tdnq4rwm6RqLcVZjBW1VjtneWlGUYzCXQPADoBdwAk8LIdqByNFeJHWxwkNz4UsCVwL/Vzh+B3qAUZaYYDxTVlr6yc4xDAIuaK+g9YmPIrQc/VfeRk1lBb5ICq1QeW4olCSZzXPJimqMBpUaqignwlEDgZTyO1LKZinlmwsX9z7girm8uRDCKITYCYwBjwIHgZCUMlc4ZRBoPs62K4tUOpdneMqK4bwmefLAOOe0eTnjwPdx+l5m+OJ/I+NZRoPHRjYv8cf0UtTFieLNHVUnpe2Kcjqay2RxvRDidiHEQ4XHa4Cb5vLmUsq8lHIj0AKcD6ye6bRZPvdDQohtQoht4+Pjc/k4ZRGQUjIYTKJN2WzstaEw/niG9zf0U/fKdwiuupHQqncBlOYARguBYyCQwG01sbzGOe29FUU5PnMZGvop8DDQVHjcCXz8WD5EShkCngQuBCqFEMVcvxZgeJbX/FBKuVlKubm2tvZYPk45hY1H0yTS+bJjjx8Yo8kS581dXyDjWcbwxV8uPddQoW8281JvkPt2DfFs9wRnt3owqGEhRTlh5hIIaqSUvwA0gMKwTv7ILwEhRK0QorLwsx24GtgHPAHcWDjtJuC+42i3sgglMjnGoumyY6lsnhcOTvBfrp9gSgXpv+J7aObJu/1atxWTQfDbnUP8+Jke3DYTH7iwfaGbriintbmswokLIaopDOEUMn/Cc3hdI3CHEMKIHnB+IaX8nRBiL3CPEOIrwCvA7cfXdGUx0TR9SOjw3ca2HvLzLu0RNiWeZ/iCL5CqWVf2vNEg+Pxb11DpMPOGlTVUu6wL2GpFWRrmEgj+AbgfWCGEeA6oZfKOflZSyleBTTMcP4Q+X6AsISORFOmsNu14155tfN18F5Hmy/Cv+4sZX/um9Q3UuWfej1hRlNfvqIFASrlDCHEZcCYggAOFtQSKMieRVJZAbPoG9H0+P7dM/Ds5s4Ohy/4TxPSRSpfNpIKAosyzowYCIcSfHXboHCEEUsr/nac2KaeRXF5jMDB9L+FOX5T4/Z/lrYZ+dl38QwyOumnnmIyCVq9aOawo820uQ0PnTfnZBlwF7ABUIFCOqJgqmtfKJwZ2D4XZ8rufc7vxQfpXvh/DmddPe60Q0FblwGRUFUUVZb7NZWjoo1MfCyE8wJ3z1iLltPBU5zj/cO9OPnHtmaxr9pSO7+gLctuDL/Kg+fvEK1YSecMXZnx9nduqKooqygI5ntutBLDqRDdEOb08vGcUfzzDl363h70jekWSl3oCfPn3e/hP2+14DXGGr/oe0jR9/N9pNVJXoeYFFGWhzGWO4AEmV/8agDXAL+azUcriJqVke2+Q5TVO0jmNf7l/D+8+p5l7Xh7glorneUPqRUYu+Dyp6ul7CRsNglZVTE5RFtRc+t7/MeXnHNAnpRycp/Yop4GhUJKD4zHevqGJt29o4jO/2c1dL/ZzRU2UjydvJ9Z0CRPrZi4621Jlx6zmBRRlQc1ljuCphWiIcnpIZHJs6w2S0yQr61xUu6z8+zvX8+S+YT4x9HfIrIWBy74xY6pojdtChc18ElqtKEvbrIFACBFl5oJwAr3KtNqyUimjafrew8UKocUN5KtdVm4x/BrXxC76rvo+OWfjtNfaLYZSXSFFURbWrIFASuleyIYoi99QKEkmp9E1FsNtM1Hv1stBOHzbqdv5XYKrbiSy7C3TXicEtHgdpS0qFUVZWHPOzxNC1KGvIwBAStk/Ly1SFqVwIksooS847x6LsarOhRACQyZGy5MfI+tsZviiL8342qZKOzazcSGbqyjKFHPZj+DtQoguoAd4Cn23sofmuV3KIpLNawyF9NXDqWyePn+clXV6h7Jx679iiQ0ycPk30CzTO5keu5kqtfm8opxUc0nP+DL6PgKdUspl6CuLn5vXVimLytTVw70TcTQJq+pcVPQ+TFXnPYyffQuJhgumvc5sEjSrEhKKctLNJRBkpZR+wCCEMEgpn0DfjF5RGI+miaVypcedY/o21WvdSZqf/SeS1esYO+fvZ3xti9eh9h1WlFPAXOYIQkIIF/AMcLcQYgx9PYGyxCUzeXyRVNmx7rEoVXYz67d/BkM2zsDl30Yapw/91LgtuFQJCUU5JczaIxBCfE8IcQlwA3pZiY8Df0DfgP5tC9M85VSlaZKBYGLaRjNdYzH+1vUE7sEnGbngc6S906uRqFRRRTm1HOmWrAt9VXEjcC/wcynlHQvSKuWUNxxOTttoJpHJYQt182e2HxNtuZzA6sMrmKtUUUU5Fc3aI5BSfltKeRFwGRAAfiKE2CeE+LwQ4owFa6FyygknsgTj0/cm6vGF+Kb5v8ibnAxe+h/6Vf8wjR6bShVVlFPMUSeLpZR9UsqvSSk3AX8CvAt9E3plCcrkNAZDiRmfa3zlG6w39HLoon8nN8NGM26bSe05rCinoLmsIzALId4mhLgbff1AJ/DueW+ZcsqRUp8X0KZvPYxt6AXeOHY3j1ivJX/G9NXDRoNKFVWUU9WRag1dA7wPeAvwEnAP8CEpZXyB2qacYsaiaRLp/LTjhnSYxsc/Rp9WR/95X2B6JSFo9qqqoopyqjrS/8zPAC8Aq6WUb5NS3q2CwNIVS+cYi6RnfK75+c9jT4/xL+aPsfmM1mnPe51mPHZVVVRRTlVHKjp3xUI2RDl15fIaA4GZ5wU83b+l8uBv+c/sjbRuuHTaAjGLyUCTRw0JKcqpTPXVlaMaCCb56fO9fPfxrrLj5mg/zc99lkO2tfxQvoNr1zaUPS8EtFbZMajVw4pySlOBQDmisWgKfyzN/TuHeXSvj7HiSmItR+sTH0Mi+VDiw1y0sh6vo3wFca3bisOiVg8ryqlOBQJlVvHCvMC23iDJbB4JbNk/BkDdzu/iHNvOg62fpDtTzVvOLp8itluM1LlVqqiiLAbzFgiEEK1CiCcKi9D2CCE+VjheJYR4VAjRVfjuna82KMcvl9foD+glJJ7qHMfrMLOuqYIt+33YRl+m7pVv41/xLv61fy1nNbg5s36yxHRxSEitHlaUxWE+ewQ54B+llKvRy1j/rRBiDfBpYIuUchWwpfBYOcX0BxLk8pJ4Ose2vgBvWFnDdWsbSEYCND72UbKuZm6v+BsmYhnef2F72UW/qdKO1aRWDyvKYjFvA7hSyhFgpPBzVAixD2hGL2J3eeG0O4AngX+ar3Yox84XSREvrBfYeshPNi+59IxallU7WGW9HVvKx4ErfsndfwiyocXDhpbK0msr7Ca10YyiLDILMpMnhOgANgEvAvWFIIGUcqSwBaZyioiksmXrBZ7uGqe+wsqZ9W6qDtzDZrGVW/PvY7SvmnByiA9c2FE612QUNFeqVFFFWWzmfbK4sJfBr4CPSykjx/C6Dwkhtgkhto2Pj89fA5WSdC5ftl4gnMyycyDEpatqsYW6aXrhi/hqLuK27Fv41Y4hzu+o4syGybmBZq8dk1o9rCiLzrz+rxVCmNGDwN1Syl8XDvuEEI2F5xuBsZleK6X8oZRys5Ryc21t7Xw2U0HfX6DfX15H6NnuCTQJly930/r436KZHfiv+TatVS4A3n9he+ncKpeFCptaPawoi9G8DQ0JffbwdmCflPIbU566H7gJ+Grh+33z1QZl7oZCSVKH7S/wVOc4rVUOLuy8FXtwPz3X/S95ZwMfutTGUDDJshonAFazgUa10YyiLFrzOUdwCfABYLcQYmfh2GfQA8AvhBA3A/3Ae+axDcocjEfThBLl+wsMh5LsG4lw61ldVB/4GWNn30Ks9XIANrRUliaIhYBWr0OtHlaURWw+s4aeBWa7Olw1X5+rHJtoKjtt32HQF44tN4zwzsGvE687F9/mT8z4+voKG3aLShVVlMVMrf9fwtK5fGnR2FR5TfLsvgHuctwGRjMDV/4XGKaP/7tsJmrV6mFFWfRUIFii8pqkzz/zJjO7h8J8JP0jlpsO0nvlT8i6mqadYzQIWtRGM4pyWlC5fkuQlJL+QGLa5vNFyZfv4k9MTzCy/m+Its08iqc2mlGU04fqESxBw+EUsVRuxue00df4gP9bdNo3kD5v5nmBKpdFbTSjKKcRdUu3iGTzGv3+mTeImauJWJpALDPjc4ZMlNbHPkwMO11v/DYYpt8n2FSqqKKcdlQgWER+/lI/V3/zKcLJ7NFPnkEklWUkND1DCACp0fLk3+NJDfKv1k/Q2rZs2il6VVGVKqoopxsVCBaR3YNhMjmN3olj3zo6kckdsTdRu+s2PP2P8G/ZP6VxwzUzlpBuqrRjM6tUUUU53ahAsIh0j8cA6PUfWyBI5/L0TkxPEy1yDT5F/bZbedp6Gb82v41r19ZPO8djN6uqoopymlKBYJGQUtLt0wNB3zHME+TyGr0TCfLazFHAEuml9YmPEHGv4q/DN3HDppZpewlYTAaaVaqoopy2VCBYJMaiaaJpPdNnrj2CvCbp9cfJ5GZOEzVkYrQ/+pcAfMH+aQxWJ29eP30D+rYqB0Y1L6Aopy0VCBaJ7jG9N2AxGebUI9AKQSCZmTkIIDVanvo41tBBXt78Te7rt/G2s5umbTbf4FElJBTldKcCwSJRDASXrKim7wg9glg6x+d+s5vXhsIkCruMzaR++zfw9D3CyIVf4L8HWrGZDbzt7PIVxB67mRqXKiGhKKc7FQgWie6xGG6rifOWVTERyxBNzZxC+uieUe56sZ9nuydmfa/K7t9Qt/M7BM58L33L38+z3RNcu6aBiimLxNS8gKIsHSoQLBLdYzFW1LlYVq3vATDT8JCUki379X1+xqPpac8DOEZfpvnpTxJrvIjhi7/Ci70B8prksjMmN/8RAtqr1byAoiwVKhAsEt3jMVbWuWg/QiAYDCbZ3hcEYDw2PRBYIr20P/ZXZN0t9F/1A6TRwvMH/dS4rKyqc5XOU+sFFGVpUYFgEQgns4xH04VA4ADKM4eklAwEEhwcizES1lcOH94jMKZCtD/85yA1eq/9CXlbJYlMjh39QS5eUV1aQOZ1qvUCirLUqKJzi0BxonhlrQunVd8DoDhhrAeBJOFklt1DYQCqnZayHoHIpWh/9C+wRAfoedPdZDx6+YiXe4PkNMklK2sAsFuMNFeqeQFFWWpUj2AROFgMBIXhm45qB73+RKmcdLH20GtDYZxWI5vbvZM9Ai1P65Mfw+HbzuDl3yTReEHpfZ/rnqDKYeGsBjdGg6CtyjFjaQlFUU5vKhAssP9+6iC/2j54TK/pHo9hMRlordKHhdqrnfT54/T6E0SSk+Wkdw+FWdfkoa7CRjSVI5XJ0bj1S3h6H2Lkwi8QXv620rmpbJ7t/UEuWlGtB4FqBxaT+uegKEuR+p+/wO58oY9f7TjGQDAWY3mNs5TF0+a144ukmZgyD+CPpRkOp1jX5CltH+nZ9i1q9v6U8fUfwr/u5rL33N4XJJPTuGRFNfUVNlxWNUqoKEuVCgQLSNMkY9HUrKmds+kei5WGhTI5DXth9W9xYhgozQ+sa/ZQ67Jyk/Fhztj7HYKrbmT0/M9Me8/nDk7gsZu5eGWN2ndYUZY4FQgWUCCRIZuXTMyQ2jmbVDbPQDDByjoXyUyeg+Ox0mrfkXCydN5rQ2GcFiPLapyc7X+QL5nvoKvqMgbf+HUQ5X/NubzGtt4gb1hZU8pCUhRl6VKBYAH5IvodfDCRJZufpQbQYQ6Nx5FSz+0/OB4jl5c0evQdwg7vEaxt8uDte4g1L3+GZ/NruaPx8zPuMtYzESeZzXPd2no1OawoigoEC6kYCAD8s2wXebjiHgR2k7G0n4DTasJjNzMSShbeS58fuNG5k7bHP0KidiOfNn+akVlKEu33RQE4b1nVcf4miqKcTlQgWECj4ckhobnME2iaZEdfEIPQewRTNXpsjIRTpLJ5vr2li6sM27nx0OdJ1qyn9/r/xemunHF1Meirkpsr7TR61JoBRVFUIFhQU3sER5snSOf0+YBD4zHq3LZpqZ2NHhsDwQT/8sAeqoe28APrd0jVrKHnTXeiWdzUuq0zBpsGj43dgyHOafeemF9KUZRFb94CgRDif4QQY0KI16YcqxJCPCqE6Cp8X1JXI18kRbGO25F6BOFklu6xGKmsxlAoOa03ANDosRNMZOnwPcIPrN8iU7OGnuvvRLNUAFDrtjIRSyOn7E9Z5bKQzWsMh1Oc21Z5Yn85RVEWrfnsEfwUuP6wY58GtkgpVwFbCo+XjNFIiuW1ehroTMM2UkqGQ0n6/Qk0rfg4RXOlbdq5Z9a7+SPzs3zX/D1SdefQ86afoVknL+61LivZvCytOq6wm2jy2EpF6c5tV/MDiqLo5i0QSCmfBgKHHb4BuKPw8x3AO+br80+kdC7Pj585NOuWj3Pli6TpqHbgtpqm9Qj0oaB42SRyMJElmc3PWP/n6vj9fM34fRJNFxV6Au6y54trA8ajaewWI61evXzE9r4gdrORsxrd095TUZSlaaHnCOqllCMAhe91C/z5x+Xpzgm+8vt9PH9w9s1e5sIXSVFXYaPGbS3rEfhjabp8MZKZ8h3FhgpZQWVDQ1JSv+1Wmp//PNG2q+m99idI8/S1AMVAEEpm6ah2YCiMSe3oD7KxtRKzUU0PKYqiO2WvBkKIDwkhtgkhto2Pj5/UtgwE9Nr/g8HkUc6cXTqXJxDP0FBho9ZlZSKaJpPT6JmIMxxKkclppLLlgWC4EAhKPQItR/Mz/0Tdzu/qu4td/QOkafqwEehDQwA5TcNUuOgnMjn2DEc4V00UK4oyxUIHAp8QohGg8H1sthOllD+UUm6WUm6ura2d7bQFUQwAxTv04zEW0XsADRU2at1WRiMpusaixFJ60bgfPH2Iz/32tbLXDIWSmI2CGrcVQyZKxyM3U9V5D2Mb/46hN3xtxsViRZUOMzazAV9ksufx6mCYvCZVIFAUpcxCB4L7gZsKP98E3LfAn39cBoJ6j2DodfQIiqmjXqcZs1EwEUujTZly6PJFOeCLEktPVhMdDiVp8tixxodZ/sC7cdpRllwAABGcSURBVA09zdAl/45v8yf0/SRnYTDA8loXTZX2Uq8CKE0Ub1IZQ4qiTDGf6aM/B14AzhRCDAohbga+ClwjhOgCrik8PuWdiB5BsRxEOqvhspmJp/OlyWcpZem9D4xGS68ZCiW51NnHivvejiU+TO91dxBY/f4jfo4Q0FHtLG0yMzUQ7OgLsrLORaVD7UCmKMqkeas9LKV83yxPXTVfnzlfBl9njyAYz5Sqg3odFirtZgBCyQx1bhuBeIZ0ISjsH9XH8POa5NLog3wx9VM0VwMH3/xz0t4zjvg5QkBHjRNnoaR0k8fO/kJgyWuS7f1BrlvTcFy/g6Iopy9VhP4owsks0VQOt9WEL6pP6s51A5doKstoOEUqqzEeTWM2Ctw2E97CHXkokaXObSv1BgyC/9/enUfHVZ53HP8+s2o0M1osWbYsybZsC2Nj8ILBNpjEZQkmcbAhpGUpkMIpaZoACSkJTZuWcHpK0nJCaQskZCeBQFgKhFICB0ygLLaxsbHxKm+yLGuXLM1IGs1o3v5xr2RJlmTX1mjsO8/nHB3NvbrWvPe88v3Nfd/3vi/ba9uRnhgFb/09/+x5mv25i4iueIyerJHb9XtDoP+6ApPyAjS0x+iK9/Cd/9pMa0eci2edFgO1lFJj6JQdNXSq6B0xtHBqPsZAbb8ZP4cTjSXY0xBhX2MHXXHrk35ztJuCoB8RIS/bviPosJ4Z6A2CBZPzidTupvz311C292keTXyetRecWAgATLIfRPvqExt4fsNB7rrsDC4/S+8IlFIDaRAcQ2//wOJpBdZ2a8ewx0ZjCfY2RtnTECUaGzgUtCkSY1zQuhPovSNo6bCe+q1p7cTncfGlvE087/o2vtY9PDP9fn6QuI7i/NCI5XO5oHyIEIAjzx+8sb2eryybzu0XzzieU1ZKZRhtGjqG3v6B3iAYqp+grStOQ3uMjkEX//6aot19q4z13hG02HcEjc0t/GvWL1mx41U2mumsmfMAGyO5hPyN5GQNX0W9IZDtG/qY6eNDuARuvmAq37p8pq49oJQakgbBMVS3dBLye5hVnIPIkWacZNLQ0tFNc7S7r/lnOMYYu2nIuhPwul2E/B5aO+IE6jdwf8NfU2JqaTj7Nm7ZvIx5rWGaIp2U5AWGvXh73EJ5YZAsr3vY952Ym8WG716mo4SUUiPSIDiG6pYOSvMD+DwuisJ+qpo7OHS4k+Zo94DnAEYS7e4hlkj2NQ0BFAUMnzn0I6ZXPsPB5Dh+XP4QSxetoqJxKztq24klkpxdkjPk7/N7XUwtCB5Xp7WGgFLqWLSP4BiqW6xP5q0dVmdvZV2ExvajQ6Ar3sPvN9XQNMSsos1RqwmoIGhN+5Bdu4bHu7/JysjTHJyyiuWx79NVsgSAmRPDHGztpDESG3KyuYDPzbTC4wsBpZQ6HnpHMAxjDO1dcaqaO6iYEOJAcyeFIR876yJHHbuhqoWHV1dS3x6jprWTL396+oCf94ZDsa+DSf97DwXbn6TePZE7vP/A/BlXE9m+te+if+bEI3cBg9chyAl4KMs/MoGcUkqNBg2CfowxRGIJDnfGaetM0BLtpqO7h6KwNQxzfDiL93Y3kTQGlwhJY/jP1ZW8vrWOkrwA5YVBNlS1HPV7myOdXO9+gxVvP4c33k7DnL/kX2JX8/r2NooGzTBaUWR18CYNA+4ICsM+XVpSKZUSGR8EiZ4kkViCts4E7bH4gCafenvNgAn2lM5FYT+JpKEl2k1ByM/WmjZe31rHirOL+YsLy3n1k1p+8s4eatu6mJhjhUew5l3+9KN/ZJJ3J+35i9h34T8RGzeT7PUH6Iy3sKcxSsjv6RsdlOV1M60wRGVDhOLcAGKvV9y/f0EppUZTxgVBMmnoiPcQ6UoQicXp7B6+x7d3orgi+6JeZAdCfXuMgpCfdfua8biEG5dMwedxscCezO2jqhZWlbQxce395Bx4kybPBO42d3LTirv6JovLD1gX9q01bUeNDjq/3Fo9LBzwMHlcdt+UEUoplQqOv8L0JA3R7gSd3T1EYtb3fsv4jqi+3QqCCTm9TUNHgmBWMazb18yckty+cfwleQEWhhq5YONPqPjgLZK+MIfO/w53Vy3mQFuSm/pd7POC1rMEtW1dR60Wdu15Zdx6UTmTx2XrAjJKqZRzfBBsO9R23Bf+weraYgR97r6ndnv7Curbuzh0uJMDLZ0sn1MMgL9lJ+M3PcLvEi/QFfdSf86XaZr7FXqy8qnbvpGC0MCmnbzAke3Bo4MKw36Kc7P0ATCl1JhwfBCcjN6lJXsFfO6+9YbX7bM6hS8N7WXKa98lp+p1ku4sNpddzy27LuRvSpcyJyuXlmg31S2dfOqMgYvr5NtPF8ORIHC7hNJxAXKyvCil1FjJyCBo74rTGOmmvDA44nH17TGKcwcuBVmU46e5LUJF3f/wSuAFZq/eRcKfR938r9N01pfolBxaKj9gQ1ULc0py+fHbu0kkk6yaN2nA78kNHLnYT8oLEPS7Kc3P1ucDlFJjLiOD4PH39/PWznqeuHXxsBdeYwx1bV3MKzuympe3bT+3mydZWv8HCmml3ldGzbn30XLGF0l6rVAJYj0LsKHKWgTm3d1N3LRkCqX5AxeY97hd5GR5aOtKMH9yHlMKRg4lpZRKFUcHwSNvVbKnIcrNS6YO2L+pupWueJJttW3MLR162ca2rgSxRJLSYA95u54jf9ezhGre5QxcvNkzl1/3fIYrrries0qOniJ6weQ8frOmikfe2s208UGumlcy5HuMC/rI8ro1BJRSaeXoINjbEOXVLbX8+aIpuO2ncZsisb5lIz+uPjxkEEhPDNn1Gg96n2TFxvV4k110h8qoW3AXTyU+zYNrrbH/3ygeOkTmT87nN2uqaO+Kc9+VZ+EZNPJHxBqJNLc0j/jxTliklFIp4uggWDaziGfWV7Ojrp3ZxdbUDb1LRoazPGw60MqNi6cA4IpHCR18h5x9rxLe/xqeeITDrmwOTb6S+Jwv0jHhPBAXgd2NwHbOnZLfFy6DzSgKMbUgm4sqxjNt/MD1BEJZHiblZeH3uPnhn83DnOiQJqWUGiWODoKlFYW4BNbvb+kLgi0HDxP0ubl81gQ2bPqQ8MfrKaj5I8FD7+HqiRHzhPnvxEJejC9ixqLP8fkFUwf8ztJxVlv/Ent9gqG4RPiP6xYM2OfzuJiYmzWgk9gKEh0iqpRKL0cHQW7Ay6ziHD7c38yNi6fg6aij+MDLPBTcxqJ9mwn7DsJaiOVMoXnWjbzjOo+71wYoLcjhrstmDjmqqCw/mx/dcG7fMpDH4nLB+JCfwpBfJ4tTSp2SHB0ENO7ituA7xPe9z/Sn9pMd2c+9QFcyRNekxdy773JcMy7hqksuoidpeODJDZQVunjgmrkjDuMsyT/25G8ikB/0MSHsP6qPQCmlTiXODoJX7mblgdU0ucLU+BZQV76K+7eP55YvXEnFxDw2vriF1vpurgLe293IwdZOvr38zJMayy9i3YkU5fjxe4ZfPUwppU4Vzg6Cy+5jR3OC655rZFZ2LjluD3u8DUwrygXgnNJcHn9/Py0d3TyzvpqSvMCIbf8j6Q2A8WH/iMtHKqXUqcbZQVB8DvHkYc6dmuS9ykZyA15mF+f0jfaxho7u5xfv7mVvY5Q7L6kYdiTQcESs5wEKQ359KlgpdVrKiCvXuZPziXb3UHO4i7NLcvv2Tx8fIuhzs3pHA0VhP8sGzQc0Eq9HmJDr58yJYSblBTQElFKnrbRcvURkuYjsEJFKEbkn1e83ryyv75P+nH5B4HZJ3/bVC0qPq1M3lOVhckE2MyeEKQpnaUewUuq0N+ZNQyLiBh4GLgOqgXUi8pIxZmuq3jPo9zBrYpjdDVGmD3rA65JZE4jEElw6q2jYf+/zuMjL9pKf7dNP/kopx0lHH8H5QKUxZg+AiDwFrARSFgQAty6dRlM0dlQfwJJpBUN2ELtdQm62l7yAV1cIU0o5WjqucCXAgX7b1cCiwQeJyG3AbQCTJ08+6TedURRiBqERj/G4hZyAl9yAl6DPrQvDKKUyQjqCYKir61ET7hhjHgMeA1i4cGHKJuQJ+NyEszzkZHkJ+HTYp1Iq86QjCKqBsn7bpUDNWL25z+Mi6HcT9nsJ+t3a2auUynjpCIJ1QIWIlAMHgWuB61P1ZgGfmyyvm6DPTbbPo529Sik1yJgHgTEmISJfA/4AuIGfG2M+SdX7DR4lpJRSaqC0DIcxxrwCvJKO91ZKKTWQtpMopVSG0yBQSqkMp0GglFIZToNAKaUynAaBUkplOA0CpZTKcBoESimV4TQIlFIqw4kxKZvPbdSISAOw/wT/eSHQOIrFOR3oOWcGPWfnO9nznWKMOebSi6dFEJwMEfnQGLMw3eUYS3rOmUHP2fnG6ny1aUgppTKcBoFSSmW4TAiCx9JdgDTQc84Mes7ONybn6/g+AqWUUiPLhDsCpZRSI3B0EIjIchHZISKVInJPussz2kSkTERWi8g2EflERO60948TkddFZJf9PT/dZR1tIuIWkY9E5GV7u1xE1tjn/LSI+NJdxtEkInki8qyIbLfre4nT61lEvmH/XW8Rkd+KSJbT6llEfi4i9SKypd++IetVLP9uX88+FpEFo1UOxwaBiLiBh4ErgNnAdSIyO72lGnUJ4JvGmFnAYuCr9jneA7xhjKkA3rC3neZOYFu/7R8AD9rn3ALcmpZSpc5DwKvGmDOBuVjn7th6FpES4A5goTFmDtZqhtfivHr+JbB80L7h6vUKoML+ug14dLQK4dggAM4HKo0xe4wx3cBTwMo0l2lUGWMOGWM22K/bsS4OJVjn+Sv7sF8Bq9JTwtQQkVLgc8BP7W0BLgaetQ9x1DmLSA7wKeBnAMaYbmNMKw6vZ6wVFAMi4gGygUM4rJ6NMW8DzYN2D1evK4HHjeUDIE9EikejHE4OghLgQL/tanufI4nIVGA+sAaYYIw5BFZYAEXpK1lK/BvwLSBpbxcArcaYhL3ttLqeBjQAv7Cbw34qIkEcXM/GmIPAA0AVVgAcBtbj7HruNVy9puya5uQgkCH2OXKIlIiEgOeArxtj2tJdnlQSkRVAvTFmff/dQxzqpLr2AAuAR40x84EoDmoGGordLr4SKAcmAUGsppHBnFTPx5Kyv3MnB0E1UNZvuxSoSVNZUkZEvFgh8IQx5nl7d13vLaP9vT5d5UuBC4ErRWQfVnPfxVh3CHl2EwI4r66rgWpjzBp7+1msYHByPV8K7DXGNBhj4sDzwAU4u557DVevKbumOTkI1gEV9igDH1ZH00tpLtOostvGfwZsM8b8sN+PXgJutl/fDLw41mVLFWPM3xpjSo0xU7Hq9E1jzA3AauAa+zCnnXMtcEBEZtq7LgG24uB6xmoSWiwi2fbfee85O7ae+xmuXl8CbrJHDy0GDvc2IZ00Y4xjv4DPAjuB3cDfpbs8KTi/pVi3hh8DG+2vz2K1mb8B7LK/j0t3WVN0/suAl+3X04C1QCXwDOBPd/lG+VznAR/adf0CkO/0ega+B2wHtgC/BvxOq2fgt1h9IHGsT/y3DlevWE1DD9vXs81YI6pGpRz6ZLFSSmU4JzcNKaWUOg4aBEopleE0CJRSKsNpECilVIbTIFBKqQznOfYhSmUOEekdugcwEejBmt4BoMMYc0FaCqZUCunwUaWGISL3AhFjzAPpLotSqaRNQ0odJxGJ2N+XicgfReR3IrJTRL4vIjeIyFoR2Swi0+3jxovIcyKyzv66ML1noNTQNAiUOjFzsdZEOBu4ETjDGHM+1tTYt9vHPIQ1d/55wBfsnyl1ytE+AqVOzDpjz/MiIruB1+z9m4E/sV9fCsy2psoBIEdEwsZaO0KpU4YGgVInJtbvdbLfdpIj/69cwBJjTOdYFkyp/y9tGlIqdV4Dvta7ISLz0lgWpYalQaBU6twBLLQXGt8K/FW6C6TUUHT4qFJKZTi9I1BKqQynQaCUUhlOg0AppTKcBoFSSmU4DQKllMpwGgRKKZXhNAiUUirDaRAopVSG+z/NforlKgd3ZQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Show the generated data\n", "simulated_values = problem.evaluate(y1[:2])\n", "\n", "plt.figure()\n", "plt.xlabel('Time')\n", "plt.ylabel('Values')\n", "plt.plot(times, noisy_values)\n", "plt.fill_between(times, simulated_values - sigma, simulated_values + sigma, alpha=0.2)\n", "plt.plot(times, simulated_values)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can now estimate profile likelihood confidence intervals for the carrying capacity parameter by fixing the other parameters at their maximum likelihood estimates. In classical statistics, it is assumed that the log-likelihood near the maximum likelihood estimates is well approximated by a normal distribution. When this assumption breaks down, because the likelihood distribution is skewed, some prefer to use profile likelihood approaches to construct confidence intervals (others prefer Bayesian approaches which don't rely on such assumptions).\n", "\n", "To start, let's plot the profile log-likelihood." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAEKCAYAAADenhiQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl8VNX9//HXhyQkBBIghD2syiIIIiAIuAtVKm5ViwsgiCDUtbXfVqp+v11sf7a11bqgIioqqFBQFNzBXRZlibLvWyBA2MMWspzfH3OjqQ1Z4E7uZPJ+Ph7zyMyZOzPvuWg+ufece4455xAREfFTtaADiIhI9FFxERER36m4iIiI71RcRETEdyouIiLiOxUXERHxnYqLiIj4TsVFRER8p+IiIiK+iw06QFBSU1Ndy5Ytg44hIlKpLFy4cJdzrn5p21XZ4tKyZUsWLFgQdAwRkUrFzDaVZTudFhMREd+puIiIiO9UXERExHcqLiIi4jsVFxER8Z2Ki4iI+E7FRUREfFdlr3MRqQh5+QXsOXyMPYeOceBIHtlHc8k+mseR3HyOHMsnJ6+AvPwC8p2joMCBGTFmxMYY1WOqkRBXjYS4GGrGx5KUEEtSQhx1E+NIqVmdWvGxmFnQX1GkWCouIidh/5FcNu46xOY9h9m85zAZew+Tuf8o2/cfZceBo+w9nBu2z64eU436SfE0qp1Ao+QEmtRJoFlKIs1SEmnh/YyL0ckJCYaKi0gZHM3NZ9X2bFZkHmDl9mxWbc9mbdZBsrJz/mO7ejWr07hOAml1E+nesi6pteKpV7M6KTXjSa4RS3JCHLUSYkmsHkNCbAwJcTHExoSOVj7+eDbOOS66uC95BQUcyyvgaG4BR3PzOZiTR/bRPA4cyWXfkVx2H8xhz6Fj7MzOYfv+o6zIPMCsFTvIySv4PktsNaN5vUROrV+L9o2Tad8oidMaJ9MiJZFq1XTEI+Gl4iLyI8451mUdZMHGvaRv2ce3GftZvSOb/AIHQM3qMbRtlMQFbetzSoNatE6tSYt6NUmrW4Oa8Sf+v9RDDz0EQL9+/YipFkN8bAxJCeXLnZWdw5a9h9mw6zDrsw6yLusga3YeZNaKHXjxSYqPpWPTZDqn1aFr8zp0bV6XBsnl+CCRMlBxkSqvoMCxakc2c9ftZu763SzYuOf701m1a8TROa02F7VvTaemtTmtcTLN6kbmX/5mRoPkBBokJ9CtRcp/PHc0N581Ow6yPHM/S7buZ8nWA0z4aiPjPg8d6TRLqUGPlvXo2TqFXq3r0SwlMYivIFFExUWqpN0Hc/hizS4+W53F56uz2H3oGAAt6iXS97SGnNUyhW4t69I6tWZUdJonxMXQKa02ndJqM/CsUFtOXj5Ltx5g8ea9fLNxDx+v3MG0RRkANE9J5Nw2qZzbpj59Tq1HUkJcgOmlMlJxkSpj7c6DfLR8B7NW7GDR5r04Byk1q3Nem1T6nJpK71NTaVqnRtAxK0x8bAzdWtSlW4u63HpuawoKHGt2HmTuul18uXYX0xdvZdL8zcTFGD1apXBhuwb8pEMjmtfTUY2UzpxzQWcIRPfu3Z2m3I9+q3dkM/O7TN5bksmanQcB6NS0Nn1Pa8iF7etzepPaEXOK64ILLgDg008/DTRHoWN5BSzavJdPVu7k45U7v99/pzVO5pKODbmsU2PaNEwKOKVUNDNb6JzrXup2Ki4SbbbuO8Lb6dt4K30rK7dnYwY9Wqbw006N+UnHhjSuHZlHJ6tWrQKgXbt2AScp3ubdh/lw+XY+WLadBZtCR37tGiZx+RmNueKMpjqiqSJUXEqh4hJdjhzL54Nl25myYAtz1u0GoGvzOlzZpSn9OzWiQXmGXUmpdh44yntLtzPj220s2LQXgO4t6nJ116YM6NyE2jXURxOtVFxKoeISHVZuP8Cr8zfz5uKtZB/No1lKDa7t2oyrz6x8f0nPmDEDgMsvvzzgJOWzdd8R3krfyhuLtrJ250HiY6vR//RGDDyrOWe3TomKARHyAxWXUqi4VF7H8gp4b2kmL8/dxMJNe6keW43LOjXm592b0bNVSsT0oZRXpPW5lJdzjiVb9zNlwRbeSt9G9tE8WqXW5MYezbm2Wxp1a1YPOqL4QMWlFCoulc/ugzlMmr+ZifM2sTM7h1apNbmpZ3Ou6Rodv7gqe3Ep6sixfN5flsmkeZtZ4P0BcMUZTRjWpyUdm9QOOp6chLIWFw1Floi3Pusg47/cwLSFGeTkFXBBu/r8tXdLzm9Tv9IepUS7GtVjuPrMNK4+M40VmQeYOG8TbyzaytSFGfRolcItfVrRr0NDYvTvF7VUXCRiLd26n6c+Wcv7y7YTF1ONa7o2Zfg5rTm1Qa2go0k5nNY4mT9f3YnfXNqeKd9s4aW5Gxk1cSGtUmsy/JxWXNstjYS4mKBjis9UXCTiLNy0l8dnr+Gz1VkkxcfyiwtOYWjvVtRPig86mpyE2jXiGHFea245pxXvL93OuM/X8cD0pTw2azW3nNOKwWe30EwAUUR9LhIxFm/ey6Oz1vD56ixSalZn+DmtGNyrBclV5BfOli1bAGjWrFnASSqGc475G/Yw9tN1fL46i6SEWIb1bskt57SiTmLl70OLVurQL4WKS+RYkXmARz5YxeyVO6mbGMdt55/CkF4tSKyuA+uq4ruMfTz1yVo+WLaDpPhYhvVpyfBzWlM7sWr8YVGZVNriYma/B0YAWV7T75xz73rPjQGGA/nAXc65D7z2S4F/ATHAeOfcw6V9jopL8LbsOcw/P1rN9PStJMXHctv5p3Bz75bUOolp6yuzyZMnAzBw4MCAkwRnReYBHp+9hveWbicpIZaR54ZOo53MUgbir8peXA465x75UXsH4DWgB9AEmAW09Z5eDfQDMoBvgBucc8tL+hwVl+DsP5zLk5+s4aU5m6hWDYb1acWo806p8n+lRtNQ5JO1fNsBHp21mo+W76BezercfuGp3HR2c+Jj1fEftGgcinwl8LpzLgfYYGZrCRUagLXOufUAZva6t22JxUUqXm5+ARPnbeKxWWs4cDSX67ql8at+7WhUW1OzyH/q0CSZ54Z0Z9Hmvfz9/VX8ceZyXpyzgd9c0p4BnRvrqv9KIFIX2L7DzL4zsxfMrK7X1hTYUmSbDK/teO3/xcxGmtkCM1uQlZVV3CYSJp+vzqL/v77gDzOW06lpbd6581z+du0ZKixSoq7N6/LayLN5+ZYe1Kwey52vLeaqp77im417go4mpQjkyMXMZgGNinnqfuBp4E+A837+A7gFKO5PFUfxBbLYc33OuXHAOAidFit3cCm3jL2H+eOM5Xy4fAct6iXy3JDu9D2tgf7ylHI5r219+pyaypuLt/LIB6u47pm5XNa5MWP6tyetbuWaQ66qCKS4OOf6lmU7M3sOmOk9zACKjtFMA7Z594/XLgHJyctn/BcbeOLjNRjG/1zSjlvPbaVz5nLCYqoZ13ZL46edGvHsZ+t59vN1fLR8B6POP4XR559Cjer6byuSRGKHfmPnXKZ3/5dAT+fc9WbWEXiVHzr0ZwNtCB3RrAYuBrYS6tC/0Tm3rKTPUYd++Mxfv5vfvbmEdVmHuLRjIx68vEOVWuHxRO3atQuA1NTUgJNUDtv2HeHh91by9rfbaFqnBg8OOI1LOjbSUXGYVeYO/b+ZWRdCp7Y2ArcBOOeWmdkUQh31ecDtzrl8ADO7A/iA0FDkF0orLBIe+w/n8v/eW8Hr32whrW4NXhx2Fhe2axB0rEpDRaV8mtSpweM3nMmNPZvz+7eXMWriIs5vW58/XtmRFvVqBh2vyou4I5eKoiMXf72/dDsPTF/K3sPHuPWcVtzdt40ugiynCRMmADB06NBAc1RGefkFvDx3E//8aDW5+QXcedGpjDzvFKrHRuqYpcqr0l7nUlFUXPyx+2AO//v2Mt75LpOOTZL56zWdOb2pplQ/EbrO5eRt33+UP81czjtLMmnToBYPX9OJbi1Sgo4VVcpaXFTW5YS9tySTfo9+zofLtvPrn7Rl+u19VFgkUI1qJ/DUTV15cehZHD6Wz7XPzOXB6UvJPpobdLQqR+ctpNz2H8nl928v483FWzm9aTL/uO5s2jVKCjqWyPcubN+AD395Ho98uIoJczYye8UOHr6mM+e1rR90tCpDRy5SLnPW7uKSRz/n7W+3cffFbXjzF31UWCQi1YyP5f8u78i00b2pUT2GIS98zW+mfsv+IzqKqQgqLlImOXn5/OXdFdw4fj6J8TG8Mbo3v+zXlrgY/Sckka1r87q8c9e5jL7gFKYuzODSxz7nyzW7go4V9dShL6Vau/Mgd722mOWZB7ipZ3Puv+w0jQQLg8OHDwOQmKgrzsMlfcs+7p2SzrqsQwzp1YL7+rfXf8vlVJmvc5EIMnVhBg9OX0pCXDWeG9Kdfh0aBh0paqmohF+XZnV4565z+dv7q3jhqw18uWYXj13fhc5pdYKOFnV0TkOKdSgnj19NTufX//6Wzmm1ee/u81RYwmzs2LGMHTs26BhRLyEuhv+9vAOvjujJkdx8fjZ2Dk99spb8gqp5FidcdFpM/suaHdmMnrSI9VkHueviNtx5URtiqmlKjXDTdS4Vb//hXO6fvoSZ32XSo1UKjw3sQhNNVVQiXeciJ+SNRRlc8eRX7Ducy8Rbe3JP37YqLBK1aifG8cQNZ/KP685g2db99P/XF3ywbHvQsaKCiosAcCyvgAemL+FXU76lU1pt3r3rHHqformuJPqZGdd0S2PmXefSPCWR215ZyIPTl3I0Nz/oaJWaiouwff9RBo6by8R5m7nt/Na8emtPGiRrES+pWlql1mTa6N6MOLcVr8zbxDVPz2HjrkNBx6q0VFyquG827mHAE1+yans2Y2/qypj+pxGra1ekiqoeW437L+vA+CHdydh7hAFPfMk732UGHatSUod+Ffbq/M3839tLSaubyLjB3WjTUFfaixTK2HuYO19bzOLN+xjWpyVj+p+mWZZRh76UoLB/5XdvLqH3KalMv72PCovIj6TVTWTyyF4M7d2SF7/ayA3PzWP7/qNBx6o0VFyqmL2HjjHkhfnf96+8MPQsateICzqWAI888giPPPJI0DGkiOqx1fj9FR15/IYzWZF5gAFPfMH89buDjlUpqLhUIeuyDnL12K9YtGkfjw48gzH9T9Mw4wgyc+ZMZs6cGXQMKcYVZzThrdv7kJwQx03j5/PiVxuoql0KZaXiUkXMWbuLq5/6iuyjebw2sidXn5kWdCSRSqVNwySm39GHC9rV5w8zlnPvlG81XLkEKi5VwL8XbGHIC1/TqHYC02/vo5X5RE5QckIc4wZ355d92/LG4q0MfHau+mGOQ8Ulijnn+OdHq/mfqd9xdut6TB3dm2YpmhxR5GRUq2bc3bcNzw7uxtqdB7niyS9ZvHlv0LEijopLlMrNL+Def3/L47PXcF23NF4cdhbJCeq4j2Q1atSgRg3Na1VZXNKxEW/8og/xcdUY+Ow8pi/eGnSkiKLrXKLQwZw8Rk9cyBdrdvHLvm256+JTMVPHvUg47D10jFETFzJ/wx7uvOhUftm3LdWieKCMrnOporKyc7h+3FzmrNvNX6/pxN1926iwiIRR3ZrVeWV4TwZ2b8YTH6/ljtcWceSYOvq1WFgU2bz7MINfmM/OAzk8N6QbF7XX+iuVyZ/+9CcAHnzwwYCTSHlVj63Gw9d04tQGtfjLeyvYtm8ezw3pTv2k+KCjBUZHLlFiReYBrnlmDvuP5DJpRE8Vlkpo9uzZzJ49O+gYcoLMjBHnteaZQd1Yuf0AV4/9irU7s4OOFRgVlyiwYOMeBj47lxgz/n1bL7o2rxt0JJEq65KOjZg8shdHcwu4euwc5q6rmlf0q7hUcp+vzmLQ8/NJrRXP1NG9NEeYSAQ4o1kdpt/em4bJCdz8wte8/e22oCNVOBWXSuz9pdu59aUFtEqtxZRRvUirq2tYRCJFWt1Epo3qTZfmdbjrtcWM+3xdlZoyRsWlknpzcQa3v7qI05sm8/qIs0mtVXU7DqNFvXr1qFevXtAxxEe1E+N4+ZYeXNa5MX95dyV/nLmcgoKqUWA0WqwSmvzNZu57Ywm9WtfjuSHdqRmvf8ZoMG3atKAjSBgkxMXwxPVn0iApnhe/2sjug8d45Lozon5tGP1WqmRembuRB99axvlt6/Ps4G4kxMUEHUlESlGtmvG/AzrQICmBv76/kj2HjvHM4G7UiuI/DAMpnWZ2nZktM7MCM+v+o+fGmNlaM1tlZpcUab/Ua1trZvcVaW9lZvPNbI2ZTTaz6hX5XSrS819u4MG3ltH3tIaMG6LCEm3GjBnDmDFjgo4hYWJmjL7gFP5+bWfmrt/NTc/NY++hY0HHCpugjsuWAj8DPi/aaGYdgOuBjsClwFgzizGzGOApoD/QAbjB2xbgr8Cjzrk2wF5geMV8hYo1/ov1/Gnmcvqf3oixN3UlPlaFJdrMnTuXuXPnBh1Dwuy67s14ZlA3VmzP5udRPKtyIMXFObfCObeqmKeuBF53zuU45zYAa4Ee3m2tc269c+4Y8DpwpYXmNbkImOq9/iXgqvB/g4o1/ov1PPTOCi7r1JjHbzgz6s/VikS7fh0a8tKwHmTuP8o1T89hw65DQUfyXaT9lmoKbCnyOMNrO157PWCfcy7vR+1Ro2hheez6LsTFRNo/mYiciF6n1OO1EWdzJDefnz87l1Xbo+tq/rD9pjKzWWa2tJjblSW9rJg2dwLtx8s00swWmNmCrKyskr9ABHjxqw0qLCJRrFNabSaPPBsDBo6by3cZ+4KO5Juw/bZyzvV1zp1ezO2tEl6WATQr8jgN2FZC+y6gjpnF/qj9eJnGOee6O+e6169f/0S+VoWZNH8Tf5ixnEs6NlRhqSLS0tJIS9Py01VNm4ZJTB3Vm1rxsdz43Hy+2bgn6Ei+iLTfWG8D15tZvJm1AtoAXwPfAG28kWHVCXX6v+1Cl7t+Alzrvf5moKTiVSlM+WYL97+5lIvbN+CJG7qqsFQREydOZOLEiUHHkAA0r5fIv0f1okFyPEOe/5o563YFHemkBTUU+WozywB6Ae+Y2QcAzrllwBRgOfA+cLtzLt/rU7kD+ABYAUzxtgX4LfArM1tLqA/m+Yr9Nv56K30rv33jO85tk8pTN3VV571IFdG4dg1eH3k2zVJqMOzFb/hsdeSfui+JVqKMIB8t38GoiQvp3qIuE4b1oEZ1DTeuSu655x4AHnvssYCTSJB2H8xh0PNfs27nQZ4e1JWLT4us5TO0EmUl8+WaXdw+aRGnN63N80PPUmGpgtLT00lPTw86hgSsXq14XhvRk/aNkxg1cSGzlu8IOtIJKXHuATObQQmjr5xzV/ieqApauGkvI15eQOv6NXlp2FlRPSWEiJSuTmJo6eQhz89n9KSFPHVjV37SsVHQscqltCOXR4B/ABuAI8Bz3u0goavs5SSt2p7NLRO+oWFyPK8M70mdxKidvUZEyqF2jTheubUnHZrU5vZXF/HBsu1BRyqXEouLc+4z59xnwJnOuYHOuRne7UbgnIqJGL027z7M4OfnkxBXjVeG96zS622LyH9LTojjleE96NikNne8uqhSnSIra59LfTNrXfjAGyYc2ReKRLis7BwGvzCfY/kFvDK8J81StNBXVde2bVvatm0bdAyJMMkJcbw8vAcdGifzi0mL+GTlzqAjlUmZRouZ2aXAOGC919QSGOmc+zB80cIryNFi2UdzuX7cPNZnHWLSiJ5a815ESrX/cC43PT+P1TsO8tyQ7pzfNpi/730dLeace5/QBY13e7d2lbmwBCknL59RExeyans2Ywd1VWERkTKpnRjHxOE9ObV+LUa+vCDiL7QsU3ExszjgNuBB7zbCa5NyKChw3DvlW75au5u/XtOZC9s1CDqSRJCRI0cycuTIoGNIBAuNIutB85REbn1pAQsieKqYsva5PA10A8Z6t25em5TDX95dwczvMrmvf3uu6aY5pOQ/rV69mtWrVwcdQyJcvVrxTBrRk4bJCQx78ZuIneyyrMXlLOfczc65j73bMOCscAaLNi98uYHxX25gaO+W3HZe69JfICJyHA2SEph0a09qJ8Yx+PmvI3K6/rIWl3wzO6XwgTdyLD88kaLPe0sy+dM7oRmOHxzQgdAaZyIiJ65JnRq8NuJsEuKqMej5+WyMsAXHylpc/gf4xMw+NbPPgI+Be8MXK3os3LSHuyenc2azOvzr+jOJqabCIiL+aJaSyMThPckvcNw0fj7b9h0JOtL3yjpabDah0WJ3ebd2zrlPwhksGmzcdYhbX1pAk9oJjL/5LBLiNF+YHF+XLl3o0qVL0DGkkmnTMImXb+nBgSO5DHp+PrsP5gQdCSj7dS5xwGjgPK/pU+BZ51xu+KKFV7ivc9l3+Bg/GzuHPYeP8eYv+tAqtWbYPktE5JuNexg0fj5tGybx6oieJCWEZ0Cv37Mia7RYORzLK+C2VxaSsfcI4wZ3V2ERkbA7q2UKTw/qyorMA4x8eSFHc4PtFtdoMZ855xjzxhLmb9jD36/rTI9WKUFHkkpi0KBBDBo0KOgYUold1L4hj1x3BnPX7+bO1xaTl18QWBaNFvPZM5+tZ9qiDO6+uA1XdmkadBypRDIyMsjIyAg6hlRyV53ZlN9f3oGPlu/g/jeXEtSCkGVdOKRwtNh6wIAWwLCwpaqkPli2nb99sJIBnRtzT982QccRkSpqaJ9W7D50jCc+Xkv9pHh+fUm7Cs9QpuLinJttZm2AdoSKy0rnXGQMSYgQS7fu557X0+mcVodHrjtD17KISKB+1a8tWdk5PPnJWlJrVWdon1YV+vnlWfKwG6HZkGOBM8wM59zLYUlVyWRl5zDy5QXUSYzjucHdNORYRAJnZjx01ensPnSMP8xcTv2kBC7r3LjCPr9MxcXMXgFOAdL5oa/FAVW+uBTOcrzn8DGmjupNg+SEoCNJJdWrV6+gI0iUiY2pxhM3nMlN4+fzyynppNaqTs/W9Srks8t6ncsKoIMLqmcoDPy4zsU5x2+nfceUBRk8eeOZDOjcxKd0IiL+2XvoGNc8M4dd2TlMHd2btg2TTvi9/L7OZSnQ6ITTRKkJczYyZUEGd110qgqLiESsujWr89KwHsTHxTD0ha/Zvv9o2D+zxNNiZjaD0OmvJGC5mX0NfN+R75y7IrzxItecdbt46J0V9OvQkHv6amlaOXnXXHMNANOmTQs4iUSjZimJTBh2Fg9MX0pBBZyEKq3P5ZGwJ6iEMvYe5o5XF9MqtSaPDuxCNU1GKT7YvXt30BEkynVsUps3RveukNGsJRYX59xnYU9QyRw5ls9trywkN7+AcYO7USu+PAPuRESCVVGXSZR2WuxL59w5ZpZN6PTY908BzjmXHNZ0ESY0tct3LM88wPM3d6d1/VpBRxIRiUilHbmc4/088aEFUcTM6HNqKu0bJ3NR+4ZBxxERiVilHbmUOOuic26Pv3Ei33XdmwUdQaLUxRdfHHQEEd+U1mGwkNDpsOJO0jlAi8GL+OTBBx8MOoKIb0o7LVaxk9GIiEhUKNNFlBYyyMwe9B43N7MeJ/qhZnadmS0zswIz616kvaWZHTGzdO/2TJHnupnZEjNba2aPmzfkwcxSzOwjM1vj/ax7orlEgtS/f3/69+8fdAwRX5T1Cv2xQC/gRu9xNvDUSXzuUuBnwOfFPLfOOdfFu40q0v40MBJo490u9drvA2Y759oAs73HIpXOkSNHOHLkSNAxRHxR1uLS0zl3O3AUwDm3F6h+oh/qnFvhnFtV1u3NrDGQ7Jyb681v9jJwlff0lcBL3v2XirSLiEhAylpccs0sBu9aFzOrD4Rr/cxWZrbYzD4zs3O9tqZA0SX6Mrw2gIbOuUwA72eDMOUSEZEyKuvl5Y8DbwINzOzPwLVAiUNbzGwWxU92eb9z7q3jvCwTaO6c221m3YDpZtaR449WKxczG0no1BrNmzcv78tFRKSMyroS5SQzWwhcTOgX/VXOuRWlvKZvecN4q1vmePcXmtk6oC2hI5W0IpumAdu8+zvMrLFzLtM7fbazhPcfB4yD0JT75c0nEk4DBgwIOoKIb8q6WNhw59zzwMoibQ8753ztPPdOt+1xzuWbWWtCHffrnXN7zCzbzM4G5gNDgCe8l70N3Aw87P083lGRSET79a9/HXQEEd+Utc/lWjO7qfCBmY0F6p/oh5rZ1WaWQWgE2jtm9oH31HnAd2b2LTAVGFVkFoDRwHhgLbAOeM9rfxjoZ2ZrgH7eYxERCVBZV6KsQegI4QWgP6Gji3vCnC2s/FiJUsRPF1xwAQCffvppoDlESlLWlSjLM7fYrcB04Cvgj2aWUhXnFhMRkdKVZ26xwp+XeTfNLSYiIsXS3GIiIuK70k6LXeSc+9jMflbc8865N8ITS0REKrPSToudD3wMXF7Mcw5QcRHxyc9//vOgI4j4pkyjxaKRRouJiJSfX6PFflXS8865f5Y3mIgU7/DhwwAkJiYGnETk5JV2WiypQlKICD/96U8BXeci0aG00WJ/qKggIiISPco6/cv3zGxROIKIiEj0KHdxofjp70VERL53IsXlHd9TiIhIVCnrYmHfc849EI4gIlXd0KFDg44g4puyrueSzX+v/LgfWADc65xb73cwkapGxUWiSVmPXP5JaOXHVwn1uVxPaAnjVYSm4b8gHOFEqpJdu3YBkJqaGnASkZNX1uJyqXOuZ5HH48xsnnPuj2b2u3AEE6lqrr32WkDXuUh0KGuHfoGZ/dzMqnm3opMgVc35Y0RE5LjKWlxuAgYDO73bYGCQt0LlHWHKJiIilVSZTot5HfbFzYwM8KV/cUREJBqU6cjFzNLM7E0z22lmO8xsmpmlhTuciIhUTmXt0H+R0Eix67zHg7y2fuEIJVIVjR49OugIIr4pa3Gp75x7scjjCWZ2TzgCiVRVAwcODDqCiG/K2qG/y8wGmVmMdxsE7A5nMJGqZsuWLWzZsiXoGCK+KOuRyy3Ak8CjhIYezwGGhSuUSFU0ePBgQNe5SHQo05GLc26zc+4K51x951wD59xVwM/CnE1ERCqpE5kVuVCJSyCLiEjVdTLFReu6iIhIsU6muGjaFxERKVaJHfrHmWofQkctNcL4MjfNAAAOUUlEQVSSSKSKuvfee4OOIOKbEouLcy6pooKIVHWXX368GZZEKp+TOS0mIj5atWoVq1atCjqGiC/KvcyxiITHbbfdBug6F4kOgRy5mNnfzWylmX3nTYhZp8hzY8xsrZmtMrNLirRf6rWtNbP7irS3MrP5ZrbGzCabWfWK/j4iIvKfgjot9hFwunOuM7AaGANgZh0ILaHcEbgUGFs45QzwFNAf6ADc4G0L8FfgUedcG2AvMLxCv4mIiPyXQIqLc+5D51ye93AeUDh9/5XA6865HOfcBmAt0MO7rXXOrXfOHQNeB640MwMuAqZ6r38JuKqivoeIiBQvEjr0bwHe8+43BYrO3JfhtR2vvR6wr0ihKmwvlpmNNLMFZrYgKyvLp/giIvJjYevQN7NZQKNinrrfOfeWt839QB4wqfBlxWzvKL4IuhK2L5ZzbhwwDqB79+66CFQiygMPPBB0BBHfhK24OOf6lvS8md0MDAAuds4V/qLPAJoV2SwN2ObdL659F1DHzGK9o5ei24tUKn37lvi/jEilEtRosUuB3wJXOOcOF3nqbeB6M4s3s1ZAG+Br4BugjTcyrDqhTv+3vaL0CXCt9/qbgbcq6nuI+Ck9PZ309PSgY4j4IqjrXJ4E4oGPQn3yzHPOjXLOLTOzKcByQqfLbnfO5QOY2R3AB0AM8IJzbpn3Xr8FXjezh4DFwPMV+1VE/HHPPaHFXXWdi0SDQIqLc+7UEp77M/DnYtrfBd4tpn09odFkIiISISJhtJiIiEQZFRcREfGdiouIiPhOE1eKRIi//OUvQUcQ8Y2Ki0iE6N27d9ARRHyj02IiEWLOnDnMmTMn6BgivtCRi0iE+N3vfgfoOheJDjpyERER36m4iIiI71RcRETEdyouIiLiO3Xoi0SIxx57LOgIIr5RcRGJEF26dAk6gohvdFpMJELMmjWLWbNmBR1DxBc6chGJEA899BCgFSklOujIRUREfKfiIiIivlNxERER36m4iIiI79ShLxIhnn322aAjiPhGxUUkQrRr1y7oCCK+0WkxkQgxY8YMZsyYEXQMEV/oyEUkQvzjH/8A4PLLLw84icjJ05GLiIj4TsVFRER8p+IiIiK+U3ERERHfqUNfJEK88sorQUcQ8Y2Ki0iEaNasWdARRHyj02IiEWLy5MlMnjw56BgivtCRi0iEePrppwEYOHBgwElETl4gRy5m9nczW2lm35nZm2ZWx2tvaWZHzCzduz1T5DXdzGyJma01s8fNzLz2FDP7yMzWeD/rBvGdRETkB0GdFvsION051xlYDYwp8tw651wX7zaqSPvTwEigjXe71Gu/D5jtnGsDzPYei4hIgAIpLs65D51zed7DeUBaSdubWWMg2Tk31znngJeBq7ynrwRe8u6/VKRdREQCEgkd+rcA7xV53MrMFpvZZ2Z2rtfWFMgosk2G1wbQ0DmXCeD9bBDuwCIiUrKwdeib2SygUTFP3e+ce8vb5n4gD5jkPZcJNHfO7TazbsB0M+sIWDHv404g00hCp9Zo3rx5eV8uElZTp04NOoKIb8JWXJxzfUt63sxuBgYAF3ununDO5QA53v2FZrYOaEvoSKXoqbM0YJt3f4eZNXbOZXqnz3aWkGkcMA6ge/fu5S5OIuGUmpoadAQR3wQ1WuxS4LfAFc65w0Xa65tZjHe/NaGO+/Xe6a5sMzvbGyU2BHjLe9nbwM3e/ZuLtItUKhMmTGDChAlBxxDxRVDXuTwJxAMfeSOK53kjw84D/mhmeUA+MMo5t8d7zWhgAlCDUB9NYT/Nw8AUMxsObAauq6gvIeKnwsIydOjQQHOI+CGQ4uKcO/U47dOAacd5bgFwejHtu4GLfQ0oIiInJRJGi4mISJRRcREREd+puIiIiO80caVIhHj33XeDjiDiGxUXkQiRmJgYdAQR3+i0mEiEGDt2LGPHjg06hogvVFxEIsSUKVOYMmVK0DFEfKHiIiIivlNxERER36m4iIiI71RcRETEdxqKLBIhPv3006AjiPhGRy4iIuI7FRcREfGdiouIiPhOxUVERHyn4iIiIr5TcREREd+puIiIiO9UXERExHcqLiIi4jtzzgWdIRBmlgVsOsGXpwK7fIzjF+UqH+UqH+Uqn2jN1cI5V7+0japscTkZZrbAOdc96Bw/plzlo1zlo1zlU9Vz6bSYiIj4TsVFRER8p+JyYsYFHeA4lKt8lKt8lKt8qnQu9bmIiIjvdOQiIiK+U3EphZklmNnXZvatmS0zsz947a3MbL6ZrTGzyWZWPUJyTTCzDWaW7t26VGQuL0OMmS02s5ne40D3VQm5At9XXo6NZrbEy7DAa0sxs4+8ffaRmdWNkFy/N7OtRfbZTwPIVcfMpprZSjNbYWa9ImR/FZcr0P1lZu2KfHa6mR0ws3sqYn+puJQuB7jIOXcG0AW41MzOBv4KPOqcawPsBYZHSC6A/3HOdfFu6RWcC+BuYEWRx0Hvq0I/zgXB76tCF3oZCoeI3gfM9vbZbO9xJOSC0L9l4T57N4BM/wLed861B84g9G8aCfuruFwQ4P5yzq0q/GygG3AYeJMK2F8qLqVwIQe9h3HezQEXAVO99peAqyIkV6DMLA24DBjvPTYC3lfF5aoEriS0ryCgfRaJzCwZOA94HsA5d8w5t4+A91cJuSLJxcA659wmKmB/qbiUgXc6JR3YCXwErAP2OefyvE0ygKZB53LOzfee+rOZfWdmj5pZfAXHegz4DVDgPa5HBOyrYnIVCnJfFXLAh2a20MxGem0NnXOZAN7PBhGSC+AOb5+9EMDpp9ZAFvCid4pzvJnVJPj9dbxcEOz+Kup64DXvftj3l4pLGTjn8r3DyjSgB3BacZtVbKr/zmVmpwNjgPbAWUAK8NuKymNmA4CdzrmFRZuL2bRC99VxckGA++pH+jjnugL9gdvN7LyAcvxYcbmeBk4hdCo2E/hHBWeKBboCTzvnzgQOEdwpw6KOlyvo/QWA1895BfDvivpMFZdy8A5zPwXOBuqYWaz3VBqwLQJyXeqcy/ROmeUALxIqhhWlD3CFmW0EXid0Ouwxgt9X/5XLzCYGvK++55zb5v3cSeh8eA9gh5k1BvB+7oyEXM65Hd4fNQXAc1T8PssAMoocpU8l9Es96P1VbK4I2F+F+gOLnHM7vMdh318qLqUws/pmVse7XwPoS6ij7hPgWm+zm4G3IiDXyiL/wRih86hLKyqTc26Mcy7NOdeS0CH4x865mwh4Xx0n16Ag91UhM6tpZkmF94GfeDneJrSvIJj/vorNVbjPPFdTwfvMObcd2GJm7bymi4HlBLy/jpcr6P1VxA38cEoMKmB/xZa+SZXXGHjJzGIIFeMpzrmZZrYceN3MHgIW43XkRUCuj82sPqHTUenAqArOVZzfEuy+Op5JEbCvGgJvhuobscCrzrn3zewbYIqZDQc2A9dFSK5XLDRk2wEbgdsqOBfAnYT+7aoD64FheP8PBLi/jpfr8aD3l5klAv1+9NkPE+b9pSv0RUTEdzotJiIivlNxERER36m4iIiI71RcRETEdyouIiLiOxUXqRLMrJGZvW5m68xsuZm9a2Ztw/RZ7xZeg1QZmdkfzayvd/8ebyirSLloKLJEPe8iyTnAS865Z7y2LkCSc+6LMr7evKusC9tinHP54cocKbxZDbo753YFnUUqFx25SFVwIZBbWFgAnHPpzrkvzKyWmc02s0UWWrvkSgAza2mhNTnGAouAZmZ20Purfj7wgJm9Wfh+ZtbPzN7w7m80s9Qi7/Gchdbc+dCbTQEzO8ubzHCumf3dzIq9ctvMfuPl+tbMHvbaRpjZN17btMIjCwutT/OMmX1hZqu9OdUKv8sX3ndcZGa9S3n/CWZ2rZndBTQBPjGzT8xsuJk9WuS1I8zsn378A0kUcs7ppltU34C7CK2pUdxzsUCydz8VWEvoiv2WhGZQPrvItg74uXffgJVAfe/xq8Dl3v2N3nu1BPKALl77FGCQd38p0Nu7/zCwtJhs/QkdcSV6j1O8n/WKbPMQcKd3fwLwPqE/GtsQmu8qAUgEErxt2gALSnn/CcC1Rb+Ld78moRnB47zHc4BOQf/76haZNx25SFVnwF/M7DtgFqHlABp6z21yzs0rsm0+MA1C6+kArwCDvP6VXsB7xbz/BvfDImQLgZbe9knOuTle+6vHydYXeNE5d9j7zD1e++nekcgS4CagY5HXTHHOFTjn1hCagqQ9obV+nvO2/zfQoZT3L5Zz7hDwMTDAzNoTKjJLSnqNVF2aW0yqgmX8MHHmj90E1Ae6OedyvT6GBO+5Qz/a9qj7z36WF4EZwFHg3+6HNWuKyilyPx+oQfHLEBTHKH55ggnAVc65b81sKHBBked+vL0DfgnsILQ6YjUvb0nvX5LxwO8IHbW9WM7XShWiIxepCj4G4s1sRGGD1+dxPlCb0FovuWZ2IdCirG/qQlPSbwMeIPQLv6yv2wtk2w/LUl9/nE0/BG4p0qeS4rUnAZlmFkeoOBZ1nZlVM7NTCC1gtYrQd8x0oQEJg4GYUt6/qGzv8wqzzweaATfyn7PsivwHFReJet4prKuBft5Q5GXA7wkVhklAdzNbQOgX9cpyvv0kYItzbnk5XzccGGdmcwkdQewvJvf7hKZGX2ChFUd/7T31IDCf0KqoP867CviM0Cm6Uc65o8BY4GYzmwe0xTsiK+H9ixoHvGdmnxRpmwJ85RVJkWJpKLLISTCzJ4HFzrlyLSNgZrWccwe9+/cBjZ1zd59klgnATOfc1JN5nzJ8zkxCAyRmh/NzpHLTkYvICTKzhUBnYOIJvPwyM0v3hiCfS2jUV0Qzszpmtho4osIipdGRi4iI+E5HLiIi4jsVFxER8Z2Ki4iI+E7FRUREfKfiIiIivlNxERER3/1/U/OjdmEblJoAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "kappa = np.linspace(30, 70, 100)\n", "log_prob = [log_likelihood([y1[0], k, y1[2]]) for k in kappa]\n", "plt.plot(kappa, log_prob)\n", "plt.vlines(y1[1], ymin=-2700, ymax=log_likelihood(y1), linestyles='dashed')\n", "plt.xlabel('Carrying capacity')\n", "plt.ylabel('Log-likelihood')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To construct a profile likelihood $(1-\\alpha)$% confidence interval we determine the region in parameter space that satisfies:\n", "\n", "$$\\text{log } p(X|\\theta, \\hat{\\phi}) \\geq \\text{log } p(X|\\hat{\\theta}, \\hat{\\phi}) - \\frac{1}{2}\\chi(1)^2_{1-\\alpha},$$\n", "\n", "where $\\theta$ is the parameter we are seeking a confidence interval for and $\\phi$ is a vector of other parameters; the $(\\hat{\\theta},\\hat{\\phi})$ variables indicate the maximum likelihood estimates; and $\\chi(1)^2_{1-\\alpha}$ represents the $\\alpha$% critical values of a chi-squared distribution with one degree of freedom.\n", "\n", "First we obtain the threshold value of log-likelihood to construct a 95% confidence interval." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "-253.16865388426388\n" ] } ], "source": [ "import scipy.stats\n", "chi2 = scipy.stats.chi2.ppf(0.95, df=1)\n", "log_likelihood_min = log_likelihood(y1) - chi2 / 2\n", "print(log_likelihood_min)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next we find the bounds of this region, which yields the 95% confidence interval on this parameter." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Confidence interval = [48.977858109304364, 50.63527979747698]\n" ] } ], "source": [ "import scipy.optimize\n", "def log_likelihood_bounds(k):\n", " return (log_likelihood([y1[0], k, y1[2]]) - log_likelihood_min)**2\n", "res = scipy.optimize.minimize(log_likelihood_bounds, 40)\n", "kappa_min = res.x[0]\n", "res = scipy.optimize.minimize(log_likelihood_bounds, 60)\n", "kappa_max = res.x[0]\n", "\n", "print('Confidence interval = ' + str([kappa_min, kappa_max]))" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.8" } }, "nbformat": 4, "nbformat_minor": 2 }