4.5 MAP categorical distribution.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As a discrete data $\\{x_i\\}_{i=1}^{I}$  where $x_i\\in\\{1,2,...6\\}$  \n",
    "$Pr(x=k|\\lambda _1.._k) = \\prod_{k=1}^{K}{\\lambda_k^{x_k}} = \\lambda_k$  \n",
    "\n",
    "Similar with ML, we choose prior **Dirichlet distribute** (conjugate function)  \n",
    "  \n",
    "  \n",
    "  \n",
    "$\\frac{\\Gamma{\\left[\\sum_{k=1}^{N}\\alpha_k\\right]}}{\\prod_{k=1}^{N}{\\Gamma[\\alpha_k]}}\\prod_{k=1}^{N}\\lambda_k^{\\alpha_k-1}$  \n",
    "$Pr(\\lambda_1..\\lambda_6|x_1..x_n) = Pr(x_1..x_n|\\theta)Pr(\\theta)$,  *it has removed evident*  \n",
    "$Pr(\\lambda_1..\\lambda_6|x_1..x_n) = \\prod_{i=1}^{N}Cat_{x_i}(\\lambda_{1..6})Dir_{\\lambda_{1..6}}(\\alpha_{1..6})$  \n",
    "$Pr(\\lambda_1..\\lambda_6|x_1..x_n) = \\frac{\\Gamma[\\sum_{k=1}^{N}\\alpha_k]}{\\prod_{k=1}^{N}{\\Gamma[\\alpha_k]}}\\prod_{k=1}^{N}\\lambda_k^{\\alpha_k+N_k-1}$  \n",
    "\n",
    "Maximize the log probability and use the Lagrange, we have:  \n",
    "$\\hat{\\lambda} = \\underset{\\lambda_1.._6}{argmax}\\left [\\log{\\frac{\\Gamma[\\sum_{k=1}^{N}\\alpha_k]}{\\prod_{k=1}^{N}{\\Gamma[\\alpha_k]}}}+\\sum_{k=1}^{N}(\\alpha_k+N_k-1)\\log{\\lambda_k} + \\nu\\left(\\sum_{k=1}^{N}\\lambda_k-1\\right)\\right]$  \n",
    "$\\hat{\\lambda} = \\frac{N_k+\\alpha_k-1}{\\sum_{k=1}^{6}(N_k+a_k-1)}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import math\n",
    "from scipy import stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[10, 100, 1000, 1000, 10, 1]\n",
      "[2565 1451  994 1000 1493 2497]\n",
      "[1 2 3 4 5 6]\n"
     ]
    }
   ],
   "source": [
    "I = 10000 # change to recognize result\n",
    "\n",
    "# alphaPrior = np.ones(6) # it a special type of ML\n",
    "alphaPrior = [10, 100, 1000, 1000, 10, 1]\n",
    "\n",
    "original_probabilities = [0.25,0.15,0.1,0.1,0.15,0.25]\n",
    "data = np.random.choice(np.arange(7)[1:], I, p=original_probabilities)\n",
    "\n",
    "# Histogram of categorical distribution\n",
    "x = np.histogram(data, bins=[1, 2, 3, 4, 5, 6, 7])\n",
    "print(alphaPrior)\n",
    "print(x[0])\n",
    "print(x[1][:-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEaCAYAAAACBmAUAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAZ1UlEQVR4nO3de7hddX3n8ffHhIuCYpToCElMUHSMSsEnQjtWvHARrYLtYy3WCzhWxlY6VmdqqXakYnWsttX2ERWqUbwgWrw0ozg0Ixa1LZhwsRYUDZFLAgoasEUUDXznj7WCm+05J2ufs8/eO4f363nOk71+6/db67tO8uxP1vrtvVaqCkmSduY+4y5AkrRrMDAkSZ0YGJKkTgwMSVInBoYkqRMDQ5LUiYGhiZOkkjxylmOvSXLkNOuenOSqqfomeV2S982u4oFr/PUk1ye5LckhQ9je3b+vJB9M8mcDjB2o/4B1/WmSj7SvV7THu2g+9qXRMDA0FO2b74/bN4XvtW9Ee4+7rl5V9eWqevQ0695SVb8DkGRl+ya8eJ5K+Qvg5Krau6oum6d9/IIkJyb5yqj216uqrmuP985x7F/DYWBomJ5TVXsDTwDWAH/S32Ee34R3JQ8HrujS0d+XJomBoaGrqq3A54HHwd2XTF6Z5NvAt9u2lyfZlGRbknVJ9uvbzLOSbE7y/SRvT3KfdtwjklyQ5Aftuo8meWDf2CcmuTLJLUk+kGTPduxTk2yZqubeyyfAl9o/b23PmJ7S1vn4nv4PSXJ7kqVTbOs+Sf4kybVJbkryoST7JNkjyW3AIuBrSa6eppbZ/L5mlOQxwHuBX2mP6dae1UuSfC7JfyS5OMkjesb95yTr2/1eleT5M+xjVZIL2+2sB/btWXePs7YkD2r/bm5o/54+09P32UkuT3Jrkn9OctAgx6r5Y2Bo6JIsB54F9F5ueS5wGLA6ydOB/w08H3gYcC1wTt9mfp3mLOUJwHHAf92x+XbsfsBjgOXAn/aNfSHwDOARwKOY4kxnJw5v/3xgexnlwra+F/X0eQHwhaq6eYrxJ7Y/TwMOAPYG3lVVd7RnYAC/VFWPmGLsDs9lsN/XjKrqG8ArgH9pj+mBPauPB94ILAE2AW8GSLIXsB44G3hI2+/dSVZPs5uzgUtoguJNwAkzlPRh4H7AY9ttv6Pd5yHAWuC/AQ8GzgDWJdljkOPVPKkqf/yZ8w9wDXAbcCvNG9q7gfu26wp4ek/f9wNv61neG/gZsLKn/zE963+P5s15qv0+F7isr45X9Cw/C7i6ff1UYEtf3yPb138KfKR9vbKtYXFP38OA64C0yxuB509T0xeA3+tZfnR7fIt7ju+RM/wuZ/P7emT7+oPAn02z3ROBr/S1fRB4X9/v65vt698CvtzX/wzg1Cm2vQLYDuzV03b2VL9TmtC7C1gyxXbeA7ypr+0q4Cnj/jfuT+H1UQ3Tc6vq/02z7vqe1/sBl+5YqKrbkvwA2J/mTby//7XtGJI8FPhr4MnA/WnOkm+ZYV93j52Lqro4ye3AU5PcCDwSWDdN9/3a/fbWsBh4KLC14y4H/X3NxXd7Xt9OE0jQzLUc1nf5ajHN2UG//YBbqupHPW3X0pwB9lsObKuq/r+3Hfs8Icnv97TtzhD+DjV3BoZGpfe2yDfQvDEAd1/6eDD3fDNdzs8nhle0YwDe0m7r8VW1LclzgXf17av3Tap37Gxq7XUWzWWp7wLnVtVPpul3j+Pj5//7/t4sa+jy+xp0m11cD1xYVUd16HsjzVzIXj2hsWKafV4PPCjJA6vq1inWvbmq3jxgrRoB5zA0Dh8DXprk4Pba9FuAi6vqmp4+f5hkSTsf8irg4237/Wkuff0wyf7AH06x/VcmWZbkQcDre8Z2dTPNJZMD+to/QjO38iLgQzOM/xjw6nYSeG+a4/t4VW0fsI7e7e3s99XF94BlSXbv2P+zwKOSvDjJbu3PE9sJ9HuoqmtpLtO9McnuSX4VeM5UG62qG2k+FPHu9u94tyQ75o3+FnhFksPS2CvJryW5/4DHqnlgYGjk2stW/wv4JM3/TB9BM6Ha6+9pJlAvBz5Hcx0fmsnZJwA/bNs/NcUuzgb+AdgMXA0M9MW0qrqdZuL3n9pP6vxy2349zaWhAr48wybW0ly2+RLwHeAnwO/P0H9n9XT5fXVxAc1Z23eTfL/Dfv8DOLrd1w00Z1Z/Dkw3Af3bNHM924BTmTlUX0wzD/NN4CbgD9p9bgReTnPWeAvNJPyJO6tVo7FjAk9SB0nWAjdU1aCfvJJ2ec5hSB0lWQn8BjDn23lIuyIvSUkdJHkT8G/A26vqO+OuRxoHL0lJkjrxDEOS1MmCncPYd999a+XKleMuQ5J2KZdccsn3q+oX7pEGCzgwVq5cycaNG8ddhiTtUpJcO906L0lJkjoxMCRJnRgYkqRODAxJUicGhiSpEwNDktTJSAMjyTHtc4E3JTllivWvSfMs5n9N8oUkvc8AuLN9zu/lSaZ7cI0kaZ6M7HsYSRYBpwNHAVuADUnWVdWVPd0uA9ZU1e1Jfhd4G81jIgF+XFUHj6peSdI9jfIM41BgU1Vtrqqf0jzE/rjeDlX1xfZZBAAXActGWJ8kaQaj/Kb3/tzzOcVbaB62Mp2X0TyVa4c9k2ykedTlW6vqM/0DkpwEnASwYsWKORX7jvXfmtP4UXn1UY8adwnSvca9/X1hIm8NkuRFwBrgKT3ND6+qrUkOAC5I8vWqurp3XFWdCZwJsGbNGm/DK0lDNMpLUluB5T3Ly5jiIfZJjqR5DvOxVXXHjvaq2tr+uRn4R3yIjSSN1CgDYwNwYJJV7UPojwfu8WmnJIcAZ9CExU097UuS7NG+3hd4EtA7WS5JmmcjuyRVVduTnAycDywC1lbVFUlOAzZW1Trg7cDewN8lAbiuqo4FHgOckeQumpB7a9+nqyRJ82ykcxhVdR5wXl/bG3peHznNuH8GHj+/1UmSZuI3vSVJnRgYkqRODAxJUicGhiSpEwNDktSJgSFJ6sTAkCR1YmBIkjoxMCRJnRgYkqRODAxJUicGhiSpEwNDktSJgSFJ6sTAkCR1YmBIkjoxMCRJnRgYkqRODAxJUicGhiSpEwNDktSJgSFJ6sTAkCR1YmBIkjoxMCRJnRgYkqRODAxJUicGhiSpEwNDktSJgSFJ6sTAkCR1YmBIkjoZaWAkOSbJVUk2JTllivWvSXJlkn9N8oUkD+9Zd0KSb7c/J4yybknSCAMjySLgdOCZwGrgBUlW93W7DFhTVQcB5wJva8c+CDgVOAw4FDg1yZJR1S5JGu0ZxqHApqraXFU/Bc4BjuvtUFVfrKrb28WLgGXt62cA66tqW1XdAqwHjhlR3ZIkRhsY+wPX9yxvadum8zLg84OMTXJSko1JNt58881zLFeS1GsiJ72TvAhYA7x9kHFVdWZVramqNUuXLp2f4iTpXmqUgbEVWN6zvKxtu4ckRwKvB46tqjsGGStJmj+jDIwNwIFJViXZHTgeWNfbIckhwBk0YXFTz6rzgaOTLGknu49u2yRJI7J4VDuqqu1JTqZ5o18ErK2qK5KcBmysqnU0l6D2Bv4uCcB1VXVsVW1L8iaa0AE4raq2jap2SdIIAwOgqs4Dzutre0PP6yNnGLsWWDt/1UmSZjKRk96SpMljYEiSOjEwJEmdGBiSpE4MDElSJwaGJKkTA0OS1ImBIUnqxMCQJHViYEiSOjEwJEmdGBiSpE4MDElSJwaGJKkTA0OS1ImBIUnqxMCQJHViYEiSOjEwJEmdDBQYSS5PcnKSJfNVkCRpMg16hvE54LXADUk+luSIeahJkjSBBgqMqno98HDgN4BFwOeSfCfJG5KsmI8CJUmTYeA5jGp8vqqeD+wHnAm8Dtic5Pwkxwy7SEnS+M160jvJLwNvBU4BbgDeCFwNnJvknUOpTpI0MRYP0jnJQ4CXAC8FHgGsA55XVet7+nwYWA/8wfDKlCSN20CBAWwBNgHvB86qqu9P0ecKYMNcC5MkTZZBA+OIqvryTB2q6t+Bp82+JEnSJBp0DmNbkoP6G5MclGT1kGqSJE2gQQPjTOBxU7SvbtdJkhaoQQPjIOCrU7RvAB4/93IkSZNq0MC4E9hnivYlQOZejiRpUg0aGBcCr0+yaEdDksXA64EvDbMwSdJkGTQwXgs8GdiU5MPtdy6+Dfwq8Ic7G5zkmCRXJdmU5JQp1h+e5NIk25M8r2/dne3NDy9Psm7AuiVJczTQx2qr6qr2U1InAwe3zR8F3l1VN8w0tj0rOR04iub7HBuSrKuqK3u6XQecCPzPKTbx46o6eIp2SdIIDPo9DKrqRppLUIM6FNhUVZsBkpwDHAfcHRhVdU277q5ZbF+SNI8GDowk96M5u3gIfZe0qupTMwzdH7i+Z3kLcNgAu94zyUZgO/DWqvrMAGMlSXM06L2kjgQ+Bjx4itVFc8vz+fLwqtqa5ADggiRfr6qr++o7CTgJYMUK77be6x3rvzXuEjp59VGPGncJGiL/3S0sg056/zXNQ5SWVdV9+n52FhZbgeU9y8vatk6qamv752bgH4FDpuhzZlWtqao1S5cu7bppSVIHgwbGSuBNO5vgnsYG4MAkq5LsDhxPc7fbnUqyJMke7et9gSfRM/chSZp/gwbGPwGPns2Oqmo7zaerzge+AXyiqq5IclqSYwGSPDHJFuA3gTOSXNEOfwywMcnXgC/SzGEYGJI0QoNOer8X+Isk+wFfB37Wu7KqLp1pcFWdB5zX1/aGntcbaC5V9Y/7Z7z1iCSN1aCBcW7751Q3GpzvSW9J0hgNGhir5qUKSdLEG/Sb3tfOVyGSpMk26KQ3SZ6Z5LNJrkyyvG37nSRHDL88SdKkGCgwkrwQ+ATNDQdXAbu1qxbR3JhQkrRAzeZutS+vqlfT3KJjh4v4+c0IJUkL0KCBcSDwL1O03wY8YO7lSJIm1aCBcQMw1U1XDgeunqJdkrRADBoYZwJ/k+RJ7fLyJCcAbwPeM9TKJEkTZdCP1b4tyT7AemBPmtt03AH8RVWdPg/1SZImxGweoPT6JG8GVtOcoVxZVbcNvTJJ0kQZODAAqup2YOOQa5EkTbBBH6A04+3Iq+rYuZUjSZpUg55h/KBveTfgl2gejDTT41klSbu4QSe9XzpVe5K/BP59KBVJkibSwPeSmsYZwCuHtC1J0gQaVmDM6il8kqRdx6CT3n/T3wQ8DHgmsHZYRUmSJs+gk979j0m9C7gZeDUGhiQtaINOej9tvgqRJE22Yc1hSJIWuEHnML4IVJe+VfX0WVUkSZpIg85hfAN4IfBd4OK27VDgPwFnA3cOrzRJ0iQZNDDuAM4CXlVVd59pJHknkKp61RBrkyRNkEHnMF4CvKs3LFrvBl48nJIkSZNo0MAIv/jRWqZpkyQtIINekloLvC/JgcBFbdsvA68FPjDMwiRJk2XQwHgtcBPwKuAtbduNwFuBvxxiXZKkCTPoF/fuonl+99uSPKBt8y61knQvMKsv7iVZQ3P/qDvb5b2SzOrpfZKkXcOgX9x7KPD3NN+9KOBAYDPwV8BPaC5VSZIWoEHPMN4BfA94MHB7T/vfAUcPqyhJ0uQZ9DLSEcARVXVLkt72q4EVQ6tKkjRxBj3DuC/w0ynal9JckppRkmOSXJVkU5JTplh/eJJLk2xP8ry+dSck+Xb7c8KAdUuS5mjQwPgScGLPciVZBPwR8IWZBrb9TqeZLF8NvCDJ6r5u17XbP7tv7IOAU4HDaOZPTk2yZMDaJUlzMJvvYVyY5InAHjTfvXgssA/wpJ2MPRTYVFWbAZKcAxwHXLmjQ1Vd0667q2/sM4D1VbWtXb8eOAb42ID1S5JmadDvYVyZ5PHA79LciHBPmgnv06vqxp0M3x+4vmd5C80ZQxdTjd2/v1OSk4CTAFascEploXvH+m+Nu4ROXn3Uozr1W2jHo4Wnc2Ak2Q34CvCSqjp1/kqavao6EzgTYM2aNZ2e2yFJ6qbzHEZV/QxYRccHKE1hK7C8Z3lZ2zbfYyVJQzDopPdZwMtnua8NwIFJViXZHTgeWNdx7PnA0UmWtJPdR7dtkqQRGXTSey/ghUmOAi4BftS7sqr++3QDq2p7kpNp3ugXAWur6ookpwEbq2pdO5n+aWAJ8Jwkb6yqx1bVtiRvogkdgNN2TIBLkkajU2AkOQi4AngMcGnbfEBft51eqqqq84Dz+tre0PN6A83lpqnGrqW5vbokaQy6nmFcBjysqp4GkORzwO90+GSUJGmB6DqHkb7lJ9N861uSdC8xq9ub84sBIkla4LoGRvGLcxR+z0GS7kW6zmEE+EiSO9rlPYG/TdJ7i3Oq6thhFidJmhxdA+OsvuWPDLsQSdJk6xQYVfXS+S5EkjTZZjvpLUm6lzEwJEmdGBiSpE4MDElSJwaGJKkTA0OS1ImBIUnqxMCQJHViYEiSOjEwJEmdGBiSpE4MDElSJwaGJKkTA0OS1ImBIUnqxMCQJHViYEiSOjEwJEmdGBiSpE4MDElSJwaGJKkTA0OS1ImBIUnqxMCQJHViYEiSOhlpYCQ5JslVSTYlOWWK9Xsk+Xi7/uIkK9v2lUl+nOTy9ue9o6xbkgSLR7WjJIuA04GjgC3AhiTrqurKnm4vA26pqkcmOR74c+C32nVXV9XBo6pXknRPozzDOBTYVFWbq+qnwDnAcX19jgPOal+fCxyRJCOsUZI0jVEGxv7A9T3LW9q2KftU1Xbgh8CD23WrklyW5MIkT55qB0lOSrIxycabb755uNVL0r3crjLpfSOwoqoOAV4DnJ3kAf2dqurMqlpTVWuWLl068iIlaSEbZWBsBZb3LC9r26bsk2QxsA/wg6q6o6p+AFBVlwBXA4+a94olSXcbZWBsAA5MsirJ7sDxwLq+PuuAE9rXzwMuqKpKsrSdNCfJAcCBwOYR1S1JYoSfkqqq7UlOBs4HFgFrq+qKJKcBG6tqHfB+4MNJNgHbaEIF4HDgtCQ/A+4CXlFV20ZVuyRphIEBUFXnAef1tb2h5/VPgN+cYtwngU/Oe4GSpGntKpPekqQxMzAkSZ0YGJKkTgwMSVInBoYkqRMDQ5LUiYEhSerEwJAkdWJgSJI6MTAkSZ0YGJKkTgwMSVInBoYkqRMDQ5LUiYEhSerEwJAkdWJgSJI6MTAkSZ0YGJKkTgwMSVInBoYkqRMDQ5LUiYEhSerEwJAkdWJgSJI6MTAkSZ0YGJKkTgwMSVInBoYkqRMDQ5LUiYEhSerEwJAkdWJgSJI6GWlgJDkmyVVJNiU5ZYr1eyT5eLv+4iQre9b9cdt+VZJnjLJuSdIIAyPJIuB04JnAauAFSVb3dXsZcEtVPRJ4B/Dn7djVwPHAY4FjgHe325MkjcgozzAOBTZV1eaq+ilwDnBcX5/jgLPa1+cCRyRJ235OVd1RVd8BNrXbkySNyOIR7mt/4Pqe5S3AYdP1qartSX4IPLhtv6hv7P79O0hyEnBSu3hbkquGU/rQ7At8f5gbfM0wNza4hXY8sPCOaaEdDyy8Y5q043n4dCtGGRjzrqrOBM4cdx3TSbKxqtaMu45hWWjHAwvvmBba8cDCO6Zd6XhGeUlqK7C8Z3lZ2zZlnySLgX2AH3QcK0maR6MMjA3AgUlWJdmdZhJ7XV+fdcAJ7evnARdUVbXtx7efoloFHAh8dUR1S5IY4SWpdk7iZOB8YBGwtqquSHIasLGq1gHvBz6cZBOwjSZUaPt9ArgS2A68sqruHFXtQzSxl8tmaaEdDyy8Y1poxwML75h2meNJ8x94SZJm5je9JUmdGBiSpE4MjBFIsjbJTUn+bdy1DEOS5Um+mOTKJFckedW4a5qLJHsm+WqSr7XH88Zx1zQMSRYluSzJZ8ddyzAkuSbJ15NcnmTjuOsZhiQPTHJukm8m+UaSXxl3TTNxDmMEkhwO3AZ8qKoeN+565irJw4CHVdWlSe4PXAI8t6quHHNps9LeTWCvqrotyW7AV4BXVdVFOxk60ZK8BlgDPKCqnj3ueuYqyTXAmqoa6pfcxinJWcCXq+p97adH71dVt465rGl5hjECVfUlmk99LQhVdWNVXdq+/g/gG0zxzftdRTVuaxd3a3926f9JJVkG/BrwvnHXoqkl2Qc4nObToVTVTyc5LMDA0By1dxQ+BLh4zKXMSXv55nLgJmB9Ve3SxwO8E3gtcNeY6ximAv4hySXtbYB2dauAm4EPtJcO35dkr3EXNRMDQ7OWZG/gk8AfVNW/j7ueuaiqO6vqYJq7CByaZJe9dJjk2cBNVXXJuGsZsl+tqifQ3PH6le2l3l3ZYuAJwHuq6hDgR8AvPPZhkhgYmpX2Wv8ngY9W1afGXc+wtJcEvkhzG/1d1ZOAY9tr/ucAT0/ykfGWNHdVtbX98ybg0+z6d6zeAmzpOZs9lyZAJpaBoYG1k8TvB75RVX817nrmKsnSJA9sX98XOAr45liLmoOq+uOqWlZVK2nulnBBVb1ozGXNSZK92g9Y0F62ORrYpT91WFXfBa5P8ui26Qiau1lMrAV1t9pJleRjwFOBfZNsAU6tqvePt6o5eRLwYuDr7XV/gNdV1XnjK2lOHgac1T6U6z7AJ6pqQXwUdQF5KPDp5v8qLAbOrqr/O96ShuL3gY+2n5DaDLx0zPXMyI/VSpI68ZKUJKkTA0OS1ImBIUnqxMCQJHViYEiSOjEwJEmdGBiSpE4MDElSJwaGNALt7UduTHJqT9tBSX6S5DfHWZvUld/0lkYkyTOA/wM8Bbgc2Ah8taom+nYQ0g4GhjRCSd4JHAtcCDwZOLjn4U3SRDMwpBFKsgfwNeBA4L8sgAc16V7EOQxptFYCy2meHnfAeEuRBuMZhjQi7UOnLgK+RfNI21OBX6qq68ZamNSRgSGNSJK3Ar8NHAT8EPg8sCfw9KpaSM/e1gLlJSlpBJI8BfgfwEuq6tZq/qd2IrAa+KNx1iZ15RmGJKkTzzAkSZ0YGJKkTgwMSVInBoYkqRMDQ5LUiYEhSerEwJAkdWJgSJI6+f/5jsthOtUCHgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plt.style.use('ggplot')\n",
    "# plt.hist(data, bins=6, color = 'b')\n",
    "plt.bar(np.array(x[1][:-1]), np.array(original_probabilities), align='center', alpha=0.5)\n",
    "plt.xlabel('x', size = 14)\n",
    "plt.ylabel('Frequency', size = 14)\n",
    "plt.title('Probability of roll the dice')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def MapCategoricalNorm(N, k, alpha):\n",
    "    sumLow = 0\n",
    "    for idx in range(k):\n",
    "        sumLow += N[idx] + alpha[idx] - 1\n",
    "    lambda_ = (N + alpha - 1)/sumLow\n",
    "    return lambda_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Maximum likelihood: [0.21246389 0.12794057 0.16450681 0.16500206 0.12397854 0.20610813]\n",
      "Original: [0.25, 0.15, 0.1, 0.1, 0.15, 0.25]\n"
     ]
    }
   ],
   "source": [
    "lambdaList = MapCategoricalNorm(x[0], 6, alphaPrior)\n",
    "\n",
    "print('MAP:', lambdaList)\n",
    "print('Original:', original_probabilities)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEJCAYAAACaFuz/AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAUuUlEQVR4nO3df7BfdX3n8eeLUOgK6qJkf5QAiTVU41ahe4HusLpdBQn9Ae0sKjhsseJk7MpOF+p06I+FNt3OaNsVu1NsSSUttTpIce2kbiyTBewfWm1ukEUTBEKkkGjHSGjVomDkvX98T5gvlxvuubn3e+/9fu7zMfOdfM/nfM7J+8K5r3zu5/s556aqkCS166jFLkCSNFoGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wz6JSzJw0nePctjKsnF81zHryf54nyeU9LCMehHLMlJSTYl2ZvkqST7kvxRklU9Dj8T+MAs/8p/Dfzl7CuVRiPJn3QDkJum2ffebt8nprT/yyTfSfJIkufkVDcIqu71RJIvJtkwyq9jnBn0I5RkDTAJ/BvgcuDlwGXAq4DtSVYf5rhjAKpqf1U9MZu/s6r+vqqenEvd0gg8Crw5yXGHGpIcDfws8Mg0/S9nMGD5DnD+Yc65kcHA5tXAXwA3JnnLPNbcDIN+tG4AngbOrao7quqRqroLOLdrvwEgyaeS/EGS302yH/h01/6sqZskpyX5626kc3+SH0/yrSRvG+rzzNRNktXd9n9Ksq0b+exKct5Q/xVJbkry5STfTvJgkl+abhQlzcG9wIPAm4fafoJBkH9qmv5vB/4U+BBwxWHO+c1uYLO7qn6tO/9Pz1fBLfGbeUSSvARYD9wwdVTebX8AuCDJCV3zZUCA1zIY5Uw931HAx4GDwI8CbwOuA47tUc5vAf8LeA2wHbglyfHdvqOAfQy+AV8J/CrwK8DP9fxSpb5uYhDgh7wd+GPgWc9hSfJa4KXAXwF/BvxkkpU9zv8d4Pvmp9S2GPSjs5ZBcN93mP27uv1ru+0vV9UvVtWXqmq6Y84Dfgj42aq6p6r+BrgKOLpHLddX1V9W1YMMQvwlwOkAVfXdqrq2qrZX1cNVdSvwh8Cl/b5MqbePABNJ1ib5VwwGQn8yTb8rgI921+aXgb9lMJUzrSRHdz/V/jBwx7xX3YA+IaGFsWOG/a8AvlJV+4batjOYAprJvUPvv9L9+S8ONSR5J/AO4FTgnzEYFf1dj/NKvVXV40k+zmAk/w/Ap6rqkSTP9EnyIuBNwOuHDv0QcDXwu1NO+VtJfp3BT7VPAb8D3Diq+seZQT86uxn8SLqOwZTLVOu6/bu77X8aYS3fPfSmqqr7xjoKoPvw6v3Au4HPAN8A3gX8zAjr0fK1GbgZ+BZw7TT73wq8APj08D8AwIok51TVp4fa3sdgOugJ4Kvlo3gPy6mbEamqx4Dbgf+S5AXD+7rtdwGfrKoDPU/5JeAHkvzAUNsEc/9/+O+Bz1XV71fV3VW1G/jBOZ5TOpw7GIy+T2SwUmaqK4DfZzC1OPz6Pzz3Q9nHug9iv2LIPz+DfrSuZPBT0/9N8vokJyf5MWAbg/n5K2dxrm3A/cDNSV6T5EcZjGgOMuXDrFl6APiRJBd0c6f/HfgPczifdFhdIL8aWDN1GXCSVzMYvPxRVX1x+MVg+ubNSV648FWPP4N+hKrqIQYX7k4GF+oeBh9I3Qec2X3Q1PdcTzOYTjmWwYdTNzNYTVMMVhscqRuBW7u6tgOrgf85h/NJz6uqvllV35hm1zuAB6vq3mn2fYJBXrlI4AjEn3jGV5LXAPcAE1U104e5kpYpg36MJPkZBh/aPshg5P0+BlNAZzhHKelwnLoZLy9k8EHVLuDDDKaAzjfkj0yS9d0dxruTXDPN/qu7O4nvTXJHklOH9n0vyT3da8vCVi7NjiN6LUtJVjD4IPo8YC+DzycurapdQ33+I4MVSU8k+Xngx6rqLd2+b1XV8dOcWlpyHNFruToL2F1Ve6rqKeAW4KLhDlV119DjKz4L9HniqLTkLLkbpk488cRavXr1Ypehhu3YsePrwM8zeKLiIXuBs5/nsCuATw5tf3+SSQbLW99TVX8x3UHdo3M3ABx33HH/9hWveMUcKpcOb8eOHV+vqmmfCbTkgn716tVMTk4udhlqWJJZPd4hyWUMlskO319walXtS/Iy4M4kX+iW0z5LVW0CNgFMTEyU17ZG5fmua6dutFztA04e2l7VtT1LknMZPNHzwuEbfA49c6iq9jB4zO4ZoyxWmguDXsvVdmBtkjXdL3q5BHjW6pkkZzC4oezCqvraUPsJSY7t3p8InMNgJZS0JC25qRtpIVTVwSRXMnge0Qpgc1XtTLIRmKyqLQyehng88OfdA7YeqaoLGTy3/8YkTzMYLL1neLWOtNQY9Fq2qmorsHVK27VD7889zHGfYfDsc2ksOHUjSY0z6CWpcb2C3lvFJWl8zThH390qfgNDt4on2TLlw6fPM3iC4qFbxX8beEu379tVdfr8li1J6qvPiN5bxSVpjPVZdXMSI75VfPg28VNOOeWwJ75+2wM9yu3nqvNOm7dzaX74/1cajXldXnmkt4pPvU18PmuSpOWuz9SNt4pL0hjrE/TeKi5JY2zGqRtvFZek8dZrjt5bxSVpfHlnrCQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcfP6y8ElHZnrtz0wb+e66rzT5u1caoMjeklqnEEvSY0z6CWpcQa9JDXOoJekxrnqRpJGaCmsqHJEL0mNM+i1rCVZn+T+JLuTXDPN/quT7Epyb5I7kpw6tO/yJA92r8sXtnKpP4Ney1aSFcANwAXAOuDSJOumdPs8MFFVrwZuA367O/YlwHXA2cBZwHVJTlio2qXZMOi1nJ0F7K6qPVX1FHALcNFwh6q6q6qe6DY/C6zq3p8PbKuqA1X1OLANWL9AdUuzYtBrOTsJeHRoe2/XdjhXAJ+czbFJNiSZTDK5f//+OZYrHRmDXuohyWXABPA7szmuqjZV1URVTaxcuXI0xUkzMOi1nO0DTh7aXtW1PUuSc4FfBS6sqidnc6y0FBj0Ws62A2uTrElyDHAJsGW4Q5IzgBsZhPzXhnbdDrwxyQndh7Bv7NqkJccbprRsVdXBJFcyCOgVwOaq2plkIzBZVVsYTNUcD/x5EoBHqurCqjqQ5DcZ/GMBsLGqDizClyHNqFfQJ1kP/B6Db4YPVtV7puy/GngHcBDYD7y9qv6u23c58Gtd1/9RVTfPU+3SnFXVVmDrlLZrh96f+zzHbgY2j646aX7MOHXjWmNJGm995uhdayxJY6xP0LvWWJLG2LyuunGtsSQtPX2C3rXGkjTG+gS9a40laYzNuLzStcaSNN56raN3rbEkjS8fgSBJjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJalyvXzwiSYvp+m0PzNu5rjrvtHk717gw6CXNmUG8tDl1I0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQa9lK8n6JPcn2Z3kmmn2vy7J3UkOJrl4yr7vJbmne21ZuKql2fPOWC1LSVYANwDnAXuB7Um2VNWuoW6PAG8D3j3NKb5dVaePuk5pPhj0Wq7OAnZX1R6AJLcAFwHPBH1VPdzte3oxCpTmi1M3Wq5OAh4d2t7btfX1/Ukmk3w2yU8frlOSDV2/yf379x9hqdLcGPTSkTm1qiaAtwLvT/KD03Wqqk1VNVFVEytXrlzYCqWOQa/lah9w8tD2qq6tl6ra1/25B/gUcMZ8FifNp15B7+oENWg7sDbJmiTHAJcAva7PJCckObZ7fyJwDkNz+9JSM+OHsa5OUIuq6mCSK4HbgRXA5qramWQjMFlVW5KcCXwcOAH4qSS/UVWvAl4J3Nh9SHsU8J4p3w/SktJn1Y2rE9SkqtoKbJ3Sdu3Q++0MpnSmHvcZ4IdHXqA0T/pM3SzI6gRJ0mgsxDr6U6tqX5KXAXcm+UJVPTTcIckGYAPAKaecsgAlTc9fh/b8/O8jjac+I/qRr05wCZokjU6foHd1giSNsRmDvqoOAodWJ9wH3HpodUKSCwGSnJlkL/AmBqsRdnaHvxKYTPL/gLtwdYIkLbhec/SuTpCk8eWdsZLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TG+cvBpWXAB9Itb47oJalxBr0kNc6pG0nL2nKY1nJEL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6LWsJVmf5P4ku5NcM83+1yW5O8nBJBdP2Xd5kge71+ULV7U0Owa9lq0kK4AbgAuAdcClSdZN6fYI8DbgI1OOfQlwHXA2cBZwXZITRl2zdCQMei1nZwG7q2pPVT0F3AJcNNyhqh6uqnuBp6ccez6wraoOVNXjwDZg/UIULc2WQa/l7CTg0aHtvV3bvB2bZEOSySST+/fvP+JCpbnoFfTOY0pHpqo2VdVEVU2sXLlyscvRMjVj0DuPqYbtA04e2l7VtY36WGlB9RnRO4+pVm0H1iZZk+QY4BJgS89jbwfemOSEbvDyxq5NWnKO7tFnurnIs3uev/c8JrAB4JRTTul56vFz/bYH5u1cV5132oKfvzVVdTDJlQwCegWwuap2JtkITFbVliRnAh8HTgB+KslvVNWrqupAkt9k8I8FwMaqOrAoX4g0gz5BP3JVtQnYBDAxMVGLXI6WkaraCmyd0nbt0PvtDKZlpjt2M7B5pAVK86DP1I3zmJI0xvoEvfOYkjTGZgz6qjoIHJrHvA+49dA8ZpILAZKcmWQv8CbgxiQ7u2MPAIfmMbfjPKYkLbhec/TOY0rS+PLOWElqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g17LWpL1Se5PsjvJNdPsPzbJR7v9n0uyumtfneTbSe7pXn+44MVLPR292AVIiyXJCuAG4DxgL7A9yZaq2jXU7Qrg8ap6eZJLgPcCb+n2PVRVpy9kzdKR6DWid9SjRp0F7K6qPVX1FHALcNGUPhcBN3fvbwPekCQLWKM0ZzMG/dCo5wJgHXBpknVTuj0z6gGuZzDqOeShqjq9e71znuqW5sNJwKND23u7tmn7VNVB4B+Bl3b71iT5fJK/TvLa6f6CJBuSTCaZ3L9///xWL/XUZ0TvqEd6rq8Cp1TVGcDVwEeSvGhqp6raVFUTVTWxcuXKBS9Sgn5B76hHrdoHnDy0vaprm7ZPkqOBFwOPVdWTVfUYQFXtAB4CTht5xdIRGPWqG0c9Wsq2A2uTrElyDHAJsGVKny3A5d37i4E7q6qSrOymNUnyMmAtsGeB6pZmpc+qm9mMevZOGfUU8CQMRj1JDo16JudauDRXVXUwyZXA7cAKYHNV7UyyEZisqi3ATcCHkuwGDjD4xwDgdcDGJN8FngbeWVUHFv6rkGbWJ+ifGfUwCPRLgLdO6XNo1PM3TBn1AAeq6nuOerQUVdVWYOuUtmuH3n8HeNM0x30M+NjIC5TmwYxB76hHksZbrxumHPVI0vjyEQiS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktS4XkGfZH2S+5PsTnLNNPuPTfLRbv/nkqwe2vfLXfv9Sc6fx9qlOfPa1nIwY9AnWQHcAFwArAMuTbJuSrcrgMer6uXA9cB7u2PXAZcArwLWAx/ozictOq9tLRd9RvRnAburak9VPQXcAlw0pc9FwM3d+9uANyRJ135LVT1ZVV8Gdnfnk5YCr20tC0f36HMS8OjQ9l7g7MP1qaqDSf4ReGnX/tkpx5409S9IsgHY0G1+K8n9vaqfg6vndviJwNdHeP4ZzeH8M9Y+x/P3ssj1n9r96bX9bF4bc7TIuXLq4Xb0CfqRq6pNwKbFrqOvJJNVNbHYdRyJca4dxq/+cbq2x+2/7VTWf3h9pm72AScPba/q2qbtk+Ro4MXAYz2PlRaL17aWhT5Bvx1Ym2RNkmMYfAC1ZUqfLcDl3fuLgTurqrr2S7qVC2uAtcDfzk/p0px5bWtZmHHqppuXvBK4HVgBbK6qnUk2ApNVtQW4CfhQkt3AAQbfMHT9bgV2AQeBd1XV90b0tSyksfhR/DDGuXaYx/q9tp/Da2Nxjaz+DAYnkqRWeWesJDXOoJekxhn0PSU5OcldSXYl2ZnkFxa7piORZEWSzyf5xGLXMltJ/nmS25J8Kcl9Sf7dYtfUghauba/r57ck1tGPiYPAL1bV3UleCOxIsq2qdi12YbP0C8B9wIsWu5Aj8HvAX1XVxd0qmRcsdkGNaOHa9rp+Ho7oe6qqr1bV3d37bzK4qJ5zJ+RSlmQV8BPABxe7ltlK8mLgdQxWwVBVT1XVPyxqUY0Y92vb63pmBv0R6J5geAbwuUUuZbbeD/wS8PQi13Ek1gD7gT/ufkT/YJLjFruo1ozptf1+vK6fl0E/S0mOBz4G/Leq+sZi19NXkp8EvlZVOxa7liN0NPAjwB9U1RnAPwHPeaywjtw4Xtte1/0Y9LOQ5PsYfCN8uKr+92LXM0vnABcmeZjBUxpfn+TPFrekWdkL7K2qQyPN2xh8g2gejPG17XXdg0HfU/do2puA+6rqfYtdz2xV1S9X1aqqWs3g7s47q+qyRS6rt6r6e+DRJD/UNb2BwV2pmqNxvra9rvtx1U1/5wD/GfhCknu6tl+pqq2LV9Ky81+BD3crE/YAP7fI9bTCa3txjfy69hEIktQ4p24kqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWrc/wdkYIeuHB9YkgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axs = plt.subplots(1, 2)\n",
    "axs[0].bar(np.array(x[1][:-1]), np.array(original_probabilities), align='center', alpha=0.5)\n",
    "axs[0].set_title('Original', size = 14)\n",
    "rang = np.arange(0, 0.3, 0.05)\n",
    "axs[0].set_yticks(rang)\n",
    "\n",
    "axs[1].bar(np.array(x[1][:-1]), np.array(lambdaList), align='center', alpha=0.5)\n",
    "axs[1].set_title('MAP', size = 14)\n",
    "axs[1].set_yticks(rang)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "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.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}