{ "cells": [ { "cell_type": "markdown", "id": "6a017bc6-02c7-4313-bd81-bd811ef4965e", "metadata": {}, "source": [ "# Geostatistics revisited\n", "\n", "To attempt to cover geostatistics in any depth in roughly one hour would be impossible. Thus, I plan to introduce a few concepts you may find useful below and we'll see how much we can cover. These may be fairly basic geostatistical examples, but hopefully the provide some general sense of how Python can be used for statistical analyses. The concepts we will cover (briefly) are:\n", "\n", "- Least squares regressions (with and without weighting)\n", "- Linear correlation\n", "- Goodness-of-fit tests\n", "- Monte Carlo inversion" ] }, { "cell_type": "markdown", "id": "420a65c5-2660-4d19-9cec-8557cde5dc3a", "metadata": {}, "source": [ "## Some data\n", "\n", "For this exercise we can work with some hastily compiled climate data from a [global climate database](https://www.weatherbase.com/weather/countryall.php3). I've taken annual average, minimum, and maximum for 10 cities and put the data into [a csv file](data/global-climate-data.csv) that we can load and explore using the geostatistical concepts below. The 10 cities are largely from the northern hemisphere, and the goal here is to test the hypothesis that average annual temperatueres vary strongly as a function of latitude of the cities (a bold hypothesis, I know).\n", "\n", "Let's start by loading the data using pandas." ] }, { "cell_type": "code", "execution_count": 1, "id": "d13ec712-16e0-4ba1-994e-0ff7c54449da", "metadata": {}, "outputs": [], "source": [ "import pandas as pd\n", "\n", "# Read data file\n", "fp = \"data/global-climate-data.csv\"\n", "data = pd.read_csv(fp)" ] }, { "cell_type": "code", "execution_count": 2, "id": "2180cebd-d165-41a3-b34d-f493d8f8440d", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
CityLatitudeAverage temperatureMin temperatureMax temperature
0Helsinki60.195.01.08.0
1Paris48.4412.48.816.0
2Athens37.5418.514.322.0
3Cairo31.2421.016.027.0
4Chennai13.0028.524.232.9
5Kuala Lumpur3.0727.023.031.0
6Quito-0.0914.010.019.0
7Guayaquil-2.0926.022.030.0
8Ivalo68.361.0-5.02.0
9Kathmandu27.4218.013.022.0
\n", "
" ], "text/plain": [ " City Latitude Average temperature Min temperature \\\n", "0 Helsinki 60.19 5.0 1.0 \n", "1 Paris 48.44 12.4 8.8 \n", "2 Athens 37.54 18.5 14.3 \n", "3 Cairo 31.24 21.0 16.0 \n", "4 Chennai 13.00 28.5 24.2 \n", "5 Kuala Lumpur 3.07 27.0 23.0 \n", "6 Quito -0.09 14.0 10.0 \n", "7 Guayaquil -2.09 26.0 22.0 \n", "8 Ivalo 68.36 1.0 -5.0 \n", "9 Kathmandu 27.42 18.0 13.0 \n", "\n", " Max temperature \n", "0 8.0 \n", "1 16.0 \n", "2 22.0 \n", "3 27.0 \n", "4 32.9 \n", "5 31.0 \n", "6 19.0 \n", "7 30.0 \n", "8 2.0 \n", "9 22.0 " ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "data" ] }, { "cell_type": "code", "execution_count": 3, "id": "3af78b7d-c88b-441a-8584-35903b66589d", "metadata": {}, "outputs": [], "source": [ "# Calculate \"deviation\" by dividing difference between min and max temperatures by 2\n", "data[\"Deviation\"] = (data[\"Max temperature\"] - data[\"Min temperature\"]) / 2.0" ] }, { "cell_type": "code", "execution_count": 4, "id": "e5d501ee-9c8f-4666-aa2d-ffc1be510a22", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
CityLatitudeAverage temperatureMin temperatureMax temperatureDeviation
0Helsinki60.195.01.08.03.50
1Paris48.4412.48.816.03.60
2Athens37.5418.514.322.03.85
3Cairo31.2421.016.027.05.50
4Chennai13.0028.524.232.94.35
5Kuala Lumpur3.0727.023.031.04.00
6Quito-0.0914.010.019.04.50
7Guayaquil-2.0926.022.030.04.00
8Ivalo68.361.0-5.02.03.50
9Kathmandu27.4218.013.022.04.50
\n", "
" ], "text/plain": [ " City Latitude Average temperature Min temperature \\\n", "0 Helsinki 60.19 5.0 1.0 \n", "1 Paris 48.44 12.4 8.8 \n", "2 Athens 37.54 18.5 14.3 \n", "3 Cairo 31.24 21.0 16.0 \n", "4 Chennai 13.00 28.5 24.2 \n", "5 Kuala Lumpur 3.07 27.0 23.0 \n", "6 Quito -0.09 14.0 10.0 \n", "7 Guayaquil -2.09 26.0 22.0 \n", "8 Ivalo 68.36 1.0 -5.0 \n", "9 Kathmandu 27.42 18.0 13.0 \n", "\n", " Max temperature Deviation \n", "0 8.0 3.50 \n", "1 16.0 3.60 \n", "2 22.0 3.85 \n", "3 27.0 5.50 \n", "4 32.9 4.35 \n", "5 31.0 4.00 \n", "6 19.0 4.50 \n", "7 30.0 4.00 \n", "8 2.0 3.50 \n", "9 22.0 4.50 " ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Check DataFrame contents\n", "data" ] }, { "cell_type": "code", "execution_count": 5, "id": "24bf8978-3a9a-4845-beb9-1835e7267d36", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAs0AAAHgCAYAAABelVD0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAijElEQVR4nO3df5ClaVkf/O81uyA0gwFZmFpZpo8aQiQkLu4Urq74zixoVpOAJphIHRUjSZu8mpIApcS2SlC68ouoiSbRjhDWeMIEiT9WwACyOyAmJcwiwpLVYHinJ8iGlVJ+DI3grvf7R59ZetuZvnu2++lzpvfzqeo657nPc85z7VVdXd+95z73U621AAAAF3do1gUAAMC8E5oBAKBDaAYAgA6hGQAAOoRmAADoEJoBAKDjylkXsBNXXXVVG41Gsy7jAfvUpz6VRzziEbMu40DS22Ho63D0dhj6Ogx9HY7eDmMv+nr77bd/tLX22K3jl0VoHo1GOX369KzLeMBOnTqV48ePz7qMA0lvh6Gvw9HbYejrMPR1OHo7jL3oa1WtXWjc8gwAAOgQmgEAoENoBgCADqEZAAA6hGYAAOgQmgEAoENoBgCADqEZAAA6hGYAAOgQmgEAoENoBgCADqEZAAA6hGYAAOgQmgEAoENoBgCADqEZtjGZTDIajXLo0KGMRqNMJpNZlwQAzMCVsy4A5tVkMsnS0lLW19eTJGtra1laWkqSjMfjWZYGAOwzM81wEcvLy/cF5vPW19ezvLw8o4oAgFkRmuEizp49e0njAMDBJTTDRRw9evSSxgGAg0tohotYWVnJwsLC/cYWFhaysrIyo4oAgFkRmuEixuNxVldXs7i4mKrK4uJiVldXfQkQAB6E7J4B2xiPx0IyAGCmGQAAeoRmAADoEJoBAKBDaAYAgA6hGQAAOoRmAADoEJoBAKBDaAYAgA6hGQAAOoRmAADoEJoBAKBDaAYAgA6hGQAAOoTmA2IymWQ0GuXQoUMZjUaZTCazLgkA4MC4ctYFsHuTySRLS0tZX19PkqytrWVpaSlJMh6PZ1kaAMCBYKb5AFheXr4vMJ+3vr6e5eXlGVUEAHCwCM0HwNmzZy9pHACASyM0HwBHjx69pHEAAC6N0HwArKysZGFh4X5jCwsLWVlZmVFFAAAHi9B8AIzH46yurmZxcTFVlcXFxayurvoSIADAHrF7xgExHo+FZACAgZhpBgCADqEZAAA6hGYAAOgQmgEAoENoBgCADqEZAAA6hGYAAOgQmi9gMplkNBrl0KFDGY1GmUwmsy4JAIAZcnOTLSaTSZaWlrK+vp4kWVtby9LSUpK4eQgAwIOUmeYtlpeX7wvM562vr2d5eXlGFQEAMGtC8xZnz569pHEAAA4+oXmLo0ePXtI4AAAHn9C8xcrKShYWFu43trCwkJWVlRlVBADArAnNW4zH46yurmZxcTFVlcXFxayurvoSIADAg5jdMy5gPB4LyQAA3MdMMwAAdAjNAADQITQDAECH0AwAAB1CMwAAdAwWmqvqCVV1W1XdWVXvr6rvnY6/tKp+v6reM/35hqFqAACAvTDklnP3JHlRa+3dVfXIJLdX1Vumr/1Ya+0VA14bAAD2zGChubV2V5K7ps8/WVV3Jnn8UNcDAICh7Mua5qoaJXlqkt+cDn1PVb23ql5VVY/ejxoAAOCBqtbasBeoOpzkbUlWWmu/UFVHknw0SUvyI0mubq195wXet5RkKUmOHDly3cmTJwetc0jnzp3L4cOHZ13GgaS3w9DX4ejtMPR1GPo6HL0dxl709cSJE7e31o5tHR80NFfVQ5K8PsmbWms/eoHXR0le31p7ynafc+zYsXb69OlhitwHp06dyvHjx2ddxoGkt8PQ1+Ho7TD0dRj6Ohy9HcZe9LWqLhiah9w9o5K8MsmdmwNzVV296bRvSnLHUDUAAMBeGHL3jBuSfFuS91XVe6ZjP5DkuVV1bTaWZ5xJ8l0D1gAAALs25O4Z70hSF3jpjUNdEwAAhuCOgAAA0CE0AwBAh9AMAAAdQjMAAHQIzQAA0CE0AwBAh9AMAAAdQjMAAHQIzQAA0CE0AwBAh9AMAAAdQjMAAHQIzQAA0CE0AwBAh9AMAAAdQjMAAHQIzQAA0CE0AwBAh9AMAAAdQjMAAHQIzQAA0CE0AwBAh9AMAAAdQjMAAHQIzQAA0CE0AwBAh9AMAAAdQjMAAHQIzQAA0CE0AwBAh9AMAAAdQjMAAHQIzQAA0CE0AwBAh9AMAAAdQjMAAHQIzQAA0CE0AwBAh9AMAAAdQjMAAHQIzQAA0CE0AwBAh9AMAAAdQjMAAHQIzQAA0CE0AwBAh9AMAAAdQjMAAHQIzQAA0CE0AwBAh9AMMEcmk0lGo1EOHTqU0WiUyWQy65IASHLlrAsAYMNkMsnS0lLW19eTJGtra1laWkqSjMfjWZYG8KBnphlgTiwvL98XmM9bX1/P8vLyjCoC4DyhGWBOnD179pLGAdg/QjPAnDh69OgljQOwf4RmgDmxsrKShYWF+40tLCxkZWVlRhUBcJ7QDDAnxuNxVldXs7i4mKrK4uJiVldXfQkQYA7YPQNgjozHYyEZYA6ZaQYAgA6hGQAAOoRmAADoEJoBAKBDaAYAgI7BQnNVPaGqbquqO6vq/VX1vdPxL6iqt1TVB6aPjx6qBgAA2AtDzjTfk+RFrbUvTXJ9ku+uqicneUmSt7bWnpjkrdNjAACYW4OF5tbaXa21d0+ffzLJnUken+TZSW6ennZzkm8cqgYAANgL+7KmuapGSZ6a5DeTHGmt3ZVsBOskj9uPGgAA4IGq1tqwF6g6nORtSVZaa79QVR9rrT1q0+t/1Fr7M+uaq2opyVKSHDly5LqTJ08OWueQzp07l8OHD8+6jANJb4ehr8PR22Ho6zD0dTh6O4y96OuJEydub60d2zo+6G20q+ohSf5rkklr7Remwx+pqqtba3dV1dVJ7r7Qe1trq0lWk+TYsWPt+PHjQ5Y6qFOnTuVyrn+e6e0w9HU4ejsMfR2Gvg5Hb4cxZF+H3D2jkrwyyZ2ttR/d9NItSZ43ff68JL88VA0AALAXhpxpviHJtyV5X1W9Zzr2A0n+WZLXVtXzk5xN8s0D1gAAALs2WGhurb0jSV3k5WcMdV0AANhr7ggIAAAdQjMAAHQIzQAA0CE0AwBAh9AMAAAdQjMAAHQIzQAA0CE0AwBAh9AMAAAdQjMAAHQIzQAA0CE0AwBAh9AMAAAdQjMA+2IymWQ0GuXQoUMZjUaZTCazLglgx66cdQEAHHyTySRLS0tZX19PkqytrWVpaSlJMh6PZ1kawI6YaQZgcMvLy/cF5vPW19ezvLw8o4oALo3QDMDgzp49e0njAPNGaAZgcEePHr2kcYB5IzQDMLiVlZUsLCzcb2xhYSErKyszqgjg0gjNAAxuPB5ndXU1i4uLqaosLi5mdXXVlwCBy4bdMwDYF+PxWEgGLltmmgGm7CMMwMWYaQaIfYQB2J6ZZoDYRxiA7QnNALGPMADbE5oBYh9hALYnNAPEPsIAbE9oBoh9hAHYnt0zAKbsIwzAxZhpBgCAjm1nmqvqK5N8a5KnJ7k6yaeT3JHkDUl+rrX28cErBACAGbvoTHNV/WqSv5fkTUluykZofnKSH0zysCS/XFXP2o8iAQBglrabaf621tpHt4ydS/Lu6c+/qqqrBqsMAADmxHZrmh9VVTdsHayqp1fVlyTJBUI1AAAcONuF5h9P8skLjH96+hoAADwobBeaR621924dbK2dTjIarCIAAJgz24Xmh23z2sP3uhAAAJhX24Xmd1XV3986WFXPT3L7cCUBAMB82W73jBck+cWqGudzIflYkocm+aaB6wIAgLlx0dDcWvtIkq+qqhNJnjIdfkNr7dZ9qQwAAObEtncETJLW2m1JbtuHWgAAYC5td0fAb66qX6qqX6yqv7OfRQEAwDzZbqb5+5M8bfr8XUn+y/DlAADA/NkuNP9ckp+dPv/5fagFAADm0nZfBPzxqnpEkmqtndvHmgAAYK5cNDRXVbXWPrXdm6fntL0vCwAA5sd2Nze5rar+UVUd3TxYVQ+tqhur6uYkzxu2PAAAmL3t1jTflOQ7k7ymqr4oyceycWvtK5K8OcmPtdbeM3SBAAAwa9utaf7jJP8uyb+rqockuSrJp1trH9un2gAAYC50b26SJK21P0ly18C1AADAXNpuTTMAABChGQAAunYUmqtqsaqeOX3+8Kp65LBlAQDA/OiG5qr6+0lel+Snp0PXJPmlAWsCAIC5spOZ5u9OckOSTyRJa+0DSR43ZFEAADBPdhKaP9Na++z5g6q6Mom7AAIA8KCxk9D8tqr6gSQPr6qvTfLzSX5l2LIAAGB+7CQ0f3+SP0jyviTfleSNSX5wyKIAAGCebHtzk6o6lOS9rbWnJPkP+1MSAADMl21nmltrf5rkt6vq6D7VAwAAc2cnt9G+Osn7q+qdST51frC19qzBqgIAgDmyk9D8ssGrAACAOdYNza21t+1HIQAAMK92ckfAT1bVJ6Y/f1xV91bVJ3bwvldV1d1VdcemsZdW1e9X1XumP9+w2/8AAAAY2k5mmh+5+biqvjHJ03bw2a9O8pNJfnbL+I+11l6xw/oAAGDmdrJP8/201n4pyY07OO/tSf7wAdQEAABzpTvTXFV/c9PhoSTHsrvbaH9PVX17ktNJXtRa+6NdfBYAAAyuWts+/1bVf9x0eE+SM0n+Q2vt7u6HV42SvH56c5RU1ZEkH81G6P6RJFe31r7zIu9dSrKUJEeOHLnu5MmTvcvNrXPnzuXw4cOzLuNA0tth6Otw9HYY+joMfR2O3g5jL/p64sSJ21trx7aO7yQ039Ba+43e2EXeO8qm0LzT17Y6duxYO336dO+0uXXq1KkcP3581mUcSHo7DH0djt4OQ1+Hoa/D0dth7EVfq+qCoXkna5p/YodjOyni6k2H35TkjoudCwAA8+Kia5qr6iuTfFWSx1bVCze99PlJruh9cFW9JsnxJFdV1YeS/FCS41V1bTaWZ5xJ8l0PtHAAANgv230R8KFJDk/P2bzt3CeSPKf3wa21515g+JWXVB0AAMyBi4bm6Z0A31ZVr26tre1jTQAAMFe6W84lWa+qf5nkLyV52PnB1lp3r2YAADgIdvJFwEmS30nyRUlelo21yO8asCYAAJgrOwnNj2mtvTLJn7TW3jbdV/n6gesCAIC5sZPlGX8yfbyrqv5akg8nuWa4kgAAYL7sJDS/vKr+XJIXZWN/5s9P8o8HrQoAAObItsszquqKJE9srX28tXZHa+1Ea+261tot+1TfZW0ymWQ0GuXGG2/MaDTKZDKZdUkAADwA24bm1tq9SZ61T7UcKJPJJEtLS1lbW0trLWtra1laWhKcAQAuQzv5IuB/r6qfrKqnV9WXn/8ZvLLL3PLyctbX1+83tr6+nuXl5RlVBADAA7WTNc1fNX384U1jLYl9mrdx9uzZSxoHAGB+dUNza+3EfhRy0Bw9ejRra3/2RopHjx6dQTUAAOxGd3lGVR2pqldW1a9Oj59cVc8fvrTL28rKShYWFu43trCwkJWVlRlVBADAA7WTNc2vTvKmJF84Pf5fSV4wUD0Hxng8zurqahYXF1NVWVxczOrqasbj8axLAwDgEu0kNF/VWnttkj9NktbaPUnuHbSqA2I8HufMmTO59dZbc+bMGYEZAOAytZPQ/Kmqekw2vvyXqro+yccHrQoAAObITnbPeGGSW5J8SVX9RpLHJnnOoFUBAMAc2cnuGe+uqv8nyZOSVJLfba39yeCVAQDAnOiG5qp6WJL/N8lXZ2OJxq9X1U+11v546OIAAGAe7GR5xs8m+WSSn5gePzfJf0ryzUMVBQAA82QnoflJrbUv23R8W1X99lAFAQDAvNnJ7hm/Nd0xI0lSVV+R5DeGKwkAAObLTmaavyLJt1fV2enx0SR3VtX7krTW2l8ZrDoAAJgDOwnNNw1eBQAAzLGdbDm3VlWPTvKEzee31t49ZGEAADAvdrLl3I8k+Y4k/zvTuwJOH28criwAAJgfO1me8beTfElr7bNDFwMAAPNoJ7tn3JHkUQPXAQAAc2snM83/NBvbzt2R5DPnB1trzxqsKgAAmCM7Cc03J/nnSd6X5E+HLQcAAObPTpZnfLS19m9aa7e11t52/mfwygCAXZtMJhmNRjl06FBGo1Emk8msS4LL0k5mmm+vqn+a5Jbcf3mGLecAYI5NJpMsLS1lfX09SbK2tpalpaUkyXg8nmVpcNnZSWh+6vTx+k1jtpwDgDm3vLx8X2A+b319PcvLy0IzXKKd3NzkxH4UAgDsrbNnz17SOHBx3TXNVXWkql5ZVb86PX5yVT1/+NIAgN04evToJY0DF7eTLwK+Osmbknzh9Ph/JXnBQPUAAHtkZWUlCwsL9xtbWFjIysrKjCqCy9dFQ3NVnV+6cVVr7bWZbjfXWrsnyb37UBsAsAvj8Tirq6tZXFxMVWVxcTGrq6vWM8MDsN2a5ncm+fIkn6qqx2Tjy3+pquuTfHwfagMAdmk8HgvJsAe2C801fXxhNrab+5Kq+o0kj03ynKELAwCAebFdaH5sVb1w+vwXk7wxG0H6M0memeS9A9cGAABzYbvQfEWSw/ncjPN5Cxc4FwAADqztQvNdrbUf3rdKAABgTm235dzWGWYAAHhQ2i40P2PfqgAAgDl20dDcWvvD/SwEAADm1U7uCAgAAA9qQjMAAHQIzQAA0CE0AwBAh9AMAAAdQjMAAHQIzQAA0CE0AwBAh9AMAAAdQjMAAHQIzQAA0CE0AwBAh9AMAAAdQjMAAHQIzQAA0CE0AwBAh9AMAAAdQjMAAHQIzQAA0DFYaK6qV1XV3VV1x6axL6iqt1TVB6aPjx7q+gAAsFeGnGl+dZKbtoy9JMlbW2tPTPLW6TEAAMy1wUJza+3tSf5wy/Czk9w8fX5zkm8c6voAALBX9ntN85HW2l1JMn183D5fHwAALlm11ob78KpRkte31p4yPf5Ya+1Rm17/o9baBdc1V9VSkqUkOXLkyHUnT54crM6hnTt3LocPH551GQeS3g5DX4ejt8PQ12Ho63D0dhh70dcTJ07c3lo7tnX8yl196qX7SFVd3Vq7q6quTnL3xU5sra0mWU2SY8eOtePHj+9TiXvv1KlTuZzrn2d6Owx9HY7eDkNfh6Gvw9HbYQzZ1/1ennFLkudNnz8vyS/v8/UBAOCSDbnl3GuS/I8kT6qqD1XV85P8syRfW1UfSPK102MAAJhrgy3PaK099yIvPWOoawIAwBDcERAAADqEZgAA6BCaAQCgQ2gGAIAOoRkAADqEZgAA6BCaAQCgQ2gGAIAOoRkAADqEZgAA6BCaAQCgQ2gGAIAOoRkAADqEZgAA6BCaAQCgQ2gGAIAOoRkAADqEZgAA6BCaAQCgQ2gGAIAOoRkAADqEZgAA6BCaAQCgQ2gGAIAOoRkAADqEZgAA6BCaAQCgQ2gGAIAOoRkAADqEZgAA6BCaAQCgQ2gGAIAOoRkAADqEZgAA6BCaAQCgQ2gGAIAOoRkAADqEZgAA6BCaAQCgQ2gGAIAOoRkAADqEZgAA6BCaAQCgQ2gGAIAOoRkAADqEZgAA6BCaAQCgQ2gGAIAOoRkAADqEZgAA6BCaAQCgQ2gGAIAOoRkAADqEZgAA6BCaAQCgQ2gGAIAOoRkAADqEZgAA6BCaAQCgQ2gGAIAOoRkAADqEZgAA6BCaAQCg48pZXLSqziT5ZJJ7k9zTWjs2izoAAGAnZjnTfKK1dq3ADAA82Pzar/1aRqNRDh06lNFolMlkMuuS6JjJTDMAwIPVZDLJK17xinzmM59JkqytrWVpaSlJMh6PZ1ka25jVTHNL8uaqur2qlmZUAwDAvlteXr4vMJ+3vr6e5eXlGVXETlRrbf8vWvWFrbUPV9XjkrwlyT9qrb19yzlLSZaS5MiRI9edPHly3+vcK+fOncvhw4dnXcaBpLfD0Nfh6O0w9HUY+jqMG2+8MRfKX1WVW2+9dQYVHRx78Tt74sSJ2y+0fHgmofl+BVS9NMm51torLnbOsWPH2unTp/evqD126tSpHD9+fNZlHEh6Owx9HY7eDkNfh6GvwxiNRllbW/sz44uLizlz5sz+F3SA7MXvbFVdMDTv+/KMqnpEVT3y/PMkX5fkjv2uAwBgFlZWVvJ5n/d59xtbWFjIysrKjCpiJ2axpvlIkndU1W8neWeSN7TW/tsM6gAA2Hfj8TgvfvGLs7i4mKrK4uJiVldXfQlwzu377hmttQ8m+bL9vi4AwLx45jOfmZe//OWzLoNL4I6AAADQITQDAECH0AwAAB1CMwAAdAjNAADQITQDAECH0AwAAB1CMwAAdAjNAADQITQDAECH0AwAAB1CMwAAdAjNAADQITQDAECH0AwAAB1CMwAAdAjNAADQITQDAECH0AwAAB1CMwAAdAjNAADQITQDAECH0AwAAB1CMwAAdAjNAADQITQDAECH0AwAAB1CMwAAdAjNAADQITQDAECH0AwAAB1CMwAAdAjNAADQITQDAECH0AwAAB1CMwAAdAjNAADQITQDAECH0AwAAB1CMwAAdAjNAADQITQDAECH0AwAAB1CMwAAdAjNAADQITQDAECH0AwAAB1CMwAAdAjNAADQITQDADAXJpNJRqNRDh06lNFolMlkMuuS7nPlrAsAAIDJZJKlpaWsr68nSdbW1rK0tJQkGY/HsywtiZlmAADmwPLy8n2B+bz19fUsLy/PqKL7E5oBAJi5s2fPXtL4fhOaAQCYuaNHj17S+H4TmgEAmLmVlZUsLCzcb2xhYSErKyszquj+hGYAAGZuPB5ndXU1i4uLqaosLi5mdXV1Lr4EmNg9AwCAOTEej+cmJG9lphkAADqEZgAA6BCaAQCgQ2gGAIAOoRkAADqEZgAA6JhJaK6qm6rqd6vq96rqJbOoAQAAdmrfQ3NVXZHk3yb5+iRPTvLcqnryftcBAAA7NYuZ5qcl+b3W2gdba59NcjLJs2dQBwAA7Ei11vb3glXPSXJTa+3vTY+/LclXtNa+Z8t5S0mWkuTIkSPXnTx5cl/r3Evnzp3L4cOHZ13GgaS3w9DX4ejtMPR1GPo6HL0dxl709cSJE7e31o5tHZ/FbbTrAmN/Jrm31laTrCbJsWPH2vHjxwcuazinTp3K5Vz/PNPbYejrcPR2GPo6DH0djt4OY8i+zmJ5xoeSPGHT8TVJPjyDOgAAYEdmEZrfleSJVfVFVfXQJN+S5JYZ1AEAADuy78szWmv3VNX3JHlTkiuSvKq19v79rgMAAHZqFmua01p7Y5I3zuLaAABwqfZ994wHoqr+IMnarOvYhauSfHTWRRxQejsMfR2O3g5DX4ehr8PR22HsRV8XW2uP3Tp4WYTmy11Vnb7Q1iXsnt4OQ1+Ho7fD0Ndh6Otw9HYYQ/Z1JrfRBgCAy4nQDAAAHULz/liddQEHmN4OQ1+Ho7fD0Ndh6Otw9HYYg/XVmmYAAOgw0wwAAB1C88Cq6qaq+t2q+r2qesms67mcVdWrquruqrpj09gXVNVbquoD08dHz7LGy1FVPaGqbquqO6vq/VX1vdNxvd2FqnpYVb2zqn572teXTcf1dQ9U1RVV9VtV9frpsb7ugao6U1Xvq6r3VNXp6Zje7lJVPaqqXldVvzP9W/uV+rp7VfWk6e/q+Z9PVNULhuqt0Dygqroiyb9N8vVJnpzkuVX15NlWdVl7dZKbtoy9JMlbW2tPTPLW6TGX5p4kL2qtfWmS65N89/T3VG935zNJbmytfVmSa5PcVFXXR1/3yvcmuXPTsb7unROttWs3bdult7v3r5P8t9baX0zyZdn43dXXXWqt/e70d/XaJNclWU/yixmot0LzsJ6W5Pdaax9srX02yckkz55xTZet1trbk/zhluFnJ7l5+vzmJN+4nzUdBK21u1pr754+/2Q2/pg/Pnq7K23DuenhQ6Y/Lfq6a1V1TZK/luRnNg3r63D0dheq6vOTfE2SVyZJa+2zrbWPRV/32jOS/O/W2loG6q3QPKzHJ/k/m44/NB1j7xxprd2VbIS/JI+bcT2XtaoaJXlqkt+M3u7adAnBe5LcneQtrTV93Rs/nuT7kvzppjF93RstyZur6vaqWpqO6e3ufHGSP0jyH6dLin6mqh4Rfd1r35LkNdPng/RWaB5WXWDMdiXMpao6nOS/JnlBa+0Ts67nIGit3Tv9Z8Nrkjytqp4y45Iue1X115Pc3Vq7fda1HFA3tNa+PBvLCr+7qr5m1gUdAFcm+fIk/7619tQkn4qlGHuqqh6a5FlJfn7I6wjNw/pQkidsOr4myYdnVMtB9ZGqujpJpo93z7iey1JVPSQbgXnSWvuF6bDe7pHpP8WeysaafH3dnRuSPKuqzmRjyduNVfVz0dc90Vr78PTx7mysDX1a9Ha3PpTkQ9N/aUqS12UjROvr3vn6JO9urX1kejxIb4XmYb0ryROr6oum/xf0LUlumXFNB80tSZ43ff68JL88w1ouS1VV2Vhrd2dr7Uc3vaS3u1BVj62qR02fPzzJM5P8TvR1V1pr/6S1dk1rbZSNv6m3tta+Nfq6a1X1iKp65PnnSb4uyR3R211prf3fJP+nqp40HXpGkv8Zfd1Lz83nlmYkA/XWzU0GVlXfkI31d1ckeVVrbWW2FV2+quo1SY4nuSrJR5L8UJJfSvLaJEeTnE3yza21rV8WZBtV9dVJfj3J+/K5NaI/kI11zXr7AFXVX8nGF1CuyMYExWtbaz9cVY+Jvu6Jqjqe5MWttb+ur7tXVV+cjdnlZGNJwX9ura3o7e5V1bXZ+OLqQ5N8MMnfzfTvQvR1V6pqIRvfH/vi1trHp2OD/M4KzQAA0GF5BgAAdAjNAADQITQDAECH0AwAAB1CMwAAdAjNAADQITQDDKSqzl3Cucer6qs2Hf+Dqvr26fPvqKovfADXP1NVV13ie1433a/3fE2nq+pfbHr9VFWd3nR8rKpOTZ//5ap69aXWCXA5EJoB5sPxJPeF5tbaT7XWfnZ6+B1JLjk0X6qq+ktJrmitfXA69A+TPD3JFVX1Fzed+riq+vqt72+tvS/JNVV1dOhaAfab0Aywj6rqb1TVb1bVb1XVr1XVkaoaJfkHSf5xVb2nqp5eVS+tqhdX1XOSHEsymb728M0zyFtmeh9TVW+efvZPJ6lN1/3Wqnrn9DN+uqquuEB549z/drOHkrRs3CmyNo3/yyQ/eJH/xF/Jxu2tAQ4UoRlgf70jyfWttacmOZnk+1prZ5L8VJIfa61d21r79fMnt9Zel+R0kvH0tU9v89k/lOQd08++JRu3kE1VfWmSv5PkhtbatUnuzUZA3uqGJLdvOv6ZJP89yaHW2p2bxv9Hks9U1YkLfMbpbMxOAxwoV866AIAHmWuS/JequjrJQ5P8f3v42V+T5G8mSWvtDVX1R9PxZyS5Lsm7qipJHp7k7gu8/+okf3D+oLX2piRvusi1Xp6N2ebv3zJ+d/ZhKQnAfjPTDLC/fiLJT7bW/nKS70rysAfwGffkc3+/t76/XeD8SnLzdKb62tbak1prL73AeZ/eaT2ttVun516/5aWHTT8H4EARmgH2159L8vvT58/bNP7JJI+8yHu2vnYmGzPHSfK3No2/PdNlF9Mv6j16Ov7WJM+pqsdNX/uCqlq8wHXuTPLnd/RfsWElyfdtGfsLSe64hM8AuCwIzQDDWaiqD236eWGSlyb5+ar69SQf3XTuryT5pvNfBNzyOa9O8lPnvwiY5GVJ/vX0M+7ddN7LknxNVb07ydclOZskrbX/mY2lFG+uqvcmeUs2lmJs9YZs7OKxI621N2bTco6pE9PPAThQqrUL/UseAA8200B+Wza+MHhv7/wLvP/zkrwtyVe31u7Z6/oAZkloBuA+VfVXk9zZWjv7AN77xCSPb62d2vPCAGZMaAYAgA5rmgEAoENoBgCADqEZAAA6hGYAAOgQmgEAoOP/B+PHOKIHacAgAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Create a scatter plot of the average temperatures as a function of latitude\n", "import matplotlib.pyplot as plt\n", "\n", "# Create figure, axes\n", "fig, ax = plt.subplots(1, 1, figsize=(12,8))\n", "\n", "# Plot temperature data\n", "ax.plot(data[\"Latitude\"], data[\"Average temperature\"], 'ko')\n", "\n", "# Add axis labels\n", "ax.set_xlabel(\"Latitude (°N)\")\n", "ax.set_ylabel(\"Temperature (°C)\")\n", "\n", "# Enable plot grid\n", "ax.grid()" ] }, { "cell_type": "markdown", "id": "0870b458-1902-45d7-b551-40f41f37bd16", "metadata": {}, "source": [ "## Least squares regressions\n", "\n", "Least squares regressions can be used to fit a function to a scattered set of data points on an x-y plot. Typically we are fitting a line to the data, so we will use that example here as well.\n", "\n", "### Unweighted least squares regressions\n", "\n", "We can start with the idea that we want a line that fits a scattered set of x and y points. As you may recall, the equation for a line is\n", "\n", "$$\n", "\\large\n", "y = A + B x.\n", "$$\n", "\n", "Thus, to plot the line we would need to be able to find a slope $B$ and y-intercept $A$. We can find those values from our scatted x and y points using the following:\n", "\n", "$$\n", "\\large\n", "A = \\frac{\\sum{x^{2}} \\sum{y} - \\sum{x} \\sum{xy}}{\\Delta}\n", "$$\n", "\n", "where $x$ is the $i$th data point plotted on the $x$-axis, $y$ is the $i$th data point plotted on the $y$-axis, and $\\Delta$ is defined below.\n", "\n", "The line slope can be found using\n", "\n", "$$\n", "\\large\n", "B = \\frac{N \\sum{xy} - \\sum{x} \\sum {y}}{\\Delta}\n", "$$\n", "\n", "where $N$ is the number of values in the regression.\n", "\n", "And the value of $\\Delta$ is\n", "\n", "$$\n", "\\large\n", "\\Delta = N \\sum{x^{2}} - \\left( \\sum{x} \\right)^{2}\n", "$$\n", "\n", "Let's test this out for a set of x and y points below." ] }, { "cell_type": "code", "execution_count": 6, "id": "ba81cbda-f23c-4323-b2fb-0e2fb63424b1", "metadata": {}, "outputs": [], "source": [ "def calculate_delta(x_values):\n", " \"\"\"Calculates divisor term for slope and intercept equations\"\"\"\n", " x_squared_sum = (x_values * x_values).sum()\n", " x_sum_squared = x_values.sum()**2\n", " delta = len(x_values) * x_squared_sum - x_sum_squared\n", " return delta\n", "\n", "def calculate_a(x_values, y_values):\n", " \"\"\"Calcuates y-intercept for regression line\"\"\"\n", " left_side = (x_values * x_values).sum() * y_values.sum()\n", " right_side = x_values.sum() * (x_values * y_values).sum()\n", " delta = calculate_delta(x_values)\n", " a = (left_side - right_side) / delta\n", " return a\n", " \n", "def calculate_b(x_values, y_values):\n", " \"\"\"Calculates slope for regression line\"\"\"\n", " left_side = len(x_values) * (x_values * y_values).sum()\n", " right_side = x_values.sum() * y_values.sum()\n", " delta = calculate_delta(x_values)\n", " b = (left_side - right_side) / delta\n", " return b" ] }, { "cell_type": "code", "execution_count": 7, "id": "f9a51184-efaf-4929-9f72-b24e40483e25", "metadata": {}, "outputs": [], "source": [ "# Calculate slope and intercept\n", "a = calculate_a(data[\"Latitude\"].values, data[\"Average temperature\"].values)\n", "b = calculate_b(data[\"Latitude\"].values, data[\"Average temperature\"].values)" ] }, { "cell_type": "code", "execution_count": 8, "id": "dd0b5283-c7a0-470c-bab1-eb5772e17399", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAs0AAAHgCAYAAABelVD0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABO2UlEQVR4nO3dd3hVVbrH8e+iakBFBRmKEFRQUREwIopKEAv2rmh07LGPqGOZQb22OHadsccyWDIwig17JVhRQbBgVzrWURwxI4Ku+8c+ocNJSE5OyvfzPOfJOfu0N+/N6O9u135XiDEiSZIkafkaZbsASZIkqbYzNEuSJElpGJolSZKkNAzNkiRJUhqGZkmSJCkNQ7MkSZKURpNsF1ARrVu3jrm5udkuo1b6+eefadGiRbbLqNPsYdXZw6qzh9XDPladPaw6e1h12ezh+PHjv4sxtlnyeJ0Izbm5uYwbNy7bZdRKpaWl5OfnZ7uMOs0eVp09rDp7WD3sY9XZw6qzh1WXzR6GEKYu67jLMyRJkqQ0DM2SJElSGoZmSZIkKQ1DsyRJkpSGoVmSJElKw9AsSZIkpWFoliRJktIwNEuSJElpGJolSZKkNAzNkiRJUhqGZkmSJCkNQ7MkSZKUhqFZkiRJSsPQLEmSJKVhaJYkSZLSMDRLkiRJaRiaJUmSpDQMzZIkSVIahmYpg/Lz88nPz892GZIkqYoMzZIkSVIahmZJkiQpDUOzJEmSlIahWZIkSUrD0CxJkiSlYWiWJEmS0jA0S5IkSWkYmiVJkqQ0DM2SJElSGoZmSZIkKQ1DsyRJkpSGoVmSJElKw9AsSZIkpWFoliRJktIwNEuSJElpGJolSZKkNAzNkiRJUhqGZkmSJCkNQ7MkSZKUhqFZK5Sfn09+fn62y5AkScoqQ7MkSZKUhqFZkiRJSsPQLEmSJKVhaJYkSZLSMDRLkiRJaRiaJUmSpDQMzZIkSVIaGQvNIYRVQghvhhDeCSFMCiFclDq+VgjhuRDCp6mfa2aqBkmSJKk6ZPJM81xghxjj5kBPYFAIoS9wLvBCjLEr8ELqsSRJklRrZSw0x8Sc1MOmqVsE9gbuTh2/G9gnUzVIkiRJ1SGja5pDCI1DCBOBb4DnYoxvAG1jjF8CpH6uk8kaJEmSpKpqkskPjzH+BvQMIbQCHg4hbFrR94YQCoFCgLZt21JaWpqRGrNlyJAhAFx//fVV+pw5c+ZktDezZ88GqHf9X1Qme9gQ+geZ/ztsCOxh9bCPVWcPq84eVl1t7GFGQ3O5GOPsEEIpMAj4OoTQLsb4ZQihHclZ6GW9pxgoBsjLy4v5+fk1UWqNadWqFQBV/b1KS0ur/BkrUl111maZ7GFD6B9k/u+wIbCH1cM+Vp09rDp7WHW1sYeZnJ7RJnWGmRDCqsCOwEfAKOCI1MuOAB7NVA2SJElSdcjkmeZ2wN0hhMYk4fz+GOPjIYTXgftDCMcA04ADM1iDJEmSVGUZC80xxneBXss4/h9gYKa+V5IkSapu7ggoSZIkpWFoliRJktIwNEuSJElpGJolSZKkNAzNkiRJUhqGZkmSJCkNQ7MkSZKUhqFZkiRJSsPQLGVISUkJY8eOZcyYMeTm5lJSUpLtkiRJ0koyNEsZUFJSQmFhIXPnzgVg6tSpFBYWGpwlSaqjDM1SBgwdOpSysrLFjpWVlTF06NAsVSRJkqrC0CxlwLRp0yp1XJIk1W6GZikDOnXqVKnjkiSpdjM0SxlQVFRETk7OYsdycnIoKirKUkWSJKkqDM1SBhQUFFBcXEzz5s0B6Ny5M8XFxRQUFGS5MkmStDKaZLsAqb4qKCjg9ttvB6C0tDS7xUiSpCrxTLMkSZKUhqFZkiRJSsPQLEmSJKVhaJYkSZLSMDRLkiRJaRiaJUmSpDQMzZIkSVIahmZJkiQpDUPz8syfDz//nO0qJEmSVAsYmpfn8cehfXv405/gww+zXU1WlJSUMHbsWMaMGUNubi4lJSXZLkmSJCkrDM3Ls956sOeecNtt0L07DBgA998Pv/6a7cpqRElJCYWFhcydOxeAqVOnUlhYaHCWJEkNkqF5eXr0gPvugxkz4PLLYcoUOPhg6NQJzj8fpk/PdoUZNXToUMrKyhY7VlZWxtChQ7NUkSRJUvYYmtNp0wbOOQc++wyeeAK23BKKiiA3F/beG555Bn7/PdtVVrtp06ZV6rgkSVJ9ZmiuqMaNYbfd4LHH4IsvkiD9+uswaBB06wZXXw3/+U+2q6w2nTp1qtRxSZKk+szQvDJyc+Gyy5IlGv/6V3LB4FlnQYcO8Mc/wtixEGO2q6ySoqIicnJyFjuWk5NDUVFRliqSJEnKHkNzVTRvDoccAi+9BO+9B8ccA488AltvDb17Q3ExzJmT7SpXSkFBAcXFxTRv3hyAzp07U1xcTEFBQZYrkyRJqnmG5uqy6aZw000wcybcemuyzvn445Ozz6eeCh98kO0KK62goIC+ffvSv39/pkyZYmCWJEkNlqG5uq22WhKWJ06EV1+FvfZKzjhvsgnk58O//91gxtZJkiTVF4bmTAkBttkG7r03GVt3xRUwbRoMHgydOnHM5Mms88sv2a5SUh2Tn59Pfn5+tsuQpAbH0FwT2rSBs89OxtY9+ST06cOh06Yx/I03krF1Tz9dL8fWSZIk1ReG5prUqBHsuislBx/M6s2a0QTIffxxSnbdFbp2hauugu++y3aVkiRJWoKhuYaVb09dllrXPPX33yls1oySpk2Ts9EdO8LhhyczoOv42DpJkqT6wtBcw5a5PfWvvzL0l1/g/ffh2GPh0UeT9dC9etXpsXWSJEn1haG5hq1we+pNNoEbb4RZs+C225Injj8+2TzllFNg0qQarFSSJEnlDM01rELbU7dsCYWFMGECvPYa7LMP3H57Mgu6f38YMcKxdZIkSTXI0FzDKrU9dQjJ7oL33JNsmnLllcn4ukMOgXXXhaFDaf7VVzVUuSRJUsNlaK5hK709devWcNZZ8Omn8NRT0LcvXH45fQsKkg1UnnrKsXWSJEkZYmjOgiptT92oEQwalFwsOHky0w49FN58E3bbLRlbd+WVjq2TJEmqZobmuqxTJyYfc0yy0+CIEcmSjXPOgQ4dkrF1r73m2DpJkqRqYGiuD5o1g4MPhtLSZGxdYSGMGgX9+kHPnskkDsfWSZIkrTRDc32zySZwww3JhYPFxclyjhNOSMbWnXxyEqolSZJUKRkLzSGEdUMIo0MIH4YQJoUQTksdvzCEMDOEMDF12y1TNTRoLVvCccfB228nuwvusw/ceSdsthlsv71j6yRJkiohk2ea5wNnxhg3BvoCJ4cQuqeeuy7G2DN1ezKDNSiEZNLGPfck4+quuio5C10+tu6vf4WpU7NdpSRJUq2WsdAcY/wyxvh26v5PwIdAh0x9nyqgdWv485+TsXVPP53MgL7iCujSBfbcE558En77LdtVSpIk1TpNauJLQgi5QC/gDaAfcEoI4Y/AOJKz0T8s4z2FQCFA27ZtKS0trYlSa8zs2bMBqvx7zZkzZ+U+o3lzGDKE5occQrsnnqD944/T7PHH+V+7dszac0++2nVX5rVqVW111mYr3cMKaAj9g8z2sKGoaA8byt/UyvJvsersYdXZw6qrjT0MMcMjyUIILYExQFGM8aEQQlvgOyAClwDtYoxHr+gz8vLy4rhx4zJaZ03Lz88Hqv4vvtLS0gWfVSW//gqPPAK33JJM4WjWDA48kFMmTeL91VendMyYqn9HLVVtPVyG6vq/c22XyR42FBXtYUP5m1pZ/i1WnT2sOntYddnsYQhhfIwxb8njGZ2eEUJoCjwIlMQYHwKIMX4dY/wtxvg7cDvQJ5M1qIKaNYODDoLRo2HSJDj+eHjsMW6cOJE7xo+HW2+Fn37KdpWSJElZkcnpGQG4E/gwxnjtIsfbLfKyfQFnoNU23bvDP/4Bs2ZxVbdu/B4CnHhismmKY+skSVIDlMkzzf2Aw4Edlhgvd2UI4b0QwrvAAOD0DNagqmjRgifataOwd28YOxb23Xfh2LrttoPhw2Hu3GxXKUmSlHEZuxAwxvgKEJbxlCPm6poQYKutktu118KwYcna50MPhTZt4JhjkuUcubnZrlSSJCkj3BFQlbP22nDmmfDJJ/DMM7DNNnDllbDeerDHHo6tkyRJ9ZKhWSunUSPYeedk4saUKXDeeTB+POy+O2ywAVx+OXzzTbarlCRJqhaGZlXduuvCxRfDtGlw//3JZil/+UtyvKAAXnkFMjzaUJIkKZMMzao+TZvCgQfCiy/CBx/ACSfAE08kFw1uvnmyDtqxdZIkqQ4yNCszNt4Y/v53mDkTbr8dmjSBk06C9u2Tn++9l+0Ka0RpaambUEiSVA8YmpVZLVrAsccm653HjoX994d//hN69EjOQP/rX46tkyRJtZ6hWTWjfGzdsGEwYwZcfTV89VWy5nnddZM10FOmZLtKSZKkZTI0q+aVj637+GN49lno12/h2Lrdd0/WQTu2TpIk1SKGZmVPo0aw007w8MMwdSqcfz68/XYy73n99eFvf3NsnSRJqhUMzaodOnaEiy5KxtY98EBy1vmvf02OH3qoY+skSVJWGZpVuzRtCgcckIyt+/DDZNLGk08mFw326AE33wz//W+2q5QkSQ2MoVm110YbwfXXJ2Pr7rgDmjWDk0+GDh3gxBPh3XezXaEkSWogDM2q/Vq0gGOOgXHj4I03kjPRw4YlG6Zsuy2UlDi2TpIkZZShWXVHCNCnTzLneeZMuOaa5ELBww5L1j6fey5MnpztKiVJUj1kaFbdtNZacMYZ8NFHydi67bZLZj+vv34ytu7xxx1bJ0mSqo2hWXVb+di6hx5KNkc5/3yYMAH23NOxdZIkqdoYmlV/lI+tmzoVRo5MQnP52LpDDoGXX3ZsnSRJWimGZq1QaWkppaWl2S6jcpo2hf33hxdeSJZvnHwyPPUUbL89bLYZ3HSTY+skSVKlGJpVv224IVx3HcyaBXfeCausAqecAu3bwwknwDvvZLtCSZJUBxia1TDk5MDRRydj6958Ew46CO6+G3r2pNcpp8B998Evv2S7SkmSVEsZmtXwbLkl3HVXMrbu2mtp+uOPcPjhsO66cM458MUX2a5QkiTVMoZmNVxrrQWnn86bd98Nzz2XrHm+5hrYYAPYbTd47DHH1kmSJMDQLCVj63bcER58MJm8ccEFyVrnvfaC9daDoiL4+utsV6laJD8/n/z8/GyXIUmqQYZmaVEdOsCFFyYzn0eOhK5d4bzzkqUbgwfDmDGOrZMkqQEyNEvLUj627vnnF46te+YZyM+HTTeFG2+EH3/MdpVSneYZe0l1iaFZSqd8bN3MmckFhDk5cOqpyVnp44+HiROzXaEkScowQ7NUUTk5cNRR8NZbye3gg+Gee6BXL9hmG7j3XsfWSZJUTxmapZWRl5dsljJrVnIW+j//gT/+Mdmy++yz4fPPs12hJEmqRobmLKmT21NraWuuCUOGJOuen38+WfN87bXJ2LpBg2DUKMfWSZJUDxiapeoQAgwcmEzcmDo1mcDx3nuw997QpUsytu6rr7JdpSRJWkmGZqm6degA//d/ydi6Bx9MLiQsH1t38MGOrZMkqQ4yNEuZ0rQp7Ldfstvgxx8nEzeee27h2LobbnBsnSRJdYShWaoJ3bola51nzEjG1rVoAX/6E7RvD4WFMGFCtiuUJEkrYGiWalL52Lo330zG1g0eDPfdB717w9ZbJyPsHFsnSVKtY2iWsqV8bN3MmXD99fDDD3DEEcnYurPOcmydJEm1iKFZyrY114TTToMPP4QXXoABA5LZz+Vj6x59FObPz3aVkiQ1aIZmqbYIAXbYAR54AKZNg4sugvffh332ScbWXXqpY+sauJKSEsaOHcuYMWPIzc2lpKQk2yVJUoNhaJZqo/bt4YILkrF1Dz0EG28M55+fjK076CAoLXVsXQNTUlJCYWEhc+fOBWDq1KkUFhYanCWphhiapdqsSRPYd1949ln45JNk4sbzzydLODbZxLF1DcjQoUMpKytb7FhZWRlDhw7NUkWS1LAYmqW6omtXuOaa5MLBf/4TVltt4di6446Dt9/OdoXKoGnTplXquCSpehmapbpm1VXhyCPhjTdg3Dg45BAoKYEttoC+feHuu+F//8t2lapmnTp1qtRxSVL1MjRLddkWW8Add8CsWcnYutmzk0DdsSP8+c/w2WdZLlDVpaioiJycnMWO5eTkUFRUlKWKJKlhMTRL9UGrVgvH1r34YjKF4+9/T5Z07LILPPKIY+vquIKCAoqLi2nevDkAnTt3pri4mIKCgixXJkkNQ5NsFyCpGoWQXCQ4YAB8+WVyFvq225KLCTt2TLbsPvZYaNcu25VqJRQUFHD77bcDUFpamt1iJKmB8UyzVF+1a5eMqZsyBR5+GLp3T8bYdeqUjK0bPdqxdZIkVZChWarvmjRJNkh55plkbN1ppyVj63bYIQnS//hHshZakiQtl6FZaki6doWrr07G1g0bBmuskYToDh2SZRuOrZMkaZkyFppDCOuGEEaHED4MIUwKIZyWOr5WCOG5EMKnqZ9rZqoGScux6qpwxBEwdiyMHw+HHgrDhyfTOLbayrF1kiQtIZNnmucDZ8YYNwb6AieHELoD5wIvxBi7Ai+kHkvKlt694fbbk7PPf/87/Pe/ydi6Dh3gzDPh00+zXaEkSVmXsdAcY/wyxvh26v5PwIdAB2Bv4O7Uy+4G9slUDZIqoVWrZIfBDz5ILhLcccdkvXO3brDzzsnFhI6tkyQ1UCHWwNXzIYRc4CVgU2BajLHVIs/9EGNcaolGCKEQKARo27btFiNGjMh4nXXRnDlzaNmyZbbLqNPs4fI1+89/aPfEE7R7/HFW+fZb5rZuzaw99uDL3Xfn19atF7yuofVwyJAhAFx//fXV9pkV7WEmvjtbstlHLZ89rDp7WHXZ7OGAAQPGxxjzljye8dAcQmgJjAGKYowPhRBmVyQ0LyovLy+OGzcuo3XWVaWlpeTn52e7jDrNHlbA/PnwxBNw883w7LMLJ3KceCIMGEDpmDENqoflv2t1zkqu6N9hJr47W7LZRy2fPaw6e1h12exhCGGZoTmj0zNCCE2BB4GSGONDqcNfhxDapZ5vB3yTyRokVYMmTWDvvZOxdZ9+CkOGJDsPDhwI3bvTYeRIx9ZJkuq1TE7PCMCdwIcxxmsXeWoUcETq/hHAo5mqQVIGbLABXHUVzJiRTNlo1YquN90E7dsnY+vGj892hZIkVbtMnmnuBxwO7BBCmJi67QZcDuwUQvgU2Cn1WFJds+qq8Mc/wuuvM664GA47LBlbl5cHffokc6AdWydJqicyOT3jlRhjiDH2iDH2TN2ejDH+J8Y4MMbYNfXz+0zVIKlmzOnaFYqLYdasZOLGnDlw1FHJ2Lozzkh2IpQkqQ5zR0BJ1WeNNeDUU2HSpGRs3U47wQ03wIYbJvcdWydJqqMMzZKqXwiQnw///jdMnw6XXAIffwz77Qe5uXDRRclZaUmS6ghDs6TM+sMf4Lzz4Isv4NFHYdNN4cILoVMnOOAAeOEFqIF58ZIkVYWhWVLNaNIE9toLnn46GVt3+ukLdx7ceGO4/nr44YdsVylJ0jIZmiXVvPKxdTNnwj33wJprJiG6Qwc45hhwMyNJUi1jaJaUPausAocfDq+/DhMmJPdHjIAtt0zG1v3zn1BWlu0qJUkyNEuqJXr2hNtuSy4QvOEG+PlnOPpox9ZJkmoFQ7Ok2mWNNeCUU+D996G0FHbZZeHYuh13hIcecmydJKnGrTA0hxC2DiHcFEJ4N4TwbQhhWgjhyRDCySGENWqqSEkNUAjQv3+yXGP6dLj00uRs8/77Q+fOyQSOmTOzXaVWUklJCWPHjmXMmDHk5uZSUlKS7ZIkaYWWG5pDCE8BxwLPAIOAdkB34DxgFeDREMJeNVGkpAbuD3+AoUNh8uRkbF2PHnDxxUl43n9/x9bVMSUlJRQWFjJ37lwApk6dSmFhocFZUq22ojPNh8cYj4kxjooxzooxzo8xzokxvh1jvCbGmA+8VkN1ShI0bpyMrXvqqWRs3RlnwJgxybKNjTaC665zbF0dMHToUMqWuMCzrKyMoUOHZqkiSUpvRaG5VQih35IHQwjbhRDWB4gxfpexyrRC+fn55OfnZ7sMKXvWXx+uvBJmzIB774W1105CdIcOyQWEb72V7Qq1HNOmTavUcUmqDVYUmq8HflrG8f+lnpOk7FtlFTjsMHjttYVj6+6/PxlZt+WWcNddjq2rZTp16lSp45JUG6woNOfGGN9d8mCMcRyQm7GKJGlllY+tmzkTbrwxCcvHHJOcfT79dPj442xXKKCoqIicnJzFjuXk5FBUVJSliiQpvRWF5lVW8Nyq1V2IJFWbNdaAk09OxtaNGQODBsFNNyXrnnfcER58EObNW6mPdupD1RUUFFBcXEzz5s0B6Ny5M8XFxRQUFGS5MklavhWF5rdCCMcteTCEcAwwPnMlSVI1CQG23x6GD0/G1hUVJRcQHnDASo2tc+pD9SkoKKBv377079+fKVOmGJgl1XorCs1DgKNCCKUhhGtStzEkY+hOq5HqJKm6tG0Lf/0rfPEFjBqVLOUoH1u3337w/PPw++8r/AinPkhSw7Xc0Bxj/DrGuA1wETAldbsoxrh1jPGrmilPkqpZ48aw557w5JPw2Wdw5pnw8suw007J8o1rr4Xvv1/mW536IEkNV9pttGOMo2OMN6RuL9ZEUZJUI9ZbD664Ilm6ce+90KZNEqI7dICjjlpqbJ1THySp4VrRjoAHhhAeCSE8HEI4uCaLkqQaVT627tVXYeJEOOIIeOCBZGxdXh7ceSeUlTn1QZIasBWdaT4H2A/YHzi7ZsqRpCzbfHO49VaYNSuZuPHLL3DssdChAwVvvUXxhRc69UGSGqAmK3juPuCe1P0HaqAWSao9Vl8dTjoJTjwRXnkFbr4Zbr6Zgnnz2KhVK0atvz4XTZwITZtmu1JJUg1YbmiOMV4fQmgBhBjjnBqsSZJqjxBgu+2S29dfw1130eGii7jogw+SyRvHHZfcOnbMdqWSpAxa0ZrmEGP8eUWBOYQQMlOWJNVCbdvCX/7CoVttxV823RR69YJLLoHcXNh3X3juubRj6yRJddOK1jSPDiGcGkJY7LLwEEKzEMIOIYS7gSMyW54k1T6/h8Dra68NTzwBn38Of/5zsoRj551hww1XOLZOklQ3rSg0DwJ+A4aHEGaFED4IIXwBfAocAlwXYxxWAzVKUu3VpQtcfjnMmAElJcnZ6PKxdUceCW++CTFmu0pJUhWtaHOTX2KMN8cY+wGdgYFA7xhj5xjjcTHGiTVVpCTVes2bw6GHJmec33knmfP84IOw1VbJ2Lo77oCff852lZKklZR2cxOAGOO8GOOXMcbZGa5Hkuq+Hj2SaRuzZiU/f/01uViwQwc47TT46KNsVyhJqqQKhWZJ0kpYbbVkZN277yZbde++O9xyC2y8MeywQ7KByrx52a5SklQBhmZJyrQQYNttkzXPM2bA3/4GkyfDQQclY+suuIDm336b7SolSStQodAcQugcQtgxdX/VEMJqmS1LkuqpddaBc8+Fzz6Dxx+H3r3h0kvpO3hwMrbu2WcdWydJtVDa0BxCOA4YCdyWOtQReCSDNUlS/de4cbJc4/HH4fPPmTZ4MLz6KuyySzK27ppr4D//yXaVkqSUipxpPhnoB/wXIMb4KbBOJouSpAalSxcmH3ccTJ+eLOH4wx+S2c8dOsARR8Abbzi2TpKyrCKheW6M8dfyByGEJoD/9Jak6lY+tu7ll5OLB48+Gh56CPr2hS22gDvuYJXffst2lZLUIFUkNI8JIfwVWDWEsBPwAPBYZsuSpAZus80Wjq275RaYPx+OO46Rr7/OqZ99Bh9+mO0KJalBqUhoPgf4FngPOB54Ejgvk0VJklJWWw1OOCHZMOWVV3h97bXZc9Ys6N4dBgyA++9P5kBLkjJqhaE5hNAIeC/GeHuM8cAY4wGp+y7PkKSaFAL060fRxhtzYN++ydbdU6bAwQcnY+vOPz9ZEy1JyogVhuYY4+/AOyGETjVUjyQpjR+bNYNzzknG1j3xRLJNd1ER5ObCPvvAM884tk6SqllFlme0AyaFEF4IIYwqv2W6MElSGo0bw267wWOPwRdfJEH6tddg0CDo1g2uvtqxdZJUTSoSmi8C9gAuBq5Z5CZJqi1yc+Gyy5IlGv/6F7RvD2edtXBs3dixjq1rYPLz88nPz892GVK9kTY0xxjHLOtWE8VJkiqpeXM45BB46SV47z045hh4+GHYeutk98Hbb4eff852lZJU51RkR8CfQgj/Td1+CSH8FkL4b00UJ0mqgk03hZtugpkzk7F1v/8OhYXJWehTT4UPPsh2hZJUZ1TkTPNqMcbVU7dVgP2BGzNfmiSpWpSPrZs4Mdmqe889obgYNtkE8vMdWydJFVCRNc2LiTE+AuxQ/aVIkjIqBNhmG7jvPpgxA664AqZNS8bWdeoE552XPJYkLaUiyzP2W+R2QAjhctxGW5LqtjZt4Oyzk7F1Tz4JffokFxJ26QJ77w1PP+3YOklaRJMKvGbPRe7PB6YAe2ekGknSCpWWllbvBzZqBLvumtymTk2WbdxxB4waBeuvD8cfD0cdBa1bV+/3SlIdU5HlGXfEGI9K3Y6LMRYBXTNdmCSphnXunGySMn06DB+ejKs7+2zo2BH++Ed4/XXH1klqsCoSmm+o4LHFhBDuCiF8E0J4f5FjF4YQZoYQJqZuu1WmWElSDWjWDAYPhjFjkrF1xx4LjzySrIfu3Ts5Gz1nTrarlKQatdzQHELYOoRwJtAmhHDGIrcLgcYV+OxhwKBlHL8uxtgzdXtypaqWJNWMTTeFG2+EWbPg1luTM83HH5+chXZsnaQGZEVnmpsBLUnWPa+2yO2/wAHpPjjG+BLwfTXUKEnKtpYtk7A8YUKyVfdeey0cW9e/P/z7346tk1SvLfdCwNSuf2NCCMNijFOr8TtPCSH8ERgHnBlj/GFZLwohFAKFAG3btq3+i1/quNmzZwMwZ84ce1NF9rDqGloPy//3V52/c53r4THH0HT//fnDU0/R/rHHWHXwYH5dc02+3G03Zu25J3Pbtk37EfYxs1a2v/aw6uxh1dXGHlZkekZZCOEqYBNglfKDMcaVmdV8C3AJyci6S4BrgKOX9cIYYzFQDJCXlxfz8/NX4uvqr1atWgHQsmVL7E3VlJaW2sMqamg9LP/fX3X+znW2h/vsk+w2+OyzNLv5ZjoPH07n4cNh993hpJNg552TCR3LYB8za2X7aw+rzh5WXW3sYUUuBCwBPgK6ABeRjJx7a2W+LMb4dYzxtxjj78DtQJ+V+RxJUi3SqBEMGpSMqZs8Gf7yF3jjjWSMXdeucOWV8N132a5SkqqkIqF57RjjncC8GOOYGOPRQN+V+bIQQrtFHu4LvL+810qS6qBOneDSS5OxdSNGwLrrwjnnJGPrDj88WQ/t2DpJdVBFlmfMS/38MoSwOzAL6JjuTSGE4UA+0DqEMAP4PyA/hNCTZHnGFOD4ypcsSar1mjVLtuc++GCYNCmZvHH33ckW3ptvDieeSOnjjycXGEpSHVCRM82XhhDWAM4E/gzcAZye7k0xxkNijO1ijE1jjB1jjHfGGA+PMW4WY+wRY9wrxvhlFeuXJNV2m2wCN9yQjK277TYIAU44Adq3h1NOSUK1JNVyKwzNIYTGQNcY448xxvdjjANijFvEGEfVUH2SpPqiZUsoLIS3306WaeyzT7Jl96abJmPrRoxwbJ2kWmuFoTnG+BuwVw3VIklqCEKArbeGe+6BGTOSCwVnzIBDDknWQA8dClOrc9KpJFVdRZZnvBZCuDGEsF0IoXf5LeOVSZLqv9at4ayz4NNP4emnoW9fuPxyWG892HNPeOop+P33bFcpSRUKzduQzGi+mGSu8jXA1ZksSitWUlLC2LFjGTNmDIMHD6akpCTbJUlS1TRqBLvsAo8+unBs3VtvwW67wQYbwBVXwLffZrtKSQ1Y2tCcWse85G1lNjZRNSgpKaGwsJC5c+cC8PXXX1NYWGhwllR/lI+tmzYt2Z67Uyc499xkbN1hhzm2TlJWpA3NIYS2IYQ7QwhPpR53DyEck/nStCxDhw6lrKxssWNlZWUMHTo0SxVJUoY0awYHHQSlpcmEjeOPh8ceg379oGfPZIzdTz9lu0pJDURFlmcMA54B2qcefwIMyVA9SmPatGmVOi5J9UL37vCPf8DMmVBcnCznOPFE6NABTj4Z3nevLEmZVZHQ3DrGeD/wO0CMcT7wW0ar0nJ16tSpUsclqV5p2RKOOy4ZW/f667DvvnDnnbDZZrD99qzzwguQWr4mSdWpIqH55xDC2iS7+BFC6Av8mNGqtFxFRUXk5OQsdiwnJ4eioqIsVSRJWRBCMmnj7ruTs89XXQWzZtH90kuTsXV//StMmZLtKiXVIxUJzWcAo4D1QwivAvcAp2a0Ki1XQUEBxcXFNG/eHIC2bdtSXFxMQUFBliuTpCxZe23485/hk09454orYJttkmkb5WPrnnwSfvM/kEqqmopMz3gb6E8yeu54YJMY47uZLkzLV1BQQN++fenfvz8jRowwMEsSQKNG/NCnDzzySHKWeehQGDcOdt/dsXWSqqwi0zNWAf4EXAJcBJycOiZJDVJpaSmlpaXZLkMrsu66cMklydi6+++H3NyFY+sKCuDVVx1bJ6lSKrI84x6SzU1uAG4EugP3ZrIoSZKqRdOmcOCBMHr0wrF1jz8O224Lm28Ot9zi2DpJFVKR0LxhjPGYGOPo1K0Q6JbpwiRJqlblY+tmzYLbb4cmTeCkk6B9++Tne+9lu0JJtVhFQvOE1MQMAEIIWwGvZq4kSZIyqEULOPZYGD8exo6F/feHu+6CHj1gu+3gX/9ybJ2kpVQkNG8FvBZCmBJCmAK8DvQPIbwXQvCCQElS3RQCbLUVDBuWjK27+mr46qtkzfO668Jf/uLYOkkLVCQ0DwK6kEzQ6J+6vxuwB7Bn5kqTJKmGrL02nHkmfPwxPPNMslX3lVcmY+t23x2eeMKxdVIDV5GRc1OB/wJrAGuX32KMU1PPSZJUPzRqBDvvDA8/nJxlPu+8ZPfBPfaA9deHv/0Nvvkm21VKyoKKjJy7BHgX+AdwTep2dYbrkiQpu9ZdFy6+eOHYuvXWS3Ya7NgRDj0UXnnFsXVSA1KR5RkHAevHGPNjjANStx0yXZgkSbVC+di6F1+EDz6AE09Mdhncbrvk4sGbb4b//jfbVUrKsIqE5veBVhmuQ5Kk2m/jjeHvf08uHLzjDmjWDE4+GTp0SML0u7Xj+viSkhLGjh3LmDFjyM3NpaSkJNslSXVeRULz30jGzj0TQhhVfst0YZIk1VotWsAxxyTbdL/xRjK2btiwZMOUbbeFkpKsja0rKSmhsLCQuanvnzp1KoWFhQZnqYoqEprvBq4ALmfhmuZrMlmUJEl1QgjQp8/CsXXXXANffw2HHZasfT73XJg8uUZLGjp0KGVlZYsdKysrY+jQoTVah1TfVCQ0fxdj/EdqN8Ax5beMVyZJUl2y1lpwxhnJ2Lpnn03WPF91VTJ1Y/fdk+27a2Bs3bRp0yp1XFLFVCQ0jw8h/C2EsHUIoXf5LeOVSZJUFzVqBDvtBA89BFOnwvnnw4QJsOeeNTK2rlOnTpU6LqliKhKaewF9gctw5JwkSRXXsSNcdFESnh94IAnN5WPrDjkEXn652sfWFRUVkZOTs9ixnJwcioqKqvV7pIamIpubDFjGzZFzkiRVVNOmcMAB8MIL8OGHcNJJ8NRTsP32sNlmcNNN1Ta2rqCggOLiYpo3bw5A586dKS4upqCgoFo+X2qoKrK5SdsQwp0hhKdSj7uHEI7JfGmSJNVDG20E11+/cGzdKqvAKadA+/ZwwgnwzjtV/oqCggL69u1L//79mTJlioFZqgYVWZ4xDHgGaJ96/AkwJEP1SJLUMCw6tu7NN5MNVO6+G3r2hH794L774Jdfsl2lpJTlhuYQQpPU3dYxxvuB3wFijPOBzF/+K0lSQ7HllvDPfyZnn6+9Fr79Fg4/PNnK+5xz4Isvsl2h1OCt6Ezzm6mfP4cQ1gYiQAihL/BjpguTJKnBWWstOP10+OgjeO65ZGzdNdfABhvAbrvBY4/VyNg6SUtbUWgOqZ9nAKOA9UMIrwL3AKdmujBJkhqsRo1gxx0Xjq274AKYOBH22gvWWw8uuyzZREVSjVlRaG4TQjgDyAceBq4EngJuB3bMfGmSJIkOHeDCC5PwPHIkdO0KQ4cmSzcGD4aXXqr2sXWSlrai0NwYaAmsBrQAmqSO5aSOSZKkmtK0Key/Pzz/fLJ84+ST4ZlnoH//ah9bJ2lpTVbw3JcxxotrrBJJklQxG24I110HRUUwYgTccksytu6cc6CgAE48MdsVSvVORdY0S5Kk2ignB44+Gt56Kxlbd9BBcM890KsXN06YwE5ff+3YOqmarCg0D6yxKiRJUtVsuSXcddeCsXWrz5vH0I8+SrbsPvts+PzzbFco1WnLDc0xxu9rshBJklQNUmPr/rjllpzRo0ey5vnaa5OxdbvuCqNGObZOWgkV2RFQkiTVNSHw9pprwoMPJpM3LrwQ3n0X9t47GVtXVARffZXtKqU6w9AsSVJ916ED/N//wZQpSYju1g3OO2/h2LoxYxxbJ6VhaJYkqaFo2hT22y/ZbfDjj+HUU5Oxdfn5sOmmcOONNJ4zJ9tVSrWSoVmSpIaoW7dkrfPMmckFhC1awKmnss2BB0JhYbIDoaQFDM2SJDVkOTlw1FHJyLq33uKbAQPgvvugVy/Yemu4917H1kkYmiVJUrm8PD4+++zk7PN118H338Mf/5iMrTvrLMfWqUEzNEuSpMWtuSYMGZJs1/3CC8ma5+uuS8bWDRoEjz4K8+dnu0qpRhmaJUnSsoUAO+wAI0cuHFv33nuwzz7J2LpLL3VsnRoMQ7MkSUqvfGzd1Knw0EOw4YZw/vnJ2LqDD4bSUsfWqV7LWGgOIdwVQvgmhPD+IsfWCiE8F0L4NPVzzUx9vyRJyoAmTWDffReOrfvTn5L7AwbAJpvADTfAjz9mu0qp2mXyTPMwYNASx84FXogxdgVeSD2WJEl1UbducM01yYWD//wnrLZaEqLbt0/G1k2YkO0KpWqTsdAcY3wJ+H6Jw3sDd6fu3w3sk6nvlyRJNWTVVeHII+GNN2DcODjkkGRsXe/e0Lcv3HOPY+tU54WYwfVHIYRc4PEY46apx7NjjK0Wef6HGOMyl2iEEAqBQoC2bdtuMWLEiIzVWRcNGTIEgEsvvZSWLVtmt5g6bs6cOfawiuxh1dnD6mEfFyr/98T1119fqfdVVw+b/PQTbZ95hg6jRpEzfTrzVl+drwYNYtZee/G/Dh2q/Pm1mX+HVZfNHg4YMGB8jDFvyeNNslFMRcQYi4FigLy8vJifn5/dgmqZVq1aAdCyZUvsTdWUlpbawyqyh1VnD6uHfVyo/N8Tle1HtfZwzz2TNc6jR9P0lltY96GHWPf++2HnneGkk2D33ZM10vWMf4dVVxt7WNPTM74OIbQDSP38poa/X5Ik1aTysXUPPJBM3rjoIpg0KRlb16ULXHIJfPlltquU0qrp0DwKOCJ1/wjg0Rr+fkmSlC3t28MFF8CUKfDww7DxxsnjTp3goINg9GjH1qnWyuTIueHA68CGIYQZIYRjgMuBnUIInwI7pR5LkqSGpEmT5Ezzs8/CJ58kEzeefz45I73JJvCPf8Ds2dmuUlpMJqdnHBJjbBdjbBpj7BhjvDPG+J8Y48AYY9fUzyWna0iSpIaka9eFY+uGDUvG1p12WrKZynHHwdtvZ7tCCXBHwDqrtLSU0tLSbJchSVL1WHVVOOKIZGzd+PFw6KFQUgJbbJGMrbv7bvjf/7JdpRowQ7MkSapdeveG22+HWbPg739Pdhg88kjo2BH+/Gf47LNsV6gGyNAsSZJqp1atkvXOH3wAL74IAwcmIbpr12Rs3SOPwPz52a5SDYShWZIk1W4hwIABcP/9MG0aXHwxfPgh7LtvMrbu4osdW6eMMzRLkqS6o107OP98mDw5GVvXvTv83/8lY+sOPNCxdcoYQ7MkSap7ysfWPfNMMrbutNOSJRw77JAE6b//3bF1qlaGZkmSVLd17QpXXw0zZiRTNtZYA4YMSTZTOfbYZBqHVEWGZkmSVD+suir88Y8wdmwSlAsKYPhwyMuDrbZK5kA7tk4rydAsSVI91ODn+ZePrZs5M9lh8L//haOOSjZNOfNM+PTTbFeoOsbQLEmS6q9WreDUU5OxdaNHw047JSG6W7dkbN3DDzu2ThViaJYkSfVfCJCfD//+dzK27pJL4KOPYL/9IDc3GVs3a1a2q1QtZmiWJEkNS7t2cN558MUXyQYpm266cGzdAQckUziqMLZuyJAh5OfnV1u5qh0MzZIkqWFq0gT23huefjpZ43z66ckSjoEDYeON4frr4Ycfsl2laglDsyRJ0gYbwFVXLRxbt+aaSYju0AGOOcaxdTI0S5IkLVA+tu711+Htt+Gww2DEiGRsXZ8+8M9/QllZtqtUFhiaJUmSlqVXLyguTi4QvOEGmDMHjj4aOnaEM85IdiJUg2FoliRJWpE11oBTToFJk6C0NBlbd8MNsOGGyf2HHnJsXQNgaJYkSaqIEKB//2Rs3fTpcOml8PHHsP/+0LkzXHSRY+vqMUOzJElSZf3hDzB0KEyeDI8+Cj16JKG5UyeumjKF3j/8UKWxdap9DM2SJEkrq3Fj2GsveOqpZGzdGWeQN2cO1777Lmy0kWPr6hFDsyRJUnVYf3248kp26d6doo02grXXXji27uijYdy4bFeoKjA0S5IkVaNfGzXiubZt4bXXYMIEOPxwuP9+2HLL5HbXXY6tq4MMzZIkSZnSsyfcdhvMnAk33piE5WOOSc4+n356ciGh6gRDsyRJUqatsQacfDK8/z6MGQODBsFNNyXrnnfcER58EObNy3aVWgFDsyRJUk0JAbbfHoYPT8bWFRUlFxAecADk5sKFFyZnpVXrGJolSZKyoW1b+Otf4YsvYNQo2HxzuPjiZObz/vvD88/D779nu0qlGJolSZKyqXFj2HNPePJJ+OwzOPNMeOmlZLfBjTaC665zbF0tYGiWJEmqLdZbD664Ilm6ce+90KYNnHEGtG+fjK17661sV9hgGZolSZJqm1VWgcMOg1dfhYkT4YgjkrF1ffpAXp5j67LA0CxJklSbbb453HorzJqVTNz45ZeFY+uGDHFsXQ0xNEuSJNUFq68OJ50E772XrHkeNAhuvjlZ9zxwIIwc6di6DDI0S5Ik1SUhwHbbLRxbd9llyQWEBx6YTN74v/+DGTOyXWW9Y2iWJEmqq9q2hb/8JRlb99hj0KsXXHJJMvN5v/3gueccW1dNDM2SJEl1XePGsMce8MQT8Pnn8Oc/w8svw847J8s3rr0Wvv8+21XWaYZmSZKk+qRLF7j88mSJxn33wTrrJLOfO3SAo46CN9+EGLNdZZ1jaJYkSaqPmjeHggJ45RV45x048sjkYsGttkrG1t15p2PrKsHQLEmSVN/16AG33AIzZyZj6379FY49Ntk0ZcgQ+OijbFdY6xmaJUmSGorysXXvvpused5tt2Rs3cYbww47OLZuBQzNkiRJDU0IsO228K9/JWuf//a3ZAJH+di6Cy5wbN0SDM2SJEkN2TrrwLnnJlM3Hn88GVt36aXJ2Lp994Vnn3VsHYZmSZIkQTK2bvfdF46tO+us5CLCXXaBDTeEa66B//wn21VmjaFZkiRJi+vSJVmyMWMGlJTAH/6QzH7u2DGZwvHGGw1ubJ2hWZIkScvWvDkcemhy0eC77yZznh98EPr2TcbW3XEH/PxztqusEYZmSZIkpbfZZsmkjVmzkp/z5sFxxyWbppx2Gnz4YbYrzChDsyRJkiputdXgxBOTDVNeeSVZB33rrdC9OwwYAA88UC/H1hmaJUmSVHkhQL9+yZrn6dOTNdBTpsBBB0GnTsnYuunTs11ltclKaA4hTAkhvBdCmBhCGJeNGiRJklRNysfWffZZMn1jiy0Wjq3bZ596MbYum2eaB8QYe8YY87JYgyRJUrUpKSnhgw8+YMyYMeTm5lJSUpLtkmpW48bJLoOPP55slnLOOfDaa8nYum7d4Oqr6+zYOpdnSJIkVYOSkhIKCwuZl1rPO3XqVAoLCxtecC6XmwuXXZYs0fjXv6B9+2T2c4cOcMQRMHZsnRpbl63QHIFnQwjjQwiFWapBkiSp2gwdOpSysrLFjpWVlTF06NAsVVRLNG8OhxwCL72UjK075hh46CHYeutkGcftt9eJsXUhZiHhhxDaxxhnhRDWAZ4DTo0xvrTEawqBQoC2bdtuMWLEiBqvsy6YM2cOLVu2zHYZdZo9rDp7WHX2sHrYx6qzhytvhx12YFm5KoTAiy++mIWKaq/GZWW0ff552j/6KC2/+IL5LVrw1c47M2vvvSnr3Dmrf4cDBgwYv6zlw1kJzYsVEMKFwJwY49XLe01eXl4cN87rBZeltLSU/Pz8bJdRp9nDqrOHVWcPq4d9rDp7uPJyc3OZOnXqUsc7d+7MlClTar6guiDGZM3zLbcko+p+/RX692dS//5scsEFyRrpGhZCWGZorvHlGSGEFiGE1crvAzsD79d0HZIkSdWpqKiInJycxY7l5ORQVFSUpYrqgPKxdffdl2zZffnlMG0aXf75T2hUuy69y0Y1bYFXQgjvAG8CT8QYn85CHQ3Gsv7zxq233so999yThWqWXY8kSXVdQUEBxcXFNG3aFEjOMBcXF1NQUJDlyuqINm2SaRuffcY7V12VBOpapElNf2GM8Qtg85r+Xi3uhBNOyHYJkiTVOwUFBVx11VW0atWK0tLSbJdTNzVqxNy2bbNdxVJq13lv1ZgLL7yQq69OlpHn5+dzzjnn0KdPH7p168bLL78MwG+//cZZZ53FlltuSY8ePbjtttuW+pxzzjmHm2++ebHPveaaa5gzZw4DBw6kd+/ebLbZZjz66KNLvbe0tJQ99thjweNTTjmFYcOGATB+/Hj69+/PFltswS677MKXX34JwD/+8Q+6d+9Ojx49GDx4cLX1Q5IkaUUMzQJg/vz5vPnmm1x//fVcdNFFANx5552sscYavPXWW7z11lvcfvvtTJ48ebH3DR48mH//+98LHt9///0ceOCBrLLKKjz88MO8/fbbjB49mjPPPHOZVxQvy7x58zj11FMZOXIk48eP5+ijj14wrufyyy9nwoQJvPvuu9x6660AjBs3jmOPPbY62iBJkqpRfn5+vbmwtMaXZ6h22m+//QDYYostFlzh++yzz/Luu+8ycuRIAH788Uc+/fRTunTpsuB9vXr14ptvvmHWrFl8++23rLnmmnTq1Il58+bx17/+lZdeeolGjRoxc+ZMvv76a/7whz+kreXjjz/m/fffZ6eddgKSM97t2rUDoEePHhQUFLDPPvuwzz77AJCXl8cdd9xRXa2QJElaiqFZADRv3hyAxo0bM3/+fABijNxwww3ssssuK3zvAQccwMiRI/nqq68WLJkoKSnh22+/Zfz48TRt2pTc3Fx++eWXxd7XpEkTfl9kH/ry52OMbLLJJrz++utLfdcTTzzBSy+9xKhRo7jkkkuYNGkSTZr4ZyxJkjLL5Rlarl122YVbbrllwXagn3zyCT8vY8eewYMHM2LECEaOHMkBBxwAJGel11lnHZo2bcro0aOXO7fygw8+YO7cufz444+88MILAGy44YZ8++23C0LzvHnzmDRpEr///jvTp09nwIABXHnllcyePZs5c+Zk6teXJElawFN0DUBZWRkdO3Zc8PiMM86o0PuOPfZYpkyZQu/evYkx0qZNGx555JGlXrfJJpvw008/0aFDhwXLKAoKCthzzz3Jy8ujZ8+ebLTRRku9b9111+Wggw6iR48edO3alV69egHQrFkzRo4cyZ/+9Cd+/PFH5s+fz5AhQ+jWrRuHHXYYP/74IzFGTj/9dFq1asW4ceO49dZbXaIhSZIyxtDcACy6BGJZFh2J07p16wVrmhs1asRll13GZZddlvY73nvvvcUet27depnLK4DFzg5feeWVXHnllUu9pmfPnrz00ktLHX/llVeWOuaaZkmSlGkuz5AkSZLSMDRLkiRJaRiaG4CioiI22WQTevToQc+ePXnjjTeAZHbiuHHjslzd8v3tb39jgw02YMMNN+SZZ55Z5mvOP//8Bb/XzjvvzKxZsyr1fkmSpIpwTXM99/rrr/P444/z9ttv07x5c7777jt+/fXXbJeV1gcffMCIESOYNGkSs2bNYscdd+STTz6hcePGi73urLPO4pJLLgGS3QIvvvhibr311gq/X5IkqSI801zPffnll7Ru3XrBHObWrVvTvn37pV43fPhwNttsMzbddFPOOeecBcdbtmzJmWeeSe/evRk4cCDffvstAJ9//jmDBg1iiy22YLvttuOjjz6q1rofffRRBg8eTPPmzenSpQsbbLABb7755lKvW3311Rfc//nnnwkhVOr9kiRJFWForud23nlnpk+fTrdu3TjppJMYM2bMUq+ZNWsW55xzDi+++CITJ07krbfeWjBa7ueff6Z37968/fbb9O/ff8EW24WFhdxwww2MHz+eq6++mpNOOmmpzx09ejQ9e/Zc6rbNNtukrXvmzJmsu+66Cx537NiRmTNnLvO1Q4cOZd1116WkpISLL7640u+XJElKx9Bcz7Vs2ZLx48dTXFxMmzZtOPjggxk2bNhir3nrrbfIz8+nTZs2NGnShIKCggXj3ho1asTBBx8MwGGHHcYrr7zCnDlzeO211zjwwAPp2bMnxx9/PF9++eVS3z1gwAAmTpy41O21115LW3eMcalj5WeRl1RUVMT06dMpKCjgxhtvrPT7JUmS0nFNcwPQuHFj8vPzyc/PZ7PNNuPuu+/myCOPXPD8sgLm8oQQ+P3332nVqhUTJ05c4WtHjx7N6aefvtTxnJycpYLzww8/vOAs9h133EHHjh2ZPn36gudnzJixzGUlizr00EPZfffdueiii1bq/ZIkScvjmeZ67uOPP+bTTz9d8HjixIl07tx5sddstdVWjBkzhu+++47ffvuN4cOH079/fyDZGGXkyJEA/Otf/2Lbbbdl9dVXp0uXLjzwwANAErrfeeedpb67Mmea99133wXP5+XlsddeezFixAjmzp3L5MmT+fTTT+nTp89S71v0dxs1atSCnQcr+n5JkqSK8ExzPTdnzhxOPfVUZs+eTZMmTdhggw0oLi5e7DXt2rXjb3/7GwMGDCDGyG677cbee+8NQIsWLZg0aRJbbLEFa6yxBv/+978BKCkp4cQTT+TSSy9l3rx5DB48mM0337za6t5kk0046KCD6N69O02aNOGmm25aMPni2GOP5YQTTiAvL49zzz2Xjz/+mEaNGtG5c2duvfXWtO+XJEmqrFCZ/zSfLXl5ebE2zxPOptLSUvLz8zP2+S1btlxs2+v6KNM9bAjsYdXZw+phH6vOHlZdz549adWqFaWlpdkuJevK/5Yq24ts/h2GEMbHGPOWPO7yDEmSJCkNQ7NWqL6fZZYkSaoIQ7MkSZKUhqG5AWjZsmXa11x//fWUlZXVQDVLmz17NjfffHO1fNa1115L9+7d6dGjBwMHDmTq1KkLnmvcuPGCDVb22muvZb7/9NNPX/Cabt260apVqwXPDRo0iFatWrHHHntUS62SJKnuMDQLWLnQ/Ntvv1XLd1dnaO7Vqxfjxo3j3Xff5YADDuDss89e8Nyqq666YKzdqFGjlvn+6667bsFrTj31VPbbb78Fz5111lnce++91VKnJEmqWwzNDUj5lagHHHAAG220EQUFBcQY+cc//sGsWbMYMGAAAwYMAODZZ59l6623pnfv3hx44IEL1jbn5uZy8cUXs+222/LAAw/w9NNP07t3bzbffHMGDhwIJFtvH3300Wy55Zb06tWLRx99FIBhw4ax9957M2jQIDbccMMFm5mce+65fP755/Ts2ZOzzjqrSr/jgAEDyMnJAaBv377MmDFjpT9r+PDhHHLIIQseDxw4kNVWW61K9UmSpLrJOc0NzIQJE5g0aRLt27enX79+vPrqq/zpT3/i2muvZfTo0bRu3ZrvvvuOSy+9lOeff54WLVpwxRVXcO2113LBBRcAsMoqq/DKK6/w7bff0rt3b1566SW6dOnC999/DyTbWu+www7cddddzJ49mz59+rDjjjsC8Oabb/L++++Tk5PDlltuye67787ll1/O+++/v9wdBrfbbjt++umnpY5fffXVCz53We6880523XXXBY9/+eUX8vLyaNKkCeeeey777LPPct87depUJk+ezA477JCupZIkqQEwNDcwffr0oWPHjkAyR3LKlClsu+22i71m7NixfPDBB/Tr1w+AX3/9la233nrB8wcffPCC122//fZ06dIFgLXWWgtIzlKPGjWKq6++GkjC6rRp0wDYaaedWHvttQHYb7/9eOWVV1YYXgFefvnlSv+e9913H+PGjWPMmDELjk2bNo327dvzxRdfsMMOO7DZZpux/vrrL/P9I0aM4IADDnBDFEmSBBiaG5zmzZsvuN+4cWPmz5+/1GtijOy0004MHz58mZ/RokWLBa8LISzz/Q8++CAbbrjhYsffeOONpV6/rPcvqbJnmp9//nmKiooYM2bMYr9v+/btAVhvvfXIz89nwoQJKwzNN910U9raJElSw+CaZgGw2mqrLQimffv25dVXX+Wzzz4DoKysjE8++WSp92y99daMGTOGyZMnAyxYnrHLLrtwww03UL7b5IQJExa857nnnuP777/nf//7H4888gj9+vVb7LuX5eWXX15wcd6it2UF5gkTJnD88cczatQo1llnnQXHf/jhB+bOnQvAd999x6uvvkr37t2X+X0ff/wxP/zww2Jn1yVJUsNmaBYAhYWF7LrrrgwYMIA2bdowbNgwDjnkEHr06EHfvn356KOPlnpPmzZtKC4uZr/99mPzzTdfsGzj/PPPZ968efTo0YNNN92U888/f8F7tt12Ww4//HB69uzJ/vvvT15eHmuvvTb9+vVj0003rfKFgGeddRZz5szhwAMPXGy03IcffkheXh6bb745AwYM4Nxzz10Qmu+6667FpmkMHz6cwYMHL3UWfLvttuPAAw/khRdeoGPHjjzzzDNVqlWSJNUdofxsYG2Wl5cXx40bl+0yaqVs7s1eWcOGDWPcuHHceOON2S5lMXWph7WVPaw6e1g97GPV2cOq69mzJ61ataK0tDTbpWRd+d9SZXuRzb/DEML4GGPeksc90yxJkiSl4YWAqjFHHnkkRx55ZLbLkCRJqjTPNEuSJElpGJolSZKkNAzNkiRJ1ej666/3IsB6yNAsSZIkpWFoliRJktIwNEuSJKnalZSUMHbsWMaMGUNubi4lJSXZLqlKDM2SJEmqViUlJRQWFjJ37lwApk6dSmFhYZ0OzoZmSZIkVauhQ4dSVla22LGysjKGDh2apYqqztAsSZKkajVt2rRKHa8LDM2SJEmqVp06darU8brA0CxJkqRqVVRURE5OzmLHcnJyKCoqylJFVWdoliRJUrUqKCiguLiY5s2bA9C5c2eKi4spKCjIcmUrr0m2C5AkSVL9U1BQwO233w5QL3ZI9EyzJEmSlIahWZIkSUojK6E5hDAohPBxCOGzEMK52ahBkiRJqqgaD80hhMbATcCuQHfgkBBC95quQ5IkSaqobJxp7gN8FmP8Isb4KzAC2DsLdUiSJEkVko3Q3AGYvsjjGaljkiRJUq2UjZFzYRnH4lIvCqEQKARo27ZtvRhVkglz5syxN1VkD6vOHladPawe9rHq7GHV2cOFZs+eDVR+5Fxt7GE2QvMMYN1FHncEZi35ohhjMVAMkJeXF/Pz82ukuLqmtLQUe1M19rDq7GHV2cPqYR+rzh5WnT1cqFWrVgCV7kdt7GE2lme8BXQNIXQJITQDBgOjslCHJEmSVCE1fqY5xjg/hHAK8AzQGLgrxjippuuQJEmSKior22jHGJ8EnszGd0uSJEmV5Y6AkiRJUhqGZkmSJCkNQ7MkSZKUhqFZkiRJSsPQLEmSJKVhaJYkSZLSMDRLkiRJaRiaJUmSpDQMzZIkSVIahmZJkiQpDUOzJEmSlIahWZIkSUrD0CxJkiSlYWiWJEmS0jA0S5IkSWkYmiVJkqQ0DM2SJElSGoZmSZIkKQ1DsyRJkpSGoVmSJElKw9AsSZIkpWFoliRJktIwNEuSJElpNMl2AZIkSaqfSktLs11CtfFMsyRJkpSGoVmSJElKw9AsSZIkpWFoliRJktIwNEuSJElpGJolSZKkNAzNkiRJUhqGZkmSJCkNQ7MkSZKUhqFZkiRJSsPQLEmSJKVhaJYkSZLSMDRLkiRJaRiaJUmSpDQMzZIkSVIahmZJkiQpDUOzJEmSlIahWZIkSUrD0CxJkiSlEWKM2a4hrRDCt8DUbNdRS7UGvst2EXWcPaw6e1h19rB62Meqs4dVZw+rLps97BxjbLPkwToRmrV8IYRxMca8bNdRl9nDqrOHVWcPq4d9rDp7WHX2sOpqYw9dniFJkiSlYWiWJEmS0jA0133F2S6gHrCHVWcPq84eVg/7WHX2sOrsYdXVuh66plmSJElKwzPNkiRJUhqG5joqhDAohPBxCOGzEMK52a6nrggh3BVC+CaE8P4ix9YKITwXQvg09XPNbNZY24UQ1g0hjA4hfBhCmBRCOC113D5WUAhhlRDCmyGEd1I9vCh13B5WUgihcQhhQgjh8dRje1gJIYQpIYT3QggTQwjjUsfsYSWEEFqFEEaGED5K/XNxa3tYOSGEDVN/g+W3/4YQhtS2Phqa66AQQmPgJmBXoDtwSAihe3arqjOGAYOWOHYu8EKMsSvwQuqxlm8+cGaMcWOgL3By6u/PPlbcXGCHGOPmQE9gUAihL/ZwZZwGfLjIY3tYeQNijD0XGe9lDyvn78DTMcaNgM1J/h7tYSXEGD9O/Q32BLYAyoCHqWV9NDTXTX2Az2KMX8QYfwVGAHtnuaY6Icb4EvD9Eof3Bu5O3b8b2Kcma6prYoxfxhjfTt3/ieRfEB2wjxUWE3NSD5umbhF7WCkhhI7A7sAdixy2h1VnDysohLA6sD1wJ0CM8dcY42zsYVUMBD6PMU6llvXR0Fw3dQCmL/J4RuqYVk7bGOOXkARCYJ0s11NnhBBygV7AG9jHSkktK5gIfAM8F2O0h5V3PXA28Psix+xh5UTg2RDC+BBCYeqYPay49YBvgX+mlgndEUJogT2sisHA8NT9WtVHQ3PdFJZxzDEoqlEhhJbAg8CQGON/s11PXRNj/C31nyI7An1CCJtmuaQ6JYSwB/BNjHF8tmup4/rFGHuTLPc7OYSwfbYLqmOaAL2BW2KMvYCfcSnGSgshNAP2Ah7Idi3LYmium2YA6y7yuCMwK0u11AdfhxDaAaR+fpPlemq9EEJTksBcEmN8KHXYPq6E1H/KLSVZa28PK64fsFcIYQrJErUdQgj3YQ8rJcY4K/XzG5I1pH2wh5UxA5iR+i9FACNJQrQ9XDm7Am/HGL9OPa5VfTQ0101vAV1DCF1S/1/ZYGBUlmuqy0YBR6TuHwE8msVaar0QQiBZv/dhjPHaRZ6yjxUUQmgTQmiVur8qsCPwEfawwmKMf4kxdowx5pL8M/DFGONh2MMKCyG0CCGsVn4f2Bl4H3tYYTHGr4DpIYQNU4cGAh9gD1fWISxcmgG1rI9ublJHhRB2I1nP1xi4K8ZYlN2K6oYQwnAgH2gNfA38H/AIcD/QCZgGHBhjXPJiQaWEELYFXgbeY+Fa0r+SrGu2jxUQQuhBclFLY5KTF/fHGC8OIayNPay0EEI+8OcY4x72sOJCCOuRnF2GZJnBv2KMRfawckIIPUkuRm0GfAEcRep/19jDCgsh5JBcr7VejPHH1LFa9bdoaJYkSZLScHmGJEmSlIahWZIkSUrD0CxJkiSlYWiWJEmS0jA0S5IkSWkYmiVJkqQ0DM2SlCEhhDmVeG1+CGGbRR6fEEL4Y+r+kSGE9ivx/VNCCK0r+Z6Rqfm95TWNCyFcucjzpSGEcYs8zgshlKbubxZCGFbZOiWpLjA0S1LtkA8sCM0xxltjjPekHh4JVDo0V1YIYROgcYzxi9ShE4HtgMYhhI0Week6IYRdl3x/jPE9oGMIoVOma5WkmmZolqQaFELYM4TwRghhQgjh+RBC2xBCLnACcHoIYWIIYbsQwoUhhD+HEA4A8oCS1HOrLnoGeYkzvWuHEJ5NffZtQFjkew8LIbyZ+ozbQgiNl1FeAYtvU9sIiCQ7P4ZFjl8FnLecX/Exkm2tJaleMTRLUs16BegbY+wFjADOjjFOAW4Frosx9owxvlz+4hjjSGAcUJB67n8r+Oz/A15JffYokq1nCSFsDBwM9Isx9gR+IwnIS+oHjF/k8R3Aa0CjGOOHixx/HZgbQhiwjM8YR3J2WpLqlSbZLkCSGpiOwL9DCO2AZsDkavzs7YH9AGKMT4QQfkgdHwhsAbwVQgBYFfhmGe9vB3xb/iDG+AzwzHK+61KSs83nLHH8G2pgKYkk1TTPNEtSzboBuDHGuBlwPLDKSnzGfBb+83vJ98dlvD4Ad6fOVPeMMW4YY7xwGa/7X0XriTG+mHpt3yWeWiX1OZJUrxiaJalmrQHMTN0/YpHjPwGrLec9Sz43heTMMcD+ixx/idSyi9SFemumjr8AHBBCWCf13FohhM7L+J4PgQ0q9FskioCzlzjWDXi/Ep8hSXWCoVmSMicnhDBjkdsZwIXAAyGEl4HvFnntY8C+5RcCLvE5w4Bbyy8EBC4C/p76jN8Wed1FwPYhhLeBnYFpADHGD0iWUjwbQngXeI5kKcaSniCZ4lEhMcYnWWQ5R8qA1OdIUr0SYlzWf8mTJDU0qUA+muSCwd/SvX4Z728OjAG2jTHOr+76JCmbDM2SpAVCCLsAH8YYp63Ee7sCHWKMpdVemCRlmaFZkiRJSsM1zZIkSVIahmZJkiQpDUOzJEmSlIahWZIkSUrD0CxJkiSl8f/fnUj7L4dJJgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Plot results\n", "fig, ax = plt.subplots(1, 1, figsize=(12,8))\n", "\n", "# Use ax.errorbar rather than ax.plot to show temperature deviations\n", "ax.errorbar(data[\"Latitude\"], data[\"Average temperature\"], fmt='ko', yerr=data[\"Deviation\"])\n", "\n", "# Calculate regression line points\n", "x1 = -5.0 # Minimum latitude\n", "x2 = 70.0 # Maximum latitude\n", "y1 = a + b * x1\n", "y2 = a + b * x2\n", "\n", "# Plot regression line in red\n", "ax.plot([x1, x2], [y1, y2], 'r-')\n", "\n", "# Add axis labels\n", "ax.set_xlabel(\"Latitude (°N)\")\n", "ax.set_ylabel(\"Temperature (°C)\")\n", "\n", "# Add text with regression line info\n", "ax.text(1, 2, f\"Line values:\\nSlope = {b:.2f}\\nIntercept = {a:.2f}\")\n", "\n", "# Enable plot grid\n", "ax.grid()" ] }, { "cell_type": "markdown", "id": "33869df2-76fa-4716-9f92-c97d5514c4ea", "metadata": {}, "source": [ "### Weighted least squares regressions\n", "\n", "Let's now imagine a case where the values on the y-axis have some nominal uncertainties $\\sigma$. We can modify our regression to include the uncertainties as weights for the point that are inversely proportional to the uncertainty:\n", "\n", "$$\n", "\\large\n", "w_{i} = 1 / \\sigma_{i}^{2}\n", "$$\n", "\n", "where $i$ is the ith weight for the ith uncertainty.\n", "\n", "We can thus modify our slope and intercept calculations above as follows:\n", "\n", "$$\n", "\\large\n", "A_{\\mathrm{w}} = \\frac{\\sum{wx^{2}} \\sum{wy} - \\sum{wx} \\sum{wxy}}{\\Delta_{\\mathrm{w}}},\n", "$$\n", "\n", "$$\n", "\\large\n", "B_{\\mathrm{w}} = \\frac{\\sum{w} \\sum{wxy} - \\sum{wx} \\sum {wy}}{\\Delta_{\\mathrm{w}}},\n", "$$\n", "\n", "and \n", "\n", "$$\n", "\\large\n", "\\Delta_{\\mathrm{w}} = \\sum{w} \\sum{wx^{2}} - \\left( \\sum{wx} \\right)^{2}.\n", "$$\n", "\n", "Let's consider another example below." ] }, { "cell_type": "code", "execution_count": 9, "id": "a4a6a807-969c-431d-94ce-a37bf6a3551c", "metadata": {}, "outputs": [], "source": [ "def calculate_deltaw(x_values, w_values):\n", " \"\"\"Calculates divisor term for weighted slope and intercept equations\"\"\"\n", " x_squared_sum = (w_values * x_values * x_values).sum()\n", " x_sum_squared = (w_values * x_values).sum()**2\n", " delta = w_values.sum() * x_squared_sum - x_sum_squared\n", " return delta\n", "\n", "def calculate_aw(x_values, y_values, w_values):\n", " \"\"\"Calcuates y-intercept for weighted regression line\"\"\"\n", " left_side = (w_values * x_values * x_values).sum() * (w_values * y_values).sum()\n", " right_side = (w_values * x_values).sum() * (w_values * x_values * y_values).sum()\n", " delta = calculate_deltaw(x_values, w_values)\n", " a = (left_side - right_side) / delta\n", " return a\n", " \n", "def calculate_bw(x_values, y_values, w_values):\n", " \"\"\"Calculates slope for weighted regression line\"\"\"\n", " left_side = w_values.sum() * (w_values * x_values * y_values).sum()\n", " right_side = (w_values * x_values).sum() * (w_values * y_values).sum()\n", " delta = calculate_deltaw(x_values, w_values)\n", " b = (left_side - right_side) / delta\n", " return b" ] }, { "cell_type": "code", "execution_count": 10, "id": "50e0190d-b7a1-4d18-b7a1-44cadc5ab6de", "metadata": {}, "outputs": [], "source": [ "# Modify the deviation for Quito to give it more weight in the regression line calculation\n", "data.loc[(6, \"Deviation\")] = 0.5" ] }, { "cell_type": "code", "execution_count": 11, "id": "ed50c31f-d71c-44b6-a193-85a7a060cbfa", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
CityLatitudeAverage temperatureMin temperatureMax temperatureDeviationWeight
0Helsinki60.195.01.08.03.500.081633
1Paris48.4412.48.816.03.600.077160
2Athens37.5418.514.322.03.850.067465
3Cairo31.2421.016.027.05.500.033058
4Chennai13.0028.524.232.94.350.052847
5Kuala Lumpur3.0727.023.031.04.000.062500
6Quito-0.0914.010.019.00.504.000000
7Guayaquil-2.0926.022.030.04.000.062500
8Ivalo68.361.0-5.02.03.500.081633
9Kathmandu27.4218.013.022.04.500.049383
\n", "
" ], "text/plain": [ " City Latitude Average temperature Min temperature \\\n", "0 Helsinki 60.19 5.0 1.0 \n", "1 Paris 48.44 12.4 8.8 \n", "2 Athens 37.54 18.5 14.3 \n", "3 Cairo 31.24 21.0 16.0 \n", "4 Chennai 13.00 28.5 24.2 \n", "5 Kuala Lumpur 3.07 27.0 23.0 \n", "6 Quito -0.09 14.0 10.0 \n", "7 Guayaquil -2.09 26.0 22.0 \n", "8 Ivalo 68.36 1.0 -5.0 \n", "9 Kathmandu 27.42 18.0 13.0 \n", "\n", " Max temperature Deviation Weight \n", "0 8.0 3.50 0.081633 \n", "1 16.0 3.60 0.077160 \n", "2 22.0 3.85 0.067465 \n", "3 27.0 5.50 0.033058 \n", "4 32.9 4.35 0.052847 \n", "5 31.0 4.00 0.062500 \n", "6 19.0 0.50 4.000000 \n", "7 30.0 4.00 0.062500 \n", "8 2.0 3.50 0.081633 \n", "9 22.0 4.50 0.049383 " ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Calculate weights as 1 over the deviation squared\n", "data[\"Weight\"] = 1 / data[\"Deviation\"]**2\n", "\n", "# Check DataFrame contents\n", "data" ] }, { "cell_type": "code", "execution_count": 12, "id": "19ccb6dd-423a-41ae-8029-2e879306e2fb", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAs0AAAHgCAYAAABelVD0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAA7EElEQVR4nO3dfXjcdZ3v/+e7N7SEAgUKtaW0KSoohVJKQLAgKajgzQqyoHDiqj9xo+vNLsrZg8esu7gaZRVdjjcsm6oL7smCUEVY9SAsNC13BVpAoYAgpS3QCmW1QAmEtnx+f8wk5GaS76STyUyS5+O65srMZ74z8877iuxrv/3M+xspJSRJkiT1b1ylC5AkSZKqnaFZkiRJymBoliRJkjIYmiVJkqQMhmZJkiQpg6FZkiRJyjCh0gUUY9q0aam2trbSZVSlF198kd12263SZYxo9rB09rB09nBo2MfS2cPS2cPSVbKHq1evfjaltG/v9RERmmtra1m1alWly6hKbW1t1NfXV7qMEc0els4els4eDg37WDp7WDp7WLpK9jAi1hdad3uGJEmSlMHQLEmSJGUwNEuSJEkZDM2SJElSBkOzJEmSlMHQLEmSJGUwNEuSJEkZDM2SJElSBkOzJEmSlMHQLEmSJGUwNEuSJEkZDM2SJElSBkOzJEmSlMHQLEmSJGUwNEuSJEkZDM2SJElSBkOzJEmSlMHQLEmSJGUwNEtlVF9fT319faXLkCRJJTI0S5IkSRkMzZIkSVIGQ7MkSZKUwdAsSZIkZTA0S5IkSRkMzZIkSVIGQ7MkSZKUwdAsSZIkZTA0S5IkSRkMzZIkSVIGQ7MkSZKUwdAsSZIkZTA0S5IkSRkMzZIkSVIGQ7MkSZKUwdAsSZIkZTA0S5IkSRkMzZIkSVIGQ7MkSZKUwdCsAdXX11NfX1/pMiRJkirK0CxJkiRlMDRLkiRJGQzNkiRJUgZDsyRJkpTB0CxJkiRlMDRLkiRJGQzNkiRJUoayheaImBwRd0XEbyJiTUR8Ob++d0TcGBGP5n/uVa4aJEmSpKFQzjPNHcCJKaXDgQXAKRFxDPAF4KaU0huBm/KPJUmSpKpVttCccrbmH07M3xJwKnB5fv1y4LRy1SBJkiQNhbLuaY6I8RFxH/AMcGNK6U5gekppE0D+537lrEGSJEkq1YRyvnlKaQewICKmAtdExKHFvjYiGoFGgOnTp9PW1laWGivl3HPPBeDiiy8u6X22bt1a1t5s2bIFYNT1v7ty9nAs9A/K/3c4FtjDoWEfS2cPS2cPS1eNPSxraO6UUtoSEW3AKcDTETEjpbQpImaQOwtd6DUtQAtAXV1dqq+vH45Sh83UqVMBKPX3amtrK/k9BjJUdVazcvZwLPQPyv93OBbYw6FhH0tnD0tnD0tXjT0s5/SMffNnmImIXYG3Aw8D1wEfyR/2EeDactUgSZIkDYVynmmeAVweEePJhfOrUkq/iIg7gKsi4hxgA3BmGWuQJEmSSla20JxS+i1wRIH1/wZOKtfnSpIkSUPNKwJKkiRJGQzNkiRJUgZDsyRJkpTB0CxJkiRlMDRLkiRJGQzNkiRJUgZDsyRJkpTB0CxJkiRlMDRLZdLa2srKlStZvnw5tbW1tLa2VrokSZK0kwzNUhm0trbS2NhIR0cHAOvXr6exsdHgLEnSCGVolsqgqamJ9vb2Hmvt7e00NTVVqCJJklQKQ7NUBhs2bBjUuiRJqm6GZqkMZs+ePah1SZJU3QzNUhk0NzdTU1PTY62mpobm5uYKVSRJkkphaJbKoKGhgZaWFiZNmgTAnDlzaGlpoaGhocKVSZKknTGh0gVIo1VDQwNLliwBoK2trbLFSJKkknimWZIkScpgaJYkSZIyGJolSZKkDIZmSZIkKYOhWZIkScpgaJYkSZIyGJolSZKkDIZmSZIkKYOhWZIkScpgaFa/WltbWblyJcuXL6e2tpbW1tZKlyRJklQRhmYV1NraSmNjIx0dHQCsX7+exsZGg7MkSRqTDM0qqKmpifb29h5r7e3tNDU1VagiSZKkyjE0q6ANGzYMal2SJGk0MzSroNmzZw9qXZIkaTQzNKug5uZmampqeqzV1NTQ3NxcoYokSZIqx9CsghoaGmhpaWHSpEkAzJkzh5aWFhoaGipcmSRJ0vCbUOkCVL0aGhpYsmQJAG1tbZUtRpIkqYI80yxJkiRlMDRL0ghSX19PfX19pcuQpDHH0CxJkiRlMDRXgJenliRJGlkMzcPMy1NLkiSNPIbmYeblqSVJkkYeQ/Mw8/LUkiRJI4+heZh5eWpJkqSRx9A8zLw8tSRJ0shjaB5mXp5akiRp5PEy2hXg5aklSZJGFs80S5IkSRkMzZIkSVIGQ7MkSZKUoWyhOSIOiIhlEfFQRKyJiL/Jr18QEU9FxH3527vLVYMkSZI0FMr5RcDtwHkppXsiYndgdUTcmH/un1NKF5XxsyVJkqQhU7bQnFLaBGzK338hIh4C9i/X50mSJEnlMiwj5yKiFjgCuBNYBHwmIj4MrCJ3NvpPBV7TCDQCTJ8+fdSNZtuyZQtQ+si5rVu3lrU3Q1VnNStnD8dC/6D8f4djQbE9HCt/UzvLv8XS2cPS2cPSVWMPyx6aI2IK8FPg3JTS8xHxL8BXgJT/+S3gY71fl1JqAVoA6urqUn19fblLHVZTp04FoNTfq62treT3GMhQ1VnNytnDsdA/KP/f4VhQbA/Hyt/UzvJvsXT2sHT2sHTV2MOyTs+IiInkAnNrSulnACmlp1NKO1JKrwJLgKPLWYMkSZJUqnJOzwjgh8BDKaVvd1uf0e2w9wMPlKsGSZIkaSiUc3vGIuAvgPsj4r782heBsyNiAbntGeuAT5SxBkmSJKlk5ZyecSsQBZ76Vbk+U5IkSSoHrwgoSZIkZTA0S5IkSRkMzZIkSVIGQ7MkSZKUYViuCCiNVdV2NSNJkrRzPNMsSZIkZTA0S5IkSRkMzZIkSVIGQ7MkSZKUwdAsSZIkZTA0S5IkSRkMzZIkSVIGQ7MkSZKUwdAsSZIkZTA0S5IkSRkMzZIkSVIGQ7MkSZKUYUKlC1B1a2trq3QJkiRJFeeZZkmSJCmDoVmSJEnKYGiWJEmSMhiaJUmSpAyGZkkapPr6eurr6ytdhiRpGBmaJUmSpAyGZklSRXjGXtJIYmiWJEmSMhiaJUmSpAyGZkmSJCmDl9GuEC9PLUmSNHJ4plmSJEnKYGiWJEmSMhiaJUmSpAyGZkmSJCmDoVmSJEnKYGiWJEmSMhiaJUmSpAyGZkkaIVpbW1m5ciXLly+ntraW1tbWSpckSWOGoVmSRoDW1lYaGxvp6OgAYP369TQ2NhqcJWmYGJolaQRoamqivb29x1p7eztNTU0VqkiSxhZDsySNABs2bBjUuiRpaBmaJWkEmD179qDWJUlDy9AsSSNAc3MzNTU1PdZqampobm6uUEWSNLYYmiVpBGhoaKClpYVJkyYBMGfOHFpaWmhoaKhwZZI0NkyodAGSpOI0NDSwZMkSANra2ipbjCSNMZ5pliRJkjIYmiVJkqQMhmZJkiQpQ9lCc0QcEBHLIuKhiFgTEX+TX987Im6MiEfzP/cqVw2SJEnSUCjnmebtwHkppTcDxwCfjohDgC8AN6WU3gjclH8sSZIkVa2yheaU0qaU0j35+y8ADwH7A6cCl+cPuxw4rVw1SJIkSUNhWEbORUQtcARwJzA9pbQJcsE6Ivbr5zWNQCPA9OnTHa/Uj61bt9qbEtnD0o21Hm7ZsgUY2rFvxfawHJ9dKZXso/pnD0tnD0tXjT0se2iOiCnAT4FzU0rPR0RRr0sptQAtAHV1dam+vr5sNY5kbW1t2JvS2MPSjbUeTp06FWBIf+die1iOz66USvZR/bOHpbOHpavGHpZ1ekZETCQXmFtTSj/LLz8dETPyz88AnilnDZIkSVKpyjk9I4AfAg+llL7d7anrgI/k738EuLZcNUiSJElDoZzbMxYBfwHcHxH35de+CFwIXBUR5wAbgDPLWIMkSZJUsrKF5pTSrUB/G5hPKtfnSpIkSUPNKwJKkiRJGQzNkiRJUgZDsyRJkpTB0CxJkiRlMDRLkiRJGQzNkiRJUgZDsyRJkpTB0CxJkiRlGDA0R8SxEfH9iPhtRGyOiA0R8auI+HRE7DlcRUqSRpfW1lZWrlzJ8uXLqa2tpbW1tdIlSdKA+g3NEfH/gI8DvwZOAWYAhwB/B0wGro2I9w1HkZKk0aO1tZXGxkY6OjoAWL9+PY2NjQZnSVVtoDPNf5FSOieldF1KaWNKaXtKaWtK6Z6U0rdSSvXA7cNUpyRplGhqaqK9vb3HWnt7O01NTRWqSJKyDRSap0bEot6LEXF8RLweIKX0bNkqkySNShs2bBjUuiRVg4FC88XACwXWX8o/J0nSoM2ePXtQ65JUDQYKzbUppd/2XkwprQJqy1aRJGlUa25upqampsdaTU0Nzc3NFapIkrINFJonD/DcrkNdiCSNBE59KF1DQwMtLS1MmjQJgDlz5tDS0kJDQ0OFK5Ok/k0Y4Lm7I+IvU0pLui9GxDnA6vKWJUnVp7+pD4CBb5AaGhpYsiT3f17a2toqW4wkFWGg0HwucE1ENPBaSK4DdgHeX+a6JKnqDDT1wdAsSaNbv6E5pfQ08NaIWAwcml/+ZUrp5mGpTJKqjFMfJGnsGuhMMwAppWXAsmGoRZKq2uzZs1m/fn3BdUnS6DbQFQHPjIifR8Q1EfHB4SxKkqqRUx8kaewa6Ezz+cDR+ft3Az8pfzmSVL069y2fc845dHR0MGfOHJqbm93PLEljwECh+f8CP87fv3oYapGkqufUB0kamwb6IuDFEbEbECmlrcNYkyRJklRV+g3NEREppRcHenH+mDT0ZUmSJEnVY6ArAi6LiM9GRI+vhUfELhFxYkRcDnykvOVJkiRJlTfQnuZTgI8BV0TEXGALuUtrjwduAP45pXRfuQuUJEmSKm2gPc0vA5cAl0TERGAa8FJKacsw1SZJkiRVhcyLmwCklLYBm8pciyRJklSVBtrTLEmSJAlDsyRJkpSpqNAcEXMi4u35+7tGxO7lLUuSJEmqHpmhOSL+ElgK/Gt+aRbw8zLWJEmSJFWVYs40fxpYBDwPkFJ6FNivnEVJkiRJ1aSY0NyRUnql80FETAC8CqAkSZLGjGJC8/KI+CKwa0S8A7ga+M/yliVJkiRVj2JC8/nAZuB+4BPAr4C/K2dRkiRJUjUZ8OImETEO+G1K6VBgyfCUJEmSJFWXAc80p5ReBX4TEbOHqR5JkiSp6hRzGe0ZwJqIuAt4sXMxpfS+slUlSZIkVZFiQvOXy16FJEkaUvX19QC0tbVVtA5ptMgMzSml5cNRiCRJklStMkNzRLzAa3OZdwEmAi+mlPYoZ2GSJElStSjmTPPu3R9HxGnA0eUqSJIkSao2xcxp7iGl9HPgxKEvRZIkSapOxWzPOL3bw3FAHV5GW5IkSWNIMdMz/qzb/e3AOuDUslQjSRqQkxAkqTKKCc0/SCnd1n0hIhYBz5SnJEmSJKm6FLOn+btFrvUQET+KiGci4oFuaxdExFMRcV/+9u7BFDus/vu/4f774cUXs4+VJEnSqNbvmeaIOBZ4K7BvRHy+21N7AOOLeO/LgO8BP+61/s8ppYsGWefw+9Wv4MMfzt2fPh1e/3o48MC+P1/3OoiobK2SJEkqq4G2Z+wCTMkf033s3PPAGVlvnFJaERG1JVVXSYsXwxVXwNq18NhjuZ/Ll0NrK6Ru34PcddeeIbr7/blzYdKkyv0OkiRJGhL9hub8lQCXR8RlKaX1Q/iZn4mIDwOrgPNSSn8qdFBENAKNANOnT6/Ml19e97rc7a1vfa2uV15h8tNPs+vGjUzeuJFdN25k102bmHz//ex6ww2Mf/nlrmNTBB3TpvHyzJm8NGMGL+2/Py/PmMFLM2fy8syZbNtjj5LPUm/dutUvBpXIHpZurPVwy5YtwNB+KW+s9RDsY7ntbH/tYensYemqsYfFfBGwPSK+CcwDJncuppR2ZlbzvwBfITey7ivAt4CPFTowpdQCtADU1dWl+vr6nfi4YZYSPP1019npWLuWyY89xuS1a5n6m9/A9df3PH6PPfqene78OXs2TJyY+ZFtbW2MiN5UMXtYurHWw6lTpwIM6e881noI9rHcdra/9rB09rB01djDYkJzK/AT4L3AJ4GPAJt35sNSSk933o+IJcAvduZ9qlZEwbPTXdrb4fHHX9vu0flzzRr4xS/glVdeO3b8+Fxw7h2mO+/vuefw/V6SJEljXDGheZ+U0g8j4m+6bdlYvjMfFhEzUkqb8g/fDzww0PGjTk0NzJuXu/X26qvw1FO5EN09UD/2GPzsZ/Dssz2P33tveP3rOWTKFLjxxp7Bev/9c6FbkiRJQ6KY0Lwt/3NTRLwH2AjMynpRRFwB1APTIuJJ4B+A+ohYQG57xjrgE4MveZQaNw4OOCB3O+GEvs8/91zBs9S7r1kDK1bAjh2vHbvLLlBbW3jbx4EHwm67DduvJUn9qbb9ipI0kGJC81cjYk/gPHLzmfcAPpf1opTS2QWWfzi48tRlzz1hwYLcrZs729qoP+44eOKJvts+HnsMbr8dnn++53sVGqHXed8RepIkSX0MGJojYjzwxpTSL4DngMXDUpUydW6Ov+CCC2DChNx4u7lz+x6YEvzxj4W3faxY0f8IvUJnqWtrYfLkvp8hSZI0yg0YmlNKOyLifcA/D1M9GmoRsM8+udtRR/V9vqMD1q8vfJb6pptyX17s/l777184UL/+9bnP8Cy1JEkahYrZnnF7RHyP3ASNrmtKp5TuKVtVGj6TJsFBB+VuvaUEzzxTOFBffz1s2tTz+CEYoSdJklSNignNnbPT/rHbWgJ2Zk6zhkBraysrV66ko6ODs846i29961s0NDQM/QdF5PY/T59e3Ai9zkD94IPwy1/mzmJ36hyh199ZakfoSZKkKpYZmlNK7mOuIq2trTQ2NtKRD6RPP/00jY2NAOUJzgPJGqG3cWPhs9QDjNAr+OVER+hJkqQKywzNETEd+BowM6X0rog4BDg2peQkjApoamqivfs+Y6C9vZ2mpqbhD80DGTcOZs3K3QqN0Hv++cJfTly1Cn76U9i+/bVjHaEnSZIqrJjtGZcB/wY05R8/Qm5/s6G5AjZs2DCo9aq1xx4FR+gBucDcfYRe92B9xx25mdXdTZ/e/7YPR+hJkqQhUExonpZSuioi/jdASml7ROzIepHKY/bs2axfv77g+qiRNULvT38qvO3DEXqSJKlMignNL0bEPuS+/EdEHENuZrMqoLm5mcbGxh5bNGpqamhubq5gVcMoIrf/ee+9Bx6h1ztQr13b7wi9BXvvDUce6Qg9SZLUr2JC8+eB64DXR8RtwL7AGWWtSv3q3Ld8zjnn0NHRwfTp08s3PWMk2pkRevfeW3iE3u679/1SYudPR+hJkjSmFDM9456IOAE4GAjgdymlbWWvTP1qaGhgyZIlQO6KgJ1XB1SGfkbo3dfWluth5wi93mepHaEnSdKYV8z0jMnAp4DjyG3RuCUiLk0pvVzu4tS/tra2Hj81BIoZoVdo20fWCL3egdoReiOe/7uTpLGnmO0ZPwZeAL6bf3w28O/AmeUqSqo63Ufove1tfZ93hJ4kSaNaMaH54JTS4d0eL4uI35SrIGlEKmaEXqGz1I7QkyRpRCgmNN8bEceklFYCRMRbgNvKW5Y0inQfoXfSST2f62+E3tq12SP0egdqR+hJklQ2xYTmtwAfjojOq2fMBh6KiPuBlFKaX7bqpNGuDCP0+t32MW2aZ6klSdpJxYTmU8pehaTCihmhVyhQ9zdCr79tH47QkyRpQMWMnFsfEXsBB3Q/PqV0TzkLk5Sh+wi9Y4/t+3yhEXpr1xY/Qq/7/alTh+3XkiSpGhUzcu4rwEeBx8hfFTD/88TylSWpZEM9Qq+/s9T77z88v48kSRVUzPaMDwCvTym9Uu5iJA2TwYzQ6x6oC43QmziRo6dPh0MP7Ruo586FKVOG7/eSJKlMignNDwBTgWfKW4qkqjHIEXpbV66k5plnHKEnVYnW1lZWrlxJR0cHtbW1NDc309DQUOmypBGtmND8dXJj5x4AujZBppTeV7aqJFWvAiP0HmxrY7/Oy7n/8Y+Ft330N0Jv7tzCgdoRetJOaW1tpbGxkY789xbWr19PY2MjgMFZKkExofly4J+A+4FXy1uOpBGvc4ReXV3f5155BdatK3z1xJtvhhdffO3YQiP0ut93hJ5UUFNTE+3dx1EC7e3tNDU1GZqlEhQTmp9NKX2n7JVIGv122WXnRuj9+te5Ly525wg9qaANGzYMal1ScYoJzasj4uvAdfTcnuHIOUlDp5gReuvW9Q3UhUbojRuXC86FArUj9DTKzZ49m/Xr1xdcl7TzignNR+R/HtNtzZFzkoZXTQ0cckju1luhEXqd93/+c9i8uefxvUfodb8/a1ZubrU0QjU3N9PY2Nhji0ZNTQ3Nzc0VrEoa+Yq5uMni4ShEknZaMSP0Hn+871nqfkboUVtb+Cy1I/Q0AnTuWz7nnHPo6Ohgzpw5Ts+QhkAxFzeZDnwNmJlSeldEHAIcm1L6Ydmrk6ShsMcecPjhuVtv27fDk0/2DdSPPVZ4hN5++/W/7WPGDL+cqKrQ0NDAkiVLAGhra6tsMdIoUcz2jMuAfwOa8o8fAX4CGJoljXwTJuTOLNfWdo3Q66H3CL3O+7feCldckdsa0skRepI0avUbmiNiQkppOzAtpXRVRPxvgJTS9ojYMWwVSlIlZY3QW7++8Fnq3iP0IDdCr7+z1N3nV0uSqs5AZ5rvAhYCL0bEPuS+/EdEHAM8N8DrJGls2GUXeOMbc7feUsp9AbFQoL7hhj4j9I6rqcm9T6EvJ86Z4wg9SaqwgUJz58a8z5MbN/f6iLgN2Bc4o9yFSdKIFpHb/7zffkWN0PvDihXM6uiAhx5yhJ4kVaGBQvO+EfH5/P1rgF+RC9IdwNuB35a5NkkavXqN0Pv94Yczq/NS5K++Cps2FT5LPdAIvUKB2hF6kjQkBgrN44EpvHbGuVNN+cqRJDFuXG7/8/77DzxCr3egvuce+NnPihuh13lzhJ4kFWWg0LwppfSPw1aJJKk4OzNCb+1aWLkStmzpebwj9CSpKMXsaZYkjRRZI/T+9KfC2z4GGqFXKFDPnesIPUljykChucB/bSVJI9pee+XG5w00Qq93oF67FpYtK26EXuf9adM8Sy1pVOk3NKeU/jichUiSKmxnRuitXVtwhB67797/lxMdoSdpBCrmioCSpLEua4TeSy8V/nLiww/Dr35V3Ai9zvuO0JNUhQzNkqTS7bprjxF6PfQeodc9WBcaobfXXv1/OdERepIqxNAsSSqvrBF6L7zQN0xnjdArFKgPPHDYfiVJY4+hWZJUWbvvnj1Cr9CXE++8s88IvbfutRe86U2FQ/XrXpcL8JK0EwzNkqTq1X2E3okn9n2++wi9tWt59tZbmfnSS4VH6E2e3P+XEx2hJymDoVmSNHL1GqH3SFsbMzsvR17KCL3ewdoRetKYZ2iWJI1OxYzQKxSoHaEnqQBDsyRp7Ok+Qu+YY/o+X2iE3tq1A4/QKxSoHaEnjRqGZkmSeitmhF6hs9TFjtDrvO8IPWnEKFtojogfAe8FnkkpHZpf2xv4CVALrAM+kFL6U7lqkCRpyHUfoXf88X2ff+GF3Fnq3oF6Z0boTZkybL+WpIGV80zzZcD3gB93W/sCcFNK6cKI+EL+8fllrEGSpOG1++4wf37u1tuOHfDEE323fTz2GNx1V24aSHf77df/lxMdoScNq7KF5pTSioio7bV8KlCfv3850IahWZI0Vowfnz1Cr9C2j8GO0KutHaZfSBo7IqVUvjfPheZfdNuesSWlNLXb839KKe3Vz2sbgUaA6dOnH3nllVeWrc6RbOvWrUzxn+9KYg9LZw9LZw+HxmjuY2zbxuRnnmHyU0+x66ZN7LpxI5M3bWLXp55i140bGf/yyz2Of3riRJ7cZRdmHn88L8+YwUszZ/LyzJm8NHMm2/bcs98ReqO5h8PFHpaukj1cvHjx6pRSXe/1qv0iYEqpBWgBqKurS/WdczfVQ1tbG/amNPawdPawdPZwaIzZPhYYoXfP97/PzJdfZsb998P11/c8fsqUfr+cuPzxxzlhLPZwCI3Zv8MhVI09HO7Q/HREzEgpbYqIGcAzw/z5kiSNPgVG6P3TTTcBufDBSy/BunV9t308/DD8v/8H3c5Sv80RelJBwx2arwM+AlyY/3ntMH++JEljz667wpvfnLv11n2E3tq1rL/pJmpffTUXrK+9Fp7pdX5rr736D9SO0NMoVs6Rc1eQ+9LftIh4EvgHcmH5qog4B9gAnFmuz5ckSUXoNUJv3Zw51Hb/Z/H+Rujdey9ccw1s2/basYVG6HW/7z5fjWDlnJ5xdj9PnVSuz5QkSUMsa4Tek0/2DdQDjdDr7yy1I/RU5ar2i4CSJKnKjR8Pc+bkbgON0OsdqG+7rfAIvblzCwfq2trcFhOpggzNkiSpPPbaC448Mnfr7ZVXYMOGvmep166FtjbYurXn8fvv3/9Z6mnT+h2hJw0VQ7MkSRp+u+wCb3hD7tZbSvDss4W3fdx4Izz1VM/jp0zpP1DPnp37LKlEhmZJklRdImDffXO3/Ai9HnqP0OsM1L/7XZ8RehQaodf9/l4Fr7Em9WFoliRJI0vWCL0//KHwWWpH6KkEhmZJkjR6jBsHM2fmbscf3/f5rVsLfzmxvxF6c+YUDtSO0BtzDM2SJGnsmDJl8CP01q4d1Ai9XTZvzp3xdoTeqGJoliRJgiEbofdWcITeKGRoliRpFGpra6t0CaNPkSP0Hrn+eg6aMOG1YD3QCL1C2z723dcRelXI0CxJklSqbiP0Nk6axEHdL0U+0Ai9//ovuPzynu/lCL2qZGiWJEkqp6wRei+/DI8/3jdQZ43QKxSsHaFXNoZmSZKkIXTuuecyderU4rfITJ5c/Ai97sH6uuuKH6F34IFwwAGO0CuBoVmSJKlaOUKvahiaJUmSRqpiRuj1DtRZI/QKBeoZM8b8CD1DsyRJ0mjUfYTe4sV9ny80Qm/tWrj9drjyyq4RekD/I/QOPDC3PgZG6BmaJUmSxqKBRuht2wbr1xc+S11ohN7Mmf1v+xglI/QMzZIkSepp4sSuEXp9dI7QKxSoixmh1/3+nDkjZoSeoVmSJEnF6z5C7y1v6ft8oRF6a9f2P0LvgAP6nJ2e8txz0H3WdRUwNEuSJGnoFDNCr9BZ6m4j9Oa97nXwl385zIUPzNAsSZKk4dF9hN5xx/V9futWePxxHm5r44jhr25AhmZJkiRVhylT4LDDeO6//7vSlfQxtgfuSZIkSUUwNEuSJEkZDM2SJElSBkOzJEmSlMHQLEmSJGUwNEuSJEkZDM2SJElSBkOzJEmSlMHQLEmSJGUwNEuSJEkZDM2SJElSBkOzJEmSlMHQLEmSJGUwNEuSJEkZDM2SJElSBkOzJEmSlMHQLEmSJGUwNEuSJEkZDM2SJElSBkOzJEmSlMHQLEmSJGUwNEuSJEkZDM2SJElSBkOzJEmSlMHQLEmSJGWYUIkPjYh1wAvADmB7SqmuEnVIkiRJxajkmebFKaUFBmZJkjRatLa28uCDD7J8+XJqa2tpbW2tdEkaIm7PkCRJGgKtra00Njaybds2ANavX09jY6PBeZSoVGhOwA0RsToiGitUgyRJ0pBpamqivb29x1p7eztNTU0VqkhDKVJKw/+hETNTShsjYj/gRuCzKaUVvY5pBBoBpk+ffuSVV1457HWOBFu3bmXKlCmVLmNEs4els4els4dDwz6Wzh7uvBNPPJFCuSoiuPnmmytQ0chVyb/DxYsXry60fbgioblHAREXAFtTShf1d0xdXV1atWrV8BU1grS1tVFfX1/pMkY0e1g6e1g6ezg07GPp7OHOq62tZf369X3W58yZw7p164a/oBGskn+HEVEwNA/79oyI2C0idu+8D7wTeGC465AkSRpKzc3N1NTU9Firqamhubm5QhVpKFViT/N04NaI+A1wF/DLlNL1FahjzCj0zxuXXnopP/7xjytQTeF6JEka6RoaGmhpaWHixIlA7gxzS0sLDQ0NFa5MQ2HY5zSnlNYChw/356qnT37yk5UuQZKkUaehoYFvfvObTJ06lba2tkqXoyHkyLkx6oILLuCii3LbyOvr6zn//PM5+uijOeigg7jlllsA2LFjB3/7t3/LUUcdxfz58/nXf/3XPu9z/vnnc8kll/R4329961ts3bqVk046iYULF3LYYYdx7bXX9nltW1sb733ve7sef+Yzn+Gyyy4DYPXq1ZxwwgkceeSRnHzyyWzatAmA73znOxxyyCHMnz+fs846a8j6IUmSNBBDswDYvn07d911FxdffDFf/vKXAfjhD3/Innvuyd13383dd9/NkiVLePzxx3u87qyzzuInP/lJ1+OrrrqKM888k8mTJ3PNNddwzz33sGzZMs4777yC3yguZNu2bXz2s59l6dKlrF69mo997GNd43ouvPBC7r33Xn77299y6aWXArBq1So+/vGPD0UbJEnSEKqvrx81XyytyGW0VX1OP/10AI488siub/jecMMN/Pa3v2Xp0qUAPPfcczz66KPMnTu363VHHHEEzzzzDBs3bmTz5s3stddezJ49m23btvHFL36RFStWMG7cOJ566imefvppXve612XW8rvf/Y4HHniAd7zjHUDujPeMGTMAmD9/Pg0NDZx22mmcdtppANTV1fGDH/xgqFohSZLUh6FZAEyaNAmA8ePHs337dgBSSnz3u9/l5JNPHvC1Z5xxBkuXLuUPf/hD15aJ1tZWNm/ezOrVq5k4cSK1tbW8/PLLPV43YcIEXn311a7Hnc+nlJg3bx533HFHn8/65S9/yYoVK7juuuv4yle+wpo1a5gwwT9jSZJUXm7PUL9OPvlk/uVf/qXrcqCPPPIIL774Yp/jzjrrLK688kqWLl3KGWecAeTOSu+3335MnDiRZcuW9Tu38sEHH6Sjo4PnnnuOm266CYCDDz6YzZs3d4Xmbdu2sWbNGl599VWeeOIJFi9ezDe+8Q22bNnC1q1by/XrS5IkdfEU3RjQ3t7OrFmzuh5//vOfL+p1H//4x1m3bh0LFy4kpcS+++7Lz3/+8z7HzZs3jxdeeIH999+/axtFQ0MDf/Znf0ZdXR0LFizgTW96U5/XHXDAAXzgAx9g/vz5vPGNb+SII44AYJdddmHp0qX89V//Nc899xzbt2/n3HPP5aCDDuJDH/oQzz33HCklPve5zzF16lRWrVrFpZde6hYNSZJUNobmMaD7FohCuo/EmTZtWtee5nHjxvG1r32Nr33ta5mfcf/99/d4PG3atILbK4AeZ4e/8Y1v8I1vfKPPMQsWLGDFihV91m+99dY+a+5pliRJ5eb2DEmSJCmDoVmSJEnKYGgeA5qbm5k3bx7z589nwYIF3HnnnUBuduKqVasqXF3/vv71r/OGN7yBgw8+mF//+tcFj7n66quZN28e48aN6/O7FPN6SZKkYrineZS74447+MUvfsE999zDpEmTePbZZ3nllVcqXVamBx98kCuvvJI1a9awceNG3v72t/PII48wfvz4Hscdeuih/OxnP+MTn/jETr1ekiSpGJ5pHuU2bdrEtGnTuuYwT5s2jZkzZ/Y57oorruCwww7j0EMP5fzzz+9anzJlCueddx4LFy7kpJNOYvPmzQA89thjnHLKKRx55JEcf/zxPPzww0Na97XXXstZZ53FpEmTmDt3Lm94wxu46667+hz35je/mYMPPninXy9JklQMQ/Mo9853vpMnnniCgw46iE996lMsX768zzEbN27k/PPP5+abb+a+++7j7rvv7hot9+KLL7Jw4ULuueceTjjhhK5LbDc2NvLd736X1atXc9FFF/GpT32qz/suW7aMBQsW9Lm99a1vzaz7qaee4oADDuh6PGvWLJ566qmif+9SXy9JktSd2zNGuSlTprB69WpuueUWli1bxgc/+EEuvPBCPvrRj3Ydc/fdd1NfX8++++4L5GYsr1ixgtNOO41x48bxwQ9+EIAPfehDnH766WzdupXbb7+dM888s+s9Ojo6+nz24sWLue+++3aq7pRSn7WIGLbXS5IkdWdoHgPGjx9PfX099fX1HHbYYVx++eU9QnOhgNmfiODVV19l6tSpmYF42bJlfO5zn+uzXlNTw+23395j7Zprruk6i/2DH/yAWbNm8cQTT3Q9/+STTxbcVtKfUl8vSZLUndszRrnf/e53PProo12P77vvPubMmdPjmLe85S0sX76cZ599lh07dnDFFVdwwgknALkLoyxduhSA//iP/+C4445jjz32YO7cuVx99dVALnT/5je/6fPZnWeae996B2aA97///V3P19XV8b73vY8rr7ySjo4OHn/8cR599FGOPvroon/vUl8vSZLUnWeaR7mtW7fy2c9+li1btjBhwgTe8IY30NLS0uOYGTNm8PWvf53FixeTUuLd7343p556KgC77bYba9as4cgjj2TPPffkJz/5CQCtra381V/9FV/96lfZtm0bZ511FocffviQ1T1v3jw+8IEPcMghhzBhwgS+//3vd02++PjHP84nP/lJ6urquOaaa/jsZz/L5s2bec973sOCBQv49a9/PeDrJUmSBisG80/zlVJXV5eqeZ5wJbW1tVFfX1+2958yZUqPy16PRuXu4VhgD0tnD4eGfSydPSzdggULmDp1Km1tbZUupeI6/5YG24tK/h1GxOqUUl3vdbdnSJIkSRkMzRrQaD/LLEmSVAxDsyRJkpTB0DwGTJkyJfOYiy++mPb29mGopq8tW7ZwySWXDMl7rVixgoULFzJhwoSuqR/dPf/88+y///585jOf6fc9rrrqKg455BDmzZvH//gf/wPoe6GWyZMnd10ARpIkjX6GZgE7F5p37NgxJJ89lKF59uzZXHbZZV1ht7cvfelLXeP0Cnn00Uf5+te/zm233caaNWu4+OKLgZ7j826++WZqamp45zvfOSQ1S5Kk6mdoHkM6v4l6xhln8KY3vYmGhgZSSnznO99h48aNLF68mMWLFwNwww03cOyxx7Jw4ULOPPPMrr3NtbW1/OM//iPHHXccV199Nddffz0LFy7k8MMP56STTgJyl97+2Mc+xlFHHcURRxzBtddeC8Bll13GqaeeyimnnMLBBx/cdTGTL3zhCzz22GMsWLCAv/3bvy3pd6ytrWX+/PmMG9f3T3v16tU8/fTTA4bdJUuW8OlPf5q99toLgP3226/PMUuXLuVd73oXNTU1JdUqSZJGDuc0jzH33nsva9asYebMmSxatIjbbruNv/7rv+bb3/42y5YtY9q0aTz77LN89atf5b/+67/Ybbfd+Kd/+ie+/e1v8/d///cATJ48mVtvvZXNmzezcOFCVqxYwdy5c/njH/8IQHNzMyeeeCI/+tGP2LJlC0cffTRvf/vbAbjrrrt44IEHqKmp4aijjuI973kPF154IQ888EC/Vxg8/vjjeeGFF/qsX3TRRV3vm+XVV1/lvPPO49///d+56aab+j3ukUceAWDRokXs2LGDCy64gFNOOaXHMVdeeSWf//zni/pcSZI0Ohiax5ijjz6aWbNmAbk5kuvWreO4447rcczKlSt58MEHWbRoEQCvvPIKxx57bNfzH/zgB7uOe9vb3sbcuXMB2HvvvYHcWerrrruOiy66CICXX36ZDRs2APCOd7yDffbZB4DTTz+dW2+9ldNOO23Amm+55ZZSfmUALrnkEt797ndzwAEHDHjc9u3befTRR2lra+PJJ5/k+OOP54EHHmDq1KkAbNq0ifvvv5+TTz655JokSdLIYWgeYyZNmtR1f/z48Wzfvr3PMSkl3vGOd3DFFVcUfI/ddtut67iIKPj6n/70pxx88ME91u+8884+xxd6fW9Dcab5jjvu4JZbbuGSSy5h69atvPLKK0yZMoULL7ywx3GzZs3imGOOYeLEicydO5eDDz6YRx99lKOOOgrIfUnw/e9/PxMnTizqcyVJ0ujgnmYBsPvuu3cF02OOOYbbbruN3//+9wC0t7d3bVvo7thjj2X58uU8/vjjAF3bM04++WS++93v0nm1yXvvvbfrNTfeeCN//OMfeemll/j5z3/OokWLenx2IbfcckvXl/C634oNzJC77PeGDRtYt24dF110ER/+8If7BGaA0047jWXLlgHw7LPP8sgjj3DggQd2PX/FFVdw9tlnF/25kiRpdDA0C4DGxkbe9a53sXjxYvbdd18uu+wyzj77bObPn88xxxzDww8/3Oc1++67Ly0tLZx++ukcfvjhXds2vvSlL7Ft2zbmz5/PoYceype+9KWu1xx33HH8xV/8BQsWLODP//zPqaurY5999mHRokUceuihJX8R8O6772bWrFlcffXVfOITn2DevHmZr/nRj37EddddB+QC/z777MMhhxzC4sWL+eY3v9m1nWTdunU88cQTA07fkCRJo1N0ng2sZnV1dWnVqlWVLqMqVfLa7IN12WWXsWrVKr73ve9VupQeRlIPq5U9LJ09HBr2sXT2sHQLFixg6tSptLW1VbqUiuv8WxpsLyr5dxgRq1NKdb3XPdMsSZIkZfCLgBo2H/3oR/noRz9a6TIkSZIGzTPNkiRJUgZDsyRJkpTB0CxJkjSELr74Yr8EOAoZmiVJkqQMhmZJkiQpg6FZkiRJQ661tZWVK1eyfPlyamtraW1trXRJJTE0S5IkaUi1trbS2NhIR0cHAOvXr6exsXFEB2dDsyRJkoZUU1MT7e3tPdba29tpamqqUEWlMzRLkiRpSG3YsGFQ6yOBoVmSJElDavbs2YNaHwkMzZIkSRpSzc3N1NTU9Firqamhubm5QhWVztAsSZKkIdXQ0EBLSwuTJk0CYM6cObS0tNDQ0FDhynbehEoXIEmSpNGnoaGBJUuWAIyKKyR6plmSJEnKYGiWJEmSMlQkNEfEKRHxu4j4fUR8oRI1SJIkScUa9tAcEeOB7wPvAg4Bzo6IQ4a7DkmSJKlYlTjTfDTw+5TS2pTSK8CVwKkVqEOSJEkqSiVC8/7AE90eP5lfkyRJkqpSJUbORYG11OegiEagEWD69OmjYlRJOWzdutXelMgels4els4eDg37WDp7WDp7+JotW7YAgx85V409rERofhI4oNvjWcDG3gellFqAFoC6urpUX18/LMWNNG1tbdib0tjD0tnD0tnDoWEfS2cPS2cPXzN16lSAQfejGntYie0ZdwNvjIi5EbELcBZwXQXqkCRJkooy7GeaU0rbI+IzwK+B8cCPUkprhrsOSZIkqVgVuYx2SulXwK8q8dmSJEnSYHlFQEmSJCmDoVmSJEnKYGiWJEmSMhiaJUmSpAyGZkmSJCmDoVmSJEnKYGiWJEmSMhiaJUmSpAyGZkmSJCmDoVmSJEnKYGiWJEmSMhiaJUmSpAyGZkmSJCmDoVmSJEnKYGiWJEmSMhiaJUmSpAyGZkmSJCmDoVmSJEnKYGiWJEmSMhiaJUmSpAyGZkmSJCmDoVmSJEnKYGiWJEmSMkyodAGSJEkandra2ipdwpDxTLMkSZKUwdAsSZIkZTA0S5IkSRkMzZIkSVIGQ7MkSZKUwdAsSZIkZTA0S5IkSRkMzZIkSVIGQ7MkSZKUwdAsSZIkZTA0S5IkSRkMzZIkSVIGQ7MkSZKUwdAsSZIkZTA0S5IkSRkMzZIkSVIGQ7MkSZKUwdAsSZIkZTA0S5IkSRkipVTpGjJFxGZgfaXrqFLTgGcrXcQIZw9LZw9LZw+Hhn0snT0snT0sXSV7OCeltG/vxRERmtW/iFiVUqqrdB0jmT0snT0snT0cGvaxdPawdPawdNXYQ7dnSJIkSRkMzZIkSVIGQ/PI11LpAkYBe1g6e1g6ezg07GPp7GHp7GHpqq6H7mmWJEmSMnimWZIkScpgaB6hIuKUiPhdRPw+Ir5Q6XpGioj4UUQ8ExEPdFvbOyJujIhH8z/3qmSN1S4iDoiIZRHxUESsiYi/ya/bxyJFxOSIuCsifpPv4Zfz6/ZwkCJifETcGxG/yD+2h4MQEesi4v6IuC8iVuXX7OEgRMTUiFgaEQ/n/7t4rD0cnIg4OP832Hl7PiLOrbY+GppHoIgYD3wfeBdwCHB2RBxS2apGjMuAU3qtfQG4KaX0RuCm/GP1bztwXkrpzcAxwKfzf3/2sXgdwIkppcOBBcApEXEM9nBn/A3wULfH9nDwFqeUFnQb72UPB+f/ANenlN4EHE7u79EeDkJK6Xf5v8EFwJFAO3ANVdZHQ/PIdDTw+5TS2pTSK8CVwKkVrmlESCmtAP7Ya/lU4PL8/cuB04azppEmpbQppXRP/v4L5P4PxP7Yx6KlnK35hxPzt4Q9HJSImAW8B/hBt2V7WDp7WKSI2AN4G/BDgJTSKymlLdjDUpwEPJZSWk+V9dHQPDLtDzzR7fGT+TXtnOkppU2QC4TAfhWuZ8SIiFrgCOBO7OOg5LcV3Ac8A9yYUrKHg3cx8L+AV7ut2cPBScANEbE6Ihrza/aweAcCm4F/y28T+kFE7IY9LMVZwBX5+1XVR0PzyBQF1hyDomEVEVOAnwLnppSer3Q9I01KaUf+nyJnAUdHxKEVLmlEiYj3As+klFZXupYRblFKaSG57X6fjoi3VbqgEWYCsBD4l5TSEcCLuBVjp0XELsD7gKsrXUshhuaR6UnggG6PZwEbK1TLaPB0RMwAyP98psL1VL2ImEguMLemlH6WX7aPOyH/T7lt5Pba28PiLQLeFxHryG1ROzEi/i/2cFBSShvzP58ht4f0aOzhYDwJPJn/lyKApeRCtD3cOe8C7kkpPZ1/XFV9NDSPTHcDb4yIufn/r+ws4LoK1zSSXQd8JH//I8C1Fayl6kVEkNu/91BK6dvdnrKPRYqIfSNiav7+rsDbgYexh0VLKf3vlNKslFItuf8G3pxS+hD2sGgRsVtE7N55H3gn8AD2sGgppT8AT0TEwfmlk4AHsYc762xe25oBVdZHL24yQkXEu8nt5xsP/Cil1FzZikaGiLgCqAemAU8D/wD8HLgKmA1sAM5MKfX+sqDyIuI44Bbgfl7bS/pFcvua7WMRImI+uS+1jCd38uKqlNI/RsQ+2MNBi4h64H+mlN5rD4sXEQeSO7sMuW0G/5FSaraHgxMRC8h9GXUXYC3w/5H/3zX2sGgRUUPu+1oHppSey69V1d+ioVmSJEnK4PYMSZIkKYOhWZIkScpgaJYkSZIyGJolSZKkDIZmSZIkKYOhWZIkScpgaJakMomIrYM4tj4i3trt8Scj4sP5+x+NiJk78fnrImLaIF+zND+/t7OmVRHxjW7Pt0XEqm6P6yKiLX//sIi4bLB1StJIYGiWpOpQD3SF5pTSpSmlH+cffhQYdGgerIiYB4xPKa3NL/0VcDwwPiLe1O3Q/SLiXb1fn1K6H5gVEbPLXaskDTdDsyQNo4j4s4i4MyLujYj/iojpEVELfBL4XETcFxHHR8QFEfE/I+IMoA5ozT+3a/czyL3O9O4TETfk3/tfgej2uR+KiLvy7/GvETG+QHkN9LxM7TggkbvyY3Rb/ybwd/38iv9J7rLWkjSqGJolaXjdChyTUjoCuBL4XymldcClwD+nlBaklG7pPDiltBRYBTTkn3tpgPf+B+DW/HtfR+7Ss0TEm4EPAotSSguAHeQCcm+LgNXdHv8AuB0Yl1J6qNv6HUBHRCwu8B6ryJ2dlqRRZUKlC5CkMWYW8JOImAHsAjw+hO/9NuB0gJTSLyPiT/n1k4AjgbsjAmBX4JkCr58BbO58kFL6NfDrfj7rq+TONp/fa/0ZhmEriSQNN880S9Lw+i7wvZTSYcAngMk78R7bee2/371fnwocH8Dl+TPVC1JKB6eULihw3EvF1pNSujl/7DG9npqcfx9JGlUMzZI0vPYEnsrf/0i39ReA3ft5Te/n1pE7cwzw593WV5DfdpH/ot5e+fWbgDMiYr/8c3tHxJwCn/MQ8IaifoucZuB/9Vo7CHhgEO8hSSOCoVmSyqcmIp7sdvs8cAFwdUTcAjzb7dj/BN7f+UXAXu9zGXBp5xcBgS8D/yf/Hju6Hfdl4G0RcQ/wTmADQErpQXJbKW6IiN8CN5LbitHbL8lN8ShKSulXdNvOkbc4/z6SNKpESoX+JU+SNNbkA/kycl8Y3JF1fIHXTwKWA8ellLYPdX2SVEmGZklSl4g4GXgopbRhJ177RmD/lFLbkBcmSRVmaJYkSZIyuKdZkiRJymBoliRJkjIYmiVJkqQMhmZJkiQpg6FZkiRJyvD/Aynt1qvP4UinAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Calculate weighted slope and intercept\n", "aw = calculate_aw(data[\"Latitude\"].values, data[\"Average temperature\"].values, data[\"Weight\"].values)\n", "bw = calculate_bw(data[\"Latitude\"].values, data[\"Average temperature\"].values, data[\"Weight\"].values)\n", "\n", "# Create figure and axes\n", "fig, ax = plt.subplots(1, 1, figsize=(12,8))\n", "\n", "# Plot average temperatures and deviations as a function of latitude\n", "ax.errorbar(data[\"Latitude\"], data[\"Average temperature\"], fmt='ko', yerr=data[\"Deviation\"])\n", "\n", "# Calculate regression line points\n", "x1w = -5.0\n", "x2w = 70.0\n", "y1w = aw + bw * x1\n", "y2w = aw + bw * x2\n", "\n", "# Plot weighted regression line\n", "ax.plot([x1w, x2w], [y1w, y2w], 'r-')\n", "\n", "# Add axis labels\n", "ax.set_xlabel(\"Latitude (°N)\")\n", "ax.set_ylabel(\"Temperature (°C)\")\n", "\n", "# Add regression line info\n", "ax.text(1, 2, f\"Line values:\\nSlope = {bw:.2f}\\nIntercept = {aw:.2f}\")\n", "\n", "# Enable plot gridlines\n", "ax.grid()" ] }, { "cell_type": "markdown", "id": "d3a5a3bd-bcc2-4c47-bfce-7a81f57c8993", "metadata": {}, "source": [ "### More information\n", "\n", "Additional introductory information about least-squares regressions can be found on the [Introduction to Quantitative Geology course page](https://introqg-site.readthedocs.io/en/latest/notebooks/L2/least-squares.html)." ] }, { "cell_type": "markdown", "id": "715202a1-2103-4f6a-917e-c98af4baf102", "metadata": {}, "source": [ "## Linear correlation\n", "\n", "The linear correlation coefficient $r$ measures how well a given set of points fit to a line, with values ranging from 1 for a perfect positive correlation to -1 for a perfect negative correlation. For a given set of x and y points, $r$ can be calculated as follows:\n", "\n", "$$\n", "\\large\n", "r = \\frac{\\sum{\\left(x_{i} - \\bar{x} \\right)\\left(y_{i} - \\bar{y} \\right)}}{\\sqrt{\\sum{\\left(x_{i} - \\bar{x} \\right)^2} \\sum{\\left(y_{i} - \\bar{y} \\right)^2}}},\n", "$$\n", "\n", "where $\\bar{x}$ and $\\bar{y}$ are the means of the values on the x and y axes, respectively." ] }, { "cell_type": "code", "execution_count": 13, "id": "2950c9d8-8797-4c28-8519-ec54c3a0d625", "metadata": {}, "outputs": [], "source": [ "# Import NumPy for square root function\n", "import numpy as np" ] }, { "cell_type": "code", "execution_count": 14, "id": "8e47bee8-1af9-4afb-bf19-d674d221055d", "metadata": {}, "outputs": [], "source": [ "# Rename plot x and y values for simplicity\n", "x = data[\"Latitude\"].values\n", "y = data[\"Average temperature\"].values" ] }, { "cell_type": "code", "execution_count": 15, "id": "45ead911-18a1-4cbe-bbc0-22cf8e46a548", "metadata": {}, "outputs": [], "source": [ "# Calculate correlation coefficient\n", "top = ((x - x.mean()) * (y - y.mean())).sum()\n", "bottom = np.sqrt(((x - x.mean())**2).sum() * ((y - y.mean())**2).sum())\n", "r = top / bottom" ] }, { "cell_type": "code", "execution_count": 16, "id": "8b7a82b0-9725-493c-aebd-c7084fe458a9", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Correlation coefficient r: -0.8181\n" ] } ], "source": [ "print(f\"Correlation coefficient r: {r:.4f}\")" ] }, { "cell_type": "code", "execution_count": 17, "id": "5ccb4be0-df98-4480-98a2-b57a95870011", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAs0AAAHgCAYAAABelVD0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABXbElEQVR4nO3dd3yV5fnH8c/NEI2oqCBlyFABBYQAAXEScIB7Kxj3iLt11NFSrahprbvijKOgRlFxUSeu4EQJggruAQg4sE7kJ4Levz+eEDYnkJycjM/79TqvnPOcdeVqar99vJ/rDjFGJEmSJK1cvUwXIEmSJFV3hmZJkiQpBUOzJEmSlIKhWZIkSUrB0CxJkiSlYGiWJEmSUmiQ6QLKo2nTprFdu3aZLqNa+vnnn1l33XUzXUaNZg8rzh5WnD2sHPax4uxhxdnDistkDydOnPhNjLHZssdrRGhu164dJSUlmS6jWiouLiY3NzfTZdRo9rDi7GHF2cPKYR8rzh5WnD2suEz2MIQwfUXHXZ4hSZIkpWBoliRJklIwNEuSJEkpGJolSZKkFAzNkiRJUgqGZkmSJCkFQ7MkSZKUgqFZkiRJSsHQLEmSJKVgaJYkSZJSMDRLkiRJKRiaJUmSpBQMzZIkSVIKhmZJkiQpBUOzJEmSlIKhWZIkSUrB0CxJkiSlYGiWJEmSUjA0S2mUm5tLbm5upsuQJEkVZGiWJEmSUjA0S5IkSSkYmiVJkqQUDM2SJElSCoZmSZIkKQVDsyRJkpSCoVmSJElKwdAsSZIkpWBoliRJklIwNEuSJEkpGJolSZKkFAzNkiRJUgqGZkmSJCkFQ7MkSZKUgqFZkiRJSsHQLEmSJKVgaJYkSZJSMDRLkiRJKRiaJUmSpBQMzVql3NxccnNzM12GJElSRhmaJUmSpBQMzZIkSVIKhmZJkiQpBUOzJEmSlIKhWZIkSUrB0CxJkiSlYGiWJEmSUkhbaA4hrB1CeCOE8FYIYWoIYVjp8Y1CCM+EED4q/blhumqQJEmSKkM6zzTPBwbEGLsD2cCgEEJf4HzguRhjB+C50seSJElStZW20BwTc0sfNiy9RWBfYGTp8ZHAfumqQZIkSaoMaV3THEKoH0KYDHwNPBNjfB1oHmP8AqD05ybprEGSJEmqqAbp/PAY429AdgihCfBwCKFred8bQsgH8gGaN29OcXFxWmrMlDPOOAOAa6+9tkKfM3fu3LT25vvvvweodf1fUjp7WBf6B+n/O6wL7GHlsI8VZw8rzh5WXHXsYVpD8yIxxu9DCMXAIOCrEEKLGOMXIYQWJGehV/SeQqAQICcnJ+bm5lZFqVWmSZMmAFT09youLq7wZ6xKZdVZnaWzh3Whf5D+v8O6wB5WDvtYcfaw4uxhxVXHHqZzekaz0jPMhBDWAXYB3gfGAEeVvuwo4NF01SBJkiRVhnSeaW4BjAwh1CcJ5/fHGB8LIbwG3B9COA6YARycxhokSZKkCktbaI4xvg30WMHx/wE7p+t7JUmSpMrmjoCSJElSCoZmSZIkKQVDsyRJkpSCoVmSJElKwdAsSZIkpWBoliRJklIwNEuSJEkpGJolSZKkFAzNUpoUFRUxfvx4xo0bR7t27SgqKsp0SZIkaQ0ZmqU0KCoqIj8/n/nz5wMwffp08vPzDc6SJNVQhmYpDYYOHcq8efOWOjZv3jyGDh2aoYokSVJFGJqlNJgxY8ZqHZckSdWboVlKgzZt2qzWcUmSVL0ZmqU0KCgoICsra6ljWVlZFBQUZKgiSZJUEYZmKQ3y8vIoLCykUaNGALRt25bCwkLy8vIyXJkkSVoTDTJdgFRb5eXlceuttwJQXFyc2WIkSVKFeKZZkiRJSsHQLEmSJKVgaJYkSZJSMDRLkiRJKRiaJUmSpBQMzZIkSVIKhmZJkiQpBUOzJEmSlIKheWUWLoSff850FZIkSaoGDM0r89hj0LIl/PGP8N57ma4mI4qKihg/fjzjxo2jXbt2FBUVZbokSZKkjDA0r8xmm8Hee8Mtt0DnztC/P9x/P/z6a6YrqxJFRUXk5+czf/58AKZPn05+fr7BWZIk1UmG5pXp1g3uvhtmzoTLLoNp0+DQQ6FNG7jgAvj880xXmFZDhw5l3rx5Sx2bN28eQ4cOzVBFkiRJmWNoTqVZMzjvPPj4Y3j8cejdGwoKoF072HdfePpp+P33TFdZ6WbMmLFaxyVJkmozQ3N51a8Pe+wB//0vfPppEqRfew0GDYKOHeHKK+F//8t0lZWmTZs2q3VckiSpNjM0r4l27eAf/0iWaNxzT3LB4DnnQKtWcOSRMH48xJjpKiukoKCArKyspY5lZWVRUFCQoYokSZIyx9BcEY0awZAh8OKL8M47cNxx8MgjsO220LMnFBbC3LmZrnKN5OXlUVhYSKNGjQBo27YthYWF5OXlZbgySZKkqmdorixdu8INN8CsWXDzzck65xNPTM4+n346vPtupitcbXl5efTt25d+/foxbdo0A7MkSaqzDM2Vbb31krA8eTK88grss09yxrlLF8jNhfvuqzNj6yRJkmoLQ3O6hADbbQd33ZWMrfvXv2DGDBg8GNq04bjPPmOTX37JdJWSapjc3Fxyc3MzXYYk1TmG5qrQrBmce24ytu6JJ6BPHw6bMYN7X389GVv31FO1cmydJElSbWForkr16sHuu1N06KGsv9ZaNADaPfYYRbvvDh06wBVXwDffZLpKSZIkLcPQXMUWbU89r3Rd8/Tffyd/rbUoatgwORvdujUccUQyA7qGj62TJEmqLQzNVWyF21P/+itDf/kFpkyB44+HRx9N1kP36FGjx9ZJkiTVFobmKrbK7am7dIHrr4fZs+GWW5InTjwx2TzltNNg6tQqrFSSJEmLGJqrWLm2p27cGPLzYdIkePVV2G8/uPXWZBZ0v34wapRj6yRJkqqQobmKrdb21CEkuwveeWeyacrllyfj64YMgU03haFDafTll1VUuSRJUt1laK5ia7w9ddOmcM458NFH8OST0LcvXHYZffPykg1UnnzSsXWSJElpYmjOgAptT12vHgwalFws+NlnzDjsMHjjDdhjj2Rs3eWXO7ZOkiSpkhmaa7I2bfjsuOOSnQZHjUqWbJx3HrRqlYyte/VVx9ZJkiRVAkNzbbDWWnDooVBcnIyty8+HMWNg++0hOzuZxOHYOkmSpDVmaK5tunSB4cOTCwcLC5PlHCedlIytO/XUJFRLkiRptaQtNIcQNg0hvBBCeC+EMDWE8KfS4xeFEGaFECaX3vZIVw11WuPGcMIJ8Oabye6C++0Ht98OW28NO+3k2DpJkqTVkM4zzQuBs2OMWwF9gVNDCJ1Ln7smxphdensijTUohGTSxp13JuPqrrgiOQu9aGzdX/8K06dnukpJkqRqLW2hOcb4RYzxzdL7PwHvAa3S9X0qh6ZN4c9/TsbWPfVUMgP6X/+C9u1h773hiSfgt98yXaUkSVK106AqviSE0A7oAbwObA+cFkI4EighORv93Qrekw/kAzRv3pzi4uKqKLXKfP/99wAV/r3mzp27Zp/RqBGccQaNhgyhxeOP0/Kxx1jrscf4vxYtmL333ny5++4saNKk0uqszta4h+VQF/oH6e1hXVHeHtaVv6k15d9ixdnDirOHFVcdexhimkeShRAaA+OAghjjQyGE5sA3QAQuAVrEGI9d1Wfk5OTEkpKStNZZ1XJzc4GK/w9fcXFx2WdVyK+/wiOPwE03JVM41loLDj6Y06ZOZcr661M8blzFv6OaqrQerkBl/edc3aWzh3VFeXtYV/6m1pR/ixVnDyvOHlZcJnsYQpgYY8xZ9nhap2eEEBoCDwJFMcaHAGKMX8UYf4sx/g7cCvRJZw0qp7XWgkMOgRdegKlT4cQT4b//5frJk7lt4kS4+Wb46adMVylJkpQR6ZyeEYDbgfdijFcvcbzFEi/bH3AGWnXTuTNcdx3Mns0VHTvyewhw8snJpimOrZMkSXVQOs80bw8cAQxYZrzc5SGEd0IIbwP9gTPTWIMqYt11ebxFC/J79oTx42H//RePrdtxR7j3Xpg/P9NVSpIkpV3aLgSMMb4MhBU85Yi5miYE2Gab5Hb11TBiRLL2+bDDoFkzOO64ZDlHu3aZrlSSJCkt3BFQq2fjjeHss+HDD+Hpp2G77eDyy2GzzWCvvRxbJ0mSaiVDs9ZMvXqw227JxI1p0+Bvf4OJE2HPPWGLLeCyy+DrrzNdpSRJUqUwNKviNt0ULr4YZsyA++9PNkv5y1+S43l58PLLkObRhpIkSelkaFbladgQDj4Ynn8e3n0XTjoJHn88uWiwe/dkHbRj6yRJUg1kaFZ6bLUV/PvfMGsW3HorNGgAp5wCLVsmP995J9MVVoni4mI3oZAkqRYwNCu91l0Xjj8+We88fjwceCD85z/QrVtyBvqeexxbJ0mSqj1Ds6rGorF1I0bAzJlw5ZXw5ZfJmudNN03WQE+blukqJUmSVsjQrKq3aGzdBx/A2LGw/faLx9btuWeyDtqxdZIkqRoxNCtz6tWDXXeFhx+G6dPhggvgzTeTec+bbw7//Kdj6yRJUrVgaFb10Lo1DBuWjK174IHkrPNf/5ocP+wwx9ZJkqSMMjSremnYEA46KBlb9957yaSNJ55ILhrs1g1uvBF+/DHTVUqSpDrG0Kzqa8st4dprk7F1t90Ga60Fp54KrVrBySfD229nukJJklRHGJpV/a27Lhx3HJSUwOuvJ2eiR4xINkzZYQcoKnJsnSRJSitDs2qOEKBPn2TO86xZcNVVyYWChx+erH0+/3z47LNMVylJkmohQ7Nqpo02grPOgvffT8bW7bhjMvt5882TsXWPPebYOkmSVGkMzarZFo2te+ihZHOUCy6ASZNg770dWydJkiqNoVm1x6KxddOnw+jRSWheNLZuyBB46SXH1kmSpDViaNYqFRcXU1xcnOkyVk/DhnDggfDcc8nyjVNPhSefhJ12gq23hhtucGydJElaLYZm1W6dOsE118Ds2XD77bD22nDaadCyJZx0Erz1VqYrlCRJNYChWXVDVhYce2wytu6NN+CQQ2DkSMjOpsdpp8Hdd8Mvv2S6SkmSVE0ZmlX39O4Nd9yRjK27+moa/vADHHEEbLopnHcefPpppiuUJEnVjKFZdddGG8GZZ/LGyJHwzDPJmuerroIttoA99oD//texdZIkCTA0S8nYul12gQcfTCZvXHhhstZ5n31gs82goAC++irTVaoayc3NJTc3N9NlSJKqkKFZWlKrVnDRRcnM59GjoUMH+NvfkqUbgwfDuHGOrZMkqQ4yNEsrsmhs3bPPLh5b9/TTkJsLXbvC9dfDDz9kukqpRvOMvaSaxNAspbJobN2sWckFhFlZcPrpyVnpE0+EyZMzXaEkSUozQ7NUXllZcMwxMGFCcjv0ULjzTujRA7bbDu66y7F1kiTVUoZmaU3k5CSbpcyenZyF/t//4Mgjky27zz0XPvkk0xVKkqRKZGjOkBq5PbWWt+GGcMYZybrnZ59N1jxffXUytm7QIBgzxrF1kiTVAoZmqTKEADvvnEzcmD49mcDxzjuw777Qvn0ytu7LLzNdpSRJWkOGZqmytWoFf/97MrbuwQeTCwkXja079FDH1kmSVAMZmqV0adgQDjgg2W3wgw+SiRvPPLN4bN3w4Y6tkySphjA0S1WhY8dkrfPMmcnYunXXhT/+EVq2hPx8mDQp0xVKkqRVMDRLVWnR2Lo33kjG1g0eDHffDT17wrbbJiPsHFsnSVK1Y2iWMmXR2LpZs+Daa+G77+Coo5Kxdeec49g6SZKqEUOzlGkbbgh/+hO89x489xz075/Mfl40tu7RR2HhwkxXKUlSnWZolqqLEGDAAHjgAZgxA4YNgylTYL/9krF1l17q2Lo6rqioiPHjxzNu3DjatWtHUVFRpkuSpDrD0CxVRy1bwoUXJmPrHnoIttoKLrggGVt3yCFQXOzYujqmqKiI/Px85s+fD8D06dPJz883OEtSFTE0S9VZgwaw//4wdix8+GEycePZZ5MlHF26OLauDhk6dCjz5s1b6ti8efMYOnRohiqSpLrF0CzVFB06wFVXJRcO/uc/sN56i8fWnXACvPlmpitUGs2YMWO1jkuSKpehWapp1lkHjj4aXn8dSkpgyBAoKoJevaBvXxg5Ev7v/zJdpSpZmzZtVuu4JKlyGZqlmqxXL7jtNpg9Oxlb9/33SaBu3Rr+/Gf4+OMMF6jKUlBQQFZW1lLHsrKyKCgoyFBFklS3GJql2qBJk8Vj655/PpnC8e9/J0s6Bg6ERx5xbF0Nl5eXR2FhIY0aNQKgbdu2FBYWkpeXl+HKJKluaJDpAiRVohCSiwT794cvvkjOQt9yS3IxYevWyZbdxx8PLVpkulKtgby8PG699VYAiouLM1uMJNUxnmmWaqsWLZIxddOmwcMPQ+fOyRi7Nm2SsXUvvODYOkmSysnQLNV2DRokG6Q8/XQytu5Pf0rG1g0YkATp665L1kJLkqSVMjRLdUmHDnDllcnYuhEjYIMNkhDdqlWybMOxdZIkrVDaQnMIYdMQwgshhPdCCFNDCH8qPb5RCOGZEMJHpT83TFcNklZinXXgqKNg/HiYOBEOOwzuvTeZxrHNNo6tkyRpGek807wQODvGuBXQFzg1hNAZOB94LsbYAXiu9LGkTOnZE269NTn7/O9/w48/JmPrWrWCs8+Gjz7KdIWSJGVc2kJzjPGLGOObpfd/At4DWgH7AiNLXzYS2C9dNUhaDU2aJDsMvvtucpHgLrsk6507doTddksuJnRsnSSpjgqxCq6eDyG0A14EugIzYoxNlnjuuxjjcks0Qgj5QD5A8+bNe40aNSrtddZEc+fOpXHjxpkuo0azhyu31v/+R4vHH6fFY4+x9pw5zG/alNl77cUXe+7Jr02blr2urvXwjDPOAODaa6+ttM8sbw/T8d2Zksk+auXsYcXZw4rLZA/79+8/McaYs+zxtIfmEEJjYBxQEGN8KITwfXlC85JycnJiSUlJWuusqYqLi8nNzc10GTWaPSyHhQvh8cfhxhth7NjFEzlOPhn696d43Lg61cNFv2tlzkou799hOr47UzLZR62cPaw4e1hxmexhCGGFoTmt0zNCCA2BB4GiGONDpYe/CiG0KH2+BfB1OmuQVAkaNIB9903G1n30EZxxRrLz4M47Q+fOtBo92rF1kqRaLZ3TMwJwO/BejPHqJZ4aAxxVev8o4NF01SApDbbYAq64AmbOTKZsNGlChxtugJYtk7F1EydmukJJkipdOs80bw8cAQwIIUwuve0BXAbsGkL4CNi19LGkmmaddeDII+G11ygpLITDD0/G1uXkQJ8+yRxox9ZJkmqJdE7PeDnGGGKM3WKM2aW3J2KM/4sx7hxj7FD689t01SCpaszt0AEKC2H27GTixty5cMwxydi6s85KdiKUJKkGc0dASZVngw3g9NNh6tRkbN2uu8Lw4dCpU3LfsXWSpBrK0Cyp8oUAublw333w+edwySXwwQdwwAHQrh0MG5aclZYkqYYwNEtKrz/8Af72N/j0U3j0UejaFS66CNq0gYMOgueegyqYFy9JUkUYmiVVjQYNYJ994KmnkrF1Z565eOfBrbaCa6+F777LdJWSJK2QoVlS1Vs0tm7WLLjzTthwwyREt2oFxx0HbmYkSapmDM2SMmftteGII+C112DSpOT+qFHQu3cytu4//4F58zJdpSRJhmZJ1UR2NtxyS3KB4PDh8PPPcOyxjq2TJFULhmZJ1csGG8Bpp8GUKVBcDAMHLh5bt8su8NBDjq2TJFW5VYbmEMK2IYQbQghvhxDmhBBmhBCeCCGcGkLYoKqKlFQHhQD9+iXLNT7/HC69NDnbfOCB0LZtMoFj1qxMV6k1VFRUxPjx4xk3bhzt2rWjqKgo0yVJ0iqtNDSHEJ4EjgeeBgYBLYDOwN+AtYFHQwj7VEWRkuq4P/wBhg6Fzz5LxtZ16wYXX5yE5wMPdGxdDVNUVER+fj7z588HYPr06eTn5xucJVVrqzrTfESM8bgY45gY4+wY48IY49wY45sxxqtijLnAq1VUpyRB/frJ2Lonn0zG1p11Fowblyzb2HJLuOYax9bVAEOHDmXeMhd4zps3j6FDh2aoIklKbVWhuUkIYftlD4YQdgwhbA4QY/wmbZVJ0qpsvjlcfjnMnAl33QUbb5yE6FatkgsIJ0zIdIVaiRkzZqzWcUmqDlYVmq8FflrB8f8rfU6SMm/tteHww+HVVxePrbv//mRkXe/ecMcdjq2rZtq0abNaxyWpOlhVaG4XY3x72YMxxhKgXdoqkqQ1tWhs3axZcP31SVg+7rjk7POZZ8IHH2S6QgEFBQVkZWUtdSwrK4uCgoIMVSRJqa0qNK+9iufWqexCJKnSbLABnHpqMrZu3DgYNAhuuCFZ97zLLvDgg7BgwRp9tFMfKi4vL4/CwkIaNWoEQNu2bSksLCQvLy/DlUnSyq0qNE8IIZyw7MEQwnHAxPSVJEmVJATYaSe4995kbF1BQXIB4UEHrdHYOqc+VJ68vDz69u1Lv379mDZtmoFZUrW3qtB8BnBMCKE4hHBV6W0cyRi6P1VJdZJUWZo3h7/+FT79FMaMSZZyLBpbd8AB8Oyz8Pvvq/wIpz5IUt210tAcY/wqxrgdMAyYVnobFmPcNsb4ZdWUJ0mVrH592HtveOIJ+PhjOPtseOkl2HXXZPnG1VfDt9+u8K1OfZCkuivlNtoxxhdijMNLb89XRVGSVCU22wz+9a9k6cZdd0GzZkmIbtUKjjlmubF1Tn2QpLprVTsCHhxCeCSE8HAI4dCqLEqSqtSisXWvvAKTJ8NRR8EDDyRj63Jy4PbbYd48pz5IUh22qjPN5wEHAAcC51ZNOZKUYd27w803w+zZycSNX36B44+HVq3ImzCBwosucuqDJNVBDVbx3N3AnaX3H6iCWiSp+lh/fTjlFDj5ZHj5ZbjxRrjxRvIWLGDLJk0Ys/nmDJs8GRo2zHSlkqQqsNLQHGO8NoSwLhBijHOrsCZJqj5CgB13TG5ffQV33EGrYcMY9u67yeSNE05Ibq1bZ7pSSVIarWpNc4gx/ryqwBxCCOkpS5KqoebN4S9/4bBttuEvXbtCjx5wySXQrh3svz8880zKsXWSpJppVWuaXwghnB5CWOqy8BDCWiGEASGEkcBR6S1Pkqqf30PgtY03hscfh08+gT//OVnCsdtu0KnTKsfWSZJqplWF5kHAb8C9IYTZIYR3QwifAh8BQ4BrYowjqqBGSaq+2reHyy6DmTOhqCg5G71obN3RR8Mbb0CMma5SklRBq9rc5JcY440xxu2BtsDOQM8YY9sY4wkxxslVVaQkVXuNGsFhhyVnnN96K5nz/OCDsM02ydi6226Dn3/OdJWSpDWUcnMTgBjjghjjFzHG79NcjyTVfN26JdM2Zs9Ofv76a3KxYKtW8Kc/wfvvZ7pCSdJqKldoliStgfXWS0bWvf12slX3nnvCTTfBVlvBgAHJBioLFmS6SklSORiaJSndQoAddkjWPM+cCf/8J3z2GRxySDK27sILaTRnTqarlCStQrlCcwihbQhhl9L764QQ1ktvWZJUS22yCZx/Pnz8MTz2GPTsCZdeSt/Bg5OxdWPHOrZOkqqhlKE5hHACMBq4pfRQa+CRNNYkSbVf/frJco3HHoNPPmHG4MHwyiswcGAytu6qq+B//8t0lZKkUuU503wqsD3wI0CM8SNgk3QWJUl1Svv2fHbCCfD558kSjj/8IZn93KoVHHUUvP66Y+skKcPKE5rnxxh/XfQghNAA8J/eklTZFo2te+ml5OLBY4+Fhx6Cvn2hVy+47TbW/u23TFcpSXVSeULzuBDCX4F1Qgi7Ag8A/01vWZJUx2299eKxdTfdBAsXwgknMPq11zj944/hvfcyXaEk1SnlCc3nAXOAd4ATgSeAv6WzKElSqfXWg5NOSjZMefllXtt4Y/aePRs6d4b+/eH++5M50JKktFplaA4h1APeiTHeGmM8OMZ4UOl9l2dIUlUKAbbfnoKttuLgvn2TrbunTYNDD03G1l1wQbImWpKUFqsMzTHG34G3QghtqqgeSVIKP6y1Fpx3XjK27vHHk226CwqgXTvYbz94+mnH1klSJSvP8owWwNQQwnMhhDGLbukuTJKUQv36sMce8N//wqefJkH61Vdh0CDo2BGuvNKxdZJUScoTmocBewEXA1ctcZMkVRft2sE//pEs0bjnHmjZEs45Z/HYuvHjHVtXx+Tm5pKbm5vpMqRaI2VojjGOW9GtKoqTJK2mRo1gyBB48UV45x047jh4+GHYdttk98Fbb4Wff850lZJU45RnR8CfQgg/lt5+CSH8FkL4sSqKkyRVQNeucMMNMGtWMrbu998hPz85C3366fDuu5muUJJqjPKcaV4vxrh+6W1t4EDg+vSXJkmqFIvG1k2enGzVvffeUFgIXbpAbq5j6ySpHMqzpnkpMcZHgAGVX4okKa1CgO22g7vvhpkz4V//ghkzkrF1bdrA3/6WPJYkLac8yzMOWOJ2UAjhMtxGW5JqtmbN4Nxzk7F1TzwBffokFxK2bw/77gtPPeXYOklaQoNyvGbvJe4vBKYB+6alGknSKhUXF1fuB9arB7vvntymT0+Wbdx2G4wZA5tvDieeCMccA02bVu73SlINU57lGbfFGI8pvZ0QYywAOqS7MElSFWvbNtkk5fPP4d57k3F1554LrVvDkUfCa685tk5SnVWe0Dy8nMeWEkK4I4TwdQhhyhLHLgohzAohTC697bE6xUqSqsBaa8HgwTBuXDK27vjj4ZFHkvXQPXsmZ6Pnzs10lZJUpVYamkMI24YQzgaahRDOWuJ2EVC/HJ89Ahi0guPXxBizS29PrFHVkqSq0bUrXH89zJ4NN9+cnGk+8cTkLLRj6yTVIas607wW0Jhk3fN6S9x+BA5K9cExxheBbyuhRklSpjVunITlSZOSrbr32Wfx2Lp+/eC++xxbJ6lWW+mFgKW7/o0LIYyIMU6vxO88LYRwJFACnB1j/G5FLwoh5AP5AM2bN6/8i19qiblz59qbCrKHFVfXevj9998DlXtRXo3r4XHH0fDAA/nDk0/S8r//ZZ3Bg/l1ww35Yo89mL333sxv3jzlR9jH9FrT/trDirOHFVcdexhiios6QgjNgHOBLsDai47HGFPOag4htAMeizF2LX3cHPiGZGTdJUCLGOOxqT4nJycnlpSUpHpZnVRcXExubm6my6jR7GHF1bUeLvpdK/Mf6DW6h7//DmPHwo03wuOPJ8f23BNOOQV22y2Z0LEC9jG91rS/9rDi7GHFZbKHIYSJMcacZY+X50LAIuB9oD0wjGTk3IQ1KSLG+FWM8bcY4+/ArUCfNfkcSVI1Uq8eDBqUjKn77DP4y1/g9deTMXYdOsDll8M332S6SkmqkPKE5o1jjLcDC2KM40rPDPddky8LIbRY4uH+wJSVvVaSVAO1aQOXXpqMrRs1CjbdFM47Lxlbd8QRyXpox9ZJqoHKs7nJgtKfX4QQ9gRmA61TvSmEcC+QCzQNIcwE/g7khhCySZZnTANOXP2SJUnV3lprJdtzH3ooTJ2aTN4YOTLZwrt7dzj5ZIofeyy5wFCSaoDynGm+NISwAXA28GfgNuDMVG+KMQ6JMbaIMTaMMbaOMd4eYzwixrh1jLFbjHGfGOMXFaxfklTddekCw4cnY+tuuQVCgJNOgpYt4bTTklAtSdXcKkNzCKE+0CHG+EOMcUqMsX+MsVeMcUwV1aeVyM3N9SIDSTVL48aQnw9vvpks09hvv2TL7q5dk7F1o0Y5tk5StbXK0Bxj/A3Yp4pqkSTVBSHAttvCnXfCzJnJhYIzZ8KQIcka6KFDYXplTjqVpIorz/KMV0MI14cQdgwh9Fx0S3tlkqTar2lTOOcc+OgjeOop6NsXLrsMNtsM9t4bnnwyGWknSRlWntC8HcmM5ouBq0pvV6azKK1aUVER48ePZ9y4cQwePJiioqJMlyRJFVOvHgwcCI8+unhs3YQJsMcesMUW8K9/wZw5ma5SUh2WMjSXrmNe9pZyYxOlR1FREfn5+cyfPx+Ar776ivz8fIOzpNpj0di6GTOS7bnbtIHzz0/G1h1+uGPrJGVEytAcQmgeQrg9hPBk6ePOIYTj0l+aVmTo0KHMmzdvqWPz5s1j6NChGapIktJkrbXgkEOguDiZsHHiifDf/8L220N2djLG7qefMl2lpDqiPMszRgBPAy1LH38InJGmepTCjBkzVuu4JNUKnTvDddfBrFlQWJgs5zj5ZGjVCk49Faa4V5ak9CpPaG4aY7wf+B0gxrgQ+C2tVWml2rRps1rHJalWadwYTjghGVv32muw//5w++2w9daw005s8txzULp8TZIqU3lC888hhI1JdvEjhNAX+CGtVWmlCgoKyMrKWupYVlYWBQUFGapIkjIghGTSxsiRydnnK66A2bPpfOmlydi6v/4Vpk3LdJWSapHyhOazgDHA5iGEV4A7gdPTWpVWKi8vj8LCQho1agRA8+bNKSwsJC8vL8OVSVKGbLwx/PnP8OGHvPWvf8F22yXTNhaNrXviCfjNf0EqqWLKMz3jTaAfyei5E4EuMca3012YVi4vL4++ffvSr18/Ro0aZWCWJIB69fiuTx945JHkLPPQoVBSAnvu6dg6SRVWnukZawN/BC4BhgGnlh5TBhUXF1NcXJzpMqQ6yf/+1QCbbgqXXJKMrbv/fmjXbvHYurw8eOUVx9ZJWi3lWZ5xJ8nmJsOB64HOwF3pLEqSpErRsCEcfDC88MLisXWPPQY77ADdu8NNNzm2TlK5lCc0d4oxHhdjfKH0lg90THdhkiRVqkVj62bPhltvhQYN4JRToGXL5Oc772S6QknVWHlC86TSiRkAhBC2AV5JX0mSJKXRuuvC8cfDxIkwfjwceCDccQd06wY77gj33OPYOknLKU9o3gZ4NYQwLYQwDXgN6BdCeCeE4AWBkqSaKQTYZhsYMSIZW3fllfDll8ma5003hb/8xbF1ksqUJzQPAtqTTNDoV3p/D2AvYO/0lSZJUhXZeGM4+2z44AN4+ulkq+7LL0/G1u25Jzz+uGPrpDquPCPnpgM/AhsAGy+6xRinlz4nSVLtUK8e7LYbPPxwcpb5b39Ldh/cay/YfHP45z/h668zXaWkDCjPyLlLgLeB64CrSm9XprkuSZIya9NN4eKLF4+t22yzZKfB1q3hsMPg5ZcdWyfVIeVZnnEIsHmMMTfG2L/0NiDdhUmSVC0sGlv3/PPw7rtw8snJLoM77phcPHjjjfDjj5muUlKalSc0TwGapLkOSZKqv622gn//O7lw8LbbYK214NRToVWrJEy/XT2ujy8qKmL8+PGMGzeOdu3aUVRUlOmSpBqvPKH5nyRj554OIYxZdEt3YZIkVVvrrgvHHZds0/3668nYuhEjkg1TdtgBiooyNrauqKiI/Px85pd+//Tp08nPzzc4SxVUntA8EvgXcBmL1zRflc6iJEmqEUKAPn0Wj6276ir46is4/PBk7fP558Nnn1VpSUOHDmXevHlLHZs3bx5Dhw6t0jqk2qY8ofmbGON1pbsBjlt0S3tlkiTVJBttBGedlYytGzs2WfN8xRXJ1I0990y2766CsXUzZsxYreOSyqc8oXliCOGfIYRtQwg9F93SXpkkSTVRvXqw667w0EMwfTpccAFMmgR7710lY+vatGmzWscllU95QnMPoC/wDxw5J0lS+bVuDcOGJeH5gQeS0LxobN2QIfDSS5U+tq6goICsrKyljmVlZVFQUFCp3yPVNeXZ3KT/Cm6OnJMkqbwaNoSDDoLnnoP33oNTToEnn4SddoKtt4Ybbqi0sXV5eXkUFhbSqFEjANq2bUthYSF5eXmV8vlSXVWezU2ahxBuDyE8Wfq4cwjhuPSXJklSLbTllnDttYvH1q29Npx2GrRsCSedBG+9VeGvyMvLo2/fvvTr149p06YZmKVKUJ7lGSOAp4GWpY8/BM5IUz2SJNUNS46te+ONZAOVkSMhOxu23x7uvht++SXTVUoqtdLQHEJoUHq3aYzxfuB3gBjjQiD9l/9KklRX9O4N//lPcvb56qthzhw44ohkK+/zzoNPP810hVKdt6ozzW+U/vw5hLAxEAFCCH2BH9JdmCRJdc5GG8GZZ8L778MzzyRj6666CrbYAvbYA/773yoZWydpeasKzaH051nAGGDzEMIrwJ3A6ekuTJKkOqtePdhll8Vj6y68ECZPhn32gc02g3/8I9lERVKVWVVobhZCOAvIBR4GLgeeBG4Fdkl/aZIkiVat4KKLkvA8ejR06ABDhyZLNwYPhhdfrPSxdZKWt6rQXB9oDKwHrAs0KD2WVXpMkiRVlYYN4cAD4dlnk+Ubp54KTz8N/fpV+tg6SctrsIrnvogxXlxllUiSpPLp1AmuuQYKCmDUKLjppmRs3XnnQV4enHxypiuUap3yrGmWJEnVUVYWHHssTJiQjK075BC4807o0YPrJ01i16++cmydVElWFZp3rrIqJElSxfTuDXfcUTa2bv0FCxj6/vvJlt3nnguffJLpCqUabaWhOcb4bVUWIkmSKkHp2Loje/fmrG7dkjXPV1+djK3bfXcYM8axddIaKM+OgJIkqaYJgTc33BAefDCZvHHRRfD227DvvsnYuoIC+PLLTFcp1RiGZkmSartWreDvf4dp05IQ3bEj/O1vi8fWjRvn2DopBUOzJEl1RcOGcMAByW6DH3wAp5+ejK3LzYWuXeH666k/d26mq5SqJUOzJEl1UceOyVrnWbOSCwjXXRdOP53tDj4Y8vOTHQgllTE0S5JUl2VlwTHHJCPrJkzg6/794e67oUcP2HZbuOsux9ZJGJolSdIiOTl8cO65ydnna66Bb7+FI49Mxtadc45j61SnGZolSdLSNtwQzjgj2a77ueeSNc/XXJOMrRs0CB59FBYuzHSVUpUyNEuSpBULAQYMgNGjF4+te+cd2G+/ZGzdpZc6tk51hqFZkiSltmhs3fTp8NBD0KkTXHBBMrbu0EOhuNixdarV0haaQwh3hBC+DiFMWeLYRiGEZ0IIH5X+3DBd3y9JktKgQQPYf//FY+v++Mfkfv/+0KULDB8OP/yQ6SqlSpfOM80jgEHLHDsfeC7G2AF4rvSxJEmqiTp2hKuuSi4c/M9/YL31khDdsmUytm7SpExXKFWatIXmGOOLwLfLHN4XGFl6fySwX7q+X5IkVZF11oGjj4bXX4eSEhgyJBlb17Mn9O0Ld97p2DrVeCGmcf1RCKEd8FiMsWvp4+9jjE2WeP67GOMKl2iEEPKBfIDmzZv3GjVqVNrqrMnmzp1L48aNM11GjWYPK84eVpw9rBz2cbEzzjgDgGuvvXa13ldZPWzw0080f/ppWo0ZQ9bnn7Ng/fX5ctAgZu+zD//XqlWFP7868++w4jLZw/79+0+MMeYse7zahuYl5eTkxJKSkrTVWZMVFxeTm5ub6TJqNHtYcfaw4uxh5bCPiy3qQ3Fx8Wq9r9J7GCO88ALcdBM88kgyqm633eCUU2DPPZM10rWMf4cVl8kehhBWGJqrenrGVyGEFqUFtQC+ruLvlyRJVWnR2LoHHkgmbwwbBlOnJmPr2reHSy6BL77IdJVSSlUdmscAR5XePwp4tIq/X5IkZUrLlnDhhTBtGjz8MGy1VfK4TRs45JDkjLRj61RNpXPk3L3Aa0CnEMLMEMJxwGXAriGEj4BdSx9LkqS6pEGD5Ezz2LHw4YfJxI1nn03OSHfpAtddB99/n+kqpaWkc3rGkBhjixhjwxhj6xjj7THG/8UYd44xdij9uex0DUmSVJd06LB4bN2IEcnYuj/9KdlM5YQT4M03M12hBLgjoCRJqg7WWQeOOioZWzdxIhx2GBQVQa9eydi6kSPh//4v01WqDjM0S5Kk6qVnT7j1Vpg9G/7972SHwaOPhtat4c9/ho8/znSFqoMMzZIkqXpq0iRZ7/zuu/D887DzzkmI7tAhGVu3aISdVAUMzZIkqXoLAfr3h/vvhxkz4OKL4b33YP/9k7F1F1/s2DqlnaFZkiTVHC1awAUXwGefJWPrOneGv/89GVt38MGOrVPaGJolSVLNs2hs3dNPJ2Pr/vSnZAnHgAFJkP73vx1bp0plaJYkSTVbhw5w5ZUwc2YyZWODDeCMM5LNVI4/PpnGIVWQoVmSJNUO66wDRx4J48cnQTkvD+69F3JyYJttkjnQjq3TGjI0S5JUCxUXF1NcXJzpMjJn0di6WbOSHQZ//BGOOSbZNOXss+GjjzJdoWoYQ7MkSaq9mjSB009Pxta98ALsumsSojt2TMbWPfywY+tULoZmSZJU+4UAublw333J2LpLLoH334cDDoB27ZKxdbNnZ7pKVWOGZkmSVLe0aAF/+xt8+mmyQUrXrovH1h10UDKFowJj68444wxyc3MrrVxVD4ZmSZJUNzVoAPvuC089laxxPvPMZAnHzjvDVlvBtdfCd99lukpVE4ZmSZKkLbaAK65YPLZuww2TEN2qFRx3nGPrZGiWJEkqs2hs3WuvwZtvwuGHw6hRydi6Pn3gP/+BefMyXaUywNAsSZK0Ij16QGFhcoHg8OEwdy4ceyy0bg1nnZXsRKg6w9AsSZK0KhtsAKedBlOnQnFxMrZu+HDo1Cm5/9BDjq2rAwzNkiRJ5REC9OuXjK37/HO49FL44AM48EBo2xaGDXNsXS1maJYkSVpdf/gDDB0Kn30Gjz4K3boloblNG66YNo2e331XobF1qn4MzZIkSWuqfn3YZx948slkbN1ZZ5Ezdy5Xv/02bLmlY+tqEUOzJElSZdh8c7j8cgZ27kzBllvCxhsvHlt37LFQUpLpClUBhmZJkqRK9Gu9ejzTvDm8+ipMmgRHHAH33w+9eye3O+5wbF0NZGiWJElKl+xsuOUWmDULrr8+CcvHHZecfT7zzORCQtUIhmZJkqR022ADOPVUmDIFxo2DQYPghhuSdc+77AIPPggLFmS6Sq2CoVmSJKmqhAA77QT33puMrSsoSC4gPOggaNcOLrooOSutasfQLEmSlAnNm8Nf/wqffgpjxkD37nDxxcnM5wMPhGefhd9/z3SVKmVoliRJyqT69WHvveGJJ+Djj+Hss+HFF5PdBrfcEq65xrF11YChWZIkqbrYbDP417+SpRt33QXNmsFZZ0HLlsnYugkTMl1hnWVoliRJqm7WXhsOPxxeeQUmT4ajjkrG1vXpAzk5jq3LAEOzJElSdda9O9x8M8yenUzc+OWXxWPrzjjDsXVVxNAsSZJUE6y/PpxyCrzzTrLmedAguPHGZN3zzjvD6NGOrUsjQ7MkSVJNEgLsuOPisXX/+EdyAeHBByeTN/7+d5g5M9NV1jqGZkmSpJqqeXP4y1+SsXX//S/06AGXXJLMfD7gAHjmGcfWVRJDsyRJUk1Xvz7stRc8/jh88gn8+c/w0kuw227J8o2rr4Zvv810lTWaoVmSJKk2ad8eLrssWaJx992wySbJ7OdWreCYY+CNNyDGTFdZ4xiaJUmSaqNGjSAvD15+Gd56C44+OrlYcJttkrF1t9/u2LrVYGiWJEmq7bp1g5tuglmzkrF1v/4Kxx+fbJpyxhnw/vuZrrDaMzRLkiTVFYvG1r39drLmeY89krF1W20FAwY4tm4VDM2SJEl1TQiwww5wzz3J2ud//jOZwLFobN2FFzq2bhmGZkmSpLpsk03g/POTqRuPPZaMrbv00mRs3f77w9ixjq3D0CxJkiRIxtbtuefisXXnnJNcRDhwIHTqBFddBf/7X6arzBhDsyRJkpbWvn2yZGPmTCgqgj/8IZn93Lp1MoXj9dfr3Ng6Q7MkSZJWrFEjOOyw5KLBt99O5jw/+CD07ZuMrbvtNvj550xXWSUMzZIkSUpt662TSRuzZyc/FyyAE05INk3505/gvfcyXWFaGZolSZJUfuutByefnGyY8vLLyTrom2+Gzp2hf3944IFaObbO0CxJkqTVFwJsv32y5vnzz5M10NOmwSGHQJs2ydi6zz/PdJWVJiOhOYQwLYTwTghhcgihJBM1SJIkqZIsGlv38cfJ9I1evRaPrdtvv1oxti6TZ5r7xxizY4w5GaxBkiSp0hQVFfHuu+8ybtw42rVrR1FRUaZLqlr16ye7DD72WLJZynnnwauvJmPrOnaEK6+ssWPrXJ4hSZJUCYqKisjPz2dB6Xre6dOnk5+fX/eC8yLt2sE//pEs0bjnHmjZMpn93KoVHHUUjB9fo8bWZSo0R2BsCGFiCCE/QzVIkiRVmqFDhzJv3ryljs2bN4+hQ4dmqKJqolEjGDIEXnwxGVt33HHw0EOw7bbJMo5bb60RY+tCzEDCDyG0jDHODiFsAjwDnB5jfHGZ1+QD+QDNmzfvNWrUqCqvsyaYO3cujRs3znQZNZo9rDh7WHH2sHLYx4qzh2tuwIABrChXhRB4/vnnM1BR9VV/3jyaP/ssLR99lMaffsrCddfly912Y/a++zKvbduM/h32799/4oqWD2ckNC9VQAgXAXNjjFeu7DU5OTmxpMTrBVekuLiY3NzcTJdRo9nDirOHFWcPK4d9rDh7uObatWvH9OnTlzvetm1bpk2bVvUF1QQxJmueb7opGVX366/Qrx9T+/Wjy4UXJmukq1gIYYWhucqXZ4QQ1g0hrLfoPrAbMKWq65AkSapMBQUFZGVlLXUsKyuLgoKCDFVUAywaW3f33cmW3ZddBjNm0P4//4F61evSu0xU0xx4OYTwFvAG8HiM8akM1FFnrOhfb9x8883ceeedGahmxfVIklTT5eXlUVhYSMOGDYHkDHNhYSF5eXkZrqyGaNYsmbbx8ce8dcUVSaCuRhpU9RfGGD8Fulf192ppJ510UqZLkCSp1snLy+OKK66gSZMmFBcXZ7qcmqlePeY3b57pKpZTvc57q8pcdNFFXHllsow8NzeX8847jz59+tCxY0deeuklAH777TfOOeccevfuTbdu3bjllluW+5zzzjuPG2+8canPveqqq5g7dy4777wzPXv2ZOutt+bRRx9d7r3FxcXstddeZY9PO+00RowYAcDEiRPp168fvXr1YuDAgXzxxRcAXHfddXTu3Jlu3boxePDgSuuHJEnSqhiaBcDChQt54403uPbaaxk2bBgAt99+OxtssAETJkxgwoQJ3HrrrXz22WdLvW/w4MHcd999ZY/vv/9+Dj74YNZee20efvhh3nzzTV544QXOPvvsFV5RvCILFizg9NNPZ/To0UycOJFjjz22bFzPZZddxqRJk3j77be5+eabASgpKeH444+vjDZIkqRKlJubW2suLK3y5Rmqng444AAAevXqVXaF79ixY3n77bcZPXo0AD/88AMfffQR7du3L3tfjx49+Prrr5k9ezZz5sxhww03pE2bNixYsIC//vWvvPjii9SrV49Zs2bx1Vdf8Yc//CFlLR988AFTpkxh1113BZIz3i1atACgW7du5OXlsd9++7HffvsBkJOTw2233VZZrZAkSVqOoVkANGrUCID69euzcOFCAGKMDB8+nIEDB67yvQcddBCjR4/myy+/LFsyUVRUxJw5c5g4cSINGzakXbt2/PLLL0u9r0GDBvy+xD70i56PMdKlSxdee+215b7r8ccf58UXX2TMmDFccsklTJ06lQYN/DOWJEnp5fIMrdTAgQO56aabyrYD/fDDD/l5BTv2DB48mFGjRjF69GgOOuggIDkrvckmm9CwYUNeeOGFlc6tfPfdd5k/fz4//PADzz33HACdOnVizpw5ZaF5wYIFTJ06ld9//53PP/+c/v37c/nll/P9998zd+7cdP36kiRJZTxFVwfMmzeP1q1blz0+66yzyvW+448/nmnTptGzZ09ijDRr1oxHHnlkudd16dKFn376iVatWpUto8jLy2PvvfcmJyeH7Oxsttxyy+Xet+mmm3LIIYfQrVs3OnToQI8ePQBYa621GD16NH/84x/54YcfWLhwIWeccQYdO3bk8MMP54cffiDGyJlnnkmTJk0oKSnh5ptvdomGJElKG0NzHbDkEogVWXIkTtOmTcvWNNerV49//OMf/OMf/0j5He+8885Sj5s2bbrC5RXAUmeHL7/8ci6//PLlXpOdnc2LL7643PGXX355uWOuaZYkSenm8gxJkiQpBUOzJEmSlIKhuQ4oKCigS5cudOvWjezsbF5//XUgmZ1YUlKS4epW7p///CdbbLEFnTp14umnn17hay644IKy32u33XZj9uzZq/V+SZKk8nBNcy332muv8dhjj/Hmm2/SqFEjvvnmG3799ddMl5XSu+++y6hRo5g6dSqzZ89ml1124cMPP6R+/fpLve6cc87hkksuAZLdAi+++GJuvvnmcr9fkiSpPDzTXMt98cUXNG3atGwOc9OmTWnZsuVyr7v33nvZeuut6dq1K+edd17Z8caNG3P22WfTs2dPdt55Z+bMmQPAJ598wqBBg+jVqxc77rgj77//fqXW/eijjzJ48GAaNWpE+/bt2WKLLXjjjTeWe936669fdv/nn38mhLBa75ckSSoPQ3Mtt9tuu/H555/TsWNHTjnlFMaNG7fca2bPns15553H888/z+TJk5kwYULZaLmff/6Znj178uabb9KvX7+yLbbz8/MZPnw4EydO5Morr+SUU05Z7nNfeOEFsrOzl7ttt912KeueNWsWm266adnj1q1bM2vWrBW+dujQoWy66aYUFRVx8cUXr/b7JUmSUjE013KNGzdm4sSJFBYW0qxZMw499FBGjBix1GsmTJhAbm4uzZo1o0GDBuTl5ZWNe6tXrx6HHnooAIcffjgvv/wyc+fO5dVXX+Xggw8mOzubE088kS+++GK57+7fvz+TJ09e7vbqq6+mrDvGuNyxRWeRl1VQUMDnn39OXl4e119//Wq/X5IkKRXXNNcB9evXJzc3l9zcXLbeemtGjhzJ0UcfXfb8igLmyoQQ+P3332nSpAmTJ09e5WtfeOEFzjzzzOWOZ2VlLRecH3744bKz2LfddhutW7fm888/L3t+5syZK1xWsqTDDjuMPffck2HDhq3R+yVJklbGM8213AcffMBHH31U9njy5Mm0bdt2qddss802jBs3jm+++YbffvuNe++9l379+gHJxiijR48G4J577mGHHXZg/fXXp3379jzwwANAErrfeuut5b57dc4077///mXP5+TksM8++zBq1Cjmz5/PZ599xkcffUSfPn2We9+Sv9uYMWPKdh4s7/slSZLKwzPNtdzcuXM5/fTT+f7772nQoAFbbLEFhYWFS72mRYsW/POf/6R///7EGNljjz3Yd999AVh33XWZOnUqvXr1YoMNNuC+++4DoKioiJNPPplLL72UBQsWMHjwYLp3715pdXfp0oVDDjmEzp0706BBA2644YayyRfHH388J510Ejk5OZx//vl88MEH1KtXj7Zt23LzzTenfL8kSdLqCqvzr+YzJScnJ1bnecKZVFxcTG5ubto+v3Hjxktte10bpbuHdYE9rDh7WDnsY8XZw4rLzs6mSZMmFBcXZ7qUjFv0t7S6vcjk32EIYWKMMWfZ4y7PkCRJklIwNGuVavtZZkmSpPIwNEuSJEkpGJrrgMaNG6d8zbXXXsu8efOqoJrlff/999x4442V8llXX301nTt3plu3buy8885Mnz697Ln69euXbbCyzz77rPD9Z555ZtlrOnbsSJMmTcqeGzRoEE2aNGGvvfaqlFolSVLNYWgWsGah+bfffquU767M0NyjRw9KSkp4++23Oeiggzj33HPLnltnnXXKxtqNGTNmhe+/5ppryl5z+umnc8ABB5Q9d84553DXXXdVSp2SJKlmMTTXIYuuRD3ooIPYcsstycvLI8bIddddx+zZs+nfvz/9+/cHYOzYsWy77bb07NmTgw8+uGxtc7t27bj44ovZYYcdeOCBB3jqqafo2bMn3bt3Z+eddwaSrbePPfZYevfuTY8ePXj00UcBGDFiBPvuuy+DBg2iU6dOZZuZnH/++XzyySdkZ2dzzjnnVOh37N+/P1lZWQD07duXmTNnrvFn3XvvvQwZMqTs8c4778x6661XofokSVLN5JzmOmbSpElMnTqVli1bsv322/PKK6/wxz/+kauvvpoXXniBpk2b8s0333DppZfy7LPPsu666/Kvf/2Lq6++mgsvvBCAtddem5dffpk5c+bQs2dPXnzxRdq3b8+3334LJNtaDxgwgDvuuIPvv/+ePn36sMsuuwDwxhtvMGXKFLKysujduzd77rknl112GVOmTFnpDoM77rgjP/3003LHr7zyyrLPXZHbb7+d3XffvezxL7/8Qk5ODg0aNOD8889nv/32W+l7p0+fzmeffcaAAQNStVSSJNUBhuY6pk+fPrRu3RpI5khOmzaNHXbYYanXjB8/nnfffZftt98egF9//ZVtt9227PlDDz207HU77bQT7du3B2CjjTYCkrPUY8aM4corrwSSsDpjxgwAdt11VzbeeGMADjjgAF5++eVVhleAl156abV/z7vvvpuSkhLGjRtXdmzGjBm0bNmSTz/9lAEDBrD11luz+eabr/D9o0aN4qCDDnJDFEmSBBia65xGjRqV3a9fvz4LFy5c7jUxRnbddVfuvffeFX7GuuuuW/a6EMIK3//ggw/SqVOnpY6//vrry71+Re9f1uqeaX722WcpKChg3LhxS/2+LVu2BGCzzTYjNzeXSZMmrTI033DDDSlrkyRJdYNrmgXAeuutVxZM+/btyyuvvMLHH38MwLx58/jwww+Xe8+2227LuHHj+OyzzwDKlmcMHDiQ4cOHs2i3yUmTJpW955lnnuHbb7/l//7v/3jkkUfYfvvtl/ruFXnppZfKLs5b8raiwDxp0iROPPFExowZwyabbFJ2/LvvvmP+/PkAfPPNN7zyyit07tx5hd/3wQcf8N133y11dl2SJNVthmYBkJ+fz+67707//v1p1qwZI0aMYMiQIXTr1o2+ffvy/vvvL/eeZs2aUVhYyAEHHED37t3Llm1ccMEFLFiwgG7dutG1a1cuuOCCsvfssMMOHHHEEWRnZ3PggQeSk5PDxhtvzPbbb0/Xrl0rfCHgOeecw9y5czn44IOXGi333nvvkZOTQ/fu3enfvz/nn39+WWi+4447lpqmce+99zJ48ODlzoLvuOOOHHzwwTz33HO0bt2ap59+ukK1SpKkmiMsOhtYneXk5MSSkpJMl1EtZXJv9tU1YsQISkpKuP766zNdylJqUg+rK3tYcfawctjHirOHFZednU2TJk0oLi7OdCkZt+hvaXV7kcm/wxDCxBhjzrLHPdMsSZIkpeCFgKoyRx99NEcffXSmy5AkSVptnmmWJEmSUjA0S5IkSSkYmiVJkirRtdde60WAtZChWZIkSUrB0CxJkiSlYGiWJElSpSsqKmL8+PGMGzeOdu3aUVRUlOmSKsTQLEmSpEpVVFREfn4+8+fPB2D69Onk5+fX6OBsaJYkSVKlGjp0KPPmzVvq2Lx58xg6dGiGKqo4Q7MkSZIq1YwZM1breE1gaJYkSVKlatOmzWodrwkMzZIkSapUBQUFZGVlLXUsKyuLgoKCDFVUcYZmSZIkVaq8vDwKCwtp1KgRAG3btqWwsJC8vLwMV7bmGmS6AEmSJNU+eXl53HrrrQC1YodEzzRLkiRJKRiaJUmSpBQyEppDCINCCB+EED4OIZyfiRokSZKk8qry0BxCqA/cAOwOdAaGhBA6V3UdkiRJUnll4kxzH+DjGOOnMcZfgVHAvhmoQ5IkSSqXTITmVsDnSzyeWXpMkiRJqpYyMXIurOBYXO5FIeQD+QDNmzevFaNK0mHu3Ln2poLsYcXZw4qzh5XDPlacPaw4e7jY999/D6z+yLnq2MNMhOaZwKZLPG4NzF72RTHGQqAQICcnJ+bm5lZJcTVNcXEx9qZi7GHF2cOKs4eVwz5WnD2sOHu4WJMmTQBWux/VsYeZWJ4xAegQQmgfQlgLGAyMyUAdkiRJUrlU+ZnmGOPCEMJpwNNAfeCOGOPUqq5DkiRJKq+MbKMdY3wCeCIT3y1JkiStLncElCRJklIwNEuSJEkpGJolSZKkFAzNkiRJUgqGZkmSJCkFQ7MkSZKUgqFZkiRJSsHQLEmSJKVgaJYkSZJSMDRLkiRJKRia64Avv/ySwYMHs/nmm9O5c2f22GMPPvzww0r/nmnTptG1a9eUr7nnnnvKHpeUlPDHP/6x0mupiCFDhtCtWzeuueYa3n//fbKzs+nRoweffPIJ22233Srfe+GFF/Lss8+u0fdOnjyZJ55Iz+7yI0eOpEOHDnTo0IGRI0eu8DUzZsygf//+9OjRg27dupXVMnnyZLbddlu6dOlCt27duO+++9JSoyRJ1VmDTBeg9Ioxsv/++3PUUUcxatQoIAlBX331FR07dkz5/t9++4369euv9PHqWhSaDzvsMABycnLIyclZ48+rbF9++SWvvvoq06dPB+Cyyy5j3333ZdiwYQC8+uqrq3z/xRdfvMbfPXnyZEpKSthjjz3W+DNW5Ntvv2XYsGGUlJQQQqBXr17ss88+bLjhhku97tJLL+WQQw7h5JNP5t1332WPPfZg2rRpZGVlceedd9KhQwdmz55Nr169GDhwIE2aNKnUOiVJqs4801zLvfDCCzRs2JCTTjqp7Fh2djY77rgjMUbOOeccunbtytZbb112BrG4uJj+/ftz2GGHsfXWWy/3+LfffuOcc86hd+/edOvWjVtuuWW57502bRo77rgjPXv2pGfPnmVh8/zzz+ell14iOzuba665huLiYvbaay8gCXf77bcf3bp1o2/fvrz99tsAXHTRRRx77LHk5uay2Wabcd11163wd33qqafo2bMn3bt3Z+edd17lZ/78888ce+yx9O7dmxNOOIFHH30UgN12242vv/6a7Oxshg0bxrXXXsttt91G//79AWjcuHHZ911++eVsvfXWdO/enfPPPx+Ao48+mtGjRwMwceJE+vXrVxYyv/jiCwByc3M577zz6NOnDx07duSll17i119/5cILL+S+++4jOzu7Us/mPv300+y6665stNFGbLjhhuy666489dRTy70uhMCPP/4IwA8//EDLli0B6NixIx06dACgZcuWbLLJJsyZM6fS6pMkqSbwTHMtN2XKFHr16rXC5x566CEmT57MW2+9xTfffEPv3r3ZaaedAHjjjTeYMmUK7du3p7i4eKnHhYWFbLDBBkyYMIH58+ez/fbbs9tuuxFCKPvsTTbZhGeeeYa1116bjz76iCFDhlBSUsJll13GlVdeyWOPPQYkAX2Rv//97/To0YNHHnmE559/niOPPJLJkycD8P777/PCCy/w008/0alTJ04++WQaNmxY9t45c+Zwwgkn8OKLL9K+fXu+/fbbVX5mQUEBAwYM4I477uCxxx7jrLPOYpdddmHMmDHstddeZd8bY6Rx48b8+c9/Xqp3Tz75JI888givv/46WVlZZd+3yIIFCzj99NN59NFHadasGffddx9Dhw7ljjvuAGDhwoW88cYbPPHEEwwbNoxnn32Wiy++mJKSEq6//vrl/rP64IMPOPTQQ1f4n2NxcfEqz/rOmjWLTTfdtOxx69atmTVr1nKvu+iii9htt90YPnw4P//88wqXmbzxxhv8+uuvbL755iv9PkmSaiNDcx328ssvM2TIEOrXr0/z5s3p168fEyZMYP3116dPnz60b9++7LVLPh47dixvv/122RnVH374gY8++mip5R4LFizgtNNOY/LkydSvX79ca6hffvllHnzwQQAGDBjA//73P3744QcA9txzTxo1akSjRo3YZJNN+Oqrr2jdunXZe8ePH89OO+1UVuNGG220ys8cO3YsY8aM4corr2Tu3Ln8+uuvzJgxg3XWWadcvXv22Wc55phjyMrKWur7Fvnggw+YMmUKu+66K5Asa2nRokXZ8wcccAAAvXr1Ytq0aSm/r1OnTmVBfnXFGJc7tuT/wVnk3nvv5eijj+bss8/mtdde44gjjmDKlCnUq5f8C6kvvviCI444gpEjR5YdkySprjA013JdunQpC7fLWlGYWmTddddd6eMYI8OHD2fgwIFLvWbJ8HfNNdfQvHlz3nrrLX7//XfWXnvtlLWuKtw1atSo7Fj9+vVZuHDhcu9dURBc2WfGGHnwwQfp1KkTxcXF5ObmLvc7pKp1Rd+35PNdunThtddeW+Hzi36fFf0uK7I6Z5pff/11TjzxRCBZY926deulzujPnDmz7Pdd0u233162bGPbbbfll19+4ZtvvmGTTTbhxx9/ZM899+TSSy+lb9++KeuVJKm28XRRLTdgwADmz5/PrbfeWnZswoQJjBs3jp122on77ruP3377jTlz5vDiiy/Sp0+flJ85cOBAbrrpJhYsWADAhx9+yM8//7zUa3744QdatGhBvXr1uOuuu/jtt98AWG+99fjpp59W+Lk77bQTRUVFQBIEmzZtyvrrr1+u33Pbbbdl3LhxfPbZZwBlyyVW9pkDBw5k+PDhZaF60qRJ5fqeRXbbbTfuuOMO5s2bt9T3LdKpUyfmzJlTFpoXLFjA1KlTV/mZq+rNojPNK7otuzRjm222KXtun332YeDAgYwdO5bvvvuO7777jrFjxy73f3gA2rRpw3PPPQfAe++9xy+//EKzZs349ddf2X///TnyyCM5+OCDy9UfSZJqG0NzLRdC4OGHH+aZZ55h8803p0uXLlx00UW0bNmS/fffn27dutG9e3cGDBjA5Zdfzh/+8IeUn3n88cfTuXNnevbsSdeuXTnxxBOXO1t6yimnMHLkSPr27cuHH35Ydqa6W7duNGjQgO7du3PNNdcs9Z6LLrqIkpISunXrxvnnn7/S0Wgr0qxZMwoLCznggAPo3r172VnZlX3mBRdcwIIFC+jWrRvHHHMMF1xwQbm/C2DQoEHss88+5OTkkJ2dzZVXXrnU82uttRajR4/mvPPOo3v37mRnZ6ecvNG/f3/efffdSr8QcKONNuKCCy6gd+/e9O7dmwsvvLBsOcmFF17ImDFjALjqqqu49dZb6d69O0OGDGHEiBGEELj//vt58cUXGTFiBNnZ2WRnZ6/xUhFJkmqqsKp/RV9d5OTkxJKSkkyXUS0tubRAa8YeVpw9rDh7WDnsY8XZw4qzh4st6sOSywTLI5M9DCFMjDEuNw/XM82SJElSCoZmSZIkKQVDsyRJkpSCoVmSJElKwTnNkiRJSovVvQCwOvNMsyRJkpSCoVmSJElKwdAsSZIkpWBoliRJklIwNEuSJEkpGJolSZKkFAzNkiRJUgqGZkmSJCkFQ7MkSZKUgqFZkiRJSsHQLEmSJKVgaJYkSZJSMDRLkiRJKRiaJUmSpBQMzZIkSVIKhmZJkiQpBUOzJEmSlIKhWZIkSUrB0CxJkiSlEGKMma4hpRDCHGB6puuoppoC32S6iBrOHlacPaw4e1g57GPF2cOKs4cVl8keto0xNlv2YI0IzVq5EEJJjDEn03XUZPaw4uxhxdnDymEfK84eVpw9rLjq2EOXZ0iSJEkpGJolSZKkFAzNNV9hpguoBexhxdnDirOHlcM+Vpw9rDh7WHHVroeuaZYkSZJS8EyzJEmSlIKhuYYKIQwKIXwQQvg4hHB+puupKUIId4QQvg4hTFni2EYhhGdCCB+V/twwkzVWdyGETUMIL4QQ3gshTA0h/Kn0uH0spxDC2iGEN0IIb5X2cFjpcXu4mkII9UMIk0IIj5U+toerIYQwLYTwTghhcgihpPSYPVwNIYQmIYTRIYT3S/+5uK09XD0hhE6lf4OLbj+GEM6obn00NNdAIYT6wA3A7kBnYEgIoXNmq6oxRgCDljl2PvBcjLED8FzpY63cQuDsGONWQF/g1NK/P/tYfvOBATHG7kA2MCiE0Bd7uCb+BLy3xGN7uPr6xxizlxjvZQ9Xz7+Bp2KMWwLdSf4e7eFqiDF+UPo3mA30AuYBD1PN+mhorpn6AB/HGD+NMf4KjAL2zXBNNUKM8UXg22UO7wuMLL0/EtivKmuqaWKMX8QY3yy9/xPJ/0C0wj6WW0zMLX3YsPQWsYerJYTQGtgTuG2Jw/aw4uxhOYUQ1gd2Am4HiDH+GmP8HntYETsDn8QYp1PN+mhorplaAZ8v8Xhm6TGtmeYxxi8gCYTAJhmup8YIIbQDegCvYx9XS+mygsnA18AzMUZ7uPquBc4Ffl/imD1cPREYG0KYGELILz1mD8tvM2AO8J/SZUK3hRDWxR5WxGDg3tL71aqPhuaaKazgmGNQVKVCCI2BB4EzYow/ZrqemibG+Fvpv4psDfQJIXTNcEk1SghhL+DrGOPETNdSw20fY+xJstzv1BDCTpkuqIZpAPQEboox9gB+xqUYayyEsBawD/BApmtZEUNzzTQT2HSJx62B2RmqpTb4KoTQAqD059cZrqfaCyE0JAnMRTHGh0oP28c1UPqvcotJ1trbw/LbHtgnhDCNZInagBDC3djD1RJjnF3682uSNaR9sIerYyYws/TfFAGMJgnR9nDN7A68GWP8qvRxteqjoblmmgB0CCG0L/1/ZYOBMRmuqSYbAxxVev8o4NEM1lLthRACyfq992KMVy/xlH0spxBCsxBCk9L76wC7AO9jD8stxviXGGPrGGM7kn8GPh9jPBx7WG4hhHVDCOstug/sBkzBHpZbjPFL4PMQQqfSQzsD72IP19QQFi/NgGrWRzc3qaFCCHuQrOerD9wRYyzIbEU1QwjhXiAXaAp8BfwdeAS4H2gDzAAOjjEue7GgSoUQdgBeAt5h8VrSv5Ksa7aP5RBC6EZyUUt9kpMX98cYLw4hbIw9XG0hhFzgzzHGvexh+YUQNiM5uwzJMoN7YowF9nD1hBCySS5GXQv4FDiG0v9eYw/LLYSQRXK91mYxxh9Kj1Wrv0VDsyRJkpSCyzMkSZKkFAzNkiRJUgqGZkmSJCkFQ7MkSZKUgqFZkiRJSsHQLEmSJKVgaJakNAkhzF2N1+aGELZb4vFJIYQjS+8fHUJouQbfPy2E0HQ13zO6dH7voppKQgiXL/F8cQihZInHOSGE4tL7W4cQRqxunZJUExiaJal6yAXKQnOM8eYY452lD48GVjs0r64QQhegfozx09JDJwM7AvVDCFsu8dJNQgi7L/v+GOM7QOsQQpt01ypJVc3QLElVKISwdwjh9RDCpBDCsyGE5iGEdsBJwJkhhMkhhB1DCBeFEP4cQjgIyAGKSp9bZ8kzyMuc6d04hDC29LNvAcIS33t4COGN0s+4JYRQfwXl5bH0NrX1gEiy82NY4vgVwN9W8iv+l2Rba0mqVQzNklS1Xgb6xhh7AKOAc2OM04CbgWtijNkxxpcWvTjGOBooAfJKn/u/VXz234GXSz97DMnWs4QQtgIOBbaPMWYDv5EE5GVtD0xc4vFtwKtAvRjje0scfw2YH0Lov4LPKCE5Oy1JtUqDTBcgSXVMa+C+EEILYC3gs0r87J2AAwBijI+HEL4rPb4z0AuYEEIAWAf4egXvbwHMWfQgxvg08PRKvutSkrPN5y1z/GuqYCmJJFU1zzRLUtUaDlwfY9waOBFYew0+YyGL//m97PvjCl4fgJGlZ6qzY4ydYowXreB1/1feemKMz5e+tu8yT61d+jmSVKsYmiWpam0AzCq9f9QSx38C1lvJe5Z9bhrJmWOAA5c4/iKlyy5KL9TbsPT4c8BBIYRNSp/bKITQdgXf8x6wRbl+i0QBcO4yxzoCU1bjMySpRjA0S1L6ZIUQZi5xOwu4CHgghPAS8M0Sr/0vsP+iCwGX+ZwRwM2LLgQEhgH/Lv2M35Z43TBgpxDCm8BuwAyAGOO7JEspxoYQ3gaeIVmKsazHSaZ4lEuM8QmWWM5Rqn/p50hSrRJiXNG/yZMk1TWlgfwFkgsGf0v1+hW8vxEwDtghxriwsuuTpEwyNEuSyoQQBgLvxRhnrMF7OwCtYozFlV6YJGWYoVmSJElKwTXNkiRJUgqGZkmSJCkFQ7MkSZKUgqFZkiRJSsHQLEmSJKXw/0jw/IHAFZiwAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Create figure and axes\n", "fig, ax = plt.subplots(1, 1, figsize=(12,8))\n", "\n", "# Plot average temperatures as a function of latitude (with errorbars)\n", "ax.errorbar(data[\"Latitude\"], data[\"Average temperature\"], fmt='ko', yerr=data[\"Deviation\"])\n", "\n", "# Calculate regression line points\n", "# NOTE: Using the unweighted regression line here!\n", "x1 = -5.0\n", "x2 = 70.0\n", "y1 = a + b * x1\n", "y2 = a + b * x2\n", "\n", "# Plot unweighted regression line\n", "ax.plot([x1, x2], [y1, y2], 'r-')\n", "\n", "# Add axis labels\n", "ax.set_xlabel(\"Latitude (°N)\")\n", "ax.set_ylabel(\"Temperature (°C)\")\n", "\n", "# Add text with statistical values\n", "ax.text(1, 2, f\"Line values:\\nSlope = {b:.2f}\\nIntercept = {a:.2f}\")\n", "ax.text(1, -2, f\"Correlation coefficient = {r:.2f}\")\n", "\n", "# Enable gridlines\n", "ax.grid()" ] }, { "cell_type": "markdown", "id": "f147b6b3-5fe7-4d90-910c-4336d063b6af", "metadata": {}, "source": [ "## Goodness-of-fit tests\n", "\n", "Goodness-of-fit tests provide a way to quantify how well predictions fit a given set of observations that have known uncertainties. One common measure is the *weighted sum of the squared errors*, also often refered to simply as chi squared:\n", "\n", "$$\n", "\\large\n", "\\chi^{2} = \\sum \\frac{(O_{i} - E_{i})^{2}}{\\sigma_{i}^{2}}\n", "$$\n", "\n", "where $O_{i}$ is the $i$th observed value, $E_{i}$ is the $i$th expected value, and $\\sigma_{i}$ is the $i$th standard deviation. For a given observation $i$, this equation will be less than 1 if the prediction falls within the uncertainties of the measurement or greater than 1 otherwise. This is fine, but it can then be confusing when the values are summed. Another option involves dividing by the number of values used in the calculation, as shown below.\n", "\n", "$$\n", "\\large\n", "\\chi^{2} = \\frac{1}{N} \\sum \\frac{(O_{i} - E_{i})^{2}}{\\sigma_{i}^{2}}\n", "$$\n", "\n", "where $N$ is the number of values in the summation. This results in a $\\chi^{2}$ value that should be less than 1 if the predictions are within the uncertainties *on average*, and greater than 1 otherwise." ] }, { "cell_type": "code", "execution_count": null, "id": "e6e9fd35-d1aa-4634-b514-af6826b473fe", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.10" } }, "nbformat": 4, "nbformat_minor": 5 }