hw_1/ml2019fall_hw1/bayes_decision_rule/run.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bayes Decision Rule\n",
    "*Complete and hand in this completed worksheet (including its outputs and any supporting code outside of the worksheet) with your assignment submission. Please check the pdf file for more details.*\n",
    "\n",
    "In this exercise you will:\n",
    "    \n",
    "- implement the calculation of **likelihood** of each features given particular class \n",
    "- implement the calculation of **posterior** of each class given particular feature\n",
    "- implement the calculation of **minimal total risk** of bayes decision rule"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# some basic imports\n",
    "import scipy.io as sio\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "data = sio.loadmat('data.mat')\n",
    "x1_train, x1_test, x2_train, x2_test = data['x1_train'], data['x1_test'], data['x2_train'], data['x2_test']\n",
    "all_x = np.concatenate([x1_train, x1_test, x2_train, x2_test], 1)\n",
    "data_range = [np.min(all_x), np.max(all_x)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "from get_x_distribution import get_x_distribution\n",
    "\n",
    "train_x = get_x_distribution(x1_train, x2_train, data_range)\n",
    "test_x = get_x_distribution(x1_test, x2_test, data_range)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 1 likelihood\n",
    "Your code for this section will be written inside **bayes_decision_rule/likehood.py**. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEKCAYAAADn+anLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAFH9JREFUeJzt3X+w3XWd3/Hnyxg2LTCuhqiUGzeBZdalqIh30f6yriINuE3olE4j3S5W1gxTU3a7ZWbDZKUMakGd7tYqwmYXRvaHoquue1fiAK46nbYDTZAETZAlZEBuoUsMbJCuKJF3/zgn9OTmfHPvTc75nnPD8zFzJ+f7PZ/7Pe98zz33db+f7/fz+aaqkCSpn5eMugBJ0vgyJCRJjQwJSVIjQ0KS1MiQkCQ1MiQkSY1aDYkkq5I8kGRXkg19nn9Pkj1JtnW/frXN+iRJB3tpWy+UZBFwPfBOYBrYkmSqqnbOaPq5qlrfVl2SpGZtHkmcA+yqqt1V9WPgVmBNi68vSZqn1o4kgFOAR3uWp4E392n3z5O8FfhL4N9X1aMzGyRZB6wDOP7449/02te+dgjlStKx65577vl+VS2brV2bIZE+62bOCfLnwGer6kdJLgNuAd5+yDdVbQI2AUxOTtbWrVsHXaskHdOSPDKXdm12N00Dy3uWJ4DHehtU1d6q+lF38feAN7VUmySpjzZDYgtwepKVSY4D1gJTvQ2SnNyzuBq4v8X6JEkztNbdVFX7k6wHbgcWATdX1Y4k1wBbq2oKuDzJamA/8CTwnrbqkyQdKgt9qnDPSUgalueee47p6WmeffbZUZdyxJYsWcLExASLFy8+aH2Se6pqcrbvb/PEtSQtKNPT05x44omsWLGCpN+1N+Otqti7dy/T09OsXLnyiLbhtByS1ODZZ59l6dKlCzIgAJKwdOnSozoSMiQk6TAWakAccLT1GxKSpEaek5CkOVqx4baBbu/h69410O0Ng0cSkqRGhoQkjbmpqSkuuuiig9bdcMMNXH755UN/bUNCksbcxo0bufrqqw9ad9ppp7Fz58w7LQyeISFJY2z79u08//zznHnmmTzyyCPccMMNQGegXxJ2797NpZdeesiRxqAYEpI0xrZt28ab3tSZ6/TOO+/kwQcfBGDnzp284Q1v4NRTT+Wmm24a2usbEpI0xp5//nmeeeYZfvKTn/ClL32JH/zgB/zwhz/k05/+NBdffPHQX99LYCVpjkZxyeoFF1zAxz/+cc466yw+/OEPc9111zE5Ocm6des4++yzh/76hoQkjbFXvepVbNu27YXl1atXH/T83r172bhxI/feey/XXnstV1555UBf35CQpAVs6dKl3HjjjUPbvuckJEmNDAlJUiNDQpLUyJCQJDUyJCRJjQwJSVIjQ0KS1MhxEpI0V1e/bMDb2zfY7Q2BRxKSpEaGhCSNOW86JElq5E2HJEl9zXbToS9/+cu8733vY82aNdxxxx0Df31PXEvSGJvtpkMXXnghF154IU899RRXXHEF55133kBf3yMJSRpjc73p0Ic+9CHe//73D/z1PZKQpLkawSWrs910qKrYsGED559//lBuQmRISNIYm+2mQ5/4xCf42te+xr59+9i1axeXXXbZQF/fkJCkBezyyy8f6qWwnpOQJDUyJCRJjQwJSVIjQ0KSDqOqRl3CUTna+lsNiSSrkjyQZFeSDYdpd1GSSjLZZn2S1GvJkiXs3bt3wQZFVbF3716WLFlyxNto7eqmJIuA64F3AtPAliRTVbVzRrsTgcuBu9uqTZL6mZiYYHp6mj179oy6lCO2ZMkSJiYmjvj727wE9hxgV1XtBkhyK7AGmDlD1QeBjwJXtFibJB1i8eLFrFy5ctRljFSb3U2nAI/2LE93170gyRuB5VX1lcNtKMm6JFuTbF3ICS9J467NkEifdS909CV5CfA7wH+YbUNVtamqJqtqctmyZQMsUZLUq82QmAaW9yxPAI/1LJ8InAl8M8nDwFuAKU9eS9LotBkSW4DTk6xMchywFpg68GRV7auqk6pqRVWtAO4CVlfV1hZrlCT1aC0kqmo/sB64Hbgf+HxV7UhyTZLVh/9uSdIotDrBX1VtBjbPWHdVQ9u3tVGTJKmZI64lSY0MCUlSI0NCktTIkJAkNTIkJEmNDAlJUiNDQpLUyJCQJDUyJCRJjQwJSVIjQ0KS1MiQkCQ1MiQkSY0MCUlSI0NCktTIkJAkNTIkJEmNDAlJUiNDQpLUyJCQJDUyJCRJjQwJSVIjQ0KS1Oiloy5A0nhaseG2vusfvu5dLVeiUfJIQpLUyJCQJDUyJCRJjQwJSVIjQ0KS1MirmyQNTb8rpLw6amHxSEKS1MiQkCQ1MiQkSY0MCUlSo1ZDIsmqJA8k2ZVkQ5/nL0vy7STbkvz3JGe0WZ8k6WCthUSSRcD1wPnAGcC7+4TAZ6rqdVV1FvBR4Lfbqk+SdKg2jyTOAXZV1e6q+jFwK7Cmt0FVPd2zeDxQLdYnSZqhzXESpwCP9ixPA2+e2SjJ+4HfAI4D3t5vQ0nWAesAXvOa1wy8UElSR5tHEumz7pAjhaq6vqpOA34T+K1+G6qqTVU1WVWTy5YtG3CZkqQD2gyJaWB5z/IE8Nhh2t8KXDjUiiRJh9VmSGwBTk+yMslxwFpgqrdBktN7Ft8FPNhifZKkGVo7J1FV+5OsB24HFgE3V9WOJNcAW6tqClif5FzgOeAp4JK26pMkHarVCf6qajOweca6q3oe/1qb9UiSDs8R15KkRoaEJKnRvEMiyfHd0dOSpGPcrCGR5CVJLk5yW5IngO8CjyfZkeRjM65IkiQdQ+ZyJPEN4DTgSuDVVbW8ql4J/CPgLuC6JL88xBolSSMyl6ubzq2q52aurKongS8CX0yyeOCVSZJGbtaQOBAQST4O/DydqTS205mxdVtvG0nSsWU+4yTuB74CLKYz1fcfJbmxqj45lMokSSM355Coqht7Fjcn+SSdqTYMCUk6Rs17xHWSy4CfBU4Enp6luSRpATuSwXSb6XQ9TQDXDrYcSdI4mXNIJPl8kp+vqu9V1U3APwU+PLzSJEmjNp/upj8CPpckwD3ACcDzQ6lKkjQW5nPiegqYSvJ64Cw6RyGbD/9dkqSFbN4nrqvqPuC+IdQiSRozs4ZEktfMcVt/XVVe7SRJx5C5HEncQmeUdQ7TpoBPA38wgJokSWNiLtNy/GIbhUiSxs98LoH9L8MsRJI0fuYzmO6ZJH+e5HiAJOcl+R9DqkuSNAbmcwnsbyW5GPhmkh8B/xfYMLTKJEkjN+eQSPIO4H10wuFk4NKqemBYhUmSRm8+3U0bgQ9U1duAi+iMvn77UKqSJI2FuYyTSHW8EAhV9e0k5wNfAL5+oM0wC5UktW9O97hO8u/6DKrbS+f+1rcAlwy+NEnSqM3lnMQq4L3AZ5OcCjwFLAEWAXcAv3PgNqaSpGPLXAbTPQt8CvhUksXAScAPq+qvh12cJGm0Zu1uSnJJku8neRL4feAZA0KSXhzmck7iA8A7gdcC3wP+01ArkiSNjbmck3i6qu7tPv5AkruHWZAkaXzMJSROTrKOzn2tvwssHm5JkqRxMZeQ+I/A64F/BbwOOCHJZmA7cF9VfXaI9UmSRmguVzdt6l1OMkEnNF4HXAAYEpJ0jDqS25dOA9N4f2tJOubNZ+4mSdKLTKshkWRVkgeS7EpyyDTjSX4jyc4k9yX5iyQ/02Z9kqSDtRYSSRYB1wPnA2cA705yxoxm9wKTVfV6OpMHfrSt+iRJh2rzSOIcYFdV7a6qHwO3Amt6G1TVN6rqb7qLdwETLdYnSZqhzZA4BXi0Z3m6u67JpcBX+z2RZF2SrUm27tmzZ4AlSpJ6tRkS6bOu7z0okvwyMAl8rN/zVbWpqiaranLZsmUDLFGS1Gvel8AehWlgec/yBPDYzEZJzqVzF7x/XFU/aqk2SVIfbR5JbAFOT7IyyXHAWmCqt0GSNwK/C6yuqidarE2S1EdrIVFV+4H1wO105oH6fFXtSHJNktXdZh8DTgD+JMm2JFMNm5MktaDN7iaqajMzRmpX1VU9j89tsx5J0uE54lqS1MiQkCQ1MiQkSY0MCUlSI0NCktTIkJAkNTIkJEmNDAlJUiNDQpLUyJCQJDUyJCRJjQwJSVIjQ0KS1KjVWWAltWPFhtv6rn/4une1XIkWOkNC0tgy7EbP7iZJUiNDQpLUyJCQJDUyJCRJjQwJSVIjQ0KS1MiQkCQ1MiQkSY0MCUlSI0dcSy8mV7+sYf2+duvQguGRhCSpkSEhSWpkd5O0UNhVpBEwJKS5ONpf0P6C1wJld5MkqZEhIUlqZEhIkhoZEpKkRoaEJKlRqyGRZFWSB5LsSrKhz/NvTfKtJPuTXNRmbZKkQ7UWEkkWAdcD5wNnAO9OcsaMZt8D3gN8pq26JEnN2hwncQ6wq6p2AyS5FVgD7DzQoKoe7j73fIt1SZIatNnddArwaM/ydHfdvCVZl2Rrkq179uwZSHGSpEO1GRLps66OZENVtamqJqtqctmyZUdZliSpSZvdTdPA8p7lCeCxFl9f0iD0m2LE6UWOWW0eSWwBTk+yMslxwFpgqsXXlyTNU2shUVX7gfXA7cD9wOerakeSa5KsBkjyC0mmgX8B/G6SHW3VJ0k6VKuzwFbVZmDzjHVX9TzeQqcbShocZ2CVjpgjriVJjQwJSVIjQ0KS1MiQkCQ1MiQkSY0MCUlSI0NCktSo1XES0rw5xkEaKUNCUrsM/gXF7iZJUiNDQpLUyJCQJDUyJCRJjTxxreHxBKW04HkkIUlqZEhIkhoZEpKkRoaEJKmRISFJamRISJIaGRKSpEaOk5C08DgGpzUeSUiSGhkSkqRGhoQkqZEhIUlqZEhIkhoZEpKkRoaEJKmR4ySkMbNiw2191z+8pOVCjnWOtZgTQ0KSjlS/oDnGQsaQUDP/0tIxwCOzo+M5CUlSI48kjlUeBUgaAI8kJEmNWg2JJKuSPJBkV5INfZ7/qSSf6z5/d5IVbdYnSTpYa91NSRYB1wPvBKaBLUmmqmpnT7NLgaeq6meTrAU+AvzLtmociEF089hVJL04LIDfF22ekzgH2FVVuwGS3AqsAXpDYg1wdffxF4BPJklVVWtV+gtaR6nf1TReSbOwvZivkEpbv3+TXASsqqpf7S7/a+DNVbW+p813um2mu8sPddt8f8a21gHruos/BzzQwn9hNicB35+11WgthBrBOgfNOgfrWKnzZ6pq2WwbafNIIn3WzUyoubShqjYBmwZR1KAk2VpVk6Ou43AWQo1gnYNmnYP1YquzzRPX08DynuUJ4LGmNkleCrwMeLKV6iRJh2gzJLYApydZmeQ4YC0wNaPNFHBJ9/FFwNdbPR8hSTpIa91NVbU/yXrgdmARcHNV7UhyDbC1qqaAm4A/TLKLzhHE2rbqG4Cx6v5qsBBqBOscNOscrBdVna2duJYkLTyOuJYkNTIkJEmNDIkjlOTqJP87ybbu1wUN7Q47FcmQa/xYku8muS/Jnyb56YZ2Dyf5dvf/sbXF+sZ+mpYky5N8I8n9SXYk+bU+bd6WZF/Pz8JVbdfZreOw72M6/mt3f96X5OwR1PhzPftpW5Knk/z6jDYj2Z9Jbk7yRHe81oF1r0hyZ5IHu/++vOF7L+m2eTDJJf3aDLHG4X7Oq8qvI/iiMzL8ilnaLAIeAk4FjgO2A2e0WON5wEu7jz8CfKSh3cPASS3vv1n3DfBvgRu7j9cCnxvB+3wycHb38YnAX/ap823AV9qubb7vI3AB8FU645HeAtw94noXAf+HzqCuke9P4K3A2cB3etZ9FNjQfbyh32cIeAWwu/vvy7uPX95ijUP9nHskMVwvTEVSVT8GDkxF0oqquqOq9ncX76IzNmVczGXfrAFu6T7+AvCOJP0GXA5NVT1eVd/qPv4BcD9wSps1DNAa4A+q4y7gp5OcPMJ63gE8VFWPjLCGF1TVf+PQcVm9P4O3ABf2+dZ/AtxZVU9W1VPAncCqtmoc9ufckDg667uHeDc3HIaeAjzaszzN6H7BvJfOX5H9FHBHknu6U560YS775oU23Q/BPmBpK9X10e3ueiNwd5+n/16S7Um+muTvtlrY/zfb+zhOP4/QOTr8bMNz47A/AV5VVY9D5w8G4JV92ozTfh3459ybDh1Gkq8Br+7z1EbgBuCDdHb8B4H/TOcNOmgTfb53oNccH67GqvqzbpuNwH7gjxs28w+q6rEkrwTuTPLd7l8swzSwaVrakOQE4IvAr1fV0zOe/hadLpNnuuemvgyc3naNzP4+jtP+PA5YDVzZ5+lx2Z9zNRb7dVifc0PiMKrq3Lm0S/J7wFf6PDWXqUiOymw1dk+i/RLwjup2TPbZxmPdf59I8qd0uoKGHRLzmaZlepTTtCRZTCcg/riqvjTz+d7QqKrNST6V5KSaMTHlsM3hfRz6z+M8nA98q6r+auYT47I/u/4qyclV9Xi3a+6JPm2m6ZxHOWAC+GYLtb1gmJ9zu5uO0Iy+3H8GfKdPs7lMRTI0SVYBvwmsrqq/aWhzfJITDzymcxKs3/9l0BbENC3dcyA3AfdX1W83tHn1gXMlSc6h87na216Vc34fp4Bf6V7l9BZg34GulBF4Nw1dTeOwP3v0/gxeAvxZnza3A+cleXm32/m87rpWDP1zPowz8C+GL+APgW8D99H5QTq5u/7vAJt72l1A54qYh+h0AbVZ4y46faXbul83zqyRztVF27tfO9qssd++Aa7p/rADLAH+pPv/+F/AqSN4n/8hna6D+3r24wXAZcBl3Tbru/tuO50Th39/BHX2fR9n1Bk6N/56qPuzO9l2nd06/jadX/ov61k38v1JJ7QeB56jc3RwKZ1zYH8BPNj99xXdtpPA7/d873u7P6e7gH/Tco1D/Zw7LYckqZHdTZKkRoaEJKmRISFJamRISJIaGRKSpEaGhCSpkSEhSWpkSEgDluQXuhM/LumOdN2R5MxR1yUdCQfTSUOQ5EN0Roz/LWC6qq4dcUnSETEkpCHozke1BXiWzrQSPxlxSdIRsbtJGo5XACfQuZvdkhHXIh0xjySkIUgyRedueyvpTP64fsQlSUfE+0lIA5bkV4D9VfWZJIuA/5nk7VX19VHXJs2XRxKSpEaek5AkNTIkJEmNDAlJUiNDQpLUyJCQJDUyJCRJjQwJSVKj/wcPE3j1iySAhwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from likelihood import likelihood\n",
    "\n",
    "l = likelihood(train_x)\n",
    "width = 0.35\n",
    "p1 = plt.bar(np.arange(data_range[0], data_range[1] + 1), l.T[:,0], width)\n",
    "p2 = plt.bar(np.arange(data_range[0], data_range[1] + 1) + width, l.T[:,1], width)\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('$P(x|\\omega)$')\n",
    "plt.legend((p1[0], p2[0]), ('$\\omega_1$', '$\\omega_2$'))\n",
    "plt.axis([data_range[0] - 1, data_range[1] + 1, 0, 0.5])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x1_test: err 3/all 100\n",
      "x2_test: err 61/all 200\n"
     ]
    }
   ],
   "source": [
    "#TODO\n",
    "#compute the number of all the misclassified x using maximum likelihood decision rule\n",
    "\n",
    "# begin answer\n",
    "# under this condition, solution of the equation is that p(w_i|x) = p(x|w_i)\n",
    "x1_test_idx = np.squeeze(x1_test - data_range[0])\n",
    "likelyhood_test_x1 = l[:, x1_test_idx]\n",
    "pred_x1 = np.argmax(likelyhood_test_x1, axis=0)\n",
    "err_x1 = np.sum(pred_x1)\n",
    "\n",
    "x2_test_idx = np.squeeze(x2_test - data_range[0])\n",
    "likelyhood_test_x2 = l[:, x2_test_idx]\n",
    "pred_x2 = 1-np.argmax(likelyhood_test_x2, axis=0)\n",
    "err_x2 = np.sum(pred_x2)\n",
    "\n",
    "print('x1_test: err {}/all {}'.format(err_x1, x1_test.shape[1]))\n",
    "print('x2_test: err {}/all {}'.format(err_x2, x2_test.shape[1]))\n",
    "# end answer\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 2 posterior\n",
    "Your code for this section will be written inside **bayes_decision_rule/posterior.py**. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEKCAYAAADn+anLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAFUtJREFUeJzt3X+wXOVdx/H3tyF4FZC2SdoiNzGhRivSUuiVolWL0mISnQRnGAVaizYlw9iITq1jmCjN0FZoO/6gLT+MgJRqoWgR75RUoLadOlYwoSQpJEUuEWQNmjSlabFNIeXrH3uCm5t97t172T27N7xfM3ey55xnz/nm7N188pwfz4nMRJKkdl7U7wIkSYPLkJAkFRkSkqQiQ0KSVGRISJKKDAlJUlFtIRERN0TEroh4oLD8LRGxtfr5UkScXFdtkqT26uxJ3AgsmWD5fwBvzMzXAO8F1tdRlCSp7Ii6NpSZX4yIhRMs/1LL5D3AcK9rkiRNrLaQmKKVwGdKCyNiFbAK4Kijjnrdq171qrrqkqTDwn333fe1zJw3WbuBC4mI+HmaIfEzpTaZuZ7qcNTIyEhu2rSppuok6fAQEY910m6gQiIiXgNcByzNzD39rkeSXugG5hLYiFgA3Ab8emb+e7/rkSTV2JOIiJuBM4C5EdEA3gPMBsjMa4FLgTnA1REBsD8zR+qqT5J0qDqvbjpvkuXvAN5RUzmSNKlnnnmGRqPBvn37+l3KtA0NDTE8PMzs2bOn9f6BOichSYOk0WhwzDHHsHDhQqojHDNKZrJnzx4ajQaLFi2a1joG5pyEJA2affv2MWfOnBkZEAARwZw5c55XT8iQkKQJzNSAOOD51m9ISJKKPCchSR1auOaOrq7v0St+qavr6wV7EpKkIkNCkgbc6Ogo55xzzkHzrrnmGi6++OKeb9uQkKQBt3btWtatW3fQvFe+8pVs27at59s2JCRpgG3ZsoVnn32Wk046iccee4xrrrkGaN7oFxHs2LGDlStXHtLT6BZDQpIG2ObNm3nd614HwN13383DDz8MwLZt2zj55JM54YQTuP7663u2fUNCkgbYs88+y1NPPcX3vvc9brvtNr71rW/xne98hxtvvJHzzz+/59v3ElhJ6lA/LlldtmwZV155Ja997Wt5//vfzxVXXMHIyAirVq3i1FNP7fn2DQlJGmAvf/nL2bx583PTy5cvP2j5nj17WLt2Lffffz+XX345l1xySVe3b0hI0gw2Z84crr322p6t33MSkqQiQ0KSVGRISJKKDAlJUpEhIUkqMiQkSUWGhCSpyPskJKlT647t8vr2dnd9PWBPQpJUZEhI0oDzoUOSpCIfOiRJamuyhw7dfvvtXHjhhaxYsYK77rqr69v3xLUkDbDJHjp09tlnc/bZZ/Pkk0/y7ne/m7POOqur26+tJxERN0TEroh4oLA8IuLDETEWEVsjovcDpUvSgOv0oUPve9/7eOc739n17dfZk7gR+ChwU2H5UmBx9fN64JrqT0kaDH24ZHWyhw5lJmvWrGHp0qU9eQhRbSGRmV+MiIUTNFkB3JSZCdwTES+OiOMy84laCpSkATTZQ4c+8pGP8NnPfpa9e/cyNjbGRRdd1NXtD9I5ieOBx1umG9U8Q0KSCi6++OKeXgo7SCERbeZl24YRq4BVAAsWLOhuFaU7KqfSzWy3juf7/sNpHYNQQz/WIY238/7283/olHrXMYFBCokGML9lehjY2a5hZq4H1gOMjIy0DRJp4DzfwJT6YJDukxgF3lZd5XQ6sNfzEZLUX7X1JCLiZuAMYG5ENID3ALMBMvNaYAOwDBgDvg38Zl21SVJJZhLR7mj4zNC8Fmj66ry66bxJlifQ/Yt8JWmahoaG2LNnD3PmzJmRQZGZ7Nmzh6GhoWmvY5DOSUjSQBkeHqbRaLB79+7ebOAbu9rP37u9a+sYGhpieHh4ioX9P0NCkgpmz57NokWLereBdacX5k/lyrsurGMCg3TiWpI0YAwJSVKRISFJKjIkJElFhoQkqcirmyRpul4AY3oZEtJM8QL4B0mDx8NNkqQiQ0KSVGRISJKKDAlJUpEhIUkqMiQkSUWGhCSpyJCQJBUZEpKkIkNCklRkSEiSigwJSVKRISFJKjIkJElFhoQkqciQkCQVGRKSpCJDQpJUZEhIkopqDYmIWBIRD0XEWESsabN8QUR8PiLuj4itEbGszvokSQerLSQiYhZwFbAUOBE4LyJOHNfsD4FbM/MU4Fzg6rrqkyQdqs6exGnAWGbuyMyngVuAFePaJPCD1etjgZ011idJGqfOkDgeeLxlulHNa7UOeGtENIANwG+3W1FErIqITRGxaffu3b2oVZJEvSERbebluOnzgBszcxhYBnw8Ig6pMTPXZ+ZIZo7MmzevB6VKkqDekGgA81umhzn0cNJK4FaAzPxXYAiYW0t1kqRD1BkSG4HFEbEoIo6keWJ6dFyb/wTOBIiIH6cZEh5PkqQ+qS0kMnM/sBq4E9hO8yqmByPisohYXjX7PeDCiNgC3Az8RmaOPyQlSarJEXVuLDM30Dwh3Trv0pbX24A31FmTJKnMO64lSUWGhCSpyJCQJBUZEpKkIkNCklRU69VNkjQw1h1bmL+33joGnD0JSVKRISFJKjIkJElFhoQkqciQkCQVeXWTpJ5ZuOaOQ+Y9esUv9aESTZc9CUlSkSEhSSoyJCRJRYaEJKloyiEREUdFxKxeFCNJGiyTXt0UES+i+TzqtwA/CXwX+L6I2E3zKXPrM/PhnlYp6QWp3dVR4BVSdeqkJ/F54JXAJcArMnN+Zr4M+FngHuCKiHhrD2uUJPVJJ/dJvCkznxk/MzO/DnwK+FREzO56ZZKkvpu0J3EgICLizyMiJmojSTq8TOXE9VPAaEQcBRARZ0XEv/SmLEnSIOh4WI7M/MOIOB/4QkR8F/hfYE3PKpMk9V3HIRERZwIX0gyH44CVmflQrwqTJPXfVAb4Wwtcmpn/HBGvBj4ZEe/KzM/1qDZJfeTlp4KpHW76hZbXX4mIpTSvbvrpXhQmqQd8rrOmqJOb6RZMsHhly/JvZOY3u1OWJGkQdNKT+BiQQLvLXw/MT+BG4KaJVhQRS4ArgVnAdZl5RZs2vwqsq9a5JTPP76BGSVIPTBoSmfnz3dhQNd7TVcCbgQawMSJGM3NbS5vFNO/sfkNmPhkRL+vGtiVJ01PnKLCnAWOZuSMznwZuAVaMa3MhcFVmPgmQmbtqrE+SNM5ULoG9EvhxqsNAwCcyc/MUtnU88HjLdAN4/bg2P1pt619oHpJal5n/2KaWVcAqgAULJjplIkl6PqbSk9gOfIjmOYVdwF9HxOopvL90TqPVEcBi4AzgPOC6iHjxIW/KXJ+ZI5k5Mm/evCmUIEmaiqlcAntty+SGiPgosBH4aIeraADzW6aHgZ1t2txTjQX1HxHxEM3Q2NhpnZKk7pnKzXQARMRFwI8AxwBTueR1I7A4IhYB/0XzGRXjr1y6nWYP4saImEvz8NOOqdYoSeqO6Zy43kDz0NMwcHmnb8rM/cBq4M7q/bdm5oMRcVlELK+a3QnsiYhtNJ9j8fuZuWcaNUqSumAqJ65vBd6TmduB6yPir4D7gU93uo7M3EAzZFrnXdryOoF3VT+SpD6byuGmv6Y5XlMA9wFHA8/2pCpJ0kDoZFiOyKZRms+TOBk4meahqg2tbXpbqiSpbh094zoifvvAGE2ZuSUzbwI+AZwUER8DLuhlkZKk/ujkcNMS4O3AzRFxAvAkMETzZre7gD+b4k11kqQZopOxm/YBVwNXR8RsYC7wncz8Rq+LkyT116SHmyLigoj4WkR8HbgOeMqAkKQXhk7OSfwRzZFbXwX8J/DHPa1IkjQwOjkn8c3MvL96/UcRcW8vC5IkDY5OQuK4atTV7cBXgdm9LUmSNCg6CYn3AK8B3gK8Gjg6IjbQHC58a2be3MP6JEl91MnVTetbpyNimGZovBpYBhgSknSYmvIosJnZoDmk94bJ2kqSZrY6H18qSZphDAlJUpEhIUkqMiQkSUWGhCSpyJCQJBUZEpKkIkNCklRkSEiSigwJSVKRISFJKjIkJElFhoQkqciQkCQVGRKSpKJaQyIilkTEQxExFhFrJmh3TkRkRIzUWZ8k6WC1hUREzAKuApYCJwLnRcSJbdodA1wM3FtXbZKk9ursSZwGjGXmjsx8GrgFWNGm3XuBDwL7aqxNktRGnSFxPPB4y3SjmveciDgFmJ+Zn55oRRGxKiI2RcSm3bt3d79SSRJQb0hEm3n53MKIFwF/BvzeZCvKzPWZOZKZI/PmzetiiZKkVnWGRAOY3zI9DOxsmT4GOAn4QkQ8CpwOjHryWpL6p86Q2AgsjohFEXEkcC4wemBhZu7NzLmZuTAzFwL3AMszc1ONNUqSWtQWEpm5H1gN3AlsB27NzAcj4rKIWF5XHZKkzh1R58YycwOwYdy8Swttz6ijJulwtHDNHW3nPzpUcyGa8WoNCUmHgXXHFubvrbcO1cJhOSRJRYaEJKnIkJAkFRkSkqQiQ0KSVGRISJKKvARW0szjZbi1sSchSSoyJCRJRYaEJKnIkJAkFRkSkqQir26SVC+vTJpR7ElIkooMCUlSkSEhSSoyJCRJRYaEJKnIkJAkFXkJrKTD2sI1d7Sd/+hQzYXMUPYkJElFhoQkqciQkCQVGRKSpCJDQpJUZEhIkopqDYmIWBIRD0XEWESsabP8XRGxLSK2RsQ/RcQP11mfJOlgtYVERMwCrgKWAicC50XEieOa3Q+MZOZrgL8DPlhXfZKkQ9XZkzgNGMvMHZn5NHALsKK1QWZ+PjO/XU3eAwzXWJ8kaZw6Q+J44PGW6UY1r2Ql8JmeViRJmlCdw3JEm3nZtmHEW4ER4I2F5auAVQALFizoVn2SpHHq7Ek0gPkt08PAzvGNIuJNwFpgeWZ+t92KMnN9Zo5k5si8efN6Uqwkqd6Q2AgsjohFEXEkcC4w2togIk4B/oJmQOyqsTZJUhu1hURm7gdWA3cC24FbM/PBiLgsIpZXzT4EHA38bURsjojRwuokSTWodajwzNwAbBg379KW12+qsx5J0sS841qSVGRISJKKDAlJUpEhIUkqMiQkSUWGhCSpyJCQJBUZEpKkIkNCklRkSEiSigwJSVKRISFJKjIkJElFhoQkqciQkCQVGRKSpCJDQpJUVOuT6SRNbuGaO9rOf3So5kIk7ElIkiZgT0KSJvFC7t3Zk5AkFRkSkqQiQ0KSVGRISJKKDAlJUpEhIUkqMiQkSUWGhCSpqNaQiIglEfFQRIxFxJo2y78vIj5ZLb83IhbWWZ8k6WC13XEdEbOAq4A3Aw1gY0SMZua2lmYrgScz80ci4lzgA8Cv1VWj1A0v5LtzdfipsydxGjCWmTsy82ngFmDFuDYrgI9Vr/8OODMiosYaJUktIjPr2VDEOcCSzHxHNf3rwOszc3VLmweqNo1q+pGqzdfGrWsVsKqa/DHgoRr+CpOZC3xt0lb9NxPqnAk1gnV2m3V212R1/nBmzptsJXUO8NeuRzA+oTppQ2auB9Z3o6huiYhNmTnS7zomMxPqnAk1gnV2m3V2V7fqrPNwUwOY3zI9DOwstYmII4Bjga/XUp0k6RB1hsRGYHFELIqII4FzgdFxbUaBC6rX5wCfy7qOh0mSDlHb4abM3B8Rq4E7gVnADZn5YERcBmzKzFHgeuDjETFGswdxbl31dcFAHf6awEyocybUCNbZbdbZXV2ps7YT15Kkmcc7riVJRYaEJKnIkJimiFgXEf8VEZurn2WFdhMORdLjGj8UEV+NiK0R8fcR8eJCu0cj4ivV32NTjfUN/DAtETE/Ij4fEdsj4sGI+J02bc6IiL0tvwuX1l1nVceEn2M0fbjan1sj4tQ+1PhjLftpc0R8MyJ+d1ybvuzPiLghInZV92sdmPfSiLg7Ih6u/nxJ4b0XVG0ejogL2rXpcZ29+65npj/T+AHWAe+epM0s4BHgBOBIYAtwYo01ngUcUb3+APCBQrtHgbk1779J9w3wW8C11etzgU/24XM+Dji1en0M8O9t6jwD+HTdtU31cwSWAZ+heT/S6cC9fa53FvDfNG/q6vv+BH4OOBV4oGXeB4E11es17b5DwEuBHdWfL6lev6TmOnv2Xbcn0VudDEXSM5l5V2burybvoXlvyqCYEcO0ZOYTmfnl6vW3gO3A8XXW0EUrgJuy6R7gxRFxXB/rORN4JDMf62MNz8nML3LofVmtv4MfA85u89ZfBO7OzK9n5pPA3cCSOuvs5XfdkHh+VlfduxsK3dDjgcdbphv07x+Yt9P8X2Q7CdwVEfdVQ57UoZN981yb6guwF5hTS3VtVIe7TgHubbP4pyJiS0R8JiJ+otbC/t9kn+Mg/T5Cs3d4c2HZIOxPgJdn5hPQ/A8D8LI2bQZtv3b1u17nsBwzTkR8FnhFm0VrgWuA99Lc6e8F/oTmh3PQKtq8t6vXHE9UY2b+Q9VmLbAf+JvCat6QmTsj4mXA3RHx1ep/K73UtWFa6hARRwOfAn43M785bvGXaR4yeao6N3U7sLjuGpn8cxyk/XkksBy4pM3iQdmfnRqk/dr177ohMYHMfFMn7SLiL4FPt1nUyVAkz8tkNVYn0X4ZODOrg5Jt1rGz+nNXRPw9zUNBvQ6JqQzT0ujnMC0RMZtmQPxNZt42fnlraGTmhoi4OiLm5riBKXutg8+x57+PU7AU+HJm/s/4BYOyPyv/ExHHZeYT1aG5XW3aNGieRzlgGPhCDbUdpFffdQ83TdO4Y7m/AjzQplknQ5H0TEQsAf4AWJ6Z3y60OSoijjnwmuYJsHZ/l26bEcO0VOdArge2Z+afFtq84sC5kog4jeb3ak99VXb8OY4Cb6uucjod2HvgUEofnEfhUNMg7M8Wrb+DFwD/0KbNncBZEfGS6rDzWdW82vT0u96rM/CH+w/wceArwFaav0jHVfN/CNjQ0m4ZzStiHqF5CKjOGsdoHivdXP1cO75GmlcXbal+Hqyzxnb7Bris+kUHGAL+tvp7/BtwQh8+55+heehga8t+XAZcBFxUtVld7bstNE8a/nQf6mz7OY6rM2g++OuR6nd3pO46qzp+gOY/+se2zOv7/qQZWk8Az9DsHaykeQ7sn4CHqz9fWrUdAa5ree/bq9/TMeA3+1Bnz77rDsshSSrycJMkqciQkCQVGRKSpCJDQpJUZEhIkooMCUlSkSEhSSoyJKQui4ifrAZ+HKrucn0wIk7qd13SdHgzndQDEfE+mneMfz/QyMzL+1ySNC2GhNQD1XhUG4F9NIeV+F6fS5KmxcNNUm+8FDia5tPshvpcizRt9iSkHoiIUZpP21tEc/DH1X0uSZoWnychdVlEvA3Yn5mfiIhZwJci4hcy83P9rk2aKnsSkqQiz0lIkooMCUlSkSEhSSoyJCRJRYaEJKnIkJAkFRkSkqSi/wMJdsR/ByjZdgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from posterior import posterior\n",
    "\n",
    "p = posterior(train_x)\n",
    "width = 0.35\n",
    "p1 = plt.bar(np.arange(data_range[0], data_range[1] + 1), p.T[:,0], width)\n",
    "p2 = plt.bar(np.arange(data_range[0], data_range[1] + 1) + width, p.T[:,1], width)\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('$P(\\omega|x)$')\n",
    "plt.legend((p1[0], p2[0]), ('$\\omega_1$', '$\\omega_2$'))\n",
    "plt.axis([data_range[0] - 1, data_range[1] + 1, 0, 1.2])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x1_test: err 15/all 100\n",
      "x2_test: err 32/all 200\n"
     ]
    }
   ],
   "source": [
    "#TODO\n",
    "#compute the number of all the misclassified x using optimal bayes decision rule\n",
    "\n",
    "# begin answer\n",
    "# compute prior\n",
    "total = np.sum(train_x)\n",
    "prior_prob = np.sum(train_x, axis=1)/total\n",
    "# computer p_x\n",
    "p_x = np.sum(train_x, axis=0)/total\n",
    "\n",
    "# compute prosterior prob of test_x1\n",
    "x1_test_idx = np.squeeze(x1_test - data_range[0])\n",
    "likelyhood_test_x1 = l[:, x1_test_idx]\n",
    "p_x_test_x1 = p_x[x1_test_idx]\n",
    "posterior_test_x1 = (likelyhood_test_x1.T* prior_prob).T / p_x_test_x1\n",
    "\n",
    "pred_test_x1 = np.argmax(posterior_test_x1, axis=0)\n",
    "err_x1 = np.sum(pred_test_x1)\n",
    "\n",
    "# compute prosterior prob of test_x2\n",
    "x2_test_idx = np.squeeze(x2_test - data_range[0])\n",
    "likelyhood_test_x2 = l[:, x2_test_idx]\n",
    "p_x_test_x2 = p_x[x2_test_idx]\n",
    "posterior_test_x2 = (likelyhood_test_x2.T* prior_prob).T / p_x_test_x2\n",
    "\n",
    "pred_test_x2 = np.argmax(posterior_test_x2, axis=0)\n",
    "err_x2 = np.sum(1-pred_test_x2)\n",
    "\n",
    "print('x1_test: err {}/all {}'.format(err_x1, x1_test.shape[1]))\n",
    "print('x2_test: err {}/all {}'.format(err_x2, x2_test.shape[1]))\n",
    "# end answer\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 3 risk\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ming total risk: 41.66259277457047\n"
     ]
    }
   ],
   "source": [
    "risk = np.array([[0, 1], [2, 0]])\n",
    "#TODO\n",
    "#get the minimal risk using optimal bayes decision rule and risk weights\n",
    "\n",
    "# begin answer\n",
    "err_x1_post = np.sum(posterior_test_x1[1, pred_test_x1==1])\n",
    "err_x2_post = np.sum(posterior_test_x2[0, (1-pred_test_x2)==1])\n",
    "cls_table = np.array([[0, err_x2_post],[err_x1_post, 0]])\n",
    "total_risk = np.multiply(risk, cls_table)\n",
    "print('ming total risk:', np.sum(total_risk))\n",
    "# end answer"
   ]
  }
 ],
 "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}