{ "cells": [ { "cell_type": "markdown", "id": "rocky-imaging", "metadata": {}, "source": [ "# Results for the wine quality dataset\n", "\n", "\n", "\n", "This notebook produces the results of EiV and non-EiV models for \n", "the wine quality dataset (taken from https://archive.ics.uci.edu/ml/datasets/wine+quality) as presented in the preprint \"Errors-in-Variables for deep learning: rethinking aleatoric uncertainty\" submitted to NeurIPS 2021.\n", "\n", "\n", "This notebook produces Figures 5a and 6 of the preprint. \n", "\n", "How to use this notebook: \n", "\n", "+ This notebook assumes that the corresponding trained networks exist in `saved_networks`. To achieve this, either run the training scripts described in the `README` or load the pre-trained networks from the link in the `README` into the `saved_networks` folder. \n", "\n", "+ To run this notebook, click \"Run\" in the menu above. \n", "\n", "+ To run this notebook with a GPU, set `use_gpu` to `True` in cell [2] (default is `False`)\n", "\n", "+ Plots will be displayed inline and, in addition, saved to `saved_images`\n" ] }, { "cell_type": "code", "execution_count": 31, "id": "attractive-punch", "metadata": {}, "outputs": [], "source": [ "import random\n", "import os\n", "\n", "import numpy as np\n", "import torch\n", "import torch.nn as nn\n", "from torch.utils.data import DataLoader, TensorDataset\n", "import matplotlib.pyplot as plt\n", "import matplotlib\n", "from tqdm.notebook import tqdm\n", "\n", "from EIVArchitectures import Networks\n", "from generate_wine_data import test_x, test_y, train_x, train_y\n", "from EIVTrainingRoutines import train_and_store\n", "\n", "%matplotlib inline" ] }, { "cell_type": "markdown", "id": "92190849", "metadata": {}, "source": [ "## Fix relevant hyperparameters" ] }, { "cell_type": "markdown", "id": "7f0337bf", "metadata": {}, "source": [ "### Values that can be changed" ] }, { "cell_type": "code", "execution_count": 32, "id": "2a8a4fbc", "metadata": {}, "outputs": [], "source": [ "# the Deming factor to use for the scatter plots (Figure 5)\n", "# Choose one of 0.01,0.1, 0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.6 and 2.0\n", "deming = 1.0\n", "\n", "# Switch to True if GPU should be used\n", "use_gpu = False\n", "\n", "# Uncertainty coverage factor (1.96 taken from the standard normal)\n", "k=1.96" ] }, { "cell_type": "code", "execution_count": 33, "id": "copyrighted-taxation", "metadata": {}, "outputs": [], "source": [ "# graphics\n", "fontsize=15\n", "matplotlib.rcParams.update({'font.size': fontsize})" ] }, { "cell_type": "code", "execution_count": 34, "id": "3f5515e7", "metadata": {}, "outputs": [], "source": [ "# Set device\n", "if not use_gpu:\n", " device = torch.device('cpu')\n", "else:\n", " device = torch.device('cuda:0' if torch.cuda.is_available() else 'cpu')" ] }, { "cell_type": "markdown", "id": "d6ee9d43", "metadata": {}, "source": [ "### Values to keep fixed\n", "The following values assume the settings from the training scripts. To change the following values, these scripts must be adapted and rerun." ] }, { "cell_type": "code", "execution_count": 35, "id": "564904dc", "metadata": {}, "outputs": [], "source": [ "from train_eiv_wine import dim, init_std_y_list, precision_prior_zeta, deming_factor_list, seed_list\n", "init_std_y = init_std_y_list[0]" ] }, { "cell_type": "code", "execution_count": 36, "id": "auburn-immune", "metadata": {}, "outputs": [], "source": [ "# function to fix seeds (for reproducability)\n", "def set_seeds(seed):\n", " torch.backends.cudnn.benchmark = False \n", " random.seed(seed)\n", " np.random.seed(seed)\n", " torch.manual_seed(seed)" ] }, { "cell_type": "markdown", "id": "extraordinary-suspension", "metadata": {}, "source": [ "## Comparison EiV and non-EiV for one deming factor\n", "Produces Figure 6 of the preprint" ] }, { "cell_type": "code", "execution_count": 37, "id": "daily-draft", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "FNNEIV(\n", " (main): Sequential(\n", " (0): EIVInput()\n", " (1): Linear(in_features=11, out_features=200, bias=True)\n", " (2): LeakyReLU(negative_slope=0.01)\n", " (3): EIVDropout()\n", " (4): Linear(in_features=200, out_features=100, bias=True)\n", " (5): LeakyReLU(negative_slope=0.01)\n", " (6): EIVDropout()\n", " (7): Linear(in_features=100, out_features=50, bias=True)\n", " (8): LeakyReLU(negative_slope=0.01)\n", " (9): EIVDropout()\n", " (10): Linear(in_features=50, out_features=1, bias=True)\n", " )\n", ")" ] }, "execution_count": 37, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Load EiV model\n", "net = Networks.FNNEIV(p=0.5, init_std_y=init_std_y,\n", " precision_prior_zeta=precision_prior_zeta, deming=deming,\n", " h=[dim, 200,100,50,1])\n", "saved_file = os.path.join('saved_networks', \n", " 'eiv_wine_init_std_y_%.3f_deming_factor_%.3f_seed_%i.pkl'\n", " % (init_std_y, deming, 0))\n", "train_loss, test_loss, stored_std_x, stored_std_y, state_dict, extra_list = train_and_store.open_stored_training(saved_file, \n", " net=net, extra_keys=['rmse'], device=device)\n", "rmse = extra_list[0]\n", "net.to(device)" ] }, { "cell_type": "code", "execution_count": 38, "id": "swiss-paris", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "FNNBer(\n", " (main): Sequential(\n", " (0): Linear(in_features=11, out_features=200, bias=True)\n", " (1): LeakyReLU(negative_slope=0.01)\n", " (2): Dropout(p=0.5, inplace=False)\n", " (3): Linear(in_features=200, out_features=100, bias=True)\n", " (4): LeakyReLU(negative_slope=0.01)\n", " (5): Dropout(p=0.5, inplace=False)\n", " (6): Linear(in_features=100, out_features=50, bias=True)\n", " (7): LeakyReLU(negative_slope=0.01)\n", " (8): Dropout(p=0.5, inplace=False)\n", " (9): Linear(in_features=50, out_features=1, bias=True)\n", " )\n", ")" ] }, "execution_count": 38, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Load non-EiV model\n", "ber_net = Networks.FNNBer(p=0.5, init_std_y=init_std_y,\n", " h=[dim, 200,100,50,1])\n", "ber_saved_file = os.path.join('saved_networks', \n", " 'noneiv_wine_init_std_y_%.3f_seed_%i.pkl'\n", " % (init_std_y, 0))\n", "ber_train_loss, ber_test_loss, ber_stored_std_x,\\\n", " ber_stored_std_y, ber_state_dict, ber_extra_list\\\n", " = train_and_store.open_stored_training(ber_saved_file, net=ber_net, extra_keys=['rmse'], device=device)\n", "ber_rmse = ber_extra_list[0]\n", "ber_net.to(device)" ] }, { "cell_type": "code", "execution_count": 39, "id": "91a110fd", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "std_y for EiV: 0.593\n", "std_y for non-EiV: 0.587\n" ] } ], "source": [ "# print std_y for EiV and non-EiV\n", "print('std_y for EiV: %.3f' % (stored_std_y[-1]))\n", "print('std_y for non-EiV: %.3f' % (ber_stored_std_y[-1]))" ] }, { "cell_type": "code", "execution_count": 40, "id": "brave-fiction", "metadata": {}, "outputs": [], "source": [ "# Collect predictions for test data\n", "\n", "set_seeds(0)\n", "\n", "# Save network state\n", "ber_net_train_state = ber_net.training\n", "net_train_state = net.training\n", "net_noise_state = net.noise_is_on\n", "\n", "# Switch on Dropout and input noise\n", "ber_net.train()\n", "net.train()\n", "net.noise_on()\n", "\n", "# Collect predictions and compute errors and uncertainties\n", "pred, _ = [t.cpu().detach().numpy()\n", " for t in net.predict(test_x.to(device), number_of_draws=1000,\n", " take_average_of_prediction=False)]\n", "ber_pred, _ = [t.cpu().detach().numpy()\n", " for t in ber_net.predict(test_x.to(device), number_of_draws=1000,\n", " take_average_of_prediction=False)]\n", "val_y = test_y.detach().cpu().numpy()\n", "err = np.mean(pred, axis=1).flatten()-val_y\n", "unc = np.std(pred, axis=1).flatten()\n", "ber_err = np.mean(ber_pred, axis=1).flatten()-val_y\n", "ber_unc = np.std(ber_pred, axis=1).flatten()\n", "\n", "# Restore networks\n", "if net_train_state:\n", " net.train()\n", "else:\n", " net.eval()\n", "if net_noise_state:\n", " net.noise_on()\n", "else:\n", " net.noise_off()\n", "if ber_net_train_state:\n", " ber_net.train()\n", "else:\n", " ber_net.eval()" ] }, { "cell_type": "code", "execution_count": 41, "id": "descending-newark", "metadata": {}, "outputs": [], "source": [ "def diagonal(x_array, eps=0.1):\n", " min_x = np.min(x_array)\n", " max_x = np.max(x_array)\n", " x = np.linspace(min_x-eps, max_x+eps)\n", " return x,x" ] }, { "cell_type": "markdown", "id": "e770c860", "metadata": {}, "source": [ "Figure 6 in the preprint" ] }, { "cell_type": "code", "execution_count": 42, "id": "mobile-consumption", "metadata": { "scrolled": true }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAARYAAAEOCAYAAABMyDi5AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAnrUlEQVR4nO3de5xVZb348c+XAQRSUW6pKTMSnjxW51SOFmlqICJTinIxE1HUIvCUmKXnFMcO1cFz1F+Wv2OaVN7HGdO8VWapeUkFETKTY5goAnnhIoiLu8x8zx/P2rBmsfbea+9Ze699+b5fr/Xas+7fmWG+POt5nvU8oqoYY0ySeqQdgDGm9lhiMcYkzhKLMSZxlliMMYmzxGKMSZwlFmNM4nqmHUA5DBo0SJuamtIOw5ias2jRorWqOji8vS4SS1NTEwsXLkw7DGNqjogsj9puj0LGmMRZYjHGJK6siUVEJonI/SLyuohsFJFFIvLFGOf1F5EbRWS9iGwQkVYRGViOmI0xhSt3ieUiYCPwdeBk4FHgdhH5Wp7z7gCOA74ETAWOAO4tVZDGmO4pd+XtSaq6NrD+BxE5AJdw/ifqBBEZAYwBjlXVJ/xtrwPPiMjxqvpwqYM2xhSmrCWWUFLJeA4YkuO0scCqTFLxr7MAWObvM8ZUmEqovP008GKO/YcCSyK2/9XfZ4ypMKkmFhEZBYwDfpzjsH2BdyK2r/f3GWMS4nkejz/+eLevk1piEZEm4HbgPlW9Kc/hUaNRSZbtmetPE5GFIrJwzZo1RcdpTL3wPI+xY8fS0tLC6tWru3WtVBKLiAwAfgusAM7Mc/h6YJ+I7fsQXZIBQFXnqmqzqjYPHrxbj2NjTEAmqcyfP5+bb76ZIUNyVXvmV/bEIiL9gF8DvYHPqeqmPKcsIbouJVvdizGmAMGk0t7ezsSJE7t9zXJ3kOsJ3AkcAoxV1Tjlrd8C+4nI0YHrNAPD/H3GmG644YYbEk0qAFLOwbRFZC7wZWAmsCC0+zlV3SYiS4HHVfW8wHkPAv8AfBPoBC4HVqvqZ+Lct7m5We0lRGOiqSrPPfccn/jEJwo+V0QWqWpzeHu5H4VO8D+vBuaFlv39fT2BhtB5pwOPAzcAtwCLgFNLHawxtcrzPL7whS+wdOlSRKSopJJLWXveqmpTMceo6jvAOf5ijOmGYJ3KGWecwfDhwxO/RyV0kDP1qLUVmpqgRw/32dqadkR1IVxRO27cuJLcpy4GejIVprUVpk2DzZvd+vLlbh1g8uT04qpxwaTS1taWWEVtFCuxmPKbNWtXUsnYvNltNyWjqjQ0NNDW1sakSZNKei8rsZjyW7GisO2mWzzPo0ePHuy99948+uij9OhR+vKElVhM+Q0dWth2U7TM48+pp56KqpYlqYAlFpOGOXOgX7+u2/r1c9tNYoJ1Kl/+8pcRkbLd2xKLKb/Jk2HuXGhsBBH3OXeuVdwmKNz6U+o6lTCrYzHpmDzZEkkJTZ06NfFu+oWwEosxNej73/8+d911VypJBSyxGFMzPM/juuuuQ1U57LDDOOWUU1KLxRKLMTUgU6fyta99jb/85S9ph2OJxZhqF66o/ed//ue0Q7LEYkw18zyPlpaWVCtqo1hiMaaKPfPMMyxcuLCikgpYc7MxVUlVERGOP/54Xn31Vfbff//8J5WRlViMqTKe5zFy5EjuvvtugIpLKmCJxZiqkqlT+eMf/0hHR0fa4WRlj0LGVIlMUpk3b17F1amEWYnFmCqwdevWqkkqYInFmKqwxx57MGLEiKpIKmCPQsZUNM/zWLVqFcOHD+eKK65IO5zYrMRiTIXK9KgdOXIkW7ZsSTucgliJxZgKFO6m37dv37RDKoiVWIypMKWYS7ncLLEYU2Fmz55d1UkF7FHImIrzve99j5aWFkaNGpV2KEWzEosxFcDzPC688EI8z+N973tfVScVsMRiTOoydSrXXHMN8+fPTzucRGRNLCIyUUT6lDMYY+pNeDyV0aNHpx1SInKVWH4BrBKRm0XkRBFpKFdQxtSDanr3p1C5EsuRwE+B44AHgLdE5FoR+Uw5AjOm1q1Zs4aVK1fWXFIBEFXNf5DI0cDpwARgCPAG0A60q+qikkaYgObmZl24cGHaYRgDwJYtW+jTpw8iwtatW+nTp3prHERkkao2h7fHqrxV1SdV9avAB4ATgAeBc4AFIvI3EfluotEaU6M8z2P06NF84xvfAKjqpJJLQa1Cqtqpqo+o6peBg4DrgA8C/16K4IypJcEetSNGjEg7nJIqqIOcuFmlj2XXY9FA4G+4xyJjTBbBpNLW1lb2uZTLLVZiEZERuGQyCXg/8DpwE9Cmqn8qWXTG1ABVZdy4cXWTVCBHYhGRj+OSyWnAUOBt4C5cMvljecIzpvqJCDNnzmTGjBl1kVQgd4llEeAB9wJtwEOqWrmj9xpTYTzPY968eZxwwgmMGzcu7XDKKlfl7SRgiKqeraoPJpVURGS4iFwvIs+LSIeIPBbjnCYR0YjF6nZMRcp0fjv55JN544030g6n7LKWWFT1lyW654eBFmA+0LvAc78JPBVYX5tUUMYkJdyj9oADDkg7pLLLVceyGhijqs+JyBogZ086VR0S856/UtX7/HvcBQyKGyzwkqrWxltapibVwiBNSchVx/JjYFXg6/xddGNQ1c4krmNMJWpra6v7pAIxu/SX7OZ+iUVVj8tzXBOwDPfoMwBYjatQnqWqeUcZti79plxUlcWLF/PRj3407VDKoltd+vNcuKeIlPohchuu1HQeMAq4HpiBdcwzFcDzPMaPH8/ixYsRkbpJKrnkGo9lu4gcEVjvISJ/EJFDQoceDqwsVYAAqvqmqn5VVe9X1cdUdTZwEXCyiHws6hwRmSYiC0Vk4Zo1a0oZnqljmTqV+++/n5dffjntcCpGrhJLT0AC64IbQmGvUgZUgLv8z09E7VTVuararKrNgwcPLmNYpl6EK2pPPfXUtEOqGNU8NKWGPo0pG2v9ya2aE0vmN1nx48GY2tPQ0MBee+1lSSWLsk//ISL9cB3kwI3vsreIZH4zD6jqZhFZCjyuquf558zGPYI9BbwLHANcDNytqn8pZ/ymvnmeR2dnJ/379+eBBx7AvfBvwvIllq+JyJv+15mf4EwRWRU4Zv8C7zkEuDO0LbN+MPCaH1dwjN0luF63XwL6AiuAK4E5Bd7bmKJlHn8AnnjiCXr0qOYCf2nlSiwrgKND25bjSgtRx8aiqq/RtVI46pim0Ho71rRsUhSuU7Gkkluud4WayhiHMRXLKmoLZ2nXmDymTZtmSaVAlliMyWPOnDn88pe/tKRSAEssxkTwPI+rrrqKzs5Ohg0bVncDNXWXJRZjQjLjqVxyySX86U82pHMxyt6PxZhKFq6obW7e7cVdE0NRJRYRGSoilpRMTbHWn+QUnFj8yeGXAf+UfDjGpOf555/nz3/+syWVBBRb6rB+zKZmdHZ20qNHD44++miWLVuGvQ3ffVZ5a+qa53l89rOf5ZZbbgGwpJKQYhKL4rr2b0s4FmPKKlOn8tRTT9GvX7+0w6kpBT8K+YNhH1yCWIwpG6uoLS17FDJ1Z9u2bZZUSswSi6k7e+yxB2PGjLGkUkLWF8XUDc/zWLlyJYcddhiXXnpp2uHUNCuxmLqQqVMZOXIkGzduTDucmherxCIiA1R1XamDMaYUghW1bW1t7LnnnmmHVPPilljeFJFfiMhYEbFSjqka4aQyadKktEOqC3GTxHTcWLW/BlaKyGUi8qHShWVMMv7rv/7LkkoKCpq7WUSGAVOBKcBQYD5wA3CHqlbsg6vN3Vy/tm7dyjPPPMOxxx6bdig1KZG5m1X1VVX9jqoeDIwGOoC5wFsicpOIRM5KaEw5eZ7HjBkzWLduHX369LGkkoJi3m7uJyJTge/gRvF/Efgh8I/AsyJycaIRGlOAzCBNP/3pT1mwYEHa4dSt2IlFRI4RkRuBt4CrgZeAT6nqR1X1UlX9JPAt4N9KE6oxuWWSyrx582hvb+fEE09MO6S6FSuxiMgrwKPAcOACYH9V/Yqqhv9LeATYN9kQjckvnFSsR2264va8/SXwM1X9W66DVHUR1unOpGDDhg2sXr3akkqFiJtYFgNvR+0QkQHA51X1lsSiMiamTZs20bdvXw488EBeeOEFevfunXZIhvilixuBD2bZd7C/35iy8jyPMWPGMH36dABLKhUkbmLJNRTlQODdBGIxJrZgj9rRo0enHY4JyfooJCLjgOAsTZeKyJrQYX2AzwDPliA2YyJZN/3Kl6uOZQjw0cD6B4H9QsdsB34P/GfCcRkTSVWZMGGCJZUKF6tLv4g8CsxQ1SWlDyl51qW/tjzyyCOsX7/eWn8qQLYu/bFahVT1s8mHZEx8nufxhz/8gXHjxjFq1Ki0wzF5xB5BTkQOAD4PHIirWwlSVf3XJAMzJiNTp7JgwQJefvllGhsb0w7J5BG35+2pwKvAj4HzgEkRi0lbays0NUGPHu6ztTXtiLotWFHb2tpqSaVKxC2xXIarpJ1qI8lVqNZWmDYNNm9268uXu3WAyZPTi6sbrPWnesXtx3IQ8P8tqVSwWbN2JZWMzZvd9ip17733WlKpUnFLLE8DHwIeLmEspjtWrChsexWYMmUKRxxxBIceemjaoZgCxS2xXARME5GzReQAf0yWLkspgzQxDB1a2PYK5Xken//853n2Wdfn0pJKdYqbWP6C6yx3I7AS8CKWWERkuIhcLyLPi0iHiDwW87z+InKjiKwXkQ0i0ioiA+Pet+bNmQPh+Yf79XPbq0SmTuXBBx9kRRWXtEz8R6FzcZPBJ+HDQAtuvNxC3hq7A/c49iWgE7gcuBf3SoHJVNDOmuUef4YOdUmlSipuwxW1EyZMSDsk0w0FDaadyA1FevgTyyMidwGDVPW4POeMwNXzHKuqT/jbjgSeAUaras66H+t5W9k2btzIiSeeaBW1VSiRwbSTkEkqBRoLrMokFf86C4Bl/j5TxXr37s1+++1nSaWG5Hq7eQGu38qLIvIseR6FVPXIpIMLOBSIek/pr/4+U4U8z2Pbtm0MGjSIO++8E5Fco3OYapKrxPK/wJbA1/mWUtoXeCdi+3pqeYzdGuxJm5GpUxkzZgwdHR2WVGpM1hKLqp4T+HpqWaLJLarEJFm2IyLTgGkAQ6usyRWoyZ60GcGK2vb2dhoaGtIOySSsWga+Xg/sE7F9H6JLMqjqXFVtVtXmwYMHly6yUqnBnrSwe1KxoQ9qUyFvNzcBZwL/wO5vN6OqpyUX1m6WEN2sfCiuybn21GBPWoCvfvWrO1t/LKnUrrhvNx+OG6l/sr8cAjQDE4FPAYNKFaDvt8B+InJ0IKZmYJi/r/bUSE/asMsuu4x77rnHWn9qXNxHoStxcwt9BFevcZ6qDsNNsarAFXFv6L8CMFFEJgIfAAZn1jOvBojIUhH5eeYcVZ0H/A64RUTGi8gpQCvwZL4+LFWrBnrSZniex2WXXUZHRwcf+MAHOOmkk9IOyZSaquZdgHXAGFxS6QQ+Hdh3LvDnONfxj2/CJaOopck/5jXgptB5++BeKXgHNyvA7bjOdXnvefjhh2tVuu021cZGVRH3edttaUdUsHfffVePOuoobWho0KeffjrtcEzCgIUa8TcXt45Fge2qqiKyGmjE9YQF9+7QIQUkstfIPZ0IqtoUse0d4Bx/qQ+TJ1d1C1C4m/6IESPSDsmUSdxHoRfZNWHZPODrInKIiDQClwCvlCI4U71skKb6FrfEMhdXSgH4Nm40uUxP2E24SlxjdlqyZAmLFy+2pFKn4o7Sf2vg67+KyD8CI4C+wHxVXV2i+EyV2bFjBz179uSII45g2bJl7Ltv7XaMNtnFbW4+Kzj2iapuVNWHVPV+YIeInFWyCE3V8DyPkSNHcu211wJYUqljNim8SYTnebS0tPD0008zZMiQtMMxKbNJ4U23ZZLKvHnzrJu+AWxS+PrV2prIaHM7duywpGJ2Y5PC16ME35zu2bMnEyZMYObMmZZUzE42KXy9aW2Fs8+Gjo7d9zU2wmuvxbqM53ksXbqUj3/848nGZ6pK0UNTikgf3Ds9TSWIy2RTikGeMiWVqKQCsd+cznR+GzVqFBs2bOh+XKbm5O3HoqpbRWQf3DtCphxKNchT1BgvQTHenA6Pp9K/f//i4zE1K26rUCv19I5OKcQpgWSOOfPMwgZ5Ov986NkTRNzn+edHH5erRBLjzWkbpMnEFvVmYngBvg68DiwEvgf8C3B+YJkR5zppLam/3Xzbbar9+qnCrqVfv65vK0cdE15Edr/2jBnRx86YsfuxjY3RxzY0xHpzevbs2drQ0KB33nln8T8LU1PI8nZz3MTSmWfpiHOdtJbUE0u2P+jGxvzHZDs+o6EhfhKKk+By2L59uz755JNF/ABMrcqWWGI9CqlqjzyLjYacS5xhJvNVnPbuHf2okq0iVnX3x63Jk2HuXNf6I+I+587NWW/jeR7nnnsub731Fr169eKoo47KHacxVM9g2tUt3zCTra2u7iWXXr1cHUu4jibXCPdRdTKTJ7sm5c5O95knqbS0tHDLLbewaNGi3PEZExA7sYjIEBG5XEQeEZG/iciH/e0z/SlQTTa5hpnM1wScsWmTax1S3dVKlDk3m24MvB3upv+5z32u6GuZ+hP37eYjgZeBCbhhIz8I7OHv3h/4RimCqxnZHkHAdVbL1QScTaaV6NprYc89o48pcuBte/fHdFfcEssPgUdxU398ha4vJS4ASjm9anXLNCFPmeLWb70VWlrgrLNcs3K+kkouy5e7z5/8JN7A2zE73W3evJkNGzZYUjHFi6rRDS+4qVZP8L9uwLUEfcJfPxbYGuc6aS2ptQrddptqr167t9bka/2Ju2SaiQcO7Lp94MDdW3pitAh5nqfbt29XVdX33nuvnD8pU6XoTqsQsAHINp3gMGBVN3Jb7Zo5E957r+s2zf9uVmwdHXDOOfD22123e97ux+aZWdHzPE488UTOPvtswL1caEyx4iaW+4DvisiwwDYVkUHAN4G7E4+sFoT/4JPW0LB74gLYvn33FqEcTd7BHrXjx49PPk5Td+Imln/DDeb0IvCEv+0nwEu4x6TvJB+ayUkkd/1MOJFkqcj1DjyQlpYW66ZvEhW3g9x63FSq/wIsBx4GluESzlGqGlH2NgwcmP+YYojAyJHuM5twIsnS5H36gAHW+mMSF/tBWlW3Az/3FxPH1VfD1KmwY0ey1731Vld/k62+JqqXbqYjXGjUuFkHH8y5b77JhAkTko3R1LVYiUVERgEHqepNEfumAstV9dFkQ6sRuUoVxWj0p3fKVX9zww3RPWr9mRU9z+M3v/kNp59+Op9ONjpjgPh1LHOA92fZNwi4LJlwakhmpLaoytViZfqmZBs+AVziydNNf+zYsUyZMoWlS5cmF5sxAXETy4dxQyZEeQ44LJlwakTcbvqFyrwwmKurfo4xVcLTng4fPjzZ+IzxxU0sO4ABWfaVqIayCuUaqCkJTz3lPrN11R84sGtpJdDT1hs6lLHNzdb6Y8ojqtdceAF+heu63zu0vTfwDPDrONdJaylLz9sZM5LtVZttmTGjqIGj7gTtCXrnBReU/mdh6gbdHOjpn3D9WFYAVwIX+Z/LgXeAj8S5TlpLyRPLbbeVJ6lkuvFn7tnY6O7b2Lh7F35/4KjOwLmvZBssypgiZUsscfux/AU4AngKmAJc7n8+CRypqouTKD1VrVmzku2qn0um3iYzrsr06fD3v7vHr+B4tytW4AFj2dWjcZi/3ZhSK6Qfy0vAF0sYS/Uq9x9rU5OrpH3qKbjuul3bOzp2rnsHHsjYlSuZD3wpeG6RQykYUwh70ywJAwaU/r2goMxAT1u2RO72rr+escOHMx9oA3ZW08YYid+YJMROLCIyERgPHIibs7kLVbUxWcopS6vTJqCls5P5r7xC+wUXMPG++7o9P7MxhYrb83Y27kXD53EvIm4vYUzVp5yllTz6AB8UYWamSfnqq9MOydShuCWW84D/VtVvlzKYqtTa6rrtl6vyNgvPXw4Abpo+HayfiklR3A5yewGPJHFDETnMH5B7s4i8ISLfE5Gc04eISJOIaMTSnkRM3VLOFqEsMq0/xwPvTZvmxsE1JkVxSyztwIl0M7mIyL64IRdeBMbhBuX+AS7B/XuMS3wT1+SdsbY78XRba+uucWdTkkkq84H2QYPodf31qcZjDMRPLI8Al/sjxj2E6xTXhao+EOM604G+wHhVfRd4SET2BmaLyBX+tlxeUtX5MWMundZWN2xBynUrXZIKMHHt2l1N0VZJa1IUN7Hc4X82AWdH7FfcINv5jAV+F0og7bgOd8fiXh2obJkXDEvxLlCBLiKQVDIbM03RYMnFpCZuHcvBeZZh2U/t4lBgSXCDqq4ANvv78rlRRDpE5E0RuUpE+sa8b3KiBqVOyWXA/QSSSkZgkGxj0hCrxKKqSVUk7EvEYxSw3t+XzTbgx8Dvce8sHQf8K66OZlxCscWTcpd4D7iib18uHTyYwStX0pKt4ti67psUxe3Hkne8FVV9MeY9o/4SJMv2zLXfBL4a2PSYiKwCrhWRj6nqn3e7oMg0YBrA0CS7sZe7l23AzjqVLVsYfeutHHPMMa5OJaoC2brumxTFfRRaDLyQZ4ljPbBPxPb+RJdkcrnL//xE1E5VnauqzaraPHhwtimRCtTaGj1nTxkEK2rbBg1ySQVyzwttTEriVt5+NmLbAOAEf5kZ8zpLCNWliMhBwPsI1b3EoKHP0ps1y83ZU2Zdkgow6Uc/2rUzyyDZVnFr0hS3juXxLLvuEZH/BE4Dfh3jUr8FLhaRvXTXlCFfwM1NlO0e2WTqLBcVeF7xUqq3WIabwKkNmDRjxu5Jwx8k25hKkcTbzY8SfybEnwAXAHeLyOW41qTZwFXBJmgRWQo8rqrn+euzcb1/n8JV3h4DXAzc7Y8VUx5Dh5a1Q9x7QC/cKFuvAHv37g1HHVW2+xtTrLh1LLl8jpj1I+omPhuF6/PyK+C7wA+B/wgd2pOu/WKW4Pq53Ag8AJyBG8HujG7EXbio+owS8YCRwBX++t4QPXWqMRUobqvQLyI298bVlxwCxH450W89GpnnmKbQejuuH1i6gvUZJSy5eEALMI+IyitrRjZVIG6JZXDEsgfwR+AkVb28NOFVoMyQkCWSqaidR6hHbYY1I5sqELfyNqpVqL4NHJh4f5YO3HPlfKB9zz2Z2NnZtZevNSObKpFEHUt9uvrqxKdPbQDOwi+pTJniJihrbHT3aWzcNWGZMRXOxrztjl69EunX4gH/C3yKwMDXN9/sWoBK+NhlTKlYiSWuwKyCNDW5YRMSSipjcb0MuzxY2YuEpopZYokSTiLnn++GIli+3I0Wt3x5IvUrwR61NxAxV621AJkqZY9CYeHxVpYv7zp3T0J2G6Qp6iBrATJVykosYWUab+U68iQVawEyVcwSS1iZHj8yg/d2SSoNDdYCZGqCJZawEj5+eMCZwArcD/6TwZ39+rmWoM5O1xJkScVUMUssYSV6/MjUqbTjZn3rYuBAK6GYmmKJJWzyZPeHnqBwRe1J4QP23NOSiqkplliinHZaYpcKD9IUWVGb8txExiTNEkuUB+JMkRTPe/7SBkzKdlBDnJlTjKke1o8lSgItQxtxP9wBwNPkmXSpo6Pb9zOmkliJJSjT47abczFnHn9OI+ZMbo2N3bqfMZXGSiwZCc1wGBykqR03r0lO1hHO1CArsWQk0OM2nFR2VtQOHLir49uMGTYUgql5VmLJSKBe5UyyJJW1a7t9bWOqiZVYMhLocTsb+AWhJuV167p9XWOqjSWWjDlz3MBNcfmjx3m4qQMAPg6MDx9nbyibOmSJJWPyZNh77/jHi+xs/fkykHXiaquYNXXIEktQAY8tXmdnlx61h0UdNHCgVcyaumSJJSjmY0uXbvoNDdE9avv1cwNuG1OHLLEEzZkTa+T9J4Bn8bvp33zzrg5uma751oxs6pxoN3uZVoPm5mZduHBhvINzJBZlV4e3lcBB0O1eusZUMxFZpKrN4e1WYgnL0r3eA0YDD/rrB5UrHmOqkCWWsDlzoHfvLpsydSqP+V/vlPC4LcbUCkssYZMnww037Ewa4fFUdlbU9upllbPGZGGJJcrkybB2LZs3bWLsUUcxv6GB9gsuYFLwHZ8bb7TKWWOysHeFcujTpw8f+9jHuPDCC5k4caKVUIyJyRJLBM/zWLduHY2NjVxzzTVph2NM1bFHoRDP8xg7diwjR45k27ZtaYdjTFWyEktAJqnMnz+f9vZ29thjj7RDMqYqWYnFF04qEydGjqdvjInBEovvW9/6liUVYxJij0K+OXPmcPLJJ3PCCSekHYoxVc9KLL7+/ftbUjEmIWVPLCJymIg8IiKbReQNEfmeiOSdIUNE+ovIjSKyXkQ2iEiriFifemMqUFkfhURkX+Bh3IBr44APAj/AJbh/z3P6HcCHgC8BncDlwL3AZ0oUrjGmSOWuY5kO9AXGq+q7wEMisjcwW0Su8LftRkRGAGOAY1X1CX/b68AzInK8qj5cpviNMTGU+1FoLPC7UAJpxyWbY/OctyqTVABUdQGwzN9njKkg5U4shwJLghtUdQWw2d8X+zzfX/OcZ4xJQbkTy77AOxHb1/v7kj7PGJOCNJqbo8ZylCzbiz5PRKaJyEIRWbhmzZoCQzTGdEe5E8t6YJ+I7f2JLpHkO2+fbOep6lxVbVbV5sGDBxcSozGmm8rdKrSEUJ2IiBwEvI/oOpTgeVHNyofimpxzWrRo0VoRWQ4MAmwi5cpiv5PKU8jvJHKQ6HInlt8CF4vIXqqaGT72C8AW4PE8510qIker6pMAItIMDPP35aSqg/1zFkaNKG7SY7+TypPE76Tcj0I/AbYBd4vI8SIyDTeX+lXBJmgRWSoiP8+sq+o84HfALSIyXkROAVqBJ60PizGVp6yJRVXXA6OABuBXwHeBHwL/ETq0p39M0Om4Us0NwC3AIuDUUsZrjClOXUxYliEi01R1btpxmF3sd1J5kvid1FViMcaUhw2bYIxJXM0nlmKHaTClIyLDReR6EXleRDpE5LG0Y6p3IjJJRO4XkddFZKOILBKRLxZ7vZoeQa6bwzSY0vkw0IKbYLJ3nmNNeVyEe6n367g+LC3A7SIySFX/p9CL1XQdi4h8C7gEaMw0Z4vIJbgm7v2yDdNgSktEeqhqp//1XcAgVT0u3ajqm59A1oa23Q6MUNWDC71erT8KFTtMgymhTFIxlSOcVHzPAUOKuV6tJ5Zih2kwxsCncdUIBavpOhZsuAVjiiIio3D1kucWc36tl1ig+GEajKlLItIE3A7cp6o3FXONWk8sxQ7TYExdEpEBuBd7VwBnFnudWk8sxQ7TYEzdEZF+wK9xXQA+p6qbir1WrSeW3wJjRGSvwLY4wzQYU1dEpCdwJ3AIMFZVV3fnerVeefsT4ALcMA2X48ZvmU1omAZTXv7/jC3+6geAvUUkM2H2A6q6OZ3I6tq1uN/JTGCAiHwqsO85Vd1WyMVquoMcuC79wDXACFy9ys+A2arakWZc9cyvHFyWZffBqvpa+aIxACLyGllGg6OI30nNJxZjTPnVeh2LMSYFlliMMYmzxGKMSZwlFmNM4iyxGGMSZ4nFGJM4Syx1zp/j+pQiz71JRBaW8561QEReExGNWHYEjpktImv9ryf6+w/Pcr1mf/9p5foe8qn1nrcmv2nAYmJMVRvh+7hBs8p5z1pxOxAe8jHYqexnuLm3wL2/4+Hm1loUca3TgY3+cRXBEkudEpG+qrqlO9dQ1VeSiqcOvamq87PtVNW/A3/3v94qIvcCp4nIJRro1SoiApyGG+KgYl6FsEehEhCRx/yxXIPbjvOLqx/x15syxVd/xPoNIvJ3EfmuiPQInftPIvIrEXnHH0F9gYiMDuwf4F9jlYhsFZGnReSToWuoiFwkIj8SkTXAC/7o+IcDZweK41P9488SkSdFZJ2IrBeRR/35soPX7PIoJCJT/Wt8VEQeEpFNIrJERMYHfzZR9xSRK0XkVf8PJXiPc0Rku4gMKvw30ZV/r5kicpmIrBGR1SLyYxHZI3Tcx2TXzA7rRaRVRN4f2B/7d9eNWHc+CvnagKG4Ud2CjgIO8vdXDEss6bsCV4ydCNwGfMf/GgARORR4CtgfmI6bVvYe3D8m/D+Kh4HRwMXAKcAa4GER2S90r4v960zBvZx5Pm74iAdw71KNAH7jH9uEm8p2EnAG7n/PJ0RkWIzv6Xbgfj/Wl4F2ETnQ35ftnj8DDmb3sYinAr/KMiZrMb4BHIAba+RK4Cu4F+8AEJHBwGNAP9z3/TU/podEJDyjQM7fXR4iIj1DS65paR7CjZ5/emj76cA64Pcx71seqmpLwgvuH+ZdoW3H4Z6hP+KvN/nrt4SO+zPQHlhvw/1R981yr/OA7cAhgW09gVeAKwPbFPeWavj8hcBNeb6fHv41lwDfCWy/CVgYWJ/q3+fcwLaBwA5ger57Ak8CNwfWhwGdwOcT+r0o8ERo273A/MD6f+NeVt07sO1I/9wvFvK7yxHHa/754eWxwDGzgbWh864D3gIa/PUGf/36tP/NhxcrsaQv/D/Ni8CBgfWRwB2avT7keFyF3rLM/3z+9seB5tCxvyEmEflHEblHRFYBHcB7wIeAf4hx+s7vSVXfBlbT9XvK5ufABBHZ01+fCqwCHswRZ0Pwf/1CYvOFf95HAr/XwLAaqroAlwyOLuRaeUojtwFHhJav5Im9DXg/u0p1x/nrFfUYBPYoVAneCa1vB/oE1gcCb+Y4fxDwKdwffnA5B/9xKWBVnIDEDYz1e//8i4DP4P7hPx+KLZt3Quvh7ymbX+BKKKf5dS1n4UoFO3Kc8wqB71vckAzdiW1/on9Oq4ABca/lxxH8fYQrulep6sLQ8lKe2P+IK71mHodOB94AnshzXtlZq1BpbGX3Gf7C/yjjehv3jz2bdbhHixkR+8KD88QdI2ME7n/e0aq6cwhPEekf8/yiqOomEWnHlVSW48YHuSnPaScBwcrXN7oZxptEz6XzfqKberN5A5eMMwoaKCmKqqqI3AGcIyIXAuNxj44VN0+TJZbS+DtwTGjb6KgDY3gE9z/4LFXdmmX/CcAKLW44wajSRKZvys4/BhH5NK5uoZA/rkLumfFz3NSrs3F1H3/NdSFVfSGBeIKeAWaIyF6q6gGIyBG47/3JuBdR1e24hJ+0NlwF9JW4/6wq7jEILLGUyj3AeSLyQ1y9xmeBMUVe67vAs7gWmR/gSjAfB95W1RtwLTfTgcdE5P8Br+Ien44E3lLVH+a5/hLcuMBj/Gsvw/1hbwR+KiJX4Eovs4HXi/we8t7Tr4tBVZ8Rkf/F1Wfkq3Mohatwpb/fiRvOdE9che4LwC8TvM/+0nX4x4w/+UkpkqouEpG/+TG+oqrPJhhTYqyOpQRU9TfAt3FNj/fgivQXFnmtl3B/ZGtxTbL3+Ndd7u/fiktcD+GS0O+Bq3GDIi+IcYv/BP6Kq994FjhJVVfhmpn3A+7zY58OLC3me4hzz9D+e3EDnrcndL/YVHUN7ue5FVca+DGubmN0rj/4IpwBzItY4kxp2o6bG6vsP5+4bGhKU3FEZAHwkqpOSTsWUxx7FDIVw+/ZOxJX6fkvKYdjusESi6kkz+KacL9VqXUHJh57FDLGJM4qb40xibPEYoxJnCUWY0ziLLEYYxJnicUYkzhLLMaYxP0fC67O4jqz1CAAAAAASUVORK5CYII=\n", "text/plain": [ "<Figure size 432x288 with 1 Axes>" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAARUAAAEOCAYAAACn/4O6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAl6UlEQVR4nO3de5xcVZXo8d/qTiekI680CeK9dEfkoaAipEV5jVEuAtGRi90IGjLBVzSMGBIGvE4Qo077uaikGoEA4RlMJygIo2IAEcJrdGCCXD46GZTxYxIUyKMBgSTS0L3uH/tU+nTlnFOnqk7VOVW1vp/P+XTXqapTm9C9eu+191lbVBVjjElKS9oNMMY0FgsqxphEWVAxxiTKgooxJlEWVIwxibKgYoxJ1Li0G1BN++yzj06bNi3tZhjTcB5//PGtqjol6LmGDirTpk1j7dq1aTfDmIYjIhvCnrPhjzEmURZUjDGJsqBijEmUBRVjTKIsqBhjEpX5oCIivSLyKxEZFJG/icjvReQiERmfdtuMMbvKfFABOoA1wOeAU4AbgEXAkjQbZUzdGxiAadOgpcV9HRhI5LKZX6eiqtcUnFojInsA/ygi56oVhDGmdAMDMHcubN/uHm/YwN2f/SwnjozQOnt2RZeuh55KkEHAhj/GlGvRotGAAuSAU157jWvOO6/iS2e+p5InIq3ABOBI4MvAVdZLMaZMGzfu/DYHLAR6gc+/8ELFl66nnso273gYeBC4IN3mGFPHOjuBsQFlJdDW1VXxpespqBwDHA+cD5wKXBH0IhGZKyJrRWTtli1batk+Y+pHXx9/2W03LsIXUNrboa+v4ktLPY4gROQfgOXAgar6x7DXdXd3q91QaEyIgQGevOACDn3uOddD6euDWbNivVVEHlfV7qDn6ianUuA33te3AqFBxRizq/7+fiZMmMC8efM4PGYQKUU9DX/8jvW+/inVVhiTpjLWmeRyORYsWMCaNWuo1igl8z0VEbkb+CXwn8AwLqCcD/wwauhjTEMLWGfC3Lnu+5DeRy6XY+HChfT29jIwMICIVKdtqprpA/gW8DvgVeAl3NDnXKCt2HunT5+uxjSkri5V2PXo6gp8+ZIlSxTQ3t5eHRoaqvjjgbUa8nuX+Z6Kqn4N+Fra7TAmU3zrTOKcHxkZobe3l5UrV9LW1lbFhtXB8McYE6Cz0w15gs77bNq0iX333Zfzzz+fkZERWlqqn0at10StMc2trw/a28eeK1hn0t/fz8EHH8y6desAahJQwIKKMfVp1ixYtgy6ukDEfV22bGeStr+/nwULFvDhD3+Ygw46qKZNs+GPMfVq1qzAmZ58QKlVDqWQ9VSMaSA/+9nPUg0oYEHFmIZy8sknk8vlUgsoYEHFmIawfPlyNm/eTFtbG+edd15qAQUsqBhT93K5HGeffTbf+9730m4KYEHFmLrmX3rfl0DZgiRYUDGmTvX39+8MKGnmUApZUDGmDu3YsYNrrrkmcwEFbJ2KMXVHVZk4cSIPPfQQe+21V6YCClhPxZi6ksvl+OQnP8kbb7zBlClTdg0oVdrLpxQWVIypE/mk7BtvvBFcYClfY2XDBlcIIV9jpcaBxYKKMXUgH1B6enpYtWpV8JCnYC8fwD1etKg2jfRYUDGmFioYllxxxRU7Z3lCAwqUXGOlWiyoGFNtFQ5LjjjiCM4+++ziszwFtVSKnq8SCyrGVFuZw5InnngCgGOPPZYbb7yx+CxPjBortWBBxZhqK2NYksvlOPLII7nzzjvjf06RGiu1YutUjKm2mKUf8/xL70866aTSPiukxkotWU/FmGorYVjiDyhZWykblwUVY6ot5rDkySefrPuAAnWwl7KInA7MBqYDewK/B76nqquKvdf2Ujb1ZvXq1Zx44omZDyhReynXQ09lIW4jsQXAx4A1wEoROTfVVhmTkMsvv5yHH34YgJkzZ2Y+oBRTD4nav1fVrb7H94vIW3DB5vKU2mRMIvI5lDlz5nD88cen3ZxEZL6nUhBQ8p4Apta6LcYkyV8P5dprrx37ZAZuDCxXPfRUghwDrEu7EcaUK3IbjYEB+MxnYGjIPd6wwT2G1KeL48h8T6WQiJwAnApcmXZbjCmHqvLYY4+Fz/LMnz8aUPKGhtz5OlBXQUVEpgErgZ+o6k0hr5krImtFZO2WLVtq2Txjitq+fTsiws033xw+bTw4GPzmwcG6GBbVTVARkcnAXcBG4Kyw16nqMlXtVtXuKVOm1Kx9xhSTy+U44ogj2Lx5M+PGjStvlicD9VKKqYugIiLtwJ3AeOAjqrot5SYZU5L8LM+73vUu9t577+gXd3QEnxfJRL2UYjIfVERkHHArcBBwiqpuTrlJxoxVZEgSq8CS32WXQeFr2tpc7yRIjeulFKWqmT6AZYACXwbeX3BMiHrv9OnT1ZiqWrFCtb1d1f3Ku6O93Z1X1ZtuukkB7enp0aGhodKu29WlKuK+5h/7Pyd/dHVV4T8sGrBWw35nw57IygGs94JK0DEt6r0WVEzVFflF37Rpk1544YWlBZQwRQJYLUUFlczf+1MJu/fHVF1LS+Cw5KfAya+9xvjx45P9vIEBl0PZuNGVTujrS2XtSr3f+2NM7ZQ6ZRtQE6UfbyHVlVVYSjVrFqxfDyMj7msGF8NZUDEmr5xasgW1Uvpxd772HnUUX/rSl6rd4kyyoGJMXtxasv7ezKJFMGcOdHWNCSgrH3mk7u82LpcFFWPy4tSSDerNLF/O5gsvZPGee7ql900cUMCCijGj4mxxEdKbmXrxxfx60iRW3nYbbQcdlLlVrrVkQcWYvLBasjNnjg53CgpY54DvAAwO8o5nn6UNMrt8vlYsqBiTF1RLds4cWL58dLjjk8NVCvsPYKTwWhlcPl8r9VpPxZjqKNziYtq0XYc7jAaUXtxt84F/nbO2fL5GrKdiTJSAwJAPKD24gBKakq3xdqNZYUHFmCgBgWEScDqwCl9AERn7ohS2G80KCyrGRPElb//snZoL/JCCHsoXv5j6dqNZYTkVY6J4gSF37rksevFFfgW8BxjTL+nogKVLa9+2jLKeimluMe71yW3ezMIXX+QjRx3FYRMnjn2yvd3VPzGjwm5fboTDSh+YSDFKCSxZskQB7e3tdeULguqcNCHquZ5KJYcFFRMoquCRrxbKXXfdNTagmJ2igkro8EdEekVkt1r1mIypCf+9O2G8aeQTTzyRpUuX1vVm6WmIyqn8CNgkIstF5GQRaa1Vo4xJRFC+JOjenQI3TJ7MM888Q2trK/PmzbOAUqKooHIUcC0wA1gNPC8iS0WkMTZ8NY0trDZKVA8FVw/ls4ODLDnssPDkbR3svZOqsHGR/wCOA64AngOGgWeA7wLT47w/rcNyKk2oWL6ktTX0uZxX+7gXdCgseZuhOrFpIqlELa5ncwKuB7PVCzB/AL5RynVqdVhQaTLz5rlZmbCA4g8CBeeWhAWUwor1Gapon6aooFLSOhVVHVHV+1T188D+wFXA24CLKu0xGVORgQG4+urwvXHy8qtdfRt2vQb8gNGbAwMzKPl7gOIUcmpyJQUVcWaIyNW47UfPAZ4GvlWNxhmzi7B8xqJFxQPK+PHw6qsweza89BLgShZMAO4n5s2BcQo5NblYQUVEjhaRy4C/APcBM4GbgG5VfbuqLq5aC93nHygi14jIkyIyLCIPVPPzTEZFFaYu1lPo6HDvGRx0X4eH6Qc+huup7EVEQPHfHBhWyKlJbx4MFDYuAo4ALgH+hMudbAaWAseHvadaB27Hg2dw25/+F/BAnPdZTqXBROUzwp4TCUzehiZl/QndsFWztqq2vEQtrmf4V2A5cDLQGvbaah9Ai+/72yyoNKmwJGw+cBQmYEVc8rbgvUUDShPO5pQqKqhEDX9OB6aq6hxVvVtVhxPsIJVEVXep1meaUFQ+I6gU5A9+MHr3sPfepXjbaODlUFq9NZ35r01etiAJoUFFVX+sqq/VsjHGRJo5M/p81O59Xi7k/cDn8AJKe7urP6sKb7zhvmZ01796EnXvz2YROcL7fov3OPSoXZNN01q9urTzPo8eeCAsW8aRXV1cK0Kb9UiqJqpI05XAJt/3RebrskFE5uKKc9Fp03yNpcw1IrlcjoULF3LrrbfSu3598u0yY4UlW7J6YInaxhRnRiVshqe1NTSxmq+H0tPTY+ULEkREorbicpIiMg6X0H220muZJpVff5K/ezi//iRv0SLXG5k82S1gGxoa+/7h4dHX+4Yz+R5KT08Pq1atsruNayUs2gBDwHt9j1twCw8PKnjd+4DhsOskfWA9lcYT1gPp6Ai8Tyf0yPdwurp0HaiA9o4f76aNm3Q9SbVQ5pTyOMbW9xVcGYTdkw1rxYlIu1c0qhf4H8CU/GMRaS/2fpNxYTmRwcGitU/G8JU3eAfwS2Dl0NDoVqRnneWmm/PL+62EQVXUSzX9qbjVtH75x28F1te0NSZZnZ1F65zE0trKZdu383bgJOBDYa/bsAE+/WkXYPJDKf+Qy2aEKlIX1fRVdb2qSsixPu32mQrNnLnrZlylGj+e/uFhzgNWxHn966/vmptp4v2Pk1QXQcU0sIGB0QVoFehndKXsDZVcyEoYVKzY8OdcEXnO+z7/p2S+iGzyvWa/5JtlmkaMmrHF9AMLhoai66HEZWubKhYVVDbiykj6bQD+LuS1xsSTL0C9cWPFPRQF1lGkwFKQtraxORWwEgYJCQ0qqjqthu0wzaJwTUoFXsFNRV6Nu6U+9qxDV9do8MgHt85Od86StBWrl9kf0ygSGO4A5IDLgF8Bb6HM5OCsWRZEqsAStaa2EkiE5oCFQDcwpZwL+CvGmcRZUDG1VWEiNB9QeoBVVJCUtenjqrGgYmorqMZrzHtyVlJGQIla/2LTx1VhQcXU3sSJo99PmgQTJsR620zcXjCxA0q++ltYYLHp46ooK6iISKd3d7Ix8eVnfgYHR89t2+a2zYjwY2AHruL9tyhhyNPXF751h4hNH1dJyUHF26j9T8C7k2+OaUj5G/fOOqvkmZ8cbg1Kf6mf2dHhZnbChjiqNvNTJeUOfyq8UcM0Df9ePSXyJ2X/qZQ3trXBZZe578OGOF1dJbfHxGM5FVNdZa5LKXuWp6MDbrxxtBdim3/VXDl5EcUt17dK+6a4MmZYBoFvU2JAmTQpODeTDy62crZmRCu89yLLuru7de3atWk3o7nk7+vZsMHtpTNc3nZRfwQ6KaGHIuK25jA1ISKPq2p30HM2g2OSMTAA8+ePndkpMaDkcL2UbwFvK/XzbXo4MyynYioXNFVconwO5SnczYElsRxJplhQMZWr8CbBwqRsa9CLWrwf1a4umDdv7PamYZuCWQ3aVNjwx1Sugvqy/cSc5RkZGe2RxEmyRm37YUnaqiraUxGR3UTkWhF5fy0aZOpMhX/9pwKfIOYsTyk3AQb1nuwmwpooGlRU9W/AmcBu1W+OqQvnnAPjxrnhx1lnlXWJ9d7XTwG3UMIsT7Ep6vyQJ6z3ZDcRVl3cnMr9wAer2RBTJ845B666quypYnA5lEOAx7zHJS3PjprlibN612aJqi5uTuVK4DoRmQSsxm3cPmaBi6quS7htJouWLavo7fmkbC9wRNQL581zVfb9Q5hiszzFEsY2S1QbYVsX+g/cLJ//GPYdI1R521PgUOA+YDvwLPBNoLXY+2zb0yqIuwVpwLHE/SHSXnBbkUZtX6oab9N2P5HiW6KaRJDABu2pDX1EZG/cDpbrgFNx66IuxQ3dLkqrXU2rzFWy9zM6yxNZ9d7fmyi1hmzYToddXbB+fQmtNZWIFVRU9cFqNyTCF4GJwMdV9WXgXhHZA1gsIt/xzplamTvX5VRK9EHgemA2RZKy/gJOperr27VSvw15aq6kxW8i8j4ROV9E+ryv76tWw3xOAe4pCB634ALNB2rw+cZv6VKX72gNXKK2i2XAf+OSsZ8hxizP4GD5RalnzXI5nzgL40z1hI2L/AeQT9COAEPAc97XYeDnQHuc65RzAJuBxQHntwEXRL3Xcio1sGKFaktLZA7lH8vJweTzKiaTiMipxO2pfAc4GjgD2E1V98OtWznTO39JEgEuxN7ASwHnX/SeM2latCjw7mD/LE+unOvaepK6FTeo9ABfUdVbVXUEQFVHVPVW4P8Ap1ergZ6g+gwSdF5E5orIWhFZu2XLlio3q05Vck9M4XsDEqP+gFL23sa2nqRuxQ0qewLPhDz3DLBHMs0J9CKu5nGhPQnowajqMlXtVtXuKVPK2mqqsfkXiKmWtrFW0HsLvA7cRoyAErV1hiVX61rcoPIkME9k7E+C93ie93y1PAW8veBz98fleZ6q4uc2pnLuiYlZuPoNXBC5hxg9lEmTdi3zCK4cpCVX61rcdSr/DNwFPCUid+BW1E4FTgOm4WZoquUu4AIR2V1VX/HOnYHbtSHNqe76FJarCDt/zjlw9dXB21z49AM/wWXt3xSnHa++6gLIxInwwgtW5rGBxOqpqOr9wHuAJ3D5kz7czaW/AY5U1TXVaiBwNa4e7u0i8r9EZC6wGFiitkaldGG5iqDzAwOxA8oCYB9KzJ8MDsKOHW7Dr/XrLaA0irBpofwBTAAWAYcXe221Dtwy/ftxvZPncBUHbZl+OVasUG1vHzt9294evIS9o6Po1G8u7tJ7mz5uKFQypayqr3lBZa9qBLU4VHWdqn5IVSeq6n6q+jVVLf822WYWtUDMP7Ozzz5Fy0NejeuhVDTLAzZ93GDi5lQeBaZjOYzGEHRPTWGltBj1Zj+Au4fi+1QQUMCmjxtM3KByIbBSRIYIL31QfpFSk74S6sw+BBwPvAMo/S6gAjZ93HDiTik/irs7+PvA08DLwCsFh6m1JAs7x6wzm8P1UMr6pJaW+EWrTd2K21P5dFVbYUpXTmHn/EZfhTv1DQy4X/Iiszz+qvdnlNPmkRF3Q6JpbGEZ3PxBBmZ/yj0aevanq6u0mZSoWZ+wa5VTYKnYYRoCEbM/sbY9FZFtwExNt65KyRp629OWluCeRdj2n2HFoLu6XM8l4ufgv3H5k/9NhbM8tjVpw4ja9jRuTuUx3OyPyYpSFrFB9EraIrMvBwIPUGFAgaLDK9MY4gaVC3H3/nxJRA4QkUki0u4/qtlIE6Cvb9d7Z6JmUqKCUNC1gMuAH3vfH0uFAcU0DZv9qVelVjkLChzjx7t7cGbPdvfgdHTsfKofOI/RoJIY23q04cXNqZxNcE2TnVR1eUJtSkxD51TiKJztmTkTVq92jydPhpdfhtdfH319ezvMmUP/VVcls1I2iBWhbghROZXUZ2iqeTT07E8xxe7xCZnxye29dzKzPGGHSOlbb5jModLZH190OhSXsN0fuEFVnxeRA4FNOlqWIDOauqcSNduzfn3o7NH5wEaq0EPJ6+hwdyb7V++KwBe/aGtY6kjFsz8i8iYR+RHwO+A63F3Cb/Ge/jbw9SQaahJUrG5KQeL2Je/r94i5WXo58jmdwtsBVF2JBcu3NIS4idolwDHACcDujN3+djVwcsLtMpWKmu0ZGHAJWk8Otw5lA+5/bNxl1iXJJ5JfeCH4edXo6nOmbsQNKh/HFb5eg9uWw28D0JVoq0zlwqacZ850y/m9u5D7cUvvj2O06xlLRwe0FfRnWlrcjFLhZ65YMVqEKWpNjJVAaAhxg8pEIOxe+N3ZNdCYtPmnnMFt/rV9uzvnDT/6qaAeymWXwY03jp3SvvlmuOGG6Gnuvr7wotdWAqExhGVw/Qfegkrv+1bcpmJHeo9vBlbHuU6tj6af/Ymo3PbDSu/lqcS8ebtuph5Wfc5kEglsJnYR8HER+SXwOdwP5EwR+QGuZq0larMgXwpBxFW+jyi09BHgm1RxlifK0qWuLq2VQGhMYdGm8MCt1H4Yt7XLCG7I82/AsXGvUeujKXoqMe8yzh+rQF9OYr1Jrf/7bE1LphDRU4md6FfVfwOOF5GJeFuRqlV7S1dhTZUicrik7GICupZtbWNX10aZNCl2EytSTs0Yk7q4w5+dVHWHqj5rASUD5s8vOaD04jZxGqO1FSZMiPeZInDNNfHbWIlyNj4zqSs5qJiMOOecWMWpIcbexsPDY9athOrocLmQWvUSSt34zGRC5oOKiJwhIreLyHMiot7Njc3Jn4i9Kl7J6b8Cl+ILKCKuZxJXR8doJmXr1toOO0qtGWMyIfNBBff7MA24M+V2pMu/OXpMitvF/tf4eigjI6VVX3vllfSWz5daM8ZkQj0ElTNU9UjcOq3mVcIWGuCGPOfhAsv+FAx5wv7StwT8OAwNpZfDKLVmjMmEzAcVVbWiplBSHiGfQ/kLIUudw3oAYT2YNHMYs2a5Jf4jI7bfcp3IfFAxnph5hHxA6cHdbTxmzUB+eXxYDyC/pL/MzzYGLKjUj5A6sn6XMTrLE1i+QNUles85J3j/H8thmASUVKQpkQ8U2RPYr9jrVPWpgve9CVcL99OqelPE9ecCcwE6OzunbyghsZl5/vKQ7e2wbduYp38K/BC4iRKX3o8fD7vv7soSTJ7szr3wwtiAY4xPVJGmNILK54Bri71OVcfcyho3qPg1dOW3cePc+hLgD8DBSV67vd0SoiZSEvv+JEZVr1NVKXbUul11xwsoOeBQ3KbpO3V1je5ZXA5btWoqYDmVehC0EXtr686k7GnA0fnXtra6WZKlS93XcgOLrVo1Zcp8UBGRQ0WkF7frJkC3iPSKyAdSbFbt+Be9qe68qS53wAHBS+/zN9zlxUjwBrIZH1OusNuXs3LgbqrVgOOBYu9tiNIHAWUNHs4XWDrgAB1qaXHnW1td8aMgheUD5s0bfdzRodrWNvYzrGCSKYKI0gepB41qHqkElaTrfwTUMhkBHchXbCv8jHI+32qWmBJFBZWaz/7UUs1nf4Lqm1Q6k+Kb5bkKOB54Z+Fr8p8ByX++MQEyNaVcSzUPKsU28CqHtwo2n5SdCwRWM8knZJP+fGMCZGpKuaFVo/5HV9eYpfdXRH221R8xGWBBJUlVqP+RO+64MffyhK6U7ewM/xzV0anoUgVNZxsTwYJKkhK+d2Z4eJjVmzbR8973sqqz0xVY6ugI3rCrr89tFBa2p06+vmspQSFkOtsCi4kUlsFthKOeZ3+GhoZUVXX79u07v4/8jBUr3FRwsSr4XV3xGxFWpb+Ua5iGhM3+ZJT/BkHfjXz9e+3FLZMn84vf/IY99tgj3rXCksSFROJXfmtpcWGkkmuYhmSJ2iwqHFoMDsLgIP2qLHjxRfZfv56Jd9wR/3pxk7Gl5HesRqwpgwWVtASUh8zh29t4eJi2r5ew8WOcX/RS8ztWX8WUwYJKWgp6FtcSsI1GKVPBQQGgrc0ldsut72o1Yk0ZYu9QaBLW2TkmB3IC8CVgCb5p41KGGflf9KCKbpWYNcuCiCmJ9VTS4vUsfonbmPoA4HJ8AaWcYYYViTYZYEElLbNmkTvtNE4EbgA3TKlkqGJMRtjwJyW5XI6FAwP09vYyZ+VKl/8wpgFYTyUFuVyOhQsX0tvby8qVK2mzgGIaiAWVGtuwYQNf/epXLaCYhmXDnxrr6urikUce4fDDD7eAYhqS9VRqpL+/n5tvvhmA7u5uCyimYVlQqYFcLseCBQv4+c9/TiPfa2UMWFCpOn9SdsWKFUhYaQJjGoQFlSrKB5Senh5LypqmYUGlil555RV6enpYtWqVBRTTNDI9+yMiewDnA6cAhwA7gF8DX1HVP6TZtihbt25ln3324eKLL2Z4eJjW1ta0m2RMzWS9p9IJfB64B3cD7xeA/YBHRWT/NBsWpr+/n0MOOYSnn34awAKKaTqZ7qkAfwLepqo78idE5GFgI/AZ4BtpNSxIf38/CxYsoLe3l2nTpqXdHGNSkemgoqrbAs69ICIbgKkpNCmUP6BYUtY0s6wPf3YhIlOAA4F1abcl74477rCAYoyn7oIKcCnwKnBL2g3JO+WUU7jkkkssoBhDCtueisieuGRrJFV9KuC984ArgR5VDawKLSJzcbuD0tnZOX1DnArzZVqxYgUzZ85kcr4SvjFNImvV9E8H/ivGMYaIfAxXHO0rYQEFQFWXqWq3qnZPmTKlCs13crkcs2fP5tJLL63aZxhTj2oeVFT1OlWVYof/PSJyDG64c7WqfrfWbS7kX3q/ePHitJtjTKZkPqciIocBdwJ3A1+uyYdG7B9sBZaMiZbpKWURmYoLJq8C3weO8t2Q97KqJjsDNDAA8+e7jb3y8vsHA6+eeipXXHGFBRRjImR621MRmQGsCXn6QVWdEfX+krY9ze8YWLDBV552diIbNvD888/T0dFhAcU0tahEbaZ7Kqr6AFCbWgEBOwbm5YDfbtzIdSMjvPnNb65Jc4ypV5nPqdRMyG6AOdzOgS+3tzM8PFzTJhlTjyyo5AXsBpgPKD2traxautSGPMbEYEElr2Av4u/j7W08fjyrrr+etjlzUmuaMfUk0zmVmirYi/iwKVOYffDBXH///dZDMaYEFlT8Zs3id4cfzjvf+U5OwG2abowpjQ1/fHK5HO9+97u59957026KMXXLgorHX6R6xowZaTfHmLplQQVbem9Mkpo+qKxdu9YCijEJavpEbXd3N7fffjsf/ehHLaAYk4CmDyoAp512WtpNMKZhNP3wxxiTLAsqxphEWVAxxiTKgooxJlEWVIwxibKgYoxJlAUVY0yiLKgYYxKV6cLXlRKRLUD1tigctQ+wtQafU2/s3yVcvf/bdKlq4G59DR1UakVE1oZVFm9m9u8SrpH/bWz4Y4xJlAUVY0yiLKgkY1naDcgo+3cJ17D/NpZTMcYkynoqxphEWVBJiIjsISLfEJHHROSvIvK8iNwhIgen3bZaE5FDReQ+EdkuIs+KyDdFpDXtdqVJRE4XkZ+KyF9E5FUReVxEPpl2u6rBgkpyOoHPA/cAvcAXgP2AR0Vk/zQbVksisjfwS0CBU4FvAucD30izXRmwEHgVWAB8DFgDrBSRc1NtVRVYTiUhIjIJGFHVHb5zk4GNwHdVtSl+qUTkq8CFuMVRL3vnLgQWA2/On2s2IrKPqm4tOLcSOFpV35pSs6rCeioJUdVt/oDinXsBt6J3ajqtSsUpwD0FweMWYCLwgXSalL7CgOJ5ggb82bCgUkUiMgU4EFiXdltq6O3AU/4TqroR2O49Z0YdQwP+bFjh6+q6FDeOviXthtTQ3sBLAedf9J4zgIicgMs5fSbttiTNgkoEEdkTl2yNpKpPFZ4TkXnAWUCPqg5WoXlZFpSok5DzTUdEpgErgZ+o6k3ptiZ5FlSinQ5cG+N1MuaByMeAy4GvqOod1WhYhr0I7BVwfk+CezBNxUve34VL4J+VcnOqwnIqEVT1OlWVYof/PSJyDG64c7WqfjedlqfqKQpyJ96U+iQKci3NRkTagTuB8cBHVHVbyk2qCgsqCRKRw3A/NHcDX065OWm5CzhJRHb3nTsD2AE8mE6T0ici44BbgYOAU1R1c8pNqhpbp5IQEZkKPI7LG/wD8Dff0y+rasNl+YN4i9/WAb8DLgEOAJYA/ap6UZptS5OILMMtjpwPPFbw9BOq+lrtW1UdFlQSIiIzcKskgzyoqjNq1piUicihwBXA0bg8ynXAYlUdTrNdaRKR9UBXyNNvVdX1tWtNdVlQMcYkynIqxphEWVAxxiTKgooxJlEWVIwxibKgYoxJlAUVY0yiLKgY4xGR9SKiAccbvtcsFpGt3ve93vPTQ67X7T3/iVr9N2SB3VBozFgrcTeD+vkXc10H/Mz7/k7gFeBM3GrqQmfiSl/cmXAbM80WvzUREdlNVf8W93zMa04srHhXr7xVr7ep6j+V8J6bcRXtpqnvl0lEBFf17yFVbci7kcPY8KdOichxIvKgV7F+UESu9d/EJyJne13vo0TkARHZAVwQdt57z4dE5FER+ZuIbBKRpSLyJt81Z3jvPcmrDP8qbjl+pf8tKiLzReTbIrJFRDaLyJUiMqHgde/xVel/UUQGRGRf3/PT8sMNEbnG29Xgz94uB4n8rPuHP55VuKLnxxS89Fhgf+/5pmJBpQ6JyLHAfcDzuMr95wEzgRsDXr4K1/2eydhu+Jjz3v06dwNbgR7g68CngNsCrnk98CSuKvz1Ff8HOecDb8HVGPkubjeC+fknvdKcDwDtXrvOxfUQ7hWR8QXX+g5u2NELrAAu9r6PQ0RkXMERtb3Ivbh/szMLzp8JvAD8IubnNg5VtaPODuBhYE3BuQ/hxv7v9B6f7T2eX/C6sPO3AE8Drb5zn/Bee7T3eIb3OJfwf4/ihgn+c/8K/Lvv8f/F3Zy4h+/cUd57P+k9nuY9vrngWv8PuCVGO9Z77y88HvC9ZjGwteB9V+ECfKv3uNV7fE3aPytpHNZTqTNeoZ+jgR/5/5oCjwCvA4UzET8PuVTh+aOAO3TsncQ/Bt4Ajot5TX87WwvaV0zhX/R1wP8saN8v1FelX1UfwwWCwvZFXqtIL2QF8N6C4wtF2r4K2JfR3QJmeI+bbugDNvypR3vj/hIuxQWR/PEa0IYbx/ttCrlO4fn9Cs95AWYQmBzzmn5/9LfPq8sa5aWCx0PAblHt87WlsH2h1/La4f93+2Ph9VR1bcHx+yJtfxj4M6NDoDOBZ4GHiryvIdmUcv15CdclXwysDnj+2YLHYdN7heefo2APGu+veAcuNxDnmn5/D/gTrYXtKtUu7fPsS/B0bphncb2PvIqLI6mqisgPgU+LyHnAx4HlqjpS6bXrkQWVOqOq20Tk34FDVPWbCV76UeA0Efln3xDo47ifkUfKaOdvE2wbuPbNE5HdVfUVABF5Ly6PErt9qjoErE24beCGOufjksyTadKhD1hQqVcXAveJyAhuduYV3LTmR4BFqvqHMq75L7gd8/5VRK7C5SAuwe02+Otkml2RJcA84B4RuQR4Ey55+1tc7icp+4nI+wPO/8YLSIFU9XER+YPXxj+q6n8k2Ka6YjmVOqSqjwB/B0wBfoBb4Xkh8Azx8h1B1/xP3JalU4HbcUFmFfGnYqtKVbcAH8TV/l0FXInLZZwY9ctehk8Bvw444mxPegtuu5Zm2jxuF7ai1hiTKOupGGMSZUHFGJMoCyrGmERZUDHGJMqCijEmURZUjDGJsqBijEmUBRVjTKIsqBhjEvX/AXKN0IQrYqKxAAAAAElFTkSuQmCC\n", "text/plain": [ "<Figure size 432x288 with 1 Axes>" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# plot and save figures\n", "plt.figure(1)\n", "plt.clf()\n", "plt.scatter(ber_unc, unc, color='r')\n", "plt.plot(*diagonal(unc), color='k', linestyle='dashed')\n", "plt.xlabel('uncertainty - non-EiV')\n", "plt.ylabel('uncertainty - EiV')\n", "plt.axis('square')\n", "plt.gca().set_aspect('equal','box')\n", "plt.tight_layout()\n", "plt.savefig(os.path.join('saved_images', 'wine_unc_scatter.pdf'))\n", "\n", "plt.figure(2)\n", "plt.clf()\n", "plt.scatter(ber_err, err, color='r')\n", "plt.plot(*diagonal(ber_err), color='k', linestyle='dashed')\n", "plt.xlabel('error - non-EiV')\n", "plt.ylabel('error - EiV')\n", "plt.axis('square')\n", "plt.gca().set_aspect('equal','box')\n", "plt.tight_layout()\n", "plt.savefig(os.path.join('saved_images', 'wine_error_scatter.pdf'))" ] }, { "cell_type": "code", "execution_count": 83, "id": "4c5e40a4", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.95228493\n", "1.0843716\n", "1.6982132\n", "2.5368009\n" ] } ], "source": [ "last_std_y = stored_std_y[-1].item()\n", "ber_last_std_y = ber_stored_std_y[-1].item()\n", "print(np.sqrt(np.mean(np.abs(err)**2/(unc**2 + last_std_y**2))))\n", "print(np.sqrt(np.mean(np.abs(ber_err)**2/(ber_unc**2 + ber_last_std_y**2))))\n", "print(np.sqrt(np.mean(np.abs(err)**2/np.abs(unc)**2)))\n", "print(np.sqrt(np.mean(np.abs(ber_err)**2/np.abs(ber_unc)**2)))\n", "#print(np.corrcoef(np.abs(err),unc)[0,1])\n", "# print(np.corrcoef(np.abs(ber_err),ber_unc)[0,1])\n" ] }, { "cell_type": "markdown", "id": "statewide-aurora", "metadata": {}, "source": [ "## Deming vs RMSE\n", "Produces Figure 5a of the preprint" ] }, { "cell_type": "code", "execution_count": 84, "id": "smart-polls", "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "b23ed94b293e4a7da02cb3ab583524a1", "version_major": 2, "version_minor": 0 }, "text/plain": [ " 0%| | 0/11 [00:00<?, ?it/s]" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Cycle through seeds and deming factors\n", "# and collect the RMSE\n", "\n", "set_seeds(0)\n", "\n", "# This deming factor will be considered for error/unc quotient and coverage\n", "deming_for_quotient = 1.0\n", "\n", "rmse_list = []\n", "\n", "# quotient error vs uncertainty/uncertainty+noise\n", "unc_q = []\n", "full_q = []\n", "ber_unc_q = []\n", "ber_full_q = []\n", "ber_rmse_fixed_seed_list = []\n", "\n", "# coverage of noisy labels by uncertainty/uncertainty+noise\n", "unc_cov = []\n", "full_cov = []\n", "ber_unc_cov = []\n", "ber_full_cov = []\n", "def compute_coverage(errors, uncertainties, k=1.96):\n", " return np.mean(np.abs(errors) <= k*uncertainties)\n", "\n", "\n", "for i, deming_scale_test in enumerate(tqdm(deming_factor_list)):\n", " rmse_fixed_seed_list = []\n", " for seed in seed_list:\n", " net = Networks.FNNEIV(p=0.5, init_std_y=init_std_y, \n", " precision_prior_zeta=precision_prior_zeta,\n", " h=[dim, 200,100,50,1],\n", " deming=deming_scale_test)\n", " saved_file = os.path.join('saved_networks', \n", " 'eiv_wine_init_std_y_%.3f_deming_factor_%.3f_seed_%i.pkl'\n", " % (init_std_y, deming_scale_test, seed))\n", " net.to(device)\n", " train_and_store.open_stored_training(saved_file, net=net, device=device)\n", " net_train_state = net.training\n", " net_noise_state = net.noise_is_on\n", " net.train()\n", " net.noise_on()\n", " pred, _ = [t.cpu().detach().numpy()\n", " for t in net.predict(test_x.to(device), number_of_draws=300,\n", " take_average_of_prediction=False)]\n", " val_y = test_y.detach().cpu().numpy()\n", " if deming_scale_test == deming_for_quotient:\n", " err = np.abs(np.mean(pred, axis=1).flatten()-val_y.flatten())\n", " unc = np.std(pred, axis=1).flatten()\n", " last_std_y = stored_std_y[-1].item()\n", " unc_q.append(np.sqrt(np.mean(err**2/unc**2)))\n", " full_q.append(np.sqrt(np.mean(err**2/(unc**2+last_std_y**2))))\n", " unc_cov.append(compute_coverage(err, unc))\n", " full_cov.append(compute_coverage(err, np.sqrt(unc**2 + last_std_y**2)))\n", " r = np.sqrt(np.mean((np.mean(pred, axis=1).flatten()-val_y.flatten())**2))\n", " rmse_fixed_seed_list.append(r)\n", " if i==0:\n", " # non-EiV\n", " ber_net = Networks.FNNBer(p=0.5, init_std_y=init_std_y, h=[dim, 200,100,50,1])\n", " ber_saved_file = os.path.join('saved_networks',\n", " 'noneiv_wine_init_std_y_%.3f_seed_%i.pkl' % (init_std_y, seed))\n", " ber_net.to(device)\n", " train_and_store.open_stored_training(ber_saved_file, net=ber_net, device=device)\n", " ber_net_train_state = ber_net.training\n", " ber_pred, _ = [t.cpu().detach().numpy()\n", " for t in ber_net.predict(test_x.to(device), number_of_draws=300,\n", " take_average_of_prediction=False)]\n", " ber_r = np.sqrt(np.mean((np.mean(ber_pred, axis=1).flatten()-val_y.flatten())**2))\n", " ber_err = np.abs(np.mean(ber_pred, axis=1).flatten()-val_y.flatten())\n", " ber_unc = np.std(ber_pred, axis=1).flatten()\n", " ber_last_std_y = ber_stored_std_y[-1].item()\n", " ber_rmse_fixed_seed_list.append(ber_r)\n", " ber_unc_q.append(np.sqrt(np.mean(ber_err**2/ber_unc**2)))\n", " ber_full_q.append(np.sqrt(np.mean(ber_err**2/(ber_unc**2+ber_last_std_y**2))))\n", " ber_unc_cov.append(compute_coverage(ber_err, ber_unc))\n", " ber_full_cov.append(compute_coverage(ber_err, np.sqrt(ber_unc**2 + ber_last_std_y**2)))\n", " rmse_list.append(np.array(rmse_fixed_seed_list))\n", "\n", "# convert list to an array \n", "rmse_list = np.stack(rmse_list, axis=0)\n", "ber_rmse_list = np.stack(ber_rmse_fixed_seed_list, axis=0)\n", "\n", "# restore settings of last loaded networks\n", "if net_train_state:\n", " net.train()\n", "else:\n", " net.eval()\n", "if net_noise_state:\n", " net.noise_on()\n", "else:\n", " net.noise_off()\n", "if ber_net_train_state:\n", " ber_net.train()\n", "else:\n", " ber_net.eval()" ] }, { "cell_type": "code", "execution_count": 85, "id": "amino-expression", "metadata": { "scrolled": true }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ0AAAEOCAYAAABSLcpPAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAx80lEQVR4nO3deXxU9b3/8dcnIaxGdgStGKF1r0VMq1YFFVGpbVG0Qhdva2vVVmvtdbdqKdq69Ad661LxWrGtG7fetLeoKIKCYl0KaF1RURCRHQMBEsj2+f3xncg4mUkmITlnkryfj0cek/M933PmM+cx8Mk5383cHRERkSjkxR2AiIh0HEo6IiISGSUdERGJjJKOiIhERklHREQio6QjIiKR6RR3AG1Nv379vKioKO4wRERy2sKFC9e7e//UciWdJioqKmLBggVxhyEiktPM7MN05Xq8JiIikVHSERGRyCjpiIhIZJR0REQkMko6IiISGSUdERGJjJKOiIhERklHREQ+q6QEhg2DXr3Ca0lJi51ag0NFRGSHkhI491woL4f8fHjvvbANMG7cTp9edzoiIrLDpEngDlTA9i2wbVvYnjSpRU6vpCMiIjssWwpeBp0dunkoKyiAZcta5PSRJx0zO8DM5phZuZmtNLNJZpbfyDETzcwz/FyZUrevmU01s9VmVmFmi83sP5L2F2U4z8Ot9ZlFRNqEze9Dn+1QWQUOlFsor6qCFproONI2HTPrDcwG3gLGAkOByYTkd3UDh94DPJFSdgpwOTAz6fy7As8CW4CfAeuBA4DOac55CfB80vb67D+JiEg7s3Y+PHcKfHM7TMuHLbVQC+BgBtde2yJvE3VHgvOAbsA4dy8DnkokiolmdnOirB53XwGsSC4zs2uAxe7+alLxVUAXoNjdKxJlz2SI5R13f7H5H0VEpJ1Yej+89COorYSxJ8GR34OLLoNPPoE+feC221qkEwFEn3TGAE+mJJeHgZuAkcCMbE5iZn2A0cD1KbvOAm5NSjgiIpKJ18Jrv4I3E/+V7nMBDL8F8jrBt77bKm8ZdZvOfsDi5AJ3Xw6UJ/Zl63SggJCwADCzvYEBwEYze9zMKs1snZlNMbN0j9emmVmNma1K1OnW5E8jItJWVVfA898OCcfy4NDboPi2kHBaUdR3Or2BjWnKSxP7sjUBWOTu7yaVDUy83kxIRicBXwJ+C1QDlyX2bwfuAGYBZcAxhLahoYR2JhGR9q1iDTw7Fja8BJ0K4ajpsPuYSN46jsGhnqbMMpTXr2g2iPAo7vKUXXV3bW+6+48Tvz9tZoXAVWY20d3L3X0VcEHScXPNbA1wp5kNS2kjqnvPc4BzAAYPHpxNmCIiuWnj6zD361C+HHrsBSMfhV4HRfb2UT9eKwV6pSnvSfo7oHTOICSp6SnlnyReUzsOPE3oXDC0gXM+kngdnm6nu9/t7sXuXty/f70lv0VE2oaVM2HWkSHh9D0MTngp0oQD0SedxaS03ZjZnkAPUtp6GjABmO/uH6WUvw9Upqmf6GgeOv9l4CmvIiLtyzu3wbyvQ/VmGDweRj0D3XaLPIyok85M4MTEI68644EKYF5jB5tZEXA48FDqPnevBJ4CjkvZNYrQUWFJA6c+PfG6sLEYRETalNpq+NcFsPDC0FvtoGvgyAehUzx9p6Ju07kLuBAoMbObgCHARGBKcjdqM1sCzHP3H6UcP4HQKeAR0psEzDezaYTEdDBwBXCdu29PnHsiUEgYGFoGjAAuBUrc/bUW+IwiIrmhqgzmj4dVT0BeZzjsj7D392INKdKk4+6lZjYKuJ0wJmcjcAsh8aTGlW5qnAnAHHdfl+H8L5vZN4AbgO8Aa4HfJLbrLCbMRnA2YaDqcuB3iXoiIu3DlmXhcdqmN6FLPzj6bzDgqLijwtzVjNEUxcXFvmDBgrjDEBHJbP2LoUv0trWw635wzGOwy5BIQzCzhe5enFqu9XRERNqTZQ/Diz+A2u0w8Hg46q/QuVfcUX1KSxuIiLQH7vD6dfDPb4eE8/lz4ZjHcyrhgO50RETavprt8NLZsOx+wGD4ZNj3ojA7dI5R0hERacu2rYPnToV1z0OnHvDVh+Bz34g7qoyUdERE2qpNb4ceals+gO6fg5EzoPewuKNqkJKOiEhbtHo2PHc6VG2CPsUw8h/QbVDcUTVKHQlERNqa96bCMyeFhLPnODh+XptIOKCkIyLSdtTWwML/hH+dB14DB1wRukR36h53ZFnT4zURkbagagv88zvw8QywTvCVu2HoWXFH1WRKOiIiuW7rRzDvG7Dx39C5NxxdArsdE3dUzaKkIyKSyzYsgGe/CRWroPALYdG1XfeJO6pmU9IREclVH5XAP78HNRUwYGS4w+nSJ+6odoo6EoiI5Bp3ePNGeO60kHCGnAXHzmrzCQd0pyMikltqKkPvtA+mhe1hN8L+l+XklDbNoaQjIpIrtm8Idzdr50F+N/jq/WEcTjuipCMikgvK3g1T2mx+Lwz0HPEP6FtvOZo2T0lHRCRua+bCc+OgsjTMnTZyRphLrR1SRwIRkTi9Pw2eOSEknD2+Acc/124TDijpiIjEw2vh1SvgpR9CbRXs959w9N+gYJe4I2tVkScdMzvAzOaYWbmZrTSzSWaW38gxE83MM/xcmVK3r5lNNbPVZlZhZovN7D9S6vQ0s2lmVmpmm8zsATPr2xqfV0SknupymP8teOsmsHz4ytSw8Fpeg/8VtguRtumYWW9gNvAWMBYYCkwmJL+rGzj0HuCJlLJTgMuBmUnn3xV4FtgC/AxYDxwAdE45djqwL3A2UAvcBPwdOLrJH0pEpCnKV4YZBj5ZCAU94ehHYODxcUcVmag7EpwHdAPGuXsZ8FQiUUw0s5sTZfW4+wpgRXKZmV0DLHb3V5OKrwK6AMXuXpEoeybluCOAE4GR7v5souxj4CUzO97dZ+/shxQRSav01TCHWvkK2GVImNKm5/5xRxWpqB+vjQGeTEkuDxMS0chsT2JmfYDRwEMpu84C/piUcDLFsKYu4QC4+8vA0sQ+EZGWt2IGPHVUSDj9j4ITXupwCQeiTzr7AYuTC9x9OVCe2Jet04ECQsICwMz2BgYAG83scTOrNLN1ZjbFzJIfr9WLIeHtJsYgItI4d3h7Cjw7Fqq3QtH34LjZ0LVf3JHFIuqk0xvYmKa8NLEvWxOARe7+blLZwMTrzcDHwEnAb4GfANe3QgwiIg2rrQpT2rxyMeBw8HVwxJ8hv0vckcUmjsGhnqbMMpTXr2g2iPAo7vKUXXUJ9E13/3Hi96fNrBC4yswmunt5c2Iws3OAcwAGDx6cTZgi0tFVbgw91FbPhvyucPifYK8z4o4qdlHf6ZQCvdKU9yT93Uc6ZxASxPSU8k8Sr8+klD9N6FwwtJEYemWKwd3vdvdidy/u379/lmGKSIe1+X2YdURIOF0HwKi5SjgJUSedxaS0m5jZnkAP0rezpDMBmO/uH6WUvw9UpqlfNzVrbaYYEjK19YiIZG/tfJh1GJQthp4HwYkvQ7/D4o4qZ0SddGYCJyYeedUZD1QA8xo72MyKgMOp32sNd68EngKOS9k1itBRYUlSDAPN7Kik8xYDQ0ga8yMi0mRL74enR4XZogeNgROehx57xR1VTok66dwFbAdKzOz4RFvJRGBKcjdqM1tiZn9Mc/wEoBp4JMP5JwGHJGYbOMHMLgGuAH7r7tsB3P0F4Engz2Y2zsxOAR4g3D1pjI6INJ3Xwr+vgRfOhNpK2OdnMPIfULBr3JHlnEg7Erh7qZmNAm4HZhDaUG4hJJ7UuNLNBzEBmOPu6zKc/2Uz+wZwA/AdYC3wm8R26nluAe4lJN5HgQub/olEpMOrroAXfwDL/wcsDw79PexzftxR5Sxzz6rTmCQUFxf7ggUL4g5DRHJBxZow/mbDS9CpEI76H9j9pLijyglmttDd6y0IpPV0RESaY+MbYdG1rR+GdpuRj0Kvg+KOKucp6YiINNXKJ2D+GVC9GfoeDiP+Dt12izuqNkHr6YiINMU7t8O8k0PCGTweRj2thNMESjoiIg0pKYE99oBu3WDPHvCHn4XeagddA0c+CJ26xR1hm6KkIyKSSUkJnHsurF4F+dtgY3no81p2IRw8KfRWkybRFRMRyWTSJKitgR4e5rXPBzrtCnc1OpZdMlDSERHJZOkHULspJJtaYGsedOkOy5bFHFjbpaQjIpLOxtehTwVU10INIeHUAlVVUFQUc3Btl5KOiEiqdc/DUyNgbDXkF8BWg1oPC7KZwbXXxh1hm6WkIyKS7OPH4OnRULURThsH9/wFBg6CLl1g0CCYOhXGjYs7yjZLg0NFROosvT/Mo+Y1MPRs+PJdkJcPp4+PO7J2Q3c6IiIA7/w+zBLtNXDAFfCVu0PCkRalOx0R6djc4bVr4c3rw/Yh/w/2vzjemNoxJR0R6bhqa2DB+bBkKlg+HPZHGPL9uKNq15R0RKRjqtkeHqct/yvkd4Ujp8Pnvhl3VO2eko6IdDxVW+C5U2H17LC658gZMGBE3FF1CEo6ItKxbFsPc78Gn/wLug6AY5+E3sPijqrDUNIRkY5j60fwzAlQthh67A3HzYLCz8cdVYcSeZdpMzvAzOaYWbmZrTSzSWbWYL9EM5toZp7h58qkevdlqLNfUp2iDHUebs3PLSIx27QYnjoyJJxeX4QTnlfCiUGkdzpm1huYDbwFjAWGApMJye/qBg69B3gipewU4HJgZkr5YuCslLJlac55CfB80vb6Bt5fRNqyDf+CuWNg+wbof2Row+ncO+6oOqSoH6+dB3QDxrl7GfCUme0KTDSzmxNl9bj7CmBFcpmZXQMsdvdXU6pvdfcXs4jlnSzriUhbtno2PHsKVG+F3b8GR/0VOnWPO6oOK+rHa2OAJ1OSy8OERDQy25OYWR9gNPBQy4YnIu3K8kdg7skh4RR9D0b8XQknZlEnnf0Ij78+5e7LgfLEvmydTlhSKV07zAFmVmZm281svpllSmbTzKzGzFaZ2RQz05qzIu3Je1Nh/hlQWwn7/hyO+BPkFcQdVYcX9eO13sDGNOWliX3ZmgAscvd3U8pfAV4itBn1By4mPMI7yt1fTtTZDtwBzALKgGMIbUNDCe1MItKWucObv4XXEs3EB18PB14VliSQ2MXRZdrTlFmG8voVzQYRHsVdXu/E7v+VUvcxQgK6itDxAHdfBVyQVG2uma0B7jSzYWnaiDCzc4BzAAYPHpxNmCISB6+FRRfDO7cCBl++E75wXtxRSZKoH6+VAr3SlPck/R1QOmcQktT0xiq6ewXwODC8kaqPJF7T1nP3u9292N2L+/fvn2WYIhKp2ip44fsh4eQVwFHTlXByUNR3OotJabsxsz2BHqS09TRgAjDf3T9qwvs2dhflKa8i0pZUl4f2m5WPQacecPTfYNDouKOSNKK+05kJnGhmhUll44EKYF5jB5tZEXA4WfZaS3QOGAMsbKTq6YnXxuqJSK6pLA2zDKx8DLr0heOeVsLJYVHf6dwFXAiUmNlNwBBgIjAluRu1mS0B5rn7j1KOnwBUs+NxGEnH9AQeBe4HlgD9gF8AexAeydXVmwgUEgaGlgEjgEuBEnd/rSU+pIhEpGIVPHMibHwdun8Ojp0FPfePOyppQKRJx91LzWwUcDswg9COcwsh8aTGlW5qnAnAHHdfl2bfdmAdYWaDAcA24AVgpLsvSKq3mDAbwdmE8UHLgd8Bv2nWhxKR6EyevOP3c8bC0yfA1qWw674h4fRQR59cZ+5qxmiK4uJiX7BgQeMVRaTl7b574rUKrsqHbWugz5fhmMeha794Y5PPMLOF7l6cWt5om46ZfdXMemRRr6+Zfae5AYqINKikBNatg7WrYe16eHYN7DYKRs1RwmlDsulI8BxwYN2GmeUnRvKndi/+PPCXlgxORAQICefcc4Fq2MXD4Iv7OkPpj6GgsLGjJYdkk3TSDePV0F4Ric6vJ0LVZtiF8L9PrUGnQrj+hpgDk6aKfD0dEZEmWfscLHmT0FcI2G6wzaCgAJYtizMyaQYlHRHJTTXb4JVLYfZI6FcLtfmwNS8kHYCqKigqijVEabpsk066Lm7q9iYireOTRfBEMbz9/8JEnT8/Dbr0hmoPE3q6h/Jrr407UmmibMfp3GdmW1PK/mJm5UnbjfZwExFpUG0VvHkDvHEdeDUU7gNH/Bn6HQZFJTB+PNTUQH4+TJ0K48bFHbE0UTZJ509pyt7MUPflDOUiIg3b9FaYsPOTxDi4fS6EYTfsWHRt3Di48cYd9ZVw2qRGk467nxVFICLSQXktLL4V/n0V1G6H7oPh8Gkw8Lj6dS++OPLwpGXFsZ6OiEiw5QN48SxY+2zYHnIWDL8FOveMNy5pNdnMSHCAmU1IUz7GzBaa2VYz+8DMLmqVCEWk/XGHJXfD4weHhNN1NxjxDzj8XiWcdi6bO51rgL7Aw3UFZnYo8H/AKuAPwD7AZDP72N3/2hqBikg7Ub4SXjobVs0M24O/BcV3aiqbDiKbpHMYkDrs9+dADXC0uy8HMLN7EuVKOiJSnzt8+DAsOD+sgdO5NxTfAXtNCN2fpUPIJukMBN5LKfsaYb2b5Ull/0v6nm4i0tFtWw8LfgrLE3+TDhoDh90D3XePNy6JXDZJZyPQu27DzPYH+gDPptTbCnRvschEpH1YMQNe/nFYhqDTLjB8Cgw9W3c3HVQ2Secl4Hwzm+Hu1cC5hNkI/pFSb3/g4xaOT0TaqqoyWHgRfDAtbA8YAYffB7vsHWdUErNsOxK8AKwys43AUGC6u7+RUu87wPyWDU9E2qTVT4eu0OXLIa8LfOm3sN9FYJrusaPLZnDoG2Y2DPgh0BNYRErbjZn1B/4N3N8KMYpIW1FdDq9eAe/eFrb7HBqmsel5QLxxSc7I6s8Od3/f3X/p7he4+73uXpOyf527X+jujU6Dkxj3M8fMys1spZlNMrP8Ro6ZaGae4efKpHr3ZaizX8r5eprZNDMrNbNNZvaAmfXN5lqISAbrX4SZh4SEY53gi7+GE15QwpHPiHRGAjPrDcwG3gLGEh7VTSYkv6sbOPQe4ImUslOAy4GZKeWLgdSpe5albE8H9gXOBmqBm4C/A0c3+iFEBCZP3vH7RT+DN34Nb90YprTpeUC4u+lzaHzxSc5qNOmY2QdNOaG7D2lg93lAN2Ccu5cBT5nZrsBEM7s5UZbunCuAFSlxXQMsdvdXU6pvdfcXMwVgZkcAJwIj3f3ZRNnHwEtmdry7z274E4p0cCUlcMUVYbbnrnnwyW1w4IeAwf6XwsGTIL9r3FFKjsrmTqcI2Aw8CqzcyfcbAzyZklweJtxpjARmZHMSM+sDjAaub2YMa+oSDoC7v2xmSxP7lHREMikpgXPPDQmni0N+Ddz2IZw/AH7xvzDgqLgjlByXTdL5FXAGMB54DngIeMTdP2nG++0HPJ1c4O7LE+vy7EeWSQc4HSggaWqeJAeYWRnQBfgX8Et3n5cSw+I0x72d2CcimUyaBF4DuzjUtcR6F3h8N7hBCUca12hHAne/zt2/CAwDngcuI3SfftzMzjSzwia8X2/CYNNUpSQNQM3CBGCRu7+bUv4KcDHwDeC7hH8WT5nZV1ohBpGOZ9lSqC0L/7IcKDfo1hs+XN7YkSJA9stV4+5vuPvV7v55QoP728BvgTVm9qsmvGe6Za4tQ3n9imaDCI/iHkoT43+5+x/cfZ67PwIcRxiwetXOxGBm55jZAjNbsG7dumzCFGl/tq2DPpVQVRO632zJg2qDqiooKoo7OmkjmjtSaxHhMdlzQGfCbATZKAV6pSnvSfq7j3TOICSI6Y1VdPcK4HFgeBYx9MoUg7vf7e7F7l7cv3//LMMUaUe2rYU5x8E3t0FePmwxqPUwiacZXHtt3BFKG5F10rHgODP7b2A1YSBoJfB1wqOsbCwmpd3EzPYEepC+nSWdCcB8d/8oy/rw2TuYejEkZGrrEenY6hLOpjfg+P1g6t0h8QDk58PUqVo6WrKWTZfpI4BvExrvC4HHCONbHnf3yia+30zgUjMrdPfNibLxQAUwL/Nhn8ZSBBwO/DSbNzOzboQeaQtTYrjGzI5y9/mJesXAEOqP+RHp2CrWwNPHwaa3wvib456GbrvBx6U76ijhSBOYe8NNKWZWS+gyPYOwcNvWhuq7++MNnKs3YWDoG4Ru0kOAKcCt7n51Ur0lhKUTfpRy/BXAdcDu7r4uZV9PQrfu+4ElQD/gF8AhwJHuviCp7hOEhecuYcfg0LXu3ujg0OLiYl+wYEFj1UTavorV4Q6n7G3oeSAcNyckHJEsmNlCdy9OLc92RoJCwoSe3ya0p2Ti7OhIWX+ne6mZjQJuJySxjcAtwMQ0caU7zwRgTmrCSdgOrCPMbDAA2EaYqHRkcsJJOs8twL2ER4yPAhc28LlEOpaKVYmEsxh6HgSj5kDXAXFHJe1ANnc6ezXlhO7+4U5FlON0pyPtXvlKmHMsbH4Xeh0Mx82GrupAI03T7DudbJOImR1LGMMzpunhiUhOKP84kXDeg15fSiScfnFHJe1IVo/XzKwXcBKwJ/AB8A93r0rs+xZh4s3hQOpgTRFpK8pXwOxjYcsS6D0sJJwumnxdWlY2vde+CMwCklsQF5nZacCDhN5kbxG6TTc6dkZEctDWj8Idzpb3ofchiYTTJ+6opB3KZpzOb4Ey4AigO2Eg6CeEec0OAr7v7l9094fcvbbVIhWR1rF1Ocw5JiScPocq4UiryubxWjHwc3d/KbH9jpn9BHgPOMfdtVqoSFu19cPwSG3rUuhTDMfNgs6aglBaTzZJZzfqL4JWt/3vlgxGRCK0ZVl4pLZ1GfT5ciLh9Io5KGnvsh2nk6lfdXVLBSIiEdqyFGYfA+XLoe9hcOyT0Lln3FFJB5Bt0nnSzNIlmDmp5e6uEWQiuWzLB4mE8xH0OwKOfQIKdo07Kukgskk6v271KEQkGpvfD50GyldAv6/CsTOVcCRS2QwOVdIRaQ82Lwl3OBUfQ/8j4ZiZUNCUNRhFdl5z19MRkbak7D2YPTKRcI5WwpHYKOmItHdl78CckVCxEgaMgGMeV8KR2GTbkUBE2qJNi0O36G2rYcAxcMyj0KlH3FFJB6Y7HZH2atPbodPAttWw27FKOJITlHRE2qNNbyUSzhrYbRSMVMKR3KDHayLtzcY3wgJs29fBwNEw4v+gU7e4oxIBdKcj0r5sfD0p4ZyghCM5R3c6Im3Z5Mk7fv/haHh6FGxfD4NOghF/g/yu8cUmkoaSjkhbVpd0BlVB0Q2wfQMMGgMjSpRwJCdF/njNzA4wszlmVm5mK81skpnlN3LMRDPzDD9XZjjmlMT+BSnlRRnO83BLfk6RVldSAuvWwdrVsHY9PLcBdj9ZdziS0yK90zGz3sBswkqjY4GhwGRC8ru6gUPvAZ5IKTuFsEz2zDTv0xWYAqxp4JyXAM8nba9vOHqRHFJSAueeC14DPRw2Avd1hqPOhPwucUcnklHUj9fOA7oB49y9DHjKzHYFJprZzYmyetx9BbAiuczMrgEWu/uraQ65FPgYeJ+wumk677j7i837GCIxmzQJqrfDLg4G1AKdCuH6G+D08XFHJ5JR1I/XxgBPpiSXhwmJaGS2JzGzPsBo4KE0+wYDlwE/37lQRXLU5iWw5E3wzSHhVAPleVBQAMuWxRycSMOiTjr7AYuTC9x9OVCe2Jet04ECQsJKNRn4H3df1Mg5pplZjZmtMrMpZqZ+pZLbqivgtV/BYwdB3+pwd7PNQsIBqKqCoqI4IxRpVNSP13oTnj6nKk3sy9YEYJG7v5tcaGbHAicC+zRw7HbgDmAWUAYcQ2gbGkpoZxLJPR8/Dgt/FhZgA/jJsTD537CllE8X9jWDa6+NLUSRbMTRZTrd0teWobx+RbNBhEdxl6eUdwJ+D1zv7qszvrn7KuCCpKK5ZrYGuNPMhqVrIzKzc4BzAAYPHpxNmCItY+tyWHgRrPhb2O55EHz5DzDgKNinBMaPh5oayM+HqVNh3LhYwxVpTNSP10qBXmnKe5L+DiidMwhJanpK+Y8T5/6TmfUys15AZyA/sV3QwDkfSbwOT7fT3e9292J3L+7fv3+WYYrshJpKePNGeHT/kHA67QKHTIYxi0LCgZBg+veHgQPDqxKOtAFR3+ksJqXtxsz2BHqQ0tbTgAnAfHf/KKV8X+BzQLq7nFLgTOD+DOf0lFeR+Kx5Bv71UyhL/JMYfAYMnwLd96hf9+KLo41NZCdFnXRmApeaWaG7b06UjQcqgHmNHWxmRcDhwE/T7L4d+HtK2RXA3sC5wNsNnPr0xOvCxmIQaTUVq2DRJfDhg2G78AtQfAcMGp35GCUdaWOiTjp3ARcCJWZ2EzAEmAhMSe5GbWZLgHnu/qOU4ycQOog+klKOuy8BliSXmdkPgH7uPjepbCJQSBgYWgaMIIzrKXH313bq04k0R201vHcnvHYNVJWF2QQO/CXsf6kGekq7E2nScfdSMxtFuCuZQWjHuYWQeFLjSjc1zgRgjruv24kwFhNmIzibMD5oOfA74Dc7cU6R5ln3Aiz4KZS+GrZ3/zoU/x522TvWsERai7mrGaMpiouLfcGCBY1XFGnItvXw7yvg/T+G7R57waG/h899M964RFqImS109+LUcs0yLRIlrw2J5tUroPITyCsIj9EO/CV06h53dCKtTklHJCqfvAL/+glseCls7zYKim+Hnk2ZjEOkbVPSEWltlZtCJ4H37gh3Ot0GwfBbQldos7ijE4mUko5Ia3GHZQ/CKxfDtjVg+bDvL+DgiVCwa9zRicRCSac9S17KWOM5Wk+667zpLfjX+bB2btju99UwfU3vgyMPTySXqPdaE7Wp3mu7777j95Uro3vfjpbskq/zh+/CG9fB4ing1dClHwy7GYZ8HyzyhXpFYqPeaxKdjpR06paMrqmBHgZXF8GXNgAGnz8XvvRb6NIn7ihFcob+9BJprk+XjK6G7g7Uwh82wFt7wwkvwlfuUsIRSaGkI9IcW5bClefD9k/CpEoFhDk08nvAo4XQ7ysxByiSm5R0pGXVPW5avTq8lpTEHVHLqSyFJXfDU0fDP4bAitWQVxvmJq802JIHXQth2YdxRyqSs9SmIy2n7nFTTU3YrqkJ29B213qpqYRVM2HpX+DjGVBbGcrzu8HuXWBDDWzauqO+lowWaZCSjjSdO2zfAFuXhsdMWz4Iv1/2IGwvh54e/vqvAarL4VdXwqmntJ3eW+5h1oClf4Hl08NnBcDCLAJ7nwl7joOCpxJtOlt2HKslo0UapKQj6VVXwNZlIaEkJ5a67erN9Y9ZRWjbgLC2ayfAy+GDd+F/+8OAkbDbMTDgGOh1UO4loS0fwNL7Ydn9sPm9HeU9DwqJpug70P1zO8rr7t60ZLRI1pR02qvkrrz5+WE7+T/D2hqo+Dh9QtnyAWxLtwBrkk6FsMuQxM/e0GNvKJoCq0thY1loLcx36FIAAzxMbrnib+EHoHOflk9CzemqXVkKy/8a7mrWzd9R3nVgSDJ7nwm9vpR5uppx4+DGGz+7LSIZaXBoE7WJwaF1bSufbIA8D3cdu3SHy4+CL3tILOUfQm1V5nNYJ+hRFBJKcmKp+71zn/r/Ede974YNO8r69g1//Z8wDNbMDT9rn4HyFZ89tiWSULaDYTO203SHPU+FojNh4CjI099kIs2VaXCokk4TtYmkM2wYLP8AajeHhANhvdU+wA1J9boO3JFEdhmSSCqJ37vtAXnp1tFrREnJZx83TZ9e/69/93Bn1dJJqKGk02A7zXFJ7TSFTf7IIlKfkk4LaRNJp2ch5CUat2sTP126Q6XD239N3LUUtd76LU2dfqelklC6921qO42ItAhNg9NRlL0HvStgA6EHWUXiP+Yu3WDI52CPk+OMLj2zHe1DQ3+YOQk11Cb0zHs72rAK8uGOc2GfN2Hd8zveJ9t2GhFpNZEnHTM7ALgNOALYCNwD/Nrdaxo4ZiLwqwy7r3L3G1ILzewU4G9AvWxrZj2BW4FTCE3ejwIXuvsG2rKKNfDMiXBKDdxXAGXVhMxD2+rK29Qk9DJwr0GBQxfAquHKu+GHwBFqpxHJJZH+CzSz3sBs4C1gLDAUmEz4j//qBg69B3gipewU4HJgZpr36QpMAdZkON90YF/gbMLDp5uAvwNHZ/VBclHVZpj7tfCf80nFcOQv4Dvfj6crb0tP8tlYEvrlw+BVIeHUqe0EswbB795UO41IDom0TcfMrgQuA/Zy97JE2WXARGBgXVmW53oMGOLu+6fZdw1wAvA+cFDynY6ZHQH8Exjp7s8myr4CvASMdvfZDb1vTrbp1FbB3K/D6lmwy1A44Z/QdUB8SxtErVcv6NoZNq4PHSeqDPoPgIoK2Lgx5uBEOqZcadMZAzyZklweJtxpjARmZHMSM+sDjAauT7NvMCGxjQQuzBDDmrqEA+DuL5vZ0sS+BpPOzrj11ltZt25dC5/VOW2PRxne63W2VHdn6qIT+OSl/wLg8s07BnDe9MtftvD75o7zu3ShZ9lmulXWlTgVmzaxqbCQO9rx5xZpbf379+eiiy5q0XNGnXT2A55OLnD35WZWntiXVdIBTieMfX84zb7JwP+4+yJL31C8H7A4TfnbiX2tZt26dey1114tes4vd/4bh3R+nSrvwqzKn1O4exF1D5Py83d0eW7p980lr596KiMeeCA8dkvIy8/n9VNPbdefW6S1ffhhy09eG/U8JL0JnQdSlSb2ZWsCsMjd300uNLNjgROBhv68bakYYndgwTMc0vkJaj2Pp7adw/raorhDisWHw4fz7He/i+flYYDn5fHsd7/Lh8OHxx2aiKSIoytPukYky1Bev6LZIMKjs8tTyjsBvweud/dG5nBpWgxmdg5wDsDgwYOzCbPV7Z2/kK92ng7As9vPZEXNQTFHFK8Phw9nW2HhZ7ZFJPdEfadTCvRKU96T9Hcf6ZxBSBDTU8p/nDj3n8ysl5n1AjoD+YntuqkoM8XQK1MM7n63uxe7e3H//v2zDLP1DMx7l2O73ouZ8/L2sbxb/dW4QxIRyUrUdzqLSWk3MbM9gR6kb2dJZwIw390/SinfF/gckO4upxQ4E7g/8T7pukbvR+g2nVO++NRTn/7++ujR9M77mBO73Uknq+bNypG8WjUmxuhyy2ujR8cdgog0IuqkMxO41MwK3b2ua9V4oAKY19jBZlYEHA78NM3u26mfNK4A9gbOJXQUqIvhGjM7yt3nJ85bDAwhzZifuB2clHQ+OOFQxnT9PV2sgqXVh/DPygnsmFxNXlfSEcl5USeduwjdmEvM7CbCf/QTgSnJ3ajNbAkwz91/lHL8BMLUlY+kntjdlwBLksvM7AdAP3efm1TvBTN7EvizmV3CjsGh8xsboxOrbs6YrrexS95GVtV8nqe3/RBv5Omo/vIXkVwTadJx91IzG0W4K5lBaEO5hZB4UuNKN8XxBGCOu+/sYJcJife9l6RpcHbynK2nk9P1vK30yC+jtHYQsyp+Sg2dGz1Mf/mLSK6JvPeau78FHNdInaIM5cOa+F4/yFC+ETgr8ZOz9lq0iK6bN5PXqRa7FypO6c7jB13IdnrEHZqISLPk2HrBUmevRYsY8cAD5BXUYt3AN0DevdBv4bK4QxMRaTYlnRxV/NijFLAtJBwDKo1aNw597LG4QxMRaTYlnRzUwz6h94aVdMqvDqNVyw2vhtr8fAo3tO3VF0SkY9PiIjlmUP47HN/lv8nr79R+ArbNqFtpKK+mhs19+8YboIjITtCdTlQmT+bIl176zGDPz3IOLpjFyV1vpVveZtaP3YNKeuDVHiaydAczFp6cgyt/iohkSXc6UZk8maM2byY/P79eV+ZObGNklz8ztGAhAK9UnsSCg8Yy+Luvcvx//zd5tbXUahJLEWkHlHRi1tPWMLrrH+iTv4pK78rcbT9gWc0hgCaxFJH2R0knQnl5edTU1Hy6RsX+he8ydo8ZdM3fztrtfXlg+Wmsr+wD7FjDoqam5tPfW2NtCxGRTFpjguNIl6tuD5q9XHXy0tErPoLXfwVv/iZs73kaHD4NCgobPq49LzktIu1KrixXLd1rYd7JsOpJsDz40o2w/yWQfpVTuPjiaOMTEWlFSjpRKCmBdevAa2Cgw/89CUf3gyMfhoGjGj5WSUdE2hF1mW5tJSVw7rlg1bCLwyZgWieo/E3jCUdEpJ1R0mltkyZBbXVYpq5uQeyCXeGmO2MOTEQkeko6rW3ZMujcFbYbVBhU5EFB51AuItLBKOm0tqIiqKoKSacq0VmgqiqUi4h0MEo6re3aa0PPNP/sdDZce23ckYmIRE6911rbuHHhdfx4qKmB/HyYOnVHuYhIB6KkE4Vx4yB5ZK8Sjoh0UJE/XjOzA8xsjpmVm9lKM5tkZvmNHDPRzDzDz5VJ9X5tZq+bWZmZbTazBWY2PuVcRRnO83BrfWYREQkivdMxs97AbOAtYCwwFJhMSH5XN3DoPcATKWWnAJcDM5PKdgXuS5y/BjgdeNjMatz9kZTjLwGeT9pe34SPIiIizRD147XzgG7AOHcvA54ys12BiWZ2c6KsHndfAaxILjOza4DF7v5qUr1fpBw6y8wOBP4DSE0677j7izv1aUREpEmiTjpjgCdTksvDwE3ASGBGNicxsz7AaOD6LKpvADo3Mc6Wp+lsREQiTzr7AU8nF7j7cjMrT+zLKukQHpsVEBJWPWbWCdgFOBk4AZiQptq0RPJaCzwE/NLdK7J8/6ZT0hERiTzp9AY2pikvTezL1gRgkbu/m7rDzA4HXkhsVgMXuPvfk6psB+4AZgFlwDGEtqGhhHYmERFpJXF0mU63gI9lKK9f0WwQ4VHc5RmqvA58GehFuNO53czK3P0hAHdfBVyQVH+uma0B7jSzYcltREnveQ5wDsDgwYOzCVNERNKIust0KSEZpOpJ+jugdM4gJKnp6Xa6+1Z3X+DusxMdC/5CaDNqSF0ng7TrQbv73e5e7O7FrbGSnohIRxF10llMaLv5lJntSZiDeXGW55gAzHf3j7KsvwjY08wKGqjjKa8iItIKok46M4ETzSx5XebxQAUwr7GDzawIOJzQ8J+tI4EV7l7VQJ3TE68Lm3BeERFpoqjbdO4CLgRKzOwmYAgwEZiS3I3azJYA89z9RynHTyB0Dkgdc4OZ7QVMAx4EPiD0Xjs1ccxPkupNBAoJA0PLgBHApUCJu7/WEh9SRETSizTpuHupmY0Cbid0j94I3EJIPKlxpZsaZwIwx93Xpdm3EVhJmNlgYGL7LeBkd388qd5iwmwEZxMGqi4Hfgf8phkfSUREmsDc1YzRFMXFxb5gwYK4wxARyWlmttDdi+uVK+k0jZmtAz5s4mH90NxuTaHr1TS6XtnTtWqanblee7l7ve6+SjoRMLMF6TK+pKfr1TS6XtnTtWqa1rheWjlUREQio6QjIiKRUdKJxt1xB9DG6Ho1ja5X9nStmqbFr5fadEREJDK60xERkcgo6ewkMzvAzOaYWbmZrTSzSWaWbmBr6nE9zWyamZWa2SYze8DM+kYRc1yac63MrMjMPM1P2rWU2gsz+7yZTTWzf5tZjZnNzfK4Dve9guZdrw783fqWmf3DzD42sy1mttDMvp3FcS3y3YpjaYN2w8x6A7MJMx+MJazJM5mQzK9u5PDpwL6EmRFqCTNh/x04upXCjdVOXisIs0g8n7Td3sdaHAh8DXiRpq1826G+V0mae72g4323/hNYCvyC8Fm/BjxoZv3c/bYGjmuZ75a766eZP8CVhOUadk0quwwoTy5Lc9wRhBmtRySVfSVRdnzcnyvHrlVR4rp8Pe7PEPH1ykv6/RFgbhbHdLjv1U5er4763eqXpuxBYGkDx7TYd0uP13bOGOBJT5qslLCEdjfCQnMNHbfG3Z+tK3D3lwl/fYxpjUBzQHOvVYfk7rXNOKwjfq+AZl+vDsnd093JvQIMaOCwFvtuKensnP1IWQfI3ZcT/nrfL+0RGY5LeLuR49qy5l6rOtMSz+pXmdkUM+vWGkG2cR3xe9US9N2CrxIefWfSYt8ttensnN6kX/G0NLGvOccN2emoclNzr9V24A5gFmEpimMIS5UPJbQNyQ4d8Xu1M/TdAhIz/48FfthAtRb7binp7Lx0A50sQ3lLHNeWNfkzu/sq4IKkorlmtga408yGufurLRtim9cRv1fNou/WpwtjPgj8n7vf10j1Fvlu6fHazikFeqUp70n6vwoaO65XI8e1Zc29VunULeI3fCfiaY864veqpXWY75aZ9SGs5rwc+F4j1Vvsu6Wks3MWk/I808z2BHqQ/vlnxuMSMj03bQ+ae63S8ZRXCTri96qldYjvlpl1Bx4ldC8/2d23NnJIi323lHR2zkzgRDMrTCobD1QA8xo5bqCZHVVXYGbFhGejM1sj0BzQ3GuVzumJ14UtEVg70hG/Vy2t3X+3zKwT8FfgC8AYd1+bxWEt9t3S3Gs7ITHg8S3gDcJAqSHAFOBWd786qd4SYJ67/yip7AlgH8LAtLqBVmvdvV0O4mvutTKziUAhYfBeGTACuBR43N1Pi/IzRCnxl+jXEpsXA7sCv0psP+7u5fpe7dCc69WBv1t3Az8Gfg68nLL7FXff3qrfrbgHKrX1H+AA4GnCX+yrgOuA/JQ6y4D7Usp6AdMIz0PLCI159QZttaef5lwrYAKwANgEVAJLgElAl7g/TytfqyLCI550P0X6Xu389erA361lcX63dKcjIiKRUZuOiIhERklHREQio6QjIiKRUdIREZHIKOmIiEhklHRERCQySjoiIhIZJR2RNszMfm5mryfWui8zsxfM7IC44xLJREsbiLRRZvZj4DeEZb9fJ0zeOByoiTMukYZoRgKRNsrMHiSsU3+0h7VhRHKeHq+JtF13AYOBj8zseTM738zy4w5KpCFKOiJtkJl1Af6TMFN3MTADmAw8FGdcIo1Rm45I2zQFwN2vSGy/mlgn5Toz293dV8YXmkhmSjoibUxixdWfACNTdi1KvPYBlHQkJ+nxmkjbMxaoAl5IKd+NsLjWx5FHJJIlJR2RtmcvoNTdq1PKTwL+6e6lMcQkkhUlHZG2pwzon1gCHAAzOxQ4DZgaW1QiWdA4HZE2xswOBl4BZgO3AnsD1wOz3H1CjKGJNEpJR6QNMrNvA5MIj9o+Au4FbkrzyE0kpyjpiIhIZNSmIyIikVHSERGRyCjpiIhIZJR0REQkMko6IiISGSUdERGJjJKOiIhERklHREQio6QjIiKR+f+/mcI+L1mPDQAAAABJRU5ErkJggg==\n", "text/plain": [ "<Figure size 432x288 with 1 Axes>" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# plot and save figure\n", "plt.figure()\n", "plt.plot(deming_factor_list, np.mean(rmse_list,axis=1), linewidth=2, color='orange')\n", "plt.errorbar(deming_factor_list, np.mean(rmse_list,axis=1), np.std(rmse_list,axis=1)/np.sqrt(rmse_list.shape[1]), color='r', linewidth=3,alpha=0.9,ecolor='r',fmt='o',linestyle='None')\n", "plt.fill_between(deming_factor_list, np.mean(ber_rmse_list)+ np.std(ber_rmse_list)/np.sqrt(len(ber_rmse_list)), np.mean(ber_rmse_list)- np.std(ber_rmse_list)/np.sqrt(len(ber_rmse_list)), color='k', alpha=0.4)\n", "plt.xlabel(r'$\\delta$')\n", "plt.ylabel('RMSE')\n", "plt.tight_layout()\n", "plt.savefig(os.path.join('saved_images','wine_deming_rmse.pdf'))" ] }, { "cell_type": "code", "execution_count": 88, "id": "d502a44f", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "EiV: 1.716 (0.009)\n", "non-EiV: 2.545 (0.026)\n", "EiV: 0.960 (0.009)\n", "non-EiV: 1.085 (0.002)\n" ] } ], "source": [ "# Error vs. Uncertainty\n", "print('EiV: %.3f (%.3f)' % (np.mean(unc_q), np.std(unc_q)/np.sqrt(len(unc_q))))\n", "print('non-EiV: %.3f (%.3f)' % (np.mean(ber_unc_q), np.std(ber_unc_q)/np.sqrt(len(ber_unc_q))))\n", "print('EiV: %.3f (%.3f)' % (np.mean(full_q), np.std(unc_q)/np.sqrt(len(full_q))))\n", "print('non-EiV: %.3f (%.3f)' % (np.mean(ber_full_q), np.std(ber_full_q)/np.sqrt(len(ber_full_q))))" ] }, { "cell_type": "code", "execution_count": 89, "id": "df5e5040", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Uncertainty only \n", "--------\n", "EiV: 0.842 (0.002)\n", "non-EiV: 0.674 (0.003)\n", "Uncertainty + Noise \n", "--------\n", "EiV: 0.958 (0.002)\n", "non-EiV: 0.931 (0.001)\n" ] } ], "source": [ "# Coverage\n", "print('Uncertainty only \\n--------')\n", "print('EiV: %.3f (%.3f)' % (np.mean(unc_cov), np.std(unc_cov)/np.sqrt(len(unc_cov))))\n", "print('non-EiV: %.3f (%.3f)' % (np.mean(ber_unc_cov), np.std(ber_unc_cov)/np.sqrt(len(ber_unc_cov))))\n", "print('Uncertainty + Noise \\n--------')\n", "print('EiV: %.3f (%.3f)' % (np.mean(full_cov), np.std(unc_cov)/np.sqrt(len(full_cov))))\n", "print('non-EiV: %.3f (%.3f)' % (np.mean(ber_full_cov), np.std(ber_full_cov)/np.sqrt(len(ber_full_cov))))" ] }, { "cell_type": "code", "execution_count": null, "id": "607b92ac", "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.8.5" } }, "nbformat": 4, "nbformat_minor": 5 }