sgd.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Linear regression and stochastic gradient descent"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### LR:\n",
    "$\\Large \\frac{1}{\\ell}\\sum_{i=1}^\\ell{{((w_0 + w_1x_{i1} + w_2x_{i2} +  w_3x_{i3}) - y_i)}^2} \\rightarrow \\min_{w_0, w_1, w_2, w_3},$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### GD:\n",
    "$\\Large w_0 \\leftarrow w_0 - \\frac{2\\eta}{\\ell} \\sum_{i=1}^\\ell{{((w_0 + w_1x_{i1} + w_2x_{i2} +  w_3x_{i3}) - y_i)}}$\n",
    "$\\Large w_j \\leftarrow w_j - \\frac{2\\eta}{\\ell} \\sum_{i=1}^\\ell{{x_{ij}((w_0 + w_1x_{i1} + w_2x_{i2} +  w_3x_{i3}) - y_i)}},\\ j \\in \\{1,2,3\\}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\eta$ - GD step."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### SGD:\n",
    "$\\Large w_0 \\leftarrow w_0 - \\frac{2\\eta}{\\ell} {((w_0 + w_1x_{k1} + w_2x_{k2} +  w_3x_{k3}) - y_k)}$\n",
    "$\\Large w_j \\leftarrow w_j - \\frac{2\\eta}{\\ell} {x_{kj}((w_0 + w_1x_{k1} + w_2x_{k2} +  w_3x_{k3}) - y_k)},\\ j \\in \\{1,2,3\\},$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$k$ - random index, $k \\in \\{1, \\ldots, \\ell\\}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Normal equation\n",
    "$\\Large w = {(X^TX)}^{-1}X^Ty$ ([*Generalized inverse*](https://en.wikipedia.org/wiki/Generalized_inverse))\n",
    "\n",
    "\n",
    "$\\Large X^TXw = X^Ty$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Data: [faculty.marshall.usc.edu](http://www-bcf.usc.edu/~gareth/ISL/data.html).**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>TV</th>\n",
       "      <th>Radio</th>\n",
       "      <th>Newspaper</th>\n",
       "      <th>Sales</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>230.1</td>\n",
       "      <td>37.8</td>\n",
       "      <td>69.2</td>\n",
       "      <td>22.1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>44.5</td>\n",
       "      <td>39.3</td>\n",
       "      <td>45.1</td>\n",
       "      <td>10.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>17.2</td>\n",
       "      <td>45.9</td>\n",
       "      <td>69.3</td>\n",
       "      <td>9.3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>151.5</td>\n",
       "      <td>41.3</td>\n",
       "      <td>58.5</td>\n",
       "      <td>18.5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>180.8</td>\n",
       "      <td>10.8</td>\n",
       "      <td>58.4</td>\n",
       "      <td>12.9</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      TV  Radio  Newspaper  Sales\n",
       "1  230.1   37.8       69.2   22.1\n",
       "2   44.5   39.3       45.1   10.4\n",
       "3   17.2   45.9       69.3    9.3\n",
       "4  151.5   41.3       58.5   18.5\n",
       "5  180.8   10.8       58.4   12.9"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "adver_data = pd.read_csv('advertising.csv')\n",
    "\n",
    "adver_data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>TV</th>\n",
       "      <th>Radio</th>\n",
       "      <th>Newspaper</th>\n",
       "      <th>Sales</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>count</th>\n",
       "      <td>200.000000</td>\n",
       "      <td>200.000000</td>\n",
       "      <td>200.000000</td>\n",
       "      <td>200.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mean</th>\n",
       "      <td>147.042500</td>\n",
       "      <td>23.264000</td>\n",
       "      <td>30.554000</td>\n",
       "      <td>14.022500</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>std</th>\n",
       "      <td>85.854236</td>\n",
       "      <td>14.846809</td>\n",
       "      <td>21.778621</td>\n",
       "      <td>5.217457</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>min</th>\n",
       "      <td>0.700000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.300000</td>\n",
       "      <td>1.600000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25%</th>\n",
       "      <td>74.375000</td>\n",
       "      <td>9.975000</td>\n",
       "      <td>12.750000</td>\n",
       "      <td>10.375000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>50%</th>\n",
       "      <td>149.750000</td>\n",
       "      <td>22.900000</td>\n",
       "      <td>25.750000</td>\n",
       "      <td>12.900000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>75%</th>\n",
       "      <td>218.825000</td>\n",
       "      <td>36.525000</td>\n",
       "      <td>45.100000</td>\n",
       "      <td>17.400000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>max</th>\n",
       "      <td>296.400000</td>\n",
       "      <td>49.600000</td>\n",
       "      <td>114.000000</td>\n",
       "      <td>27.000000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "               TV       Radio   Newspaper       Sales\n",
       "count  200.000000  200.000000  200.000000  200.000000\n",
       "mean   147.042500   23.264000   30.554000   14.022500\n",
       "std     85.854236   14.846809   21.778621    5.217457\n",
       "min      0.700000    0.000000    0.300000    1.600000\n",
       "25%     74.375000    9.975000   12.750000   10.375000\n",
       "50%    149.750000   22.900000   25.750000   12.900000\n",
       "75%    218.825000   36.525000   45.100000   17.400000\n",
       "max    296.400000   49.600000  114.000000   27.000000"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "adver_data.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[230.1  37.8  69.2]\n",
      " [ 44.5  39.3  45.1]\n",
      " [ 17.2  45.9  69.3]\n",
      " [151.5  41.3  58.5]\n",
      " [180.8  10.8  58.4]]\n"
     ]
    }
   ],
   "source": [
    "X = adver_data[['TV', 'Radio', 'Newspaper']].values\n",
    "y = adver_data[['Sales']].values\n",
    "\n",
    "print(X[:5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[147.0425  23.264   30.554 ] [85.63933176 14.80964564 21.72410606]\n"
     ]
    }
   ],
   "source": [
    "# scaling\n",
    "\n",
    "means, stds = X.mean(axis=0), X.std(axis=0)\n",
    "print(means, stds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.96985227  0.98152247  1.77894547]\n",
      " [-1.19737623  1.08280781  0.66957876]\n",
      " [-1.51615499  1.52846331  1.78354865]]\n"
     ]
    }
   ],
   "source": [
    "X = (X - means) / stds\n",
    "print(X[:3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1.]\n",
      " [1.]\n",
      " [1.]\n",
      " [1.]\n",
      " [1.]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "# for w0\n",
    "ones = np.ones((X.shape[0], 1))\n",
    "print(ones[:5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.96985227  0.98152247  1.77894547  1.        ]\n",
      " [-1.19737623  1.08280781  0.66957876  1.        ]\n",
      " [-1.51615499  1.52846331  1.78354865  1.        ]\n",
      " [ 0.05204968  1.21785493  1.28640506  1.        ]\n",
      " [ 0.3941822  -0.84161366  1.28180188  1.        ]]\n"
     ]
    }
   ],
   "source": [
    "X = np.hstack((X, np.ones((X.shape[0], 1))))\n",
    "print(X[:5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [],
   "source": [
    "# MSE function\n",
    "\n",
    "def mserror(y, y_pred):\n",
    "    err = np.mean((y - y_pred)**2)\n",
    "    return err"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "MSE of forecasting Sales if the median is always predicted from the original sample:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "28.346\n"
     ]
    }
   ],
   "source": [
    "print(round(mserror(y, np.median(y)), 3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 3.91925365]\n",
      " [ 2.79206274]\n",
      " [-0.02253861]\n",
      " [14.0225    ]]\n"
     ]
    }
   ],
   "source": [
    "# Normal equation function\n",
    "\n",
    "def normal_equation(X, y):\n",
    "    return np.dot(np.linalg.pinv(X), y)\n",
    "\n",
    "\n",
    "norm_eq_weights = normal_equation(X, y)\n",
    "print(norm_eq_weights)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "14.022\n"
     ]
    }
   ],
   "source": [
    "y_prediction = np.mean(np.dot(X, norm_eq_weights))\n",
    "print(np.round(y_prediction, 3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.784\n"
     ]
    }
   ],
   "source": [
    "def linear_prediction(X, w):\n",
    "    return np.dot(X, w)\n",
    "\n",
    "\n",
    "answer3 = mserror(y, linear_prediction(X, norm_eq_weights))\n",
    "print(round(answer3, 3)) \n",
    "\n",
    "# much better than median prediction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 3.91925365]\n",
      " [ 2.79206274]\n",
      " [ 0.02253861]\n",
      " [14.0225    ]]\n"
     ]
    }
   ],
   "source": [
    "# SGD step function\n",
    "# train_ind - index of the training sample object\n",
    "\n",
    "def stochastic_gradient_step(X, y, w, train_ind, eta=0.01):\n",
    "    \n",
    "    grad = np.ndarray(shape=(4,1))\n",
    "    \n",
    "    # computation the gradient for each parameter\n",
    "    for i in range(4):\n",
    "        grad[i,0] = (2.0/200) * X[train_ind,i] * (linear_prediction(X[train_ind,:4], w) - y[train_ind])\n",
    "    \n",
    "    return w - eta * grad\n",
    "\n",
    "grad = np.ndarray(shape=(4,1))\n",
    "print(grad)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {},
   "outputs": [],
   "source": [
    "def stochastic_gradient_descent(X, y, w_init, eta=1e-2, max_iter=1e4, min_weight_dist=1e-8, seed=42, verbose=False):\n",
    "    \n",
    "    weight_dist = np.inf\n",
    "    w = w_init\n",
    "    errors = []\n",
    "    iter_num = 0\n",
    "    np.random.seed(seed)\n",
    "        \n",
    "    while weight_dist > min_weight_dist and iter_num < max_iter:\n",
    "        random_ind = np.random.randint(X.shape[0])\n",
    "        \n",
    "        weight_dist = np.linalg.norm(w - stochastic_gradient_step(X, y, w, random_ind, eta=0.01))\n",
    "        w = stochastic_gradient_step(X, y, w, random_ind, eta=0.01)\n",
    "        errors.append(mserror(y, linear_prediction(X, w)))\n",
    "        iter_num += 1\n",
    "        \n",
    "    return w, errors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$10^5$ iterations of SGD:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 10.4 s, sys: 15.6 ms, total: 10.5 s\n",
      "Wall time: 10.5 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "stoch_grad_desc_weights, stoch_errors_by_iter = stochastic_gradient_descent(X, y, np.zeros(shape=(4,1)), eta=1e-2, max_iter=1e5, min_weight_dist=1e-8, seed=42, verbose=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Error on the first 50 iterations of SGD:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'MSE')"
      ]
     },
     "execution_count": 101,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAEJCAYAAABc/7oDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXhV1dn38e9NEhLmMUCAMM+gggQQAQecqFLnqSqgaBGHvmpbizjU1j7to7V1qCMooFa09hGwap1wRhAkzEOQQeYwQ8IQCATu94+zaY8xg4RzcjL8PteVi33WXnude7XInb3X2muZuyMiIhIJVWIdgIiIVBxKKiIiEjFKKiIiEjFKKiIiEjFKKiIiEjFKKiIiEjFRSypmlmpmn5lZhpktMbM7gvJHzWyZmS00sylmVjco721m84OfBWZ2SSHtTgurl2lmbwXlZ5hZdti530arbyIiUjCL1nsqZpYCpLj7XDOrBcwBLgaaA5+6e56ZPQLg7qPMrDpwMChPARYATd09r4jvmAT8y91fMbMzgF+7++CodEhERIoVH62G3X0TsCk43mNmGUAzd/8orNpM4PKgTk5YeRJQZLYLEtVA4IaSxtiwYUNv1apVSS8XEamU5syZs93dkws6F7WkEs7MWgE9gFn5Tg0H3gir1wcYD7QEhhR1lwJcAnzi7rvDyvqa2QIgk9Bdy5Ki4mrVqhXp6ek/thsiIgKY2drCzkV9oN7MagKTgDvDE4CZ3QfkAROPlrn7LHfvCvQCRptZUhFN/wx4PezzXKClu58EPAW8VUg8I8ws3czSt23bVtJuiYhIAaKaVMwsgVBCmejuk8PKhwGDgWu9gEEdd88A9gHdCmm3AdAb+HfYNbvdfW9w/B6QYGYNC2h7rLunuXtacnKBd28iIlJC0Zz9ZcA4IMPdHwsrHwSMAi4MH0cxs9ZmFh8ctwQ6AmsKaf4K4F13PxB2fZPgOzGz3oT6tiOinRIRkSJFc0ylHzAEWGRm84Oye4G/AYnA1CAHzHT3kUB/4B4zOwQcAW519+0AZvYecJO7ZwbtXA08nO/7LgduMbM8YD9wdUF3QSIiEj1Rm1JcHqSlpbkG6kVEjo2ZzXH3tILO6Y16ERGJGCUVERGJGCWVEjhw6DC/f2cJq7fvi3UoIiJlipJKCSzckM3EWesY+NfPuW3iXBZvzI51SCIiZYKSSgn0bl2f6aMGMvL0tny5fBuDn/qKoeO/YdZ3O6jMEx9ERDT76zhnf2XvP8SrM9cy/qvV7Nh3kJ4t6/HHS7rRqUntCEUpIlK2aPZXFNWplsBtZ7Zj+j0DeeiirqzdsY+bXk4nO+dQrEMTESl1SioRkpQQx9C+rRg7NI3N2Qe4+80FehQmIpWOkkqEndyiHvf8pBMfLd3ChOlrYh2OiEipUlKJghv7t+bszo353/czmL8+K9bhiIiUGiWVKDAz/nrFSTSqlcRtE+dqfEVEKg0llSipUz2BZ649ma17DvBrja+ISCWhpBJF3VPrcs9POjN16RbGfbU61uGIiESdkkqUDe/XinO7NObh95cxd92uWIcjIhJVSipRZmY8evlJNKmTxK2vzmXr7gPFXyQiUk4pqZSCOtUTGDskjez9hxjx9zkcOHQ41iGJiESFkkop6dK0No9f1Z3567MYPXmRBu5FpEKK5h71qWb2mZllmNkSM7sjKH/UzJaZ2UIzm2JmdYPy3mY2P/hZYGaXFNLuS2a2Oqxu96DczOxvZrYyaPvkaPWtpAZ1a8KvzunAlHkbGfPld7EOR0Qk4qJ5p5IH/MrdOwOnALeZWRdgKtDN3U8ElgOjg/qLgTR37w4MAsaYWXwhbd/t7t2Dn/lB2U+A9sHPCOC5qPTqON0+sB2DT0zhkQ+W8fHSLbEOR0QkoqKWVNx9k7vPDY73ABlAM3f/yN3zgmozgeZBnZyw8iTgWJ8PXQS84iEzgbpmlnLcHYmwowP33ZrW4Y5/zOPbzXtiHZKISMSUypiKmbUCegCz8p0aDrwfVq+PmS0BFgEjw5JMfn8MHnE9bmaJQVkzYH1YnQ1BWf5YRphZupmlb9u2rUT9OV7VqsYxdmhPqifGc9Mrs9m572BM4hARibSoJxUzqwlMAu50991h5fcRekQ28WiZu89y965AL2C0mSUV0ORooFNQpz4w6miTBdT9wd2Ou4919zR3T0tOTi5hr45fSp1qjB3Sky27c7n8uRnc/X8LeOqTFfxr/kbmrdvFjr25GswXkXKnsDGLiDCzBEIJZaK7Tw4rHwYMBs7yAv7ldPcMM9sHdAPS853bFBzmmtkE4NfB5w1AaljV5kBmpPoSDT1a1OOZa05m7Jer+GL5Nrbuyf3e+XrVE7js5OYM7duKFg2qxyhKEZEfL2pJxcwMGAdkuPtjYeWDCN1dnO7uOWHlrYH17p5nZi2BjsCaAtpNcfdNQfsXExrgB3gbuN3M/gH0AbLDElCZdU6XxpzTpTEA+w8eZsOuHNbuyGHdzhzS1+5kwow1jJu+mjM7NmLYqa0Y0K4hVaoUdFMmIhJ7UdtO2Mz6A9MIjY8cCYrvBf4GJAI7grKZ7j7SzIYA9wCHgvoPuftbQVvvATe5e6aZfQokE3rcNZ/Q2MveIMk8TWjmWA5wg7sXuVdwJLYTjrbN2Qd4bdZaXvtmHdv3HqRNwxoM7duSn/VpQWJ8XKzDE5FKqKjthLVHfRlPKkfl5h3m/UWbeWnGGuavz+KUNvUZOzSN2kkJsQ5NRCoZ7VFfASTGx3Fxj2a8dVs/Hr/qJNLX7OLK579mc7bWEhORskNJpRy6pEdzJtzQi/U7c7jsuRms3Kp3XUSkbFBSKacGtE/mjZv7kpt3hMue+5o5a3fGOiQRESWV8qxbszpMvuVU6teoyjUvzOLDJZtjHZKIVHJKKuVciwbVeXNkXzql1OaWV+fw5pwNsQ5JRCoxJZUKoEHNRF7/eR/6tm3AqEkL+XJ5bJafERFRUqkgqleN5/nretK+UU1unTiXpZm7i79IRCTClFQqkFpJCUy4oRc1E+MZ/tJsNmXvj3VIIlLJKKlUMCl1qjH++l7szc3jhgmz2XPgUKxDEpFKREmlAurStDbPXHsyK7bu5daJczl0+EjxF4mIRICSSgV1eodk/nRJN6at2M79UxZrGX0RKRVRXfpeYuuqXi1Yv3M/T3+2kka1E/nlOR0IrbspIhIdSioV3K/O7cCW3Qd46tOVbNi1n/+99ASSErS6sYhEh5JKBWdm/PnyE2lRvzp/nbqc1dv3MXZoTxrVKmhTTRGR46MxlUrAzPjFWe15/rqT+XbzHi56ejqLN2bHOiwRqYCUVCqRQd1SePOWvhhw+fMzeG9Rmd8YU0TKmaglFTNLNbPPzCzDzJaY2R1B+aNmtszMFprZFDOrG5T3NrP5wc8CM7ukkHYnmtm3ZrbYzMabWUJQfoaZZYe18dto9a0869q0Dv+6vT9dUmpz68S5/OHdpcxYtZ3s/XqfRUSOXzS3E04BUtx9rpnVAuYQ2lO+OfBpsBf9IwDuPsrMqgMHg/IUYAHQ1N3z8rV7PvB+8PE14Et3f87MzgB+7e6Df2yM5Wnnx0jLzTvM/VMW839hC1C2bFCdbs3qcEKzOpzcoh69WtXTbDER+YGidn6M2kC9u28CNgXHe8wsA2jm7h+FVZsJXB7UyQkrTwIKzHbu/t7RYzP7hlCSkmOUGB/Ho1ecxD0/6cTizN0s3pjN4o3ZLFifxb8Xhh6LDe/Xmvsv6EyVKkosIvLjlMrsLzNrBfQAZuU7NRx4I6xeH2A80BIYkv8uJV+bCcAQ4I6w4r5mtgDIJHTXsqSA60YAIwBatGhRgt5ULA1qJnJ6h2RO75D8n7Jd+w7y5CcrGD99NTv35fLny0+iaryG30SkeFH/l8LMagKTgDvdfXdY+X1AHjDxaJm7z3L3rkAvYLSZFTXv9VlCj76mBZ/nAi3d/STgKeCtgi5y97HunubuacnJyQVVqfTq1ajKgz/twt3ndeSt+Zn8/JV0cg4Wmt9FRP4jqkkluJuYBEx098lh5cOAwcC1XsCgjrtnAPuAboW0+yCQDPwy7Jrd7r43OH4PSDCzhhHsTqViZtx2ZjsevvQEpq3YxjUvzGLXvoOxDktEyrhozv4yYByQ4e6PhZUPAkYBF4aPo5hZazOLD45bAh2BNQW0exNwHvAzdz8SVt4k+E7MrDehvu2IQtcqlat7t+C563qydNNurhjzNZlZWk5fRAoXzTuVfoTGPAaGTfM9H3gaqAVMDcqeD+r3BxaY2XxgCnCru28HMLP3zKxpUO95oDHwdb6pw5cDi4Mxlb8BVxd0FyTH7ryuTXhleG+2ZB/gsudmsHLrnliHJCJlVNSmFJcHlXlKcUksycxm2PjZHD5yhAk39KZ7at1YhyQiMVDUlGJN6ZEfrWvTOky6pS+1khK45oWZTFuxLdYhiUgZo6Qix6Rlgxq8ObIvLepXZ/hLs3l3YWasQxKRMkRJRY5Zo9pJvHFzX3qk1uMXr8/j7zPXxjokESkjlFSkROpUS+CVG3tzVqdGPPDWYp78eIV2lxQRJRUpuaSEOJ6/rieXndycxz9ezs9fmcP89VmxDktEYkibdMlxiY+rwl+uOJE2yTUY88UqLn5mC31a1+fm09twRodGWjdMpJLRlGJNKY6Yvbl5vDF7PeOmfUdm9gE6NK7Jzwe04aLuzbR2mEgFUtSUYiUVJZWIO3T4CO8uzGTMF9+xbPMeOjauxcSf96FhzcRYhyYiEaD3VKRUJcRV4ZIezXn/jgGMGdKTdTtzuOaFmezYmxvr0EQkypRUJGrMjPO6NmHc9Wms25nDtS/OYqcWpRSp0JRUJOpObduQccN6sXr7Pq55YaZWOxapwJRUpFT0a9eQF4am8d32fVz74iyycpRYRCoiJRUpNad1SOaFoWms3LZXiUWkglJSkVJ1eodkxgzpyYotocSyZvu+WIckIhGkpCKl7syOjRgzNDQrbNCTX/LS9NUcOVJ5p7aLVCRKKhITZ3ZsxNS7TueUNg343TtL+dkLM1m3I6f4C0WkTIvmdsKpZvaZmWWY2RIzuyMof9TMlpnZQjObYmZ1g/LeYTtELjCzSwppt7WZzTKzFWb2hplVDcoTg88rg/OtotU3iYwmdZKYcH0v/nzZiSzN3M2gJ7/k7zPX6q5FpByL5p1KHvArd+8MnALcZmZdgKlAN3c/EVgOjA7qLwbS3L07MAgYc3TP+nweAR539/bALuDGoPxGYJe7twMeD+pJGWdmXNkrlQ/uOo2eLevxwFuLuW7cLL5ZvVOrHouUQ1FLKu6+yd3nBsd7gAygmbt/5O55QbWZQPOgTk5YeRLwg39RzMyAgcCbQdHLwMXB8UXBZ4LzZwX1pRxoVrcarwzvzZ8uOYHFG7O5cszXnPt4aLwle/+hWIcnIj9SqYypBI+iegCz8p0aDrwfVq+PmS0BFgEjw5LMUQ2ArLDyDUCz4LgZsB4gOJ8d1Jdywsy4pk8LZt17Nn++/ESqJ8bzu3eW0udPH/ObNxcwf30WB/OOxDpMESlC1Je+N7OawCTgTnffHVZ+H6FHZBOPlrn7LKCrmXUGXjaz9939QHhzBXyF/4hz4fGMAEYAtGjR4hh7I6WhWtU4rkxL5cq0VBZvzGbirLX8a34m/0zfAEBCnFEjMZ4aVeOpXjWOGonxtGtUk7M6NWJAh2RqJmpHB5FYieoqxWaWALwLfOjuj4WVDwNGAme5e4FTfszsM+Bud08PKzNgG9DE3fPMrC/wO3c/z8w+DI6/DsZiNgPJXkQHtUpx+bHnwCE+WLyZLbsPsO/gYXJy89h38DD7cvPYm5vHwg3ZZO8/RNW4KvRpU5+zOzdmYKdGpNavHuvQRSqcolYpjtqvdEECGAdk5Esog4BRwOnhCcXMWgPrg2TREugIrAlv0909SDaXA/8AhgH/Ck6/HXz+Ojj/aVEJRcqXWkkJXJGWWuj5vMNHSF+7i08ytvBJxlYefHsJD769hIGdGvHC0DTitFmYSKmI2p2KmfUHphEaHzn6IPxe4G9AIrAjKJvp7iPNbAhwD3AoqP+Qu78VtPUecJO7Z5pZG0IJpT4wD7jO3XPNLAn4O6Gxm53A1e7+XVEx6k6l4vpu217enLOBZz9fxYM/7cIN/VrHOiSRCkObdBVCSaVic3eunzCb9DU7mfrL02lat1qsQxKpELRJl1RKZsb/XNyNIw6//ddivfciUgqUVKRCS61fnbvOac/HGVv5YPHmWIcjUuEpqUiFN7xfa7qk1ObBt5ew+4BepBSJJiUVqfDi46rw8GUnsH1vLn/+YFmswxGp0JRUpFI4sXldrj+1NRNnrWPO2p2xDkekwlJSkUrjV+d2IKV2EqMnL9JyLyJRoqQilUaNxHgeuqgby7fs5YVpRb7CJCIlpKQilcrZXRpz/glNeHzqcka9uZDlW/bEOiSRCkUr70ml86dLTqBe9apMmruBN9LXc3qHZG4a0Jr+7Rqi3RJEjo/eqNcb9ZXWzn0HeW3WWl7+ei3b9uTSsXEtbhzQmstObq61wkSKoDfqRQpQv0ZVbh/Ynq9GnclfrjgJM/jNmwu5beJcDhw6HOvwRMolJRWp9BLj47i8Z3Pev2MADwzuwgdLNnP9hG/YoxclRY6ZkopIwMy4sX9rnriqO+lrdnH12Jls25Mb67BEyhUlFZF8Lu7RjBeGpbFq214uf34G63YUuI+ciBRASUWkAGd2bMTEm04hK+cQlz0/g4xNu4u/SESUVEQK07NlPf5vZF/izLhyzNe8PGMNm7L3xzoskTJNU4o1pViKsTFrPyNeSWdJZuhu5YRmdTinS2PO6dKYTk1q6d0WqXRisvOjmaUCrwBNCG0PPNbdnzSzR4GfAgeBVcAN7p5lZucADwNVg3N3u/unBbT7BqH96wHqAlnu3t3MWgEZwLfBuZnuPrKoGJVU5Mdyd1Zu3cvUjC1MXbqF+euzcIfm9aox5JSW/HxAG6ro3RapJEqcVMzsOnd/NTju5+7Tw87d7u5PF3FtCpDi7nPNrBYwB7gYaA586u55ZvYIgLuPMrMewJZgH/puwIfu3qyYjv0VyHb3h4Kk8q67dyvqmnBKKlJSW/cc4NOMrbyzMJPpK3dwbpfG/PXKk6iVlBDr0ESi7nhefvxl2PFT+c4NL+pCd9/k7nOD4z2E7iKauftH7p4XVJtJKMng7vPcPTMoXwIkmVliYe1b6JnDlcDrxfRBJOIa1Uri6t4tePXGPvx2cBc+WbaVi5+Zzqpte2MdmkhMFZdUrJDjgj4X3kjoLqIHMCvfqeHA+wVcchkwz92LeklgAKE7mxVhZa3NbJ6ZfWFmAwqJZYSZpZtZ+rZt235sF0QKZGYM79+aV2/sw66cQ1z89HSmLt0S67BEYqa4pOKFHBf0uUBmVhOYBNzp7rvDyu8D8oCJ+ep3BR4Bbi6m6Z/x/buUTUALd+9B6A7rNTOrnf8idx/r7mnunpacnPxjuiBSrL5tG/DOL/rTqmENfv5KOk98vJwjRyrvJBipvIpbpbiTmS0kdFfSNjgm+NymuMbNLIFQQpno7pPDyocBg4GzPGxQx8yaA1OAoe6+qoh244FLgZ5Hy4K7mtzgeI6ZrQI6ABo0kVLRrG41/m9kX+6dsognPl7BjFU7uOCEFPq1a0Db5JqaJSaVQnFJpXNJGw7GPMYBGe7+WFj5IGAUcLq754SV1wX+DYwOnxBQiLOBZe6+Iez6ZGCnux82szZAe0A7MUmpSkqI469XnESP1LqM+fI7Hnx7CQCNaiVyatsGnNquIae2bUCzutWUZKRCOqYpxWbWADgNWOfuc4qp2x+YBiwiNKUY4F7gb0AisCMom+nuI83sfmA0ED5Gcq67bzWzF4Hn3T09aPul4Lrnw77vMuAhQo/UDgMPuvs7RcWo2V8Sbet25DB91XZmrNrB16u2s33vQQAa106kR2o9ureoS4/UupzQvA7Vq2p7IykfjmdK8bvAPe6+OJgiPJfQ46S2hN47eSIaAZcWJRUpTe7Ot1v2MHPVDuavz2Le+izWBuuKxVUxuqTU5r4LOnNKmwYxjlSkaEUlleJ+NWrt7ouD4xuAqe4+NHjvZDpQrpOKSGkyMzo1qU2nJv+dP7Jjby7z12cxf30W7yzIZOj4b3j2mpM5u0vjGEYqUnLFzf4K31DiLOA9+M97J0cKvEJEfrQGNRM5q3NjfnVuR6bc2o/OTWpx86tzmDJvQ/EXi5RBxSWV9Wb2CzO7BDgZ+ADAzKoBenVYJILq1ajKxJ+fQp/W9bnrjQW8NH11rEMSOWbFJZUbga7A9cBV7p4VlJ8CTIhiXCKVUs3EeMZf34tzuzTmd+8s5cmPV1CZF32V8qfIMRV33wr8YFFGd/8M+CxaQYlUZkkJcTx77cncM3kRj3+8nKz9B3nggi5asFLKhSKTipm9XdR5d78wsuGICEB8XBX+fNmJ1E5KYPz01eQddh66qKvebZEyr7jZX32B9YSWQ5nFMaz3JSLHp0oV44HBnUmIN8Z88R3N6lVj5OltYx2WSJGKSypNgHMIrbN1DaE33l939yXRDkxEQtOQR53XiU1ZB3j4/WU0rVuNC09qGuuwRApV5EC9ux929w/cfRihwfmVwOdm9otSiU5EqFLFePSKE+nTuj6//ucCZn23o/iLRGKk2D3qzSzRzC4FXgVuI7TMyuSirxKRSEqMj2PskDRS61fj56+ks3LrnliHJFKgIpOKmb0MzCD0jsrv3b2Xu//B3TeWSnQi8h91qifw0g29qRofx/UTZrN1z4FYhyTyA8XdqQwhtHz8HcAMM9sd/Owxs93FXCsiEZZavzrjr09jx96D3PhSOjkH84q/SKQUFTemUsXdawU/tcN+arn7DzbAEpHoO7F5XZ6+pgdLMrO5fsJsNmbtj3VIIv9R7JiKiJQ9Z3VuzGNXdmfxxmzOe/xLXpu1Tm/eS5mgpCJSTl3coxkf3nkaJzavw71TFjFk3Dds2JVT/IUiUaSkIlKOpdavzsSb+vCnS05g3rpdnPf4l7w6c63uWiRmopZUzCzVzD4zswwzW2JmdwTlj5rZMjNbaGZTgm2EMbNzzGyOmS0K/hxYSLu/M7ONZjY/+Dk/7NxoM1tpZt+a2XnR6ptIWWJmXNOnBR/edRont6zH/W8t5uqxM5m/Pqv4i0Ui7Ji2Ez6mhkM7Raa4+9xgU685wMVAc+BTd88zs0cA3H2UmfUAtrh7ppl1Az5092YFtPs7YK+7/yVfeRdCy8n0BpoCHwMd3P1wYTFq50epaNydN2av588ffsvOfQc5q1Mj7jqnA92a1Yl1aFKBFLXzY9TuVNx9k7vPDY73ABlAM3f/yN2PzoOcSSjJ4O7z3D0zKF8CJJlZ4jF85UXAP9w9191XE3r7v3ck+iJSXpgZV/duwbTfnMnd53Ukfe0uBj/1FTf/PZ1lm/UWgERfqYypmFkroAehRSnDDQfeL+CSy4B57p5bSJO3B4/PxptZvaCsGaHFL4/aEJTlj2WEmaWbWfq2bduOoRci5UeNxHhuO7Md00adyZ1nt2fGyh385Mlp3P7aXNbt0GC+RE/Uk4qZ1QQmAXe6++6w8vuAPGBivvpdgUeAmwtp8jmgLdAd2AT89eilBdT9wbM9dx/r7mnunpacnHyMvREpX2onJXDn2R2YNupMbj2jLZ9kbOXsx7/g0Q+XsS9XL05K5EU1qZhZAqGEMtHdJ4eVDwMGA9d62KCOmTUHpgBD3X1VQW26+5ZgocsjwAv89xHXBiA1rGpzIDP/9SKVUd3qVbn7vE589uszGHxCCs98tooz//I5k+du4MgRzRSTyInm7C8DxgEZ7v5YWPkgYBRwobvnhJXXJbS0/mh3n15EuylhHy8BFgfHbwNXBwtgtgbaA99Eqj8iFUGTOkk8dlV3Jt96Kil1q/HLfy7g0udmMG/drliHJhVENO9U+hFaO2xgvum/TwO1gKlB2fNB/duBdsADYfUbAZjZi2Z2dKbBn4NpxwuBM4G7AII9Xv4JLAU+AG4rauaXSGV2cot6TLnlVP56xUlszNrPJc/O4Jf/nM/W3VqkUo5P1KYUlweaUiwCe3PzeOazlYybtpqq8VX4xcB23NCvNVXj9W60FCwmU4pFpHyomRjPqEGd+PCu0+jTuj7/+/4yBj3xJZ9/uzXWoUk5pKQiIgC0bliDcdf3YsL1vXDg+gmzuenl2azfqSnI8uMpqYjI95zZqREf3nkao3/Sia9X7eCnT3/F16u0hbH8OEoqIvIDVeOrcPPpbXnvjgEk10xkyLhZvDZrXazDknJASUVECtWyQQ0m3Xoq/ds35N4pi/j9O0vIO3wk1mFJGaakIiJFqp2UwLhhvbixf2smTF/D8JfT2X3gUKzDkjJKSUVEihVXxXhgcBcevvQEZqzczqXPzmDtjn2xDkvKICUVEfnRru7dgr/f2Ifte3O55oVZHDik94vl+5RUROSY9G3bgOev68nGrP28OO27WIcjZYySiogcs1PaNGBQ1yY8+/kqtmhpFwmjpCIiJTL6/E7kHXYe/fDbWIciZYiSioiUSMsGNbihXysmzd3Aog3ZsQ5HygglFREpsdsGtqN+9ar84d2lVObFaeW/lFREpMRqJyXwy3M78M2anXyweHOsw5EyQElFRI7LVWmpdGpSiz+9n6EpxqKkIiLHJz6uCvdf0IX1O/fz0ow1sQ5HYiya2wmnmtlnZpZhZkvM7I6g/FEzW2ZmC81sSrCNMGZ2jpnNCXZ1nGNmAwtpt7DrW5nZ/rBdI58v6HoRibz+7RtydudGPP3pSrbtyY11OBJD0bxTyQN+5e6dgVOA28ysCzAV6ObuJwLLgdFB/e3AT939BGAY8PdC2i3seoBV7t49+BkZ+S6JSGHuPb8zBw4d5rGpmmJcmUUtqbj7JnefGxzvATKAZu7+kbvnBdVmAs2DOvPcPTMoXwIkmVliAe0WeL2IxFab5JoM7duK179ZzxMfL+fIEc0Gq4xKZUzFzFoBPYBZ+U4NB94v4JLLgHnuXtx9dP7rW5vZPDP7wswGlDBcESmh3wzqyKU9mvHExyu4deJc9uXmFX+RVChRTypmVhOYBNzp7rvDyu8j9IhsYr76XYFHgJuLaTf/9Z/7OLwAABL9SURBVJuAFu7eA/gl8JqZ1S7guhFmlm5m6du2bSt5x0TkB5IS4vjrlSdx/wWd+WjpZi59dgbrdmg74sokqknFzBIIJZSJ7j45rHwYMBi41sPemDKz5sAUYKi7ryqi3R9c7+657r4jOJ4DrAI65L/W3ce6e5q7pyUnJ0eimyISxsy4aUAbXh7em827D3DhM18xfeX2WIclpSSas78MGAdkuPtjYeWDgFHAhe6eE1ZeF/g3MNrdpxfRbmHXJ5tZXHDcBmgPaAlVkRgZ0D6Zt2/vR6NaiQwd/w3jvlqtt+4rgWjeqfQDhgADw6b5ng88DdQCpuab+ns70A54IKx+IwAze9HM0oJ6hV1/GrDQzBYAbwIj3X1nFPsnIsVo2aAGk2/tx1mdGvGHd5fy+3eWagC/grPK/JtDWlqap6enxzoMkQrvyBHnj+9lMO6r1VzaoxmPXH4iCXF697q8MrM57p5W0Ln40g5GRCqfKlWM+y/oTL3qCfzlo+XsPnCIp685maSEuFiHJhGmXxVEpFSYGbcPbM8fLurKJ8u2Mmz8N+w5cCjWYUmEKamISKka0rcVT1zVnTlrd/GzF2ayY6+WdalIlFREpNRd1L0ZY4f2ZMWWvVwx5ms2Ze+PdUgSIUoqIhITAzs15u839mHr7lxumDBbb99XEEoqIhIzvVvX55lrT2b5lj3c9cZ8TTeuAJRURCSmTu+QzAODu/DR0i385SOtcFzeaUqxiMTc9ae2YvmWvTz7+SraN67JJT20+Hh5pTsVEYk5M+Ohi7pySpv6jJq0iLnrdsU6JCkhJRURKRMS4qrw3LU9SamTxIhX5pCZpRlh5ZGSioiUGfVqVGXcsDRyDx3mppfTyTmoGWHljcZURKRMadeoFk9d04PhL82m38Ofklq/Ok3rVCOlbtJ//6xbjaZ1qpFcK5G4KhbrkCWMkoqIlDlndGzE2CFpfJyxhczsA6zctpdpK7ax7+Dh79WLr2I0rp1E07pJpNSpRuuGNbiqVypN61aLUeSiVYq1SrFIueDu7N6fR2b2fjKz9pOZfYBNWfvZlH0g+Lyfjbv2E1fFuLh7M24+vS3tGtWMddgVklYpFpFyz8yoUz2BOtUT6Jzyg53CAdiwK4cXvvyOf8xez5tzN3BelybcckZbTkqtW8rRVl66U9GdikiFs31vLi9NX8MrX69h94E8TmlTn54t65FSpxpNgzGZlDrVqJ0UT2iTWjkWRd2pKKkoqYhUWHsOHOL1b9bx+jfrWbczh8P5loGpmRjPGR2TGXl6W7o1qxOjKMufmCQVM0sFXgGaAEeAse7+pJk9CvwUOAisAm5w9ywzOwd4GKganLvb3T8toN36wBtAK2ANcKW777LQrxtPAucDOcD17j63qBiVVEQqj8NHnK17DpCZdYBN2fvZlHWA1Tv28fb8TPbm5jGgfUNuOaMtfds00N1LMWKVVFKAFHefa2a1gDnAxUBz4FN3zzOzRwDcfZSZ9QC2uHummXUDPnT3ZgW0+2dgp7s/bGb3APWC688HfkEoqfQBnnT3PkXFqKQiItn7DzFx1lrGf7WG7XtzOal5HW45oy3ndmlCFU1XLlCZePxlZv8Cnnb3qWFllwCXu/u1+eoasB1o6u65+c59C5zh7puCxPW5u3c0szHB8ev56xUWk5KKiBx14NBhJs3dwJgvvmPdzhxqJ8XTKaU2nZrUomOTWnRqUpuOTWpRM1Hzm2I++8vMWgE9gFn5Tg0n9Cgrv8uAefkTSqDx0UQRJJZGQXkzYH1YvQ1B2feSipmNAEYAtGjR4pj6ISIVV1JCHNf2aclVaal8uGQL01dtZ9mm3Uyeu5G9YXu9pLWsx0vDeyu5FCLq/6uYWU1gEnCnu+8OK78PyAMm5qvfFXgEOPdYv6qAsh/chrn7WGAshO5UjvE7RKSCi4+rwgUnpnDBiSlA6P2YDbv28+3mPSzcmM0zn61k1JsLefqaHhp7KUBUk4qZJRBKKBPdfXJY+TBgMHCWhz1/M7PmwBRgqLuvKqTZLWaWEvb4a2tQvgFIDavXHMiMXG9EpDIyM1LrVye1fnXO7tKY6lXjePj9ZZw8vR439m8d6/DKnKgtKBmMi4wDMtz9sbDyQcAo4EJ3zwkrrwv8Gxjt7tOLaPptYFhwPAz4V1j5UAs5BcguajxFRKQkbj6tDed0acz/vpdB+pqdsQ6nzInmKsX9gCHAQDObH/ycDzwN1AKmBmXPB/VvB9oBD4TVbwRgZi+a2dFBoYeBc8xsBXB0GjLAe8B3wErgBeDWKPZNRCopM+MvV5xEs3rVuO21uWzfW9DQb+Wllx81+0tESmBJZjaXPjuDtFb1eGV4n0q1WnJRs7+0n4qISAl0bVqHP1zcjekrd/D41OWxDqfMUFIRESmhK9NSuSotlac/W8knGVtiHU6ZoKQiInIcfn9RV7qk1OauN+YzZd4GjhypvEMKoKQiInJckhLiGDOkJy0aVOeuNxYw+KmvmLZiW6zDihklFRGR45Ravzpv39afJ6/uzu4Dhxgy7huGjJvF4o3ZsQ6t1Gn2l2Z/iUgE5eYd5tWZ63jq0xVk5Rziou5NOblFPapXjaNmYjzVE+OpUTWO6lXjaZNcg6SEuFiHfMzKxIKSZZGSiohES/b+Qzz/xSomTF/NgUNHCqzTsGZVbujXmutOaUmdagmlHGHJKakUQklFRKItN+8wew/kkXPwMPsO5rEv9zA5B/PIyjnEpLkb+PzbbdSoGse1p7Tkxv6taVw7KdYhF0tJpRBKKiISa0szdzPmy1W8syCT+CpVuKRHM0ac3oa2yTVjHVqhlFQKoaQiImXF+p05vDDtO96YvZ6Dh49w/gkp3H5mOzqn1I51aD+gpFIIJRURKWu2781l/FereeXrtezNzePszo247cx29GhRL9ah/YeSSiGUVESkrMrOOcTLX69h/PTVZOUcon+7htxxdnt6taof69C09peISHlTp3oC/++s9nw1aiD3nt+JZZv3cNWYr5m6tGwvB6OkIiJShtVMjGfEaW354u4zOKFZHX7x+lzmr8+KdViFUlIRESkHaiTG8+KwXiTXSuTGl2azdse+WIdUICUVEZFyIrlWIi/d0JvD7lw/YTY79x2MdUg/EM3thFPN7DMzyzCzJWZ2R1D+qJktM7OFZjYl2EYYM2sQ1N9rZk8X0e4bYTtDrjGz+UF5KzPbH3bu+cLaEBEpr9om1+TFoWlszNrPTS/P5sChw7EO6XuieaeSB/zK3TsDpwC3mVkXYCrQzd1PBJYDo4P6B4AHgF8X1ai7X+Xu3d29OzAJmBx2etXRc+4+MsL9EREpE9Ja1eeJq7ozb30Wd/5jPofL0HL7UUsq7r7J3ecGx3uADKCZu3/k7nlBtZlA86DOPnf/ilByKZaZGXAl8HrEgxcRKePOPyGF+87vzAdLNvM//15K3uGC1xcrbfGl8SVm1groAczKd2o48EYJmx0AbHH3FWFlrc1sHrAbuN/dpxUQywhgBECLFi1K+NUiIrF304A2bMzaz4Tpa3hzzgb6tW3IaR2SOa1DQ5rXqx6TmKKeVMysJqHHVHe6++6w8vsIPSKbWMKmf8b371I2AS3cfYeZ9QTeMrOu4d8J4O5jgbEQevmxhN8tIlImPHBBF/q0rs/n327jy+Xb+GDJZgDaNKzBaR2S6d+uIX3a1KdWUumsghzVpGJmCYQSykR3nxxWPgwYDJzlJXil38zigUuBnkfL3D0XyA2O55jZKqADoFfmRaTCqlLFGNQthUHdUnB3Vm3bx5fLtzFtxTbemL2el2asIa6K0T21Lv3aNaR/u4Z0T61L1fjojH5ELakEYx7jgAx3fyysfBAwCjjd3XNK2PzZwDJ33xDWbjKw090Pm1kboD3wXYk7ICJSzpgZ7RrVpF2jmgzv35rcvMPMXZvF9JXbmbZyO09/uoK/fbKC6lXjuKZ3C+4f3CXiMUTzTqUfMARYdHTaL3Av8DcgEZgayjvMPDpTy8zWALWBqmZ2MXCuuy81sxeB59396F3H1fxwgP404CEzywMOAyPdfWfUeiciUsYlxsfRt20D+rZtwK/P60h2ziG+/m4HX63cRtO61aLynVpQUgtKiogcEy0oKSIipUJJRUREIkZJRUREIkZJRUREIkZJRUREIkZJRUREIkZJRUREIkZJRUREIqZSv/xoZtuAtcfRRENge4TCKU/U78pF/a5cfky/W7p7ckEnKnVSOV5mll7YW6UVmfpduajflcvx9luPv0REJGKUVEREJGKUVI7P2FgHECPqd+Wiflcux9VvjamIiEjE6E5FREQiRkmlBMxskJl9a2YrzeyeWMcTLWY23sy2mtnisLL6ZjbVzFYEf9aLZYzRYGapZvaZmWWY2RIzuyMor9B9N7MkM/vGzBYE/f59UN7azGYF/X7DzKrGOtZoMLM4M5tnZu8GnytLv9eY2SIzm29m6UFZif+uK6kcIzOLA54BfgJ0AX5mZpHfk7NseAkYlK/sHuATd28PfBJ8rmjygF+5e2fgFOC24P/jit73XGCgu58EdAcGmdkpwCPA40G/dwE3xjDGaLoDyAj7XFn6DXCmu3cPm0pc4r/rSirHrjew0t2/c/eDwD+Ai2IcU1S4+5dA/i2ZLwJeDo5fBi4u1aBKgbtvcve5wfEeQv/QNKOC991D9gYfE4IfBwYCbwblFa7fAGbWHLgAeDH4bFSCfhehxH/XlVSOXTNgfdjnDUFZZdHY3TdB6B9foFGM44kqM2sF9ABmUQn6HjwCmg9sBaYCq4Asd88LqlTUv+9PAL8BjgSfG1A5+g2hXxw+MrM5ZjYiKCvx3/X4KARY0VkBZZpCVwGZWU1gEnCnu+8O/fJasbn7YaC7mdUFpgCdC6pWulFFl5kNBra6+xwzO+NocQFVK1S/w/Rz90wzawRMNbNlx9OY7lSO3QYgNexzcyAzRrHEwhYzSwEI/twa43iiwswSCCWUie4+OSiuFH0HcPcs4HNCY0p1zezoL6AV8e97P+BCM1tD6HH2QEJ3LhW93wC4e2bw51ZCv0j05jj+riupHLvZQPtgZkhV4Grg7RjHVJreBoYFx8OAf8UwlqgInqePAzLc/bGwUxW672aWHNyhYGbVgLMJjSd9BlweVKtw/Xb30e7e3N1bEfrv+VN3v5YK3m8AM6thZrWOHgPnAos5jr/revmxBMzsfEK/ycQB4939jzEOKSrM7HXgDEKrlm4BHgTeAv4JtADWAVe4e/7B/HLNzPoD04BF/PcZ+72ExlUqbN/N7ERCg7JxhH7h/Ke7P2RmbQj9Bl8fmAdc5+65sYs0eoLHX79298GVod9BH6cEH+OB19z9j2bWgBL+XVdSERGRiNHjLxERiRglFRERiRglFRERiRglFRERiRglFRERiRglFamUzGxv8GcrM7smwm3fm+/zjEi2H2lmdr2ZPR3rOKRiUFKRyq4VcExJJVipuijfSyrufuoxxlSu/Ij/PaQSUVKRyu5hYECwl8RdwYKKj5rZbDNbaGY3Q+iluGCPldcIvRSJmb0VLMK35OhCfGb2MFAtaG9iUHb0rsiCthcH+1dcFdb252b2ppktM7OJVsBCY0GdR4I9T5ab2YCg/Ht3Gmb27tE1rMxsb3DNHDP72Mx6B+18Z2YXhjWfamYfWGifoAfD2rou+L75ZjbmaAIJ2n3IzGYBfSP1f4aUf1pQUiq7ewjeoAYIkkO2u/cys0Rgupl9FNTtDXRz99XB5+HuvjNY0mS2mU1y93vM7HZ3717Ad11KaJ+SkwitUjDbzL4MzvUAuhJaX2o6ofWoviqgjXh37x2s6vAgoaVUilID+NzdR5nZFOB/gHMI7QX0Mv9dYqg30A3ICeL6N7APuIrQgoOHzOxZ4FrglaDdxe7+22K+XyoZJRWR7zsXONHMjq75VAdoDxwEvglLKAD/z8wuCY5Tg3o7imi7P/B6sBLwFjP7AugF7A7a3gAQLD3fioKTytHFLecEdYpzEPggOF4E5AYJYlG+66e6+47g+ycHseYBPQklGYBq/HdhwcOEFtwU+R4lFZHvM+AX7v7h9wpDj5P25ft8NtDX3XPM7HMg6Ue0XZjwNaUOU/h/m7kF1Mnj+4+yw+M45P9di+nI0evd/UjYCrzww2XdPYj3ZXcfXUAcB4LkKPI9GlORym4PUCvs84fALcHS95hZh2D11vzqALuChNKJ0BLxRx06en0+XwJXBeM2ycBpwDcR6MMaQnugVDGzVEKPso7VORbal7waoV3+phPaRvZyC+2zcXTf8pYRiFcqMN2pSGW3EMgzswXAS8CThB4LzQ0Gy7dR8FaqHwAjzWwh8C0wM+zcWGChmc0NllA/agqhQe0FhO4EfuPum4OkdDymA6sJPd5aDMwtQRtfAX8H2hFaqTYdwMzuJ7QrYBXgEHAbsPY445UKTKsUi4hIxOjxl4iIRIySioiIRIySioiIRIySioiIRIySioiIRIySioiIRIySioiIRIySioiIRMz/BybHKdspS22sAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%pylab inline\n",
    "plot(range(50), stoch_errors_by_iter[:50])\n",
    "xlabel('Iteration number')\n",
    "ylabel('MSE')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'MSE')"
      ]
     },
     "execution_count": 102,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEGCAYAAACKB4k+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAe2klEQVR4nO3deZgddZ3v8ff3dKf3Jd2d7nQnnRUSSNhCCMgiXhxAwS04bjg6gMvFBRyXe2ce0HmeGe+Mz4PXOy5ct4uK4oyIijgyiCBiABFZOiEkgQTShoR01k46pLN2evneP+p3OiehF0j6dJ1z6vN6nvNU1a/qVH9PWfhJbb8yd0dERAQgFXcBIiKSOxQKIiIySKEgIiKDFAoiIjJIoSAiIoOK4y7geEyaNMlnzpwZdxkiInll6dKlO9y9cah5eR0KM2fOpK2tLe4yRETyipltGG6eTh+JiMgghYKIiAxSKIiIyCCFgoiIDFIoiIjIIIWCiIgMUiiIiMigRIbCmq3dfOX+Nby8/1DcpYiI5JREhsL6Hfv51pK/0LHrQNyliIjklESGwuSaUgC27zkYcyUiIrklkaHQVFMGwLbunpgrERHJLYkMhUlVJQDs3KtQEBHJlMhQKC0uoqq0mJ37dKFZRCRTIkMBoL6yhC6FgojIERIbCg1VJezQ6SMRkSMkNhSaqkvp3KNQEBHJlNhQaKwuZbtCQUTkCMkNhaoyXt7fS2//QNyliIjkjMSGwqTq9G2putgsIpKW2FCor4hCYZf6PxIRGZTYUJgYQkG3pYqIHJbYUEg/1azbUkVEDktsKDRVR/0f6bZUEZHDEhsKNeXFlBSndFuqiEiGxIaCmdFUXcr2bnWfLSKSlthQgOipZh0piIgclrVQMLNpZrbEzFab2bNm9unQXm9mD5jZ2jCsC+1mZjebWbuZrTCzhdmqLa2pukyhICKSIZtHCn3A/3D3ecC5wHVmNh+4AXjQ3ecAD4ZpgMuBOeFzLfCdLNYGQHNtGdt26/SRiEha1kLB3be4+7IwvgdYDUwFFgO3hcVuA64I44uBH3vkcWCimbVkqz6I+j/a09PHwd7+bP4ZEZG8MS7XFMxsJnAm8AQw2d23QBQcQFNYbCqwMeNrHaHt6HVda2ZtZtbW2dl5XHU1VkXvatZtqSIikayHgplVAb8EPuPu3SMtOkSbv6LB/RZ3X+TuixobG4+rtsbqEAp6gE1EBMhyKJjZBKJA+Im73xWat6VPC4Xh9tDeAUzL+HorsDmb9aVDQbeliohEsnn3kQE/AFa7+1czZt0NXB3GrwZ+ndF+VbgL6Vxgd/o0U7ZMromeatYdSCIikeIsrvsC4G+BlWa2PLR9HrgJ+LmZfQR4CXhPmHcv8BagHdgPfCiLtQHQUFlCccrYpiMFEREgi6Hg7o8y9HUCgIuHWN6B67JVz1BSKWNyTRmbX1YoiIhAwp9oBphWX86GnfviLkNEJCckPhSmTCxnqx5gExEBFAo010RdXQwMvOLuVxGRxFEo1JbRN+Ds1BvYREQUCs3htlSdQhIRUSjQUlsOwObdB2KuREQkfgqFidGRwpaXFQoiIokPhYbKEkqKU2zW6SMREYWCmdFSW8ZmHSmIiCgUAFpqy3ShWUQEhQIQXWzeqv6PREQUChD1lrqt+6AeYBORxFMoAFMmltHb7+zQy3ZEJOEUCsCU8KzCJl1sFpGEUygAU+sUCiIioFAAYFp9BQAvde2PuRIRkXgpFICq0mLqKibQsUtHCiKSbAqFoKVW71UQEVEoBK115Tp9JCKJp1AIpkwsZ5uOFEQk4RQKQVNNKXt6+th/qC/uUkREYqNQCNLPKqhjPBFJMoVCcPhZBZ1CEpHkUigEUybqSEFERKEQTK4upShlbNKzCiKSYAqFoLgoRXONXrYjIsmmUMgwdWI5HQoFEUkwhUKGqXXlOn0kIommUMgwrb6CLbsPcKhvIO5SRERioVDI0DqxnAGHbXo1p4gklEIhQ/q2VL1XQUSSSqGQYcrEMkDPKohIcikUMqSPFLaoYzwRSSiFQoayCUU0VJboZTsiklgKhaO01pXTsUvvVRCRZFIoHGV6QyUbdioURCSZshYKZnarmW03s1UZbf9sZpvMbHn4vCVj3o1m1m5mz5vZm7NV12hm1Few6eUD9PbrWQURSZ5sHin8CLhsiPavufuC8LkXwMzmA1cCp4TvfNvMirJY27Cm11fQP+BsURfaIpJAWQsFd38E6HqViy8G7nD3Hnd/EWgHzslWbSOZ3lABwIaufXH8eRGRWMVxTeF6M1sRTi/VhbapwMaMZTpC2yuY2bVm1mZmbZ2dnWNe3PT6EAq6riAiCTTeofAd4ARgAbAF+LfQbkMs60OtwN1vcfdF7r6osbFxzAtsrimjpDjFxi6Fgogkz7iGgrtvc/d+dx8AvsfhU0QdwLSMRVuBzeNZW1oqZUyrK9eRgogk0riGgpm1ZEy+E0jfmXQ3cKWZlZrZLGAO8OR41pZpen0FG3SkICIJVJytFZvZT4GLgElm1gH8E3CRmS0gOjW0HvgYgLs/a2Y/B54D+oDr3L0/W7WNZkZDJU+t34W7YzbUmS0RkcKUtVBw9/cP0fyDEZb/EvClbNXzWkyvr2BvTx9d+w7RUFUadzkiIuNGTzQPIX0H0ks6hSQiCaNQGMI0hYKIJJRCYQh6VkFEkkqhMITykiKaa8pYv1NPNYtIsigUhjGjoYKXdKQgIgmjUBjGzIZKHSmISOIoFIYxc1IlO/YeYs/B3rhLEREZNwqFYZzQWAlA+/a9MVciIjJ+FArDOLm5BoDnt+6JuRIRkfGjUBhGa105FSVFrFEoiEiCKBSGkUoZcydX88I2hYKIJIdCYQRzJ1cpFEQkURQKI5g7uZodew+xY29P3KWIiIwLhcIIdLFZRJJGoTCCk5qrAXSxWUQSQ6EwgsbqUiZVlbB6S3fcpYiIjAuFwihOatYdSCKSHAqFUZw0uYYXtu2hf8DjLkVEJOsUCqM4ubmag70DbFDneCKSACOGgpl9MGP8gqPmXZ+tonLJyS3RxWbdgSQiSTDakcLnMsb/71HzPjzGteSkOU3VmOkOJBFJhtFCwYYZH2q6IJWXFDG9voL2TvWWKiKFb7RQ8GHGh5ouWCc2VvGCjhREJAGKR5l/spmtIDoqOCGME6ZnZ7WyHHLK1Fr+8Px2Dhzqp7ykKO5yRESyZrRQmDcuVeS4+S01uMOard2cOb0u7nJERLJmxNNH7r4h8wPsBRYCk8J0IpwyJeoD6Tk92SwiBW60W1LvMbNTw3gLsIrorqN/N7PPjEN9OaG1rpzq0mJ1dyEiBW+0C82z3H1VGP8Q8IC7vx14HQm5JRXAzJg3pYZnNysURKSwjRYKvRnjFwP3Arj7HmAgW0XlotOm1vLc5m76+hP1s0UkYUYLhY1m9ikzeyfRtYT7AMysHJiQ7eJyycnN1fT0DbCha3/cpYiIZM1oofAR4BTgGuB97v5yaD8X+GEW68o581qii826riAihWzEW1LdfTvw8SHalwBLslVULpo7uZqSohQrO3bzttOnxF2OiEhWjBgKZnb3SPPd/R1jW07uKilOMa+lmmc6Xh59YRGRPDXaw2vnARuBnwJPkJD+joZzeutEfvX0JgYGnFQq0ZtCRArUaNcUmoHPA6cC3wAuBXa4+8Pu/nC2i8s1p7XWsrenj3U79G4FESlMoz3R3O/u97n71UQXl9uBh8zsU+NSXY45vbUWgFWbdsdciYhIdoz65jUzKzWzvwb+A7gOuBm461V871Yz225mqzLa6s3sATNbG4Z1od3M7GYzazezFWa28Nh/Uvac2FhFSXFK3V2ISMEarZuL24DHiJ5R+KK7n+3u/+Lum17Fun8EXHZU2w3Ag+4+B3gwTANcDswJn2uB77zqXzCOiotSzGupYYUuNotIgRrtSOFvgbnAp4HHzKw7fPaY2Yj/XHb3R4Cuo5oXA7eF8duAKzLaf+yRx4GJoa+lnHNGay0rO3YzMJCY10mISIKMdk0h5e7V4VOT8al295pj+HuT3X1LWPcWoCm0TyW6yymtI7S9gplda2ZtZtbW2dl5DCUcn1On1LLvUL+ebBaRgjTqNYVxMtT9nUP+U9zdb3H3Re6+qLGxMctlvVL6yebn1DmeiBSg8Q6FbenTQmG4PbR3ANMylmsFNo9zba/K3OYqJhQZK3UHkogUoPEOhbuBq8P41cCvM9qvCnchnQvsTp9myjWlxUXMa6lh+cZdcZciIjLmshYKZvZT4M/ASWbWYWYfAW4CLjWztUQPwt0UFr8XWEf0HMT3gE9mq66xsHB6Hc9s3K1utEWk4IzWzcUxc/f3DzPr4iGWdaJnIPLCmdMn8qPH1rNm6x5OnVobdzkiImMmVy4055WzZtQBsOwlnUISkcKiUDgGUyeW01RdyrINCgURKSwKhWNgZpw1o46lOlIQkQKjUDhGZ82oY2PXAbZ1H4y7FBGRMaNQOEbnzKoH4MkXj+7JQ0QkfykUjtH8lhoqS4p4ar1CQUQKh0LhGBUXpVgwfSJt63VdQUQKh0LhOCycXseard3s6+mLuxQRkTGhUDgOZ8+sZ8ChTbemikiBUCgch0Uz65hQZDzWviPuUkRExoRC4ThUlBRz1ow6HlmrUBCRwqBQOE5vmNvI6i3ddO7pibsUEZHjplA4Tq8/cRIAf163M+ZKRESOn0LhOJ0ypZbqsmJdVxCRgqBQOE5FKeN1sxp4XEcKIlIAFApj4JxZdazfuZ/t6gdJRPKcQmEMnH9CdF3hj7oLSUTynEJhDMxvqaGhsoSHXuiMuxQRkeOiUBgDqZTxxpObePj57Xpvs4jkNYXCGLn45Ca6D/apywsRyWsKhTFy4dxGSopS/P65bXGXIiJyzBQKY6SqtJjzTmjggdXbcPe4yxEROSYKhTF0yfzJbNi5n7907o27FBGRY6JQGEOXzGsC4IHntsdciYjIsVEojKGW2nLOaK3lt6u2xF2KiMgxUSiMsctPa2FFx246du2PuxQRkddMoTDGLj+1GYD7Vm2NuRIRkddOoTDGZjRUMq+lht8qFEQkDykUsuCtpzWzdMMuNnbpFJKI5BeFQhZcceZUAO5atinmSkREXhuFQha01lVw/gkN3LlsIwMDepBNRPKHQiFL3rOolY1dB3jixa64SxERedUUClly2SktVJUWc+fSjrhLERF51RQKWVJeUsTbTm/h3pVb2NvTF3c5IiKvikIhi96zqJUDvf3cu1JPOItIflAoZNHC6XXMnlTJnW06hSQi+SGWUDCz9Wa20syWm1lbaKs3swfMbG0Y1sVR21gyM951VitPru9i/Y59cZcjIjKqOI8U3ujuC9x9UZi+AXjQ3ecAD4bpvPeuha2kDF1wFpG8kEunjxYDt4Xx24ArYqxlzDTXlnHRSU3c8dRGevr64y5HRGREcYWCA78zs6Vmdm1om+zuWwDCsGmoL5rZtWbWZmZtnZ2d41Tu8bn6/Jns2NvDPc/ogrOI5La4QuECd18IXA5cZ2ZveLVfdPdb3H2Ruy9qbGzMXoVj6MITJ3FyczXfefgvesJZRHJaLKHg7pvDcDvwK+AcYJuZtQCEYcG8viyVMj5x0Qm0b9/L757bFnc5IiLDGvdQMLNKM6tOjwNvAlYBdwNXh8WuBn493rVl01tPa2FGQwXffqgddx0tiEhuiuNIYTLwqJk9AzwJ/Mbd7wNuAi41s7XApWG6YBQXpfj4fzuBFR27ebR9R9zliIgMqXi8/6C7rwPOGKJ9J3DxeNcznv564VS+/vsX+NaSdi6ckx/XQ0QkWXLpltSCV1pcxH+/cDaPr+vi8XU74y5HROQVFArj7AOvm8HkmlK+fN8aXVsQkZyjUBhn5SVFfPaSuTz90svc/6zuRBKR3KJQiMG7z2rlhMZKvnL/Gvr6B+IuR0RkkEIhBsVFKf7+zSfxl859/Fw9qIpIDlEoxOTNpzRz9sw6vnL/Grr2HYq7HBERQKEQGzPjX684jT0H+7jpt6vjLkdEBFAoxOqk5mo+cuEsft7WwZ/0QJuI5ACFQsw+e8lcZjdW8j9/8Qy79/fGXY6IJJxCIWZlE4r4+vsW0Lmnhy/+17NxlyMiCadQyAGnt07kkxedwF1Pb+I3K/TOBRGJj0IhR3zq4jksmDaRG+9aQceu/XGXIyIJpVDIEROKUtx85ZkMOHzmjuV6qE1EYqFQyCHTGyr40jtPpW3DLr7++7VxlyMiCTTuXWfLyBYvmMpj7Tv55pJ2WuvKufKc6XGXJCIJolDIQf9yxals7T7Ijb9aSUNVKZfOnxx3SSKSEDp9lINKilN894NncdrUWq6/fRlL1hTM66pFJMcpFHJUeUkRP7zmbE5squKjP27jF20b4y5JRBJAoZDDGqpK+dnHzuO82Q38/Z0r+NaSdr2YR0SySqGQ46pKi7n1mrNZvGAKX7n/ef757mfpH1AwiEh26EJzHigpTvG19y6gqbqU7/3xRTr39vDV9y6gbEJR3KWJSIFRKOSJVMr4wlvnM7mmjH/9zWq2dz/BN95/JlMnlsddmogUEJ0+yjMfvXA23/ybM1m9pZvLv/6I+koSkTGlUMhDbzt9Cvd++kJmN1Zx3e3L+PQdT7NLb28TkTGgUMhTMxoq+cXHz+Nzl87l3pVbuPRrj3DPis26O0lEjotCIY9NKErxdxfP4e7rX09LbRnX3/40H/j+EzyxbqfCQUSOieXz/3ksWrTI29ra4i4jJ/QPOP/x+Aa+8eBauvYd4ozWWj78+llcfmoLJcXKfhE5zMyWuvuiIecpFArLgUP9/HJZB7f+6UXWde5jck0pV503k/efM536ypK4yxORHKBQSKCBAefhtZ3c+uiL/HHtDkqLU7z9jCm8a2Err5tVTyplcZcoIjEZKRT0nEKBSqWMN57UxBtPamLttj388LH1/PrpTdy5tIOm6lIunjeZS+Y1ccGJk/QQnIgM0pFCguw/1MeDq7dz78ot/HHtDvb29FE2IcX5J0zighMnsXD6ROZPqaG0WCEhUsh0pCAAVJQU8/YzpvD2M6bQ09fP4+u6WLJmO0ue384fQvfcJcUp5rXUcPrUWk6bWsv8KTXMmVyloBBJCB0pCABbdx9k+cZdtK3fxarNu1m1qZu9PX0ATCgyTmyqZk5TFXOaqpg5qZLWunKm11dQX1mCma5PiOQTHSnIqJpry7istoXLTm0BogvV63fu47kt3aza1M2ard0se2kXdz+z+YjvlU8oorWunNa6cppry2msKqGusoRJVaXUV5ZQWz6B2vIJ1JRPoLq0WBe4RXKcQkGGlEoZsxurmN1YxdtOnzLYvv9QHy917aej6wAvde1n08sH6Ni1n45dB1jRsZuu/YcY7uAzZVBdNoHqsmIqS4qpLC2isjQarygtonxCERUlRZSXFFNSZEwoSkWf4tQR0yXFKUrS84oszM+YDsukzChOGSkzUikoSo+bkTKwjKGIRHIuFMzsMuAbQBHwfXe/KeaSJENFSTEnN9dwcnPNkPP7B5zdB3rZsbeHrn2H2H2gl90HeukOw90Hetnb08e+nj72H+pnb08f27oPsq+nn4O9/ew71MfB3oFx/lVgxuGwwAanhxtGBzw2+F0bHEbLpOdmBk5m9hxeZuj5Y/KbxnZ1Yx6eWYniBG3DK8+exkcvnD1m60vLqVAwsyLgW8ClQAfwlJnd7e7PxVuZvFpFKaO+suS4HpRzd/oGnN7+AQ71DXCof4Defqe3byBqS0/3D9Dbd9R0xncGBqL1DHh0OqzfffAFRQOh3QlDd9xhwA+3ux+5XHq+O/SHw6FoELV5WC7d7oO/53B71HDEYPA3j6WxvlI41pces3ElM9e34VivcFJV6diuMMipUADOAdrdfR2Amd0BLAYUCgliZoOngSr0ELbIuMq1TnGmAplvqO8IbYPM7FozazOzts7OznEtTkSk0OVaKAx1wu2Igy53v8XdF7n7osbGxnEqS0QkGXItFDqAaRnTrcDmYZYVEZExlmuh8BQwx8xmmVkJcCVwd8w1iYgkRk5daHb3PjO7Hrif6JbUW9392ZjLEhFJjJwKBQB3vxe4N+46RESSKNdOH4mISIwUCiIiMiive0k1s05gwzF+fRKwYwzLKTTaPiPT9hmdttHI4tw+M9x9yHv68zoUjoeZtQ3Xdaxo+4xG22d02kYjy9Xto9NHIiIySKEgIiKDkhwKt8RdQI7T9hmZts/otI1GlpPbJ7HXFERE5JWSfKQgIiJHUSiIiMigRIaCmV1mZs+bWbuZ3RB3PdlkZtPMbImZrTazZ83s06G93sweMLO1YVgX2s3Mbg7bZoWZLcxY19Vh+bVmdnVG+1lmtjJ852bLs5cem1mRmT1tZveE6Vlm9kT4nT8LnTNiZqVhuj3Mn5mxjhtD+/Nm9uaM9rzf18xsopndaWZrwn50nvafw8zss+G/rVVm9lMzK8vrfSh6DWFyPkQd7f0FmA2UAM8A8+OuK4u/twVYGMargReA+cD/Bm4I7TcAXw7jbwF+S/Rui3OBJ0J7PbAuDOvCeF2Y9yRwXvjOb4HL4/7dr3EbfQ64HbgnTP8cuDKMfxf4RBj/JPDdMH4l8LMwPj/sR6XArLB/FRXKvgbcBnw0jJcAE7X/DG6bqcCLQHnGvnNNPu9DSTxSGHzlp7sfAtKv/CxI7r7F3ZeF8T3AaqIdeTHRf+yE4RVhfDHwY488Dkw0sxbgzcAD7t7l7ruAB4DLwrwad/+zR3v3jzPWlfPMrBV4K/D9MG3AXwF3hkWO3jbpbXYncHFYfjFwh7v3uPuLQDvRfpb3+5qZ1QBvAH4A4O6H3P1ltP9kKgbKzawYqAC2kMf7UBJDYdRXfhaqcKh6JvAEMNndt0AUHEBTWGy47TNSe8cQ7fni68A/AANhugF42d37wnTm7xncBmH+7rD8a91m+WQ20An8MJxi+76ZVaL9BwB33wT8H+AlojDYDSwlj/ehJIbCqK/8LERmVgX8EviMu3ePtOgQbX4M7TnPzN4GbHf3pZnNQyzqo8wruG2ToRhYCHzH3c8E9hGdLhpOorZRuJaymOiUzxSgErh8iEXzZh9KYigk7pWfZjaBKBB+4u53heZt4dCdMNwe2ofbPiO1tw7Rng8uAN5hZuuJDsv/iujIYWI4FQBH/p7BbRDm1wJdvPZtlk86gA53fyJM30kUEtp/IpcAL7p7p7v3AncB55PH+1ASQyFRr/wM5yt/AKx2969mzLobSN8BcjXw64z2q8JdJOcCu8PpgfuBN5lZXfjX0ZuA+8O8PWZ2bvhbV2WsK6e5+43u3uruM4n2gz+4+weAJcC7w2JHb5v0Nnt3WN5D+5XhzpJZwByii6d5v6+5+1Zgo5mdFJouBp5D+0/aS8C5ZlYR6k9vn/zdh+K+eh/Hh+gOiReIrup/Ie56svxbX090uLkCWB4+byE6j/kgsDYM68PyBnwrbJuVwKKMdX2Y6AJYO/ChjPZFwKrwnW8SnpTPpw9wEYfvPppN9B9kO/ALoDS0l4Xp9jB/dsb3vxB+//Nk3D1TCPsasABoC/vQfxLdPaT953D9XwTWhN/w70R3EOXtPqRuLkREZFASTx+JiMgwFAoiIjJIoSAiIoMUCiIiMkihICIigxQKknfMbG8YzjSzvxnjdX/+qOnHxnL9Y83MrjGzb8ZdhxQOhYLks5nAawoFMysaZZEjQsHdz3+NNeWVV7E9JGEUCpLPbgIuNLPloU/7IjP7ipk9Ffry/xiAmV1k0Tslbid6oAoz+08zWxr6wb82tN1E1NvlcjP7SWhLH5VYWPcqi/r+f1/Guh+yw+8b+El4svUIYZkvm9mTZvaCmV0Y2o/4l76Z3WNmF6X/dvjOUjP7vZmdE9azzszekbH6aWZ2n0V97v9Txro+GP7ecjP7f+kACOv9X2b2BFGX1SKHxf00oD76vNYPsDcMLyI8hRymrwX+MYyXEj2FOysstw+YlbFs+gnccqInURsy1z3E33oXUXfPRcBkou4NWsK6dxP1SZMC/gy8foiaHwL+LYy/Bfh9GL8G+GbGcvcAF4VxJzzZCvwK+B0wATgDWJ7x/S1ETxinf8siYB7wX8CEsNy3gasy1vveuP931Cc3P+kOm0QKwZuA080s3edMLVEfMoeAJz3qpz7t78zsnWF8Wlhu5wjrfj3wU3fvJ+oM7mHgbKA7rLsDwMyWE53WenSIdaQ7I1walhnNIeC+ML4S6HH3XjNbedT3H3D3neHv3xVq7QPOAp4KBy7lHO60rp+og0SRV1AoSCEx4FPufv8RjdHpmH1HTV8CnOfu+83sIaI+aUZb93B6Msb7Gf6/q54hlunjyNO4mXX0unu6H5qB9PfdfSCjB054ZVfK6S6Xb3P3G4eo42AIN5FX0DUFyWd7iF4xmnY/8AmLugrHzOZa9EKYo9UCu0IgnEz02si03vT3j/II8L5w3aKR6G1kT47Bb1gPLDCzlJlNI3rT1mt1qUXvTC4nesPXn4g6qXu3mTXB4Du5Z4xBvVLgdKQg+WwF0GdmzwA/Ar5BdFplWbjY28nQr3a8D/i4ma0g6pHy8Yx5twArzGyZR91op/2K6KLsM0T/Ev8Hd98aQuV4/InoHb8ria4HLDuGdTxK1DvnicDt7t4GYGb/CPzOzFJAL3AdsOE465UCp15SRURkkE4fiYjIIIWCiIgMUiiIiMgghYKIiAxSKIiIyCCFgoiIDFIoiIjIoP8PQVPVMNi4uj8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%pylab inline\n",
    "plot(range(len(stoch_errors_by_iter)), stoch_errors_by_iter)\n",
    "xlabel('Iteration number')\n",
    "ylabel('MSE')\n",
    "\n",
    "# the algorithm converges"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 3.91069256e+00],\n",
       "       [ 2.78209808e+00],\n",
       "       [-8.10462217e-03],\n",
       "       [ 1.40190566e+01]])"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stoch_grad_desc_weights"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "MSE on the last iteration:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.7844125884067035"
      ]
     },
     "execution_count": 108,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stoch_errors_by_iter[-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}