{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Regression" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For many processes, the exact functional relation between input variables, i.e. the \"features\", and output variables, i.e. the \"targets\", is not known. In such cases, assuming a linear or low-order polynomial relation between the features and the targets may be a viable approach. The coefficients of the polynomial may then be learned from data. This approach is known as linear or polynomial regression." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Imports" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We need to import the following packages, classes, and functions." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# for handling data:\n", "import numpy as np\n", "import pandas as pd\n", "\n", "# for plotting:\n", "import matplotlib.pyplot as plt\n", "\n", "import halerium.core as hal\n", "\n", "# for graphs:\n", "from halerium.core import Graph, Entity, Variable, StaticVariable\n", "from halerium.core.regression import linear_regression, polynomial_regression, connect_via_regression\n", "\n", "# for models:\n", "from halerium.core import DataLinker, get_data_linker\n", "from halerium.core.model import MAPModel, ForwardModel, Trainer\n", "from halerium.core.model import get_posterior_model\n", "\n", "# for predictions:\n", "from halerium import Predictor\n", "\n", "# for analysing graphs:\n", "from halerium.core.utilities.print import print_child_tree" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Example data" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Consider a process with one feature \"x\" and one target \"y\". \n", "\n", "We first generate some data for this process:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEGCAYAAACKB4k+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAhxUlEQVR4nO3df5RcdZnn8fdT1dXVIY0QkgD5pUES1k1ykjjTgk5cjoCujGKiE2AYlWF3WZk/yKgzLAnoKuO67EAQXfcMOiciI+ygbEycSQZZFflxHJATaLATSMAlZ0DSSSRNDD86k1R3Vz37R92qVFXf6q5O+tatH5/XOW1Xf+ve6sc+oZ76/nq+5u6IiIgAJOIOQEREGoeSgoiIFCkpiIhIkZKCiIgUKSmIiEhRR9wBnIgZM2b4/Pnz4w5DRKSpPP3006+5+8yw55o6KcyfP5/e3t64wxARaSpm9ptqz2n4SEREipQURESkSElBRESKlBRERKRISUFERIqUFEREmszBwQzb97zOwcHMpL92Uy9JFRFpN1v69rJu8w5SiQTDuRzrVy9l5fI5k/b66imIiDSJg4MZ1m3ewdHhHG9lRjg6nGPt5h2T2mNQUhARaRL9h46QSpS/bacSCfoPHZm036GkICLSJOZOm8JwLlfWNpzLMXfalEn7HUoKIiJNYnp3mvWrl9KVSnByuoOuVIL1q5cyvTs9ab9DE80iIk1k5fI5rFgwg/5DR5g7bcqkJgRQUhARaTrTu9OTngwKNHwkIiJFkSUFM+sysyfNbLuZ7TSzrwTtf2Vme82sL/j6SMk9N5rZbjP7tZl9OKrYREQkXJTDRxngQncfNLMU8JiZ/d/guW+4+9dKLzazRcAVwGJgNvBzMzvH3bMRxigiEruDg5nI5ggmKrKk4O4ODAY/poIvH+OWVcB97p4BXjKz3cC5wBNRxSgiEreodyhPVKRzCmaWNLM+4ADwoLtvC55aY2Y7zOwuM5sWtM0B9pTc3h+0Vb7mNWbWa2a9AwMDUYYvIhKpeuxQnqhIk4K7Z919OTAXONfMlgDfBs4GlgP7gduDyy3sJUJec4O797h7z8yZoUeMiog0hXrsUJ6ouqw+cvfXgUeBi9391SBZ5IDvkB8ignzPYF7JbXOBffWIT0QkDvXYoTxRUa4+mmlmpwaPpwAfBF4ws1kll30CeC54vBW4wszSZnYWsBB4Mqr4RETiNr07zeU9c8vaLu+ZG+tkc5Srj2YBd5tZknzy2eju95vZ/zaz5eSHhl4G/gzA3Xea2UZgFzACXKuVRyLSyg4OZtjY21/WtrG3n89ddE5siSHK1Uc7gHeHtF85xj03AzdHFZOISCMpzCkc5dgQUmFOIa6koB3NIiIxaas5BRERGVs9qp5OlAriiYjEKOqqpxOlpCAiErMoq55OlIaPRESkSElBRESKlBRERKRISUFERIqUFEREpEhJQUREipQURESkSElBRESKlBRERKRISUFERIqUFEREpEhJQUREipQURESkKMozmrvM7Ekz225mO83sK0H7aWb2oJm9GHyfVnLPjWa228x+bWYfjio2EREJF2VPIQNc6O7LgOXAxWb2XuAG4CF3Xwg8FPyMmS0CrgAWAxcD3wrOdxYRkTqJLCl43mDwYyr4cmAVcHfQfjfw8eDxKuA+d8+4+0vAbuDcqOITEZHRIp1TMLOkmfUBB4AH3X0bcIa77wcIvp8eXD4H2FNye3/QVvma15hZr5n1DgwMRBm+iEjbiTQpuHvW3ZcDc4FzzWzJGJdb2EuEvOYGd+9x956ZM2dOUqQiIgJ1Wn3k7q8Dj5KfK3jVzGYBBN8PBJf1A/NKbpsL7KtHfCIitTo4mGH7ntc5OJiJO5RIRLn6aKaZnRo8ngJ8EHgB2ApcFVx2FbAleLwVuMLM0mZ2FrAQeDKq+EREJmpL315W3Pown75zGytufZitfXvjDmnSdUT42rOAu4MVRAlgo7vfb2ZPABvN7GrgFeAyAHffaWYbgV3ACHCtu2cjjE9E2sjBwQz9h44wd9oUpnenq7aNdf+6zTs4OpzjKDkA1m7ewYoFM8a9t5lElhTcfQfw7pD2g8BFVe65Gbg5qphEpD1t6dvLus07SCUSDOdyrF+9FIdRbSuXj1rbUtR/6AipRKKYEABSiQT9h44oKYiINIuwT/jXb9oBOJkRr/lT/9xpUxjO5crahnM55k6bEmn89aYyFyLS0gqf8EslE0bSytsKn/qrmd6dZv3qpXSlEpyc7qArlWD96qUt1UsA9RREpMWFfcLP5pzKFe+1fOpfuXwOKxbMqHkeohmppyAiLS3sE/5tly7ltkuXHden/undaZbNO7UlEwKopyAibaDaJ/xW/9R/PJQURKQtTO9Oj3rjD2trdxo+EpGW1Oo7j6OinoKItJywfQlj7UGQY9RTEJGWUrov4a3MCEeHc6zdvEM9hhopKYhISwnblzDeHgQ5RklBRFpK2L6EoWyWN44MqbdQAyUFEWkplfsSOhKQc7j23l+1bGXTyaSkICItoXS10crlc3h83YXc8anfI5lIMJx1zS/USKuPRKTpVVttdMqUFJ3JBJmR1q5sOpnUUxCRpjbWaqN2qWw6mZQURKSpjbXaqF0qm06myIaPzGwecA9wJpADNrj7N83sr4DPAAPBpV9w9weCe24ErgaywGfd/adRxScirWFqZ5LMSPkhjaW9gXaobDqZopxTGAGuc/dnzOxk4GkzezB47hvu/rXSi81sEXAFsBiYDfzczM7RkZwiUk1hLiGRMMg66aRhCRvVG1CNo9pFeRznfmB/8PgtM3seGGuf+SrgPnfPAC+Z2W7gXOCJqGIUkeZVOpdQ4Gb8eM37WXDGyTFG1tzqMqdgZvPJn9e8LWhaY2Y7zOwuM5sWtM0B9pTc1k9IEjGza8ys18x6BwYGKp8WkTYROpeQNPpUBO+ERJ4UzKwb2Ax83t3fBL4NnA0sJ9+TuL1wacjtPqrBfYO797h7z8yZM6MJWkQaXtjKosOZLDdt3alNaicg0qRgZinyCeFed/8RgLu/6u5Zd88B3yE/RAT5nsG8ktvnAvuijE9EmkdlKezSlUVT08nidYeHstqkdgKiXH1kwHeB59396yXts4L5BoBPAM8Fj7cC3zezr5OfaF4IPBlVfCLSPKptTiusLHrkhQPctHUnh4eOrUvRJrXjE2VPYQVwJXChmfUFXx8B1pvZs2a2A7gA+AsAd98JbAR2AT8BrtXKIxEZrxT29O40F7zrdLJePtqsTWrHJ8rVR48RPk/wwBj33AzcHFVMItJ8ChPKR6leqqIwlLS2ojehXsLEqfaRiDS0WktVaJPa5FCZCxFpaBMpVTG9O82yeacqIZwA9RREpO4ODmYm9IlevYD6UVIQkbqqtpJoPCpVUR8aPhKRuhlvJZHET0lBROqm/9ARPFe+dNRzTv+hIzFFJJWUFESkbqZ2Jslky5NCJutM7UxWuUPqTUlBROrm8FCWrlT5205XKlG2E1nipaQgInVTbYexdh43DiUFEakbHY/Z+LQkVUTqSnsOGpuSgojUnfYcNC4lBRGZNAcHM+zc9ybgLJ59it74m5CSgohMii19e7luYx8jQe26VNK4/bJlNe1WlsahiWYRmbDKU9AODmZYu2l7MSEADGed6zdpt3KzUU9BRCYkrHbRO6ZPJWkJoHy/QTJhOv2syainICI1q1a7aGpnkqznRl2fzbn2IDSZyJKCmc0zs0fM7Hkz22lmnwvaTzOzB83sxeD7tJJ7bjSz3Wb2azP7cFSxicjxKZyCVippxr43jrLmgoUkS85aTCWN2y7VHoRmE+Xw0Qhwnbs/Y2YnA0+b2YPAfwAecvdbzOwG4AZgnZktAq4AFgOzgZ+b2Tk6p1mkcYSdgnZ4KMvVdz9FV0eSjmSCa1bM531nT9fqoyYVWU/B3fe7+zPB47eA54E5wCrg7uCyu4GPB49XAfe5e8bdXwJ2A+dGFZ+ITNz07jRf+uiiUe3DWeetzAiZkRx/98uXlRCaWF3mFMxsPvBuYBtwhrvvh3ziAE4PLpsD7Cm5rT9oq3yta8ys18x6BwYGIo1bRModHMyQ7kiMWdU0lUioFHYTi3z1kZl1A5uBz7v7m2ZW9dKQNh/V4L4B2ADQ09Mz6nkRiUZh1VHSbMyqpsO5nCaXm1ikPQUzS5FPCPe6+4+C5lfNbFbw/CzgQNDeD8wruX0usC/K+ETkmMq9B5XPFVYdlSaEqekkHYn8pLIK3LWGyHoKlu8SfBd43t2/XvLUVuAq4Jbg+5aS9u+b2dfJTzQvBJ6MKj4ROaZy78GXPrqIJXNOKRasK6w6OsqxSeapnUm+8rHFXPCu/AiwCty1hiiHj1YAVwLPmllf0PYF8slgo5ldDbwCXAbg7jvNbCOwi/zKpWu18kgkWvlaRW+wdtN2MiNefNP/4j8+R3c6yUjOWb96KSsWzBi16mg45yyfd2oxCSgZtAZzb95h+Z6eHu/t7Y07DJGmc3Aww73bXuGOR3aTTBj/OsYcQVcqwePrLuQnO3/LV/5pFwZkRnKkk4YljPWrl6q+UZMxs6fdvSfsOZW5EGkhBwcz4w7jbOnbW+wZ1CKVSHDvtlf41qO7SSXg8FC+x5DJOmSdtZt3sGLBDPUUWoSSgkiLCKtJVPkJvjBhXC0hdHUkODpSPkw0lM1xxyMvVr2nsARVSaE1qPaRSAuoVpOociVRWJmKgnRHgg1/2sPNn1hSdlzmmgsW0Jmsvi9BS1Bbi3oKIi0gbHVQ2Cf4sDIVkE8It126lPPPmQnAxYvPLA5DAdzx6O5R9+SL4LmWoLYYJQWRFhD2Zh/2CX56d5r1q5dy/aYdJBPGSDbHn1+4kE+e9/ayN/bK4zLXr17K2tIlq5csYsnsU7QEtQUpKYi0gMKb/dqKOYWwN2wv/K8bZvCO6SeN+8a+cvkcViyYob0IbWDcpGBma8jvSD5Uh3hE5DiN98ZduSehcCBOrauHKnsP0ppq6SmcCTxlZs8AdwE/9Wbe3CDSwqq9cRdWJiXMRq0i0uohKTXu6iN3/6/kS058l/xZCC+a2f8ws7Mjjk1EJkHpyqSwTWpaPSSlalqSGvQMfht8jQDTgE1mtj7C2ERkAqoVtNu5700SIUWIT0olVcBORqllTuGz5AvXvQbcCVzv7sNmlgBeBNZGG6KIjCds49qKBTOKpSwyFRvS0h0J/vbK32fx7LcpIUiZWuYUZgB/5O6/KW1095yZXRJNWCJSq9LhocI+het+uB3DCStp1JmkbE+CSKlxk4K7f3mM556f3HBEZKLCNq4NZ6uvBfHQ86xE8lTmQqTJVdulXM1wUMQu7DAdESUFkSZTOaFc2LhWqFeU7kjQMc5/2TpHWarRjmaRJlKtEmrlxrXHd7/Gf/nhdoaqDCNpGapUE1lPwczuMrMDZvZcSdtfmdleM+sLvj5S8tyNZrbbzH5tZh+OKi6RZjVeJdTp3WmWBSehrVw+hwc+++/oDOkypDtMy1ClqiiHj74HXBzS/g13Xx58PQBgZouAK4DFwT3fMrPqtXpF2lBY2etkwti5743Q/QkLzjiZr11aPqx03YfO4Zc3XKST0qSqyIaP3P0XZja/xstXAfe5ewZ4ycx2A+cCT0QVn0izCZtQPpzJ8p++9xRTUh2hB+uokJ1MVBwTzWvMbEcwvDQtaJsD7Cm5pj9oE2lLYbuTp3en+dIli0ZdO5JjzIN1SoeVRMZT76TwbeBsYDmwH7g9aA9bOB06Q2Zm15hZr5n1DgwMRBKkSJy29O1lxa0P8+k7t7Hi1ofZ2re3+NyS2acwtbP6yGoqkag6nCRSi7quPnL3VwuPzew7wP3Bj/3AvJJL5wL7qrzGBmADQE9Pj6q1SksJ251cWtp67rQpZMcoUnx0JMtn7umlM5msek6zyFjq2lMws1klP34CKKxM2gpcYWZpMzuLfFXWJ+sZm0icCsNFO/e9OWoyuXRPQeWehI4EpJJWnEh2dzIjPuZwkshYIuspmNkPgA8AM8ysH7gJ+ICZLSc/NPQy8GcA7r7TzDYCu8hXYb3W3UOqtoi0ntK9B0PZHNlxjtWsnDyG/MqkN44Mc+29z/BWZqR4rc5KkImKcvXRn4Q0f3eM628Gbo4qHpFGFDZclEoa6Q7KhoAq39QrD9OZ3p3m4GCmpnOaRcaiHc0iETg4mKH/0BGmdiY5PJStuhw0rJhdV0eSv/6jJRwdzrF83qksOOPkmn7nRM5pFqlGSUFkkhWGgzznZLJOVyo/R1A66VuaNCo/3f/r0AjX/XAHHQljOJvjpo8t5lPvfUdNv1v7EuREKSmITKLS4aCCwuPCKqLHdr9WVr9o5bLZbH66n0KZoqxDdiRHYXr4i//4HBh86rzaEkO1c5pFaqEqqSKTKKwURUFhD0Fl/aKNvccSQjVf+addWkUkdaGkIDKJxjrbIN9ueG7i22tSSVOpa6kLJQWRSVS6j6CyQOnlPXOZfUoXmfG6BSGyOdcqIqkLJQWRSbZy+RzuX/N+ElZevWVjbz/73jhanHiuJt2R4PKeuXR2JJiaTtKVSmgVkdSNJppFjkNh9VC1FT6Hh7KkO5IMZcs3klUp6cXUziQjuRxrLljIJ897O9O706y7+F1aRSR1p6QgMkHVTj8rODiYYc/vDnN0ZKTsvuFcjsWzTxm1l+BLlyxiyexTRr35axWRxEFJQWQCxitYt6VvL9dt7GOkZK45nTQscey0M+0lkEampCBtabzhn2p27nuTREWl99KCdWs3bS9LCACO8eM17y/bmaxegDQqJQVpO+MN/4x139pNO8iMhNcX6j90hKQlgPJajh1J4/CQ6jtKc9DqI2krpcM/EykvvfvVt7j+h9tHJQSAlctmFUtWZH3081pOKs1EPQVpK2EF6MYrL72lby/Xb9rBUJX9BRt79/LAs79lJOf88Xvm8f1trxSHkFJJ47ZLtZxUmoeSgrSVsB3HY5WXLvQshkJ6CKUGM/nhoY29/fzkc+ez742jgLN49ilKCNJUNHwkbaXy5LLxNobt3PfGqE1oY0klEhweynL+OTM5/5zTlRCk6ainIG2n1iWh+Ynl7WRGyoeNOhJGwiCVTIyaQNahNtLsojyO8y7gEuCAuy8J2k4D/g8wn/xxnJe7+6HguRuBq8kv3fisu/80qthECktCC2cjF5JD6TkH6zbvGJUQ0h0Jbrt0aTGpPLfvDb56/y4daiMtw9wnXpyrphc2Ox8YBO4pSQrrgd+5+y1mdgMwzd3Xmdki4AfAucBs4OfAOeOd09zT0+O9vb2RxC/NrZZ9CIWeQNISZD3HH79nHht7+0klEmRGsiQSVnYuwkmpJH975e9z/jkzJ/y7RBqJmT3t7j1hz0V5RvMvzGx+RfMq4APB47uBR4F1Qft97p4BXjKz3eQTxBNRxSetq5Z9CAcHMyU7j/OfPe554hWAYyuTKlYb5XAWz37bqN+njWjSSuo90XyGu+8HCL6fHrTPAfaUXNcftI1iZteYWa+Z9Q4MDEQarDSfWvch7Nz35qidx5U6EtDZUduEtEiraJSJ5rDlHaHjWu6+AdgA+eGjKIOSxhY2bNN/6MioQ2w85yH7EMb/pzOSg02fOZdUR1JDQ9I26p0UXjWzWe6+38xmAQeC9n5gXsl1c4F9dY5Nmki1IaKpnclRh9hkss7UzmRZ2+LZp5BKGsNjHHiTThqpjiTL5p0axf8FkYZU7+GjrcBVweOrgC0l7VeYWdrMzgIWAk/WOTZpEmMNET3w3G9HXd+VGr10dHp3mtsvW0a6I8FJnUnSHaNPSrOEaXmptJ0ol6T+gPyk8gwz6wduAm4BNprZ1cArwGUA7r7TzDYCu4AR4NrxVh5J+worVZFMGDv3vckdj+wedb17eO2hyv0Kj+9+reycA80hSDuKcvXRn1R56qIq198M3BxVPNI65k6bwlC2/DPD4UyWe554ObRg3ZoLFlZ9cy9dOaRzDkQaZ6JZpGaP7X6NXMhUwM+fPzCqLd2R4A+XnFm2QW0sWl4q7U5JQZpKYT5hrAniUhe963Qu+ZvHQvcsaNOZyGhKChKL8d6Qqz0fNp9QTWcSHnrhAJmR0UdnPrb7teM6aEek1SkpSN2Nt+N4rOfD5hOq+eCiM/jn/3ewbJ4hlUiwc98bY56zLNLOVDpb6mq8HcfVnt/96lts3/M6kJ84rsVDzx9gKDv67AQwUonyf/ql5yyLtDP1FKSuxjv5LOx5d/jI//pn0h1JhnM5vvTRRaQ7EqErjUp1JpNcc/47uePR3WW9jsWz3zahg3ZE2omSgkQmbF5gvJPPwp4vvPkPZUcA+OqPd/Hljy3iq/fvIpkwDmfCh5OGczk+ed7b+cMlZ9K353WWzzuVBWecDMD61Uu1J0EkRGSls+tBpbMb11jzAlv79o56Qy6dU7h322/44j88V/W1T0538Pf/+TzmTpuSP9Ng7xt89ce78JyTyTpdqfzQ0PrVS3GoGodWH0m7iqV0trSv0nmBsIncapvECm/S86adxNTO5KjSFAVD2WzxvundaZbNO5WLl5xZPBzn8FC22PNYcevDVePQngSR0ZQUZNKNN28AozeJlfYshrI5srnq8wVhO5TD3uC373l93DhEpJxWH8mkG2/eoHAEZrUVR5mRHGZGRWFTIL9D+ZPnvX1S4hCR0ZQUZNJN706zfvVSulKjD6jZ0reXFbc+zKfv3MaKWx9ma9/eYs+iVFdHkjuveg/Xfegc0iUH3dx2ae0TwmPFISLhNNHcxqKeaC28fmGcf2pnkkv+5rGyc4/THcbtly3nLzf+itIphK5UgsfXXcj07vQJx6kJZZFymmiWUWo5x/hETe9Ol5WTyGRzWMWHkMyI85cb+8oSQkeCsk/0JzohrAllkdpp+KgN1XqO8WT/nqGR3KhT0QCGKtqSiQQrFsyY1FhEpDZKCm0obAw/ijIPoXMFqQSpkAnkUp1JlZwQiYuSQhuq16qcsN8D8I3LlzMlVf2fnlYIicQnlqRgZi+b2bNm1mdmvUHbaWb2oJm9GHyfFkds7eB4VuVULiM9kd/zvrNnELa8oTNpWiEkErNYVh+Z2ctAj7u/VtK2Hvidu99iZjcA09x93Vivo9VHJ6bWVTljTUrX8hph13x5y7Pc88QrxWtWLpvF1e9/p1YIidRBs6w+WgV8IHh8N/AoMGZSkBNTy6qcsUpW1HpQTeXvOTiYYWNvf9k1P9v1Kjd9bLESgkjM4ppTcOBnZva0mV0TtJ3h7vsBgu+nh91oZteYWa+Z9Q4MDNQp3NY21tBQtUnp0oNqJrqCqV4T3SIycXH1FFa4+z4zOx140MxeqPVGd98AbID88FFUAbaL8fYrVJuULhxUczx1hVR+QqRxxdJTcPd9wfcDwD8A5wKvmtksgOD7gThiaydj7Vco9B6A0MnixbPfxpHhkbLXOzI8UtMbu8pPiDSuuvcUzGwqkHD3t4LH/x74b8BW4CrgluD7lnrH1m6qVTO9d9srfKvitLLH111YNll8cDCDmUHJOqL8z7WpVj5bROIVR0/hDOAxM9sOPAn82N1/Qj4ZfMjMXgQ+FPwsEQobxhnK5rjjkRdH9R4Als07tfjm3X/oCF0d5bvQujqSE5oXKJyFoIQg0jjqnhTc/V/cfVnwtdjdbw7aD7r7Re6+MPj+u3rH1ogmuj9gIteHDeOsuWABncnyN/uwSWDNC4i0pkZakioVJlq07niK3FUO4wDc8ejusmvC3uwLCUXnHIu0FpXOblAHBzPFoyQLSstJn+j1YxnvDOXK36t5AZHm0iyb16RELUdansj1Y5nIJLDKUou0FiWFBjXRMfvJHuPXm71Ie1KV1AY10bX8WvsvIpNBcwoNbqJj9pVHYGqsX0QqaU6hiU10GKfyCMyojtoUkdak4aMWU6+jNkWkNSkptBhVIBWRE6Gk0GK001hEToSSQh0cz1GWx0urkETkRGiiOWLHU3riRKkCqYgcLyWFCI11lGXUb9TafCYix0PDRxHSpK+INBslhQhp0ldEmo2SQoQ06Ssizabh5hTM7GLgm0ASuNPdm/oENk36ikgzaaikYGZJ4A7yx3H2A0+Z2VZ33xVvZCdGk74i0iwabfjoXGB3cGTnEHAfsCrmmERE2kajJYU5wJ6Sn/uDtiIzu8bMes2sd2BgoK7BiYi0ukZLChbSVlbb2903uHuPu/fMnDmzTmGJiLSHRksK/cC8kp/nAvtiikVEpO00WlJ4ClhoZmeZWSdwBbA15phERNpGQ60+cvcRM1sD/JT8ktS73H1nzGGJiLSNhkoKAO7+APBA3HGIiLSjRhs+EhGRGLVtUqjnGQciIs2i4YaP6iGOMw5ERJpB2/UUdLC9iEh1bZcUdMaBiEh1bZcUdMaBiEh1bZcUdMaBiEh1bTnRrDMORETCtWVSAJ1xICISpu2Gj0REpDolBRERKVJSEBGRIiUFEREpUlIQEZEic/fxr2pQZjYA/CbmMGYAr8UcQ6PQ3+IY/S2O0d/imEb5W7zD3UPPM27qpNAIzKzX3XvijqMR6G9xjP4Wx+hvcUwz/C00fCQiIkVKCiIiUqSkcOI2xB1AA9Hf4hj9LY7R3+KYhv9baE5BRESK1FMQEZEiJQURESlSUjhOZjbPzB4xs+fNbKeZfS7umOJmZkkz+5WZ3R93LHEys1PNbJOZvRD8+3hf3DHFxcz+Ivjv4zkz+4GZdcUdU72Y2V1mdsDMnitpO83MHjSzF4Pv0+KMMYySwvEbAa5z938LvBe41swWxRxT3D4HPB93EA3gm8BP3P1dwDLa9G9iZnOAzwI97r4ESAJXxBtVXX0PuLii7QbgIXdfCDwU/NxQlBSOk7vvd/dngsdvkf8Pf068UcXHzOYCHwXujDuWOJnZ24Dzge8CuPuQu78ea1Dx6gCmmFkHcBKwL+Z46sbdfwH8rqJ5FXB38Phu4OP1jKkWSgqTwMzmA+8GtsUcSpz+J7AWyI1zXat7JzAA/F0wlHanmU2NO6g4uPte4GvAK8B+4A13/1m8UcXuDHffD/kPlsDpMcczipLCCTKzbmAz8Hl3fzPueOJgZpcAB9z96bhjaQAdwO8B33b3dwOHacAhgnoIxstXAWcBs4GpZvbpeKOS8SgpnAAzS5FPCPe6+4/ijidGK4CVZvYycB9woZn9fbwhxaYf6Hf3Qq9xE/kk0Y4+CLzk7gPuPgz8CPiDmGOK26tmNgsg+H4g5nhGUVI4TmZm5MeNn3f3r8cdT5zc/UZ3n+vu88lPJD7s7m35idDdfwvsMbN/EzRdBOyKMaQ4vQK818xOCv57uYg2nXQvsRW4Knh8FbAlxlhCdcQdQBNbAVwJPGtmfUHbF9z9gfhCkgbx58C9ZtYJ/AvwH2OOJxbuvs3MNgHPkF+t9yuaoMzDZDGzHwAfAGaYWT9wE3ALsNHMriafNC+LL8JwKnMhIiJFGj4SEZEiJQURESlSUhARkSIlBRERKVJSEBGRIiUFEREpUlIQEZEiJQWRSWRm7zGzHWbWZWZTg7MElsQdl0ittHlNZJKZ2X8HuoAp5Osg/XXMIYnUTElBZJIF5S2eAo4Cf+Du2ZhDEqmZho9EJt9pQDdwMvkeg0jTUE9BZJKZ2VbyJcTPAma5+5qYQxKpmaqkikwiM/tTYMTdv29mSeCXZnahuz8cd2witVBPQUREijSnICIiRUoKIiJSpKQgIiJFSgoiIlKkpCAiIkVKCiIiUqSkICIiRf8fIoa4GUkM4EIAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "n_data = 100\n", "\n", "x_data = np.random.normal(loc=1, scale=2, size=(n_data,)) + 5\n", "y_data = np.random.normal(loc=0, scale=9, size=(n_data,)) - 12 + 4 * x_data + 3 * x_data**2\n", "\n", "data = pd.DataFrame()\n", "data[\"x\"] = x_data\n", "data[\"y\"] = y_data\n", "\n", "data.plot.scatter(\"x\", \"y\");\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can extract simple statistical properties such as the mean and standard deviation:" ] }, { "cell_type": "code", "execution_count": 3, "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", "
xy
count100.000000100.000000
mean6.078618134.753586
std1.82545877.458972
min1.470327-18.618306
25%4.84391278.938936
50%6.003939119.147073
75%7.227908184.538179
max10.549351367.353488
\n", "
" ], "text/plain": [ " x y\n", "count 100.000000 100.000000\n", "mean 6.078618 134.753586\n", "std 1.825458 77.458972\n", "min 1.470327 -18.618306\n", "25% 4.843912 78.938936\n", "50% 6.003939 119.147073\n", "75% 7.227908 184.538179\n", "max 10.549351 367.353488" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "display(data.describe())\n", "\n", "x_data_mean = data['x'].mean()\n", "x_data_std = data['x'].std()\n", "\n", "y_data_mean = data['y'].mean()\n", "y_data_std = data['y'].std()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the following, we assume we do not know the exact functional relation that generated that data." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Linear regression model by hand" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Before we discuss how to quickly bulid a regression model using convenience functions, we create a linear regresssion model 'by hand'.\n", "\n", "To build the linear model, we create a `Graph` containing the variable \"x\" representing the feature and the variable \"y\" representing the target. We then add model parameters \"slope\" and \"intercept\" to our graph, and use them to connect the feature and target variable:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "graph = Graph(\"graph\")\n", "with graph:\n", " x = Variable(\"x\", shape=(), mean=x_data_mean, variance=x_data_std**2)\n", " y = Variable(\"y\", shape=(), variance=y_data_std**2)\n", " \n", " slope = StaticVariable(\"slope\", mean=0, variance=1e4)\n", " intercept = StaticVariable(\"intercept\", mean=0, variance=1e4)\n", " \n", " y.mean = slope * x + intercept \n", " \n", "hal.show(graph)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We now train a `MAPModel` with the data." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "inferred slope = 40.32779298524231\n", "inferred intercept = -109.72532491801661\n" ] } ], "source": [ "model = MAPModel(graph=graph, \n", " data={graph.x: data[\"x\"], graph.y: data[\"y\"]})\n", "model.solve()\n", "inferred_slope = model.get_means(graph.slope)\n", "inferred_intercept = model.get_means(graph.intercept)\n", "\n", "print(\"inferred slope =\",inferred_slope)\n", "print(\"inferred intercept =\",inferred_intercept)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We extract the posterior graph from the trained model and build a `ForwardModel` with it to compute predictions." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEGCAYAAABhMDI9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAtAklEQVR4nO3de3jcZZ338fc3yUxmmqRQ6cgirZk+2scWWOhZsOwDaIoosLC9RIkXWmyAIkdRUQ4eVgRkd12xq4+2LFPavVbDaVHA9UEJgidOPaRKoYAcJiVLlSnSkjaHmUnu54/fTDqTTNJJmslkMp/XdeWayW9+M3MzNL/v3Pf3vr+3OecQERFJqyh2A0REZGJRYBARkSwKDCIikkWBQUREsigwiIhIlqpiN+BgTZ8+3YXD4WI3Q0SkpGzevHmXcy6U67GSDwzhcJhNmzYVuxkiIiXFzNqGekxDSSIikkWBQUREsigwiIhIlpLPMeSSSCRob2+nu7u72E0pS4FAgBkzZuDz+YrdFBEZhUkZGNrb26mrqyMcDmNmxW5OWXHO8eabb9Le3s6sWbOK3RwRGYVJOZTU3d3NYYcdpqBQBGbGYYcdpt6aSIHFYrBxo3c71iZlYAAUFIpIn71IYTU3Q309LFvm3TY3j+3rT9rAICIyGcVi0NQEXV2wZ49329Q0tj0HBYYC2L17Nz/4wQ+K3QwRmYSiUfD7s4/5fN7xsaLAUABDBYbe3t4itEZEJpNwGOLx7GOJhHd8rCgwpIxlIueaa67h5ZdfZt68eSxevJhTTjmFT37yk/zt3/4t0WiUY445pv/cb3/72/zjP/4jAC+//DKnnXYaCxcu5O/+7u94/vnnD74xIjKphEIQiUAwCFOnereRiHd8rEzK6aoj1dzsjdH5/V4kjkSgsXH0r3fLLbewbds2tm7dymOPPcbpp5/Otm3bmDVrFtFh+nsXXXQRa9asYfbs2Tz11FNccskl/OpXvxp9Q0RkUmpshIYGb/goHB7boAAKDFmJnK4u71hTk/ehj9WHvWTJkgPO6d+7dy+PP/4455xzTv+xnp6esWmAiEw6odDYB4S0sg8M6UROOijA/kTOWH3oNTU1/ferqqro6+vr/z0937+vr49DDz2UrVu3js2bioiMUtnnGAqRyKmrq6OjoyPnY4cffjhvvPEGb775Jj09PfzsZz8DYOrUqcyaNYt77rkH8FYQ/+EPfxh9I0SkpBRywdpIlX1gKEQi57DDDmPp0qUcc8wxXH311VmP+Xw+vva1r/H+97+fM844gzlz5vQ/9qMf/YhIJMJxxx3H0Ucfzf333z/6RohIySj0grWRMudccVtwkBYtWuQGbtSzfft25s6dO6LXicUKl8gpR6P5fyBSjmIxLxhkDmcHg9DWVthrkZltds4tyvVY2ecY0gqZyBERGcp45DlHquyHkkREimk8FqyNlAKDiEgRhULeFPlMTU3FHcFQYBARKaJYzJvwkikSKe7sJAUGEZEiGo+ieCOlwCAiUkTKMZSR2tpaAF5//XU+9rGPFbk1Y+exxx7jjDPOKHYzRCaN8SiKN1Karlpg73rXu7j33ntH/fxkMklVVeH/N/X29lJZWVnw9xGRwQpdFG+k1GNIK9B69Mwy2+vXr2f58uWcdtppzJ49my996Us5n3PyySdz3XXXcdJJJ7F69Wo2b97MSSedxMKFC/nwhz/Mzp07Adi4cSPHHnssJ5xwAldffXVWOe/M958zZw4rVqzg2GOP5WMf+xidnZ0AhMNhbrjhBk488UTuuecefvnLX3LCCSewYMECzjnnHPbu3QvAQw89xJw5czjxxBO57777xvTzERFPKASLF48gKBSwhoYCA4zrevStW7dy11138cwzz3DXXXfx2muv5Txv9+7d/PrXv+aKK67g8ssv595772Xz5s2sXLmS66+/HoDPfOYzrFmzhieeeGLYb/svvPACF110EX/84x+ZOnVq1iZCgUCA3/3udzQ0NHDjjTfS0tLCli1bWLRoEd/5znfo7u7mwgsv5MEHH+S3v/0tf/7zn8f2AxGRkSvwNUuBYTw2UM3woQ99iEMOOYRAIMBRRx1FW1tbzvM+8YlPAN5Ffdu2bSxbtox58+Zx44030t7ezu7du+no6OADH/gAAJ/85CeHfM+ZM2eydOlSAM477zx+97vfDXqfJ598kueee46lS5cyb948NmzYQFtbG88//zyzZs1i9uzZmBnnnXfemHwOIjJK43DNUo5hnNejV1dX99+vrKwkmUzmPC9dqts5x9FHH80TTzyR9fhbb72V93ua2ZC/Z77PsmXLaB7wzWPr1q2Dni8iRTQO1yz1GCbiXLEM73vf+4jFYv2BIZFI8OyzzzJt2jTq6up48sknAbjzzjuHfI0dO3b0P7+5uZkTTzxx0DnHH388v//973nppZcA6Ozs5MUXX2TOnDm8+uqrvPzyy/3PF5EiGodrlgLDBJkrdsEFFzCwSiyA3+/n3nvv5ctf/jLHHXcc8+bN4/HHHwcgEolw0UUXccIJJ+Cc45BDDgG8KbIf/ehH+19j7ty5bNiwgWOPPZa//vWvfPaznx30PqFQiPXr19PY2Mixxx7L8ccfz/PPP08gEOC2227j9NNP58QTT6S+vr7/OZs2beKCCy4Y649CRIYzDtcsld1OK8G623v37u1fL3HLLbewc+dOVq9enXVONBrljDPOYNu2bePaNpXdFimwg7xmqex2Pkqw7vZ///d/861vfYtkMkl9fT3r168vdpNEZLwU8JqlwFDCPvGJT/TPKhpKOBwe996CiJS2SZtjKPUhslKmz16ktE3KwBAIBHjzzTd1gSoC5xxvvvkmgUCg2E0RKS0FXMk8UpNyKGnGjBm0t7cTmwAfcDkKBALMmDGj2M0QKR3Nzd4itYoK6OvzZhk1NhatOQWdlWRmAeA3QDVeELrXOfd1M3sHcBcQBqLAx51zb6Wecy3QBPQCVzjnfjHce+SalSQiUjJiMZgxI3ttgt8P7e0FnRAz3KykQg8l9QAfdM4dB8wDTjOz44FrgEecc7OBR1K/Y2ZHAecCRwOnAT8wM5X8FJHJq7V18IK1eNw7XiQFDQzOszf1qy/144CzgA2p4xuAs1P3zwLudM71OOdeBV4ClhSyjSIikq3gyWczqzSzrcAbwMPOuaeAw51zOwFSt+9MnX4kkFlutD11bOBrXmRmm8xsk/IIIlLS5s/3ah1l8vm840VS8MDgnOt1zs0DZgBLzGzwpgH75arWNigJ4py7zTm3yDm3KFRii9JEpPSN6QSiUAg2bIBAAGpqvNsNG4q64Hbcpqs653YDj+HlDv5iZkcApG7fSJ3WDszMeNoM4PXxaqOIyIEUZCuExkbYsQMefdS7LeKMJChwYDCzkJkdmrofBBqA54EHgBWp01YA96fuPwCca2bVZjYLmA08Xcg2ikj5yPVNfyTf/gu6FcKIt3ArnEKvYzgC2JCaWVQB3O2c+5mZPQHcbWZNwA7gHADn3LNmdjfwHJAELnXO9Ra4jSJSBtJLBfx+b9JPJOIdH3hsuC/r47x9S9FMyuqqIiKZYjFv2Cfzgh4MgnPQ3Z19rK1t6Iv8UK8z3HMmqmKuYxARKbr0N/1MFRUwcKv09Lf/oaS3QpgZiHFyzUZmBmLF2L6l4CZlSQwRkUy5Nj3r6/N6DJny2QitkWbOtSZ68VNpcYwIUNxk8VhTj0FEJr2hNj1bt26EG6Glss/W1UXVvj3YmGafJw71GESkLDQ2QkPD4E3Pch0bUplknxUYRGRSyrXzZa5Nz0a0EVquMal8xp9KjIaSRGTSKcgiNBh6TGoS9RZA01VFZJIZlymlubojJWa46aoaShKRSWVc0gAjGn8qPRpKEpFJpUzSAAWlwCAik8pQaQCYMFsqT3gKDCIyaaQL4jU0eDmFlhbvFgqUjJ6kFBhEZFIYOBOppcUrVgoFrIg6SSkwiEjJG64cdq46SQeqiVTuFBhEpOQNd/HPTEZPJ8YiNnJIPKZk9DAUGESk5NXW7i+fPfDin05Gf9rXTBv1tLCMV/rqCbUo0TAUBQYRKWnNzbBwoVdG+1xyX/wbG2Ksr2piCl0cwh6q4ko0DEeBQURKx4B9ODNzCzVdMSIMcfGPRjElGvKmwCAipSFHAaR0bmE6Mb7BVwnSlf2cigpobdWqtxFSYBCRiW+IaUezamOc1ekNH32WtdiAp7l9++Dss725q7feCtXVXkIis/jdgF6IqFaSiJSCaJQ4fvwZPYI4Pqa/1krEmqga2FNIMfCCyIoVUFXldS/icVi92tugobnZCzjp45GId7zMqccgIhPertowya7soaBkV4Ldu6Eq6M/9pEyJhBcgOjqgpweuugq2b9fKtyEoMIjIhPfq3hCXBSN0EmQPU+kkyKWBCK8eOn9Q7sClfobl88HTT2vl2xAUGERkwguH4U4aqaeNBlqop427rJEZ87Mr5iV9Qb5R+U2urF5DJ0HiwakQCAwOAIkELFmihPQQlGMQkQkvvUitqSnEi74QiUTGxmkZmzlXhcNcSohoFDprlzN9b9S70Le0eMNEPh/9T547N/2i2ccn8T4L+dIObiIy7mIxaG+NESbKtPnhvC/GB7Vx2lBPngS7sY2GdnATkQmjuRkeWtHMDxNNJPCT9MepWp/fbKCD2jhtqCdP8t3YRkM9BhEZN7EYLHx3jOe765mSMcXUBYPYmG7KLAcyXI9ByWcRGVvDLBiLRuE9ld6ahEw9fZoNNJEoMIjI2MlRtiJTOAwvJcP4yZ4N1NeTYFdtePzaKcNSYBCRsTHcbjkpoRCs+kqIleRYk7BXw0gThQKDiIxcruGiPLdKW7UKHg80cBY/5WPc078mQcsHJg4FBhEZmaGGi/KsYBpqaeaVvnru5ePcz9l8xNei5QMTjGYliUj+YjEvGHRlFK0LBiE9oyhdlC5zwVjmNNQcz9eMpOLQOgYRGRvp4aLMwJAeLgqFoLGRXfMaeKklStfhYY45JUToAM+3zOfLhKDAICL5C4ehszP7WFdX/3BRczOsWBEikfAu8n4/rF+f0WnQhjklQTkGERmRXmc5f4/FYOVK7zqfFo8PmJiULnqUKnqXtWGOTBjqMYhI3t5qjVKRDHJIxjqErqSP5KOtRGedSmXl4OdUVAwYKcooeldu9YlKRUF7DGY208weNbPtZvasmV2ZOv4OM3vYzP6Uup2W8ZxrzewlM3vBzD5cyPaJyMhECeMbsDithn3Ufeps5vx6LfMSG5lOjOnEWIR3v68vx0hRKASLFysoTFCFHkpKAl9wzs0FjgcuNbOjgGuAR5xzs4FHUr+Teuxc4GjgNOAHZpbjO4iIFMOM+SEu9kXoJNC/GY4BlfEuaq++mEf4EO0cSTszeJhltFHPw03Nuv6XmIIGBufcTufcltT9DmA7cCRwFrAhddoG4OzU/bOAO51zPc65V4GXgCWFbKOI5C8UgrO/18D1lf/MPqZkPWZAdbyDahJUE+dQ9jCFLpZGtF1mqRm35LOZhYH5wFPA4c65neAFD+CdqdOOBF7LeFp76tjA17rIzDaZ2aaY/sGJHJxhit4N0tzMP1xVzy1V11ND54HPB22XWYLGJTCYWS3wX8DnnHNvD3dqjmODVuA5525zzi1yzi0KqY8qMnoDVzGvXTt0kIjFSJ7fhHV1Ud3TgeH9cXZQM/wey5qOWnIKPivJzHx4QeFHzrn7Uof/YmZHOOd2mtkRwBup4+3AzIynzwBeL3QbRcpSZtG71IIzd/HF9E6po9IlsVtvhQUL+mcOvdUapSLu55CMfRQ6qOWL/u/j4t2s5ioS+AjQBRi+ugAVSW2XWYoKGhjMzIAIsN05952Mhx4AVgC3pG7vzzj+YzP7DvAuYDbwdCHbKFK2cq1CBqo6OwAvSPRNqaPCJbFIhOhhDcwZMCOpil4++O2PsscfYs7nlvOeyigv94b57ndh+YKopqOWqEL3GJYCnwKeMbOtqWPX4QWEu82sCdgBnAPgnHvWzO4GnsOb0XSpc663wG0UmbwG7mec8fuu2jDTuuMMNe3PgMpUkEie38R7b7qVCpL9w0Zx/Fzsi/Cv54YIhWD58hDRaCgjFigglKq8iuiZ2edzHN4DbHbObR3rRo2EiuiJDCFd0M7vh54eOPNMePBBqK4m2RVnpYsAsCbRRJIq6ujImeQDeJs66qrjWE9P/7FOAry3agf/+h+hfLZrlglmuCJ6+QaGHwOLgAdTh04HNgJzgHucc/88Rm0dMQUGkWyxGLS3xph3dj2WWexugE6C1NMGQJgo89jCaq4iSSV17M0KEt1Ug99PIN7Rf2wPU2mghWeDi1Fx1NIzFns+HwYscM59wTn3BbwgEQL+D3D+mLRSRA5aepLR55dHebvLP+y5CXyEibKLEJtYzO2sop42Tq/+FRexJmuHtS9UraYvnsx6vo+EtxJas1EnnXwDw7shK+uUAOqdc11AT+6niMh4ypxktG3f4NIVA6Uv7Jn2BUN89YHFLFqzijmBNs6uaWFOoI0jvr6Ky4LZ23GuJMIuQpqNOgnlm3z+MfCkmaVnD50JNJtZDV6iWESKJZVQbn8rjN8foqYrRpgoV3Irq/kcQboH5Q4ccIc10eEPQY9X5BS8maWnnurdz0wmA9Tf3MiDNBAmSpQwuwgRCGg26mSU9w5uZrYQOBFvssLvnHMTYmBfOQYpaxkJZheP88OeJs7vixDHj584V1fdyj99MUbtd2+C7u6sp3YS4Jyq+/no9fNZcnrogDNLB27Odt113v7NCgql6aCTz6kXORGY7Zy7w8xCQG2qnlFRKTBI2cq1TSbZ5QOS/iBV7W3w6KOwYkVWcHDAXmqopI/kmghTVx14atHA2a9Sug46+WxmXwe+DFybOuQD/nNsmicio5JeoJZh4JBRVXWVV+bi/PNxA3oMBtSxjyl0UXNlfoXuVC27POSbfP4H4O+BfQDOudeBukI1SkTyUFs7aHhokI4OuOEG6OrqDxqOwQXIzK+pRbJfvoEh7rwxJweQSjqLSLE0N8PChd72aACBwJCnusy9NoF9TKGH7J5GRVJTi2S/fAPD3Wa2FjjUzC4EWoB/L1yzRATIXRI7R/E7entJBGrzekkfSYy+/T0Hv19TiyTLAaerpgrh3YW3yvlt4H3A15xzDxe4bSLlLbOkRTzuXbwbG70hn6rsP90+fzW9+xL4crxMHB+9VJHAh48ElSSpJmOxWkWFtwezSMoBewypIaSfOuceds5d7Zz7ooKCSIFl9gr27PFum1IJ4i1bvNxBBpfs5UvVq+kkSCdBHN501E6CrGAD9bTRQAtn8VO6Buy8ht+v/IJkyXeB25Nmttg5t7GgrRERz9q1WdNQAW8BQWsrXHXVoNM7b7yV27+2imaWEyZKB7XUsbd/IRrALkJMJ4Z/4IpoLV2WAfLNMZwCPGFmL5vZH83sGTP7YyEbJlK2YjG46abBx1NJ5PiAxHEiUEfdSQuIRLySFi9OXcyO4FyOv2wxHdXZeYNdhFjJ/tIWvdVB5RdkkHx7DB8paCtExBOLwc9/7vUOBk5Fve46ds2cz5Su7NCQ6E6ypzZMY6OXKshcgHbJJTB/vld1O+0uGnmEBuZUR/lJa5jpcxUUJFtePQbnXBvetpsJ9k+Dzm/JtIjkJ10a9fLLB+UQCARg1Spe3RsaVMzu0uoIr+71Lu4DF6DNnQt33OHVQkrXQwoEvJ7FJXcsVlCQnPLqMZjZ5cDXgb8AfanDDji2QO0SKS+ZyeZMNTXQ19c/3BMG7mRAMbueEH9z39BlKjJ7ErW1sHevSlrI8PLdqOcl4P3OuTcL36SRUa0kKVnbt8PTT8OSJd7VetkybwZSigPw+bDvfc+rVpeydi1cfPHglwsG989oFTmQsdio5zW8rTxFZCxcfjkcdRScf753u2aNt1YhgwGWSJC84qqsBW4LFkBdjoI0mTNac62LE8lXvoHhFeAxM7vWzD6f/ilkw0Qmre3b4fvfzz62bh3ccAOuunpQ8m5f3MdbrdH+38NhSCbJyefzehT19V4HpL7eS12IjES+gWEH8DDgxyuel/4RkZF6+umch9+uDvHTr7XSQ3XW8YE7rYVC3pBRrvJIiYQ30zXXujiRfOWVfHbOfaPQDRGZdNKbFwzM+C5ZMuhUB1z9+QR/9r3Gz1nNaq7qL2FxsS/Cv87PzhSnE8pr18LNN2dvnvPtb2fPdE3vyaxks+Rr2OSzmX3XOfc5M3uQHNNTnXN/X8jG5UPJZ5mQmptxTU309UJFvIu+6iAVFWDp7PDll/cPJzmgF6My9ScWx8dlfI8Xgwt41YX5p3WhYRPKmZvnwKC9ewgGoa1NgUGyjXoHNzNb6JzbbGYn5XrcOffrMWrjqCkwyIQTi5GcUU9VvGvQQ/07qoVCXq6hpYW+L1xNRaIn67xOAvzoph0c8t4Qhx7qLVLL98I+cAtOzVSSXIYLDMMOJTnnNqduix4ARErFW61RquO5/7j2xX30tUaZdmqI2PS5xKbuZU5VFQwIDH1UcvtXojztvGjg98P69fld4HOtgBYZiWEDg5k9wzArnJ1zWuAmk18+Gx1nnNO+u5ZjGNxbAC+R/AJhHmqGlSvh8Iow27v6CA44r4JeXnHh/t/jca8X0NCQ34U+FFJAkNE7UPL5jHFphchENdSeCMOcM7th8J+NA3qo5mJfhGtnhlhxhjfM00aIzxBhA+f3Vz2N42Ml6/qroqZVVCiJLOMjr5XPB3wRsyeccyeMQXtGTDkGKZhY7MCZ3Fzn5BDHx0k8yqLLlnLmmfDhD2c/Pp0Y82gFYCvzBwWFXG8tcjBGnWMYgaE3nBUpVdGo1wvIvOin536mH3/rrcHnDJCedfQIy7h4bYTdfzc4UbCLEC2cmvU2vb1emSTQ7psyvsYqMKjSqkw+4fCgMhUkErBlC+6kk+i1KiqTPdgBet0GBFPDRGsSTTxFA35/aNBL19Z6K5qvv35/aaRWrxMxollJIgdrrAKDyOSTXmKcOffz1ltJXnEVVfGu/j8eB5jP5431JBLQ1ISLREj2VVDVsw/LeMkEPuYdGmX9+hBNTV7eoK8Pbr3Vq4E0ML996qmIjLu8SmKY2WVmNm24U8aoPSLjJq9Cc42N3sB+Swu0tfHWrAV0xrO/TxngzOCee6CtjeYPfI8jE218tOc+ugaMstb4E0ybH+5/2Ucf9W5XrcreR0GkmPKtlfQ3wEYzu9vMTjOzgYHgU2PcLpGCSu+Jk7PQ3MCIkbH7TZTw4D2Tgb6qapg2jRghVq6EnUkvZ7CSdf2b6vQFglSt358oGLipjshEke8Obl8BZgMR4HzgT2Z2s5m9J/X4toK1UGSMZe6JM6jQ3DARIxaDux8NcS03DkqqVbgkhMNEo1BZuf/4XTRSTxtnBlr4w/1tWoIsJSHfHgPOm9f659RPEpgG3Gtm/1ygtokURHqyUSafD9pbh44Y6XgRvaWZm/gaXXjlsTsJ0EmQRxojbIyGqK31ZhNl2kWITbaYGfPVNZDSkO/WnlcAK4BdwO3A1c65hJlVAH8CvlS4JoqMraEmG4WJDj7ZOd5qjdLUFKKmK0aEJqZkrGquwDGPLbywbi5193izii64wKt6mkh452iqqZSafGclTQeWO+faMg865/rMTKujpaSkJxt9eWWMBdZKTXI3F1x5KNNqagavR+ju5tkdtVRUwDxa6R3Qye6hmjr2AtDR4R2LROAPf4DXXvN+11RTKTVjsvK5mLTyuczlU8col+Zm3IrzIeF1HQy88STY/1UfrxrqyfyamfGXiNBEkK6sKXidBKmnLWul8tSp3iSmxYtH+d8kMg7GYs/n0b7xOjN7w8y2ZRx7h5k9bGZ/St1Oy3jsWjN7ycxeMLMP535VkZRhpxYNI5V9tkTc21c5fTyRyAoK4A057YrX9g8hpc9N5xcuqY4MKl+RSOzfG0GkFBU0MADrgdMGHLsGeMQ5Nxt4JPU7ZnYUcC5wdOo5PzCzSkRyGXZq0f5Tcq5TiEa9lWU5JCqr6aaaPUylkyAX+yLUsZc42dlqq6khftf9XPrbRtas8da2TZ3q3SqfIKWuoCufnXO/MbPwgMNnASen7m8AHgO+nDp+p3OuB3jVzF4ClgBPFLKNUqJy1TGqquovP9qcKmtdWenNElq3LmOmaDg8eOpQSqLXWMAW6thLlDC7EiGmExu8dqGvj0NPmc/i1FqE5cu1/4FMHoXuMeRyuHNuJ0Dq9p2p40cCr2Wc1546NoiZXWRmm8xsU0y7nJenXFOLOjpgyxZiMVixAmq7Y8zdt5Ha7hgrVmT0HEIh9n7u+kFrERxwI9fzAnPZxGJ2EcLng47qEJcGInQSJB7M3S3QYjWZTIoRGIaSq6xGzsy4c+4259wi59yikP4Sy1MoBLfeOvgfyFVXse3RGOcn1vIaM3mUk2mjnuWJ5v6CdAB/OmUV3QO2x+kiyL+zKutYIgGPPAKX/KaRzufa8P/aK42hhWoymRUjMPzFzI4ASN2+kTreDszMOG8G8Po4t01KyC9iC+igLvugz8d77riOtVxMgB5q6WQKXfwHn2bKju39p82YH2KVL9JfrqKTICuJsLsq+4tGIOCNWC1eDNPnqlsg5aEYgeEBvMVypG7vzzh+rplVm9ksvBIcTxehfVICYjG48KYwVSSzjruODmY+dPug7qePJEsvnd8/cykUgo9saGROoI2zprTwv/1tLPqXxv4Zq2lmmmEk5aegyWcza8ZLNE83s3bg68AtwN1m1gTsAM4BcM49a2Z3A8/hldy41DmXO0Mo5SXHWoVoFHb7QqzsjrCBT+NPBwjnco5JGkC8J2vj5MZGaGgIEY2G+l/6yCOzq2xrhpGUo0LPShpqIPZDQ5x/E3BT4VokJWeIPZe3bIHqjhjv50n8JPOv+57egS2jwmnmhd8LFpphJOVNG/XIxJW5ViE9LbWpiV3zGnj88hb+hxX4SAwbFBzZsxqS3QmqMsaGci2cHhgsRMrNRJqVJGXmQBvlvNUaJVkxuAzqmy2trEk04R8iKLiMn6dZkJ1gdhFiqZXKo104LTLZqccgRZFzhKhh/9f35pYQX14Z5vnueNY/0r7uHqYkdw8qZpfmgK/xdXYwi6dYwiu+uRySiBEmSpQw8WCIy6PeuTk6I+n0g0hZU2CQcZdrhOihFc2cW9WE+f24eJz/l4zwWqKRlURYRxMAQbrojldw2OdX4LfeQatcHPB5/oXv8sX9BxP01zIKE+UvcQiHQzkXTg9IP4iULQ0lSUHlGi5KX5SnE2MRG3kf2/lhoglL1T2yri7WJpqYToy7aGQBm6mgDwOm0MUUuul1EMfXP2QUp5KLWJMdFPDWIZxX2Uwb9bSwjFf66gm1NA+9J0O4oB+HSElQYJCCGWoMPxyGszq9i/XDLGMr8wc9N0AXF7IWgDr20k0g6/Fugtzz6Z9xpv8XnMovOJKd3D5g1TJAiBgbfF5l1EPYQ1XcK7YXIkYkouJ3IrloPwYpiFjMCwaZQzXBoFdNIkSM5Ix67yKdMnD2EEAn1ZzFA7zGTLawMGvntE6CdD7XhpseYu1auOmm/UNDZl5PIZGAn1y3kVNuXoa/a0//c+PBqV5pi8WLR72dg0ipG24/BuUYJD8jvIIOO4ZPlKqgHzICg6Wz0BmC9HAfy6mkj3+niQuJkMCHjwStl0VYOtdrx1e+AqtW7W9e+v3DYbBdYZJfzS6anexK8HZtmOloaqpILhpKkgMbxbzO4cbwd9WG6e0e8GBFBcmq7OEiA+rYxxS6uJAIC9jMmYEWnrqrjaXfy147mVndNPP+q3tDXBbMrol0aSDCq3sVDUSGosAgw8tjQ5xc0vsqDxzDb2mBdy8McWGFd7FOBOqgupq3v/ldVto6Ogmyl5pBVVOTVFHHXjbZYo45Jf+LejgMd9JIPW000EI9bdxljUoyiwxDgaHMHWiRWf+YUKb0mNABNDZ6OYWWFti8GQ47zNs8p6YrxjNd7+Ur3EBvdw+9VsmU6z6HGdTTxj9wH10DSmLX0cGSqi0jThCnA9S+YIgXpy5mXzCkJLPIASj5XMaGKEOUbdgsch5X11iMX6yNcuFNYf5aGeLMfc1EaCJBFVPpyEo49+BjBv/DLkJcwFpu4+Ksx10wiOX7vjn+M5RkFtlvuOSzegxlKu8RoqHGhPK5ujY34+rrOf6ry3i+u57GfWuJkJ462jFoFpKfBPPwdtPZygLi1dl7LViePZVctMOaSP4UGMrUiEaIMseE8t29LBV5rMtbPzCFLr7P5SQP8E/uZB4FYFfN4L0WtAJNZHwoMJSpEa/8HelX7hyRx0+COvYN+RQDrmI104nxl74Q+1ZrBZpIMSgwlKmRjhANl6TetT3Gcxs2smt7xoM5Io+lfhzQV1NLnMpBs48S+JhTHSUSgamrRtFTEZGDpsBQxvIdIRpuGcPvL29mylH1HHn+MqYcVc/vL089mI481dWDXs/q6tj7re+zoPKZQbOPaqsT/KQ1vL8tSg6IjDvNShJg6Fk7w01Ksl0xphxVn7NUxfTUqmS2b4f586GnJ+sFWn/axikfD3HanmbW0dS/onnnNyO85yvqGYgUmmYlybCG6xEMl6R+4+koCbIfTODjjaej+w/MnQt33DFozGrG/BDxONyVsfhsTqDNGz4SkaJSraQyN8TumTQ0ePffemvoJLXVhvGR/aCPBO9cEs5+Qo6NlEN4I01NTRD3hXg2oYVnIhOFegxlbqgewdq1sPDdMb61fCOHJmL4/YOT1G56iHXWlLWV5jprwk3PcXXPkSsYzSxYESk8BYYyN9S01Re/0czz3fX8ZN8yXkrW8/G+Zu65J/sC3t4aY6WL9M82MmCli9DeOnwdpUzKLYtMPAoMJeCA9YwO4vxc01ZvvDLGmqS3QvnQ1OK0tUlvc5vMC3iY3DmGMNH8/+NEZMJRYJjgRlrxehQVsmlshB2bYzzxbxvZsTnGZ07J74I/bX6YGn92d6PGn2Da/HAe/2UiMlEpMExgI614PcoK2dDczPSF9Rx15TKmL6xn2qtb8rvgh0JUrY/ggkGSNVNxwSBV65VBFil1CgwT2EgrXo+qQnauaHLVVVT92635XfAbG7G2NqoebfEqnyqDLFLyNF11AhtpPaMR1z8CaG2FigHfD3w+WLDAu+DnU6ta+2OKTCrqMUxgI61nNPD8QACuu26YN2huhrPPhn0DCtulo4mmDImUJZXEKAEj3WQmFvPWIdx88zCb8OSqdQFeNFm3TkNCIpPccCUxNJRUAkYzUnPzzblXM/e/TjohkRkYamrgvvvg1FPHotkiUqI0lDQJ5ZWEzpWQ6OvzCt6JSFlTYJiE8kpCH8yWnSIyqSkwjJORrl4+GHlf81WsSERyUI5hHDQ3e2P8QyaCC6CxEZbNi/HG01HeuSS8f3+EgTTVVEQGUI+hwEa9Gvlg3/TGG7NWM+dVG0NEBAWGghvVauSD0dwM7343fPWr4xyNRGSyUGAosFGtRh4o3wRFunvS3T34sYJGIxGZTBQYCuygJ/+MpFxqru5J2oijkYiUqwkXGMzsNDN7wcxeMrNrit2esTDqyT8jTVDk6p6ApqKKyIhMqMBgZpXA/wU+AhwFNJrZUcVt1dgYVdmhkSYochVL+uY3NRVVREZkok1XXQK85Jx7BcDM7gTOAp4raquKZTQJisZGr/bFSIoriYhkmFA9BuBI4LWM39tTx7KY2UVmtsnMNsUm80yb0SYoVBVVRA7CROsxWI5jg8q/OuduA24Dr7pqoRtVVOoBiMg4m2iBoR2YmfH7DOD1IrVl4tDqZBEZRxNtKGkjMNvMZpmZHzgXeKDIbRIRKSsTqsfgnEua2WXAL4BKYJ1z7tkiN0tEpKxMqMAA4Jz7OfDzQr/PSHdFK703FBEZnYk2lDQuRrKYuDTfUERk9Mpuz+dcWx0Hg94asIJ8kR/3NxQRObDh9nwuux7DuFc7Hfc3FBE5OGUXGMak2umEfkMRkYNTdoFh3Lc61t7KIlJiyi7HkKZZSSJSzobLMUy46arjZdwXE2v1soiUiLIbShIRkeEpMIiISBYFBhERyaLAICIiWRQYREQkiwKDiIhkUWAQEZEsCgwiIpJFgUFERLIoMIiISJbyDQyxGGzc6N2KiEi/8gwM2lFNRGRI5RcYYjFoavJ2VNuzx7ttalLPQUQkpfwCg3ZUExEZVvkFBu2oJiIyrPILDNpRTURkWOW5UU9jIzQ0aEc1EZEcyjMwgHZUExEZQvkNJYmIyLAUGEREJIsCg4iIZFFgEBGRLAoMIiKSxZxzxW7DQTGzGNBW5GZMB3YVuQ0ThT6L/fRZ7KfPYr+J8lnUO+dyTs0s+cAwEZjZJufcomK3YyLQZ7GfPov99FnsVwqfhYaSREQkiwKDiIhkUWAYG7cVuwETiD6L/fRZ7KfPYr8J/1koxyAiIlnUYxARkSwKDCIikkWB4SCY2Uwze9TMtpvZs2Z2ZbHbVGxmVmlmrWb2s2K3pZjM7FAzu9fMnk/9+zih2G0qFjO7KvX3sc3Mms0sUOw2jRczW2dmb5jZtoxj7zCzh83sT6nbacVsYy4KDAcnCXzBOTcXOB641MyOKnKbiu1KYHuxGzEBrAYecs7NAY6jTD8TMzsSuAJY5Jw7BqgEzi1uq8bVeuC0AceuAR5xzs0GHkn9PqEoMBwE59xO59yW1P0OvD/+I4vbquIxsxnA6cDtxW5LMZnZVOD/ABEA51zcObe7qI0qriogaGZVwBTg9SK3Z9w4534D/HXA4bOADan7G4Czx7NN+VBgGCNmFgbmA08VuSnF9F3gS0BfkdtRbP8LiAF3pIbVbjezmmI3qhicc/8DfBvYAewE9jjnflncVhXd4c65neB9uQTeWeT2DKLAMAbMrBb4L+Bzzrm3i92eYjCzM4A3nHObi92WCaAKWAD80Dk3H9jHBBwuGA+p8fOzgFnAu4AaMzuvuK2SA1FgOEhm5sMLCj9yzt1X7PYU0VLg780sCtwJfNDM/rO4TSqadqDdOZfuPd6LFyjKUQPwqnMu5pxLAPcBHyhym4rtL2Z2BEDq9o0it2cQBYaDYGaGN4683Tn3nWK3p5icc9c652Y458J4ycVfOefK8puhc+7PwGtm9r7UoQ8BzxWxScW0AzjezKak/l4+RJkm4jM8AKxI3V8B3F/EtuRUVewGlLilwKeAZ8xsa+rYdc65nxevSTJBXA78yMz8wCvAZ4rcnqJwzj1lZvcCW/Bm8bVSAiUhxoqZNQMnA9PNrB34OnALcLeZNeEFznOK18LcVBJDRESyaChJRESyKDCIiEgWBQYREcmiwCAiIlkUGEREJIsCg4iIZFFgEBGRLAoMImPMzBab2R/NLGBmNam9CI4pdrtE8qUFbiIFYGY3AgEgiFc36VtFbpJI3hQYRAogVQpjI9ANfMA511vkJonkTUNJIoXxDqAWqMPrOYiUDPUYRArAzB7AKz8+CzjCOXdZkZskkjdVVxUZY2b2aSDpnPuxmVUCj5vZB51zvyp220TyoR6DiIhkUY5BRESyKDCIiEgWBQYREcmiwCAiIlkUGEREJIsCg4iIZFFgEBGRLP8f0ZKvQwz35X4AAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "posterior_graph = model.get_posterior_graph(\"posterior_graph\")\n", "\n", "prediction_model = ForwardModel(graph=posterior_graph,\n", " data={posterior_graph.x: data[\"x\"]})\n", "\n", "y_linreg_prediction = prediction_model.get_means(posterior_graph.y)\n", "\n", "data[\"y_linreg\"] = y_linreg_prediction\n", "\n", "ax = data.plot.scatter(\"x\", \"y\", color='blue', label=\"true\");\n", "ax = data.plot.scatter(\"x\", \"y_linreg\", color='red', label=\"lin.reg.pred.\", ax=ax);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The plot shows that our linear regression model correctly predicts the trend seen in the true data for values of x between 3 and 9. For smaller or larger values, the predictions are significantly off due to the curvature in the true relation between x and y. To also capture that curvature in the data, we need to go beyond a linear model." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Regression model using convenience functions" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A linear regression model with just one scalar feature and one scalar target can be quickly built in the manner described in the previous section. However, building a regression model using beyond-linear polynomials, multi-dimenensional features and targets, and/or multiple features and targets that way can become very involved very quickly. To facilitate building more complex regression models, we can use the convenience function `connect_via_regression`.\n", "\n", "For example, this creates a graph for a quadratic regression model:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "result_shape = ()\n" ] } ], "source": [ "graph = Graph(\"graph\")\n", "with graph:\n", " x = Variable(\"x\", shape=(), mean=x_data_mean, variance=x_data_std**2)\n", " y = Variable(\"y\", shape=(), variance=y_data_std**2)\n", " \n", " connect_via_regression(\n", " name_prefix=\"regression\",\n", " inputs=x,\n", " outputs=y,\n", " order=2,\n", " include_cross_terms=False,\n", " inputs_location=x_data_mean,\n", " inputs_scale=x_data_std,\n", " outputs_location=y_data_mean,\n", " outputs_scale=y_data_std,\n", " )\n", "\n", "hal.show(graph)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Using a `Trainer` and a `Predictor`, we can compute the predictions from our quadratic model and compare them to the linear predctions:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "trainer = Trainer(graph=graph, data={graph.x: data[\"x\"], graph.y: data[\"y\"]})" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "graph\n", "├─regression_y\n", "│ └─location\n", "│ ├─slope\n", "│ └─intercept\n", "├─inputs\n", "├─outputs\n", "├─x\n", "└─y\n" ] } ], "source": [ "print_child_tree(graph)" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array(1.)" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "graph.regression_y.location.slope.variance.operand.value" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEGCAYAAABhMDI9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAA48ElEQVR4nO3de3yT5fn48c/dHJrQUooSHQo0zDE5KVAEQUREizpxokyn9asDiRMPKOo8jfl18+x+OhFlU/waBtskeNzUTZ2i6HSoFCgOtIgwU62iBLQV2rQ59P798aRt0iYlhaZJm+v9evWV5slzuBPoc+U+XbfSWiOEEEI0yUl3AYQQQmQWCQxCCCFiSGAQQggRQwKDEEKIGBIYhBBCxDCnuwAHql+/ftrpdKa7GEII0a2sX79+l9baEe+1bh8YnE4n69atS3cxhBCiW1FKVSZ6TZqShBBCxJDAIIQQIoYEBiGEEDG6fR9DPMFgkKqqKurr69NdFJEkm83GgAEDsFgs6S6KEFmvRwaGqqoqevfujdPpRCmV7uKIfdBas3v3bqqqqhg8eHC6iyNE1uuRTUn19fUcfPDBEhS6CaUUBx98sNTwhOgAnw/KyozHztYjAwMgQaGbkX8vIZLn8UBREUybZjx6PJ17/h4bGIQQoify+cDlAr/yUZNXhl/5cLk6t+YggSEFqqur+cMf/pDuYggheiCvFzjKA9cWwc+mwbVF6BEeY3snkcCQAokCQzgcTkNphBA9Sf6hPvzTXGDxg60GLH7qT3WRf2jnVRkkMER0ZkfOzTffzPbt2xk9ejTjxo1j6tSpXHDBBRx11FF4vV5GjhzZvO/999/Pb37zGwC2b9/OaaedxtixY5k8eTJbtmw58MIIIXqUvSYv9lxrzDZ7roW9Jm+nXaNHDlftKI/HaLOzWiEQALcbSkv3/3z33nsvmzdvZuPGjbz55ptMnz6dzZs3M3jwYLzt1PcuvfRSHn30UYYMGcL777/PFVdcwRtvvLH/BRFC9DjOQifkBKAxamNO0NjeSbI+MDR35PiNHzCel5SAI27ewY4bP378Psfn7927lzVr1nDuuec2b2toaOicAgghegxHngP3DDeu511YTBaC4SDuGW4ceZ10w0ICA16vUVNoCgoAFouxvbMCQ15eXvPvZrOZxsaWUN80dr+xsZHCwkI2btzYORcVQvRYpSNLKRlcgrfai7PQ2alBAaSPAafTaD6KFgwa2/dX79692bNnT9zXDj30UHbu3Mnu3btpaGjg73//OwAFBQUMHjyYp59+GjBmA3/wwQf7XwghRLfS0X5OR56DcYeP6/SgABIYcDiMPgW7HQoKjEe3+8BqCwcffDCTJk1i5MiR3HDDDTGvWSwWbr31Vo499ljOOOMMhg4d2vzaE088gdvtZtSoUYwYMYLnn39+/wshhOg2Uj1hraOU1jq9JThAxxxzjG69UE9FRQXDhg3r0Hl8PqP5yOnsvCYk0TH78+8mRHfn8xnBILo5226HysrU3ouUUuu11sfEey3r+xiaOBwSEIQQXa8r+jk7KuubkoQQIp2a+zkProBRy+HgigPu5zxQUmMQQog0cjjgiKuuYrN9cfO2H/jn4XA8nLYySY1BCCHS6N9bK4ygoGj+2WRfzL+3VqStTBIYhBAijV7dvLZD27uCBAYhhEijU0aO79D2riCBIUXy8/MB+PLLLznnnHPSXJrO8+abb3LGGWekuxhC9BiTfjiMUwrngab555TCeUz6YfqGbkvnc4oddthhPPPMM/t9fCgUwmxO/T9TOBzGZDKl/DpCiLb+ec3D/HvrFby6eS2njByf1qAAUmNokaIFVKPTbC9btoyZM2dy2mmnMWTIEG688ca4x5x44oksWLCAKVOmsGjRItavX8+UKVMYO3Ysp556Kjt27ACgrKyMo48+mokTJ3LDDTfEpPOOvv7QoUOZNWsWRx99NOeccw51dXUAOJ1Obr/9do4//niefvppXn31VSZOnEhxcTHnnnsue/fuBeCVV15h6NChHH/88Tz33HOd+vkIIQyTfjiM22bOSj4opHDRZwkM0KXz0Tdu3MiTTz7Jpk2bePLJJ/n888/j7lddXc1bb73F1VdfzVVXXcUzzzzD+vXrmTNnDr/61a8AuPjii3n00Ud599132/22//HHH3PppZfyn//8h4KCgphFhGw2G++88w4lJSXceeedrFq1ig0bNnDMMcfwwAMPUF9fz89//nNefPFF3n77bb766qvO/UCEEB2X4nuWBIbovNs1NcZjZy+gGuXkk0+mT58+2Gw2hg8fTmVlZdz9zjvvPMC4qW/evJlp06YxevRo7rzzTqqqqqiurmbPnj0cd9xxAFxwwQUJrzlw4EAmTZoEwIUXXsg777zT5jrvvfceH330EZMmTWL06NEsX76cyspKtmzZwuDBgxkyZAhKKS688MJO+RyEEPupC+5Z0sfQxfPRc3Nzm383mUyEQqG4+zWl6tZaM2LECN59992Y17/99tukr6mUSvg8+jrTpk3D0+qbx8aNG9scL4RIoy64Z0mNIRV5tzvRkUceic/naw4MwWCQDz/8kL59+9K7d2/ee+89AFauXJnwHJ999lnz8R6Ph+OPP77NPhMmTODf//4327ZtA6Curo6tW7cydOhQPv30U7Zv3958vBAijbrgniWBIRV5t/fDJZdcQusssQBWq5VnnnmGm266iVGjRjF69GjWrFkDgNvt5tJLL2XixIlorenTpw9gDJE9/fTTm88xbNgwli9fztFHH80333zD5Zdf3uY6DoeDZcuWUVpaytFHH82ECRPYsmULNpuNxx57jOnTp3P88cdTVFTUfMy6deu45JJLOvujEEK0pwvuWZJ2u0k3zLu9d+/e5vkS9957Lzt27GDRokUx+3i9Xs444ww2b96cjiJ2iKTdFqIDDvCeJWm3k9EN827/4x//4J577iEUClFUVMSyZcvSXSQhRFdJ4T1LAkM3dt555zWPKkrE6XR2i9qCECJzSB+DEEKIGBIYhBAiE6RwJnNHSWAQQogO8tX6KPuiDF9tJ93Em2YyT52a8uwLyUhpYFBK2ZRSa5VSHyilPlRK3RbZfpBS6jWl1CeRx75Rx/xSKbVNKfWxUurUVJZPCCE6yrPJQ9GDRUz78zSKHizCs/kAb+I+H8yebUxYq601HmfPTmvNIdU1hgbgJK31KGA0cJpSagJwM/C61noI8HrkOUqp4cD5wAjgNOAPSilJ+Rll9uzZB5StNR2iEwkK0Z35an24XnDhD/mpaajBH/Ljet51YDWH8vK2E9YCAWN7mqQ0MGjD3shTS+RHAzOA5ZHty4GzIr/PAFZqrRu01p8C24D0rVbRTSVKs9FdryNEpvBWe7GarDHbLCYL3mpvegqUIinvY1BKmZRSG4GdwGta6/eBQ7XWOwAij4dEdj8ciE43WhXZ1vqclyql1iml1vk6qbrV6W2GwF133cWRRx5JSUkJpaWl3H///Zx44onNM5x37dqFMzKN3ev1MnnyZIqLiykuLm6e3ay1Zt68eQwfPpzp06ezc+fOuNeaPXs21113HVOnTuWmm25i+/btnHbaaYwdO5bJkyezZcsWALZv386ECRMYN24ct956a/MEudby8/P5xS9+QXFxMSeffDJNn3OyKcHXr1/PqFGjmDhxIr///e877TMVIp2chU4C4dhv98FwEGehc/9POmaMkesomsVibE8XrXWX/ACFwGpgJFDd6rVvI4+/By6M2u4GftLeeceOHatb++ijj9psa8+K/6zQ9jvtus89fbT9TrtesWlFh46PZ926dXrkyJG6trZW19TU6COOOELfd999esqUKbqsrExrrbXP59NFRUVaa61ra2u13+/XWmu9detW3fS+nn32WV1SUqJDoZD+4osvdJ8+ffTTTz/d5nqzZs3S06dP16FQSGut9UknnaS3bt2qtdb6vffe01OnTtVaaz19+nS9YoXx/h555BGdl5cXt/yA/stf/qK11vq2227TV155pdZa6ylTpujLL79ca611IBDQEydO1Dt37tRaa71y5Up98cUXa621Puqoo/Sbb76ptdb6+uuv1yNGjNjnZ9bRfzch0mHFJuN+kXdnQafdL/SKFVrbbFrn5RmPKzrhnPsArNMJ7qtdNsFNa12tlHoTo+/ga6VUf631DqVUf4zaBBg1hIFRhw0AvkxluaLbDP0hI1uh63kXJYNLcOTt/6zCt99+m7PPPptevXoBcOaZZ7a7fzAYZN68eWzcuBGTycTWrVsB+Ne//kVpaSkmk4nDDjuMk046KeE5zj33XEwmE3v37mXNmjWce+65za81NDQA8O677/K3v/0NMFJ1X3/99XHPlZOT0zx57sILL2TmzJnNr8VLCQ7GKnD9+/enpqaG6upqpkyZAsBFF13Eyy+/3O77F6Lb2FSKfqAE+nnRu5zwfYfxdfdAlJZCSUnGpOVJaWBQSjmAYCQo2IES4LfAC8As4N7I4/ORQ14AViilHgAOA4YAa1NZxqY2w6agAC1thgcSGKBtumsAs9lMY2MjAPX19c3bFy5cyKGHHsoHH3xAY2MjNput3fPE05RCu7GxkcLCQjZu3HgApY+VKFV3vJTg1dXVkqpbZKR46YU6knKoaSmEer8DvjF2drmMe/oB38szKC1PqvsY+gOrlVL/Acow+hj+jhEQpimlPgGmRZ6jtf4QeAr4CHgFuFJrHU5lAVPSZgiccMIJ/PWvf8Xv97Nnzx5efPFF43pOJ+vXrweIGV1UU1ND//79ycnJ4c9//jPhcLj5PCtXriQcDrNjxw5Wr169z2sXFBQwePBgnn76acC4gX/wwQeAkV772WefBdpP1d3Y2NhcvhUrVsRN1Z0oJXhhYSF9+vRpXhDoiSee2GeZhUi1eIuedXQhtKalEKI1LYXQk6R6VNJ/tNZjtNZHa61Haq1vj2zfrbU+WWs9JPL4TdQxd2mtj9BaH6m1Tnn7gyPPgXuGG7vZTkFuAXazHfcM9wHXFoqLiznvvPMYPXo0P/nJT5g8eTIA119/PY888gjHHXccu3btat7/iiuuYPny5UyYMIGtW7c2fys/++yzGTJkCEcddRSXX355c/MMwK233soLL7wQ9/pPPPEEbrebUaNGMWLECJ5/3qiUPfjggzzwwAOMHz+eHTt2NKfqBhg9enTz73l5eXz44YeMHTuWN954g1tvvbXNNdpLCf7HP/6RK6+8kokTJ2K325uPaZ0SXIiukGjRszlzOrYQWoYv39J5EnU+dJefzuh81lrrnXt36rVVa/XOvTs7fGwyfv3rX+v77rsvJefuiNraWt3Y2Ki11trj8egzzzwz7n6JOqVTSTqfRaqsXat1nz5aQ8tPXp7xE72toMDYtz0rVmg90LZTn5i3Vg+07eyKfuKUIBM6nzOdI89xwLWE7mD9+vXMmzcPrTWFhYUsXbo03UUSIuXifdNvbDTCAb18UOiFaifBoGOf3/5L8XC+chHGikkFULiB0pSUO11koR6RMeTfTaSSx2M0FVksRvOP2w1rvvOw+HMXNFohJ8C8gW4entvOTd7nMzojotdbttuhsjJjOo6T1d5CPT02iV53D3jZRv69RKqVlhr371WrjMeSM324d7nA4ofcGrD4ce/aR3qLLOl97pGBwWazsXv3brnZdBNaa3bv3h0zRFeIAxUvi7XDAePGGY/7ld4iS3qfe2Qfw4ABA6iqqqKz0mWI1LPZbAwYMCDdxRA9RFOzkdVq3MfdbqPG4Kv14a324ix07t9QdYfDOFnrNqlu1oy0Lz2yj0EIkb0SdQPc/o8l3PLv+VhNVkKNIdwz3ICR6cBishAMB3HPcFM6MomO5I7MistQ7fUxSGAQQvQoZWXGhLWampZtuZOW0DDtspj97GY7lddUAjTXIrJhZGKT9gJDj2xKEkJkrzbdAL18NEyd32Y/U44Jb7WXcYePy6qAkIwe2fkshMheTd0AtoN85A0pwzqoHJvF2ma/zkh901NJjUEI0WM0Nf1/N8iDus4FWMnRDYTCjW32XfSjRVJTSEACgxCiR2gaiZRzaAW1F14M5gYg0gMdsoC2Q9gM5gDm1xYx84q5aS1vJpPAIITo9pqT5A1bAqdfBaZg7A5hOzz5NNT3hWonvcwOvN5uO6Ao5SQwCCG6vfKPfYSnPgDj7oV4S4HkBOGrMfSrAydevraB0ylRIRHpfBZCdGueTR7OXF1EIEFQyDXlMm+gm58FV1FJEauYxn8bi3Cs2sfiC1lM5jEIIbotX62PAb8rIqD9cV+3YOKDKzYxjH7ooiJUD0h+11myMomeEKIHapUAqfxTLwF/26GoAGh4+GXFMPqB14vKguR3nUUCgxCie4i3Dme1E0ytktppsIbg0Rdh7uZcKC/PmuR3nUUCgxAi80WGHfmUn7K8GnzKWIez+GCY/KILexAK6sEehDteh6oHYO4G0LW1cNZZRq7thQshNxfy841mpKbkd/HSsGY5GZUkhMh8Xi9/HgZzfwTWMARMsOQlzUWfl/PGFjffbgNvITirwVHXcpgCI5verFlgNrekW120yEi3migNa5aTzmchRMbbUl5B8bPD8VtattmD8J/hT/KDKy6NzZiXDLsd1q+HsWN7xGps+0M6n4UQ3dr7u/YSCttjtoXCdt6zFbbpO9CRn3ZZLLB2bVasxrY/JDAIITLe+COdBE2x24ImGHvMGKP5x26HggJCFju3me5gfu6j1GEnYC8Am61tAAgGYfx46ZBOQAKDECLjDRvkYN5ANwTt0FAAQTvzBroZNsgRs5iz+YtKrtxxCxe9PZe6jyqxvrUKPvsMli1rDh7NHc/DhsUElZgO6SwnfQxCiC7n80FVuQ8nXvqOcSZ9M674zMfaj72MP9JpBIWOXjTeqms9YDW2/SEL9QghMobHA6/M8vBI0EUQKyFrAPOy5EYDDRvk6HhAaOJwxL/xJ9qexaTGIIToMj4fjB3kY0t9Eb1oGQ2k7XZUlowGyhQyKkkI0XXamTDm9cIRJi8BrPh6Qdlh4OsFDY0yGiiTSGAQQnSeeGkrojidsC3k5E/FtQy8Dk7+GRRdC08O8bMr35mWIou2JDAIITpH82o5fmPCmd9IWxFdc3A4oPiXzzH/xyEazLDHBn4LuGYoNuxOY9lFDAkMQoiOi9dc5PXuc8KYr9bHP03z26ybEG7MhUIvIjNIYBBCdEyi5qIkMph6q71YG9uupmOxBRgz2Nlmu0gPCQxCiOS111zkcOxzwpgznE8oUB97Tg0Pn3AnjjwZkZQpJDAIIZK3r+ai0lJ2ra/kvTtXsXpZJb6S2LkJjq/34n7Vjj0I+fWQG4JHX8tlbt6ULim+SI5McBNCJM/phLq62G1+f3NzkccDP7scQgd/C3yL5ZsxLP+Do2XumtNJ6WYo+TAqTbbOkfxEGUYCgxCiQ8JaYYrzvOIzHxf93xLCV90GphAAwbCVWfcto6Sk1GhRijQ3OVwuHHstoIOSnygDSWAQQiTt23IvOSE7fWjpZPaHLKz82++Y99VDhE/wx444MgcInu6i/OMSTmm6+ZeWQklJVuYn6i5S2seglBqolFqtlKpQSn2olJof2X6QUuo1pdQnkce+Ucf8Uim1TSn1sVLq1FSWTwjRMV6cWIgdeVTXq5arK39LQ6O/zTBUAHRO26GoDgeMGydBIUMlVWNQSj0UZ3MNsE5r/Xw7h4aAX2itNyilegPrlVKvAbOB17XW9yqlbgZuBm5SSg0HzgdGAIcBq5RSP9Rah5N/S0KIVBkwxsFlFjePBudgp55dveDlIWAOA5b4x1gtMhS1u0m2xmADRgOfRH6OBg4CXEqpBxMdpLXeobXeEPl9D1ABHA7MAJZHdlsOnBX5fQawUmvdoLX+FNgGjE/+7QghUsnhgLMeLuGyQ67jktNzGHQtXPUj2JPbasfIMmrWECx7IXYdZpH5ku1j+AFwktY6BKCUegR4FZgGbErmBEopJzAGeB84VGu9A4zgoZQ6JLLb4cB7UYdVRba1PtelwKUAgwYNSvItCCHi6sh6BB4P/3jlf/jz5ZGszArqm2oKGno3QMgE89fA1EoY8xU4zHbj/NJs1G0kGxgOB/Iwmo+I/H6Y1jqslGrY18FKqXzgWeAarfV3SsVriDR2jbOtTV5wrfVjwGNgpN3ed/GFEHF5PMYENavVmLW8cCEUF8cPEj4ftz5+EUsn67h/qfkN8PDLcPonrWoIdlkus7tJNjD8P2CjUupNjP8SJwB3K6XygFXtHaiUsmAEhSe01s9FNn+tlOofqS30B3ZGtlcBA6MOHwB8mWQZhRAdET2L2W+sjaAvu4xwr96YdAjVKkg88PID3DE5HP/rGxA2wfGf5FNQ10ADCktvGzkhGY7aHSW9UE/kBj4e47/FWq31Pm/YyqgaLAe+0VpfE7X9PmB3VOfzQVrrG5VSI4AVkescBrwODGmv81kW6hFiP5WVGfmOamrivqyBxl69ydEhdi1ZyOGfzieo4zcQmLFR8LcHOfqTYraHnTz4IMws9spw1Ax2wEt7Rm7wJwPf11rfrpQapJQar7Veu49DJwEXAZuUUhsj2xYA9wJPKaVcwGfAuQBa6w+VUk8BH2GMaLpSRiQJcQBa9x9EPd+V76RvfSBmslq0Xb3AW7gHZzVs+9/55F4MwdY7aTjFeSZ/OfdxuMLRqqtCAkJ3lWxT0h+ARuAk4HZgD0bz0Lj2DtJav0PCiicnJzjmLuCuJMslhEgkuv+goQF+/GN48UXIzSXkD3CddgNuHsVFCDO92dP8x+oZCa4ZYA1DwAT3vKwJBwOxQ1I15L11PbOPvA9HHpAnlYOeIqmmJKXUBq11sVKqXGs9JrLtA631qJSXcB+kKUmIWD4fVJX7GH1WEcrvT7hfHXaKqATAiZfRbOCWXtew4XtwQWl9y2gjjP7j//eqlRtPCWAKQ9AM976UyxMb3uZD+zhkuebu54CbkoCgUspEZISQUsqBUYMQQmSQpkrCsTle/ua30ofEgSGIBSde1jGOXThYN3Ibj89Q5DRqGlvdGcJhG6O+0FQubEl+l1eXw104m5OrSmDoOZKd4PYQ8FfgEKXUXcA7wN0pK5UQosOiBxltrm2buqI1C0G8OKGXD77/KsyYAxY/jbn1bRqAtU3xSP0i8urs/PDLAvLq7MzBzS4crdfiET3APmsMSqkc4FPgRox+AQWcpbWuSHHZhBDJiHQoV33rxGp1kOf34cTLfBayiGuwU9+mo08DD+VdwDfjl8Bxdxv5jMz1bU5tN+WBamThVDfX/raU15iJEy9enOzCgc0mo1F7on0GBq11o1Lqd1rricCWLiiTECJZUR3MowMB7m5wMRs3AaxYCXCD+UF+e72P/AfvgvqWG/9jxfDr0900mkg4PMRmsvG3859jTP8xOPIcFLjB5XKw1WLUEu5YAHPnSlDoiZLtfL4N+A/wnE524kMXkc5nkbV8PmPN5agOZk3sfT5ktWOuqoTVq2HWLCry6rl3EvxpDAkDQp4lj0bdiHuGm9KRsSuwdSR7hshsndH5fB1GGoyQUqoe47+U1loXdFIZhRAd1bTMZlRgaH2vN+eaYckSuPtu5p1Yz++PTbBjhN1s57mfttQSWnM4JCBkg6QCg9a6d6oLIoTooPz8mOahuPbsgdtv5++DgkZQSDSrSIPdlIt7hptTfnBKZ5dUdDPtBgalVHF7rzel1BZCdLGmvoWcyMBCmy1ukPD1git+FOSZke2cS8PN75q47v/KcQwalpryim5lXzWG30UebcAxwAcY3zmOxkiffXzqiiaEiNuoHz0utUk4TNCWj6V+r7FLL1gyFu48ARrMJMxbbAnDw//MYe61fwYJCiKi3cCgtZ4KoJRaCVyqtd4UeT4SuD71xRMii7VOie12G+sle71gjv3TbbTmEq4NYgGWFMPVpxupLNpLZH/uZvj9y+BotBprMAsRkewEt6FNQQFAa70ZY0U3IUQqRNcKamqMR5fL2L5hg9F3EEWHwtyYu4i7J5i47McQSFRLiHjxL/DUs5F1E6xWI9gIEZHsqKQKpdTjwF8wvm9ciLFMpxAiFZYsiW0qArBYoLwcrr02ZrOvF3x463z+8NIewlMTr5fQVFO47P0cztgeldFGpi6LVpINDBcDlwPzI8//BTySkhIJke18PrgrToLhoJH02pi85m/uR7jrBLAEHyJ8UoKFlTUQyoXy2fD+fL7dvZE6XASxkJ8bxCRTl0UryQ5XrQcWRn6EEKni88FLLxm1g9ajjBYsYNfAMfTyB/hjq36E+nCCoACoNTej/30d1Bk3/ycZxuuUMDTXy1/LnfQbJkFBxEp2oZ4hwD3AcIwRSgBorb+fonIJkX2aOpvN5jZ9CNhsMHcun3odzJ9QyrunLm23D6HJfSX3cfgPr8f1DmA3WqdsNqhVDq5wO+gnA5FEHMk2Jf0R+DVGjWEqRtNSEv8thRBJiTcEFSAvDxobmzPV5ft9vFvyRLt/ffnWfILhIItOW8TcY+bCJGPQkddrzInbu1dSWoj2JRsY7Frr15VSSmtdCfxGKfU2RrAQQuyPigpYuxbGjzfu1q3SW2iAQAD18MPGMFVgr8mLzWKlvvXayxoI2XEdvpC5ZxbjLHTGpLSQVBaiI5INDPWR9NufKKXmAV8Ah6SuWEL0cFddBYsXtzyfM8eYqxBlVy8o/16Q8MKrOeZHJ+AYNAxnoROtQs0jjADj97dvhveuY4V2cM85QB2UfSQ1A7F/kp3HcA3QC7gaGAtcBMxKUZmE6NkqKmKDAsDSpXD77ejcXDTGJLX+v4BTL4LTzw9w+LJReDZ7cOQ5WHSyG4J2qO9tjDZ68VF44x6oc2CxGCNdi4pg2jTj0eNJy7sU3VhSabczmaTdFt3O8uUwe3abzd89tIznvvsBm9+YzO8m6zb9CHazncprKnHkOVjyZx9X3+olsNPZPNoIwG4HrWMHNNntyJrMoo320m4nVWNQSq1WSr3R+qdziylED+PzQVmZUUMoKzOeg9Gn0MrOXnDCCyuZU39i3KAAkKNy8FZ7AZh7kYOqteO445cO7HYoKDACwIIFkJsbe1zTmsxCJCvZPobovEg24CdAqPOLI0QP4fGgXS4aw5AT8NOYaycnB1RTvqN585qbk1aMhNlnQdD0SrujjRp1I85CZ/NzhwNuucVYRa0pzx7A3a1WY5eJzaKj9rspSSn1ltZ6SieXp8OkKUlkHJ+P0IAizAF/m5eaV1RzOKh47++s+tcyrv/uWQKW9k9pybGw/OzlbVZUi6dpOoTFYgSFplgkRLQDXsFNKXVQ1NMcjA7o73VC2YTocb4t95IbiP/HVRuw0Fju5dbQ7Swui3RAJ/or1EDYAmt+ARuugx84oL11FSJKS1vmLcioJLE/km1KWk/LcrIh4FPAlapCCZFRklnoOGqfqup8RtK2tgBgIcgrdQEWfxA1KilBauzct6+m4b1boM5BEKMWUFKS3I1e5i2IA5FsrqTBqS6IEBkp0ZoI7ewzpOSMmJcrDob3B8DoKgsPfOdmqHVby9esJlEtuqYwjH7JxfoNi2LOk5NjxB654YtUS7YpaWZ7r2utn+uc4giRQaLTVDTNSG79tT3OPrYXnwaMgHDVj+D1I5pOGOQo/xpOCFwR93JDn5+P7btD+Oyrs1lf1zaJUWOjdCKLrpFsU5ILOA5oGqI6FXgTqMH4riOBQfQ8Xm+bNBUxYz+9Xvj227b7AK4zYOnYyJOomsEm+2L2fHcFvD8Pjo1qTnp/Hls2PhhzmXDYCAZgXEKyY4uukmxg0MBwrfUOAKVUf+D3WuuLU1YyIdLN6WyTpoJgEDZsQE+ZQliZMYUaUK1G9t0/IRIUEgw9/ap+FdY3HiZQdgUMWAtV42H3MPLzIRSCX/3KGIIKxro8AGPGSFAQXSfZwOBsCgoRXwM/TEF5hMgcDofxNT167OfChYSuvhZzwN/8x6OBXX3MeA/JJb8uxIJTg6AaE5525iGHMnoZuFzDyNk2jMZGWPgoFBe37d8+5ZQUvj8hEkg2MLyplPon4MH4OzgfWJ2yUgnRBZIZbNR67Oe35V5MATMFUbusHAmuGSGsNiv1gA5bgDjZT4ErynI49sqpHOuQIaUicyU7KmmeUups4ITIpse01n9NXbGESK12Bxu1jhhRYz+9wDAC+HpB+fegOhdcM8BvAX/TKmpxRhwd51U89pqVEff/sflcMqRUZKpkawxEAkHcYKCUeldrPbHTSiVECrU72GhV4ojh88FTqx2sLZ7OW9OfI5wo01jQBjkaQmYwB+i1agGmjdMJPO+EUyQSiMyXdGDYB9u+dxEiMyQabFRV7sORIGJ4VjmYM8/H4NEuKn78YvvrFyoFj2yA3L1Q7aSuzsE6OwwYk8p3JUTn6azA0L1zd4uskmiwkRNv25215pO15Vz0+FrCl99Jhbkh4UxlgnmgGjnF72b1d8MIBo2XZKip6G46KzAI0W00DTa6aY6PYlVOXqiaS+YX0jcvr818BM8P6pm97kzCkxMEhCYhG6x8Dr4aw9vawQcfwOefGy/JUFPR3SQ783ke8ITW+ttEu3RekYTogKSGFrVViofzw7MhaFQd1L3A7yxgsVBREGTtAPjBLqNjOUA7QaEp0d3zS+G/xthSS4GxhLMMNRXdVbJLe34PKFNKPaWUOk0p1frP5KJ4BymlliqldiqlNkdtO0gp9ZpS6pPIY9+o136plNqmlPpYKXVqh9+NyC4ez/6tYRnpfVbBAIqoe34wyFUnBxk+z1gf4fhLINxOQLCEYOSaH8MDX8DmlvxJsv6B6O6SCgxa61uAIYAbmA18opS6Wyl1ROT1zQkOXQac1mrbzcDrWushwOuR5yilhmPMjxgROeYPSilTR96MyCLRQ4tqaoxHl6tllTRaFlCL2mTweo2MdBj5jB4eB08Ng38fDouPheZooSDQ+n+gBlsQ7njbzEfDnmTpXS/w6AOxq6hJf4Lo7joyXFUrpb4CvsJIvd0XeEYp9ZrW+sYEx/xLKeVstXkGcGLk9+UYOZduimxfqbVuAD5VSm0DxgPvJv1uRPaIN7TIbG5OP+rxwJw5YDIZOYeWLo2ap+B0Qjgcm88IUImGUIQsFISCBEzwq3/B3PXg0BZYORUcMG4czJwpk9VEz5FsH8PVwCxgF/A4cIPWOqiUygE+AeIGhgQObUqvobXeoZQ6JLL9cOC9qP2qItviledS4FKAQYMGdeDSoseIN7Rozx7YsAGfcxyzZkGfoA8nXrw4mTXL0ZIU1eHgrutOYKnl1Zi+A62IO77O/MRqTtCrWbzzLvprK1aCbaoFMllN9CTJ1hj6ATO11pXRG7XWjUqpMxIc01GJBgG23aj1Y8BjYCzt2UnXF92JwwELF6Ivuyz2P86117K570xmB5/jIeYTwkQOmgssC1n5djHnn+oE4DfWBBldGk2QE255/v483vzzJKzWSeTlz8W61yvVAtHjJZsS49Z2Xqvo4DW/Vkr1j9QW+gM7I9urgIFR+w0AvuzguUUW+aevmIn0poA9LRstFo744wKW8HhzwFhSDK+cfhmvbcrjpopGFkxeQK7ZSigUbHvSRgssfQv6bYOq8dhqh2G93WguAkfkR4ieTWmd2i/ckT6Gv2utR0ae3wfs1lrfq5S6GThIa32jUmoEsAKjX+EwjI7pIVrrcIJTA0aNYd26dSl9DyLz+HwwdpCPLfVF9IpaRlMrBVqzqxesLoJnh8FTRxFTH7WZbCil8Iei+ic0KJWD+fm/ECxvGWFkt0NlpVQQRM+jlFqvtT4m3mvJDlfd3wt7MDqPj1RKVSmlXMC9wDSl1CfAtMhztNYfAk8BHwGvAFfuKyiILBFneJHXC9UWB3Nw04AZTaTdUWt+NwG+dz2c99O2QQHAnGNiweQF2M128sx5mDBzyVFX8/X1X7H8hlIZYSSyXkpnPmutSxO8dHKC/e8C7kpdiUS3kyAN6oYNkLvHx7G8h5UQu3qBtxAeHA8rRtHulMtgKMDcsXOZO3Yu3movzkInjjzj7t8qy7YEBZGVJCWGyFwJ0qDuGl3CmqtW8QWz2HZwkNIp8OwIsIagzkr78/A1LJpyd3MgoM6B9yPA2RIEZISRyHYpbUoSoj0JJ6BFfFvuJZRjjd1osbB7VTmPBl1cfkaQEfPgyaMgZIK6XNpNXWENwbEvzWHmiOuB/Z84LURPJ4FBpEXcm3JUpPB4YNQMJ4HalrkKvl7wfl8//sBnPDAh2LKu8j5qCGgo3Hwy+Q98RMXHbrzepCZOC5G1pClJdLl4LUSvzPJwvtmFslrRgQAvh9x8HixlDm5+22sO7rFB7jshjDUcJPDtzwlOo/3Edk28k+HvSzDv7ocTL1/bfDidjoRrMkQmTguR1SQwiJSKl/y06aac5zdmJu8hn0eCLlTQiBQKWIKLf/QazZMTPuDJ48JgCoOCekvkxIlGWWvg7eth98jmeQjnBD0swUUQK3mNAcyr3FBSGn9NBmfnfwZCdDcSGETKJFpX2emEGXUeHsFFACs26mls1ar515F+vj1rJJgak0vq3hQo1rvgjfuaNztsPpZbXOSE/YAfAoDLhaOyBLfbgctl1BSCbbNcCJG1Uj7BLdVkgltm8vmMvoPopprmyWL4CA0owhyImpgG7OoF5d+D6lyYNTOqdhCHJQQmqw1Tjgl/qIHGzTPJK/sNgS+HoRTYbMbN/q8Lyph69zSs/prmYwP2AqxvrYJx4/Z3OQchur32JrhJjUEkp4N30Hbb8PFitlsh4KfiYFj1fdjmgEeKIbivROsabCG4yT+PK2++tXkeAnWO5uI1Xd/pBLXLSeh/A0SPbQr5g3yX76QfMjRViHgkMIh9S9Qm1I6E6yo7YUtVPtV963l4amQyWpN9NRlpOHzNz3h49s2cfcEwgJb5CHmxN/im38u8Dh6xu1nsdxHEgoUgV9rcXLHXQb99vW8hspQ0JYn2tdsm1P5X7aZ4Et2Gv+Y7D4s/d5FDkEZzaN9DTSNMYTj8pZvxVdzTodxFTcVv6uj24qTW7pD8RyLrSVOSSGifLUQHMK4zOr1EwOzjrU/KWfz5HLDU05hE2XJD8KfnoLABRn8FtwWcHP+njt3QHQ4jILlcDrZaHNLJLEQSJDBksaRaiNprE0qGv4LFLz3En4LLQJvAUt/+/pFaggrn8MfnG/lpVFL3xfZrUSUz6Wjqa8l/JETHSGDIUgnSELWsctak5St3h8d1epZcxZzPF1NvBtoZYQSABpOGAe/8lEqvi6O/gpnhcyBqrQV1ADPQpJNZiORJSows1dRCFK3pvttGaanRp7BqlfHYqlpR4atg+cblVPhavt77PqvA9fliY8hp634EDblB49EaAHMIrn4Xvrwffv7G9+G/p1CjxmAmFHuczEAToktIjSFLdbiFKMFX7qteuorFZYubn88bP4+Hf/Qw3o/XYm0kagmdFvYQ/M0DA2tgby44q8FRZ7x2LYt4kOv4utFB7SI3BdfKDDQhuprUGLJUUwtRsovSxMuEWuGriAkKAIvXLqbCV4HzyPEEWv/v0mAPgvt5KPk6nyN2mzjmy5agABDEwtBcL243FMxtv6YihEgNCQxZbB8tRM08Hhg01MfUC8sYNNTXnJ567Rdr2+Ys0sZ2x6BhuAfOwx6EgnqwBeGON6ByIZRW9mbvPYspNm3Cjz3m8PzcIH8td7aUxeEwFlyWmoIQXUaakrJcUwtRU40getSOr9ZH+adefvbQBkKXXwthK5gCzLrPTUlJKUPqfxD3nE3bS+c+TEn5+XhnTsW5M9hSM7CH2H7k6VTlO5hT42YpLZPPdtzi5ohhEgSESCcJDCLusFVGenC94EJpM6FT9xgdyBajxyB4uovyj0sYsN3Kz9+38H/HBpvP9fP3zRQOtsIE47ljzCQcdy83LlDQ0lcwYIyDQACepJTXKYmkxHayfq4EBSHSTWY+Z7nmic3KB4VeqHZis4G6rgh/KF7XMVBfwD//ZxXFJie9hhdRebCftQNgfBUU7bZT91El/Vp/648zky7ezGjpRhCia8jMZ5FQ+cc+Go9fAuPvgnAumAI0vr8AFbISf0wRWO1Bxgx2ouscLFUurty9mKG7jdd+r1yc1y/Ot/44o5pk4pkQmUk6n7OYZ5OHGW8W0XDc/xozkm01YPETOPZugsG62J012HN6YzfbWTbTjSPPQVW5jzna3by6pgLmaDdV5cmvjyl9y0JkHgkM3UC8oaIHur+v1ofrBRf1YX+bCWg2s5lb/tXYPKLIHoSHXjTz2IhnqLymktKRRnuPEy9BYmfJBbHgxNuBdyeEyDTSlJThOprx2uOBOfN8mPp5CdXmc8sde5n7U2dLeuoIb7UXq8katx9BmQL8bH0v5q3fg7fQmIBmretFY35f+kadp+8YJyFrwFgVLSLPGsQ8xnlgb1oIkVZSY8hg0fmMamqMR5crcU3A54NZ93uov7yI2nOm0HDJcP53+xSKHizCs9kTs6+z0Ekg3GrqszbmG7h7/w9FoRCOOhgXmYCWZw3St/UN3+HAvMyNttsJ5RWg7XbMy2R2shDdnQSGDNahfEYYHcnBH7mMYaVWf/MQU3/Ij+t5F77alojiyHPgnuHGbrY3Nxfd8Tp8thBKF3gwP7QwuRt+aSmqshLz6lUomZ0sRI8gTUkZrMP5jAq9xiQ0S9vmIYvJgrfaG9OkVDqylJJKE975s3DuqG+ZgFZggeJi44afzJAhSV0qRI8igSGDNeUzmjPPh+lgL+HdTtyLHQnvwWMGO7HaAwTiTE0JhoPG2sjRPB4cLhcOf6s1Epqij9zwhchK0pSU6UZ6UNcVwaxpxuNRnoS7OvIcLJtpNA+ZsBl5jIJ2CNpx9XPHdkBHd2BEs9kki6kQWU5mPmcwX62PogdjZyDbzXYqr6lsM8ooWsVnPsac6KXhu3zI3QvVTuy61TrHZWUwbZrRq90kLw+eew5OOSVF70gIkSlk5nM3FW9Iaby+gtb2fu3A9o2Dhqh7vqWg1eJn8TowGhthzJhOK78QonuSpqQu4qv1UfZFWczIoH2JN6Q0bl9B6+OcSXRad3RBBiFE1pDA0AU8mzwUPVjE1GXT4s4pSCRmSGluAXazHfcMd7u1BejAPT/ZBRmEEFlF+hhSzFfrY8DvigjoluYgq7JT9Yv2+wlan8Nb7cVZ2HYGc3t2VfjYudbLIeOdbbOdCiGyWnt9DFJjSLHyT70E/LGz1AJ+C+WfepM+hyPPwbjDxyUfFHw+uPNO+o0tYvj8afQbW0TzsmtCCLEP0vmcatVOMLVq8DcFje2p4PHAnDlQH5mb0DQc1eUyclxLH4IQYh+kxpBiY450YHnJbcwnqC+AoB3LS27GHNmBG3Sy6VWb5ibU17d9rb1cGkIIESVrA8P+jBLaHw4HLL+hFNsjleQ9uwrbI5Usv6E0+S/uHo+xxNq0acZje01C8ZIrNWk3l4YQQrTIuM5npdRpwCLABDyutb63vf33p/PZs8lYz9hqshIIB3DPcDevMZAqcVa2TO6goqLY2cl2O7Ez1faxf9Mxsm6mECJKt+l8VkqZgN8DPwKGA6VKqeGdeY2mBWr8IT81DTVxM4+mwn6tVNbR9Kqtx6nabHDHHTIUVQjRIZnW+Twe2Ka1/i+AUmolMAP4qLMusL+zidOiw+lVkYWUhRAHLKNqDMDhwOdRz6si22IopS5VSq1TSq3zJbveZcT+ziZOi/2dnSwLKQshDkCmBQYVZ1ubThCt9WNa62O01sc4Onjz29/ZxGkjs5OFEF0s05qSqoCBUc8HAF929kVKR5ZSMrhkv2YTp4WsiyCE6EKZFhjKgCFKqcHAF8D5wAWpuJAjz5H5AUEIIdIgowKD1jqklJoH/BNjuOpSrfWHaS6WEEJklYwKDABa65eAl1J9nf2aV9CtLiiEEPsn0zqfu0RHJhN3zwsKIcT+y7iZzx3V0ZnPHZ1MfMC6/IJCCLFv3Wbmc1fo6GTi7ndBIYQ4MFkXGPZnMnH3uqAQQhyYrAsMXb7UsaytLIToZrKuj6GJjEoSQmSz9voYMm64alfp8snEMntZCNFNZF1TkhBCiPZJYBBCCBFDAoMQQogYEhiEEELEkMAghBAihgQGIYQQMSQwCCGEiCGBQQghRAwJDEIIIWJIYBBCCBEjewODzwdlZcajEEKIZtkZGGRFNSGESCj7AoPPBy6XsaJaTY3x6HJJzUEIISKyLzDIimpCCNGu7AsMsqKaEEK0K/sCg6yoJoQQ7crOhXpKS6GkRFZUE0KIOLIzMICsqCaEEAlkX1OSEEKIdklgEEIIEUMCgxBCiBgSGIQQQsSQwCCEECKG0lqnuwwHRCnlAyrTXIx+wK40lyFTyGfRQj6LFvJZtMiUz6JIax13aGa3DwyZQCm1Tmt9TLrLkQnks2ghn0UL+SxadIfPQpqShBBCxJDAIIQQIoYEhs7xWLoLkEHks2ghn0UL+SxaZPxnIX0MQgghYkiNQQghRAwJDEIIIWJIYDgASqmBSqnVSqkKpdSHSqn56S5TuimlTEqpcqXU39NdlnRSShUqpZ5RSm2J/P+YmO4ypYtS6trI38dmpZRHKWVLd5m6ilJqqVJqp1Jqc9S2g5RSrymlPok89k1nGeORwHBgQsAvtNbDgAnAlUqp4WkuU7rNByrSXYgMsAh4RWs9FBhFln4mSqnDgauBY7TWIwETcH56S9WllgGntdp2M/C61noI8HrkeUaRwHAAtNY7tNYbIr/vwfjjPzy9pUofpdQAYDrweLrLkk5KqQLgBMANoLUOaK2r01qo9DIDdqWUGegFfJnm8nQZrfW/gG9abZ4BLI/8vhw4qyvLlAwJDJ1EKeUExgDvp7ko6fQgcCPQmOZypNv3AR/wx0iz2uNKqbx0FyodtNZfAPcDnwE7gBqt9avpLVXaHaq13gHGl0vgkDSXpw0JDJ1AKZUPPAtco7X+Lt3lSQel1BnATq31+nSXJQOYgWLgEa31GKCWDGwu6AqR9vMZwGDgMCBPKXVheksl9kUCwwFSSlkwgsITWuvn0l2eNJoEnKmU8gIrgZOUUn9Jb5HSpgqo0lo31R6fwQgU2agE+FRr7dNaB4HngOPSXKZ0+1op1R8g8rgzzeVpQwLDAVBKKYx25Aqt9QPpLk86aa1/qbUeoLV2YnQuvqG1zspvhlrrr4DPlVJHRjadDHyUxiKl02fABKVUr8jfy8lkaUd8lBeAWZHfZwHPp7EscZnTXYBubhJwEbBJKbUxsm2B1vql9BVJZIirgCeUUlbgv8DFaS5PWmit31dKPQNswBjFV043SAnRWZRSHuBEoJ9Sqgr4NXAv8JRSyoUROM9NXwnjk5QYQgghYkhTkhBCiBgSGIQQQsSQwCCEECKGBAYhhBAxJDAIIYSIIYFBCCFEDAkMQgghYkhgEKKTKaXGKaX+o5SyKaXyImsRjEx3uYRIlkxwEyIFlFJ3AjbAjpE36Z40F0mIpElgECIFIqkwyoB64DitdTjNRRIiadKUJERqHATkA70xag5CdBtSYxAiBZRSL2CkHx8M9Ndaz0tzkYRImmRXFaKTKaV+BoS01iuUUiZgjVLqJK31G+kumxDJkBqDEEKIGNLHIIQQIoYEBiGEEDEkMAghhIghgUEIIUQMCQxCCCFiSGAQQggRQwKDEEKIGP8fuukhm9jYTncAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "predictor = Predictor(graph=trainer(), data={graph.x: data[\"x\"]})\n", "\n", "y_prediction = predictor(graph.y)\n", "\n", "data[\"y_quadreg\"] = y_prediction\n", "\n", "ax = data.plot.scatter(\"x\", \"y\", color='blue', label=\"true\");\n", "ax = data.plot.scatter(\"x\", \"y_linreg\", color='red', label=\"lin.reg.pred.\", ax=ax);\n", "ax = data.plot.scatter(\"x\", \"y_quadreg\", color='green', label=\"quad.reg.pred.\", ax=ax);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## More on `connect_via_regression`" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The first argument `name_prefix` is used to name the entities in the graph holding the regression parameters:" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "graph\n", "├─regression_y\n", "│ └─location\n", "│ ├─slope\n", "│ └─intercept\n", "├─inputs\n", "├─outputs\n", "├─x\n", "└─y\n" ] } ], "source": [ "print_child_tree(graph)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Printing the graph's children reveals that besides regression parameters for the mean of `y` in `regression_location_y`, there are also regression parameters for the variance of `y` in `regression_log_scale_y`. Thus, the mean of `y` as well as the residual scatter of `y` is learned as a quadratic function of `x`.\n", "\n", "The argument `inputs` specifies the input variables for the regression. This can be either a single variable or a list of variables.\n", "\n", "The argument `outputs` specifies the output variables of the regression, either as a single variable, or a list of variables.\n", "\n", "The argument `order` specifies the order of the regression, i.e. the highest power of the input variables in the regression polynomial. For example, `order=1` yields linear regression.\n", "\n", "The argument `include_cross_terms` specifies whether to include cross terms when `order>1`. Cross terms are not enabled by default, since this can significantly increase the number of model parameters and thereby make the model hard to train without overfitting.\n", "\n", "The arguments `inputs_location=x_data_mean`, `inputs_scale=x_data_std`, `outputs_location=y_data_mean`, and `outputs_scale=y_data_std` allow one to directly include scaling of the data for standardization in the regression. Standardization of the data would otherwise be a required step of data preparation. Here, one just needs to provide the empirical location and scale parameters of the data.\n", "\n", "More on these and a numbe of further arguments can be found in the documentation for `connect_via_regression`.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.7" } }, "nbformat": 4, "nbformat_minor": 4 }