2018-05-25-quick-example-of-dirichlet-distribution.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Samples from Dirichlet distribution\n",
    "\n",
    "The Dirichlet distribution is a distribution over distributions!\n",
    "In Bayesian methods, it is used as a prior for categorical and multinomial distributions. The Dirichlet distribution appears in natural language processing in [Latent Dirichlet allocation](https://en.wikipedia.org/wiki/Latent_Dirichlet_allocation) and [Bayesian HMMs](http://www.aclweb.org/anthology/P07-1094).\n",
    "\n",
    "In this quick post, I'll sample from `pymc3`'s `Dirichlet` distribution using different values of concentration parameters and plot what the corresponding distributions. \n",
    "\n",
    "For more detailed information:\n",
    " - [Wikipedia article](https://en.wikipedia.org/wiki/Dirichlet_distribution).\n",
    " - [This document](http://mayagupta.org/publications/FrigyikKapilaGuptaIntroToDirichlet.pdf) seems to give a good background on the mathematics."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "import pymc3 as pm\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# helper functions you can skip over :D\n",
    "def hide_ticks(plot):\n",
    "    plot.axes.get_xaxis().set_visible(False)\n",
    "    plot.axes.get_yaxis().set_visible(False)\n",
    "\n",
    "SAVE = True\n",
    "def maybe_save_plot(filename):\n",
    "    if SAVE:\n",
    "        plt.tight_layout()\n",
    "        plt.savefig('images/' + filename, bbox_inches=\"tight\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Eh, a story might help explain. I'm randomly choosing one of three types of cookies, and I might want to weight some cookies higher than others. I don't know how exactly they should be weighted, but I have an idea of how I want them to be weighted. For example, it might be one of:\n",
    "\n",
    " - One type of cookie (I don't know which) is more frequent than the others. \n",
    " - Each type of cookie is equally likely.\n",
    " - One particular cookie type is weighted higher than the others.\n",
    " - I have no idea. \n",
    "\n",
    "By adjusting the concentration parameters `a` in the `pymc3.Dirichlet` distribution, I can add this information to my prior.\n",
    "\n",
    "To try this out, I'll sample a few examples from four Dirichlet distributions with different concentration parameters. \n",
    "\n",
    "Heads up, in Bayesian statistics, the prior involves the distribution that these samples come from, not the samples individually. If you saw my post about [Bayesian linear regression](2018-01-03-bayesian-linreg.ipynb), each of these samples is like a single line from the samples from the prior. The real prior distribution includes all of them!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Only 20 samples in chain.\n",
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (2 chains in 2 jobs)\n",
      "NUTS: [dir_skewed_stickbreaking__, dir_ten_stickbreaking__, dir_one_stickbreaking__, dir_tenth_stickbreaking__]\n",
      "100%|██████████| 520/520 [00:09<00:00, 53.02it/s]\n",
      "There were 3 divergences after tuning. Increase `target_accept` or reparameterize.\n",
      "The acceptance probability does not match the target. It is 0.939145954445564, but should be close to 0.99. Try to increase the number of tuning steps.\n",
      "There were 5 divergences after tuning. Increase `target_accept` or reparameterize.\n"
     ]
    }
   ],
   "source": [
    "# how many categories to sample from\n",
    "values = 3\n",
    "\n",
    "with pm.Model() as model:\n",
    "    # borrowing these numbers from http://mayagupta.org/publications/FrigyikKapilaGuptaIntroToDirichlet.pdf\n",
    "    pm.Dirichlet('dir_tenth', a=0.1 * np.ones(values))\n",
    "    pm.Dirichlet('dir_one', a=1 * np.ones(values))\n",
    "    pm.Dirichlet('dir_ten', a=10 * np.ones(values))\n",
    "    pm.Dirichlet('dir_skewed', a=np.array([2, 5, 15]))\n",
    "    \n",
    "    # just sample a small number.\n",
    "    trace = pm.sample(20, nuts_kwargs=dict(target_accept=.99))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For each of the four distribution, I'll grab ten samples. These samples are another distribution! I asked for three numbers, so the samples are three numbers that sum to one.\n",
    "I'll plot these ten sampled distributions as a stacked bar plot. They should all sum to one. They will also reflect the concentration parameter values:\n",
    "\n",
    " - `dir_tenth`: One type of cookie (I don't know which) is more frequent than the others. \n",
    " - `dir_one`: I have no idea. The *probabilities* are chosen from a uniform distribution.\n",
    " - `dir_ten`: Each type of cookie is equally likely.\n",
    " - `dir_skewed`: One particular cookie type is weighted higher than the others.\n",
    "\n",
    "These can be used as the values of a discrete probability distribution, for example as the probabilities in a multinomial distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABHgAAAEYCAYAAAAnPkG+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAEPtJREFUeJzt3X2MZXddx/HPF7ftClQQQbSV7a4kYkmWyB9IjUAbNRKEgkFAXMWFP1SMDzEBYX34Yw00FOIfaoRIiELDoyhigCWmiGHlSZQoOvGhBnRgAfEJ2iAgCP35xz2lw6a7zOzM3XO/c1+vpNmZnXvP+d12zrdn33Pu2RpjBAAAAIC+7jb3AgAAAADYHYEHAAAAoDmBBwAAAKA5gQcAAACgOYEHAAAAoDmBBwAAAKA5gQcAAACgOYFnTVXVK6rq+VX1yKq6Ze713JWqenpVvWvudQDL1WEeAfuP2QPs1jLnSFWdrKpX7eU2d+OO1zr3Ojg/gWfNjTHeOcZ40IU8t6o2q+r79mIdVXW4qkZVHdiL7QH97GYeAVyoVTkXAvpyDsOqEHg4J7EFWBXmETAHswfYLXOEi0ngWRNV9dCq+uuq+nRV/X6Sg9PvX1dVH93yuM2qem5V/V2Sz5xrIFXVK5McSvLmqvqfqnrO9PvXVNV7qurWqvrbqrpuy3PeUVXPq6p3T+u4uaruO335z6dfb522911bnvfrVfWpqvrXqnrMHv5rAWaw1/NoeuzV04y5tar+vqoev+Vrr6iqF1fVqWmf76uqB275+rdX1duq6pNVdUtVPWU5rxyYU4NzIWDFLekc5rlV9bFpm7dU1ffexWMuqarXVtUbqurSqrpbVZ2oqg9V1X9X1eur6j7TY2+qqmdNH185vUviZ6bPHzid79xt+vxxVfWBaV69p6oe8tVeK6tN4FkDVXVpkj9O8sok90nyB0l+6DxP+ZEkj01y7zHGF+/qAWOMpyX5SJLrxxj3HGO8qKquTHIqyfOn/Tw7yRuq6n5bnnosyTOSfGOSS6fHJMmjpl/vPW3vvdPnD09yS5L7JnlRkt+tqtr2iwdWyjLmUVVdkuTNSW7OYrb8XJJXV9XWS6WfmuTXknx9kg8muWF67j2SvC3Ja6bnPjXJS6rqwRf4EoEV1ORcCFhhSzqHeVCSn03ysDHG5UkenWTzrMd87bTfzyd5yhjjC1mc6/xgkmuTXJHkU0lePD3ldJLrpo+vTfIvufPPWtcmeecY4/aqemiS30vyU0m+IclLk7ypqi67gNfKihB41sM1SS5J8htjjP8bY/xhkr86z+N/a4xxZozxuR3u58eSvHWM8dYxxu1jjLcleX+SH9jymJePMf552vbrk3zHV9nmh8cYLxtjfCnJTUm+Ocn9d7guYHUsYx5dk+SeSW4cY3xhjPFnSd6SxYnVHd44xvjL6QTr1blz9jwuyeYY4+VjjC+OMf4myRuSPPnCXh6wojqfCwGrYRlz5EtJLkvy4Kq6ZIyxOcb40Javf12SP0nyoSTPmP5MlCTPTPIrY4yPjjE+n+RkkidNVwqdTvKI6SqdR2XxQ/Lvnp537fT1JPnJJC8dY7xvjPGlMcZNWUSkay7gtbIiBJ71cEWSj40xxpbf+/B5Hn/mAvdzVZInT5f43VpVtyZ5RBZR5g6f2PLxZ7P4Q9n5fPnxY4zPTh9+tecAq2sZ8+iKJGfGGLeftc0rt3x+rtlzVZKHnzW3fjTJN21jv0Afnc+FgNWw53NkjPHBJL+QRaD5j6p6XVVdseUh1yR5SBY/xNq636uSvHHLnPnHLGLR/adA9Jks4vEjs/ih18enq4W2Bp6rkjzrrHn1gOl17vS1siIEnvXwb0muPOutTYfO8/hxnq+d73FnkrxyjHHvLf/cY4xx4wVsC9ifljGPPp7kAXe8n3zLNj+2jeeeSXL6rLl1zzHGT2/juUAfHc6FgNW2lDkyxnjNGOMRWQSXkeSFW758c5IXJHl7VW19F8OZJI85a9YcHGPcce5zOsmTklw6/d7pJMezeKv6B7Zs44aztnH3McZrL+C1siIEnvXw3iRfTPLz0w26npjkO/dgu/+e5Fu3fP6qJNdX1aOr6muq6uB0w7Fv2ca2/jPJ7WdtD9h/ljGP3pfFT8GfM23zuiTXJ3ndNp77liTfVlVPm557SVU9rKqu3uWagNXS4VwIWG17Pkeq6kFV9T1VdVmS/03yuSz+TPRlY4wXZXGvwLdvuSn77yS5oaqumrZzv6p6wpannc7i3j53/EU275g+f9eWt3m9LMkzq+rhtXCPqnpsVV2+jNfKxSHwrIHpRlxPTPL0JJ9M8sNJ/mgPNv2CJL86XdL37DHGmSRPSPLLWQSbM0l+Mdv4PpvefnVDkndP27tmD9YHrJhlzKNpm9cneUyS/0rykiQ/Psb4p20899NJvj+Lmyt/PIu3Trwwi/fDA/tEh3MhYLUtaY5cluTGLM5fPpHFzdd/6S72/bwsbnr8p9PflvWbSd6U5Oaq+nSSv8jiL6e5w+kkl+fOwPOuJHff8nnGGO9P8hNJfjuLmzR/cHpty5yZLFl95dvqAAAAAOjGTxMAAAAAmjsw9wJYXVV1KMk/nOPLDx5jfORirgdYX+YRMAezB9gtc4SLyVu0AAAAAJrzFi0AAACA5gQeAAAAgOYEHgAAAIDmBB4AAACA5gQeAAAAgOYEHgAAAIDmBB4AAACA5gQeAAAAgOYEHgAAAIDmBB4AAACA5gQeAAAAgOYOLGWrJ+81tvvQo0cOLWUJwMW1cXyj5tz/4ROntj13YFVsHjy2lO12/X/r3HNkp8ydOy3re5nt6XrMr4rus8fxd26ODVbVsuaOK3gAAAAAmhN4AAAAAJoTeAAAAACaE3gAAAAAmhN4AAAAAJoTeAAAAACaE3gAAAAAmhN4AAAAAJoTeAAAAACaE3gAAAAAmhN4AAAAAJoTeAAAAACaE3gAAAAAmhN4AAAAAJoTeAAAAACaE3gAAAAAmqsxxtxrAAAAAGAXXMEDAAAA0JzAAwAAANCcwAMAAADQnMADAAAA0NyBZWz08IlT7twMa2bzxsfWnPvfztzZPHjsYiwF2jl65NDcS0iSbBzfmHWO7JTznfMzc3duVY7FdWP27L0Ox7/jjTkta+64ggcAAACgOYEHAAAAoDmBBwAAAKA5gQcAAACgOYEHAAAAoDmBBwAAAKA5gQcAAACgOYEHAAAAoDmBBwAAAKA5gQcAAACgOYEHAAAAoDmBBwAAAKA5gQcAAACgOYEHAAAAoDmBBwAAAKA5gQcAAACguRpjzL0GAAAAAHbBFTwAAAAAzQk8AAAAAM0JPAAAAADNHVjGRo/edNSNfWDNbBzfqDn3f/jEKXOnmc2Dx5ay3aNHDi1luyzf3HNkp8yd5R3HXZg3+8N+nD1djk3HEOtqWXPHFTwAAAAAzQk8AAAAAM0JPAAAAADNCTwAAAAAzQk8AAAAAM0JPAAAAADNCTwAAAAAzQk8AAAAAM0JPAAAAADNCTwAAAAAzQk8AAAAAM0JPAAAAADNCTwAAAAAzQk8AAAAAM0JPAAAAADNCTwAAAAAzQk8AAAAAM3VGGPuNQAAAACwC67gAQAAAGhO4AEAAABoTuABAAAAaO7AMjZ6+MQpN/bZI5sHj829BNiek7fVnLvvOHcc38zl6JFDcy/hLm0c35h1juxUx7mzrszbnVnVGbEsZs/eWtXjbd2+r1lty5o7ruABAAAAaE7gAQAAAGhO4AEAAABoTuABAAAAaE7gAQAAAGhO4AEAAABoTuABAAAAaE7gAQAAAGhO4AEAAABoTuABAAAAaE7gAQAAAGhO4AEAAABoTuABAAAAaE7gAQAAAGhO4AEAAABoTuABAAAAaK7GGHOvAQAAAIBdcAUPAAAAQHMCDwAAAEBzAg8AAABAcwIPAAAAQHMHlrHRwydOuXMzzGTz4LF5dnzytppnxwv7Ye7M9t+Ofe/okUNzL2FbNo5vzDpHdmo/zJ11Yb7uTpcZcqHMnp7mPK73+zHB8i1r7riCBwAAAKA5gQcAAACgOYEHAAAAoDmBBwAAAKA5gQcAAACgOYEHAAAAoDmBBwAAAKA5gQcAAACgOYEHAAAAoDmBBwAAAKA5gQcAAACgOYEHAAAAoDmBBwAAAKA5gQcAAACgOYEHAAAAoDmBBwAAAKC5GmPMvQYAAAAAdsEVPAAAAADNCTwAAAAAzQk8AAAAAM0dWMZGD5845cY+K2jz4LG5l8B+dvK2mnP3qzJ3HGfr6eiRQ3MvYV/YOL4x6xzZqVWZO+yMOb039tPcM3uWZ5WOt/30PUt/y5o7ruABAAAAaE7gAQAAAGhO4AEAAABoTuABAAAAaE7gAQAAAGhO4AEAAABoTuABAAAAaE7gAQAAAGhO4AEAAABoTuABAAAAaE7gAQAAAGhO4AEAAABoTuABAAAAaE7gAQAAAGhO4AEAAABoTuABAAAAaE7gAQAAAGiuxhhzrwEAAACAXXAFDwAAAEBzAg8AAABAcwIPAAAAQHMHlrLVk/dyYx9YNydvqzl3f/jEKXOHPbV58Ni2Hnf0yKElr2R9bBzfmHWO7JS508d2j2fObT/POrNnd7odX/v5e5k+ljV3XMEDAAAA0JzAAwAAANCcwAMAAADQnMADAAAA0JzAAwAAANCcwAMAAADQnMADAAAA0JzAAwAAANCcwAMAAADQnMADAAAA0JzAAwAAANCcwAMAAADQnMADAAAA0JzAAwAAANCcwAMAAADQnMADAAAA0FyNMeZeAwAAAAC74AoeAAAAgOYEHgAAAIDmBB4AAACA5gQeAAAAgOYOLGOjh0+ccudmWCOXX30iG8c3as41mDvnt3nw2NxLWIqjRw7NvQT20NxzZKfMndWzX2fduZiBe2OdZ89+O2YcE3SxrLnjCh4AAACA5gQeAAAAgOYEHgAAAIDmBB4AAACA5gQeAAAAgOYEHgAAAIDmBB4AAACA5gQeAAAAgOYEHgAAAIDmBB4AAACA5gQeAAAAgOYEHgAAAIDmBB4AAACA5gQeAAAAgOYEHgAAAIDmBB4AAACA5mqMMfcaAAAAANgFV/AAAAAANCfwAAAAADQn8AAAAAA0d2ApWz15Lzf2gXVz8raac/eHT5wyd/bQ5sFjcy+BbTp65NDcS9gzG8c3Zp0jO2Xu3MnMWB37aSZcLOs+exy/u+OY40Isa+64ggcAAACgOYEHAAAAoDmBBwAAAKA5gQcAAACgOYEHAAAAoDmBBwAAAKA5gQcAAACgOYEHAAAAoDmBBwAAAKA5gQcAAACgOYEHAAAAoDmBBwAAAKA5gQcAAACgOYEHAAAAoDmBBwAAAKA5gQcAAACgOYEHAAAAoLkaY8y9BgAAAAB2wRU8AAAAAM0JPAAAAADNCTwAAAAAzR1YxkYPnzjlxj6wRi6/+kQ2jm/UnGtYt7mzefDY3Eu4aI4eOTT3ErhI5p4jO7Vuc2dO6zTzdsJ83Btmz1datePN9zn70bLmjit4AAAAAJoTeAAAAACaE3gAAAAAmhN4AAAAAJoTeAAAAACaE3gAAAAAmhN4AAAAAJoTeAAAAACaE3gAAAAAmhN4AAAAAJoTeAAAAACaE3gAAAAAmhN4AAAAAJoTeAAAAACaE3gAAAAAmhN4AAAAAJqrMcbcawAAAABgF1zBAwAAANCcwAMAAADQnMADAAAA0JzAAwAAANCcwAMAAADQnMADAAAA0JzAAwAAANCcwAMAAADQnMADAAAA0JzAAwAAANCcwAMAAADQnMADAAAA0JzAAwAAANCcwAMAAADQnMADAAAA0JzAAwAAANCcwAMAAADQnMADAAAA0JzAAwAAANCcwAMAAADQnMADAAAA0JzAAwAAANCcwAMAAADQnMADAAAA0JzAAwAAANCcwAMAAADQnMADAAAA0JzAAwAAANCcwAMAAADQnMADAAAA0JzAAwAAANCcwAMAAADQnMADAAAA0JzAAwAAANDc/wPaTSwuHMc3sgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1152x288 with 4 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sample_count = 10\n",
    "\n",
    "fig, axs = plt.subplots(1, 4, figsize=(16, 4))\n",
    "\n",
    "names = ['dir_tenth', 'dir_one', 'dir_ten', 'dir_skewed']\n",
    "\n",
    "locations = np.arange(sample_count)\n",
    "for plot_i, name in enumerate(names):\n",
    "    examples = trace[name][:sample_count]\n",
    "    left_offset = np.zeros(sample_count)\n",
    "    for i in range(values):\n",
    "        vals = examples[:, i]\n",
    "        axs[plot_i].barh(\n",
    "            locations,\n",
    "            vals,\n",
    "            left=left_offset,\n",
    "        )\n",
    "        left_offset += vals\n",
    "    axs[plot_i].set_title(name)\n",
    "    axs[plot_i].axis('off')\n",
    "    \n",
    "maybe_save_plot('2018-05-26-dist')    \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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}