mcmc_intuition.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": true,
    "hidePrompt": true
   },
   "outputs": [],
   "source": [
    "# A seed value for reproducibility\n",
    "set.seed(32)\n",
    "a = 2.0\n",
    "b = 1.0 / 3.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppose we have a random variable θ that follows a Gamma(a,b) distribution. Let’s say a=2 and b=1/3, where b is the rate parameter.\n",
    "\n",
    "$${\\displaystyle f(x;\\alpha ,\\beta )={\\frac {\\beta ^{\\alpha }x^{\\alpha -1}e^{-\\beta x}}{\\Gamma (\\alpha )}}\\quad {\\text{ for }}x>0{\\text{ and }}\\alpha ,\\beta >0}$$\n",
    "\n",
    "\n",
    "where $$ {\\displaystyle {\\Gamma (\\alpha )}}$$  is a complete gamma function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": true,
    "hidePrompt": true
   },
   "outputs": [],
   "source": [
    "m = 100\n",
    "theta = rgamma(n=m, shape=a, rate=b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false,
    "hidePrompt": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAAAM1BMVEUAAAAAAP9NTU1oaGh8\nfHyMjIyampqnp6eysrK9vb3Hx8fQ0NDZ2dnh4eHp6enw8PD////UNI3wAAAACXBIWXMAABJ0\nAAASdAHeZh94AAAesUlEQVR4nO3di5aivBZF4Yjaanl9/6dtQUVQUEkWsJPMb4xz2q6qLTRm\n/t4tdwEQzM29A0AKCAkQICRAgJAAAUICBAgJECAkQICQAAFCAgQICRAgJECAkAABQgIECAkQ\nICRAgJAAAUICBAgJECAkQICQAAFCAgQICRAgJECAkAABQgIECAkQICRAgJAAAUICBAgJECAk\nQICQAAFCAgQICRAgJECAkAABQgIECAkQICRAgJAAAUICBAhJzznXPvX8QtO/SXZmWzjX3NLt\ndPcetX8GQxCS3k8hHYpJDv32uu1GSI+tfgppoj1LDMdM76eQPl8nyCydO3Zs9dPWJ9qzxHDM\n9N5C+vxDE+1L66+EpMYx0+u7RjpvV9dT67/bV9zjx/b/yltf+/vI6fq31a4xeVq6zfXU3/p6\nerk5Pc5vt3TLw+WyK9zq0N586/xaVTy3Wv7/4bo3/06375w3hSuqM2/sWXOT+IKQ9HpCOhX3\nRbpqLdfV/fS6mjjcf+Q5uawGHj/lDrev3v5+2tRfqzXPr9nr5SWk22hRZfLYtUPzZ1qbxBeE\npNcT0vWK4npldL6uz11jua4fy/VWUlH/9THpyrHdNZ3z5bJp9XHNoNngTev8PoV096+51aLx\nM+1N4gtC0nNN9y/c/r/8z//5eiVTf+2yv/65O19v9V3/vN4a+7su5vKP4jlZrubyMYNT65yu\nX72udLc8Vn88t/1yfh/uIxWHqpJrO/dmzrfSHz/T3iS+4CDp9YRUxlHfFXqsz3/l9VNpU105\nrG/LvwzqMbl/Oevb/x9afzx/4OX8PoRUnu+5vlY837667h/BRxwkvZ6QtvfbUvvL82vlH9Ua\nvpyqLxSPVfv67esP/G1Wzxtml7c/6rnm+X171O71hl7RGmlsEl9wkPSeS6+91DeP1Xp6+9bj\nlHsP6fb3v2WjzM8htU4NDKk50tokvuAg6fWFdDn/3R4JW7W+VV+DFJ3XSNVfy5t6y3+746Br\npOL1m30hFa0fun+vvUl8wUHS6w2pVD3L8/za+ut9pOq7y/vXv4a0/vk+UuPPdeuemKsfbNhf\n3s8B3ThIej0hLeu79I+rinPvo3buJZL7n9+vkX541O78dr5/t8fw/p7XlefXTeILDpJeT0jX\nNb46VY85lK9UKB/DK/+sn/ZsP6PzEsmq+uF98TWk1/N7yeCx1ZdA660enj/T3iS+4CDp9d20\nezzYUP5nv3yc+nZi1Vz31TWKa72yofry/QUP5TOwh88hvZzfyzcfW30J6b7VqrHHz7Q3iS8I\nSa8vpNv9o9XtPkx5x+Tezr+i8QTTsXyt3f4tkvLLxb/j6fGChY5zf2id3+s31+27To8/z5vr\nDc/1vvUzrU3iC0Iy6Xy7I4VoEJIp7vYqneOKq4HIEJIpz4cK2i8NgnWEZEr9Vovb/X5Eg5Bs\nOW/L90EU/7g+igwhAQKEBAgQEiBASIAAIQEChAQIEBIgQEiAACEBAoQECBASIEBIgAAhAQKE\nBAgQEiBASIAAIQEChAQIEBIgQEiAACEBAoQECBASIEBIgAAhAQKEBAgQEiBASIAAIQEChAQI\nEBIgQEiAACEBAoQECBASIEBIgAAhAQKEBAgQEiAQGNJu6dx6r9kVIF6+IblqcOUqG+EOATEK\nCmnjNufL5bRxO+UuAfEJCqlw5/L02S11OwTEKCgk5xp/ATIWFNK/R0iFaneAOPmHtN7u9u7v\nevK84dEG5M4/pJvqZHFW7hIQH+97N8fjbrdeVw85bOgIueNhAkCAkAABQgIECAkQICRAIPTh\n7+ej4EDGfBPYERLw5P88UrFS7gcQNf/rkiMvDAIeAm6U7dyx/2ynuN3nQo21Y8jPBItpvJBm\nngdqhAQIEBIgQEiAgGQxfb7bTkhIHyEBAty0AwQICRAgJEDAfzEdtuvq5QHrzWGsTXxBSDDD\ndzGdl42X2nx++SohIX2+i2njir/bS+1O++Lzy1cJCenzXUxF4xWrx8+ftEpISF/Yr3Xp+ots\nE9/3YeZ5oMY1EiAQcB9pf6pOcR8J8F9Mq8ajdsuPn1lMSEhfwPNIm+p5pGK95XkkZI9XNgAC\nhAQIEBIgQEiAACEBAoQECBASIEBIgAAhAQKEBAgQEiBASIAAIQEChAQIEBIgQEiAACEBAoQE\nCBASIEBIgAAhAQKEBAgQEiBASIAAIQEChAQIEBIgQEiAACEBAoQECBASIEBIgAAhAQKEBAgQ\nEiBASIAAIQEChAQIEBIgQEiAACEBAoQECBASIEBIgAAhAQKEBAgQEiBASIAAIQEChAQIEBIg\nQEiAACEBAoQECBASIEBIgAAhAQKEBAgQEiBASIAAIQEChAQIEBIgQEiAACEBAoQECBASIEBI\ngAAhAQKEBAgQEiBASIAAIQEChAQIEBIgQEiAACEBAoQECBASIEBIgAAhAQKEBAgQEiBASIAA\nIQECCYe0KAXMA79LNaRHQ59iIiTIJBnSSzx9KRESZBIMqaOb7pQICTLJhdRz9dP1ZUKCTGoh\n9T+68J4SIUEmsZA+Pkr3+k1CgkxaIX3s6O3bhASZlEL6/KzR+08QEmQSCulrRq8/REiQSSek\nnzpq/RghQSaZkH7sqPmDhASZVEL6uaPGHSVCgkyGIdU/TEiQSSSkQR09fpyQIJNGSAM7ug8Q\nEmQyDamaICTIJBHS8I6qGUKCTAoh+XRESJDKNqTrFCFBJoGQ/Dq6zhESZOIPybejkEngRc4h\nOUqCSvQhBcTguE6CSt4hURJEYg8ppARCgkzmIVESNCIPKaiD0JuGQC37kCgJCnGHFBZB4BO6\nwFPWId0tnCfRvwMJiDqkwOuSsNfq8aJXNBDSxbskQkIt5pBC79yEvDGwPY/sEZL/ORESahGH\nFPw2CP8PT3mbR+4I6X5mgfPIHCHdzyxwHpmLN6Twt4r7feRxzzzyRkjP8wubR9YI6Xl+YfPI\nWrQhCT5O62V+cEmEhBohtc4yaB4ZizUkxUcOExJkCOnlTEPmkS9Cej3XkHlkK9KQJL+WpWN+\nUEmEhBohvZ9vwDxyRUgdZxwwj0zFGZLmV1d2zg8oiZBQI6TOs/afR54Iqfu8/eeRpShDWmjO\nmJAgQ0h95+49jxwRUu/Ze88jQ4TUe/be88hQjCE9FvpYIf1aEiGhRkgfNuA7j/wQ0qct+M4j\nOxGGVK/y8UL6rSRCQo2QPm/Dcx65IaQvG/GcR2YI6ctGPOeRmfhCei7xMUP6pSRCQo2Qvm/H\nbx5ZIaTv2/GbR1aiC6mxvscN6XtJhIQaIf2yJa955ISQftqU1zwy4r8YDtu1K603h7E20WXC\nkL6VREio+S6G89I9rUbZRKfm2iYkmOG7GDau+DtWp077wm3G2ESnSUP6UhIhoea7GAp3rE8f\nXTHGJjpNG9LnkggJNd/F4FzfX2Sb6NJa2IQEMyK7Rpo6pI8lERJqAfeR9qfq1KT3kQgJNnkv\nhlXjUbvleZRNdJg8pE8lERJqAc8jbarnkYr1drrnkdqrepKQPpRESKjF9coGQoJRhDRwo8Lt\nIyFRhfSypCcKqbckQkKNkAZvVrZ9JISQhm9XtX0kxP+VDS1jbOIdIcEq38WwmyGk1/U8WUg9\nJRESat6L4Vh8fvOEYBOvCAlm+S+G4+cXBik28WK+kLpLIiTUAhbDrvG61bez/fV23xCEBLMi\netTubTFPGFJnSYSEGiF5blywfSSEkHy3TkhoiCek95U8aUgjbB8JkSyGSZ5HIiQYRkj+O0BI\nqHHTzn8PCAm1aEIa4c4+IUGGkPz3gZBQC/jMhmk/+5uQYJnvYpj8s78thCR+qzsSEvC5dpN+\n9vcYryzwmJd+HBgS4rsYpv6kVUKCaf7vkO37i2wTLUZCkn6IPxLCNZL/fhASagH3kSb97G8r\nISl/YyAS4r0Ypv3s71HeWEdIkAl4HmnKz/62E1JjVwgJtUhe2UBIsI2QBqv3hZBQI6TBCAnv\nCGm4x84QEmpxhDTOBzQSEmQIycNCs30khJA8EBJeEZKPhWT7SAgheVkoto+ERBHSSL+fiJAg\nQ0h+FoLtIyGE5IeQ0EJInhaEhAZC8kRIaIohpJ6O5g3puleEhBoh+SIkNBCStwUhoUZI3ggJ\nT4Tkr2+/kCFCCpinJDxEEFLvciUkmEFIIfOUhDtCCpknJNwRUtA8JeGGkILmCQk39kPqX6sG\nQqIk3BBS4DwloURIgfOEhBIhhc5TEi6EFD5PSLhEENKHdWojJErChZAU85QEQhLMExIISTFP\nSSAkwTwhgZAU85SUPeshfVqihAQzCEkyT0m5IyTJPCHljpA085SUOULSzBNS5oyH9HF9WgqJ\nkjJHSKp5SspaczEst6exNzEUISEOzcXgnBujpUxCoqSsNRfD+e/fGC0REtL3upgO26W6pYD1\n+nltGguJknLWsZiOxfV6aTfqJn5FSIjE+2Lar1xpNeImfhZXSJSUsZfFcN5er46W+/O1pvVI\nmxgispAoKV+txXAoH2zYHG/fkD3DREhIX+t5pOuV0e78+EYxxiYGii0kSspW63mk9X7sTQzz\nZVkSEsxoPY80/iaGiS8kSspV+5UN9xOF7Gbd6yaGiTAkSspUV0gn3QMNr5sYhpAQi8di2Lum\n5RibGC7GkCgpT/ViWDY7OoyyiaG+LUlCghmd95HG28QgcYZESVmy/MY+QkI0HouhvDZq3Lgb\nYxODRRoSJeWIkPTzhJQhbtqNME9J+TEc0tflSEgwo7UYdsvL5bQUP/qdYUiUlJ/mYtiX943K\nt8c6E88jRRwSJWWnuRhW7u9ydMvLn/DtsRdCQg5en5A9uo36mdkcQ6Kk3LyGtHZ7QhLME1Jm\n2jftjvvyjbHctBPMU1JeXh5scG5bXiFJ3ynruV6/r0RCghnth7+L8h7SZfk33iZ+F3lIlJQX\nu0/IEhIiQkijzVNSTghptHlCyklrMWyXhl79HX1IlJST5mLYmnobRfwhUVJGmouhUP4Oiu5N\nDEBIiIjZz2z4YRGaD4mS8tFcDGs3ymetZhwSJWWjuRhOxUr7TqT3TfyOkBCT9k07Qw82pBES\nJeWCkMadJ6RMmH1CNpGQKCkThDT2PCVlob0Y9uvqzX2nETfxo1+WHyHBjNZiWN3uHrlCWlLm\nIVFSFpqLYedW5zKknfs31iZ+llBIlJSD9kuEzrdXNxh41I6QEJXXlwgR0gjzlJS+5mJY3q+R\njgZ+Yx8hISod95H24leBExIlpa+1GNb31zVIP43La73+tPLiCYmSkvf+PJJbaz9EiJAuhJQ+\no69sSC0kSkodIU0zT0iJayyG/b/ys09WG/V7ktINaYhFx9cCtw9D6gvztKov39X8r7WLI6RB\nP/3+TyKkhDwuzHPhlvvyneanv2X5QfojbGIIQkJcHhfmpvGY96r8JH39JoZIMKT3fxMhJeRx\nYS7d8/bcafZf6/LbXfPIQnr7VxFSQh4XZuue7+yvtSMkRIaQJpx/+WcRUkIIacJ5QkoXIU05\n3/53EVJCniGN9lQhIT21/mGElBBCmnSekFJl87V2yYbU+pcRUkJMhvTjKzxjDKn5byOkhBDS\nxPOElCZCmnr++Y8jpIQQ0tTzhJQkQpp8vv7XEVJCCGn6+cc/j5ASQkjTzxNSgghphvmFZvsw\nxGJIv35QCCHBDEKaY34h2T4MIaRZ5heK7cMQQpplnpBSQ0jzzC8E24chhDTPPCElhpBmml8Q\nUlIIaa75BSGlhJDmmiekpBgM6edf3BB3SJcFISWEkGabJ6SUENJ88/zOpIQQ0ozzlJQOQppx\nnpDSQUhzzlNSMghp1nlKSgUhzTpPSKkgpHnnKSkRc4b0+6//7jHWjk05T0lpmDWkzq/+vrJM\nhBA8T0lJIKS55wkpCYQ0+zwlpYCQ5p+npAQQ0vzzhJQAQjIwT0nxIyQL85QUPXMhDVhTdkII\nnSek6BGSiXlKih0h2ZinpMgRko15QoocIRmZp6S4EZKVeUqKGiGZmaekmBGSmXlCihkh2Zmn\npIhZC2nIYrIWQvA8JcWLkCzNU1K0CGnG+ZC32Qveag+hwAtjt3RuvffcRPYhvX9p0FUSIRni\ne2Hc/nO4uv2XceO3CUJ6N+W/H0JBIW3c5ny5nDZu57WJuReSyfkJH7WEUFBIhTuXp89u6bUJ\nQupASHEKCulxd/fz3V5CGjQ/3RsbIRQU0r9HSIXXJgip02SfNAsh/5DW293e/V1PnjefH20Y\nEtKkj1oZnSekGPmHVD+V4Vxx9toEIXWb6pdRQ8j7wjged7v1unrIYfOxI0IaPP/jUSAkQ4y9\nsoGQSoQUH0KyOP/bYSAkQ/wvjMN2Xd1LWm8OnpsgpF4/HQdCMsT3wjgvGy+eXPltgpB6EVJs\nfC+MjSv+jtWp077g4W/5/C8HgpAM8b0wCnesTx9lT8gOez+O5RCC5384FIRkSNirv7v+8vsm\nCOkDQooL10hW578fC0IyJOA+0v5UnVLeRyKkhq8Hg5AM8b4wVo1H7ZZvL2346R3RhPQRIcUk\n4HmkTfU8UrHe6p5HIqSmb0eDkAyx9coGQmr5cjgIyRBCsjz/+XgQkiGEZHr+4wEhJEMkF4bq\neaSBn48YQQih84QUC0KyPf/pkBCSIaZu2hHSuw/HhJAMISTr8/0HhZAMISTz871HhZAMCXhC\nVv/GPkLq1HdYCMkQ3wtjlDf2EVInQopAwItW9W/sI6RuPceFkAyx9DaKob9mK5oQgue7jwwh\nGeJ7YYzxxj5C6tV5aAjJEK6RopgnJOsC7iPJ39hHSP26jg0hGeJ9YXx+Y99PmyCkAToODiEZ\nEvA8kvyNfYT0yfvRISRDLL2ygZA+ejs8hGQIIcUz/3p8CMkQQyEN7Si6EELnCckwQopo/uUI\nEZIhhBTTfPsQEZIhhBTVfOsYEZIhhBTXfPMgEZIhhBTZfOMoEZIhhBTb/PMwEZIhhBTdfH2c\nCMkQQopunpAsshPS4I5iDSF4/nGkCMkQQopwfqHZPoQIKcb5hWT7ECKkKOcXiu1DiJCinCck\nawgpzvmFYPsQIqRI5xeEZAohxTq/ICRLzIQ0vKPIQyCkpBBStPMLQjKEkOKd9zhkGAshRTxP\nSXYQUszzlGQGIUU9T0lWEFLc85RkBCFFPk9JNhBS7POUZAIhRT9PSRZYCclnNRhZyLPPU5IB\nhJTAPCXNj5BSmKek2RFSEvOUNDdCSmOekmZGSInMU9K8CCmVeUqaFSElM09JcyKkdOYpaUaE\nlNA8Jc2HkFKap6TZEFJS85Q0FyMheS0Agwt59nlKmgkhJTZPSfMgpNTmKWkWhJTcPCXNgZDS\nm6ekGRBSgvOUND1CSnGekiZHSEnOU9LUCCnNeUqaGCElOk9J0yKkVOcpaVKElOw8JU3JRkh+\nl7nxhTz7PCVNiJASnl+Q0mQIKel5SpoKIaU9T0kTIaTE57l5Nw1CSn6ekqZASOnPU9IECCmD\neUoaHyHlME9JoyOkLOYpaWyElMc8JY3MREiel3JMC3n2eUoaFyHlMs8TSqMipHzmKWlEhJTR\nPCWNh5Bymufm3WgIKa95ShoJIWU2T0njIKTc5ilpFISU3Tx3lMZASBnOU5IeIeU4T0lyhJTl\nPCWpWQjJ91KNeCHPPk9JYoSU6TwPOWgRUrbzlKRESPnOU5IQIeU7z807IULKd/5CSjqElO98\nhZQ0CCnf+TtKUiCkfOcfKEmAkPKdr1FSOELKd/6JO0rBCCnf+SZKCmQgJO/LcO6FGPt8CyWF\nIaR859u4eReEkPKdf0VKAQgp3/l3pOSNkPKd70JJnggp3/lOXCn5IaR853tQkg9Cyne+D1dK\nHggp4/lei/5vNQVuPymExHyHn66UCKmBkJjv9MOlQkgN84fkf4Pc9EKMfv77lRIhNRAS832+\npURIDYTEfL/PKRFSAyEx/8mnlAipgZCY/6w/JUJqICTmv+lLiZAaCIn577ovI0JqICTmf9B5\npURIDYTE/E86UiKkBkJi/kdvKRFSAyEx/7OXlAipYfaQAl6yH91CjH++lRIhNRAS84M0UiKk\nBkJifqA6JUJqICTmB7unREgNhMS8hyolQmogpIjnQ4Vs/JoSITUQEvOeFnxISgMhMe8/T0o1\nQmI+ZJ6U7giJ+bB5UqoQEvOh86R0ISTmFfOkNHtIIZeAnYXEfPYP4RES86L5RdYxERLzwvl8\nUyIk5qXzuaZESMyL5/NMiZCYl8/neGeJkJgfYz67lgiJ+ZHm82qJkJgfbz6jR8QJiflx5zNJ\niZCYH3s+i5RmDinoEMeykJjPICVCYt57foDrnaW0fys6ITE/0fzbIw+EJNoEIWU332qJkESb\nIKQc5xdpfsAkITE/+fz9Rh4hiTZBSBnPLxJ7spaQmJ9vPqGYCIn5eecTaYmQmJ99PoWWCIl5\nC/PR32WaN6SwQ2dpITAvmI+5JkJi3tZ8pC0REvPm5mO8YiIk5k3Ox3Yzj5CYtzsfUUyExLzt\n+UhiIiTm7c8vFuZv6RES89HMW86JkJiPbN5mToTEfJTz91t7s/1W91f+Z3bYrqu9WW8Onpsg\nJOZD50OunkyEdF42yl75bYKQmNfM+z0YYSKkjSv+jtWp075wG69NuMAbunYuSOZNzA/MyURI\nhTvWp4+u8NoEITE/wvzPOZkIqXVP7fPdNkJifvr5H27tmQiJayTmY5hfNIm3Lzmz632k/ak6\nxX0k5iOZf8nJREiXVeNRu+XZaxOExPwM873XUHPszNVhUz2PVKy33s8jERLzEc+Pd2YDN0FI\nzMc8P96ZDdwEITEf87zozAQvESIk5iOel5yZ5CVChMR8xPOSM1O8RCj0UZO5DyTzec9Lzkzx\nhCwhMR/zvOTMFC8RIiTmY56XnBnXSMznPi85M8VLhAiJ+ZjnNWf2+SVCP72jN+x9wkAg37Xf\nZYKXCAHpm+CVDUD6CAkQICRAQBKS9m4bEB9CAgRIABAgJECAkACBCd7YB6Rvgjf2Aemb4I19\nQPomeBsFkL4J3tgHpI9rJEBggjf2Aekb6Y19QF54Yx8gwMMEgAAhAQKEBAgQEiAwZ0iTfewS\n0EW6mJVnFtG2f8H+hclq/wipH/sXJqv9I6R+7F+YrPaPkPqxf2Gy2j9C6sf+hclq/wipH/sX\nJqv9I6R+7F+YrPaPkPqxf2Gy2j9C6sf+hclq/wipH/sXJqv9I6R+7F+YrPaPkPqxf2Gy2j/r\n/1ggCoQECBASIEBIgAAhAQKEBAgQEiBASIAAIQEChAQIEBIgQEiAACEBAoQECBASIEBIgMBs\nIW0KV2zs/srMET5mXWj32DGjR/GxfzaP4m5ZHzTd8ZvrH3n7FbTLmbb+1dHmErg7PnbM6FF8\n7J/No7ip9qko8xEev5n+kQdXHC/Hwln99bNHt557F/pdj9vtYjN6FOv9M3kUj+7fubzO/Kc9\nfjOFtHH76///ue08m/9qZ3bPyn1b3ReqzaP43D+TR3F927dyF5XHb6aQ1u50MfpfrMrO7ebe\nhV5uc7kvVJtH8bl/po+i0x6/mUJyrvmHPWu3/3e9Gzr3bnQ6vh4+Y0fxuX+Gj+LZrbTHj5A6\nrW/3kldz70cP0yFdGiGZPYq78lYdIY3Oub/rf7U2Vm+aRBKS3aN4Ksqbc4Q0kbO9R5ZvIgnp\nxuBRPBfVtWQCIRVWl0Cb1f2775fZo9jeI3v7t7qlrTx+sz5qdzL2eNMbe0vgpvWoncGjaDuk\n03J1qk4oj99M/8Zt9Qj+3pl8ROdS/reqfOLb4BK9uS9Ns0exvsa0eBT39aMfyuPHKxs6bcqD\ne749YWeQ7Vc21Ptn8iieno8iJvDKhsvS7AOjlXNR7Z+5/9TfPW4sWT2K9/0zeRT/uecrAIXH\nb66QztXrbmfa+A/K/Vvae9j27hGS1aPY3D9rR9E1QhIeP2P3A4E4ERIgQEiAACEBAoQECBAS\nIEBIgAAhAQKEBAgQEiBASIAAIQEChAQIEBIgQEiAACEBAoQECBASIEBIgAAhAQKEBAgQEiBA\nSIAAIQEChAQIEBIgQEiAACEBAoQECBASIEBIgAAhAQKEBAgQUhyqX8Pa/evBTf2G1mwRUhSW\n1eXUGdKSi9ACLoUouP6Quq+mMDEuhSgQknVcCjG4/xbu6/82rthWX9otXbF7fuuyXzuTv+A8\nF4QUgzqkdXmi7Kc64Vb1t7a3X3lPSXMhpCg8btqtzpedW16vf8pT55XbP7/1d7n8cTNvNhz5\nKDxqOdxPr935eurs1u37SIQ0G458FJoPNtxu5N094znttytCmg1HPgrfQ1o9/o5ZcOSj8B7S\ny7f+ueVufyKk2XDko/Aa0trtO75FSPPhyEfBudOl2cyfK46Xy+72YMPtW4fLkftI8+HIR2Hp\nXNG68rndJSpOj29t7veZDjPvaLYIKQqH5UtI5Ssb3L9T/a3rnSS3OuzLqyjMgZAAAUICBAgJ\nECAkQICQAAFCAgQICRAgJECAkAABQgIECAkQICRAgJAAAUICBAgJECAkQICQAAFCAgQICRAg\nJECAkAABQgIECAkQICRAgJAAAUICBAgJECAkQICQAAFCAgQICRD4D+S4Dh+iR49mAAAAAElF\nTkSuQmCC",
      "text/plain": [
       "Plot with title \"Histogram of theta\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "hist(theta, freq=FALSE)\n",
    "curve(dgamma(x=x, shape=a, rate=b), col=\"blue\", add=TRUE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To calculate the mean of this distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false,
    "hidePrompt": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "5.51406757147033"
      ],
      "text/latex": [
       "5.51406757147033"
      ],
      "text/markdown": [
       "5.51406757147033"
      ],
      "text/plain": [
       "[1] 5.514068"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sum(theta) / m # sample mean"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidePrompt": true
   },
   "source": [
    "Now, theoretically we know the mean of a Gamma Distribution with shape parameters alpha and beta are;\n",
    "\n",
    "\n",
    "$$ \\mathbf {E} [X]={\\frac {\\alpha }{\\beta }}$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false,
    "hidePrompt": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "6"
      ],
      "text/latex": [
       "6"
      ],
      "text/markdown": [
       "6"
      ],
      "text/plain": [
       "[1] 6"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "a / b "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We do a second approximation with larger number of samples:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false,
    "hidePrompt": true,
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "6.00148518734382"
      ],
      "text/latex": [
       "6.00148518734382"
      ],
      "text/markdown": [
       "6.00148518734382"
      ],
      "text/plain": [
       "[1] 6.001485"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "m = 1e6\n",
    "theta = rgamma(n=m, shape=a, rate=b)\n",
    "sum(theta)/m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part 2:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Say we want to calculate the integral:\n",
    "\n",
    "$$I = \\int_0^1 x e^x dx$$\n",
    "\n",
    "We can calculate the closed form solution of this integral using integration by parts:\n",
    "\n",
    "$$u = x, dv = e^x$$\n",
    "\n",
    "$$du = dx, v = e^x$$\n",
    "\n",
    "and\n",
    "\n",
    "$$I = uv - \\int v du$$\n",
    "\n",
    "$$= xe^x - \\int e^x dx$$\n",
    "\n",
    "$$= \\left. xe^x - e^x \\right|_0^1$$\n",
    "\n",
    "$$= \\left. e^x(x-1)\\right|_0^1$$\n",
    "\n",
    "$$= 0 - (-1) = 1$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false,
    "hidePrompt": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "1.03695893700551"
      ],
      "text/latex": [
       "1.03695893700551"
      ],
      "text/markdown": [
       "1.03695893700551"
      ],
      "text/plain": [
       "[1] 1.036959"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# A seed value for reproducibility\n",
    "set.seed(12345)\n",
    " \n",
    "# The first approximation using n = 100 SAMPLES\n",
    "n1 = 100;\n",
    "x = runif(n1, min = 0, max = 1)\n",
    "sum(x *exp(x))/n1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false,
    "hidePrompt": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "0.99638471037357"
      ],
      "text/latex": [
       "0.99638471037357"
      ],
      "text/markdown": [
       "0.99638471037357"
      ],
      "text/plain": [
       "[1] 0.9963847"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    " # A second approximation using n = 100 SAMPLES\n",
    "n2 = 5000;\n",
    "x = runif(n2, min = 0, max = 1)\n",
    "sum(x *exp(x))/n2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part 3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Approximating the expected value of the Beta distribution:\n",
    "\n",
    "$$\\mathbb E[x] = \\int_{p(x)} p(x)x dx,$$\n",
    "\n",
    "where $$ x \\sim p(x) = Beta(\\alpha_1,\\alpha_2).$$\n",
    "\n",
    "$$\\mathbb E[x]_{Beta(\\alpha_1,\\alpha_2)} \\approx \\frac{1}{N} \\sum_i x_i$$\n",
    "\n",
    "\n",
    "$$\\mathbb E_{Beta(2,10)}[x] = \\frac{\\alpha_1}{\\alpha_1 + \\alpha_2} = \\frac{2}{12} = 0.167$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Left as an exercise.."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.3.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}