HW/AM207_HW8.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Homework #8: Due November 4th, 2021, 11:59pm\n",
    "## Bayesian Neural Networks\n",
    "\n",
    "**AM 207: Advanced Scientific Computing**<br>\n",
    "**Instructor: Weiwei Pan**<br>\n",
    "**Fall 2021**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Name:**\n",
    "\n",
    "**Students collaborators:**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Instructions:\n",
    "\n",
    "**Submission Format:** Use this notebook as a template to complete your homework. Please intersperse text blocks (using Markdown cells) amongst `python` code and results -- format your submission for maximum readability. Your assignments will be graded for correctness as well as clarity of exposition and presentation -- a “right” answer by itself without an explanation or is presented with a difficult to follow format will receive no credit.\n",
    "\n",
    "**Code Check:** Before submitting, you must do a \"Restart and Run All\" under \"Kernel\" in the Jupyter or colab menu. Portions of your submission that contains syntactic or run-time errors will not be graded.\n",
    "\n",
    "**Libraries and packages:** Unless a problems specifically asks you to implement from scratch, you are welcomed to use any `python` library package in the standard Anaconda distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from autograd import numpy as np\n",
    "from autograd import grad\n",
    "from autograd.misc.optimizers import adam, sgd\n",
    "from autograd import scipy as sp\n",
    "import pandas as pd\n",
    "import numpy\n",
    "import matplotlib.pyplot as plt\n",
    "import sys\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Feedforward:\n",
    "    def __init__(self, architecture, random=None, weights=None):\n",
    "        self.params = {'H': architecture['width'],\n",
    "                       'L': architecture['hidden_layers'],\n",
    "                       'D_in': architecture['input_dim'],\n",
    "                       'D_out': architecture['output_dim'],\n",
    "                       'activation_type': architecture['activation_fn_type'],\n",
    "                       'activation_params': architecture['activation_fn_params']}\n",
    "\n",
    "        self.D = (  (architecture['input_dim'] * architecture['width'] + architecture['width'])\n",
    "                  + (architecture['output_dim'] * architecture['width'] + architecture['output_dim'])\n",
    "                  + (architecture['hidden_layers'] - 1) * (architecture['width']**2 + architecture['width'])\n",
    "                 )\n",
    "\n",
    "        if random is not None:\n",
    "            self.random = random\n",
    "        else:\n",
    "            self.random = np.random.RandomState(0)\n",
    "\n",
    "        self.h = architecture['activation_fn']\n",
    "\n",
    "        if weights is None:\n",
    "            self.weights = self.random.normal(0, 1, size=(1, self.D))\n",
    "        else:\n",
    "            self.weights = weights\n",
    "\n",
    "        self.objective_trace = np.empty((1, 1))\n",
    "        self.weight_trace = np.empty((1, self.D))\n",
    "\n",
    "\n",
    "    def forward(self, weights, x):\n",
    "        ''' Forward pass given weights and input '''\n",
    "        H = self.params['H']\n",
    "        D_in = self.params['D_in']\n",
    "        D_out = self.params['D_out']\n",
    "\n",
    "        assert weights.shape[1] == self.D\n",
    "\n",
    "        if len(x.shape) == 2:\n",
    "            assert x.shape[0] == D_in\n",
    "            x = x.reshape((1, D_in, -1))\n",
    "        else:\n",
    "            assert x.shape[1] == D_in\n",
    "\n",
    "        weights = weights.T\n",
    "\n",
    "\n",
    "        #input to first hidden layer\n",
    "        W = weights[:H * D_in].T.reshape((-1, H, D_in))\n",
    "        b = weights[H * D_in:H * D_in + H].T.reshape((-1, H, 1))\n",
    "        input = self.h(np.matmul(W, x) + b)\n",
    "        index = H * D_in + H\n",
    "\n",
    "        assert input.shape[1] == H\n",
    "\n",
    "        #additional hidden layers\n",
    "        for _ in range(self.params['L'] - 1):\n",
    "            before = index\n",
    "            W = weights[index:index + H * H].T.reshape((-1, H, H))\n",
    "            index += H * H\n",
    "            b = weights[index:index + H].T.reshape((-1, H, 1))\n",
    "            index += H\n",
    "            output = np.matmul(W, input) + b\n",
    "            input = self.h(output)\n",
    "\n",
    "            assert input.shape[1] == H\n",
    "\n",
    "        #output layer\n",
    "        W = weights[index:index + H * D_out].T.reshape((-1, D_out, H))\n",
    "        b = weights[index + H * D_out:].T.reshape((-1, D_out, 1))\n",
    "        output = np.matmul(W, input) + b\n",
    "        assert output.shape[1] == self.params['D_out']\n",
    "\n",
    "        return output\n",
    "\n",
    "    def make_objective(self, x_train, y_train, reg_param):\n",
    "\n",
    "        def objective(W, t):\n",
    "            squared_error = np.linalg.norm(y_train - self.forward(W, x_train), axis=1)**2\n",
    "            if reg_param is None:\n",
    "                sum_error = np.sum(squared_error)\n",
    "                return sum_error\n",
    "            else:\n",
    "                mean_error = np.mean(squared_error) + reg_param * np.linalg.norm(W)\n",
    "                return mean_error\n",
    "\n",
    "        return objective, grad(objective)\n",
    "\n",
    "    def fit(self, x_train, y_train, params, reg_param=None):\n",
    "\n",
    "        assert x_train.shape[0] == self.params['D_in']\n",
    "        assert y_train.shape[0] == self.params['D_out']\n",
    "\n",
    "        ### make objective function for training\n",
    "        self.objective, self.gradient = self.make_objective(x_train, y_train, reg_param)\n",
    "\n",
    "        ### set up optimization\n",
    "        step_size = 0.01\n",
    "        max_iteration = 5000\n",
    "        check_point = 100\n",
    "        weights_init = self.weights.reshape((1, -1))\n",
    "        mass = None\n",
    "        optimizer = 'adam'\n",
    "        random_restarts = 5\n",
    "\n",
    "        if 'step_size' in params.keys():\n",
    "            step_size = params['step_size']\n",
    "        if 'max_iteration' in params.keys():\n",
    "            max_iteration = params['max_iteration']\n",
    "        if 'check_point' in params.keys():\n",
    "            self.check_point = params['check_point']\n",
    "        if 'init' in params.keys():\n",
    "            weights_init = params['init']\n",
    "        if 'call_back' in params.keys():\n",
    "            call_back = params['call_back']\n",
    "        if 'mass' in params.keys():\n",
    "            mass = params['mass']\n",
    "        if 'optimizer' in params.keys():\n",
    "            optimizer = params['optimizer']\n",
    "        if 'random_restarts' in params.keys():\n",
    "            random_restarts = params['random_restarts']\n",
    "\n",
    "        def call_back(weights, iteration, g):\n",
    "            ''' Actions per optimization step '''\n",
    "            objective = self.objective(weights, iteration)\n",
    "            self.objective_trace = np.vstack((self.objective_trace, objective))\n",
    "            self.weight_trace = np.vstack((self.weight_trace, weights))\n",
    "            if iteration % check_point == 0:\n",
    "                print(\"Iteration {} lower bound {}; gradient mag: {}\".format(iteration, objective, np.linalg.norm(self.gradient(weights, iteration))))\n",
    "\n",
    "        ### train with random restarts\n",
    "        optimal_obj = 1e16\n",
    "        optimal_weights = self.weights\n",
    "\n",
    "        for i in range(random_restarts):\n",
    "            if optimizer == 'adam':\n",
    "                adam(self.gradient, weights_init, step_size=step_size, num_iters=max_iteration, callback=call_back)\n",
    "            local_opt = np.min(self.objective_trace[-100:])\n",
    "            if local_opt < optimal_obj:\n",
    "                opt_index = np.argmin(self.objective_trace[-100:])\n",
    "                self.weights = self.weight_trace[-100:][opt_index].reshape((1, -1))\n",
    "            weights_init = self.random.normal(0, 1, size=(1, self.D))\n",
    "\n",
    "        self.objective_trace = self.objective_trace[1:]\n",
    "        self.weight_trace = self.weight_trace[1:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem Description: Bayesian Neural Network Regression\n",
    "In this homework you will explore Bayeisan models for regression with neural networks -- essentially, you will place priors on the network parameters $\\mathbf{W}$ and infer the posterior $p(\\mathbf{W}|\\text{Data})$. These Bayesian models are called ***Bayesian neural networks***. The data for this regression problem is in `HW8_data.csv`. You are provided an implementation of a neural network as well as an example of how to use it -- but feel free to implement your own version!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "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>x</th>\n",
       "      <th>y</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>-6.0</td>\n",
       "      <td>-3.380284</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>-5.6</td>\n",
       "      <td>-2.892117</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>-5.2</td>\n",
       "      <td>-2.690059</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>-4.8</td>\n",
       "      <td>-2.040000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>-4.4</td>\n",
       "      <td>-1.399942</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     x         y\n",
       "0 -6.0 -3.380284\n",
       "1 -5.6 -2.892117\n",
       "2 -5.2 -2.690059\n",
       "3 -4.8 -2.040000\n",
       "4 -4.4 -1.399942"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#read in the data\n",
    "data = pd.read_csv('HW8_data.csv')\n",
    "x_train = data['x'].values.reshape((1, -1))\n",
    "y_train = data['y'].values.reshape((1, -1))\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's instantiate a neural network with 1-hidden layer, 5-hidden nodes and `relu` activation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "###relu activation\n",
    "activation_fn_type = 'relu'\n",
    "activation_fn = lambda x: np.maximum(np.zeros(x.shape), x)\n",
    "\n",
    "\n",
    "###neural network model design choices\n",
    "width = 5\n",
    "hidden_layers = 1\n",
    "input_dim = 1\n",
    "output_dim = 1\n",
    "\n",
    "architecture = {'width': width,\n",
    "               'hidden_layers': hidden_layers,\n",
    "               'input_dim': input_dim,\n",
    "               'output_dim': output_dim,\n",
    "               'activation_fn_type': 'relu',\n",
    "               'activation_fn_params': 'rate=1',\n",
    "               'activation_fn': activation_fn}\n",
    "\n",
    "#set random state to make the experiments replicable\n",
    "rand_state = 0\n",
    "random = np.random.RandomState(rand_state)\n",
    "\n",
    "#instantiate a Feedforward neural network object\n",
    "nn = Feedforward(architecture, random=random)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fit the neural network to the data and visualize the fit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration 0 lower bound 1188.2522119587577; gradient mag: 3294.9472005145494\n",
      "Iteration 100 lower bound 528.2005162282458; gradient mag: 2014.0223481861012\n",
      "Iteration 200 lower bound 231.71928757081068; gradient mag: 1226.9831178661013\n",
      "Iteration 300 lower bound 103.99541756891784; gradient mag: 716.7107003654073\n",
      "Iteration 400 lower bound 55.24716113156889; gradient mag: 376.23346335998286\n",
      "Iteration 500 lower bound 39.94989749994618; gradient mag: 177.9641773056736\n",
      "Iteration 600 lower bound 35.60307206137053; gradient mag: 77.19358484414678\n",
      "Iteration 700 lower bound 34.025856107542445; gradient mag: 38.67687513287844\n",
      "Iteration 800 lower bound 32.961429361488406; gradient mag: 31.129490339137103\n",
      "Iteration 900 lower bound 31.989641942996613; gradient mag: 24.752496721932584\n",
      "Iteration 1000 lower bound 31.082564536611727; gradient mag: 24.616299563320247\n",
      "Iteration 1100 lower bound 30.155032233758106; gradient mag: 24.248289796695595\n",
      "Iteration 1200 lower bound 29.210013809597285; gradient mag: 23.699589428153818\n",
      "Iteration 1300 lower bound 28.253006614697387; gradient mag: 23.095433479427733\n",
      "Iteration 1400 lower bound 27.289273636922445; gradient mag: 22.47531773786042\n",
      "Iteration 1500 lower bound 26.323662523356894; gradient mag: 21.850371675396286\n",
      "Iteration 1600 lower bound 25.36059928299808; gradient mag: 21.2251268262558\n",
      "Iteration 1700 lower bound 24.404090286894164; gradient mag: 20.602713937891433\n",
      "Iteration 1800 lower bound 23.45772579649591; gradient mag: 19.985836957824542\n",
      "Iteration 1900 lower bound 22.524686048678994; gradient mag: 19.376902849934968\n",
      "Iteration 2000 lower bound 21.60775070123689; gradient mag: 18.77801384940906\n",
      "Iteration 2100 lower bound 20.709312052421566; gradient mag: 18.190943350296987\n",
      "Iteration 2200 lower bound 19.831392143669202; gradient mag: 17.617113686328434\n",
      "Iteration 2300 lower bound 18.9756636132038; gradient mag: 17.057579413745877\n",
      "Iteration 2400 lower bound 18.14347397435544; gradient mag: 16.513017921438422\n",
      "Iteration 2500 lower bound 17.335872838126047; gradient mag: 15.983728786591483\n",
      "Iteration 2600 lower bound 16.553641480811155; gradient mag: 15.469642874398222\n",
      "Iteration 2700 lower bound 15.79979845292681; gradient mag: 14.634100216702938\n",
      "Iteration 2800 lower bound 15.089224741859162; gradient mag: 13.79947750065518\n",
      "Iteration 2900 lower bound 14.412878302812286; gradient mag: 13.252437242012403\n",
      "Iteration 3000 lower bound 13.766721311122604; gradient mag: 12.795482048078641\n",
      "Iteration 3100 lower bound 13.147000834095968; gradient mag: 12.399165758017654\n",
      "Iteration 3200 lower bound 12.550715085203827; gradient mag: 12.043477234157677\n",
      "Iteration 3300 lower bound 11.98798870979473; gradient mag: 11.671698244924366\n",
      "Iteration 3400 lower bound 11.470331425797706; gradient mag: 11.552845486196434\n",
      "Iteration 3500 lower bound 10.989896461078454; gradient mag: 11.570145170353355\n",
      "Iteration 3600 lower bound 10.545865383076633; gradient mag: 9.157301353805478\n",
      "Iteration 3700 lower bound 10.13439033091308; gradient mag: 8.739366437372619\n",
      "Iteration 3800 lower bound 9.744699797627586; gradient mag: 8.48450205329216\n",
      "Iteration 3900 lower bound 9.36943762284221; gradient mag: 8.279364759162569\n",
      "Iteration 4000 lower bound 9.023880352484412; gradient mag: 12.814088131372916\n",
      "Iteration 4100 lower bound 8.714826491071893; gradient mag: 8.636545905932076\n",
      "Iteration 4200 lower bound 8.436787581496315; gradient mag: 9.251058145534492\n",
      "Iteration 4300 lower bound 8.187373568381847; gradient mag: 7.442042690088629\n",
      "Iteration 4400 lower bound 7.964500968027264; gradient mag: 6.05345432191102\n",
      "Iteration 4500 lower bound 7.7662669301874585; gradient mag: 5.051211535334474\n",
      "Iteration 4600 lower bound 7.589279565032134; gradient mag: 4.59773932214737\n",
      "Iteration 4700 lower bound 7.4224106345582115; gradient mag: 4.45650745335735\n",
      "Iteration 4800 lower bound 7.258760999455524; gradient mag: 4.341240693065324\n",
      "Iteration 4900 lower bound 7.106644631620015; gradient mag: 4.521272792643232\n",
      "Iteration 5000 lower bound 6.97407040094256; gradient mag: 5.201564339555988\n",
      "Iteration 5100 lower bound 6.856972561289584; gradient mag: 6.1193316776037445\n",
      "Iteration 5200 lower bound 6.752919981452237; gradient mag: 7.133712145135448\n",
      "Iteration 5300 lower bound 6.660138579044939; gradient mag: 8.182158552640534\n",
      "Iteration 5400 lower bound 6.575879541247726; gradient mag: 9.166712555050204\n",
      "Iteration 5500 lower bound 6.498596062906512; gradient mag: 3.0663031457822707\n",
      "Iteration 5600 lower bound 6.426834605421498; gradient mag: 2.4331381135547065\n",
      "Iteration 5700 lower bound 6.357433884941598; gradient mag: 2.2603654201225507\n",
      "Iteration 5800 lower bound 6.286273937604813; gradient mag: 2.239535151899242\n",
      "Iteration 5900 lower bound 6.213043834805276; gradient mag: 2.225849214024523\n",
      "Iteration 6000 lower bound 6.137748546533556; gradient mag: 2.2122700356057394\n",
      "Iteration 6100 lower bound 6.060402578608479; gradient mag: 2.1983511683951003\n",
      "Iteration 6200 lower bound 5.982027488105805; gradient mag: 2.133700892111169\n",
      "Iteration 6300 lower bound 5.905113988567324; gradient mag: 2.4321267978132246\n",
      "Iteration 6400 lower bound 5.828319129576141; gradient mag: 2.827227841216459\n",
      "Iteration 6500 lower bound 5.7518306395462595; gradient mag: 3.396839218046507\n",
      "Iteration 6600 lower bound 5.672819183849378; gradient mag: 3.511172106456379\n",
      "Iteration 6700 lower bound 5.59331943219885; gradient mag: 3.936185498220977\n",
      "Iteration 6800 lower bound 5.511770386421798; gradient mag: 4.246984794757417\n",
      "Iteration 6900 lower bound 5.428607645781736; gradient mag: 14.014120586219285\n",
      "Iteration 7000 lower bound 5.344999645585854; gradient mag: 4.847156007867547\n",
      "Iteration 7100 lower bound 5.2587731740380494; gradient mag: 5.037990425758324\n",
      "Iteration 7200 lower bound 5.170245548651128; gradient mag: 4.91389612736596\n",
      "Iteration 7300 lower bound 5.081990798375702; gradient mag: 5.148618362507468\n",
      "Iteration 7400 lower bound 4.992323968370187; gradient mag: 5.062612851660892\n",
      "Iteration 7500 lower bound 4.90219687955313; gradient mag: 5.0782333056152655\n",
      "Iteration 7600 lower bound 4.811691076733211; gradient mag: 12.726071077503757\n",
      "Iteration 7700 lower bound 4.722034252465464; gradient mag: 5.627923598935833\n",
      "Iteration 7800 lower bound 4.63046615524087; gradient mag: 5.3965319790773085\n",
      "Iteration 7900 lower bound 4.5393565041325035; gradient mag: 5.2116672953127425\n",
      "Iteration 8000 lower bound 4.4509554693924125; gradient mag: 5.442647223210916\n",
      "Iteration 8100 lower bound 4.3600620974975275; gradient mag: 5.203891398270261\n",
      "Iteration 8200 lower bound 4.271237398872255; gradient mag: 5.050050070322356\n",
      "Iteration 8300 lower bound 4.184440979576098; gradient mag: 5.24679258938914\n",
      "Iteration 8400 lower bound 4.098659993997213; gradient mag: 5.578840619974285\n",
      "Iteration 8500 lower bound 4.012291494325342; gradient mag: 5.381437473766087\n",
      "Iteration 8600 lower bound 3.928451645703673; gradient mag: 5.3461552886190535\n",
      "Iteration 8700 lower bound 3.8439395365195805; gradient mag: 9.992393975865808\n",
      "Iteration 8800 lower bound 3.7623541691665165; gradient mag: 9.79791476897376\n",
      "Iteration 8900 lower bound 3.680935626235934; gradient mag: 9.233951145662296\n",
      "Iteration 9000 lower bound 3.602698412892009; gradient mag: 9.59071358903378\n",
      "Iteration 9100 lower bound 3.522741405374152; gradient mag: 4.8315358324131195\n",
      "Iteration 9200 lower bound 3.4454461281570974; gradient mag: 4.546978550730063\n",
      "Iteration 9300 lower bound 3.368233328231543; gradient mag: 8.077544206848126\n",
      "Iteration 9400 lower bound 3.2931617717908903; gradient mag: 5.198762473040688\n",
      "Iteration 9500 lower bound 3.2158141840415206; gradient mag: 3.9061861443065014\n",
      "Iteration 9600 lower bound 3.1406962333053876; gradient mag: 3.8992430901073285\n",
      "Iteration 9700 lower bound 3.06586639695715; gradient mag: 6.9120899421115265\n",
      "Iteration 9800 lower bound 2.9908245411967624; gradient mag: 3.986129059418889\n",
      "Iteration 9900 lower bound 2.9166903463356233; gradient mag: 6.789802894150458\n",
      "Iteration 10000 lower bound 2.8425619783889498; gradient mag: 3.277143244727339\n",
      "Iteration 10100 lower bound 2.768815835197241; gradient mag: 2.6436406206319267\n",
      "Iteration 10200 lower bound 2.6967620063896365; gradient mag: 1.7554182380346504\n",
      "Iteration 10300 lower bound 2.626726890192887; gradient mag: 1.9079914587407558\n",
      "Iteration 10400 lower bound 2.558830428532131; gradient mag: 1.0932264174024424\n",
      "Iteration 10500 lower bound 2.494168018543738; gradient mag: 0.773058738036761\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration 10600 lower bound 2.433316490178128; gradient mag: 0.767066193236975\n",
      "Iteration 10700 lower bound 2.376555349336703; gradient mag: 0.7553579421583901\n",
      "Iteration 10800 lower bound 2.3238459422563476; gradient mag: 0.7228024683428217\n",
      "Iteration 10900 lower bound 2.275330834330171; gradient mag: 0.6845211470546975\n",
      "Iteration 11000 lower bound 2.231249650994422; gradient mag: 0.6413826587052399\n",
      "Iteration 11100 lower bound 2.191827804434735; gradient mag: 0.5943225261135897\n",
      "Iteration 11200 lower bound 2.1572058127250573; gradient mag: 0.5444130874954253\n",
      "Iteration 11300 lower bound 2.1273989093957457; gradient mag: 0.49279021686790947\n",
      "Iteration 11400 lower bound 2.1022813532268936; gradient mag: 0.44060151668110925\n",
      "Iteration 11500 lower bound 2.0815910694099786; gradient mag: 0.38895956298850276\n",
      "Iteration 11600 lower bound 2.0649500999980592; gradient mag: 0.33889610830749123\n",
      "Iteration 11700 lower bound 2.0518957709682137; gradient mag: 0.2913192491317288\n",
      "Iteration 11800 lower bound 2.0419171399915625; gradient mag: 0.24697771788249834\n",
      "Iteration 11900 lower bound 2.0344915040292233; gradient mag: 0.2064362191304809\n",
      "Iteration 12000 lower bound 2.0291165873537573; gradient mag: 0.17006423386258412\n",
      "Iteration 12100 lower bound 2.025335369827159; gradient mag: 0.13803877718621319\n",
      "Iteration 12200 lower bound 2.022752096644936; gradient mag: 0.11035976785043246\n",
      "Iteration 12300 lower bound 2.021039532301389; gradient mag: 0.0868752864632504\n",
      "Iteration 12400 lower bound 2.0199387199884886; gradient mag: 0.06731321537538267\n",
      "Iteration 12500 lower bound 2.019253225668887; gradient mag: 0.051315577914230893\n",
      "Iteration 12600 lower bound 2.0188400594620663; gradient mag: 0.038472237005887755\n",
      "Iteration 12700 lower bound 2.018599266629262; gradient mag: 0.028351317068662885\n",
      "Iteration 12800 lower bound 2.0184637229425677; gradient mag: 0.020524595354938835\n",
      "Iteration 12900 lower bound 2.0183901201436703; gradient mag: 0.01458699820847295\n",
      "Iteration 13000 lower bound 2.0183516178077436; gradient mag: 0.010170104881073427\n",
      "Iteration 13100 lower bound 2.0183322452743586; gradient mag: 0.006950134375104933\n",
      "Iteration 13200 lower bound 2.0183228858875357; gradient mag: 0.004651251256666468\n",
      "Iteration 13300 lower bound 2.018318552295065; gradient mag: 0.0030451959053903956\n",
      "Iteration 13400 lower bound 2.018316633241322; gradient mag: 0.0019482655846468522\n",
      "Iteration 13500 lower bound 2.0183158223137134; gradient mag: 0.0012165923762785585\n",
      "Iteration 13600 lower bound 2.0183154961261143; gradient mag: 0.0007405247328018651\n",
      "Iteration 13700 lower bound 2.0183153715620183; gradient mag: 0.00043875394233396255\n",
      "Iteration 13800 lower bound 2.0183153265300566; gradient mag: 0.00025265818478568795\n",
      "Iteration 13900 lower bound 2.0183153111653813; gradient mag: 0.00014117995446784006\n",
      "Iteration 14000 lower bound 2.018315306233877; gradient mag: 7.641641274522837e-05\n",
      "Iteration 14100 lower bound 2.018315304750106; gradient mag: 3.9991531341181805e-05\n",
      "Iteration 14200 lower bound 2.0183153043331865; gradient mag: 2.019550449158989e-05\n",
      "Iteration 14300 lower bound 2.0183153042242075; gradient mag: 9.820275320292376e-06\n",
      "Iteration 14400 lower bound 2.018315304197825; gradient mag: 4.587610150174298e-06\n",
      "Iteration 14500 lower bound 2.018315304191941; gradient mag: 2.053921182451195e-06\n",
      "Iteration 14600 lower bound 2.018315304190728; gradient mag: 8.789791752595096e-07\n",
      "Iteration 14700 lower bound 2.018315304190497; gradient mag: 3.585832781695699e-07\n",
      "Iteration 14800 lower bound 2.018316152372446; gradient mag: 0.12709802814595308\n",
      "Iteration 14900 lower bound 2.018315304250162; gradient mag: 0.001028906285981916\n"
     ]
    }
   ],
   "source": [
    "###define design choices in gradient descent\n",
    "params = {'step_size':1e-3, \n",
    "          'max_iteration':15000, \n",
    "          'random_restarts':1}\n",
    "\n",
    "#fit my neural network to minimize MSE on the given data\n",
    "nn.fit(x_train, y_train, params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3de5zOdf7/8cfbcQwpSZtNrtERMYNGKEll6YAKLTUpOUypbTvXajqpnd0tdivlW7/pUpFJQkKl6CA2hwwGhVKZcUgZ5DgNg9fvj/fM7AzDzLg+1/W+rmte99vtunEd5vN5zuk1n+t9NCKCUkqpyFXFdQCllFKB0UKulFIRTgu5UkpFOC3kSikV4bSQK6VUhKvm4qSnnHKKxMXFuTi1UkpFrCVLlmwVkQaHP+6kkMfFxZGRkeHi1EopFbGMMdmlPa5NK0opFeG0kCulVIQLuJAbY2KMMV8bY5YbY741xgz3IphSSqny8aKNfB9wuYjsMcZUB/5rjJkpIgs9OLZSSqkyBFzIxS7WsqfgbvWCmy7gopRSIeJJG7kxpqoxJhPYAswWkUWlvCbZGJNhjMnIycnx4rRKKaXwqJCLyEERaQU0Ai40xrQo5TVpIpIoIokNGhwxDFIppdRx8nQcuYjsMMbMAa4EvvHy2KoS2b0bNm8ueatTB265BWrWdJ1OqbATcCE3xjQA8guKeC2gC/BswMlU5fPss/DMM7B3b+nPjxgBL74IV10V2lxKhTkvrsgbAmONMVWxTTXvisgHHhxXVSZbt8JTT0FiIvToAQ0blrwtWgT33ANXXw3XXQdvvw21arlOrVRY8GLUygqgtQdZVGX2//4f5OXZf5s3P/L5q66CK66Af/8bUlKgf394912oonPalNLfAuXe/v0wejR07Vp6ES9UowYMG2abWKZMgYcfDl1GpcKYk0WzlCrh3Xdth+brr5fv9fffD+vW2avzuDj4y1+CGk+pcKeFXLklAs8/D02b2ivy8jDGdnquX2/bzX0+266uVCWlTSvKrf/+F5YuhXvvrVh7d9WqMGECtGkD/fqBLousKjEt5MqtF1+Ek0+2nZcVVbs2zJgBDRpA9+6QleV5PKUigRZy5c4vv8C0aXDbbRAbe3zHOO00mDkT9u2zQxN/+83bjEpFAC3kyp1x4+DAARg8OLDjNGsG778PP/wAvXrZoq5UJaKFXLkhAn4/dOxoOzoDdeml8MYbMGeO/cMgugCnqjx01IpyY948WLvWTu7xSlKSbSd/7DE7LPGZZ7w7tlJhTAu5csPvh7p1oU8fb4/76KO2mP/979CkCQwc6O3xlQqEiB0+6zFtWlGht2MHTJoEN91kR554yRj4v/+Dbt0gORlmzfL2+Eodr7lz7VpCP/7o+aG1kKvQmzDBrqsyZEhwjl+9up0tev759op/+fLgnEep8ti+3fbbXHopbNsGW7Z4fgot5Cr0/H5o1cpO5gmWunXhww/tv9dcAxs3Bu9cSpVGBMaPt535b74JDz0E334LHTp4fiot5Cq0li61t2BdjRfXqBF89BHs2mWL+a5dwT+nUmCHwnbtaie6nXkmLFkCzz3nfVNiAS3kKrTGjIGYGNs+Hgrx8TB5sr0SuuEGyM8PzXlV5bR/P6SmQosW8PXXdlXPr76ChISgnlYLuQqd3FxIT7cF9aSTQnferl0hLc12fA4dqmPMVXB89RW0bm2Hv/boAatXw5132nWBgkwLuQqdKVNg587AZ3Iej4ED7S/YmDHwj3+E/vwqev32G9x+u53ctncvfPCBHZX1xz+GLIKOI1eh4/fDOefAJZe4Of/TT5ecMJSU5CaHig4i8M47duXObdvggQdg+PCgtYMfixZyFRrff2/H0T77bFAmRJSLMfaKfONGu1DX6adD585usqjI9uOPttlk1ixo2xY++cSOxHJEm1ZUaIwZA9WqwS23uM1Rowa89x6cfTZcfz2sWuU2j4os+fnwz3/azswFC2DUKPuvwyIOWshVKOTn23G03bvbZWddq1fPDkusWdMuffvLL64TqUgwf76d+/Doo/bnZvVquPvukHRmlkULuQq+Dz6ws9lCMXa8vOLi7IShnBz7B2bvXteJVLj67Te44w64+GLbWT9tmu24P/1018mKBFzIjTFnGGO+MMasNsZ8a4y5x4tgKor4/faHvls310lKuuACmDgRli2DG2+EgwddJ1LhRMT+fDRrBq+9Zjs1V62Cnj1dJzuCF1fkB4AHRKQZ0B64yxjT3IPjqmiwYQN8/LEd/hcGb0GP0L07vPSS3TLunnt0jLmy1q2zzSf9+tkZwosX203C69RxnaxUAY9aEZHNwOaC/+82xqwGTge0F0nZtnGR8F5O9s477S/uyJF26dsHHnCdSLmSn28L9lNP2QuPF16Av/wlPC9CivF0+KExJg5oDSwq5blkIBmgcePGXp5WhatDh+xolS5dbJt0OHv2WTvG/MEHwefzfp10Ff4WLbJLH69YAdddZ9+pNWrkOlW5eNbZaYypA0wB7hWRI1YnEpE0EUkUkcQGDRp4dVoVzj77DLKz3czkrKgqVeweoh06wM032xEKqnLYuRPuust+77dvh6lT7S1Cijh4VMiNMdWxRTxdRN7z4pgqCvj9UL8+XHut6yTlU6sWTJ8OZ5xhO7TWrnWdSAWTiF1QrVkzePVV+OtfbWfmdde5TlZhXoxaMcAYYLWI/CfwSCoqbN1qr2r697fjtSPFKafAzJn2/1dfbT8PVWHp6enExcVRpUoV4uLiSE9Pr9DzQZeVZRe2uuEGO7dh0SLbHn7CCaHN4RURCegGdAQEWAFkFtyuPtbHXHDBBaKi3H/+IwIi33zjOsnx+eorkZo1RTp0EMnNdZ0moowfP15iY2OloC4IILGxsTJ+/PhyPR9U+fkiI0aIxMaK1K5tf07z84N/Xo8AGVJaHS7twWDftJBHuUOHRJo3t0Uwkk2aJGKMSPfuInv3uk4TMXw+X4kiXXjz+Xzlej5oFi0SSUiwZa9HD5Hs7OCeLwiOVsh1Zqfy3sKFtq0xEjo5j6VPH7uR84cf2pE327a5ThQR1q9ff8zHy3rec7t22fbv9u3tTN4pU+zszCgaPaeFXHnP77cTJ/78Z9dJAnfHHXYj56VL7XrT2dmuE4W9ow0vLny8rOc9I2IXSGvWDF5+2Y5MWb0aevVytwJnkGghV97atcuu0XzjjWE7C67C+vSxy5Vu3myn9U+e7DpRWEtNTSU2NrbEY7GxsaSmppbreU+sX29Hn/TuDQ0a2HeJL71kN+OORqW1twT7pm3kUSwtzbZBLlrkOon31qwRSUy0n19Sksj27a4Tha3x48eLz+cTY4z4fL4jOjLLev645efbDszatW2H5ogREdWZWRaO0kZuxMHaEomJiZKRkRHy86oQaNcOfv8dli+PurevwP/Wo37mGfjDH+CNN+BPf3KdSgFkZNgt15YutUNHR48O/xnFFWSMWSIiiYc/rk0ryjsrVtidwwcNis4iDlC9OjzxhH2rXreu3dj5rrt0GVyXdu+2C561awc//2z7ND74IOqK+LFoIVfeGTPGTv7p3991kuC74AJYsgTuvx9eecXuELNggetUlc/779vOzJdesh3Ta9bYST7ReiFxFFrIlTfy8uCtt+yIgJNPdp0mNGrVgn//Gz7/3Da5dOwIKSmwf7/rZNFvwwbbmXn99fbn7auvbFPKiSe6TuaEFnLljalT7U4qkT52/Hh07myblQYMgH/8w77F/+Yb16ki1jGn7x88CC++CM2b25FEzz5r3xl16OAucDgorQc02DcdtRKFLr9c5MwzRQ4edJ3ErWnTRE49VaRGDZHnnhM5cMB1oohyzOn7S5aIXHCBCMgXMTHSpGA2aEim9ocJdGanCpoff7TNC4MG2eVgK7OePe3V+DXXwMMP26v1n35ynSpipKSkkJubW+Ixk5vL70OHQtu2/P7DD/SvUYPL8vJYB2RnZ5OcnBz6RbfCTCX/rVOeeP11W8AHDHCdJDw0aGCngY8da5tcEhLsbFcHQ30jzeHT9LtjtxobvHs3JCfTrm5dxh/WB5Gbm0tKSgoQBqsqulLaZXqwb9q0EkXy80UaNrSLEKkjZWfbZicQueYakc2bXScKa4ULav0RZLL90ycrQK4/7TQRETHGlLrgljHG7aqKIYI2raigmDnTTl0fNMh1kvDUuDHMnm076D77DFq0gEmTXKcKW/945hnur16d1cDVwDCgY61a9B45Ejj2Oi2lNcsUv1qPaqVV92Df9Io8ivTsaa/Io2gadNCsXi3Stq29Or/pJp3if7hly4q+Pl/GxMhZpXRmHuuq+1hX69ECXY9ceW7TJpGqVUWGDXOdJHLk54sMHy5SrZrI6aeLfPKJ60Tu7d4t8sAD9mfp1FNF3n7brml/FEdbp8XZOuchpIVcee8f/7A/QmvXuk4SeTIyRJo1s1+/O+8U2bPHdSI3ZswQadzYfh2GDAnoXYq2kStVUYcO2Sn5l10GZ5/tOk3kqexT/H/+2U6l79HDLnc8bx6kpUG9esd9yKSkJNLS0vD5fBhj8Pl8pKWlkZSU5GHw8KSFXB2fL7+048cr40xOrxSf4r9/f+WY4n/woN11qVkzmDED/v53WLbMfu4eSEpKIisri0OHDpGVlVUpijhoIVfHy++3V0+9erlOEvk6d4aVK+HWW/83xX/lStepvLdiBVx8sV0t8sIL7eeYkgI1arhOFvG0kKuK277dTni5+WaIiXGdJjrUrWsnVk2bZpsdEhNhxAh7BRvp9u61s1zbtLGzXMePt+uknHOO62RRQwu5qrj0dNi3T8eOB0PhFP/u3aNjiv/HH9ux8yNGwG232WVmk5Iq3TKzwaaFXFWMCLz2mr1iTEhwnSY6NWhg9wUdN842P8TH26+5RNAU/82boV8/uOoq+67tyy/t51BZljgOMU8KuTHmdWPMFmOMrt0Z7RYvtsVFOzmDyxi7QcfKldC+PSQn26v0zZtdJzu2Q4fg1VdtZ+b778PTT0NmJnTq5DpZVPPqivxN4EqPjqXCmd8PsbFw442uk1QOZ5xh25NffNGObgnnKf4rV9rRJ0OH2uGVK1bA44/bXaNUUHlSyEVkLrDdi2OpMLZnD0yYAH/+s+2cU6FRpQr89a92mN5ZZ9mvf1KS3cgjHOTmwrBhtjPz++/tqo+ffgrnnus6WaURsjZyY0yyMSbDGJORk5MTqtMqL737ri3m2qziRtOmMH8+DB9uvxctW9qrdZc++cS+S/jXv+wopjVr4JZbtDMzxEJWyEUkTUQSRSSxQYMGoTqt8pLfb4vJRRe5TlJ5VasGTzwBCxfad0Xdutlx2Xv3hjbHr7/CTTfBlVdC9erwxRfwxhtwyimhzaEAHbWiyuvbb+0U8sGD9WorHBRO8b/vvtBO8T90yI4+adrUziV46inbFt65c/DPrY5KC7kqnzFj7JVX//6uk6hCtWrBf/5jO0Hz821H46OPBm+K/6pVcOmldgRNQgIsXw5PPqmdmWHAq+GHE4AFwHnGmI3GGJ0pEk327bNjmq+9Fk491XUadbjOne1V8YAB8M9//m/6u1d+/x0ee8xe9a9aZZtQvvjCXpWrsODVqJUbRaShiFQXkUYiMsaL46owMW0abNumnZzhrG5d+65p2jQ71jwxEZ57LvAp/rNn207V1FQ75HTNGvsHQ5vXwoo2raiy+f12y7IuXVwnUWUpnOJ/zTXwyCO2KeTHHyt+nC1b7CiUrl1t0f7sMzusUAcqhCUt5OrY1q2zV2UDB0LVqq7TqPJo0MB2RBZO8U9IsGt9l2eKf+E6802b2iGOjz9uj3H55cHPrY6bFnJ1bG+8Ya/IbrvNdRJVEYdP8b/9dnuVfqwp/qtW2fb2wYPt2PDly+0Ue13hMuxpIVdHd/CgXVq1WzfbtKIiT+PGdtLQqFG2g7JFC7t6ZfG287w8Oza9VSvbLOP3w5w5dr0UFRG0kKuj++QT2LRJOzkjXZUqcPfd/5vif/PNcN558PLLMHOmXV3xmWfs1P81a+zyxFW0NEQS/W6po/P7bXtrjx6ukygvFE7xnzTJfl/vvhuuvtq2i8+aZTd80OGlEama6wAqTP3yi91T8d57dSuuaFKtGvTpY28LFsDq1XZYYa1arpOpAGghV6UbNw4OHNBdgKJZhw72piKeNq2oI4nYZpWOHXX2nlIRQAu5OtK8ebB2rXZyKhUhtJCrI/n9dsp3nz6ukyilykELuSppxw47qiEpCWrXdp1GKVUOWshVSW+/bSeIaCenUhFDC7kqye+H1q3txgVKqYighVz9z9KldvafdnIqFVG0kKv/8fvtAkk33eQ6iVKqArSQKys31y6mdMMNcNJJrtMopSpAC7myJk+GXbu0WUWpCKSFXFl+P5xzDlxyieskSqkK0kKu4Lvv7GzOQYN0L0alIpAWcmW39qpaFW691XUSpdRx0EJe2e3fbzfV7dEDTjvNdRql1HHwpJAbY640xnxnjPnBGPM3L46pQuSDD+yO6drJqSogPT2duLg4qlSpQlxcHOnp6a4jVWoBr0dujKkKjAb+BGwEFhtjpovIqkCPrULA74fTT7f7cipVDunp6SQnJ5ObmwtAdnY2ycnJACQlJbmMVml5cUV+IfCDiPwkIvuBd4BrPTiuCrYNG+Djj+G22+zOMUqVQ0pKSlERL5Sbm0tKSoqjRMqLQn46sKHY/Y0Fj5VgjEk2xmQYYzJycnI8OK0K2Btv2E0kBg50nURFkPXr11focRV8XhTy0saryREPiKSJSKKIJDZo0MCD06qAHDxoR6t06QJNmrhOoyJI48aNK/S4Cj4vCvlG4Ixi9xsBP3twXBVMn30G69drJ6eqsNTUVGJjY0s8FhsbS2pqqqNEyotCvhg4xxjTxBhTA+gHTPfguCqY/H44+WS47jrXSVSESUpKIi0tDZ/PhzEGn89HWlqadnQ6FHAPl4gcMMb8BfgEqAq8LiLfBpxMBU9ODrz/Ptx1F9Ss6TqNikBJSUlauMOIJ0MVROQj4CMvjqVC4K23ID9fdwFSKkrozM7KRsQ2q7RvDy1auE6jlPKAFvLKZsECWL1aOzmViiJayCsbvx/q1IG+fV0nUUp5RAt5ZbJrF0ycCP362WKulIoKWsgrk3fesVu6abOKUlFFC3mUKnV1Or/fdnBeeKHreEopD+lKSVGotNXpXho8mKS8PHjhBd0FSKkoo1fkUai01eluystjH8DNNzvJpJQKHi3kUejwVehigP7AewD16ztIpJQKJi3kUejwVeiuB+oBM0491UkepVRwaSGPQoevTjcYWGcM14wc6S6UUipotJBHoeKr050NXA7s7NOHpP79XUdTSgWBFvIolZSURFZWFmuHDYMqVWj1/POuIymlgkQLeTQ7cMBu53bVVXaDZaVUVNJCHs0++gh++QWGDHGdRCkVRFrIo5nfD6edBldf7TqJUiqItJBHq02b4MMPYcAAqF7ddRqlVBBpIY9WY8fCoUMwcKDrJEqpINNCHo0OHYIxY6BzZzjnHNdplFJBpoU8Gs2ZAz/9pMvVKlVJaCGPRn4/nHQS9OrlOolSKgS0kEebbdtgyhRISoJatVynUUqFgBbyaJOeDvv369hxpSqRgAq5MeYGY8y3xphDxphEr0Kp4yQCr70GiYmQkOA6jVIqRAK9Iv8G6AXM9SCLCtTixfDNN9rJqVQlE9BWbyKyGsDo1mHhwe+H2Fi48UbXSZRSIRSyNnJjTLIxJsMYk5GTkxOq01Yee/bAhAnw5z9D3bqu0yilQqjMK3JjzKfAaaU8lSIi08p7IhFJA9IAEhMTpdwJVfm8+64t5tqsolSlU2YhF5EuoQiiAuT3Q9OmcNFFrpMopUJMhx9Gg2+/hQULYNAg0P4KpSqdQIcfXm+M2Qh0AD40xnziTSxVIWPG2BUOb7nFdRKllAOBjlqZCkz1KIs6Hvv2wbhxcO21cOqprtMopRzQppVIN22anZavnZxKVVpayCOd3w+NG0MX7ZNWqrLSQh7J1q2D2bPt5hFVq7pOo5RyRAt5JHvjDTtK5bbbXCdRSjmkhTxSHTzI3tGjmVOzJlXi4oiLiyM9Pd11KqWUA1rII9QXf/sbtbdv5+W8PESE7OxskpOTtZgrVQlpIY9Q+0aPZgswvdhjubm5pKSkuIqklHJEC3kk+uUXrvj9d8YC+Yc9tX79eheJlFIOaSGPROPGUR0YU8pTjRs3DnUapZRjWsgjjQj4/Ww591w2xMaWeCo2NpbU1FRHwZRSrmghjzTz5sHatZz66KOkpaXh8/kwxuDz+UhLSyMpKcl1QqVUiAW01opywO+3G0f06UNS7dpauJVSekUeUXbsgEmT7FZutWu7TqOUChNayCPJ229DXh4MGeI6iVIqjGghjyR+P7RqBW3auE6ilAojWsgjxdKlsGyZXa5WdwFSShWjhTxS+P0QEwM33eQ6iVIqzGghjwS5uZCeDn36QL16rtMopcKMFvJIMHky7NqluwAppUqlhTwS+P1w9tnQqZPrJEqpMKSFPNx9952dzTlokHZyKqVKpYU83I0ZY7dxGzDAdRKlVJgKqJAbY0YYY9YYY1YYY6YaY07yKpgC9u+HsWOhRw847TTXaZRSYSrQK/LZQAsRiQe+B4YFHkkV+eAD2LJFOzmVUscUUCEXkVkicqDg7kKgUeCRVKFNTz/N5qpVqd69u+7JqZQ6Ki/byAcCM4/2pDEm2RiTYYzJyMnJ8fC00WnqqFE0XL4c/8GDHADdk1MpdVRGRI79AmM+BUproE0RkWkFr0kBEoFeUtYBgcTERMnIyDiOuJXHf046ift37qQJkFXscZ/PR1ZWVukfpJSKasaYJSKSePjjZa5HLiJdyjjwrUB34IryFHFVDgcP0nvnTmZTsoiD7smplDpSoKNWrgQeAXqKSK43kRSffYYP3ZNTKVU+gbaRvwycAMw2xmQaY171IJPy+9lXpw6f1KpV4mHdk1MpVZpAR62cLSJniEirgtsdXgWrtHJy4P33qTl4MC+/9pruyamUKpPu2Rlu3noL8vNh0CCSWrTQwq2UKpNO0Q8nInaBrPbtoUUL12mUUhFCC3k4WbgQVq/WmZxKqQrRQh5O/H6oUwf69nWdRCkVQbSQh4tdu+Cdd6BfP1vMlVKqnLSQh4uJE+2WboMGuU6ilIowWsjDhd9vOzjbtXOdRCkVYbSQh4MVK+Drr20np+4CpJSqIC3k4WDMGKhRA26+2XUSpVQE0kLuWl6enQTUqxfUr+86jVIqAunMTtemToXffqv0Y8fz8/PZuHEjeXl5rqMo5VxMTAyNGjWievXq5Xq9FnLX/H5o0gQuu8x1Eqc2btzICSecQFxcHEb7CVQlJiJs27aNjRs30qRJk3J9jDatuPTjj/D553bIYZXK/a3Iy8ujfv36WsRVpWeMoX79+hV6d1q5q4dj39x/PweBRo89pntyghZxpQpU9HdBC7kjb48bR/0ZM/gI2ITuyamUOn5ayB357MEHaSiCv9hjubm5pKSkOMtU2dUJo6UR4uLi2Lp1q+sYAHTu3Bkv9tjNzMzko48+8iBRSXPmzKF79+5lvu7GG28kPj6e559/3tNzz58/v+j+q6++yrhx4zw7fnlpZ6cjPXNy2Awc/mOte3JGvoMHD1K1alXXMcpNRBARqgS5nyYzM5OMjAyuvvpqz4554MCBcr3ul19+Yf78+WRnZ3t2brCFvE6dOlx00UUA3HGHm7119IrchZ9/5hrgTeDwH0PdkxO4917o3Nnb2733VijCiBEjaNu2LfHx8Tz55JNFj1933XVccMEFnH/++aSlpRU9XqdOHZ544gnatWvHggULqFOnDikpKSQkJNC+fXt+/fVXAHJycujduzdt27albdu2fPXVVwBs27aNrl270rp1a26//XaOto95RY/71FNPMXLkyKKPb9GiBVlZWWRlZdGsWTPuvPNO2rRpw4YNGxg6dCiJiYmcf/75JT7no4mLi+PJJ5+kTZs2tGzZkjVr1gCwd+9eBg4cSNu2bWndujXTpk1j//79PPHEE0ycOJFWrVoxceJEWrZsyY4dOxAR6tevX3Ql279/fz799FPy8vK47bbbaNmyJa1bt+aLL74A4M033+SGG26gR48edO3atUSmxYsX07p1a3766acSj3ft2pUtW7bQqlUr5s2bV+JdxtatW4mLiys6dq9evbjyyis555xzePjhh4uO8fHHH9OmTRsSEhK44ooryMrK4tVXX+X5558vOm7xr3dmZibt27cnPj6e66+/nt9++w2w73AeeeQRLrzwQs4991zmzZtX5te6LFrIXRg7lmrAhJiYEg/rnpzhYdasWaxdu5avv/6azMxMlixZwty5cwF4/fXXWbJkCRkZGYwaNYpt27YBtni1aNGCRYsW0bFjR/bu3Uv79u1Zvnw5nTp14rXXXgPgnnvu4b777mPx4sVMmTKFwQXzB4YPH07Hjh1ZtmwZPXv2POo7s4oe91i+++47brnlFpYtW4bP5yM1NZWMjAxWrFjBl19+yYoVK8o8ximnnMLSpUsZOnRoUQFLTU3l8ssvZ/HixXzxxRc89NBD5Ofn8/TTT9O3b18yMzPp27cvF198MV999RXffvstZ555ZlFBW7hwIe3bt2f06NEArFy5kgkTJnDrrbcWjeRYsGABY8eO5fPPPy/KMn/+fO644w6mTZvGmWeeWSLn9OnTOeuss8jMzOSSSy455ueUmZnJxIkTWblyJRMnTmTDhg3k5OQwZMgQpkyZwvLly5k0aRJxcXHccccd3HfffaUe95ZbbuHZZ59lxYoVtGzZkuHDhxc9d+DAAb7++mteeOGFEo8fL21aCbVDh+yU/M6deWTwYFJSUli/fj2NGzcmNTVVt3YDeOEFp6efNWsWs2bNonXr1gDs2bOHtWvX0qlTJ0aNGsXUqVMB2LBhA2vXrqV+/fpUrVqV3r17Fx2jRo0aRe22F1xwAbNnzwbg008/ZdWqVUWv27VrF7t372bu3Lm89957AFxzzTXUq1ev1GwVPe6x+Hw+2rdvX3T/3XffJS0tjQMHDrB582ZWrVpFfHz8MY/Rq1evoiyF+WfNmsX06dOLCnteXl6pf5guueQS5s6di20t4VsAAAxrSURBVM/nY+jQoaSlpbFp0yZOPvlk6tSpw3//+1/uvvtuAJo2bYrP5+P7778H4E9/+hMnn3xy0bFWr15NcnIys2bN4o9//OMxM5fliiuu4MQTTwSgefPmZGdn89tvv9GpU6eicd3Fz12anTt3smPHDi699FIAbr31Vm644Yai54t/3bKysgLKC1rIgy49Pb1EsfYnJdHlxx9h+HCSkpK0cIchEWHYsGHcfvvtJR6fM2cOn376KQsWLCA2NpbOnTsXXSHGxMSUaBevXr160RCyqlWrFrXlHjp0iAULFlCrVq0jzlueIWcVPW61atU4dOhQ0f3iY5Nr165d9P9169YxcuRIFi9eTL169RgwYEC5xjHXrFnziCwiwpQpUzjvvPNKvHbRokUl7nfq1InRo0ezfv16UlNTmTp1KpMnTy66sj1a89Lh2QEaNmxIXl4ey5YtK1chL/51OfzzLPycin9eIuLp8NjSvm6B0KaVIEpPTyc5OZns7GxEhOzsbLY9+yz7YmPt2ioqLHXr1o3XX3+dPXv2ALBp0ya2bNnCzp07qVevHrGxsaxZs4aFCxdW+Nhdu3bl5ZdfLrqfmZkJ2KJWOPR05syZRe2pgR43Li6OpUuXArB06VLWrVtX6sfv2rWL2rVrc+KJJ/Lrr78yc+bMCp2/uG7duvHSSy8VFeJly5YBcMIJJ5R4l3DGGWewdetW1q5dy5lnnknHjh0ZOXJkUSEv/jX5/vvvWb9+/RF/HAqddNJJfPjhhzz66KPMmTOnzIxxcXEsWbIEgMmTJ5f5+g4dOvDll18Wff22b99e6udU6MQTT6RevXpFzUVvvfVW0dV5MGghD6KUlBRyc3OL7tcDrj14kAlVq0IpV2QqPHTt2pWbbrqJDh060LJlS/r06cPu3bu58sorOXDgAPHx8Tz++OMlmiXKa9SoUWRkZBAfH0/z5s159dVXAXjyySeZO3cubdq0YdasWRXu9D7acXv37s327dtp1aoVr7zyCueee26pH5+QkEDr1q05//zzGThwIBdffHGFP7dCjz/+OPn5+cTHx9OiRQsef/xxAC677DJWrVpV1NkJ0K5du6JMl1xyCZs2baJjx44A3HnnnRw8eJCWLVvSt29f3nzzzRJXy4f7wx/+wIwZM7jrrruOuPo/3IMPPsgrr7zCRRddVK5hng0aNCAtLY1evXqRkJBA34LtGHv06MHUqVOLOjuLGzt2LA899BDx8fFkZmbyxBNPlHme42WO9falzA825hngWuAQsAUYICI/l/VxiYmJ4sW41HBXpUqVEm8P/wK8BLQCMgP4ukej1atX06xZM9cxlAobpf1OGGOWiEji4a8N9Ip8hIjEi0gr4AMgeH9yItDhV1VDgMXADp/PSR6lVHQKqJCLyK5id2sDeplZTGpqKrGxsQAkAvHAuOrVdYihUspTAY9aMcakArcAO4GjrsVqjEkGkqHyTHopHJGSkpLCkOxsco2h4+jR9NWRKkopD5XZRm6M+RQ4rZSnUkRkWrHXDQNiRKTMKWGVpY28yJ490LAh9OkDb7zhOk1Y0jZypUqqSBt5mVfkItKlnOd9G/gQKHtub2Xz7ru2mFfyXYCUUsERUBu5MeacYnd7AmsCixOl/H5o2hQKFtZRSikvBTpq5V/GmG+MMSuArsA9HmSKLqtWwYIF9mpcN06IGIcvNnW4999/v8SUeKVcCnTUSm8RaVEwBLGHiGzyKljUGDMGqleH/v1dJ4kq6enpxMXFUaVKFSe7K2khV+FEZ3YG0759MG4cXHstnHqq6zRRo7SlD7zYXSk1NZXzzjuPLl268N133wHw2muv0bZtWxISEujduze5ubnMnz+f6dOn89BDD9GqVSt+/PHHUl+nVKhoIQ+m6dNh61bt5PTY4UsfQOC7Ky1ZsoR33nmHZcuW8d5777F48WLArlK3ePFili9fTrNmzRgzZgwXXXQRPXv2ZMSIEWRmZnLWWWeV+jqlQkVXPwwmvx8aN4Yu5R34o8rjaGt1B7K70rx587j++uuLJnD17NkTgG+++YbHHnuMHTt2sGfPHrp161bqx5f3dUoFg16RB+iobbVZWTB7NgwcCBG07VckONqEskAnmpW2TOmAAQN4+eWXWblyJU8++eRRl3Yt7+uUCgYt5AEora324SFD+HrAAOjRw75o4ECnGaNR8aUPCgW6u1KnTp2YOnUqv//+O7t372bGjBkA7N69m4YNG5Kfn1+iDf7w5UuP9jqlQiGymlb+/neYMMF1iiJt167l6/z8ovsGaPL779QaOxYSEiA9Hc44w13AKFV86QOvdldq06YNffv2pVWrVvh8vqI1sZ955hnatWuHz+ejZcuWRcW7X79+DBkyhFGjRjF58uSjvk6pUAhoGdvjddxT9P1++OQT7wMdp0mlLEi/CXgLWKLL1FaITtFXqiRPp+iHlcGDw2oEyENxcWRnZx/xuE+XqVVKhZC2kQcgGG21SilVUVrIA5CUlERaWho+nw9jDD6fj7S0NN1Q+Ti5aOZTKhxV9HchsppWwlBSUpIWbg/ExMSwbds26tev7+lu5UpFGhFh27ZtxMTElPtjtJCrsNCoUSM2btxITk6O6yhKORcTE0OjRo3K/Xot5CosVK9enSZNmriOoVRE0jZypZSKcFrIlVIqwmkhV0qpCOdkZqcxJgc4ciZN+ZwCbPUwjlc0V8VororRXBUTrrkgsGw+EWlw+INOCnkgjDEZpU1RdU1zVYzmqhjNVTHhmguCk02bVpRSKsJpIVdKqQgXiYU8zXWAo9BcFaO5KkZzVUy45oIgZIu4NnKllFIlReIVuVJKqWK0kCulVISLyEJujGlljFlojMk0xmQYYy50namQMeZuY8x3xphvjTHPuc5TnDHmQWOMGGNOcZ0FwBgzwhizxhizwhgz1RhzkuM8VxZ8734wxvzNZZZCxpgzjDFfGGNWF/xM3eM6U3HGmKrGmGXGmA9cZylkjDnJGDO54GdrtTGmg+tMAMaY+wq+h98YYyYYY8q/vGEZIrKQA88Bw0WkFfBEwX3njDGXAdcC8SJyPjDScaQixpgzgD8B611nKWY20EJE4oHvgWGughhjqgKjgauA5sCNxpjmrvIUcwB4QESaAe2Bu8IkV6F7gNWuQxzmReBjEWkKJBAG+YwxpwN/BRJFpAVQFejn1fEjtZALULfg/ycCPzvMUtxQ4F8isg9ARLY4zlPc88DD2K9dWBCRWSJyoODuQqD863Z670LgBxH5SUT2A+9g/yg7JSKbRWRpwf93Y4vS6W5TWcaYRsA1gN91lkLGmLpAJ2AMgIjsF5EdblMVqQbUMsZUA2LxsG5FaiG/FxhhjNmAvep1diV3mHOBS4wxi4wxXxpj2roOBGCM6QlsEpHlrrMcw0BgpsPznw5sKHZ/I2FSMAsZY+KA1sAit0mKvIC9ODjkOkgxZwI5wBsFTT5+Y0xt16FEZBO2Vq0HNgM7RWSWV8cP2/XIjTGfAqeV8lQKcAVwn4hMMcb8GfvXt0sY5KoG1MO+BW4LvGuMOVNCMMazjFyPAl2DnaE0x8olItMKXpOCbUJID2W2w5S2LVHYvHsxxtQBpgD3isiuMMjTHdgiIkuMMZ1d5ymmGtAGuFtEFhljXgT+BjzuMpQxph72HV4TYAcwyRhzs4iM9+L4YVvIReSohdkYMw7bNgcwiRC+tSsj11DgvYLC/bUx5hB2gZygb3tztFzGmJbYH57lBVuoNQKWGmMuFJFfXOUqlu9WoDtwRSj+4B3DRuCMYvcbESZNdsaY6tgini4i77nOU+BioKcx5mogBqhrjBkvIjc7zrUR2Cgihe9aJmMLuWtdgHUikgNgjHkPuAjwpJBHatPKz8ClBf+/HFjrMEtx72PzYIw5F6iB4xXYRGSliJwqInEiEof9QW8TiiJeFmPMlcAjQE8RyXUcZzFwjjGmiTGmBrYjarrjTBj713cMsFpE/uM6TyERGSYijQp+pvoBn4dBEafg53qDMea8goeuAFY5jFRoPdDeGBNb8D29Ag87YcP2irwMQ4AXCzoN8oBkx3kKvQ68boz5BtgP3Or4KjPcvQzUBGYXvFtYKCJ3uAgiIgeMMX8BPsGOKHhdRL51keUwFwP9gZXGmMyCxx4VkY8cZgp3dwPpBX+QfwJuc5yHgmaeycBSbDPiMjycqq9T9JVSKsJFatOKUkqpAlrIlVIqwmkhV0qpCKeFXCmlIpwWcqWUinBayJVSKsJpIVdKqQj3/wHk/Mzug+ICRgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#test x-values\n",
    "x_test = np.linspace(-8, 8, 100).reshape((1, -1))\n",
    "#predict on the test x-values\n",
    "y_test_pred = nn.forward(nn.weights, x_test)\n",
    "#visualize the function learned by the neural network\n",
    "plt.scatter(x_train.flatten(), y_train.flatten(), color='black', label='data')\n",
    "plt.plot(x_test.flatten(), y_test_pred.flatten(), color='red', label='learned neural network function')\n",
    "plt.legend(loc='best')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Part I: Inference for Bayesian Neural Networks Using HMC\n",
    "\n",
    "1. (**The Bayesian Neural Network Model**) We will implement the following Bayesian model for the data:<br>\n",
    "\\begin{align}\n",
    "\\mathbf{W} &\\sim \\mathcal{N}(0, 5^2 \\mathbf{I}_{D\\times D})\\\\\n",
    "\\mu^{(n)} &= g_{\\mathbf{W}}(\\mathbf{X}^{(n)})\\\\\n",
    "Y^{(n)} &\\sim \\mathcal{N}(\\mu^{(n)}, 0.5^2)\\\\\n",
    "\\end{align}\n",
    "where $g_{\\mathbf{W}}$ is a neural network with parameters $\\mathbf{W}$ represented as a vector in $\\mathbb{R}^{D}$ with $D$ being the total number of parameters (including biases).\n",
    "<br><br>\n",
    "Implement the log of the joint distribution in `autograd`'s version of `numpy`, i.e. implement $\\log \\left[p(\\mathbf{W})\\prod_{n=1}^N p(Y^{(n)} |\\mathbf{X}^{(n)} , \\mathbf{W}) \\right]$.\n",
    "***Hint:*** you'll need to write out the log of the various Gaussian pdf's and implement their formulae using `autograd`'s numpy functions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. (**Sample from the Posterior**) Use HMC to sample from the posterior of the Bayesian neural network in Exercise 1 with 1 hidden layer and 5 hidden nodes. Start with the following settings for your sampler:\n",
    "``` python\n",
    "params = {'step_size':1e-3, \n",
    "          'leapfrog_steps':50, \n",
    "          'total_samples':10000, \n",
    "          'burn_in':.1, \n",
    "          'thinning_factor':2,\n",
    "          'position_init': nn.weights}\n",
    "```\n",
    "  Note that you should initialize with the MLE model, otherwise convergence may be slow. Please feel free to tweak these design choices as you see fit!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. (**Visualize the Posterior Predictive**) Visualize 100 samples, randomly selected, from your posterior samples of $\\mathbf{W}$, by ploting their predicted values plus a random noise $\\epsilon \\sim \\mathcal{N}(0, 0.5^2)$ at 100 equally spaced x-values between -8 and 8:\n",
    "``` python \n",
    "x_test = np.linspace(-8, 8, 100)\n",
    "y_test = nn.forward(sample, x_test.reshape((1, -1))) \n",
    "y_test += np.random.normal(0, 0.5, size=y_test.shape)\n",
    "```\n",
    "  where `sample` is a posterior sample of $\\mathbf{W}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "4. (**Model Evaluation**) Discuss the fit of your model to the data. Discuss also what the posterior predictive tell you about the aleatoric and epistemic uncertainty of the model. Are these uncertainties what you'd want (think about want kind of uncertainties you want and where)?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Part II: Inference for BNNs Using Black-Box Variational Inference with the Reparametrization Trick\n",
    "\n",
    "1. (**BBVI with the Reparametrization Trick**) Implement BBVI with the reparametrization trick for approximating an arbitrary posterior $p(w| \\text{Data})$ by an isotropic Gaussian $\\mathcal{N}(\\mu, \\Sigma)$, where $\\Sigma$ is a diagonal matrix. See Lecture #17 or the example code from [autograd's github repo](https://github.com/HIPS/autograd/blob/master/examples/black_box_svi.py). \n",
    "\n",
    "    **Unit Test:** You shoud check that your implementation is correct by approximating the posterior of the following Bayesian logistic regression model:\n",
    "\\begin{align}\n",
    "w &\\sim \\mathcal{N}(0, 1)\\\\\n",
    "Y^{(n)} &\\sim Ber(\\text{sigm}(wX^{(n)} + 10))\n",
    "\\end{align}\n",
    "  where $w$, $Y^{(n)}$, $X^{(n)}$ are a real scalar valued random variables, and where the data consists of a single observation $(Y=1, X=-20)$.\n",
    "\n",
    "  The true posterior $p(w | Y=1, X=-20)$ should look like the following (i.e. the true posterior is left-skewed):\n",
    "<img src=\"./logistic_posterior.png\" style='height:200px;'>\n",
    "  Your mean-field variational approximation should be a Gaussian with mean around -0.321 and standard deviation around 0.876."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. (**Variational Inference for BNNs**) For the Bayesian model you implemented in Exercise 1 and 2, use BBVI with the reparametrization trick to approximate the posterior of the Bayesian neural network with a mean-field Gaussian variational family (i.e. an isotropic Gaussian). Please set learning rate and maximum iteration choices as you see fit!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. (**Visualize the Posterior Predictive**) Visualize 100 samples $\\mathbf{W}^s$ from your approximate posterior of $\\mathbf{W}$ by ploting the neural network outputs with weight $\\mathbf{W}^s$ plus a random noise $\\epsilon \\sim \\mathcal{N}(0, 0.5^2)$ at 100 equally spaced x-values between -8 and 8:\n",
    "``` python \n",
    "x_test = np.linspace(-8, 8, 100)\n",
    "y_test = nn.forward(sample, x_test.reshape((1, -1))) \n",
    "y_test += np.random.normal(0, 0.5, size=y_test.shape)\n",
    "```\n",
    "  where `sample` is a sample from the approximate posterior of $\\mathbf{W}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "4.  (**Model Evaluation**) Discuss the fit of your model to the data. Discuss also what the posterior predictive tell you about the aleatoric and epistemic uncertainty of the model. Are these uncertainties what you'd want (think about want kind of uncertainties you want and where)?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Part III: Comparison \n",
    "\n",
    "1. (**Model Evaluation**) Compare the posterior predictive visualization from BBVI with the reparametrization trick to the one you obtained using HMC. Can you say whether or not your posterior approximation is good? How does approximating the posterior effect our estimation of epistemic and aleatoric uncertainty?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. (**Quality of Posterior Vs Quality of Posterior Predictive**) Based on your convergence diagnostics of your HMC sampler, do you believe that you've obtained a set of representative samples from the BNN posterior? Why or why not? \n",
    "\n",
    "  Based on your understanding of BNN posteriors, do you believe that a mean-field Gaussian is a good approximation of the BNN posterior? Why or why not?\n",
    "\n",
    "  Which set of samples, HMC or VI, do you believe better captures the BNN posterior?\n",
    "\n",
    "  Based on the empirical results from above, do you believe there is a clear positive relationship between posterior approximation and quality of posterior predictive uncertainties - i.e. is it true that the more I capture of the BNN posterior the beter the quality of my posterior predictive uncertainties?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. (**Speed and Scalability**) Compare the two inference methods, HMC and BBVI, and discuss the scalability of each with respect to the number of model parameters and number of data points (how does increasing the number of model parameters or the number of data points affect each inference method?). "
   ]
  }
 ],
 "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}