Stick breaking process example

stick-breaking.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Stick-breaking process\n",
    "\n",
    "Code accompanying [my answer](https://stats.stackexchange.com/a/469770/35989) describing [stick-breaking process](https://en.wikipedia.org/wiki/Dirichlet_process#The_stick-breaking_process). The algorithm is described by Frigyik et al (2010) in [*Introduction to the Dirichlet Distribution and Related Processes*](http://mayagupta.org/publications/FrigyikKapilaGuptaIntroToDirichlet.pdf):\n",
    "\n",
    "> **Step 1:** Simulate $u_1 \\sim \\mathcal{B}(\\alpha_1, \\sum_{i=2}^k \\alpha_i)$,  and  set $q_1=u_1$. This is the first piece of the stick.\n",
    "> The remaining piece has length $1-u_1$.  \n",
    "> **Step 2:** For $2 \\le j \\le k-1$, if $j-1$ pieces, with lengths $u_1,u_2,\\dots,u_{j-1}$, have been broken off, the length of the\n",
    "> remaining stick is $\\prod_{i=1}^{j-1} (1-u_i)$.  We simulate $u_j\n",
    " \\sim \\mathcal{B}(\\alpha_j, \\sum_{i=j+1}^k \\alpha_i)$ and set $q_j =\n",
    " u_j  \\prod_{i=1}^{j-1} (1-u_i)$. The length of the remaining part of\n",
    "> the stick is $\\prod_{i=1}^{j-1} (1-u_i) - u_j \\prod_{i=1}^{j-1} (1 -\n",
    " u_i) =\\prod_{i=1}^j (1-u_i)$.\n",
    "> \n",
    "> **Step 3:** The length of the remaining piece is $q_k$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.stats as sp\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.20677747, 0.32461341, 0.24786725, 0.22074187])"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def stick(α):\n",
    "    k = len(α)\n",
    "    u = sp.beta(a=α[:k-1], b=[np.sum(α[i:]) for i in range(1, k)]).rvs(size=k-1)\n",
    "    q = np.zeros(k)\n",
    "    q[0] = u[0]\n",
    "    for i in range(1, k-1):\n",
    "        q[i] = u[i] * np.prod(1 - u[:i])\n",
    "    q[k-1] = np.prod(1 - u[:k])\n",
    "    return q\n",
    "\n",
    "stick([10,10,10,10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjwAAAI9CAYAAAA+U6I/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3dfdBkV10n8O8JIQwSWQhksUTIBBQSQ4Hirkoha0wo38ot0egE0JRTa3RhAxsXSteokWhcJdZSxK1C3hSjkOCAQVnXRNgVsyECq7guDNkESpiQClkRRSAJCUE4+0f3SPMwz/R5Jrdv3z7P51P1VM3z3Jdzbk/37W//zrm3S601AAA9O27dHQAAWDWBBwDonsADAHRP4AEAuifwAADdE3gAgO4JPPdRKWV/KeVzpZQ7Symnr7s/rE8p5YOllHtLKa9bd19gqpwz25RS3lZKuaeUcsO6+9ILgWcY76y1nlhrvenwH0op/6GU8jellE+WUl5TSnnAdhuXUmop5a75CeDOUspvtDZcSjmplPL78+0/XEp59lHWXTzRHP45cwdtnV1KubmU8ulSyp+WUk45yrrXzV+sh9t5/w7aKaWUy0opfz//+dVSStlm3b3zx2/xmC7eQVt758fy6fmxPf0o615SSvnslrYec3h5rfWxSX65tW3Yxb7onFlKeUIp5S2llL8rpezo5nCllBNKKb9XSrllfi44c4fbX1pKOVhK+cdSyiU73HZfKeUd8/PHdTvc9qjHXGs9K8lzdrJPjk7gWYFSynck+ekkZyfZm+QxSX5hyWZPmp8ATqy1nr+D5l6W5N4kj0jyQ0leXko54yjrv3OhnRNrrde1NFJKeXiSNyW5OMlJSd6d5MCSzZ630M7jW9qZ+/Ekz0jypCRPTPI9Sf7tkm0estDWpTto6/VJ/irJw5L8bJLfK6WcfJT1D2x5/D60g7aAI/tskjck+dFj3P6GJD+c5G+OYdu/TvJTSf7oGLb9eJLLk7z4GLa9r8fMDnUdeEopTy6l/Hkp5VOllDeXUh45r4IcsVowoB9J8pu11htrrf+Q5NIk+4dupJTyoCTnJLm41npnrfWGJP81yXlDt5Xk+5PcWGt9Y631niSXJHlSKeW0FbT1I0leUmu9rdb6kSQvyWoev8cleXKSF9Va7661Xp3kYGaPKew66zpn1lrfX2v9zSQ3HsO299ZaL5+f/z53DNv/dq312iR3HMO2/6PW+oYktx/Dtsd8zBybbgNPKeWBSf4wyauSPDzJ8ZlVKK6t23yfRinlE0f5+ekdNH9Gkvcs/P6eJI8opTzsKNtcPx8Ce1MpZW9jO49L8rla6we2tHW0Cs/Xz0uoHyilXFxKOb6xrS86plrrXUk+uKStX5m39Wc7LDMf6fE7WjtJ8uFSym2llN+aV6Na2/lQrXXxRLesrX9dSvl4KeXGUspzG9uByVvzORNWrtvAk+SbM3vBvqbWem+S30ryjUn+YLsNaq0POcrPTkqWJyb55MLvh//95dus/62ZDX2dltknhf/WGES2tnO4re3auT7JE5L888yqGM9K8pMN7RxLW/8xs6G8R2Z2Av3DUspjj7GtTyY5cZtPmX+X5F8mOSXJN8z7c+UxtnO4re2O6Q1JTk9ycpIfS/LzpZRnNbYFU7fOcyasXM+B5xFJPlxr/fz89/ckuTPJ20Zo+84kD174/fC/j1gyrbVePy/LfiLJhUlOzeyNdaftHG5ru3Y+VGs9VGv9fK31YJJfTPIDDe0cS1v/q9Z6R631M7XW307yZ0m++xjbenCSO4/0KXM+lPfuWus/1lo/muR5Sb69lLK1ry3tHG5ru2P6v7XW22utn6u1viPJr6X98YOpW+c5E1au58Dz0SQPXfj9tCQls4liR7Tl6putPz+zg7ZvzGzC7WFPSvLRWuvfN25f531d5gNJji+lfM2WtlrHhFvbSbYc03z+0GPHaCs7P6Y0tnVjkseUUhYrOqt6/GDq1nnOhJXrOfC8K8k/K6U8df77czN78T55uw22XH2z9Wcnlxv/TpIfLaV8bSnloUl+LskVR1qxlHJGKeXrSin3K6WcmNkE3Y8kOXy55pnbXaY5n0fzpiS/WEp50PxYvzfJa7dp67tKKY+Y//u0zK64evPC8itKKUfsZ5LfT/KEUso5pZQ9SX4+yXtrrTcfoZ2HlFK+o5Syp5RyfCnlh5L8qyRvWVjnaJeP/k6SF8wnTH5lkhdm+8fvm0opjy+lHDefI/VfklxXa/3kfPn+UsotR9p2Pvfp/yR50byv35fZVWFXb9PW95ZSHlpmvjHJv8/C4wcbbm3nzPlrak+SE+a/7ykLt/JYcm5KKeUB8+2T5IT59mW+bNtzwHz5/efbHpfZB8g9pZT7zZcdvu3F3m22vd982+OTHDff9v4Ly28ppew/lmNmeN0Gnlrr3Ul+MMmvlVLel+RTmV299PqFF/Sq2v7jJL+a5E+TfHj+86LDy0sp1y58+nlEZpd3fyrJhzKby/M9tdbDn6oeleSdR2nu3yV5YJK/zewS6+fWWm+ct/Po+SetR8/XPTvJe0spdyW5JrOwtHhSelRmQ09HOqaPZTbv5z8l+Yck35TkmQvH9DOllGvnv94/yS8l+Vhmc2yen+QZtdb3z9f9qsxK5Qe3OaZXZjZ58mCS92V2uegrF9q6cR6iktk8oT/ObBjqfUk+k9ncpKXHNPfMJP9ifkwvTvID82NNKeVppZQ7t6z71/O2fifJZfPhOth46zxnZjYH7+58obp6d5LFe3ctex2/f77NIzP7YHX3fJ8t2756vv6zMrs1xd35wpWuj8rs/P2RbbY9b77+y5M8bf7vVyez+wNldruLd22z7bJjZmBlm8n3NCqlnJfZm/G9SZ6yePPBgfb/G0neWGt9y9KV71s7J2Q2Zv/EhbC1qrZ+OMkZtdaLVtnOvK23Jrlw6P+Xbdp6f2Yn3DfUWv/NqtuDTbTTc+Z9PTfdl3NAKeXnknys1vrKpSt/6bbfkuSCWusxXdhQSvnvmU0k//Na69nHsg++mMADAHSv2yEtAIDDBB4AoHsCDwDQPYEHAOiewAMAdE/gAQC6J/AAAN0TeACA7gk8AED3BB4AoHsCDwDQPYEHAOiewAMAdE/gAQC6J/AAAN07fsidvew5b6vL1rn03JOW7ufiAx9fuk7LfobcV8t+Ttv3Y019uvkNrx6tT0PtZ8h9tTxOLY/RkH3q+Xl5wSvOKk2d2uWGOn8l03vOjf0a77lPm/ocGLtPQ57nhzqHDRp4gHG94ikXLl3nghwcoScA0zbJwHPWdRcsXefSc18/Qk+YuivrOUvXuTltnyIAerbb31snGXjYXC0vqNv3jdARAFhg0jIA0D2BBwDonsADAHRP4AEAuifwAADdc5UWbLCDh25ddxe60XKF4Wn77m3al1shwPSo8AAA3VPhAYBdYN9Fu/stf3cfPRvv7deft3Sdk0foBwDTZkgLAOiewAMAdM+QFhvt/HvOXrrOm/PZEXoCwJQJPAAQ3ybeO0NaAED3BB4AoHuGtADW4Mp6ztJ1xr5j8xT7BEMReACA0bXMmZq5aZD2DGkBAN1T4QEgiTuX06alMnP7vhE6skMqPABA9wQeAKB7hrSYpObJbM941Wo7AkAXBB4A2GAttxNIkp/Il624J9NmSAsA6J4KD8AauCIKxiXwTJQ7nvbNlxQCjGtjA0/rmKVQAMDYfKiZno0NPEzTvouWP6Uuz70j9AQAvkDgieEjAOidwAOwwXxggzYuSwcAuifwAADdE3gAgO6ZwwMAa2D+1bgmGXhc2gzQxh2bN5f/u3EZ0gIAuifwAADdm+SQFgzJLd6BobR+rVFywkr7wc4JPAA0e8VTLmxY67Ur7wfslCEtAKB7Ag8A0D1DWgBrcP49Zy9d58357Ag9Ycpa5iDevm+EjnRAhQcA6J7AAwB0T+ABALon8AAA3TNpmaZJcck0b8534NBlS9f57hH6AcC0qfAAAN0TeACA7hnSAmBQLd83dXNePUJP2Km27wrbzO8JU+EBALqnwjNRb7/+vKXrnDxCPwCgByo8AED3BB4AoHujD2m1TIj6iXzZCD2BnXnFUy5sWOu1K+8HsHM9T8aljTk8E+WblIez76JhnuYHD926dJ2vePQgTUH3NvmGp2wmQ1oAQPdUeMjpz7x93V2AXafla1H2PPQFI/QEdgcVHgCgewIPANA9gQcA6J45PAAMquVO8efmwAg9gS8QeGBAbff6SNzvA2BcAs8Ga7mPhXtYAGyulvuIXZ57R+jJzkyxymcODwDQPRUeyE6GogBYpuXbAu4YucIj8EDayq9n56Uj9ASAVRB4IG2fRrJH4IGxtVRfb86rl64zxTkljEvgAQCS9B0MBR4AIMk0594MReCBAT3+rVc0rXdHfny1HWGtWm4Zcc2THjtCT4DDXJYOAHRPhQeA7vU8VEMbFR4AoHsCDwDQPUNaAAzK8BFTpMIDAHRP4AEAumdIC6BzLfcFSpIDzzx3xT2B9RF4AIDR7buoLYIcHKg9gQegUcv3DCWb+11D0DOBB4CN1jRk94xXDdKW0Lu5BJ41eMVTLly6zvfe9OuDtHVlPWeQ/QDAJhN41uDgoVuXrnPbCP3YqZbw9BP5shF6AgA7I/DQrKmUe+rVg7TVOplt/zWDNAdNz7lzDo3QEWAlBB4ARtf6oeby3LvinrBbCDwAnWsNF9fetPwrIa7OMFXc3qkYTo87LQMA3VPhYVAtE7L/5FEPH6EnwE4dOHTZ8pVOX30/YBVUeACA7gk8AED3BB4AoHsCDwDQPYEHAOieq7To3hS/yqP1vigADGP0s+6YX08wVXvvuWrpOjcM1FbL4312XjpQawAwTaMHnvPvcSdPAGBc6uoDaqpeAQCjE3iYpJZ5N0nykpyy4p4A0ANXaQEA3RN4AIDuCTwAQPcEHgCgeyYtA0Dabgh6zqEROsJKqPAAAN0TeACA7gk8AED3BB4AoHsCDwDQPYEHAOiewAMAdE/gAQC6J/AAAN0TeACA7gk8AED3BB4AoHsCDwDQPYEHAOiewAMAdO/4dXdg1c667oKl69y+b4SO8E8e/9Yrlq+053tW3g8Ado/uA8++i5Yf4uW5d4SeAADr0n3gAWBzvf3685auc24OjNATNp05PABA9wQeAKB7Ag8A0D2BBwDonknLAECzlqufr71phI7skAoPANC9ja3wtFyqmCQ59erVdgSAyTtw6LLlK52++n6wPio8AED3NrbCA6323nPV0nVuGKEfAKyPwAPQ6Px7zm5a7w53/oXJEXgYXctY+guNpQMwIIEHAJIcPHTr0nUuWX03WJFJBh5POgBgSJMMPAC0aZlXdHWmd3uOlg+2iQ+3DGdjA0/r5MGWF3rLPX1ObmqNobRcWZUkz8/LV9wTAHrgPjwAQPcEHgCgewIPANC9jZ3DA5usZcLmVzx6hI4A7BICT6NXPOXCpeucc+icEXrCuviKChhO64UJ+/f8xYp7wm4h8HSu6Wq2PS9dfUcAWKuWyvJtI/RjXczhAQC6p8IDABPVUpX5k0c9fISebD4VHgCgeyo8sAatEzYBGIbAAwAbrOXrkZLk7OzuC1QMaQEA3VPh2WD7Llr+33ftTSN0BGDBgUOXLV/p1OeuviOwQIUHAOieCg/N3MQQgE0l8AA0ahqqSfLdK+4HsHOGtACA7gk8AED3DGl1rqUE/8LTR+gIAKyRCg8A0L1BKzxnXXfB8pWe8aohmwQAWMqQFsSN0oDxudXHuAxpAQDdU+EBmCjTBGA4KjwAQPdUeBjU3nuuWrrO8/PyEXoCAF8g8AAMbN9Fy0+t+68ZoSPAPzGkBQB0T+ABALon8AAA3Rt9Do/vdgKANi0Xguzf8xcj9GTzqfAAAN0TeACA7gk8AED3BB4AoHtuPNjo4KFbl65zyeq7AbBWLefCl+SUEXoCOyPwAABJ2q6kzoZeSW1ICwDongoPAAysZejvthH6wRcIPAM6/56zm9a7Ip9acU8AgEWGtACA7gk8AED3DGkBQKPWqQvZ89LVdoQdU+EBALon8AAA3TOkBRN1ZT1n3V0A+BKbesm9Cg8A0D2BBwDonsADAHTPHB6y956rmta7YcX9AIBVUeEBALqnwjOgA4cua1vx1OeutiMAwBcReAA613IZcZK8JKesuCewPgIPTNTbrz9v6Tpnx+3rAVoIPABMVst3V92RAyP0hE0n8MCAzOMCmCaBByaq6VuZfSNz1/ZdtPwUvf+a5SH7hacP0RsYVuvcsqEIPDRrql6oXAAwQe7DAwB0T4UHgI3WNvQ3QkeYNIEHgI3WMhdk7HsMtXxlj6/rGZfAE5NDgWm+aQ7Vp9bvy3t+Xj7Ivlr2w+ba1DBnDg8A0D2BBwDonsADAHRP4AEAumfSMgCT1XLDU3eSpoUKDwDQPRUeADaaS+VpIfBMVNP3VinjAkATgQcANljTzXOTXX8DXYEHABhd6x3AbxmoPZOWAYDuCTwAQPcGHdLad9Hy3e2/ZsgWl2v58r3bRugHAKzTbr+aTYUHAOieScuwwVo+sd2y+m6wIrv9EzkMSeAByLDh4o6bXtywlqACY+o+8LScxG4YoR8AwPqYwwMAdG+SFZ6m0nLLVy8kyanPvY+9AQA23SQDDwCMzRSIvhnSAgC6p8IDAAzqQMu0k5GnnKjwAADdE3gAgO4Z0gKARk1DNYkrhCdo0MDT8kWdL8kpQzYJABtpivNcejZo4PG9L+MSMAGgjTk8AED3BB4AoHsmLTea4nDdFPvEcIzvAwxHhQcA6J4KT3ySBmBzuVS+Tam1rrsPAAArZUgLAOiewAMAdE/gAQC6J/AAAN0TeACA7gk8AED3BB4AoHsCDwDQPYEHAOiewAMAdE/gAQC6J/AAAN0TeACA7gk8AED3BB4AoHsCDwDQPYEHAOiewAMAdE/gAQC6J/AAAN0TeACA7gk8AED3BB4AoHsCDwDQPYEHAOiewAMAdE/gAQC6J/AAAN0TeACA7gk8AED3BB4AoHsCDwDQPYEHAOiewAMAdE/gGUEpZX8p5XOllDtLKaevuz/sTCnlilLK3aWU29bdF2jhnMMqlFJ+oZRyVymlllKOX3d/dkrgGc87a60n1lpvSpJSyhNKKW8ppfxdKaVuXbmUclIp5ffnT64Pl1Ke3dpQKeXbSil/Wkr5ZCnlliMs3ztf/ulSys2llKfvYN/7SinvmG973RGWf10p5S/ny/+ylPJ1O9j380op7y6lfKaUcsURlp897++n5/0/ZQf7vrSUcrCU8o+llEuOsPzZ88f5rlLKH5RSTjq8rNa6P8l3tbYFEzHoOWddr79SygNKKa8ppXyqlPI3pZQX7KBdxzzgMddaX5TkjNa+TI3Asz6fTfKGJD+6zfKXJbk3ySOS/FCSl5dSWp9odyV5TZKf3Gb565P8VZKHJfnZJL9XSjm5cd8fT3J5khdvXVBKOSHJm5O8LslDk/x2kjfP/97i9iS/NO/71n0/PMmbklyc5KQk705yoHG/SfLXSX4qyR8dYd9nJHllkvMye7w/neTXd7Bv2ATHfM5Z8+vvkiRfk+SUJN+W5KdKKd/Z2K5j/lKrPOZpq7Xu6p8kT07y50k+ldmb9SOTfDhJGbCN/Ulu2GbZV8/+G77obw/K7An5uIW/vTbJi3fY7tOT3LLlb49L8pkkX77wt7cnec4O931+kuu2/O3bk3xk8bFLcmuS79zhvn8pyRVb/vbjSd6x5TG6O8lpO9z365JcsuVvv5zkqoXfHzt//BcfozOT3Lbu56ufzf/ZxHPOOl9/83PKty8svzTJ7+6wXcc80DEn2ZukJjl+qOfrWD+7usJTSnlgkj9M8qokD09yfGbp9to6/589wjafOMrPTw/Utccl+Vyt9QMLf3tPhiklnpHkQ7XWO1a07/dueezeO+C+33P4l1rrXUk+uKJ9fzDzk8IA+4Z/ssHnnLW8/kopD03ylYvLM9z5yjF/6b5Xecxrt3GTjgb2zZk9Bq+ptX6+lPJbSd6Y5EXbbVBrfcgI/ToxySe3/O2TSb58hft+5Ar3PVS/P7bCfa+q37BoU88563r9nbjw+5jtHl7umDs6F+7qCk9mY5gfrrV+fv77e5LcmeRt6+tSMu/Dg7f87cFJ7jjCuvY97X3Dok0956zr9Xfnwu9jttuyfFVt93rMa7fbA89HM5tce9hpSUpmk76OqMwu89zu52cG6tcHkhxfSvmahb89KcmNA+z7xiSPKaUsJvYh9/3EUkpZ+NsTB9z3kw7/Ukp5UGbj3qvY92OSPCCz/wcY0qaec9by+qu1/kOS/7e4PMOdrxzzl+57lce8fuueRLTOnyQPTPK3SZ46//2azK5w+oaB29mfLRMIMzvJ7UnytZlNANuT5AELy383s6upHpTkqZmVFc9YWF6TnLlNe8fN9/ddmU2G3JPkhIXl70ryn+d//74kn0hy8nzZ3vm+926z7/vNt3tOkuvn/77/fNkJ8/YuzOzF+7z57yfMl5+ZLZPotuz7+Pn+fiWziXR7Mp8Yl+Tk+WNwzvzvlyV515bH+Jaj7Pv+8+2uymxS9J4k95svOyOzCaRPmz/er8uWCYIxadnPAD+bes5Z5+svsytC/2dmQfG0zMLAdy4svyXJ/m3adcwDHvN8nb3Z0EnLa+/Aun+SfGtml969b/5E+P7MUvBTB2xjf7705HP4SbP4c8vC8pOS/EFmJ8Nbkzx7YdlXZVZifNg27Z15hH1ft6Xt6zKbff/+JE9fWPa0+Yvp/kc5lq37vmJh+dcn+cv5vv93kq9fWHZeFq4AOMK+LznCvi9ZWP70JDfP931dFkJZZpdRXnmUfV9xhH3vX1j+7PnjfFdmV86cdITHVODxc59/NvGcM1++ltdfZh+eXpNZQPhokhcsLDshs3PhEa+ccszDHvOW/W9c4CnzA2CFSinnZXbPhXuTPKXObwR2H/b3w5kl8ouG6N+Wff9cko/VWl+5gn3/RpI31lrfsoJ9vzXJhff1sd1m37+Z5AeT/G2t9auH3j8MbehzTkN7K3v9LWn3W5JcUGt91pjtztvejcf8oiQvyCyQPajW+rmx+3BfCDwAQPd2+6RlAGAXEHgAgO4JPABA9wQeAKB7Ag8A0D2BBwDonsADAHRP4AEAuifwAADdE3gAgO4JPABA9wQeAKB7Ag8A0D2BBwDo3vFD7uxlz3lbXbbOpeeetHQ/Fx/4+NJ1WvYz5L567lPLfvRpc58DF7zirNLUqV1uqPNXMr3n3Ka+nqbYp019DkyxT63PgaHOYSo8AED3BB4AoHsCDwDQPYEHAOiewAMAdE/gAQC6J/AAAN0TeACA7gk8AED3BB4AoHsCDwDQPYEHAOiewAMAdE/gAQC6J/AAAN0TeACA7gk8AED3BB4AoHsCDwDQPYEHAOiewAMAdE/gAQC6J/AAAN0TeACA7gk8AED3BB4AoHsCDwDQPYEHAOiewAMAdE/gAQC6J/AAAN0TeACA7h2/7g5sirOuu2DpOpee+/oRerI7tDzeicccgDaDBh6hAACYIhUeWAMfDgDGJfBAkivrOUvXuTmvHqEnAKyCwJO2N7vkhJX3AwBYDVdpAQDdU+GZqKGGWAzVAIDAQ1wCDkD/DGkBAN3b2ApP20RjwzXA5nL7AhjOxgYeANr0Pmw9VDAcO2COOVczSb7yuuVXG2/qc6CFwAMAE9USwm7fN0JHOjB64HHVEACt3CdtM03xvb77Cs+Y6di8orbH+7R99zbsyQmM6RnyNe6NfFweb7oPPAC0meKnchiKwAMDaq0A+CTJUFQuoI3AQ/dM+mNsnnN4DkyPwAMQb1BDMjTWt019rQg8AECSvodIBR4ARtd6M8SpVQrM09tcAg8AMLrW0JvcNEh7Ag9MVMsnyZZbxQ91sgDYZAIPG63n8WYAhnPcujsAALBqAg8A0D1DWtBoU+89AYAKDwCwC6jwANBMpZNNJfCsgRMGAIzLkBYA0D0VHpq55w0Am0qFBwDonsADAHRP4AEAuifwAADdE3gAgO4JPABA9yZ5Wbob8wEAQ1LhAQC6N8kKD21UwgCgjQoPANA9FR4mad9FbU/Ny3PvinsCQA82NvC8/frzmtY7NwdW3BMAYOoMaQEA3RN4AIDuCTwAQPc2dg7P2Fom0ZpAi7llANMk8OCKKAC6J/AAwBq0VIRVg4fTfeAxFDU9TcM+p169+o4AsGuYtAwAdE/gAQC61/2QFuNqGUI859AIHQGABSo8AED3VHhoZrLxcEymxxU6MC6BBzZYS3A6OEI/gD70HMQNaQEA3VPhAZgoQ5/D6blyMUVTfLxVeACA7qnwADBZU6wUsJkEHrpnWADYVE33Nrv+3KZ97fZgKPAAAM029QazAg9MlPse0WJTh3xa3jSTtjfOod6Ae68Gb2pQGYrAE28sQJumc0WmGTBaDBWenFOZIoFnQL2fDKfIiRWWG7KaApvKZekAQPdUeAAGttvnSsAUCTwAMDChd3oEHoB4g6JN61xNcwenR+CBRt4QATbXJAOPNxYAYEiTDDxs7s3EgL75QMqmEngg7ucD0LtBA0/vt+UGADaTGw8CAN0bfUjL0MFwVNQAoI0KDwDQPYEHAOiewAMAdE/gAQC6J/AAAN0TeACA7rnT8hq4NTsAu13Le2GSHByoPRUeAKB7Ag8A0D2BBwDonsADAHRP4AEAuucqrc75slYAUOEBAHYBFR4m6eChW5vWu2S13QCgEyo8AED3VHgY3fn3nL18pT0vXX1HANg1VHgAgO4JPABA9wQeAKB7Ag8A0D2BBwDonqu0AOhey9WhV2fcu85PsU89E3gA1sCbHYxL4AGYqJY7jl+y+m7smDDXt019Xgo8DKrlhXDbCP3YLVq/ggNgtxN4AJgs1aLNNMX/N4EHBtT0tRlxggYYm8AzIG9245vipwgApkfgoXkeyCWr7QYAW4w9QbjnL3cWeCAqRQC9E3jizQ4g2dxzoatDaSHwAEA29/4y5o+22djA4z8YgLF579lcGxt4Wm1qqXNTS8sATNOmvh8OpfvAA9BiU4czgDYCD8AGUw2GNgLPGuz2siL0zmscz4HpEXgaefICm6r15qLOYfRM4AFo1HqFzqbeibbFUB/+hDDGdty6O5HIsiAAAALxSURBVAAAsGoqPNDIsCbA5ppk4PHGMk09f6kcAH0zpAUAdE/gAQC6N8khLdh7z1VN692w4n4AsBqtV+oNReChe+aEAWBICwDonsADAHRP4AEAumcOT+fcOwcAVHgAgF1g9AqPisNwXH0EAG1UeACA7gk8AED3BB4AoHuu0oI1aJl/dcnquwGwawwaeHb7JNoDhy5rWu+Fp6+4IwDAFzGkBQB0z5AWTJRbOAAMR+CBRnvvuWrpOjeM0A8Ads6QFgDQPYEHAOiewAMAdE/gAQC6J/AAAN0TeACA7gk8AED33IeHZi1fneFrMwCYIhUeAKB73Vd4pnh33Cn2aSg9HxsAm6v7wAPAuAx/M0WGtACA7qnwMLqWT3859bmr7wjQBUPptBB4AGANDP2NS+Ch6dNRkjw/L19xTwBgNQYNPMqKMK6W19wtq+8GWzgX9qtpSD4qM1OkwgOQaYaUKfYJNvV5KfAArIHJ++Pqeb7MFKtOU3y8Rw88U3wQAIC+bWyFpzXR+oQE9MyHSGizsYGHNsrmQOuVmGPOu5hin+ibwEP3NnWCHdPTe2V5itUiH9oYisADA5ri5EFgfC0ftMa+t9lu//An8GSan2oAYOzg1HNFTeBpNMW0DgC0EXgAID7Y9k7g2WBenMDYnHcYSuuVercM1N4kA48XFAAwpOPW3QEAgFUTeACA7gk8AED3BB4AoHsCDwDQvUlepQX0fcdTgLEJPGw0oQCAFoa0AIDuqfDAGri5JsC4VHgAgO4JPABA9wQeAKB75vBMlKuPAGA4KjwAQPdUeCAqagC9U+EBALon8AAA3RN4AIDuCTwAQPcEHgCgewIPANA9gQcA6J7AAwB0T+ABALon8AAA3RN4AIDuCTwAQPcEHgCgewIPANA9gQcA6J7AAwB0T+ABALon8AAA3RN4AIDuCTwAQPcEHgCgewIPANA9gQcA6F6pta67DwAAK6XCAwB0T+ABALon8AAA3RN4AIDuCTwAQPcEHgCge/8fi+U2ncgIfewAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x720 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from pylab import rcParams\n",
    "rcParams['figure.figsize'] = 10, 10\n",
    "\n",
    "alphas = [0.5, 1, 10, 100]\n",
    "n_sticks = 25\n",
    "\n",
    "for p in range(4):\n",
    "    plt.subplot(2, 2, p+1)\n",
    "    α = [alphas[p] for _ in range(5)]\n",
    "\n",
    "    for s in range(n_sticks):\n",
    "        q = stick(α)\n",
    "        plt.bar(s, q[0], 1)\n",
    "        for i in range(1, len(q)):\n",
    "            plt.bar(s, q[i], 1, bottom=np.sum(q[:i]))\n",
    "        \n",
    "    plt.axis('off')\n",
    "    plt.title(f'α = {α}')\n",
    "    \n",
    "plt.show()"
   ]
  }
 ],
 "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}