{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Causal Structures - Prediction"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"%%capture\n",
"# execute the creation & training notebook first\n",
"%run \"02-01-creation_and_training.ipynb\""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In the [previous section](./02-01-creation_and_training.ipynb) we created and trained our causal structure.\n",
"We can now make predictions. Let's create some test data first."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"test_data_a = {\"(a)\": np.linspace(4.5, 5.5, 100)}\n",
"prediction = causal_structure.predict(data=test_data_a)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Since we passed a dictionary as data to the predict method, the result is also a dictionary."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"pandas.core.frame.DataFrame"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"Index(['(a)', '(b|a)', '(c|a,b)'], dtype='object')"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"display(type(prediction))\n",
"display(prediction.keys())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If we pass a pandas data frame we get a pandas data frame in return."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
(a)
\n",
"
(b|a)
\n",
"
(c|a,b)
\n",
"
\n",
" \n",
" \n",
"
\n",
"
0
\n",
"
4.500000
\n",
"
-3.950142
\n",
"
45.574075
\n",
"
\n",
"
\n",
"
1
\n",
"
4.510101
\n",
"
-4.436720
\n",
"
45.443720
\n",
"
\n",
"
\n",
"
2
\n",
"
4.520202
\n",
"
-4.927574
\n",
"
45.272895
\n",
"
\n",
"
\n",
"
3
\n",
"
4.530303
\n",
"
-5.408997
\n",
"
45.130555
\n",
"
\n",
"
\n",
"
4
\n",
"
4.540404
\n",
"
-5.907117
\n",
"
44.996075
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" (a) (b|a) (c|a,b)\n",
"0 4.500000 -3.950142 45.574075\n",
"1 4.510101 -4.436720 45.443720\n",
"2 4.520202 -4.927574 45.272895\n",
"3 4.530303 -5.408997 45.130555\n",
"4 4.540404 -5.907117 44.996075"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"test_data_a = pd.DataFrame(data=test_data_a)\n",
"prediction = causal_structure.predict(data=test_data_a)\n",
"prediction.head()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Note that even though the column '(c|a,b)' depends on both '(a)' and '(b|a)' and we only provided data for '(a)' we get a prediction for '(c|a,b)'. This prediction is based on averaging the possible values of '(b|a)' given the provided data for '(a)'."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can plot the prediction and compare it to the training data."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Text(0, 0.5, '(c|a,b)')"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAmcAAADrCAYAAAAsXzGbAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABL1klEQVR4nO3deXxU9dX48c+ZIKuSWECgAuLGr1I1aQkZ0ArutC7gvmuf8vigrWQmSa0bSkLVuhMSfNxqF1sXrFoVreL64MYmaEDcqtZdQVACiqIkc35/3LnjncnMZBImuTPkvF+v+2KWe28OCMdzv6uoKsYYY4wxJjcE/A7AGGOMMcZ8z4ozY4wxxpgcYsWZMcYYY0wOseLMGGOMMSaHWHFmjDHGGJNDrDgzxhhjjMkh3fwOIFv69++vw4cP9zsMY0wnWbZs2VpVHeB3HNlg+cuYriddDttqirPhw4ezdOlSv8MwxnQSEXnf7xiyxfKXMV1Puhxm3ZrGGGOMMTnEijNjjDHGmBxixZkxxhhjTA6x4swYY4wxJodYcWaMMcYYk0O6XnH2zjswZQp8+63fkRhjTNt8+y2Ul8Mrr/gdiTGmA3W94uyFF+CPf4RjjoFNm/yOxhhjMrduHdx3Hxx1FHzxhd/RGGM6SNcrzs44A26+GR55BI4+2go0Y0z+GDQI/vlP+OgjOOkkaGryOyJjTAfIueJMRGpE5GMRaYgeh2X9h0yZ4rSezZvnPIFagWaMyRdjxsCNN8ITT8DFF/sdjTGmA+RccRZVq6ol0eORDvkJZ54Jt94Kjz9uBZoxJr9Mnuw8ZF5zDSxb5nc0xpgsy9XirHP893/HF2jffON3RMYYk5mrr4YddoCzz4bmZr+jMcZkUa4WZ1NFZIWI/FlEtk91kohMEZGlIrJ0zZo17ftJkyd/X6BNnAhff93emI0xpvMUFkJtLSxdCjfd5Hc0xpgs8qU4E5EnRWRlkmMScCOwK1ACfApcl+o+qnqLqpaqaumAAUk3ds/M5Mnw17/C00/D4YfDxo3tv5cxxnSWE0+Egw+Giy6CDz7wOxpjTJZ08+OHqurBmZwnIn8EHu7gcBxnnAHdusHpp8ORR8LDD0Pv3p3yo40xpl1EnFazn/wETjkF5s938pgxJq/lXLemiAz2vD0aWNlpP/yUU+Bvf3MS3KRJNgbNGJP7dt3VWR7ohRegutrvaIwxWZBzxRlwtYi8IiIrgAOAyk796aee6nRxPvUUHHYYfPllp/54Y4xps5NPdiY4XXGFM37WGJPXcq44U9XTVXUvVd1bVSeq6qedHsQZZ8Dtt8Nzz8GBB8LatZ0egjEmP4hIgYi8LCIPR993/FqNydTXw49/7PQA2PgzY/JazhVnOeOUU+D++5097A44AD77zO+IjDG5KQy8nvBZx6/VmKh3b2drp82b4dhjbe1GY/KYFWfpHHmks83TO+/AQQdZgWaMiSMiQ4DDgVv9jgWAESOccbNLl8LUqaDqd0TGmHaw4qw1Bx7ozNx85x3n9erVfkdkjMkds4DzgEjC562u1ZiVdRqTmTQJpk2DP/3J6eo0xuQdK84y4RZo774L48fDxx/7HZExxmcicgTwmaom7p+U0VqNWVunMZnf/97Z9aSqytlD2BiTV6w4y9SBBzpJ7pNPYNw4eO89vyMyxvhrX2CiiLwHzAEOFJHbVXW1qjaragT4I1DW6ZEFAvD3v8NeezkL1f7nP50egjGm/aw4a4v99oMnnoAvvnBev/mm3xEZY3yiqheq6hBVHQ6cBDytqqf5ulaj17bbwoMPOgvVnnYaNDX5EoYxpu2sOGurYNBZpPa775wCraHB74iMMbnF37UavXbaCW68ERYuhD/8wbcwjDFtY8VZexQXO2ug9ezpLLOxYIHfERljfKSq81X1iOhr/9dq9Dr5ZGdx7d//HhYt8jUUY0xmrDhrrxEj4PnnYcAAOOQQePJJvyMyxpjk/vd/YcgQp3vzq6/8jsYY0worzrbEsGFOC9puu8HhhzvjO4wxJtcUFjrrn/3nP84MTmNMTrPibEsNHOiMQfvJT5xVue+6y++IjDGmpXHj4Pzz4Y9/tAdJY3KcFWfZsP32zizOffd1xnb88Y9+R2SMMS3NmAE//Sn88pfweuKOU8aYXGHFWbZstx08+ihMmABTpsDVV/sdkTHGxOveHf75T2cy02GH2Y4nxuQoK86yqXdvp7vgxBOd7oOLLrK97YwxuWWnneChh5zCbNIk2yDdmBxkxVm2de8Od9zhtJ5dcYWzx50VaMaYXDJ6NNx+OyxeDBUVfkdjjEngW3EmIseLyKsiEhGR0oTvLhSRt0XkTRGZ4FeM7VZQ4Cz86BZoF19sBZoxJrccc4zTwn/zzc5WT8aYnNHNx5+9EjgGuNn7oYiMxNkK5cfAD4EnRWSEqjZ3fohbIBBwCjRVZ2Xuxkaor3cKN2OMyQWXXeYsTHvWWVBS4uzFaYzxnW8tZ6r6uqom25xyEjBHVb9V1XeBt/Fj4+BsCATgppvgvPPghhvghBNsfIcxJnd06wZz5jjroB1zjPMQaYzxXS6OOdsR+NDz/qPoZy2IyBQRWSoiS9esWdMpwbVZIABXXQW1tc4sqcMOgy+/9DsqY4xxDBoE99wD770HZ5wBkYjfERnT5XVocSYiT4rIyiTHpHSXJfks6YAtVb1FVUtVtXTAgAHZCbqjVFQ4A3CffRYOOgg+/9zviIwxxvGzn8F11zmzOG2DdGN816FjzlT14HZc9hEw1PN+CPBJdiLy2amnQt++cPzxzmrdTzwBP/yh31EZYwyUl8OSJXDJJc7Ys0npnqGNMR0pF7s15wIniUgPEdkZ2B1Y4nNM2XPkkTBvHnzwgfO0+p//+B2RMWYLiEiBiLwsIg8nfH6uiKiI9PcrtjYRcXY3GT3aeZBcscLviIzpsvxcSuNoEfkIGAv8S0QeA1DVV4F/AK8B84Bz8m6mZmv23x+efhrWr3cKtOXL/Y7IGNN+YSBuLyQRGQocAnzgS0Tt1asXPPCAM0Fg4kRYt87viIzpkvycrXm/qg5R1R6qOlBVJ3i+u1xVd1XV/6eqj/oVY4caPdoZf1ZQ4HRxPv203xEZY9pIRIYAhwO3JnxVC5xHivGyOe2HP4T774ePP3aW2LA1Go3pdLnYrdl1/PjHsGABDB0KP/853Hef3xEZY9pmFk4RFpviKCITgY9VNW2TeE7PNi8rg0svdWZx/vWvfkdjTJdjxZnfhg6F5593kuGJJ8Kdd/odkTEmAyJyBPCZqi7zfNYbmAZMb+36nJ9t/rvfOUMwysvhrbf8jsaYLsWKs1xQVORMEthvPzjtNPjTn/yOyBjTun2BiSLyHjAHOBD4O7AzsDz6+RDgJREZ5FeQ7VZQ4Gzr1KMHHHssbNzod0TGdBlWnOWKbbeFf/0LDj0UzjzT6VKwsR7G5CxVvTA6bnY4zpZzT6vqsaq6g6oOj37+EfBTVV3lZ6ztNmQI3HUXrFwJ//M/lpOM6SRWnOWS3r1h7lxnle7p053BuE1NfkdljOnKDj0ULr/cKdJqa/2Oxpguwc+Nz00y3bs7A3CHDHFW6l692kmKvXv7HZkxJgVVnQ/MT/L58M6OpUNccAEsXQrnngs77wxHH+13RMZs1azlLBeJOE+q11/vbKdyyCHwxRd+R2WM6apEnPFnwSCccoozy9wY02GsOMtl55wD//iH88S6777w/vt+R2SM6ap693YeFocOdXY6sd1NjOkwVpzluuOOg8cfh08/hbFjoaHB74iMMV1V//7wyCPOxICjjrIZnMZ0ECvO8sH48c5aaAUFzuv58/2OyBjTVe22mzMO9tVXYfJkm8FpTAew4ixf7LmnM85jyBCYMAHuvdfviIwxXdWECc6EpX/8A6680u9ojNnqWHGWT4YOheeeg9JSOOEEmD3b74iMMV3VeefBySfDRRc5e3EaY7LGirN884MfwBNPwKRJEAo5U9sjkbhTNKGbIfG9McZsMRFnN5Ng0NnZ5OWX/Y7ImK2GFWf5qHdvp1tz6lS47jo49VT49lsAampqqKysjBVkqkplZSU1NTU+BmyM2Sr16gUPPAD9+jlrn61b53dExmwVfCvOROR4EXlVRCIiUur5fLiIfCMiDdHjJr9izGkFBVBfD1ddBXPmwC9+gTY20tjYSF1dXaxAq6yspK6ujsbGRmtBM8Zk36BBzsPiJ5/YBAFjssTPHQJWAscANyf57h1VLenccPKQiDPuY/BgmDwZGT+e2kcfBaCuro66ujoAwuEwtbW1iIif0RpjtlZlZc6DYlWVMxY2FPI7ImPymm8tZ6r6uqq+6dfP36qcfrqzafo77yD77kvtr38d97UVZsaYDldR4SxOe+658OyzfkdjTF7L1TFnO4vIyyLyjIjs53cweeHQQ2H+fHTjRjb+5CcEPV95x6AZY0yHEIHbbnP23jzmGHjnHb8jMiZvdWhxJiJPisjKJMekNJd9CgxT1Z8AVcCdItI3xf2niMhSEVm6Zs2ajvgt5BUdNYrLDzuMVd98w7PduhH55z8Jh8NxY9CMMabDbL+904qvCkccYRMEjGmnDi3OVPVgVd0zyfFgmmu+VdXPo6+XAe8AI1Kce4uqlqpq6YABAzrmN5FHRISm4cP5y//8D9uMGoUceyy1Q4cSDoUoKiqyrk1jOoCIFERb+h+Ovr9URFZEJzQ9LiI/9DvGTrXbbs66Z++84yz58803fkdkTN7JuW5NERkgIgXR17sAuwO2w26GampquOzmm5Gnn4ajj0bOPZfajRupuegiv0MzZmsVBl73vL9GVfeOTmp6GJjuS1R+GjcObr/d2XbupJOgqcnviIzJK34upXG0iHwEjAX+JSKPRb8aB6wQkeXAvcDZqvqFX3HmIxFx1kK75x6YNg3505+c7Vasi8GYrBKRIcDhwK3uZ6q6wXNKH6Brjic44QS4/nqYOxfOOsuW2DCmDXxbSkNV7wda7PmhqvcB93V+RFuhQAAuuwx+9CNn/aF994VHHoHhw/2OzJitxSzgPGA774cicjlwBrAeOCDZhSIyBZgCMGzYsA4N0je/+Q2sWgWXXgq77ALTpvkdkTF5Iee6NU0HOO00ePxx+PRTGDMGlizxOyJj8p6IHAF8Fh0bG0dVp6nqUOAOYGqy67vMmNkZM5zlfi6+GO64w+9ojMkLVpx1FfvvDy+84Gy3sv/+cF/rjZO2R6cxae0LTBSR94A5wIEicnvCOXcCx3Z2YDlFBP74Rxg/Hn71K3jmGb8jMibnWXHWlYwcCYsXQ3ExHHcczJyZ8lTvHp3u4e7RaUWaMaCqF6rqEFUdDpwEPK2qp4nI7p7TJgJv+BJgLunRw5nBueuucNRR8PrrrV5iTFdmxVlXs8MO8PTTTnH2299CZSVEInGnqGpsj86xY8dSUVFBRUUFdXV1rFu3joqKCttI3ZjUroyu57gCOBRnNqfZfntnzGv37nDYYbB6td8RGZOz/Nxb0/ilVy+4+25nH7xZs+CDD+Bvf4M+fQBntmdtbS2qSn19PYsXLwYgFN0vr76+nnA4jKra2mnGAKo6H5gffd21uzHT2XlnZ5HaceOcXQSeftppVTPGxJGtpYuqtLRUly5d6ncY+WfWLKcFraQEHnwQhgyJfaWqBAItG1dtI3WTC0RkmaqW+h1HNnS5/HX33c76Z2eeCbfc4oxL80h88LMHQbM1SpfDrFuzq6uocNYheustGD0aXbgQIDbGLBkrzIwxW+TEE+Gii+DWW5210Dy8412BuPGuxnQVVpwZOPxwWLiQL779lub99iPyl79QWVlJXV0dJSUllJWVxZ1eUVGRdFKAze40xmTs0kth4kQIh+Guu4D48a5ugebmosbGRssppsuwMWcGAB05kmuPP56DbrmFgyZP5th99uHZ4mJebmgAoLy8HBFh8eLF1NfXx8aluS1oNTU1NDY2xj5zk2pRUZE98RpjWgoEYM4c+PnP4YwzoG9f5PDDqa2tBaCuro66ujrAhlKYrsdazgzgTAK4/Kab+Fd5OTcC+y1YwPnLl7Pf6NGEQiHq6uqYNWsWCxYsIBwOx22kbk+7xph26dULHnro++V9FiyIPfh5WWFmuhprOTMxIsJ1dXUEZs/mHeBaoOzDDxn2wAOICJFIhKqqKvr27Rtb70xE4pKpPe0aY9qkb1949FFne7kjj0Sff57Km2+OO6WystJyielSrDgzMd5JANcB7wJ/X7WK1cOH88Dkyfxx8WIaGhooKSkBYMOGDbFuS7dAcwszsKddY0yGBgyAefPQsWNZFwxy95dfxh7u3FZ4sJxiug4rzgxAXFdkOBymsLCQBx98kJ8tX85Dmzdzxs038xTw7969aYiOQ2toaIitdwbORAEve9o1xmRsl12QRx5h27FjWThoEDtdcUVcq7x3KEWmbEkOk69szJkBnC7NoqIiwuEwM2fOZP369SxfvpxIcTGjgVeA+4CLvv6aAMRa0GZGt4CqqKigvr6eYDBIJBIhHA7HjUEzxphWjRpF9zlzGL5qFXLWWRAtpmpra9s8sciW5DD5zFrOTIx3HJl3hwCA/YHrgWnAj4FTcQq0qqoqamtrY7sIBINBgC162jXGdGHHHAM1Nc4xciRccEG7WszcSUpAXPeo7W5i8oJ3Y+vOPIBrcDYEXgHcDxR5vrsQeBt4E5iQyf1GjRqlJnsikYiGQiEF4o6poE2gy0APHTky7rtQKKSRSCTuHsZ0FGCpZi8fBYCfAIcDBwIDs3XvTA7LXwmam1VPOkkVVK+5pl23iEQiGg6H43JUOBy2vGRyRroc5me35hPAnqq6N/BvnIIMERkJnITTQPNz4AYRKfAtyi5s0aJFLT6b078/RwK7i/CX115jtOe7wsLC2BpnaovUmjwgIruKyC04D4NXAicDvwGeEJFFIvIrEbHhH50tEHD2+z3hBPjd7+Dqq9t8C1uSw+Qz35KOqj6uqk3Rt4sAd1PHScAcVf1WVd/FSZplye5hOoaqUlFRwZIlSygrKyMUCsV2CTjxxBP5tKSEsap8FwjwLM7/zQBuvPFGmpqaqKioYMyYMYwdOzY2vsPGf5gcdRlwO7Crqk5Q1dNU9bjoQ+NEoBA43dcIu6pttoE77nD24Dz/fPDMBM+Em2O8bAysyRe58kQ4GXg0+npH4EPPdx9FP2tBRKaIyFIRWbpmzZoODrHrmDFjBosXLyYUCrFo0SJqa2sZM2YMZWVl9O/fn6VLl7JNSQmlkQiLgTuBWmD92rUMHjyY+vp6lixZwuLFi1m3bh2RSMQWqTU5SVVPVtVnNclfQlX9TFVnqept6e4hIgUi8rKIPBx9f42IvCEiK0TkfhEp6qDwt37dusHf/w5HH+3sA/znP2d0mTfHhEKhuElKqbafMyaXdGhxJiJPisjKJMckzznTgCbgDvejJLdK+i9JVW9R1VJVLR0wYED2fwNdkEYH0roD/AGqqqqor69nzJgxVFdXU1BQwMSJE/kcOBinMKsA/g/otnZt7LpQKMSsWbMIBALU1tbGkmMgEIgNzHW7HazL0/hJRHqKSJWI/FNE7hORShHpmeHlYeB1z/ukQzZMO3Xr5uy9OWECnHkm3Hlnq5e4s8/dCUrgdGmGQiEWL17MjBkzkl5necjkjFSD0bwH0BM4DqgD7gH+BpwH/DiT69Pc95fAQqC357MLgQs97x8DxrZ2LxtQmz2tDaRtbm5u8f0JoF+Bvg+6Z/Qz7/mRSEQjkUjcNe5nwWBQg8Fg7Hz351dXV/v1R2DyANmdEPAP4E/AAdHjFuCeDK4bAjyFM4ng4STfHw3c0dp9LH9lYONG1fHjVQMB1bvuavV076QmN38lvveqrq6O+9zykOlo6XJYJkmrBliGs2j8KTiNJUcAVcBDOE+Je7d2nyT3/TnwGjAg4fMfA8uBHsDOwH+AgtbuZ8ktu5IVUqqq06dP15KSkliCa2pq0v79+yugPwH9CHQ96GHR2ZuXXHKJlpSU6OjRozUYDLaY3VleXt5itqdb+NnMKpNOlouz5Zl8luSce4FROKvNJCvOHgJOS3HtFGApsHTYsGEd+Ce1FfnqK9X99lMtKFC9995WT890xmayvGN5yHS0LS3ODm/l+x2A0tbuk+S6t3HGljVEj5s8300D3sFZSuMXmdzPirPsSZXQvC1mJSUl2tTUFCvU3GPP7bfXZaDNoBeDDujXL+77xPMBLS8vb7FsR3l5uS3LYdLKcnH2V2CM530QuKGVa45wz0lWnEXz2P2AtPbzLX+1wZdfqo4dq9qzp+qyZa2enupBM9l5tvSG6UxbVJzly2HJLTtae4JM1qU5aNAgHT16tIZCIW1qatLgXnvp30AV9AHQ7RKKscRCbPr06S0+Kykp0enTp8fFZN0LxisbxRnO5hcrcMaMRYD3cLaVjQArW7n2CpwJS+8Bq4Cvgduj37UYspHusPzVRqtXqw4d6hyrV6c8ra0FV6aFnDHZkJXiDBgAXAs8AjztHple39GHJbfsaW3sRWIC844pU3XGmAFaDroZ9FXQ3TznT506Ne76HXbYoUVrmlugeYtBe4o1XlkqznZKd7ThPrGWM1IM2Uh3WP5qh6VLndazn/3MGY+WoK1dla0VcsnON2ZLZKs4exz47+gT5njgz8BVmV7f0Yclt+xKlYgySWDeVrD9QdeArgP9eULxNXDgwKRFWbLDCjOTKBvFWbaOhOIs5ZCNVIflr3aaM0dVRPWAA5zxaAkyHeTfWiE3ffp0myxgsi5bxdmy6K8rPJ89k+n1HX1Ycut4mXR5eguz0aNHa79+/XQn0Jej49DOTyi6ysrKdODAgVpWVhb3+ejRo617waTVkcVZ9CH0dWBqR/0M72H5awv8/e/ODM5x41Q3bGjxdaYtXqkKObcwy7QFzphMpcthbdn4fHP0109F5HDgE75f1d90Ae7aQe76ZN7tUYqKiggEAixfvpxBgwZx3HHHEQgEePHFF9nYsycHi1D/zTdcCZTirDr8JTBmzBiWLFnC4MGD437Wu+++G/e+srLStl4xnUZV9xCR/jgTA0wuO+00Zy20006Dgw6CRx+Ffv1iXyfmjFQ5pKamBtXvN0R385u7JR1AXV1dbDN1bx40JutSVW2JB87MpEJgT5z1RpcBEzO9vqMPe/LsPJl0ebqTAxJnZ1ZGx6G9ibMeWnl5eYtz3KU5EsecJW6sbro2cqhbc0sPy19Z8OCDqj16qO6xh+qHH7b4ekvHjNlkAZNt6XJYxjsEqOrDqrpeVVeq6gGqOkpV52apRjR5JNWTqPukGQ6Hqa+vp1u3bjQ0NMSdW4uzwud2wJJAgNWzZ9O9e/e4c37zm98waNAgunfvHrtnayt7G7MlRGSMiLwoIl+JyHci0iwiG/yOy7TBxInw2GPw0UdOC9rnn8e+2tK9fd3zvRLvl3i+MVui1eJMRMZleAzrjIBNbhMRZs6cmfL7/v37UzJ1Kj8FVhQUcDdQ/uGHFHjOaWxs5Nhjj2XJkiVxCXHx4sW2F6fpKNcDJwNvAb2AM4HZvkZk2m78ePRf/4L334dJk2DTpi3e29d7fjgcjtuns7Kykurq6i0q/IxJJpMxZ7/K8F73Ax9sQSxmK6CqVFVVpfx+7dq1iAgbevdmv6+/5vpu3Zjy6af8qE8fLtxpJ15fv576+nr69OlDWVmZjfEwnUZV3xaRAlVtBv4iIgv8jsm0TU1NDY2NjdTedhty0kno6adz7uDB9I2OlW1PPkk31rawsJD169fH7llbWxtXyKlnDJsxbZKqvzPfDhuz4b/EWUyXXHJJ3BiNxPXNAC0oKNBTcfblXA16iEjK85ubm2M/J/Hnmq6H7O4Q8CzQHWff4KuBSjLYvilbh+WvLddiFuW116qC3g760+JibWpqissn7iLXbbl/sve2s4Bpr3Q5LOMxZwAicoyIzBSR60TkqPaXhGZr5H3CnDlzJhs2xA/ZCQTi/7qVlJTQ3NzMHTgzOD8D5qlyERAcPZo5c+bEnT9q1CjGjx9PZWUlkUgk9pfY7UJQ6+407Xc6zjCPqcBGYChwrK8RmTbxjnmtq6sjcO65XAicClQtX07ZqFFx58+dO5dIJNKm+yd7721Jc1kLv9liqaq2xAO4AWch2l9Fj3nA/2Z6fUcf9uSZOxJX9U9c/8w9vJueA9ob9A5QBZ0Lun105ua0adNiMzjdX9210dz7hkIhDYVCtihkF4LN1jRJJM6qfHDMGFXQe0F7EL/HbzZauKzlzLRXuhzWlnXOxgN7Rm+IiNyGsy+dMXECgUDcGA2vYDDIwoULCYfDzJ4dP976a5yn3IXAdcBLwHFr13L55ZcDzmSCs88+mxtvvJHVq1ezevVqlixZQigUAqC+vp5gMEh1dbU9tZqMiEg1zv9QWzNfVZ/t6HjMllFtOaty0qJFlAP1wA+Ao4Bly5ZRVVVFUVHRFuUK9+fV1dURCoWYNWtW7L2qMmvWrNg6aZaTTFu0pVvzTcA7I3MozobBxrRQU1MTa9oXEbbffntCoRALFy5EVXnuuecAKCgoiLuuoKCA64H9cP5yvsD3M1LWrl3LZZddxueeKfLgFGX19fWAs3itzZIybfAe8H4GR6M/4ZlMeQulcDhMc3Nz7MFtNnAK8DPgKWDab37DzJkztzhXuEM5gsHv1youLCykpKSERYsWMWPGjLihF8ZkqtWWMxF5COfJshB4XUSWRN8HAZvNZFLyPil6V98WkVgC22677XjppZfYdddd+eSTT1i7di39+vVj5Tff8DPg1q+/5s9AMfBboLmVn/nZZ5/Fpsfbk6ppjare5ncMJju8Y14LCwtjY1NddwGDRozg6rffRm6+mYuBy264IW4sbKZ5w3tedXU1X3zxRewBEYit7zhmzBgqKiqor6+32ZumbVL1d7oHcCBOl2bSw+3lbO0+HX3YmI3809zcrJFIJDb2bPDgwVpcXBw3S3Ngv346Z8cdVUH/D3Rg9PO99tpLi4uLW4xjmzp1atxYDxv3sfUiu7M1/wAUed5vD1yWrfu3dlj+yp5ke/x6dyG54cgj9buCAv33ttvqhVOmtHkz82R7cIZCIQ0Ggy3ykXvYGDSTTLoclknSmg+UA8MSPu8OHATcBvxXa/dJct9rgDdwukbvdxMjMBz4BmiIHjdlcj9LbvnLTW7eZFZWVhZLqL169dLTQDeCfgp6SLduKZOgu+WTe99Mkq3JT1kuzl5O8tlL2bp/a4flr+yKRCItiiXvpKGmf/1Lvyso0JdBL5oyJa6gcwspb/HlvW/ied7t5VLlJSvMTDJbWpz1BH6DM/znE+A14D844zD+CJS0do8U9z0U6BZ9fRVwlX5fnK1s6/0sueW3xBlW7uFtHfvDKafof3r00GbQy0C7pUiE7h6cqZKt2TpkuThbAfTwvO8FvJrhtQXAy8DD0ffHA68CEaA0k3tY/sq+ZHthRiKRWMtX87x5+m1BgS73tMiHQqFYi34oFNLx48dn1ErmFn7WcmbaIl0Oa3XMmapuwllG4wYR2QboD3yjqo2tXdvKfR/3vF0EHLcl9zP5S7XlDCvXxIkTUVV69erFRXfeyeVAHTANOBg4CSgsLmbcuHGx2Z/eCQKhUCg2Y7SyspKioiIbmGuSuR14SkT+gvM/1Mk4C9JmIgy8DvSNvl8JHAPcnO0gTWaS5ZTKykpmzpwZ28oJoOTkkzn+9ttZBBwONDc3U1lZyaJFi1iyZAkDBw7kmWeeAb5f/b++vj6216+X+7n7Goi99+4gYGPOTEZSVW2deQAPAafp9y1nG3GeRJ8B9ktz3RRgKbB02LBh2S1pTadI7BZIHEe2ww47KEnWRDsOdB3o56A3HXGEbt68WXv37q3dEro8i4uLdfr06fErh9sT7FaBLLSc4RkvC/wCuBZnJZcJyc5Jcv0QnAmABxJtOfN8Nx9rOet0LXYKSHjvXYcR0J+CfgLaCHpAklb4xNyTrJUsGAzGWuyrq6vjulBteIVJJV0O6+ii60mcp8jEY5LnnGk4Y84k+r4H0C/6ehTwIdC3tZ9lyS1/eZOZW1CRkCS9A3rdYxfQZaAKOrNbN90mRZeC9x5t3bLF5K4sFWfzST2m9kBaGVML3BvNU/tbcZY7kg3a9xZIzc3NcblhCOgK0K89BZrbxZlYiKUbc+b9eTYxybQmXQ5ryyK0baaqB6f7XkR+CRwBHBQNFFX9Fvg2+nqZiLwDjMBpITNbIXeZjRkzZsS2fvKufzZo0KDY1HSvT3r2ZJ9Nm6gFKpua2A9nLaO3UvychoYGxo8fj6pNZzcxP8fpwrxLRHbGWc+sJ844sseBWlVtSHahiBwBfBbNU/u39QeLyBSc1n+GDRvWytmmLbxL98D3Wyy5C8JWVVXFnf8RTiX+NE43zmE4DRcVFRUtFstevHhxbLiEd+sm74K2qbZ6MiZjqaq2jj5wkuJrwICEzwcABdHXuwAfAz9o7X725Ll1SOxyAGfmJq20ik0CXQv6JeipaVrO3MG+XraRen4iy9s3AdsAg/EsqdHK+Vfg/H/9PWAVziYXt3u+n4+1nOWUxC7OpqamuFb5AaCvgn4Feqgnd5SXl8dNMnLfu7nCnSHu/TmWV0xr0uWwNm18nmXXA9sBT4hIg4jcFP18HLBCRJbjdBmcrapf+BWk6TwafaJ1V/iORCKEQiGWLFnS6rUPAnvjbPl0O/AnnOl2Xg0NDVRWVhIOh2NP1tXV1VRWVsb+QUQiEVvNu4tS1c2q+qlmONlJVS9U1SGqOhxnbsrTqnpaR8Zotox3odqZM2fy29/+loaGBnr37g3AjiUlHAD8G6cF7URg9OjRiEgsJ5SXl7NkyRJqamqorKxk//33p6qqyi3IUVXGjh3L2LFj4z6zvGLaokO7NdNR1d1SfH4fcF8nh2NygDdxJtuTc8GCBYwaNSqui7Nnz55s2rQJcNZ5ORCYDlwMjBLhaFXejZ5bUlLCwoULefHFF9lxxx155JFH2Lx5Mw0NDbHi7J577mHVqlW2mrdpNxE5GmfHoAHAv0SkQVUn+ByWifJ2eXp3FGhsbGTmzJmUlpayf0MDc4E7gas/+YQL6+spLi5m+fLllJSU0NDQwKZNm2LvvbMxKyoqYjM5Kyoq4vbbtLxiMpaqSS3fDusW2Hp4m//dyQLe7s7i4mLdbrvtYl0O/fv3b9nN2b27fhGdzXl4ki5O76SDxOv79++vl1xyiY9/AiYTZLlb08/D8pd/vF2TcbM6v/pKV+yyiyrodX36JB0qEQwGkw7FSDajM9V6i9bd2XWly2F+dmsak1TinpyzZs0iEAjEnnJfeukl1q9fz9SpUwFnQ3SAvffem/LycgAe/O47RuGslPwwTjNGz+g9S0pKWLp0KWVlZXHXu9auXcuGDRtiXRLGmK2Xm2+8OWbmzJkE+vRhzzfeYNHIkVRt3MhN0OJ/mMFgMG5SgGvWrFnMmjUr7rPa2lpmzJgRG0YB1t1pUrPizOQ8N3nW1NRQW1tLIBBAROI2GgaYNGkSIhJb+PFdYAxwS58+TAVexBmX1tDQQGlpacqf179/f6677rqUXQ+JRZsVccZsHWpqaigsLIyNIZNttmH08uX8aYcdOAu4A2fWiJdbYHlVVFRQUVHR4rN169ZRV1cXK9Dc7s7GxkbLIyZeqia1fDusW6Br8c66cg931lWyddIOBf0Y9FvQc0GlldmfqRarbW39JNN5sG5Nk2XJFqx180p1796qoI+C9vXkGfd7Ny94uzPdtc+8a6Gl6u40XU+6HOZ7UsrWYcmt60i3Aniy8Wfu8QPQe0EV9BHQfinOSzbmzB0rkm7lcUuwncuKM5NOe8d2JXvwKy4u1rKyMv0V6Hegb4CeuPfeKfffDAaDGgwGkz7EJdvz03RNVpyZrU6qFqxx48a1uibaWaCbQN8H3Tfhu4KCgrg10bz39m7FYk++/rPizKTSWgt3Jmsdev+Nu61dJSUleusvf6kbevXSDaA3TJqUcpB/4tpn7hqLqfKH5ZCux4ozs1VKTGbJZk2lOn4K+jZoM+i1oEU9esS+69evX9rWsWRPvjYDq/NZcWaSaa2F291rN13hltj1WFZWFtd92fz++/pR//7aFAio3nlni3/vyX5GKBTSYDAYm+Xp7eJ0XyduL2d5ZOtmxZnZ6iWO60jcrDjZ0Qf0+mg352ugV518cixBpnu6TfzOu+mxNxYbh9axrDgzqaRqoUpcLiPdpujJxofF/Ttft0513DhV0Ad+9jONRFvKvOPUEn9GMBjU8vLyFkWZW7SlarHvrD+zdO9N9llxZroEb1dGdXV12gLNOzbtIJzJApHu3bV51iyNRLsfkrWOJSb2dIN/rbuzY1lxZtJJNbartaEJ3jzS2viw5o0b9eXddlMFXTZihDZv2BC7t3dbKG/hFYlE9JJLLkm6vmJnjGdNVoTZRCd/WHFmuozEboR0LWe9e/f+visTdN4226iCLigs1IsnT27ROpYqiXmffG0cWuex4syk0loB1lrRlaqVPFkRd8nFF2vd4MHaBLoCdNdoIXbJJZe0yDlu69zAgQOT5qRevXp1aB5pLX91dGFo4llxZrqUdFPXy8vL9Yc//GFcAtx7771jT61n4Wx6/DnorYcdFtf65raMpUrk6ZK9yT4rzkwyrY05SzY2NbEIacs93CU1DgFdC/oF6C9StJyle1hMdmS7xSzV78mW+PCHFWemy0ns4nSTjzvjcvTo0RoMBlusida7d2/dDXQRqILWge5bWhq7NhmbwekPK85MKqm66bzjSdO1EiW2kntb0rwTBxJbzIeDvoQz0egS0MEDB2pZWVmrrfiJXZsdlUdamy1qD5idy4oz0yUlPgkne9/c3NwiQZaUlOg2oNdFC7QPBgzQ5hUrUv6MTJK9yT4rzkw6qQa4tza+yv3eHR/mLey8y2M0NzcnHdfaC/S2aO54ELQQdOrUqS1azlKtydgZY84SizB7wPSHFWfGJJEsISUeE0E/79ZNI927q159tWpTU+xal9syZ4NpO5cVZ6a9UhVu3pzgDuB33w8aNEjLysriWuP33nvvlLnjHNCmQEA/7tdPB3k+Ly8vT7qLCTjL+DR5cky6FvvWfi+pzk3Med4uTXvA7Fw5WZwBlwIrgAbgceCHnu8uBN4G3gQmZHI/S26mLRKTcLoCbQBow667Ov9cxo7V+nPOiT1Zu/dKTKKW0DperhRnQAHwMvBw9P0PgCeAt6K/bt/aPSx/5Q7vUhje1nRvceXtpvzBD36QMnc0P/aYRvr00beiXZ6ANjU1xXWHJt6vvLw81qrvzSvJckpbZlmma+W35YD8kS6HdcM/16jqJQAiEgKmA2eLyEjgJODHwA+BJ0VkhKo2+xeq2dqICIWFhZSUlNDQ0EA4HObaa6+lV69eNDU1xc4rKSlhv/32I7x8OZN3243TlyzhrCVLuKC5mdL58zly0iTWr19PfX094XAYVUVEUm6abrZKYeB1oG/0/QXAU6p6pYhcEH1/vl/BmbYJBAIsW7aMgoKC2GcNDQ2EQiEA6uvr487/4osvUt6r8l//YqfDDuO/7rmHBcAJQGlpKUceeSSqSjAYREQoKiqiuLiY5cuXs2TJEmpqapg7d24sN0UiEaqqqigqKqKmpgZwGlbWrVsXi6e2tpaKiooWucjl/pxwOExtbS0iQm1tLQBFRUVUV1fHzne/szzmo1RVW2ceOC1lN3peX+j57jFgbGv3sCdP0x7uAOF0uwt4m/0vPvNMjRx5pCrofM/TcOLTayJb4DH7yIGWM2AI8BRwIN+3nL0JDI6+Hgy82dp9LH/ljlTDHdwxaMlyBNEWL2/O8L7/w6mnamT33XUzaBVoSXGxNjU1xXUnet+7R3FxcYuFcb1j5xJb3YC4PT1T/f7SvTedJ10O8zuxXQ58CKwEBkQ/ux44zXPOn4DjWruXJTfTXm5ymj59etLujOnTp2tTU1Ncwv4v0PXR4wzQZs84kcTuAFvgsWPkSHF2LzAK2N9TnDUmnLOutftY/soN6YY7lJSUpF3Y2rv6vzsrPG4D9PXrNXL00aqgc0H7p3iwS3X/xBnnbtdoYkzuuDiT+9LlsECaRrUtJiJPisjKJMckAFWdpqpDgTuAqe5lSW6lKe4/RUSWisjSNWvWdMxvwmz1RARVZf369TQ0NBAMBgmFQoRCIRoaGpg3bx6lpaUUFhbGrvkrsDfOgMnbgCf79yeyZg2VlZXU1dXR2NgY+0fW2NhIXV0dlZWVqGqLc0x+EpEjgM9UdVk7r7f8lWOSDXdobm6OvZ89ezbl5eWxLk4g9n7x4sWx7sCamhpEhIULF7Jw4UKne7BvX+S++4jU1XEIzoDrAyHWfejmhnTq6+tjeSMYDAIwe/bsuHO+++47IpFIlv9kTKdLVbV15gHsBKyMvrZuTeML73poeLomBg0apNBy6nv//v1186ZNWj94sH4L+gno4UlmONk09Y6Bzy1nwBXAR8B7wCrga+B2rFsz7yVuXN7c3BzXCpZs7cRMWsPd8/YCfRW0CfShffbRpu++ixus39TUlHRGp7dbM93EBbfF3+S2dDnMz8S2u+d1OXBv9PWPgeVAD2Bn4D9AQWv3s+RmsiHVmj/p9sFramrSYtDloAoaOeUU1Y8/bnFf7/VWmG05v4sz70F8t+Y1wAXR1xcAV7d2veWv3JNqJ5BM3yfeo8VsyQ0bdOmIEaqg8/v21avCYQ2FQtrU1BQrshILtFAoFDdLvKysrMW4t8RN103uSpfDOrRbsxVXRrs4VwCH4sx4QlVfBf4BvAbMA85Rm6lpOok709KdxeT69a9/Hfd+1apVhMNhCgsL+e1vf8tyYDQwA2i6+270//0/qKuDSATVlt0Vbhen2SpdCRwiIm8Bh0TfmzyTOFMxcRa2972qMmPGjNi/a/eorKyMza5MnC0547rruO3QQ7lqxx0p27CBX8+ezdiVK9l3n31oaGigX79+JHZ3L168mMrKSqqrqwmHwyxZsiTu+9mzZ8dmltpsyzyXqmrLt8OePE22JGs5S+w+cFvMkq0btCvoqzvt5LSiHXSQTo9uom4LPGYXOdRytqWH5a/8lWw4RLLNxL2D/t21EYm2dh35ox/p09GW97mg40eOTDoZwH0/evTo2Obp5eXl2tTUFHe+O7PUy3JN7kmXw/xc58yYnKP6/YD9cDjMzJkzGTVqFA0NDZSUlLBs2bLY9wCFhYVJ1w36R2EhNUOHIuEw57/wArtOmMDpM2ciIsycORNw1hayJ1tj8peqM+Gnvr4+NonIuw6a24I1Y8YMGhsb4wb/AwSDwdiA/odxZsVdA4x67TUu3WMPuh9yCBC/ttrAgQMZO3YsL774IsXFxQAMGTIkLq599tmHYDDIrFmz4iYbeNdJMzkuVdWWb4c9eZpsSdw0vaysLDbA1n3iDQaDSVfuThx3Muucc/TdQYNUQfXYYzWyenVsnz4ve6ptO6zlzOSAdNvAuSv9e1vLE9ctS7xmb9DXQSPdumnkxhuTLq+R2JJGtHV/8+bNcZ+7EwistT43pcthvielbB2W3Ew2JU4McLdUSbYYpCtxPTN3NlUAdO6++2qke3fd0KuXTuL7ffvcn2XrnrWdFWcmV6Ran8zdFinZItfp1k3rC/pk9+6qoM/vuaf28XyXal9Od5jF9OnT47pV3cMKs9xjxZkx7ZTpMhjp9q1zx6v9GHRZdFzJ30EvmjLFnmq3gBVnJhd4x495W6wSdwxobm5OWlS5Y8e8RRugAdD63r1VQT/q3l0PTlOYea9zW+dshnjus+LMmC2QaZJLVch5k3I30Obqat0cCOhnoKfaU227WXFm/Ja4ebg7BMItyKZOnRpbJzHxSPW5t7jbcccd9Wegq4qK1H2oG5jiGrdAcwszW1sx96XLYX4upWFMzlPNfBmMZEtwzJw5k6qqqtj7JqCqsZHAyy/zNs6KpU8AtWefbZMDjMkz3uUxFixYwJgxY2KTh4qKiggEAqxatSrumqamJsLhcIvPvWpra6mrq+ODDz7gq5ISdmps5PfA8TgrHJ8BFBcXc84558Rd17dvX0QkblJTJBIhHA7H7VLSVonXtOcepo1SVW35dtiTp8m2dF2VrXVtQsuuhsSuzgDo2aDrQDcXFGhk+nTVTZt8+t3mH6zlzOSIxH1zk+UA7/vNmze3+DzxnGTj1HYH/b/o0Ij/BR3cr1/S67K5n6/tDdxx0uUw35NStg5LbqYjZJqYWhtz5g7+92654p437cwz9fZowo2MHKm6cGFn/zbzkhVnJhclmxzgrovo/tt3dxhJHEPmXeE/VXE3qrhYbxs4UBX0RdB9ablrSWu7FbTl99KWB1TTNulymK1zZkwaNTU1qGqsy9Htuky2erh39W9vF2dhYSGBgDOCIBAIMHHiRMaPHx8779JbbqGyTx/47DNOff552GcfmDoVLr8cttuuc3/Dxph2U205DKKkpISZM2cSCARYtmxZbN1EgOXLlxMMBmObmG+//fYsW7aMgoKC2PUNDQ0Eg0HKysp47rnnWNbQwC+BB4F64Hngn0DV2rUMKSmhsLAQiN/hwJvD2sKbx+rq6mLrO3rznOkgqaq2fDvsydPkgkyfVlOet2GDanm5qojq0KGqjz3WIXFuDbCWM5NDMm1lSpy16bZypdrX1215dw/vd71BLwb9CvTrQECbr7pKZ1x8cda7IW3mZ8dIl8NsQoAxWZSsRa1N5223HdTXwwsvQJ8+MGECTJkC69d3SLzGmOxI1XoeDodju4GoatwEISDW0qaqSQfyNzQ0sM8++1BRUdGiVe5r4DJgJPBEJELg/PM588YbWeEZ/O/es7GxEaceaBv3Hokxt+depg1SVW35dtiTp9nqfPON6nnnqQYCqjvsoHrrrapNTX5HlTOwljOTg1K1iicuu5FsP87x48fHjRlz11ArKyuLW2Zj9OjRLScRFBfrnSecoJHhw1VB7wIt3MJlNGzMWcdKl8Os5cyYXNWzJ1x1FSxZArvvDmeeCcEgLFvmd2TGmBRStYq7LWuhUIhgMBjbLzMUCrF48WLq6+tj49NmzJhBRUUFALNmzWLhwoUMHDgQcPbZfPHFFwFnPNv06dOdFrbly1k8eDC89ho6YwbHAi8BpdBifJgmtHolvvfG3lproOkYNiHAmFw3ahQ89xzMmQNVVVBWBuXlcOmlNmHARyLSE3gW6IGTS+9V1WoRKQZuArYF3gNOVdUNvgVqcoY7wQicwscdYA/OIPuZM2ciIsybN4/FixejqrH1yVavXh13L3dT9UAgELtnUVER9OxJ5RdfsAi4G3gB+L/SUg548klk++2pqalpsQl7uk3RM50UZbLLWs6MyQcicPLJ8PrrcNZZzri0kSNh7ly/I+vKvgUOVNVioAT4uYiMAW4FLlDVvYD7gd/5F6LJNSKSdMHqwsJCqqqqUNXY7M3Zs2cTCARirWzJ7uX+WltbS3V1dWyM2ZhwmGFr1/LyyJHs/9JLbNxxRyLXXsum1avjFqTNZExapmNpTRal6u/s6AO4FFgBNACPAz+Mfj4c+Cb6eQNwUyb3szEbpktZsEB1zz1VQfWoo1Tff9/viDodOTTmDOiN04sUBDYAEv18KPBaa9db/upaWluw2rseovcIhUJpx31FIpG4tRnd45qTT9Z3dt5ZFTSy3Xb65KhR2t+2dvJduhzmZzLr63kdcouwaHG2sq33s+RmupzvvlO98krVXr1U+/RRvfpq57MuIheKM6Ag+hD5FXBV9LMFwKTo6yrgy9buY/mr62htwerEgsw9iouLY4tZu9dMnz49dt/q6upY8eZdlqO6uvr7wuvFF1VPPlkjgYB+CXp5dNKAFWb+SJfDfOvW1PgxGH1w/gIaYzK1zTZw/vnw2mtw4IFw3nlQUgLPPON3ZF2GqjaragkwBCgTkT2BycA5IrIM2A74Ltm1IjJFRJaKyNI1a9Z0WszGX+kG2U+cODHldcuXL491RYoIhYWFrF+/PvY/83Xr1lFfX8/YsWMB4rorY0pL0Tvu4IpTT+Uh4ALgbeC+Aw9Ev0v619T4JVXV1hkHcDnwIbASGKDft5xtBF4GngH2y+Re9uRpury5c1Wj0+j11FNVP/3U74g6FDnQcuY9gGrg3ITPRgBLWrvW8lfXk9halWwvTXcBWnfJjWAw2GKxWrf1zT2HNN2VLa5bulT/veOOqqCfFRZq5M9/1kgXan33W7oc1tHJ6slo4ZV4TEo470JgRvR1D6Bf9PWoaPHWN8X9pwBLgaXDhg3ruD9BY/LFxo2qF1+s2r27at++qtdeq/rtt35H1SH8Ls6AAUBR9HUv4DngCGCH6GcB4G/A5NbuZcVZ15ZYNE2fPr3FGLRQKBS3yn+ycWuJBVqy7krvmLTq6moNlZfrTUceqR8PHqwKum7bbfWFsWNVly3rxD+Brsm34izTA9iJFOPMgPlAaWv3sORmjMe//636i184/8R33VX1wQf9jijrcqA42zvawr8i+tA5Pfp5GPh39LiS6OSAdIflL+MtmlS/b0lzC7JkhVbitkqJxVmyljP3V29rWygU0khzs95y5JH6IOjmQMDJHYcfrvrKKx3/m++icrI4A3b3vC7HWSPIfRotiL7eBfgY+EFr97PkZkwS8+apjhypsVmdH3zgd0RZ43dxls3D8pdRzXxvXve7xJYzt9AKhUIaDAbjCrTE4s8tztzz4gq6tWudyUaFhc4OJWeeqbp6dbviNKmly2F+rnN2pYisFJEVwKE4T5sA44AVIrIcuBc4W1W/8CtIY/LahAnQ0ABXXw2PPQZ77AF/+ANs2uR3ZMaYBJmuJ6b6/fpkJSUlALEdBJ555hnq6+sJBoOEQiFnYVqgsbGxxfpm7nletbW1SL9+zmSjd96BUAj++ldnl5Krr+ayadPi9tZ075VsAVuzBVJVbfl22JOnMa14913VY45xWtGGD1e9+27VPH7ixVrOTBfmtoQlm0jgXedMVVtMIiChlS3ZtXHeeMPp4gRd16ePng16wVln2V6bWyhdDvM9KWXrsORmTIaeflp1772df/5jxqguXOh3RO1ixZnp6rzFl7fA8hZJ06dPjyvWEoszb2GVrNCKK7jmz9fIPvuogn4H+iToVNDqX/3KCrN2SJfDbG9NY7qaAw6Al16Cv/0Npk2DsWNh8mS48koYMMDv6IwxGfLujelVUVHBrFmzUFXmzp1LQ0OD0xqTYPHixbE9Or1bSrmbmrfYh3PcOCpLSxn5ox+x9s9/5mhgNsBf/gJvvAHnnAPHHw/du3f8b35rl6pqy7fDnjyNaYcvv1Q97zzVbt1Ui4pUZ83Km10GsJYz08V5W7sSuyi974uLi1N2Z7q7CnjvmXjvdDsZjACdu+++GhkxQhVUd9hB9eyzVR9/XHXzZr/+aPJCuhzme1LK1mHJzZgt8Oqrqoce6qSEESOcpTdyvJvCijNjtMVemsnGkDU3N7fo9vRu75RKa3uAer+vCIU08sgjqscdp9q7t5NLdtlF9aabVL/5pvP+QPJIuhzmbs6b90pLS3Xp0qV+h2FM/lKFRx+F3/7W6aIYNw6uuw5KS/2OLCkRWaaquRlcG1n+MltCVWMzO1WVQOD7hRiam5upqqqirq4u9pm7dRSknhHqvbf3ftOnT2f9+vXfd3Wq061aVFT0/YzNb76Bhx+Ga6+FJUtg221hv/2cbeYmTXJmfpq0OczPpTSMMblEBA47DFasgBtugNdfh9Gj4Ywz4KOP/I7OGJOCtzBLHH82atQo6urqCIfDRCIRwuFwbEmN1iS73/r165k5c2bsZ7pj1eKW0ujVyxl7tmgRPPWUk0Peew9+9zsYMQL23hsuuMBZ3ufrr7fo977VStWklm+HdQsYk2Xr16tecIFqjx6qPXs6r9et8zuqGKxb05iYxDFi3iU2SkpKtKmpKe68dN2Zye6XlWUz3n/fGdc6bpzqNts4XZ+9e6uedJLqffepfvxxzg+nyKZ0Ocy6NY0x6b33Hlx8MdxxB/zgB3DJJfCb3/g+I8u6NY2J586uLCwsjLVwVVVV0bdvXzZs2BDrelRPN2gm90vbhdleGzfC88/DAw/APffA5587n++wAxx0EBx1FIwZA337Ot2i3ba+xSXS5TArzowxmXn5ZTjvPHjySdh1V7jiCjjuOKc71AdWnBnTUiQSiY0xC4fDsQLNfe8WWplKLOQyLezaZPNmZ2zaSy/Biy86Y1/Xrv3++x494JBD4OijnaEWu+wCffpkNwYfWHFmjMkOVWecyO9+BytXQlmZszXU+PGdHooVZ8Yk57ZwJZsEkPXCqiM0N8PChc7EpC+/hHffhQcfhA8++P6cXXZxJi3tu6/zetgw2Gkn2GYb/+J94QVn26sMWXFmjMmu5mb4+9+dLs6PPoLDD3da0vbaq9NCsOLMmNQ0YZZlJBLJaGZmh7eStZeq80D42mvw9tuwbBk88wx84dl6u6DAKdRGjIDddnOOESOclv5vv4U1a5wJCKpOi3/Pns7khT32gMJC5x7NzbBqlXN+UxMUFUG/fhAIOF2xn3zi/OyGBue8zz+HxYudlr5ttoH334fBgzP6LaXLYVtfJ64xpuMVFMB//ReceCLMnu0UZsXFzqysmhoYPtznAI3putyWM6/Kysq0LWcdOr4sG0Schz/vA2Ak4rSqffCBUxS99Ra8+aZTvM2f7xRTmd57jz2c7tPXX4dNm1p+L+L8PFf37jBokFO4TZgAEyc6v7pF3hay4swY0369ejnj0M480ynQZs+Gu+5yJgxccokzgWArJSI9gWeBHji59F5VrRaREuAmoCfQBPxGVZf4FqjpUrxdmm5XpreLM1mBpqo0NjbGneO9R061oHkFAk6r2K67tvxO1WnZeusteOcdJ1cNGOCMVRNxvt+0yek2Xb7cWfajqcmZjLDbbs75BQXQ2AirVzuFWWGhc4+f/ARGjuzQLlTr1jTGZM9HH8GMGfDnPzuJzJ3Z2aNH1n+U392a4vzfqo+qfiUi2wDPA2Hg90Ctqj4qIocB56nq/unuZfnLZFN7WsHyfpxaHrIxZ8aYzvXKK3DuufD4485A3ZoaOP30rE6H97s48xKR3jjF2a9xirM/q+rdInIycKSqnpLuestfJtvaM36sPePUTPvl9A4BInKuiKiI9Pd8dqGIvC0ib4rIBD/jM8a0w157ObM6n3jCWbdo8mSnG+COO5wBt1sJESkQkQbgM+AJVV0MVADXiMiHwLXAhSmunSIiS0Vk6Zo1azorZNNFJBZVmRRmycapbS0NOPnG1+JMRIYChwAfeD4bCZwE/Bj4OXCDiBT4E6ExZoscfLCzftF99zkzo047DfbcE+6+O35wbZ5S1WZVLQGGAGUisidO61mlqg4FKoE/pbj2FlUtVdXSAQMGdFrMputJLLCSvfeOMUvc5skKtM7nd8tZLXAezm73rknAHFX9VlXfBd4GyvwIzhiTBSJwzDHO1PN//MMZxHvSSc7sznvv3VqKtEZgPs4D5S+Bf0a/ugfLX8ZHNTU1cQWWW4h5x56JCEVFRXFjzGprawmHwxQVFVnXpg98K85EZCLwsaouT/hqR+BDz/uPop8lu4d1CxiTLwIBZzPkFSvgzjudVcGPP94p0u66y5kplUdEZICIFEVf9wIOBt4APgHcVXkPBN7yJUDT5XlnYboFmttC1tjYGNciVlNTEzf4P+mG5qbTdOhSGiLyJDAoyVfTgIuAQ5NdluSzpG2qqnoLcAs4A2rbGaYxpjMVFMDJJ8MJJzgtaZdeCqecAtOmwUUXOcty5IfBwG3RYRcB4B+q+rCINAJ1ItIN2ARM8TFG04W5BRZAXV1dbCZmqlmYbR2nZjqOL7M1RWQv4Cng6+hHQ3CeNsuAXwGo6hXRcx8DalR1Ybp72mwnY/JUJAIPPeSsk1ZSAjfdlNFluTRbc0tZ/jIdyWZh5qacm62pqq+o6g6qOlxVh+N0Xf5UVVcBc4GTRKSHiOwM7A7YAo7GbK0CAZg0ydmbbtYsv6MxZqtiszDzk98TAlpQ1VeBfwCvAfOAc1R165l7b4xJzt3rzhiTFTYLM3/lxPZN0dYz7/vLgcv9icYYY4zJf6lmYQI2CzPH5URxZowxxpjsq6mpidsdwC3QrDDLbTnXrWmMMcaY7LFZmPnHijNjjDHGmBxixZkxxhhjTA6x4swYY4wxJof4sghtRxCRNcD7fscB9AfW+h1EO+Vr7BZ358qVuHdS1a1ix/Acyl+QO/9928ri7lwW95ZLmcO2muIsV4jI0nxdtTxfY7e4O1e+xm0yk6//fS3uzmVxdyzr1jTGGGOMySFWnBljjDHG5BArzrLvFr8D2AL5GrvF3bnyNW6TmXz972txdy6LuwPZmDNjjDHGmBxiLWfGGGOMMTnEirMtICIFIvKyiDyc4vv9RaRBRF4VkWc6O75U0sUtIoUi8pCILI/G/Ss/YkwkIu+JyCvRP8+lSb4XEakXkbdFZIWI/NSPOJPJIPZTozGvEJEFIlLsR5yJWovbc95oEWkWkeM6Mz6zZSx/da58zWGWv/xhG59vmTDwOtA38QsRKQJuAH6uqh+IyA6dHFs6KeMGzgFeU9UjRWQA8KaI3KGq33VqhMkdoKqp1qf5BbB79AgCN0Z/zRXpYn8XGK+q60TkFzhjInIl9nRxIyIFwFXAY50XkskSy1+dL19zmOWvTmYtZ+0kIkOAw4FbU5xyCvBPVf0AQFU/66zY0skgbgW2E2dn3G2BL4CmTgpvS0wC/qaORUCRiAz2O6hMqOoCVV0XfbsIGOJnPG1UDtwH5MTfb5MZy185KS9zmOWvjmHFWfvNAs4DIim+HwFsLyLzRWSZiJzRaZGlN4v0cV8P7AF8ArwChFU11bmdSYHHo3+WU5J8vyPwoef9R9HPckFrsXv9N/BoJ8SUibRxi8iOwNHATZ0emdlSs7D81dnyNYdZ/vKBdWu2g4gcAXymqstEZP8Up3UDRgEHAb2AhSKySFX/3TlRtpRh3BOABuBAYFfgCRF5TlU3dEqQqe2rqp9Eu1eeEJE3VPVZz/eS5JpcmYrcWuwAiMgBOMntZ50eYXKtxT0LOF9Vm52GCpMPLH/5Jl9zmOUvH1jLWfvsC0wUkfeAOcCBInJ7wjkfAfNUdWO0z/tZwO+BkpnE/Suc7gxV1bdxxhP8qHPDbElVP4n++hlwP1CWcMpHwFDP+yE4T8++yyB2RGRvnK6aSar6eedGmFwGcZcCc6J/n44DbhCRozozRtMulr98kK85zPKXT1TVji04gP2Bh5N8vgfwFM4TaG9gJbCn3/FmEPeNQE309UDgY6C/z7H2AbbzvF6AM1DZe87hOM3pAowBlvj9Z9yG2IcBbwP7+B1vW+JOOP+vwHF+x21Hm/87W/7qnHjzModZ/vLvsG7NLBKRswFU9SZVfV1E5gErcMZH3KqqK30NMAVv3MClwF9F5BWcJHG+ppnt0kkGAvdHm567AXeq6ryEuB8BDsNJEl/jPEHngkxinw70w3lyA2hS/zfmzSRusxWx/NWh8jWHWf7yie0QYIwxxhiTQ2zMmTHGGGNMDrHizBhjjDEmh1hxZowxxhiTQ6w4M8YYY4zJIVacGWOMMcbkECvOjDHGGGNyiBVnJq+ISC8ReUZEClJ8311EnhURW8PPGJNTLH+ZTFlxZvLNZJztWZqTfamq3+GsbH5ip0ZljDGts/xlMmLFmck3pwIPisi2IvKUiLwkIq+IyCTPOQ9EzzPGmFxi+ctkxHYIMHlDRLoDH6jqoGizf29V3SAi/YFFwO6qqtEug1WqOsDXgI0xJsryl2kL69c2+aQ/0Bh9LcAfRGQczt5/O+Lsp7ZKVZtF5DsR2U5Vv/QnVGOMiWP5y2TMijOTT74BekZfnwoMAEap6mYRec/zHUAPYFPnhmeMMSlZ/jIZszFnJm+o6jqgQER6AoXAZ9HEdgCwk3ueiPQD1qjqZp9CNcaYOJa/TFtYy5nJN48DPwPuAB4SkaVAA/CG55wDgEc6PzRjjEnL8pfJiE0IMHlFRH4CVKnq6WnO+Sdwoaq+2XmRGWNMepa/TKasW9PkFVV9Gfi/dIs4Ag9YYjPG5BrLXyZT1nJmjDHGGJNDrOXMGGOMMSaHWHFmjDHGGJNDrDgzxhhjjMkhVpwZY4wxxuQQK86MMcYYY3LI/wckSRVYSvQYPAAAAABJRU5ErkJggg==\n",
"text/plain": [
"
"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"pl.figure(figsize=(10,3.5))\n",
"fig = pl.subplot(1,2,1)\n",
"fig.plot(prediction_mean[\"(a)\"], prediction_mean[\"(b|a)\"], color=\"red\")\n",
"fig.fill_between(prediction_mean[\"(a)\"],\n",
" (prediction_mean - prediction_std)[\"(b|a)\"],\n",
" (prediction_mean + prediction_std)[\"(b|a)\"],\n",
" color=\"red\", alpha=0.5)\n",
"fig.scatter(data[\"(a)\"], data[\"(b|a)\"], marker=\"x\", color=\"k\")\n",
"fig.set_xlabel(\"(a)\")\n",
"fig.set_ylabel(\"(a|b)\")\n",
"fig = pl.subplot(1,2,2)\n",
"fig.plot(prediction_mean[\"(a)\"], prediction_mean[\"(c|a,b)\"], color=\"red\")\n",
"fig.fill_between(prediction_mean[\"(a)\"],\n",
" (prediction_mean - prediction_std)[\"(c|a,b)\"],\n",
" (prediction_mean + prediction_std)[\"(c|a,b)\"],\n",
" color=\"red\", alpha=0.5)\n",
"fig.scatter(data[\"(a)\"], data[\"(c|a,b)\"], marker=\"x\", color=\"k\")\n",
"fig.set_xlabel(\"(a)\")\n",
"fig.set_ylabel(\"(c|a,b)\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The uncertainty margin consists of two contributions. The learned variance of the quadratic regressions and the uncertainty of the regression parameters.\n",
"\n",
"With the regression equation,\n",
"\n",
"$y(x) = a \\cdot x + b \\cdot x^2 + c + \\xi$,\n",
"\n",
"the learned variance describes the strength of the random noise $\\xi$ and the uncertainty of the regression parameters quantifies that with finite data we know $a$, $b$, and $c$ only with a finite precision."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Backwards prediction"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can also evaluate the inverse prediction, e.g. predicting '(b|a)' and '(a)' from '(c|a,b)'.\n",
"We just have to provide data for '(c|a,b)' to the predict method and the causal structure will solve the inverse model."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"