NB.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Python: 3.8.5 (default, Sep  3 2020, 21:29:08) [MSC v.1916 64 bit (AMD64)]\n",
      "Sklearn: 0.23.2\n",
      "Matplotlib: 3.3.2\n",
      "NumPy: 1.19.2\n"
     ]
    }
   ],
   "source": [
    "import sys\n",
    "import sklearn\n",
    "import matplotlib\n",
    "import numpy as np\n",
    "from datetime import datetime\n",
    "\n",
    "from scipy import optimize\n",
    "from numpy import log, exp\n",
    "    \n",
    "from scipy.stats import multivariate_normal as mnorm\n",
    "from scipy.stats import norm\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# python magic function\n",
    "%matplotlib inline\n",
    "\n",
    "print('Python: {}'.format(sys.version))\n",
    "print('Sklearn: {}'.format(sklearn.__version__))\n",
    "print('Matplotlib: {}'.format(matplotlib.__version__))\n",
    "print('NumPy: {}'.format(np.__version__))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Training Data: (60000, 28, 28)\n",
      "Training Labels: (60000,)\n",
      "Testing Data: (10000, 28, 28)\n",
      "Testing Labels: (10000,)\n"
     ]
    }
   ],
   "source": [
    "from keras.datasets import mnist\n",
    "\n",
    "(X_train, y_train), (X_test, y_test) = mnist.load_data()\n",
    "\n",
    "print('Training Data: {}'.format(X_train.shape))\n",
    "print('Training Labels: {}'.format(y_train.shape))\n",
    "print('Testing Data: {}'.format(X_test.shape))\n",
    "print('Testing Labels: {}'.format(y_test.shape))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAqYAAAK7CAYAAAAgM03aAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAzt0lEQVR4nO3dfZznc70//udr7VqXu4hcnHKdiNhI4viik6skiRJZoqJvTqjf4eiUIxJCOrmILpwo6ahzxEr8cA6LXP04js5NWrnIWiw2rF1r7ca8f3/MuH3n+H5eHzOfmZ15zsz9frvtrbm9Hvt+f56z9Zp5zHtmXpWmaQIAAIbbuOEeAAAAIhRTAACSUEwBAEhBMQUAIAXFFACAFBRTAABSUEyTKqVcXEr55nDPAfSNPQsjh/2al2LaR6WUx0opz5RSlu+19rlSyvRhHKtjPZtycSnlpV5/lhruuWCwjMI9O7GU8uNSyrxSytOllP9nuGeCwTLa9uvrSimrlFLmlFJ+O9yzjBSKaf+Mj4ijh3uI/mpTOM9ommaFXn9eG9LBYMkbTXv2xIh4R0SsExEfiIi/L6XsPpRzwRI2mvbr606PiD8M1SyjgWLaP2dGxDGllJXeGJRS1i2lNKWU8b3WppdSPtfz9iGllNtKKf9USplbSnm0lLJdz/qsUsqzpZRPv+G2q5ZSbiilzC+l3FxKWafXvTfuyZ4vpTxYStmvV3ZxKeWCUso1pZQF0f1JDMai0bRnD46Ik5umeaFpmj9ExI8i4pAB/NtANqNpv0YpZduI2CwiLhrIP8pYo5j2zz0RMT0ijunw+m0i4r8j4i0R8fOIuCwito6IDSNiakScV0pZodffPzAiTo6IVSPivoi4NCKi51sdN/Tc460RcUBEnF9K2bTXtZ+KiFMiYsWIqH0L4YieTfefpZR9O3yfILNRsWdLKStHxFoR8btey7+LiN7Xw0g3KvZrzz2WiojvRcQXI8L/93s/KKb9d0JEHFlKWa2Da//UNM1FPd8y/0VEvD0ivtE0zaKmaa6PiMXRvYFe95umaW5pmmZRRHwtIrYtpbw9IvaMiMd67vVq0zT3RsTlEfHxXtdOa5rmtqZpupqmeaXFLOdE97cF3xoR/xgRF5dS/rqD9wmyGw179vVPpi/2Wnsxuj8pwmgyGvZrRMRREXFX0zT/2cH7MaaNf/O/Qm9N09xfSrk6Ir4S/f+5kWd6vb2w535vXOv91dysXq/7Uinl+eh+arJORGxTSpnb6++Oj4hLWl3bSs9Ge901pZRLI2KfiLjtzd8NGDlGyZ59qec/J0XEK73env/m7wKMHKNhv5ZS1oruYrpVP+cnFNNOfT0i7o2Is3qtLej5z+UiYl7P22sM8HXe/vobPd9+WCUinoruDXFz0zS7tLm2v986aCKi9HtCGBlG9J5tmuaFUsrsiNgiur/FGD1v/36A80JGI3q/RsT7ImLNiHiglBIRsWxELFtKeToi/sovGrfnW/kdaJrm4ej+NsFRvdbmRMSTETG1lLJUKeUzEbHBAF9qj1LK9qWUpaP752DuappmVkRcHREblVIOKqVM6PmzdSllk77euJTy8VLKCqWUcaWUXaP752+uGuC8kNJo2LMR8dOIOL6UsnIpZeOIOCwiLh7gvJDOKNiv10bEuhExpefPCRHxXxExRSl9c4pp574REcu/Ye2wiDg2Ip6L7l9KuH2Ar/Hz6P7K8fno/pbAgRERTdPMj4hdI2L/6P7q7unoPpJiYj/ufXR0b/K50f2bkIc1TTN9gPNCZiN9z349Ih6JiJkRcXNEnNk0zf87wHkhqxG7X3t+pvXp1/9E98+D/6Xnbd5EaRq/LAYAwPDzxBQAgBQUUwAAUlBMAQBIQTEFACAFxRQAgBQUUwAAUlBMAQBIQTEFACAFxRQAgBQUUwAAUlBMAQBIQTEFACAFxRQAgBQUUwAAUlBMAQBIQTEFACAFxRQAgBQUUwAAUlBMAQBIQTEFACAFxRQAgBTGtwtLKc1QDQLZNE1ThnuG/rJnGctG2p61XxnLavvVE1MAAFJQTAEASEExBQAgBcUUAIAUFFMAAFJQTAEASEExBQAgBcUUAIAUFFMAAFJQTAEASEExBQAgBcUUAIAUFFMAAFJQTAEASEExBQAgBcUUAIAUFFMAAFJQTAEASEExBQAgBcUUAIAUFFMAAFJQTAEASEExBQAgBcUUAIAUFFMAAFJQTAEASGH8cA8AMBZstdVW1eyLX/xiNTv44INbrv/0pz+tXnPuuedWs3vvvbeaAQw3T0wBAEhBMQUAIAXFFACAFBRTAABSUEwBAEhBMQUAIIXSNE09LKUe0mdLLbVUNZs8efKgv167o2eWW265luvvfOc7q9f87d/+bTX79re/Xc0OOOCAavbKK69Us29961vV7KSTTqpmg61pmjJkLzZI7NnhNWXKlGp24403VrNJkyYN6hwvvvhiNXvLW94yqK+VyUjbs/Yrb+aDH/xgNbv00kur2Y477ljNHnzwwQHNNFhq+9UTUwAAUlBMAQBIQTEFACAFxRQAgBQUUwAAUlBMAQBIYfxwDzBc1l577Wq29NJLV7Ptttuu5fr2229fvWallVaqZvvuu281G0pPPPFENTvnnHOq2cc+9rFqNn/+/Gr2u9/9rprdfPPN1QyG2/ve975qdvnll1ezdkfDtTu2r7aPFi9eXL2m3ZFQ73//+6vZvffeW83avR757bDDDtWs3f9errjiiiUxDn209dZbV7O77757CCcZOp6YAgCQgmIKAEAKiikAACkopgAApKCYAgCQgmIKAEAKo/q4qClTplSzG2+8sZq1O9ZlpOvq6mq5fvzxx1eveemll6rZpZdeWs1mz55dzV544YVq9uCDD1YzGCzLLbdcNdtyyy2r2c9+9rNqtuaaaw5oplYeeuihlutnnHFG9ZrLLrusmt12223VrN3HgdNOO62akd9OO+1Uzd7xjndUM8dFLXnjxtWfEa633nrVbJ111qlmpZQBzTScPDEFACAFxRQAgBQUUwAAUlBMAQBIQTEFACAFxRQAgBRG9XFRjz/+eDV77rnnqlmW46LuuuuuajZ37txq9oEPfKCaLV68uOX6JZdc0ue5YDT4wQ9+UM0OOOCAIZykvdrRVSussEL1mptvvrmatTs2aPPNN+/zXIwsBx98cDW74447hnAS3qjdMXOHHXZYNWt3dN2MGTMGNNNw8sQUAIAUFFMAAFJQTAEASEExBQAgBcUUAIAUFFMAAFIY1cdFPf/889Xs2GOPrWZ77rlnNfuv//qvluvnnHNO3wfr5b777qtmu+yySzVbsGBBNdt0002r2dFHH92nuWC02GqrrVquf/jDH65eU0rp6LXaHdP061//upp9+9vfrmZPPfVUy/Xax6KIiBdeeKGa/c3f/E016/T9Jr9x4zyHyurCCy/s6LqHHnpokCfJwf9SAQBIQTEFACAFxRQAgBQUUwAAUlBMAQBIQTEFACCFUX1cVDtXXnllNbvxxhur2fz581uub7HFFtVrPvvZz1azdsfEtDsSqp3f//731ezwww/v6J6Q2ZQpU6rZDTfc0HJ90qRJ1Wuapqlm1157bTU74IADqtmOO+5YzY4//vhqVjtKZs6cOdVrfve731Wzrq6uatbuCK0tt9yymt17773VjKG1+eabt1xfffXVh3gS+mry5MkdXVf72DbSeWIKAEAKiikAACkopgAApKCYAgCQgmIKAEAKY/a38tuZN29ev6958cUXO3qtww47rJr94he/qGbtfrMWRqONNtqomh177LHVrPYbr3/+85+r18yePbua/eQnP6lmL730UjX7zW9+01E2lJZddtlq9nd/93fV7MADD1wS49CBPfbYo+V6u/9uWfLanYqw3nrrdXTPJ598stNxUvPEFACAFBRTAABSUEwBAEhBMQUAIAXFFACAFBRTAABScFzUIDnxxBOr2VZbbVXNdtxxx2q28847V7Prr7++T3PBSDJx4sRq9u1vf7ua1Y7IiYiYP39+y/WDDz64es0999xTzcbqsTtrr732cI9AH7zzne/s9zW///3vl8Ak9Nbu41e7o6T++Mc/VrPax7aRzhNTAABSUEwBAEhBMQUAIAXFFACAFBRTAABSUEwBAEjBcVGDZMGCBdXssMMOq2b33ntvNfvRj35UzW666aZq1u6om+9973st15umqV4DQ+U973lPNWt3JFQ7H/3oR1uu33zzzR3dD0abu+++e7hHSGfSpEnVbPfdd69mU6dObbm+6667djTHySefXM3mzp3b0T2z88QUAIAUFFMAAFJQTAEASEExBQAgBcUUAIAUFFMAAFJwXNQQeOSRR6rZIYccUs0uuuiianbQQQd1lC2//PIt13/6059Wr5k9e3Y1g8H0ne98p5qVUqpZu6OfHAv1P40bV38e0dXVNYSTkMUqq6wypK+3xRZbVLPaPt95552r17ztbW+rZksvvXQ1O/DAA6tZu32ycOHCanbXXXe1XF+0aFH1mvHj61XsP//zP6vZaOWJKQAAKSimAACkoJgCAJCCYgoAQAqKKQAAKSimAACk4LioYXbFFVdUs4ceeqiatTtW54Mf/GA1O/XUU1uur7POOtVrTjnllGr25JNPVjNoZc8996xmU6ZMqWZN01Szq666aiAjjSntjoRq92983333LYFpGGy1o4za/Xf7/e9/v5p99atfHfBMb7T55ptXs9pxUa+++mr1mpdffrmaPfDAA9Xsxz/+cTW75557qlm7I+ieeeaZlutPPPFE9Zpll122ms2YMaOajVaemAIAkIJiCgBACoopAAApKKYAAKSgmAIAkIJiCgBACo6LSuz++++vZvvtt181+8hHPlLNLrroopbrn//856vXvOMd76hmu+yySzWDVtodjbL00ktXs2effbaa/eIXvxjQTCPRxIkTq9mJJ57Y0T1vvPHGavYP//APHd2ToXXEEUe0XJ85c2b1mu22225JjdPS448/Xs2uvPLKlut/+MMfqtfceeedAx1p0Bx++OEt11dbbbXqNY8++uiSGmdE8sQUAIAUFFMAAFJQTAEASEExBQAgBcUUAIAUFFMAAFJwXNQINXfu3Gp2ySWXVLMLL7yw5fr48fX/Keywww7VbKeddqpm06dPr2bQX4sWLapms2fPHsJJhlbtWKjjjz++es2xxx5bzZ544olqdtZZZ1Wzl156qZqR3+mnnz7cI4wJH/zgB/t9zeWXX74EJhm5PDEFACAFxRQAgBQUUwAAUlBMAQBIQTEFACAFv5Wf2Oabb17NPv7xj1ezrbfeupq1++37mgceeKCa3XLLLf2+H3TiqquuGu4RlpgpU6ZUs9pv2H/yk5+sXjNt2rRqtu+++/Z5LmDJu+KKK4Z7hFQ8MQUAIAXFFACAFBRTAABSUEwBAEhBMQUAIAXFFACAFBwXNQTe+c53VrMvfvGL1WyfffapZmusscaAZnqj1157rZrNnj27mnV1dQ3qHIx+pZSOsr333ruaHX300QMZaUh8+ctfrmb/+I//WM0mT57ccv3SSy+tXnPwwQf3fTCARDwxBQAgBcUUAIAUFFMAAFJQTAEASEExBQAgBcUUAIAUHBfVT7Vjmg444IDqNe2OhFp33XUHOlK/3HPPPS3XTznllOo1V1111ZIahzGoaZqOsnZHpJ1zzjnV7Mc//nHL9eeee656zfvf//5qdtBBB1WzLbbYopq97W1vq2aPP/54Nbvuuutarp9//vnVa4Bc2h2Ft9FGG1WzO++8c0mMk5onpgAApKCYAgCQgmIKAEAKiikAACkopgAApKCYAgCQwpg9Lmr11VevZu9617uq2XnnnddyfeONNx7wTP1x1113VbMzzzyzmk2bNq3leldX14BngiVpqaWWqmZHHHFENdt3331brs+bN696zTve8Y6+D9ZHt99+ezW76aabqtkJJ5ww6LMAQ6vdUXjjxnlG2Jt/DQAAUlBMAQBIQTEFACAFxRQAgBQUUwAAUlBMAQBIYcQfF7XKKqtUsx/84AfVbMqUKdVs/fXXH8hI/dLuCJmzzjqrml133XXVbOHChQOaCZakO+64o5rdfffd1Wzrrbfu6PXWWGONluvtjoxr57nnnqtml112WTU7+uijO3o9YHTbdtttq9nFF188dIMk4YkpAAApKKYAAKSgmAIAkIJiCgBACoopAAApKKYAAKSQ5riobbbZppode+yx1ex973tfNfurv/qrAc3UHy+//HI1O+ecc6rZqaeeWs0WLFgwoJkgoyeeeKKa7bPPPtXs85//fDU7/vjjBzTTG5199tnV7IILLqhmDz/88KDOAYwOpZThHmHE8MQUAIAUFFMAAFJQTAEASEExBQAgBcUUAIAUFFMAAFJIc1zUxz72sY6yTj3wwAPV7Oqrr65mr776asv1s846q3rN3Llz+zwXjGWzZ8+uZieeeGJHGcBQufbaa1uuf+ITnxjiSUYuT0wBAEhBMQUAIAXFFACAFBRTAABSUEwBAEihNE1TD0uphzDKNU1ThnuG/rJnGctG2p61XxnLavvVE1MAAFJQTAEASEExBQAgBcUUAIAUFFMAAFJQTAEASEExBQAgBcUUAIAUFFMAAFJQTAEASEExBQAgBcUUAIAUFFMAAFJQTAEASEExBQAgBcUUAIAUFFMAAFJQTAEASEExBQAgBcUUAIAUStM0wz0DAAB4YgoAQA6KKQAAKSimAACkoJgCAJCCYgoAQAqKaVKllItLKd8c7jmAvrFnYeSwX/NSTPuolPJYKeWZUsryvdY+V0qZPoxjdayUsl8p5fZSyssj9X2Adkbhnv2rUsq0UsrzpZQnSin/e7hngsEyCvfrt0spD5VS5pdSZpRSDh7umUYKxbR/xkfE0cM9RH+VUpZqsfx8RHw3Ir41tNPAkBpNe/ZnEfGniFg9Ij4cEaeWUj4wpIPBkjWa9uuCiPhIREyOiE9HxNmllO2GdLARSjHtnzMj4phSykpvDEop65ZSmlLK+F5r00spn+t5+5BSym2llH8qpcwtpTxaStmuZ31WKeXZUsqn33DbVUspN/R8xXVzKWWdXvfeuCd7vpTyYCllv17ZxaWUC0op15RSFkTE//XJq2maf2+a5pcR8dSA/1Ugr1GxZ0spK0TEThFxStM0f2ma5ncR8W8R8ZkB/wtBHqNiv0ZENE3z9aZpZjRN09U0zV0RcWtEbDvQf6CxQDHtn3siYnpEHNPh9dtExH9HxFsi4ucRcVlEbB0RG0bE1Ig4r+cT0OsOjIiTI2LViLgvIi6NiOj5VscNPfd4a0QcEBHnl1I27XXtpyLilIhYMSJ+2+G8MNKNlj1b3vCfr7+9WYfvF2Q0Wvbr/1BKWbZnjt93+H6NKYpp/50QEUeWUlbr4No/NU1zUdM0r0XELyLi7RHxjaZpFjVNc31ELI7uDfS63zRNc0vTNIsi4msRsW0p5e0RsWdEPNZzr1ebprk3Ii6PiI/3unZa0zS39Xy19koHs8JoMeL3bNM08yPitoj4x1LKMqWULSNi34hYroP3CTIb8fu1he9HxO8i4roO3qcxZ/yb/xV6a5rm/lLK1RHxlYj4Qz8vf6bX2wt77vfGtd5fzc3q9bovlVKej4i1ImKdiNimlDK3198dHxGXtLoWxrJRtGcPjIjv9fy9R6P76c67+vRewAgxivZrRESUUs6M7u9sfKBpmqYv14x1imlnvh4R90bEWb3WFvT853IRMa/n7TUG+Dpvf/2Nnm8/rBLdPxM6KyJubppmlzbX2gDwf4z4Pds0zczofpLz+v1/HhH/34CmhZxG/H7tuedJEfGhiNixaZp5b/b36eZb+R1omubh6P42wVG91uZExJMRMbWUslQp5TMRscEAX2qPUsr2pZSlo/vnYO5qmmZWRFwdERuVUg4qpUzo+bN1KWWTvt64Z8ZlovuLk3E93x6cMMB5IaVRsmc3KaWsWEpZupQyNSJ2jYjvDHBeSGeU7Nd/iO6fQ92laZrnBjjnmKKYdu4bEbH8G9YOi4hjI+K5iNg0Im4f4Gv8PLq/cnw+IraK7m/lvf7zZrtGxP7R/dXd0xFxekRM7Me9D4rub2tcEBH/q+ftHw1wXshspO/Z3aL7W/gvRMT/jojdez5Zw2g00vfrqRGxdkQ8VEp5qefPVwc475hQ/MgDAAAZeGIKAEAKiikAACkopgAApKCYAgCQgmIKAEAKiikAACkopgAApKCYAgCQgmIKAEAKiikAACkopgAApKCYAgCQgmIKAEAKiikAACkopgAApKCYAgCQgmIKAEAKiikAACkopgAApKCYAgCQgmIKAEAKiikAACmMbxeWUpqhGgSyaZqmDPcM/WXPMpaNtD1rvzKW1farJ6YAAKSgmAIAkIJiCgBACoopAAApKKYAAKSgmAIAkIJiCgBACoopAAApKKYAAKSgmAIAkIJiCgBACoopAAApKKYAAKSgmAIAkIJiCgBACoopAAApKKYAAKSgmAIAkIJiCgBACoopAAApKKYAAKSgmAIAkIJiCgBACoopAAApKKYAAKSgmAIAkIJiCgBACoopAAApKKYAAKSgmAIAkIJiCgBACoopAAApKKYAAKQwfrgHIIfjjz++mp100knVbNy4+tc2O+20UzW7+eab+zQXAAy1FVdcsZqtsMIK1ezDH/5wy/XVVlutes13vvOdarZo0aJqNlp5YgoAQAqKKQAAKSimAACkoJgCAJCCYgoAQAqKKQAAKTguaow55JBDWq4fd9xx1Wu6uro6eq2maTq6DgAGw7rrrlvN2n3e23bbbavZZpttNpCR/i9rrrlmNTvqqKMG9bVGAk9MAQBIQTEFACAFxRQAgBQUUwAAUlBMAQBIQTEFACAFx0WNMeuss07L9WWWWWaIJ4G8ttlmm2o2derUarbjjjtWs0033bSjWY455piW60899VT1mu23376a/exnP6tmd911V98HgyG08cYbV7MvfelL1ezAAw+sZssuu2w1K6VUs1mzZlWz+fPnt1zfZJNNqtfst99+1ez888+vZjNmzKhmI5knpgAApKCYAgCQgmIKAEAKiikAACkopgAApKCYAgCQguOiRqGdd965mh155JH9vl+7Iyn23HPPavbMM8/0+7VgqHzyk5+sZmeffXY1W3XVVatZuyNmpk+fXs1WW221anbmmWdWs07maPda+++/f79fC/pr8uTJ1ez0009vud5uv6644ooDnumNHnrooWq22267VbMJEya0XG/3ebTdx5R22WjliSkAACkopgAApKCYAgCQgmIKAEAKiikAACkopgAApOC4qBFq++23r2YXXXRRNWt3TEdNu+NqZs6c2e/7wWAbP77+oey9731vy/Uf/ehH1WuWW265anbLLbdUs5NPPrma/fa3v61mEydOrGa//OUvW67vuuuu1Wvaueeeezq6DgbLxz72sWr2uc99bsjmeOSRR6rZLrvsUs1mzZpVzTbccMMBzYQnpgAAJKGYAgCQgmIKAEAKiikAACkopgAApKCYAgCQguOiRqhPf/rT1Wyttdbq9/2mT59ezX7605/2+34wlKZOnVrNLrzwwn7f74Ybbqhmn/zkJ6vZvHnz+v1ab3bPTo6FeuKJJ6rZT37yk37fDwbTJz7xiUG932OPPVbN7r777mp23HHHVbN2R0K1s8kmm3R0Hf+HJ6YAAKSgmAIAkIJiCgBACoopAAApKKYAAKTgt/ITW3XVVavZZz7zmWrW1dVVzebOndty/Zvf/Gaf54LhcPLJJ1ezr371q9WsaZqW6+eff371muOPP76adfqb9+187WtfG9T7HXXUUdVszpw5g/pa0F+HHXZYNTv88MNbrl9//fXVax5++OFq9uyzz/Z9sEGw+uqrD+nrjUaemAIAkIJiCgBACoopAAApKKYAAKSgmAIAkIJiCgBACo6LGmbrrrtuNbv88ssH/fXOPffclus33XTToL8W9NcJJ5xQzdodCbV48eJqdt1117VcP+6446rXLFy4sJq1s8wyy1SzXXfdtZqtvfba1ayU0nK93RFv06ZNq2Yw3J566qlqduKJJw7dIEvAtttuO9wjjHiemAIAkIJiCgBACoopAAApKKYAAKSgmAIAkIJiCgBACo6LGma77757Ndt88807uud//Md/VLOzzz67o3vCYFlppZWq2RFHHFHNmqapZrUjoSIi9t57776M1WcbbrhhNbv00kur2VZbbdXR6/3bv/1by/Uzzjijo/vBWHLUUUdVs+WXX37QX+/d7353v6+5/fbbq9kdd9wxkHFGJE9MAQBIQTEFACAFxRQAgBQUUwAAUlBMAQBIQTEFACCF0u4IllJKPaTP2h1Xc/HFF1ezdkdZtDteYr/99qtmzzzzTDXjf2qapgz3DP01EvbsW9/61mr21FNPdXTP9ddfv5q98sorLdcPPfTQ6jV77bVXNdtss82q2QorrFDN2n2sbZfts88+Ldd//etfV68Zq0banh0J+3WoLbfcci3X3/Wud1Wv+frXv17N9thjj47mGDeu/tyuq6ur3/dr97Ftp512qmaPPPJIv19rpKjtV09MAQBIQTEFACAFxRQAgBQUUwAAUlBMAQBIQTEFACCF8cM9wGix7rrrVrPLL7980F/v0UcfrWaOhCKzxYsXV7M5c+ZUs9VWW62a/elPf6pm7Y5i6kS7Y1/mzZtXzdZcc81q9uc//7maORaKkWjChAnV7D3veU81q32+bLd/Fi5cWM3a7dc77rijmu2+++7VrHakVTvjx9frVu1IuIiIs88+u5q1+1g6knliCgBACoopAAApKKYAAKSgmAIAkIJiCgBACoopAAApOC5qkBx33HHVrKura9Bf71vf+tag3xOGwty5c6vZ3nvvXc2uvvrqarbKKqtUs0ceeaTl+rRp06rXXHzxxdXs+eefr2aXXXZZNWt33E276yCrpZdeupq1O27pV7/6Vb9f66STTqpmN954YzW77bbbqlm7jxvt7rnZZptVs5p2x92ddtpp1ezxxx+vZldeeWU1W7RoUZ/mysgTUwAAUlBMAQBIQTEFACAFxRQAgBQUUwAAUlBMAQBIwXFR/TRlypSW67vuuuugv1a742wefPDBQX89GG533XVXNWt33MpQ2mGHHarZjjvuWM3aHRv36KOPDmgmWFImTJhQzdod4XTsscd29HrXXntty/Vzzz23ek27I+jafdy45pprqtm73/3uarZ48eJqdsYZZ7Rcb3fE1Ec/+tFqdumll1azf//3f69mp59+ejV74YUXqlnNfffd1+9rOuWJKQAAKSimAACkoJgCAJCCYgoAQAqKKQAAKZSmaephKfVwjHr22Wdbrq+88sod3e/OO++sZh/60Ieq2UsvvdTR69F3TdOU4Z6hv+zZJW+33XarZu1+y7fdx9o111yzms2ZM6dvgzHi9myW/brUUktVs1NOOaWaHXPMMdVswYIF1ewrX/lKNbvssstarrf7TfL3vve91ey8887r6LqHH364mn3hC1+oZjfddFPL9UmTJlWv2W677arZgQceWM322muvarb88stXs3ZmzZrVcn299dbr6H7t1ParJ6YAAKSgmAIAkIJiCgBACoopAAApKKYAAKSgmAIAkILjovrptddea7ne1dXV0f0OPvjgavYv//IvHd2TwTHSjp6JsGeHW+3jQ4TjoobCSNuzWfZru+OPzj333Gr28ssvV7PDDz+8ml1//fXVbJtttmm5fuihh1avaXe04rLLLlvNvvGNb1Sziy66qJrVjlQaagcccEA1+9SnPtXRPb/85S+3XG93fFanHBcFAEBqiikAACkopgAApKCYAgCQgmIKAEAKiikAACk4LqqFdsdEHHLIIS3XOz0uav31169mM2fO7OieDI6RdvRMxNjds0Npt912q2bXXHNNNXNc1JI30vZslv06e/bsarbaaqtVs0WLFlWzGTNmVLPll1++mm244YbVrBMnnnhiNTvttNOqWbuj3xgcjosCACA1xRQAgBQUUwAAUlBMAQBIQTEFACAFxRQAgBTGD/cAw2XKlCnVbOedd65mtWOhFi9eXL3me9/7XjV75plnqhmQT7sj3mAkevrpp6tZu+OiJk6cWM222GKLjmapHbl2yy23VK+58sorq9ljjz1WzRwJlZMnpgAApKCYAgCQgmIKAEAKiikAACkopgAApKCYAgCQwpg9LmqllVaqZmussUa/7/fkk09Ws2OOOabf9wNyuvXWW6vZuHH1r/VrR83BcNthhx2q2d57713Nttxyy2r27LPPVrMf//jH1eyFF15oud7uSEZGF09MAQBIQTEFACAFxRQAgBQUUwAAUlBMAQBIQTEFACCFMXtcFEAn7r///mr20EMPVbP111+/mm2wwQbVbM6cOX0bDDo0f/78anbJJZd0lEGnPDEFACAFxRQAgBQUUwAAUlBMAQBIQTEFACAFxRQAgBTG7HFRM2bMqGa33357Ndt+++2XxDjAKHDqqadWswsvvLCanXLKKdXsyCOPbLn+wAMP9H0wgBHCE1MAAFJQTAEASEExBQAgBcUUAIAUFFMAAFIoTdPUw1LqIYxyTdOU4Z6hv+zZ4TVp0qRq9stf/rKa7bzzztXsV7/6Vcv1Qw89tHrNggULqtloNtL2rP3KWFbbr56YAgCQgmIKAEAKiikAACkopgAApKCYAgCQgmIKAEAKjouCipF29EyEPZtZu6OkTjnllGr2hS98oeX65ptvXr3mgQce6Ptgo8hI27P2K2OZ46IAAEhNMQUAIAXFFACAFBRTAABSUEwBAEhBMQUAIAXHRUHFSDt6JsKeZWwbaXvWfmUsc1wUAACpKaYAAKSgmAIAkIJiCgBACoopAAApKKYAAKTQ9rgoAAAYKp6YAgCQgmIKAEAKiikAACkopgAApKCYAgCQgmKaVCnl4lLKN4d7DqBv7FkYOezXvBTTPiqlPFZKeaaUsnyvtc+VUqYP41gdK6XsV0q5vZTy8kh9H6CdUbhnzyilzCqlzCulzCylfG24Z4LBMgr3q8+xHVJM+2d8RBw93EP0VyllqRbLz0fEdyPiW0M7DQyp0bRn/zkiNm6aZlJEbBcRnyql7DO0k8ESNZr2q8+xHVJM++fMiDimlLLSG4NSyrqllKaUMr7X2vRSyud63j6klHJbKeWfSilzSymPllK261mfVUp5tpTy6TfcdtVSyg2llPmllJtLKev0uvfGPdnzpZQHSyn79couLqVcUEq5ppSyICI+8MZ5m6b596ZpfhkRTw34XwXyGk179sGmaRb0WuqKiA07/YeBhEbTfvU5tkOKaf/cExHTI+KYDq/fJiL+OyLeEhE/j4jLImLr6P7kMjUiziulrNDr7x8YESdHxKoRcV9EXBoR0fOtjht67vHWiDggIs4vpWza69pPRcQpEbFiRPy2w3lhpBtVe7aU8pVSyksR8URELN9zPxgtRtV+pTOKaf+dEBFHllJW6+DaPzVNc1HTNK9FxC8i4u0R8Y2maRY1TXN9RCyO//kE5DdN09zSNM2iiPhaRGxbSnl7ROwZEY/13OvVpmnujYjLI+Ljva6d1jTNbU3TdDVN80oHs8JoMWr2bNM034ruT4RbRsQlEfFiB+8TZDZq9iudUUz7qWma+yPi6oj4SgeXP9Pr7YU993vjWu+v5mb1et2XovtnVtaKiHUiYpueb1fMLaXMje6v/NZodS2MZaNtzzbd/qvntU/q+7sC+Y22/Ur/jX/zv0ILX4+IeyPirF5rr//s13IRMa/n7d7/I+7E219/o+fbD6tE98+rzIqIm5um2aXNtc0AXxtGk9G4Z8dHxAb9nhDyG437lT7yxLQDTdM8HN3fJjiq19qciHgyIqaWUpYqpXwmBv5JY49SyvallKWj++dg7mqaZlZ0fzW5USnloFLKhJ4/W5dSNunrjXtmXCa6P7mNK6UsU0qZMMB5IaWRvmdLKeNKKZ8vpaxcur0vIv42Iv5jgPNCOiN9v0b4HDsQimnnvhHdv3zQ22ERcWxEPBcRm0bE7QN8jZ9H91eOz0fEVtH9rYRommZ+ROwaEftH91d3T0fE6RExsR/3Pii6v61xQUT8r563fzTAeSGzkb5nPxYRj0TE/Ij4WUSc2/MHRqORvl99ju1QaRpPowEAGH6emAIAkIJiCgBACoopAAApKKYAAKSgmAIAkIJiCgBACoopAAApKKYAAKSgmAIAkIJiCgBACoopAAApKKYAAKSgmAIAkIJiCgBACoopAAApKKYAAKSgmAIAkIJiCgBACoopAAApKKYAAKSgmAIAkML4dmEppRmqQSCbpmnKcM/QX/YsY9lI27P2K2NZbb96YgoAQAqKKQAAKSimAACkoJgCAJCCYgoAQAqKKQAAKSimAACkoJgCAJCCYgoAQAqKKQAAKSimAACkoJgCAJCCYgoAQAqKKQAAKSimAACkoJgCAJCCYgoAQAqKKQAAKSimAACkoJgCAJCCYgoAQAqKKQAAKSimAACkoJgCAJCCYgoAQAqKKQAAKYwf7gHGgrPPPruaHXXUUdXs/vvvr2Z77rlnNZs5c2bfBgMASMQTUwAAUlBMAQBIQTEFACAFxRQAgBQUUwAAUlBMAQBIwXFRg2TdddetZlOnTq1mXV1d1WyTTTapZhtvvHE1c1wUvLmNNtqomk2YMKGa7bDDDtXs/PPPr2bt9vpQmjZtWjXbf//9q9nixYuXxDgwYO3263bbbVfNTj311Gr213/91wOaic55YgoAQAqKKQAAKSimAACkoJgCAJCCYgoAQAqKKQAAKTguapDMmTOnmt1yyy3VbK+99loS48CYsummm1azQw45pOX6Jz7xieo148bVv2Zfa621qlm7I6GapqlmQ6ndx5zvf//71exLX/pSNZs3b95ARoIBmTx5cjW76aabqtnTTz9dzdZYY42OrmPgPDEFACAFxRQAgBQUUwAAUlBMAQBIQTEFACAFxRQAgBQcFzVIFixYUM1mzpw5hJPA2HPaaadVsz322GMIJxnZDj744Gr2z//8z9XstttuWxLjwBLV7kgox0UNH09MAQBIQTEFACAFxRQAgBQUUwAAUlBMAQBIQTEFACAFx0UNkpVWWqmabbHFFkM3CIxBN9xwQzXr5LioZ599tpq1OzZp3Lj61/pdXV39niMiYrvttmu5vuOOO3Z0P6BbKWW4R6AFT0wBAEhBMQUAIAXFFACAFBRTAABSUEwBAEhBMQUAIAXHRQ2S5ZZbrpqtvfbag/56W2+9dTWbMWNGNZs5c+agzwLD7YILLqhmV155Zb/v95e//KWaPf300/2+30BMmjSp5fr9999fvWattdbq6LXa/Vvdc889Hd0Tsmqappots8wyQzgJvXliCgBACoopAAApKKYAAKSgmAIAkIJiCgBACoopAAApOC5qkDz11FPV7OKLL65mJ554Ykev1+66uXPnVrPzzjuvo9eDzF599dVqNmvWrCGcZPDttttuLddXXnnlQX+tJ554opotWrRo0F8Psnrve99bze68884hnGTs8cQUAIAUFFMAAFJQTAEASEExBQAgBcUUAIAU/Fb+EDj55JOrWae/lQ+MHvvvv381O+yww1quL7vssoM+xwknnDDo94Qlrd2pHC+++GI1mzx5cjXbYIMNBjQTnfPEFACAFBRTAABSUEwBAEhBMQUAIAXFFACAFBRTAABScFzUMBs3rv61QVdX1xBOAgzUgQceWM2+8pWvVLMNN9ywmk2YMGFAM73RfffdV83+8pe/DOprwVCYO3duNbv11lur2Z577rkEpmGgPDEFACAFxRQAgBQUUwAAUlBMAQBIQTEFACAFxRQAgBQcFzXM2h0J1TTNEE4CI9e6665bzQ466KCW6zvvvPOgz7H99ttXs8Hez/Pmzatm7Y6muuaaa6rZwoULBzQTwEB5YgoAQAqKKQAAKSimAACkoJgCAJCCYgoAQAqKKQAAKTguChgRNttss2p21VVXVbO11157SYwz7G699dZq9sMf/nAIJ4HR5y1vectwjzBmeWIKAEAKiikAACkopgAApKCYAgCQgmIKAEAKiikAACk4LgoY8UopHWWDbdy4+tf6XV1dg/pae+65ZzX70Ic+VM2uvfbaQZ0DRqO99tpruEcYszwxBQAgBcUUAIAUFFMAAFJQTAEASEExBQAgBcUUAIAUHBc1zJbE8TI77LBDNTvvvPM6uicMt/vvv7+a7bTTTtVs6tSpLdevu+666jWvvPJKn+caDJ/97Ger2ZFHHjmEk8DoctNNN1WzdkeuMXw8MQUAIAXFFACAFBRTAABSUEwBAEhBMQUAIAXFFACAFErTNPWwlHrIoHjttdeqWbv/bjq1+eabt1x/4IEHBv21Rrqmacpwz9Bf9uzINHny5Gr23HPP9ft+H/nIR6rZtdde2+/7jRQjbc/ar0vevvvuW83+9V//tZotXLiwmr3rXe+qZjNnzuzbYFT3qyemAACkoJgCAJCCYgoAQAqKKQAAKSimAACkMH64Bxjrvv/971ezz3/+84P+eocffnjL9S996UuD/lpA3+y2227DPQKMSq+++mpH15VSP+Bh4sSJnY5DH3hiCgBACoopAAApKKYAAKSgmAIAkIJiCgBACoopAAApOC5qmM2YMWO4R4AhNWHChGq26667VrMbb7yxmi1cuHBAMw2FQw89tJqdffbZQzgJjB3Tpk2rZu0+/2688cbVrN3xikcccUSf5qLOE1MAAFJQTAEASEExBQAgBcUUAIAUFFMAAFJQTAEASKE0TVMPS6mHLHF//OMfq9kGG2zQ0T3HjWv9tciGG25YveaRRx7p6LVGuqZpynDP0F9Z9uz2229fzb72ta9Vs1122aWarbfeetVs1qxZfRtsEKyyyirVbI899qhm5557bjVbccUV+z1HuyOy9tprr2p200039fu1RoqRtmez7Nex6rvf/W41a3e82+qrr17NXnnllYGMNKbU9qsnpgAApKCYAgCQgmIKAEAKiikAACkopgAApKCYAgCQwvjhHoC63//+99Vs/fXX7+ieXV1dnY4DfXbeeedVs80226yje/793/99NZs/f35H9+xEuyOtttxyy2rW7mi+dqZPn95y/YILLqheM5qPhIKh0G6/Ll68eAgnGXs8MQUAIAXFFACAFBRTAABSUEwBAEhBMQUAIAXFFACAFBwXldgPf/jDavaRj3xkCCeB4feFL3xhuEcYkGeffbaa/frXv65mRx99dMv1V155ZcAzAa1NmjSpmn30ox+tZldcccWSGGdM8cQUAIAUFFMAAFJQTAEASEExBQAgBcUUAIAUFFMAAFJwXFRiDzzwQDX7wx/+UM022WSTJTEO9NkhhxxSzY488shq9ulPf3oJTNN/jzzySDV7+eWXq9mtt95azdod/3b//ff3bTBg0Oy3337VbNGiRdWs3edfBs4TUwAAUlBMAQBIQTEFACAFxRQAgBQUUwAAUlBMAQBIoTRNUw9LqYcwyjVNU4Z7hv4aCXt24sSJ1azdMVPf/OY3q9nKK69cza688sqW6zfccEP1mmnTplWzp59+upoxvEbanh0J+3U0u+yyy6pZu2MX99prr2o2c+bMAc00ltT2qyemAACkoJgCAJCCYgoAQAqKKQAAKSimAACk4LfyoWKk/YZvhD3L2DbS9qz9yljmt/IBAEhNMQUAIAXFFACAFBRTAABSUEwBAEhBMQUAIAXFFACAFBRTAABSUEwBAEhBMQUAIAXFFACAFBRTAABSUEwBAEhBMQUAIAXFFACAFBRTAABSUEwBAEhBMQUAIAXFFACAFBRTAABSKE3TDPcMAADgiSkAADkopgAApKCYAgCQgmIKAEAKiikAACkopgAApPD/AwBnJTwFC47DAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x864 with 9 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# create figure with 3x3 subplots using matplotlib.pyplot\n",
    "fig, axs = plt.subplots(3, 3, figsize = (12, 12))\n",
    "plt.gray()\n",
    "\n",
    "# loop through subplots and add mnist images\n",
    "for i, ax in enumerate(axs.flat):\n",
    "    ax.matshow(X_train[i])\n",
    "    ax.axis('off')\n",
    "    ax.set_title('Number {}'.format(y_train[i]))\n",
    "    \n",
    "# display the figure\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(60000, 784)\n",
      "(784,)\n",
      "(10000, 784)\n",
      "(784,)\n",
      "(60000,)\n",
      "(10000,)\n",
      "Y Class Labels: {0, 1, 2, 3, 4, 5, 6, 7, 8, 9}\n"
     ]
    }
   ],
   "source": [
    "# preprocessing the images\n",
    "\n",
    "# convert each image to 1 dimensional array\n",
    "X_train = X_train.reshape(len(X_train),-1)\n",
    "X_test = X_test.reshape(len(X_test),-1)\n",
    "# y_test = y_test\n",
    "\n",
    "# normalize the data to 0 - 1\n",
    "X_train = X_train.astype(float) / 255.\n",
    "X_test = X_test.astype(float) / 255.\n",
    "\n",
    "y_classes = set(y_train)\n",
    "\n",
    "print(X_train.shape)\n",
    "print(X_train[0].shape)\n",
    "print(X_test.shape)\n",
    "print(X_test[0].shape)\n",
    "print(y_train.shape)\n",
    "print(y_test.shape)\n",
    "print(f'Y Class Labels: {y_classes}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# using an optimizer\n",
    "def negloglike(theta, x):\n",
    "    mu = theta[0]\n",
    "    var = theta[1]\n",
    "    \n",
    "    # setup lower bound for the variance\n",
    "    limit = 1e-6    \n",
    "    if var < limit:\n",
    "        var = limit\n",
    "        \n",
    "    N = x.size\n",
    "    nll = N*log(var) + np.sum((x - mu)**2/var)\n",
    "    return nll\n",
    "\n",
    "# fit to create a normal probability distribution for a feature data (one column)\n",
    "def fit_distribution(x):\n",
    "    feature_mean = np.mean(x)\n",
    "    feature_var = np.var(x)\n",
    "    \n",
    "    bnds = ((None, None), (0, None))\n",
    "    \n",
    "    params = optimize.minimize(negloglike, np.array([feature_mean, feature_var]), \n",
    "                               args=(x), bounds=bnds, tol=1e-10, method='BFGS')\n",
    "    \n",
    "    opt_mean = params.x[0]\n",
    "    opt_std = np.sqrt(params.x[1]) # convert variance to std dev\n",
    "    \n",
    "    dist = norm(opt_mean, opt_std)\n",
    "    \n",
    "    return dist\n",
    "\n",
    "# this function calculates the mle = likelihood\n",
    "def mle_prob(X, col_dist_list):\n",
    "    # feature_1 * feature_2 * .. * feature_N\n",
    "#     print(X.shape)\n",
    "    acc_prob = 1\n",
    "    for i in range(X.shape[0]):\n",
    "        dist = col_dist_list[i]\n",
    "        prob = dist.pdf(X[i])\n",
    "        if np.isnan(prob):\n",
    "            prob = 1\n",
    "        acc_prob *= prob\n",
    "    return acc_prob\n",
    "\n",
    "# this function calculates the posterior = likelihood * prior\n",
    "def map_prob(X, col_dist_list, prior):\n",
    "    # feature_1 * feature_2 * .. * feature_N * prior\n",
    "    acc_prob = 1\n",
    "    for i in range(X.shape[0]):\n",
    "        dist = col_dist_list[i]\n",
    "        prob = dist.pdf(X[i])\n",
    "        if np.isnan(prob):\n",
    "            prob = 1\n",
    "        acc_prob *= prob\n",
    "    return acc_prob * prior"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Processing Y Class: 0\n",
      "Processing Y Class: 1\n",
      "Processing Y Class: 2\n",
      "Processing Y Class: 3\n",
      "Processing Y Class: 4\n",
      "Processing Y Class: 5\n",
      "Processing Y Class: 6\n",
      "Processing Y Class: 7\n",
      "Processing Y Class: 8\n",
      "Processing Y Class: 9\n",
      "(10,)\n"
     ]
    }
   ],
   "source": [
    "X_subsets = []\n",
    "for cls in y_classes:\n",
    "    print(f'Processing Y Class: {cls}')\n",
    "#     this_x = X[y==cls]\n",
    "    X_subsets.append(X_train[y_train==cls])\n",
    "X_subsets = np.asarray(X_subsets)\n",
    "print(X_subsets.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Processing Y Class: 0\n",
      "(5923, 784)\n",
      "Processing Y Class: 1\n",
      "(6742, 784)\n",
      "Processing Y Class: 2\n",
      "(5958, 784)\n",
      "Processing Y Class: 3\n",
      "(6131, 784)\n",
      "Processing Y Class: 4\n",
      "(5842, 784)\n",
      "Processing Y Class: 5\n",
      "(5421, 784)\n",
      "Processing Y Class: 6\n",
      "(5918, 784)\n",
      "Processing Y Class: 7\n",
      "(6265, 784)\n",
      "Processing Y Class: 8\n",
      "(5851, 784)\n",
      "Processing Y Class: 9\n",
      "(5949, 784)\n",
      "Runtime: 0:01:27.098300\n"
     ]
    }
   ],
   "source": [
    "start_time = datetime.now()\n",
    "\n",
    "cls_col_dist_list = []\n",
    "for cls in y_classes:\n",
    "    # OK, since y classes are 0 and 1, so use it as the index\n",
    "    print(f'Processing Y Class: {cls}')\n",
    "    x_subset = X_subsets[cls]\n",
    "    print(x_subset.shape)\n",
    "    \n",
    "    # For each column, create a normal distribution and put in the list\n",
    "    col_dist_list = np.array([fit_distribution(x_subset[:, i]) for i in range(x_subset.shape[1])])\n",
    "    cls_col_dist_list.append(col_dist_list)\n",
    "cls_col_dist_list = np.array(cls_col_dist_list)\n",
    "cls_col_dist_list.shape\n",
    "\n",
    "end_time = datetime.now()\n",
    "\n",
    "print(f'Runtime: {str(end_time-start_time)}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.09871667, 0.11236667, 0.0993    , 0.10218333, 0.09736667,\n",
       "       0.09035   , 0.09863333, 0.10441667, 0.09751667, 0.09915   ])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "prior_list = np.array([len(X_subsets[cls]) / len(X_train) for cls in y_classes])\n",
    "prior_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# try to classify sample\n",
    "y_pred = []\n",
    "y_test_list = []\n",
    "# for i in range(len(X)):\n",
    "for i in range(500):\n",
    "    sample_index = i\n",
    "#     x_t, y_t = X[i], y[i]\n",
    "    x_t, y_t = X_test[i], y_test[i]\n",
    "#     print(X_test, y_test)\n",
    "    \n",
    "    # calculate the probability of p(y=0|x)\n",
    "    # for each class, calculate the probability\n",
    "    post_prob = [mle_prob(x_t, cls_col_dist_list[cls]) for cls in y_classes]\n",
    "#     print(post_prob)\n",
    "\n",
    "    y_pred.append(np.argmax(post_prob))\n",
    "    y_test_list.append(y_t)    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[7 2 1 0 4 1 4 9 4 7 0 6 9 0 1 3 9 7 2 4 9 6 6 5 4 0 7 4 0 1 3 1 3 0 7 2 7\n",
      " 1 2 1 1 7 4 2 3 3 3 6 4 4 6 3 5 3 6 0 4 1 9 4 7 2 7 2 7 1 2 4 3 0 7 0 2 7\n",
      " 1 7 3 2 1 7 7 6 2 7 4 4 7 3 6 1 3 6 1 3 1 4 1 2 6 9 6 0 5 4 9 9 2 1 9 4 1\n",
      " 1 1 9 7 4 4 4 7 4 3 4 7 6 7 9 0 5 8 5 6 6 3 2 1 1 0 1 6 4 6 7 3 1 7 1 8 2\n",
      " 0 2 4 4 5 3 1 5 6 0 3 4 4 6 5 4 6 3 4 5 1 4 4 7 2 3 2 1 1 4 1 8 1 8 1 0 3\n",
      " 4 2 3 0 1 1 1 0 4 0 4 1 6 4 2 3 6 1 1 1 3 9 3 2 9 4 7 9 1 9 0 3 5 7 5 7 2\n",
      " 2 7 1 2 8 1 1 7 3 1 8 7 7 7 2 2 4 1 3 3 8 4 2 3 0 2 4 1 4 1 9 3 7 7 1 1 2\n",
      " 6 8 5 7 7 9 1 3 1 2 0 3 0 1 9 9 4 1 8 2 1 2 9 1 5 9 2 6 0 1 3 4 2 9 2 0 4\n",
      " 0 0 2 1 4 4 1 2 4 0 2 7 4 3 3 0 0 3 1 9 6 5 2 5 1 7 9 1 3 4 2 0 7 1 1 2 1\n",
      " 5 3 3 4 7 2 6 3 2 1 3 5 1 3 5 1 7 1 3 0 6 1 8 5 1 7 4 0 6 7 2 4 0 6 5 6 3\n",
      " 7 2 2 8 4 5 9 1 1 4 0 7 3 7 6 1 6 2 1 9 2 8 6 1 7 3 2 3 4 4 2 3 3 4 3 4 0\n",
      " 0 3 1 7 1 3 7 9 7 1 7 2 1 4 0 9 2 0 2 9 1 4 4 1 4 4 5 9 8 4 3 7 6 0 0 3 4\n",
      " 8 0 6 4 3 3 3 3 2 3 9 1 2 3 9 0 7 6 6 6 7 8 8 2 7 5 8 7 6 1 8 4 1 2 3 3 1\n",
      " 9 7 0 4 0 4 9 7 1 0 3 2 3 7 2 9 1 4 1]\n",
      "[7 2 1 0 4 1 4 9 5 9 0 6 9 0 1 5 9 7 3 4 9 6 6 5 4 0 7 4 0 1 3 1 3 4 7 2 7\n",
      " 1 2 1 1 7 4 2 3 5 1 2 4 4 6 3 5 5 6 0 4 1 9 5 7 8 9 3 7 4 6 4 3 0 7 0 2 9\n",
      " 1 7 3 2 9 7 7 6 2 7 8 4 7 3 6 1 3 6 9 3 1 4 1 7 6 9 6 0 5 4 9 9 2 1 9 4 8\n",
      " 7 3 9 7 4 4 4 9 2 5 4 7 6 7 9 0 5 8 5 6 6 5 7 8 1 0 1 6 4 6 7 3 1 7 1 8 2\n",
      " 0 2 9 9 5 5 1 5 6 0 3 4 4 6 5 4 6 5 4 5 1 4 4 7 2 3 2 7 1 8 1 8 1 8 5 0 8\n",
      " 9 2 5 0 1 1 1 0 9 0 3 1 6 4 2 3 6 1 1 1 3 9 5 2 9 4 5 9 3 9 0 3 6 5 5 7 2\n",
      " 2 7 1 2 8 4 1 7 3 3 8 8 7 9 2 2 4 1 5 9 8 7 2 3 0 4 4 2 4 1 9 5 7 7 2 8 2\n",
      " 6 8 5 7 7 9 1 8 1 8 0 3 0 1 9 9 4 1 8 2 1 2 9 7 5 9 2 6 4 1 5 8 2 9 2 0 4\n",
      " 0 0 2 8 4 7 1 2 4 0 2 7 4 3 3 0 0 3 1 9 6 5 2 5 9 2 9 3 0 4 2 0 7 1 1 2 1\n",
      " 5 3 3 9 7 8 6 5 6 1 3 8 1 0 5 1 3 1 5 5 6 1 8 5 1 7 9 4 6 2 2 5 0 6 5 6 3\n",
      " 7 2 0 8 8 5 4 1 1 4 0 3 3 7 6 1 6 2 1 9 2 8 6 1 9 5 2 5 4 4 2 8 3 8 2 4 5\n",
      " 0 3 1 7 7 5 7 9 7 1 9 2 1 4 2 9 2 0 4 9 1 4 8 1 8 4 5 9 8 8 3 7 6 0 0 3 0\n",
      " 2 6 6 4 9 3 3 3 2 3 9 1 2 6 8 0 5 6 6 6 3 8 8 2 7 5 8 9 6 1 8 4 1 2 5 9 1\n",
      " 9 7 5 4 0 8 9 9 1 0 5 2 3 7 8 9 4 0 6]\n"
     ]
    }
   ],
   "source": [
    "y_pred = np.asarray(y_pred)\n",
    "y_test_list = np.asarray(y_test_list)\n",
    "\n",
    "print(y_pred)\n",
    "print(y_test_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.758"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum(y_pred==y_test_list)/len(y_test_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# try to classify sample\n",
    "y_pred = []\n",
    "y_test_list = []\n",
    "# for i in range(len(X)):\n",
    "for i in range(500):\n",
    "    sample_index = i\n",
    "#     x_t, y_t = X[i], y[i]\n",
    "    x_t, y_t = X_test[i], y_test[i]\n",
    "#     print(X_test, y_test)\n",
    "    \n",
    "    # calculate the probability of p(y=0|x)\n",
    "    # for each class, calculate the probability\n",
    "    post_prob = [map_prob(x_t, cls_col_dist_list[cls], prior_list[cls]) for cls in y_classes]\n",
    "# #     print(post_prob)\n",
    "\n",
    "    y_pred.append(np.argmax(post_prob))\n",
    "    y_test_list.append(y_t)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[7 2 1 0 4 1 4 9 4 7 0 6 9 0 1 3 9 7 2 4 9 6 6 5 4 0 7 4 0 1 3 1 3 0 7 2 7\n",
      " 1 2 1 1 7 4 2 3 3 3 6 4 4 6 3 5 3 6 0 4 1 9 4 7 2 7 2 7 1 2 4 3 0 7 0 2 7\n",
      " 1 7 3 2 1 7 7 6 2 7 4 4 7 3 6 1 3 6 1 3 1 4 1 2 6 9 6 0 5 4 9 9 2 1 9 4 1\n",
      " 1 1 9 7 4 4 4 7 4 3 4 7 6 7 9 0 5 8 5 6 6 3 2 1 1 0 1 6 4 6 7 3 1 7 1 8 2\n",
      " 0 2 4 4 5 3 1 5 6 0 3 4 4 6 5 4 6 3 4 5 1 4 4 7 2 3 2 1 1 4 1 8 1 8 1 0 3\n",
      " 4 2 3 0 1 1 1 0 4 0 4 1 6 4 2 3 6 1 1 1 3 9 3 2 9 4 7 9 1 9 0 3 5 7 5 7 2\n",
      " 2 7 1 2 8 1 1 7 3 1 8 7 7 7 2 2 4 1 3 3 8 4 2 3 0 2 4 1 4 1 9 3 7 7 1 1 2\n",
      " 6 8 5 7 7 9 1 3 1 2 0 3 0 1 9 9 4 1 8 2 1 2 9 1 5 9 2 6 0 1 3 4 2 9 2 0 4\n",
      " 0 0 2 1 4 4 1 2 4 0 2 7 4 3 3 0 0 3 1 9 6 5 2 5 1 7 9 1 3 4 2 0 7 1 1 2 1\n",
      " 5 3 3 4 7 2 6 3 2 1 3 5 1 3 5 1 7 1 3 0 6 1 8 5 1 7 4 0 6 7 2 4 0 6 5 6 3\n",
      " 7 2 2 8 4 5 9 1 1 4 0 7 3 7 6 1 6 2 1 9 2 8 6 1 7 3 2 3 4 4 2 3 3 4 3 4 0\n",
      " 0 3 1 7 1 3 7 9 7 1 7 2 1 4 0 9 2 0 2 9 1 4 4 1 4 4 5 9 8 4 3 7 6 0 0 3 4\n",
      " 8 0 6 4 3 3 3 3 2 3 9 1 2 3 9 0 7 6 6 6 7 8 8 2 7 5 8 7 6 1 8 4 1 2 3 3 1\n",
      " 9 7 0 4 0 4 9 7 1 0 3 2 3 7 2 9 1 4 1]\n",
      "[7 2 1 0 4 1 4 9 5 9 0 6 9 0 1 5 9 7 3 4 9 6 6 5 4 0 7 4 0 1 3 1 3 4 7 2 7\n",
      " 1 2 1 1 7 4 2 3 5 1 2 4 4 6 3 5 5 6 0 4 1 9 5 7 8 9 3 7 4 6 4 3 0 7 0 2 9\n",
      " 1 7 3 2 9 7 7 6 2 7 8 4 7 3 6 1 3 6 9 3 1 4 1 7 6 9 6 0 5 4 9 9 2 1 9 4 8\n",
      " 7 3 9 7 4 4 4 9 2 5 4 7 6 7 9 0 5 8 5 6 6 5 7 8 1 0 1 6 4 6 7 3 1 7 1 8 2\n",
      " 0 2 9 9 5 5 1 5 6 0 3 4 4 6 5 4 6 5 4 5 1 4 4 7 2 3 2 7 1 8 1 8 1 8 5 0 8\n",
      " 9 2 5 0 1 1 1 0 9 0 3 1 6 4 2 3 6 1 1 1 3 9 5 2 9 4 5 9 3 9 0 3 6 5 5 7 2\n",
      " 2 7 1 2 8 4 1 7 3 3 8 8 7 9 2 2 4 1 5 9 8 7 2 3 0 4 4 2 4 1 9 5 7 7 2 8 2\n",
      " 6 8 5 7 7 9 1 8 1 8 0 3 0 1 9 9 4 1 8 2 1 2 9 7 5 9 2 6 4 1 5 8 2 9 2 0 4\n",
      " 0 0 2 8 4 7 1 2 4 0 2 7 4 3 3 0 0 3 1 9 6 5 2 5 9 2 9 3 0 4 2 0 7 1 1 2 1\n",
      " 5 3 3 9 7 8 6 5 6 1 3 8 1 0 5 1 3 1 5 5 6 1 8 5 1 7 9 4 6 2 2 5 0 6 5 6 3\n",
      " 7 2 0 8 8 5 4 1 1 4 0 3 3 7 6 1 6 2 1 9 2 8 6 1 9 5 2 5 4 4 2 8 3 8 2 4 5\n",
      " 0 3 1 7 7 5 7 9 7 1 9 2 1 4 2 9 2 0 4 9 1 4 8 1 8 4 5 9 8 8 3 7 6 0 0 3 0\n",
      " 2 6 6 4 9 3 3 3 2 3 9 1 2 6 8 0 5 6 6 6 3 8 8 2 7 5 8 9 6 1 8 4 1 2 5 9 1\n",
      " 9 7 5 4 0 8 9 9 1 0 5 2 3 7 8 9 4 0 6]\n"
     ]
    }
   ],
   "source": [
    "y_pred = np.asarray(y_pred)\n",
    "y_test_list = np.asarray(y_test_list)\n",
    "\n",
    "print(y_pred)\n",
    "print(y_test_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.758"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum(y_pred==y_test_list)/len(y_test_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}