Mertapolis Hastings.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the given problem."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"book1_img2.png\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On solving above question we recieve this posterior distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"book1_img1.jpg\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our posterior distribution is proportional to g distribution. We're simply calculating g distribution but for numerical stability we are using log scale."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "log_gdist = function(mu, n, ybar)\n",
    "{\n",
    "    mu2 = mu^2\n",
    "    return (n*(ybar*mu - mu2/2.0) - log(1.0+mu2))\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have to write the naive implementation of metrapolis hasting algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "metrapolis_hastings = function(n, ybar, n_iter, mu_init, cand_sd) \n",
    "{\n",
    "    ## step 1, initialize\n",
    "    mu_out = numeric(n_iter)\n",
    "    accpt = 0\n",
    "    mu_now = mu_init\n",
    "    lg_now = log_gdist(mu=mu_now, n=n, ybar=ybar)\n",
    "\n",
    "    ## step 2, iterate\n",
    "    for (i in 1:n_iter) \n",
    "    {\n",
    "        ## step 2a\n",
    "        mu_cand = rnorm(n=1, mean=mu_now, sd=cand_sd) # draw a candidate\n",
    "\n",
    "        ## step 2b\n",
    "        lg_cand = log_gdist(mu=mu_cand, n=n, ybar=ybar) # evaluate log of g with the candidate\n",
    "        lalpha = lg_cand - lg_now # log of acceptance ratio\n",
    "        alpha = exp(lalpha)\n",
    "\n",
    "        ## step 2c\n",
    "        u = runif(1) # draw a uniform variable which will be less than alpha with probability min(1, alpha)\n",
    "        if (u < alpha) \n",
    "        {\n",
    "          # then accept the candidate\n",
    "          mu_now = mu_cand\n",
    "          accpt = accpt + 1 # to keep track of acceptance\n",
    "          lg_now = lg_cand\n",
    "        }\n",
    "\n",
    "        ## collect results\n",
    "        mu_out[i] = mu_now # save this iteration's value of mu\n",
    "    }\n",
    "\n",
    "    ## return a list of output\n",
    "    list(mu=mu_out, accept_rate=accpt/n_iter)\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now test it with some case study"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = c(1.2, 1.4, -0.5, 0.3, 0.9, 2.3, 1.0, 0.1, 1.3, 1.9)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "ybar = mean(y)\n",
    "n = length(y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAJYCAMAAACJuGjuAAAAMFBMVEUAAABNTU1oaGh8fHyM\njIyampqnp6eysrK9vb3Hx8fQ0NDZ2dnh4eHp6enw8PD////QFLu4AAAACXBIWXMAABJ0AAAS\ndAHeZh94AAAYiklEQVR4nO2dibaiOhAAAyIuo/j/fzsCLuiNiunuLFh13hnRS6cTqAchInEn\nAANc6grAMkEsMAGxwATEAhMQC0xALDABscAExAITEAtMQCwwAbHABMQCExALTEAsMAGxwATE\nAhMQC0xALDABscAExAITEAtMQCwwAbHABMQCExALTEAsMAGxwATEAhMQC0xALDABscAExAIT\nEAtMQCwwAbHABMQCExALTEAsMAGxwATEAhMQC0z4PbGcc49L9w+mrKNUZlM5FydTbBDLK9a/\nKsqG2ZxzI9ZCmCWW/yimTu3cIUaeBCDWh5Ui1WVxLLZhL3l1xOo2q/NSsxs/cdfV9uv+bLW/\nhBzP71bbSeSxdu15adecl+v2eC1vW7v63+m0rdzq32P6h/IexHKuuixVC/Ct/BZ8ywuxjtVF\np9WDWKvLcjNE/Lusco+sh4DrWu7f+On4/tjePrsxLW/q75mzcaNvezfIWjaIdXk579bzwao7\n7/jtZI83V2NGs6rb22uk68O2Z5W606l98OV8AJo6OfJQ3pNY/66rNk8yFslPivVHj/Hf/kTW\nnQ9Ct8/6Y4fbduezpBuOJruzLP1LdY/sher74MeHks6fbvvD2WF4ued+Ku+pj3Up5jhWoXAQ\n67J3e1luXanrHl/3x6+edhgVaC4nq909cv9U9Pjvv4eX+wpP5T2JtRlPge11paJBrMve3Ywf\nXNy6/6kb3h+HD26d6uc/n1fYtSt3E+v05+UWNy3vSaxu7L5X91IL5ifFely6vLTXntHxz5+u\nS+6vWOP7XT0x9b1YD0uPYo1HxP0yhkwR6/ZBtxsv2VYn7xGr8h6xhrf9qbFebw9fHbGq5z8O\nXbBVf934cH4tFMR62LvDKNP9s+ZjH2v4a335/KNYzds+1mDu/j6cVTSIdXmpLweT+6Gke3lV\n6J6kubx+PmK9vyq8no43Ri2PCmJdXvrT0PF2ZdZfI/avt5HPsd/zdxxrKGg1rLyvPor1XN6z\nWMfhTwvouiPW/eXaee+7WP24wLiwmnowHHHcw8j78PG/a79/GGF4J9ZTec9i9cfNhwHVckGs\n2wdD/2p1GUJqbi6tq8kA16H/rnD/R5r+42p9OF4H1D2lX3ko749YO7eMrvsviqVAZ9bB3i6k\n645YX+GG7xNPh5XV+epQLeH75wHE+oJ719vkfDWWvJA7/xDrC2631tgcVpYz1nBCrO/oNv19\nL9Xapn/d39u1Myk5AYgFJiAWmIBYYAJigQmIBSYgFpiAWGACYoEJiAUmIBaYgFhgAmKBCYgF\nJiAWmIBYYAJigQmIBSYgFpiAWGACYoEJiAUmIBaYgFhgAmKBCYgFJiAWmIBYYAJigQmIBSYg\nFpiAWGACYoEJiAUmIBaYgFhgQrBYXds/kHyzqOdmgh6hYvUPED511WSWEIAJoWKtXdOd/1kf\nz46tl/LQe9AjVKxhjqrLRFXdUqbpAD3CxTr1s6xN3gBMCD8VHvrJ/Yb5OTo6WfBMqFgHV7WH\nU1OdzdrXC5kJDRQJPontJ7ONLmX+F9BD0Dvarfv5QF2zOepVB5YC3W4wAbHABMQCE1TEYhwL\nnjESy03RSJEzLojUtTYmQvuWvgnDGrj0rYJYchDLA2LJQSwP4e37N0y87VzT/rNKUQiI5SG0\nfV096Ye+/w566ZsQsXyEtq911W64teF03Ffvb/Rb+iZELB+h7avGO2YGDu9v9Fv6JkQsH6Ib\n/bxv1FIUA2J54IglB7E8CPpY+/F2GfpYiOUhuH2ryVVh3ZmkKAXE8iAYx2qHcayq2TCOFSuo\nIBh5l4NYHhBLDmJ5QCw5iOUBseQglgfEkoNYHhBLDmJ5QCw5iOUBseQglgfEkoNYHhBLDmJ5\nQCw5iOUBseQglgfEkoNYHhBLDmJ5QCw5iOUBseQglgfEkoNYHhBLDmJ5QCw5iOUBseQglgfE\nkoNYHhBLDmJ5QCw5iOUBseQglgfEkoNYHhBLDmJ5QCw5iOUBseQglgfEkoNYHhBLDmJ5QCw5\niOUBseQglgdh+7a1c82HGaGXvgkRy0do+8YHJV+eF/n2EaSL34SI5UMkVuva7nQ6tm5rkaIY\nEMuDSKzKDU+17VxtkaIYEMuDSKzrxAFMIBApqCBEYq2vYjGBQJyggggXq9ls9253XuxaJhCI\nFFQQ4WLdZjZ2rmICgThBBRHcvsNhu22aoQvfvvVq8ZsQsXww8i4HsTwglhzE8oBYchDLg0r7\nGMeKFFQQRmK5KRopcgaxPHAqlINYHhBLDmJ5QCw5iOUhvH3/NsMMq65pmWE1VlBBhLavqye9\n85VJimJALA+h7WtdtTsMS8xij1geQttXucNt+cBtM5GCCkJ2z7vvjVqKYkAsDxyx5CCWB0Ef\na38cluhjIZaH4PatJleFNTf6xQkqCME4VjuMY1XNhnGsWEEFwci7HMTygFhyEMsDYslBLA+I\nJQexPCCWHMTygFhyEMsDYslBLA+IJQexPCCWHMTygFhyEMsDYslBLA+IJQexPCCWHMTygFhy\nEMsDYslBLA+IJQexPCCWHMTygFhyEMsDYslBLA+IJQexPCCWHMTygFhyEMsDYslBLA+IJQex\nPCCWHMTygFhyEMsDYslBLA+IJQexPCCWHMTygFhyEMsDYslBLA+IJQexPCCWHMTyIG/fx9m9\nlr4JEcsHYslBLA+h7XOPWKQoBsTyENq+fxViXUEsD8Ht6xq3Gh70zqkQsTwI2rdzbndCLMTy\nImnfceWaDrEQy4esfRtX7RELsTwI23eoP89Sv/RNiFg+xO1bIxZieeArHTmI5QGx5CCWB5X2\nMUAaKaggjMSaPSxviYtFUOW0W5sZSz4VhiSOFYNYRaTQS4xYSiBWmhjEesm/zTB1r2vaXKfu\nRayEhLavqyed15VJCjGIlZDQ9rWu2h2GpeO+cq1FCjGIlZDQ9lXucFs+uMoihRjESkj4rcmv\n3qilEINYCeGIlSYGsV5w7mPthzuT6WMhlo/g9q0mV4V1Z5JCCmIlRDCO1Q7jWFWzYRwrVlBB\nMPKeJgaxikihlxixlECsNDGIVUQKvcSIpQRipYlBrCJS6CVGLCUQK00MYhWRQi8xYimBWGli\nEKuIFHqJEUsJxEoTg1hFpNBLjFhKTNtXb47WKaKCWAl5vBHUWbiFWGpBBTFtX7dbW7iFWGpB\nBfHcvn+bWtstxFILKghP+w79k7a3pinigFgJ+du+/WrGj1BlKSKBWAl5al+3OR+u6n13tqsx\nShEPxErIQ/v+9Z33dvxdl95jrRBLLaggHsaxzger7fUHN+9/KxiaIiqIlZCHcaxmb50iKoiV\nkIdxLPsUUUGshDyOvF8WKrXT4HOKqCBWQnxiHXWfR4tYakEFcW3f/uExwLVFiuggVkJu7Zs+\noa/+8KP5wBSxQayEePtYdimiglgJ4Ua/NDE/I1Z/tBLOtfApRXQQKyGIlSbmZ8QqO4VeYsRS\nQtq+beXqD/duIZZaUEE8tG9bn07Het5ow6Fx1fa0YQIBxPIybd8wb3h/+6j7bNZhMKp16+50\nbN7fb4pYakEFMW3fyu3OwtSn3YzbR9f9k5Lb8d6a7v1IfYJNOFx9BF2FIJYSzwOkh16YGftj\nXOVym2leEwjcdHJP72cFhyQMiPk5sRq3/0Ks3XgOzGYCgceKuxefvy0hJGtAzE+JtXKHfe/I\nvFPh+nr7VrfOZAKBZ33ev31RRkjegJifEmu4w2HT76DPd5J21W03ug93McfahH/yPH8w57iF\nWEo8DjeMk5fUuzmR7VWn6u3xKqurws9qIZYSSxl59yVhuCEhCxHLmwOxErIIsV6c4b75dN5f\n9WJ+S6xNHXZ3Q+JxrFcJXnz+traIpcS0fZvQ22b+rv9wB72wih+z6/0BsdSYtq/SfMaMP0VU\n6GMl5Hnk3ThF1MIRKyHT9jXO5LfQlptQvXOHWEpM23esVl/87uvfZphh1TVtuhlWJVd4IdeM\nIXk0gwri8VQ4v7/dTX+HmOpGvw8lh/wZsZQIFat11W58kFa5s9inHK3/KbG+oXKH2/Ihm9tm\nvkzsWQGxlAhtn3Ov3qil+FQDhTXSxfyYWPtmuNlvxsO4Ux+xcroFBrE8PLRvNXavXPXZrHMf\naz+ulaSPNavMkJUQS4lp+7Zu1fVibd36c+Bq+nSat+NfBptwXpEhayGWEo9f6XTX37fMiPzX\nDuNYVbNJN46lkti9fKObRyGoIJ774PPFCkoRlbmJ3Ytl7TzyoIKYtq++HLEOeT/Rb3Z5dN4T\n4ulj7ZXvclDehPOLQ6yEPLSvmfUVjShFzNJCVkUsJf6OY7lm1o90QlNELCxkXcRSIkL7EEst\nqCAKE+urskJWRiwlJu3br/tbYVafbq+SpIjLd4ldQExIHklQQdzad7wPpa90J4UuRKxxdcRS\n4tq+rurnvzwvHHe13oxyDymiw3BDQq7taydjDKv+ySD6KaKDWAm5tq929/PfUXcgS2sTfl0O\nYiXkdpn9/KWhforoxSBWQooRK6CUkBDEUqIUsXI++iCWh1LEipU4pOmI5eEultlzPMoSK1Jf\nDrH0UkQvgVNhQkr4rjC0AMRKSAFiBccHSxLyXVBQogWDWL6YkK+vgxItl/zFCg8XnNbsbs8R\nBRVE/mLFTeyeXq3yIFYRKfQS03lXInOxRNVDrITkLZasdoiVkKzFStDtRywlchYrxTiFe7Gs\nnQexikihl9i9fKObB7GKSKGXGLGUQKx3MfNKQCwP2YqlUDGFjvisIhDLQ65iadRL4wpvThmI\n5SG4fd3audVl7miDpybnIpZVDGK9oJ9svH8yzViIulgqmx2xEhLavrZ/OFu3rYZfICaeCFMz\nMWIpEdq+agw8VvURsQJiEOtV3CWwW63UxdLa5lqSfCoHsTyEtq++zW1Yr5TFUtvkakefDwUh\nlofQ9t0nGTi6lapYeltc77T2viTE8hDcvvZm0/7Dz8W+S6G4wRX7S3mOp+RMePsOzXXpuKbz\nHieoIHIdeU+VGLGUyEws1eogVkLyEku3NrqShPwlKNEy0PnuRKmPpbyxlY8+L/+EWB6MxAp6\nwsi7FV0IIU0J+FuYWLEIqZwCeZ0KdQtRjwkZiwhKlDhGA8RKE5N35RTIRSyL7+MsYkK+SwxK\nlDJGg/C8/zbjLHTNpzlS5qT4uE42+8GzAmIp5u3qSf/w/WPhESsoUcIYDULztq7aHYal475y\nrTDF51Xy2Q9/10AsxbyVO9yWD+8n31la5/3v2EpAnpwaZEJo3i8e3700sVRi8q6cAhkcsVSu\nGkuLybtyCgj6WPtxWidpH2teDfLaD+7NO9VEaWI0CM67mlwV1t27Nd+nmFmBzPaDe/lGOVGS\nGA0E41jtMI5VNRvJONbc/LntB/diWT1RihgNchl5tyjENMZ5Fw0SJYjRALHSxORdOQWSivVF\n8pz3A2KlyfsyxTe5c9wP7unVLFHsGA0SivVV6iz3gwvPk2eDFEkn1neZ89wPLjhPpg3Sg867\nKMaF5sm1QWogVpqYvCunQCKxvk6b835ArDR5PSm+z5rvfgiZRTooEWJ9TBGQNOP9EPYLq4wb\npEIKsXLeppwKlaDzrhCTdeUCYjRArDQxeVdOAcRKE5N35RRArDQxeVdOAcTSigkZsw9KZB+j\nAWKpxXzXh8+5QRoglmLMN8E5N0gDxEoTk3flFECsNDF5V04BxEoTk3flFEAs7ZiZJeTcIA0Q\nSz1mXhE5N0gDxNKPmVVGzg3SALEMYuYUknODNECsNDF5V04BxEoTk3flFECsNDF5V04BxLKK\n+VBSzg3SALHMYvQfoIlY+ily3g8vY96plbxyxiBWmpi8K6cAYlnHvDpRqifSi9EAscxjgvr2\nSWM0QCz7GG9XK5fKWRGa94vZFn9eLO8a+VTOhtC8W8T6JubvKhlVzoTgvIfq/ZxfCimkheQc\nk3flFAjPe3g/H4VGCmEhecW4t28VE8ljNBDk3U6m0zFKISsks5jHDkNmlVOHq8KIMVO1squc\nMoiVJibvyimAWLFjXKxE4TEaIFb0mPGEmGnl1NDZ64xjfRfgMq6cEkZizR49nZ1iYTF5V04B\nToVpYoKeL4lY+ilyliTwiPW9WoilnyJnSYJPhd+q9Rti/dsMU/e6ppVM3TuXnCWhj6WYt6sn\nvfP3X0cj1tugL6J/QazWVbvxq8Ljvnr/dTRivQ+af0L8BbGqyTfQB1dZpBAXknPMQ9BctX5B\nLDd/wyBWykSlicURSxbzN2jGSPIviHXuY+2PwxJ9LLUj1ie1fkGs02pyVVh3JimkheQck3fl\nFBCMY7XDOFbVbBjH0gx6U9pviBU3Rc6SKB+xXn9vj1j6KXKWhFNhmryIFRTkWQGx9FPkLInR\nEevvKRGx9FPkLIntqXBiF2Lpp8hZEus+1k0txNJPkbMkdN7T5EUsjaCYP8DQALHSxAQEhf0s\nBbH0C8k5JjhRjNuZNUCsNDGCRN8duRBLv5CcY8SJLG8O1ACx0sQoJJp35EIs/UJyjlFL9Ekv\nxNIvJOcYzURv1UIs/UJyjlFP9OrIhVj6heQcY5FI7XnyGiBWmhizRM/P90Es/UJyjjFPdBUM\nsfQLyTkm78opgFhpYuJV7nrkimwYYqWJSVC5uGohVpqYVJVTe3bnl3mzTZGzJCWJdfvYWjDE\nShOTS+XM7EKsNDH5VO7at1c2AbHSxGRYuatg6Z6fniJFfvtBFpN35RRArDQxeVdOAcRKE5N3\n5RRArDQxeVdOAcRKE5N35RQIz8sEApKYvCunQGheJhCQxeRdOQVC8zKBgCwm78opEJqXx3HL\nYvKunAKheZlAQBaTd+UU4IiVJibvyikg6GMxgYAgJu/KKRCclwkERDF5V04BwTgWEwgIYvKu\nnAKFjLwH3ctxC/FG+4ucneYW7r/PZPLhfXG65txE82OEDVKmCLHOW0wwN7f3RyqvfrkyM8kt\n/NobeFn6ffEx5bxE82OEDVKnDLGu/wUldr5aeD/0fvKm6Ndi3f+9Lz6mnCnW7Bhhg9TR6QDZ\njmO503TPfBf3Yq+83FXf7O+7Vk9mTUq/Lz6lnJVofoywQfoYieXci20eVPxzgfAF0q0v2Gm5\n893ZY1a0rEiNI9Z3eRTXjEQ2FXnHq/6DIFpWpEIf68s8imvGIZ+avEH28IGvrgq/LfOFWJ+v\nCr/Mo7hmHMJrMvtGPw1kfQVvtLD7cQv392MmH94XQ1LOj0nXn/IRWpcvbvSDXyRUrC9u9INf\nJFSsL26bgV9E1CP2vgE4ccQCIwR9rLk3+sEvEnwSm3+jH/wignGsuTf6wS9CtxtMSCGW4lf3\nYILGTlYoI9ecS8tTVIMQq5w8RTUIscrJU1SDEKucPEU1CLHKyVNUgxCrnDxFNQixyslTVIMQ\nq5w8RTUIscrJU1SDEKucPEU1CLHKyVNUgxCrnDxFNYi7G8AExAITEAtMQCwwAbHABMQCExAL\nTEAsMAGxwATEAhMQC0xALDABscAExAITEAtMQCwwIZFYW/O8beWqNspzu+zbMqapozSoWzu3\nPnxe7xNpxDqYP7V0fC5cbZylx74tA+3QoMrcrGrIIzcriViHynpn/HPVoU9j/0w4+7aMady6\n6w+Oa+M8bZ+hdY24oBRibd3Keme0bn/+d+c2tmmitGWgGZOY56pcp5MmyY8pWvMN1Lj+ybsH\nhf/zPhChLQ/p4uRSeAx2CrEO9hvIRfofPEZbJnRxZpdp3VZcRqKrwsWIFSvJyHY4wxuzcxqP\nV0csrVwROFbmp/Yz26ZS6JsillYue7oq1jRra/m5MGrH8/5AXuudUS1RrFWMYbmBTt57X6hY\n41Xh0f6q8BRNrGO9OkZJ1KMzQ3wCrHfGZujl7qNM8hNHrH2cC8JxHOso/85ioWLFG3mPJNYx\n0jS2w8h71xTVx3pIa513nFk4ys6IItZacdKIt1RKG26pYnXD3Q3GSUYiXSFEEqu/LaSWj49y\nPxbYgFhgAmKBCYgFJiAWmIBYYAJigQmIBSYgFpiAWGACYoEJiAUmIBaYgFhgAmKBCYgFJiAW\nmIBYYAJigQmIBSYgFpiAWGACYoEJiAUmIBaYgFhgAmKBCYgFJiAWmIBYYAJigQmIBSYgFpiA\nWGACYoEJiAUmIBaYgFgSusuD9rsokwQXBWKJaMYpCiJM5VoaiCViP07SvHbxprkpBMSSUY9z\nKHMmfAaxZGz7k+A/zoR/QCwZw8x+G86Ef0AsIa3bn2rOhH9ALCEHtzpwJvwLYkmpXcWZ8C+I\nJWXvuCb0gFhSOseZ0ANiSTkfsTgT/gWxpKzk8ycvEcSSEWvi6eJALBmVa1JXIU8QC0xALDAB\nscAExAITEAtMQCwwAbHABMQCExALTEAsMAGxwATEAhMQC0xALDABscAExAITEAtMQCwwAbHA\nBMQCExALTEAsMAGxwATEAhMQC0xALDABscCE/5Vl/tSlElIkAAAAAElFTkSuQmCC",
      "text/plain": [
       "Plot with title \"Histogram of y\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "options(repr.plot.width = 5, repr.plot.height = 5)\n",
    "\n",
    "hist(y, freq=FALSE, xlim=c(-1.0, 3.0)) # histogram of the data\n",
    "curve(dt(x=x, df=1), lty=2, add=TRUE) # prior distribution for mu\n",
    "points(y, rep(0,n), pch=1) # individual data points\n",
    "points(ybar, 0, pch=19) # sample mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "List of 2\n",
      " $ mu         : num [1:1000] -0.113 1.507 1.507 1.507 1.507 ...\n",
      " $ accept_rate: num 0.122\n"
     ]
    }
   ],
   "source": [
    "set.seed(43) # set the random seed for reproducibility\n",
    "post = metrapolis_hastings(n=n, ybar=ybar, n_iter=1e3, mu_init=0.0, cand_sd=3.0)\n",
    "str(post)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message:\n",
      "\"package 'coda' was built under R version 3.4.4\""
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAJYCAMAAACJuGjuAAAAMFBMVEUAAABNTU1oaGh8fHyM\njIyampqnp6eysrK9vb3Hx8fQ0NDZ2dnh4eHp6enw8PD////QFLu4AAAACXBIWXMAABJ0AAAS\ndAHeZh94AAAVHklEQVR4nO2d7WKiMBQFg7rWWpX3f9st+FEgISQhB4LO/OgKDedeYZYipWBq\nAAFm7QbgPUEskIBYIAGxQAJigQTEAgmIBRIQCyQgFkhALJCAWCABsUACYoEExAIJiAUSEAsk\nIBZIQCyQgFggAbFAAmKBBMQCCYgFEhALJCAWSEAskIBYIAGxQAJigQTEAgmIBRIQCyQgFkhA\nLJCAWCABsUACYoEExAIJiAUSEAskIBZIQCyQgFggAbFAAmKBBMQCCYgFEhALJCAWSEAskIBY\nIAGxQAJigQTEAgmIBRIQCyQgFkhALJCAWCABsUACYoEExAIJiAUSEAskIBZIQCyQgFggAbFA\nAmKBBMQCCYgFEhALJCAWSEAskIBYIAGxQAJigQTEAgmIBRIQCyQgFkhALJCAWCABsUACYoEE\nxAIJiAUSEAskIBZIQCyQgFggAbFAAmKBBMQCCYgFEhALJCAWSFhALAObJ2Gr5xdphRKgBbFA\nAmKBBMQCCYgFEhALJCAWSEAskIBYIAGxQAJigQTEAgmIBRIQK5GUX99/EoiVyPpird+BD8RK\nZP3Nun4HPhArkfU36/od+ECsRNbfrOt34AOxEll/s67fgQ/ESmT9zbp+Bz4QK5H1N+v6HfhA\nrETW36zrd+ADsRJZf7Ou34GP5cU67Yw5nKUllmD9zbp+Bz4WFOu+Ivb3P2c8SkosyPqbdZUO\ngv8WdWmxjuZ4q+vr0ZwUJRbkU8WqA7fO0mJV5ta8vpmdosSCINb0wITsFNoV8Vwb/rWy+kab\nBrGmByZkp9CuiH9PsSpFiQVBrOmBCdkp/H4Y/Dqdzffvy9vRf/S++kabBrGmByZkp9C5v40x\n1U1RYkEQa3pgQnYSl8vpdDi0h/BHr1eIVWwHRYpVVIm5INb0wIRsMatvtGkQa3pgQrYYq0Ta\nvQeVrN8PYjlD4s5jmeL2Yog1PTAhezb2WvHeGRWxCumgeLEiSyBWIR0glhrEmh6YkC0GsQrt\nwLy+BI1MCNeCWIV2gFhyVhbLIFYKiDVZHrH+lgt/1gpiTZZHrBcnxMpYHrH+uFT71BKINSyP\nWB0uE3+cM14CsYblEavLyVzSSiDWsDxipbA1sZa/9mJVsUIKI1YifbEW725NsYLeLGIlglgh\nIxPCpSDWZHnESgGxJssjVgqINVkesVJArMnyiJUCYk2WR6wUEGuyPGKlgFiT5RErBcSyytnT\niJWhBGLZ0xNXHkk6Qiw5iBUyMiFcCmJZ5expxMpQArHsaePdyogVVAKx7GnEylACsexpxMpQ\nArHsacTKUAKx7GnEylACsexpxMpQArHsacTKUAKx7GnEylACsexpxMpQArHsacTKUKIv1vr3\npvpwsQKyESsRxAoZmRAuBbGscvY0YmUogVj2NGJlKIFY9jRiZSiBWPY0YmUo8S5iBT8E3lrQ\nnkasDCVSxAq/M2UKiWLViesPsbKQRyxv4lwQK2RkQrgUxHItOJhGrAwlEMuezitW2PC3Fyuk\nB8SaU8A3CrE8iXNBrBx5cxeZXQKx7OmsYkUdXSCWJ3EuiBU0NKqH+EVml3gzseJPrInFCuzH\nJdbIooiVyDyxohtSixXWj1Ms96KIlQhi+RdFrEQQy78oYiWCWP5FtylWUA+IFROPWDVi1YiV\nCGK5FhxMI1aGEohlTyNWhhKIZZdHrAwlEMsuX4BYxiDW7J6G6YiFWK7EuSBWjViuxLksK5al\nCWIlgVi+ao8ZiJWhBGLZM55iuZMQK6gEYtkzECtDCcSyZyBWhhKIZc9ArAwlEMuegVgZSiCW\nPWNCrKgKiPWYQizESmFNsQLvUINYNWK5EgMWQyx7WL/kR4n12tcgVkQ+Yt2nfGI9l0esiHzE\nuk+9g1gx9xpErDx8hlgRnSFWHhDLV+0xY1mxuntXxPIkBtTemFj9IdnFcnWDWB8hVu+ALUos\nM90GYn2uWHW/ecQKKIFY9gzEylACsewZiJWhBGLZMxArQwnEsmcgVoYSWcQKvGxhvIvesM2L\n1VkPiPWcShNrNN5d+93F6q0Q/+j+f8bixPr5OrQ7jMPxJ7LER4jl2ZcWLZZz2QXFuu0619Dt\n40p8hljjJRHLw9FU35f21fVcmWNUiTLEcq7cyYYcbSCWu6E0KnN5vb6YKqoEYtkzEGvYjjUx\nXeLNxbp/VJ0n1uA9f5BYOfZYzy3wbmI9MxArgd9jrPO1fZV+jPXcAoiFWH/sO58Kd7eoEvPF\n6n6WDxVr+Pk/RCz/A84zimUQ68XPsT2PVR2+Us9jzRDLFz9eO14sb5mce6waseaX2KhYzn0e\nYo02JOdtxHIlINZoQ3I2I9braAqxrPRosoiVeB6rOLGsbyFWp504RGIZ3+cpxHLMQKz5JRDL\nMQOx5pcwzz1ajVivGYg1v8Sj+08Uy31w0BWrMwKxIku8oVimP2tcLGfLHbHcHT6K+y4ejBHL\nKt7tArHG4sdrb10sx6L9Ec/FEAuxnuUR69mF94yCt0SkWI4tWr5YIxIg1hQnxCpUrL+wMdxi\nPXeyK4tVXyr/n1B4ShQv1t9G6g4fjPGKNSbBEmJ5Rg1zn8VfCQWIVV/8l/d5SiDWcC5idTl1\nrk6OKoFYw7mIlcLHiuVobUQs93cQK7IEYjnmItb8EojlmItY80uEi/W3thELsSZLIJZj7pRY\nr6tB3CBWjVgpYnVeOd8wYtWINUessdXVFytqpSDW32y5WPYPHmNv+XixTO+fYdB4h4gVVGIL\nYllVn2INLsRLEmv4DcRK4q3E6i6MWMFFJCCWVROxcpBDrMFmRSw7OZtY7j/Yi6UIsUYv7UIs\na7xcrOklgihBrPHRiGWNR6zxEoiFWDn4SLHcrSFWTrYrlj2kt/DjC2JlWmR2iQLFGhqBWJ4l\ngtimWN748W8ilqM6YnnHBNUe7PY695QdrNSSxbLPMr3aGLQ2BmItsMfanFgjbxix6ncUqzcf\nsRIXmV0CsSLF8rz7v5mI9dZi9X4pZSzcPU+K5XlPr5kjt9ZyjnVNIlbJYgUSKJb71hgjYk2P\n8Y9FrMLESlmDYWK5SyPWaAnEQqwczBfLfTXNWPz4NxHLM/YDxaqdm3U0fvybpvsCsfqTiLWQ\nWPXw856VtJJYzkGIVSMWYuWhJLHM3/ZDrP4kYiFWWGJM9RqxxsdM10Ysz1jEQqywxJjqNWKN\nj5munVOs12/lECvXIrNLrC1W94Yf6WI55weCWDlArIm2ECsNxJpoC7HSQKyJtl5ijV3Pjlhh\nJdYTyxQu1shoxAorsZpYdXfHgFj9ScRCrLDEmOr1O4n1eL1Jsf4mESvbIrNLmFFjhqM/S6yx\nKMQKLOEXy3kXtvXFGkTlF8vzIPElxDLvLtZIxtuL5Ru9iFg1YmURy3QHbF0s/5FqUHXEGiYi\nFmL9zUAsb5eDb2vFuk8j1jBRJ5anKGJNgVj3rw6xvFFvKdbo7QIQK6S2W6ze/qlwsRwn/xCr\nLlWs8XKuqDXFcqjhdCWwOmIFxk/VRizXNGJNxk/VfgexBr+bQKznDMTydjn4tuf7iNWfgVje\nLgffRqywEglixcRP1UYs1zRiTcZP1f48sQbXiyBWWvxU7Q8Uq/8SsdLip2rX3fPsiIVYo/Ej\ndxceq41YrsGIZaf9hQXVLk6sqQTECiyhE2v6nuhvKZbrl4cjU4gVEd93YrI2YjkzEMtKQyzE\nes5ALF+Tw9GIFVZCKpZPs/olltmSWJNR2cQaezOINSw3LdbjPAViIZaVNkcsdx+I5YsL7kvM\n0mL1/oLa1QxiOQYjlp1mRyDW2BRiRcQjFmI9Z5jx78bHIxZiPWcg1uhUQtQbinXaGXM4R5ZA\nLMQaX65dcH8/C3SMK4FYiDW+XLPg0RxvdX09mlNUiXcQyzt0EsQaX65ZsDK35vXN7KJKIJZe\nLPeNEbci1rMP/8pFrImseWKNXDvrqVy+WP+eYlVRJRAro1iRTfjEGs9ZVKzD1+lsvn9f3o7+\no3fEmshCrO5ynT9TM9UtqgRiIdY4l8vpdDi0h/BHr1eINZWFWGkUKJa9EGJNLRXbUvwis0sg\nFmLlALEmshBrNCSuK8RCrMAQe/u5z/u+aiLWdHJM3puKFVkCsXKLFVcZsQLjESuuMmIFxiNW\nXOXCxfr5OrRHUIfjT2QJxEKsUW67ztH5Pq4EYiHWKEdTfV/aV9dzVdwvoV+zEsVy13HPRqxc\ni7RU5vJ6fSnushnE8lUuWqze+izuBCli+SoXLRZ7rOAFJ7MQq8PvMdb52r7iGCsexBpn3/lU\nuCvtQj/E8lUuW6z659iex6oOXwnnsXI08ApArLjKhYuVXgKxECsHZYllHg88Qax+fcSy02LE\nGlsoWawkECsHiDWRhVhpIJad1bn+MeLPmHNUrhErOD5FrOHWXHTr3isO/l2uKmIFxieIZX1n\nebG6f/O7ZFXECo7fqFjrgFjh8YgVAWKFx2cRa24bWwGxwuMRKwLECo9HrAgQKzwesSJArPB4\nxIoAscLjESsCxAqPd4k1efYRsbqzEMsVn+M9INbEQgl11CBWSSBWeDxiRYBY4fGIFQFihccj\nVgSIFR6PWBEgVnB8lgsTEGtioYQ6asRi5aGQNvQg1rIU0oYexFqWQtrQg1jLUkgbehBrWQpp\nQw9iLUshbehBrGUppA09iLUshbShB7GWpZA29CDWshTShh7EWpZC2tCDWMtSSBt6EGtZCmlD\nD2ItSyFt6HH8hQliCSmkjVVALCGFtLEKiCWkkDZWAbGEFNLGKiCWkELaWAXEElJIG6uAWEIK\naWMVEEtIIW2sAmIJKaSNVUAsIYW0sQqIJaSQNlYBsYQU0sYqIJaQQtpYBcQSUkgbq4BYQgpp\nYxUQS0ghbawCYgkppI1VQCwhhbSxCoglpJA2VgGxhBTSxioglpBC2lgFxBJSSBurgFhCCmlj\nFRBLSCFtrAJiCSmkjVVALCGFtLEKiCWkkDZWAbGEFNLGKrynWIVQYk9LgVhCSuxpKRBLSIk9\nLQVigQTEAgmIBRIQCyQgFkhALJCAWCABsUBCYWL9fB2ae4abw/EnsgRilUVRYt125o99XAnE\nKouixDqa6vvSvrqeK3OMKoFYZVGUWJW5vF5fTBVVArHKoiixeg9ksZ7O4i+BWGVRlFjssd6H\nosT6PcY6X9tXHGNtnaLEqvedT4W7W1QJxCqLssSqf47teazq8MV5rG1TmFjpJRCrLBALJCAW\naDB1oWJxHmvbbEgs00VSE/JRrFiRJRCrMBALJCAWSChLLC70extKEosL/d6IksTiQr83oiSx\nuGzmjShJLC70eyNKEos91htRklhc6PdGOH8/0huQkJnaDBf6fRBc6AcSOPMOEhALJCAWSEAs\nkIBYIAGxQAJigQTEAgmIBRIQCyQgFkjYiFiwNRALJCAWSEAskIBYIAGxQAJigQTEAgmIBRIQ\nCyQgFkhALJCAWCChULFg8yRs9fwirVFSkLmRNkt964i1YORWMhFLmbmRNkt964i1YORWMhFL\nmbmRNkt964i1YORWMhFLmbmRNkt964i1YORWMhFLmbmRNkt964i1YORWMhFLmbmRNkt964i1\nYORWMjcqFnwCiAUSEAskIBZIQCyQgFggAbFAAmKBBMQCCYgFEhALJCAWSEAskIBYIAGxQAJi\ngYSlxTpWpjreMgSddq+gTubc+J/H+sgWeflnzL9r1sybOyg98/SUIGfwwmLt23uX7OYHHdug\n6tbPnBt/q8ywzXmR5/xtXqt75jVX5uV5Mxl3WmLwsmL9mOpSXyrzMzfoYv7dmv9q/3qZs+MP\n93WcL7L6Xfp2MMeMmf+atN//WLne+u9CdwncaanBy4p1NOffr9/ma27Q4d53s0o6mXPjvx93\ngsoW+d1KcDNVxkyT962fzP6R6E5LDV5WrINpduAXc8iU16ySTubM+OtzHWeL/Gcuz5fZMh8/\nrRtZc2T+qv8Qy52WGrysWJ3/bTm4mX0vc2b83lzvi2aL3Jn6q2p/aufL/Hr8KPzKk3kZLjtI\nSw3etFinZjedcYt915nFMubQHmjnzKxPzdF7dcqXiVgDrtWhzrfF2r19drGag/d/ufYud77a\nj2lf+fpErD63aj/InBW/a04KZBerOca6Nh/Ws2Wemh+Fv7KeEOuPKqdY+90wc078v/bTz33R\nTJG9rZItc2eaQ7ZbI2umzMdC7rTU4GXFun/EuOb4VHjd7a/DzDnx3VtPZ4rsnRXJlmmyZ/Y+\nFQ7TUoOXFeur3S2c28818zibvZ05J74rVqbI59LXptdsmfddSHtuLNtb7zQ7TEsNXlasbGfe\nry+vsp55f6zjbJG/R1e35njoO2Pm0TS/uDtmPJv/EGvLZ95/DxAa9tMDJ/jXeWRCJ3N2/GMd\nZ4v8cgbNy9znznweP7nTEoMXFuv+m/n5Od1ncXQyZ8c/1nG+yPPeETQz0xk0I/MpljstMXhh\nseBTQCyQgFggAbFAAmKBBMQCCYgFEhALJCAWSEAskIBYIAGxQAJigQTEAgmIBRIQCyQgFkhA\nLJCAWCABsUACYoEExAIJiAUSEAskIBZIQCyQgFggAbFAAmKBBMQCCYgFEhALJCAWSEAsD/db\n3Z2nB57/RsMdVoaHxy1OJ8ftuo8dgBZWhodWlQBfUMqGVeIBsdJhlXi4P0vm8XTM5iHUp/vc\n2655UMP5YO63E34MeQ3bPYZdD6Zqnx953huzDzhUeycQy0NXrMPrdufts+KOj1u4N/fx74q1\n7wyrHk/pOt1HntZ9MwuDWB46PwrPZn+rb/v7AxL39ydbfrcP+30Nab5+P57j8P0Ydro/S+nS\nfCfDI9Y3BGJ56FhzeDxzq32m4Y9zyP2xSfcnz+yfw+67sg/7MdiAWB561rwehfE6VL+ev/YD\nsawH/TVfj78/PC8XO/6tQSwPfrH2r4euTIjVPBf694jruvwbWBHE8jCwpju3eU7U7nS+Bon1\n+8PxuOMYC570jrHOvbmPf4ZiPY+xho+X7r36CD7r3UbysKb5GdZ+3Pv9kHfoiPVTX/adIdan\nwmfE7v75kT0WPHg+DrGqn0dUzYHSQ6zj46Dr5znEOo/1jPh+DfwgEMtD68bPrhWrOaVu/j33\nTA3/fgX6aX/q3Yc8zrxXrzPvr6/tmffP8gqxQANigQTEAgmIBRIQCyQgFkhALJCAWCABsUAC\nYoEExAIJiAUSEAskIBZIQCyQgFggAbFAAmKBBMQCCYgFEhALJCAWSEAskIBYIAGxQAJigQTE\nAgmIBRIQCyQgFkhALJDwH0KZk+CcJU3IAAAAAElFTkSuQmCC",
      "text/plain": [
       "plot without title"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "library(\"coda\")\n",
    "traceplot(as.mcmc(post$mu))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAJYCAMAAACJuGjuAAAAMFBMVEUAAABNTU1oaGh8fHyM\njIyampqnp6eysrK9vb3Hx8fQ0NDZ2dnh4eHp6enw8PD////QFLu4AAAACXBIWXMAABJ0AAAS\ndAHeZh94AAAYIUlEQVR4nO2diXqqMBBGx+Va26q8/9veyiabECB/SOSc734tV5OZAU5ZIoJl\nAAJs6wLgM0EskIBYIAGxQAJigQTEAgmIBRIQCyQgFkhALJCAWCABsUACYoEExAIJiAUSEAsk\nIBZIQCyQgFggAbFAAmKBBMQCCYgFEhALJCAWSEAskIBYIAGxQAJigQTEAgmIBRIQCyQgFkhA\nLJCAWCABsUACYoEExAIJiAUSEAskIBZIQCyQgFggAbFAAmKBBMQCCYgFEhALJCAWSEAskIBY\nIAGxQAJigQTEAgmIBRIQCyQgFkhALJCAWCABsUACYoEExAIJiAUSEAskIBZIQCyQgFggAbFA\nAmKBBMQCCYgFEhALJCAWSEAskIBYIAGxQAJigQTEAgmIBRIQCyQgFkhALJCAWCABsUACYoEE\nxAIJiAUSEAskIBZIQCyQgFggAbFAAmKBBMQCCYgFEhALJCAWSAgglkHyLFjr/kXaIAVoQSyQ\ngFggAbFAAmKBBMQCCYgFEhALJCAWSEAskIBYIAGxQAJigYTwYl2PZucfaYrZILJ3AopVXElx\nKi6quEhSLAWxvBNarItdHll2v9hVkWIpiOWd0GId7PGcfthRkWIpiOWd0GJVVxaOX2GIWMkT\nWqx/lVgHRYqlIJZ3gop1/rr+2Pff5OMyfvQeXCzM8k1Qseqr7M0OD0WKpSCWd0KOY91u1+v5\nnB/CX0a9Qqz02fPI+ysLYnkHsfIpxPINYuVTiOWbrcSKYRwLsYTEI9bK72cvqWGkGlgJu8Is\ndzpIwj2xR7Gs/dsQS8COxTLEEoJYiCUBsZ5i8Sm0dxALsSRs8CG0w4gCYiVPQLGu8YllxW/E\n8k/QqxsOJ3WKOdEru58/Ecs3QY+xbhNfzvGQYkZ0yxBLR9iD96vd1CncoyOWkj2fFSKWkF2K\nZWWSUix5wj2CWIglYddiZYglY4di1cNXiCUEsQqxMMsziIVYEhCr8QP8sWOxMsQSgliIJQGx\nEEsCYiGWBMRCLAmIhVgSECurP+ABj+xbrPrzQsTyzf7Esuang4ilArFiEWvzAvyCWIglAbEQ\nSwJiIZYExEIsCYgVk1ibF+GP3YtliCUBsQKkdAKxEEtCXtTmVXhjl2K9Dq0QSwViBUjpBGIh\nlgTEQiwJ+beGNq/CG4gVIKUTiPURYhliaUGsACmdQKyPEstiEStXHLHEBBMrQywRiBUgpQuI\n9aFibW0WYiGWBMRCLAmIhVgSEAuxJCAWYklArI8Q6zXKEJFYhlhqthBrY7MQC7EkFGJtvd30\nx97FyhBLw+7EstavuMT6pEe97kus15qLUKzPeoYwYsmTuoFYnyHWQCrE8snuxdIndQOxEEtC\ncbKKWFr2JpYh1rIusaQwGxnb3lQsQ6zExQqf1AHEWtglkhSjGyzE8gpiqZO6UIm19RmERxBL\nndQFxFrYJZIUiBUOxFIndWFzsfynRSx1Uhc2v8cuYq0LG61YhliI5Z38mmSrp7epwHvIIF0i\nSTF+UTlieQ0ZpEscKWx8CW4pVp0csaQgVuAK/M/8zsQaixyJWFuUgViro0YpVnMasZToxAqe\n1IHtxKq3k4iljLqVWVuLJfk+I2Kps06DWIu7xJECsYYzN08ePIdWd4kjBWINZ0YscdQYxAo7\nkFUkQyxxVMTyGTpIlzhSJCFW0CoQK0jUSMQKWAZiBYmKWD5DB+kSR4pIxeqkRSwhuxcrUCmI\nFSTqTsUqbgTuPXSQLiFTvO+LWL3UiOWjb1JiyYuxDLE89UWsbh7E8tIXsbp5PkWs36/z8xtY\ndr78qlJ8hFhBbvheilXl8p0voFiPo704SVJkY8NAkwNEMYmlH80qxLIPEOtih+9bPnX/OdhF\nkSL7ILHU1VRiWfJiHexWT9/soEiRjf2pI1Y3z4eI1VpU48ttjVhvO08GRSyf4YN0yQm0xUIs\nx7RWepW8WH/HWD/3fEp4jIVYzmmtUit1sbJT46zw+JCkGHu4Q0JiFetcnLbeFyYvVvZ7ycex\nDucv1TjWp4iVBRIre4nlOd2HjbzbyCUg00G3MWtELGVBDbHKtH7DB+kSLMWIWA4xIxNLOvxe\nBkcsx44fIFap1GtQXJT1M8XSjGPZx4iV6cWqsxUTHyuWNVkcNX2xmkc/ikEH6/z6MLE0KT5F\nrOJljVk9sd7W4SONuEuoFJ8mluIqh4VizS0EsTykXQdiLe8SKgViOaSzZtYZYs2rBLE8pF3H\nW7EyhViGWEv6fYZY1U/FeINIrDczMa+0+V2Kfm1EKd72dhcrsF+ItbxLzhWxxrJ2X0AsZ26H\n8a9QrE7xeuzfspgqsSZW2fALKrGyxWKNNdhUrOw2fnnf6hT2Ou5dFBOxRvrFLNbf3vA23Wh5\nik8TK9OLNTQx2C1uscQpPkqsclIu1vs6mu9NjXx8vFjvu8cq1sSbW4nV+R9ive2epFj+i3ET\nqyt0eX38WNipF1xKk4NYrm/PZrlYox8XfrpYI90Rq4yHWIt7rRNLdkHBgjcRyy+I5fr2bOaI\n9Xplt2K1RxkQayyeo1jWEmu4T7fBm2jOpclZJJbFKtaqz0VFYg3UNF8sQyy3mNXi8zyHsYn1\nZrfWEqvzXez9imWpijWVy79YncOnwUyIVXVobbgH+88Vy99srhPLt+XFP8Ry7IBYjiDWvA6p\nijXpjeKLpMO3slko1vBbiNVvEpdYPuqwzvQ7sRoiuYjVHOxCrOm0TbG8zOvIZa2BxOosmuGP\nkzvX/CFW3gGxpgpohEMs9/ZexKpk8i/W2yGMiMRqvOwiVjlT9ubaSsTqtfEvVqHW8HsBxLLu\novEs1uAGGbF6bV5iefoecvll+eH3thDrTVbEGmjfEWviE4vRvBKx3m2yphPELVZ9m7h+0gV1\nypkvVmfpRS2Wtd5zqsm18bv8nRhOYpmbWG/OTT5HrPZcbyjWwDL2JdayRT8g1nC7KbE6lZdL\nvXNGnvWbOtcpJ22xup2K6hBruk45sYi1aF57l/VmnsQauzOFU0G1Aojl2HxQrE6QaMQy67zp\nUFPZcXgdOhc00b0rVqPkXi0tsYaW08eK1ZvZoGJ13Vkt1msbsk6s8c8tZ4hliDUjpjW08C9W\nY7sxo6qOWAuOsxqZx6/gWSzWwiOP9V3EKaITq+eONV7sfsAyXVO5zhFrNevF6kuxrVjZarGs\n/DG/pnr/OXXzNMTqNkesiQhWFTbarjPyjliDYnVm1imkNQ6G/B5jVW5kC8UqV59YrHqfmSFW\n3tynWK/VoBGrP1w5Fa4qzPqxHQtCrGUp4hcr+1CxMsRyC/kaKlp4PJONimULxCrXXyCxrJFy\n92KZg1iOZ+kNsTJrLdsZ5QyLlbXFMufwsYplVrUciopY3VDxiWWbiJU5iVUVOJR0QZ1qthQr\nCyWWRStWz6O9ipV5FKtK70Os1p9/OmLZW6F2KVb3xfViVRuvmbyUHBKrXsufI9ZgRbsRyzGU\nF7Gyd2K9SrW5YlVd14k11RCx2q2HVlF0YmWfJVZzxoaTLqhTzWeI1T7Fypo7yZBiNYZ8EWtm\n60GxWi8tEKv6NZeOWP1jvWVi1XYpxcoQq9s6PrHMBgqrRkaaHx451LRKrKzu7yKWjYnV3PpO\nZ50JYk2XYzYiViP2LLHyacdOnXqqvft0SxsQ6/V/67k2mnUecYo18OpisTLzIVYmE2vh391C\nsawpliFWzGJZ1vjlUFNzzSLWCuIQa34ldRfE+lyxWq+7rgt/YtUCDIiVVWI5lRVWrAyxJlt/\niFjNDxE640zO9awS69WrcekVYhWT4cQqV8a4WCYXy9pTvsTKdiTW2+W8QiznLf6brK/rGMyr\nWC9LnMpod3E8pMvee9yYH8QabzMcy59YvSukG2L13nxblC+x3LJliFU0jkosQyyn9/10kaZ4\n23iBWPWh1UqxsqTFak93siPWGrF6v92rcRIrS0msdlDrhJ/M7wxijQZwESsLKNZLhPliDXUb\nOtcdy+9OkmItDeZdrHbb2esasVaRpljlViirfhVajYjluiFdIJZ1xepV4pAMsabecI/Xazmj\nlK5YQ6tz7H/vA3fX7ZwuVr+wRCzbr1jv/+rDimUxidX/pgNizU3hV6zlpSwRyzX0ZmJliDUa\nZSux6lfXVeJDrF4lbsn2LJbDO/GItawGL2LNP1PIhupHrPBi5edQsYq14BQ0Q6yJd7YTywuz\nxRr8mrI3scypBKdka7soU0QjliFWq+NsEGu4IWK1O84GsQabtW5SFKVYS5INi+Vw5IhY60vJ\nF3V9iyJbm3EowcyaVonV+V+vG2IhVqsxYiFWL8HMmprD7Ii1MIWDWGsOeBBrMPi+xbLJJl5K\n+WSxBv4oESuYWJlarLk1DYq1MPeAWC7fW0Os9bUEF+tt9HqGxWI5HFwg1vpa6kuQ4xGrM70m\nN2K9fWvdYCVi9V4xl73xglTzuyhTRCJW2dTT6pwowVksFwUmcg+LNd1vQar5XZQpEKv/ukex\nhoYbEKv8lbpYAxnHX0csDykQq/96LpZlXsTqJ9mHWGNzGVqszGHg0AuuYjkdZs9Ojlh7E8sQ\ny1+KsaZedgaI1csSp1jXo9n5x1uKabFWXnOXmlitA/ddiFXUcyo+arr4SjEp1tprOcd2tVZu\nEWMRqyVTOUDuv6IoxbrY5ZFl94tdPaUYFau6YHgNU2I1L9WNTixJRVGKdbDHc/phR08ppsVy\njzU3QdxiZbsSqypqvDjEml+T7Vusf5VYBz8pRmdSKVZ1/BaRWFaJZfUrOxHr/HX9se+/ycdl\n/Oh9hljjGYOIlcUglg2IZY4SzM0enVj1tYdmh4efFNuLlcUjVrZPsbLb7Xo9n/ND+MuoV8mI\nVe1rqiybifW6gHWPYilSTIi13qs4xeqmqfYFQcRymsmPFsvDKNbbVVOKZdVOsXwxArFaZ94S\nsVxArIUpyuMra70fajUi1vsg4/sw9zgTWfYjVjYoVqBNaJ94xLImnorZn1i2d7E8pZho6MOr\ndMTKumLVPzYAsebnsCp075IcxHpVFKSLLMWmYhlivSeoWL9f53x9nC+/nlJsJlb1IeFGYnWL\n2rVYj2Pj6PzkJ0UQsXqbhyxesbLGEMNexLrY4fuWT91/Dp4+hA4hVi/Je7HCrcae0Vbm754J\n7kCsg93q6Zuny2bCifUKVJ3FD40/ItaKvEtLbS0LLwOkk5+g+Nk1tcWyMbGC8UasrCfWViS9\nxZr+aE4pVralWNmYWDEQ9hjr555PLTzGGjg9Q6zyf3sWq/rmV85xwYV+nWXZuEBkZqCZ9MWy\n+iXEekPYcaxLPo51OH8tGsfaUqzXpqk94B6FWAOnp5sTVKyVKdoehRbrtUlALAcQy62ad2Jt\nBmItITqxWteiRLEau2M321fUJjGxWsOU24iVxSlWbKQrVv35WIhqXpuEfIwhhjXZKCGGcrqk\nLVaYBVqJ9RpjiGFFItYC4hXLEMsNxHKqJmqxYqimB2I5VfMS6315gbH6ZxTldEEsh2pe92cw\nxHIEsZzKQay5IJZTOVaNuCOWI0mLFWqRVmIZYjmDWE7ltMWKAsRawBuxDLFeINYCEGuacq8c\nxQeXA6Ql1uuN+gwtBIg1H8RyKgex5oJYTuW0x7GiALEWMJjCYhArVEoHSrFi2oo2SEmsDLGa\nIJafFNYVK9xajlksT7cR8A5iOZUTo1jlMkCslSlaYgVemIZYc0EsF2IWi7PClSmsufdDrAyx\nFjAoVvONLcSK7tt7iOUlRVssfQ2t3FGuP8TykgKxulRibV3HIIjlAmLNJk2xQi9LxJoNYrkQ\nq1jxnanWIJYLkYoV7lLH+SQmVutnOOIUK+A1tPNBLBcQazaI5UJ+mBzfOkQsHyk2FytwTgcQ\ny0cKxOqBWD5SIFYPxPKRolTK3r2vJNLVh1g+UiBWD8TykQKxekR68WgOYrkQ6fpDLB8pthQr\nWuJdEIgFEhALJKQmVoZYaYBYIAGxQAJigQTEAglJioVX8YNYIAGxQEJyYj1/I1b8IBZIQCyQ\ngFggAbFAQoJiRXw9LtQgFkhALJCAWCABsUACYoGEZMSqX4r1qUTQIj2xov42HVQkKFaMd6qC\nLogFEhALJCAWSEAskIBYIAGxQEKSYukLgLWkKBYkAGKBBMQCCYgFEhALJCAWSEAskIBYIAGx\nQAJigYSgYv1+nZ9XrNv58js7BWKlRUCxHkd7cZqbArHSIqBYFzt83/Kp+8/BLvNS4FViBBTr\nYLd6+maHeSkQKzECitW6jGr8mirESh62WCAh7DHWzz2f4hjr8wk53HBqnBUeH/NSIFZihB3H\nuuTjWIfz1+xxLMRKjFRG3hErMRALJCAWSNhKLMaxPpx4xLImmqQQjkR2hXiVGogFEhALJCRy\noR9ipUYiF/ohVmokcqEfYqVGKpfNQGKkcqEfJAZbLJCQyoV+kBipXOgHiZHKhX6QGImMvENq\nIBZIQCyQgFggAbFAAmKBBMQCCZGKBcmzYK37F2mLlIKYiZQZ66wjVsCQqcRELGXMRMqMddYR\nK2DIVGIiljJmImXGOuuIFTBkKjERSxkzkTJjnXXEChgylZiIpYyZSJmxzjpiBQyZSkzEUsZM\npMxYZx2xAoZMJWaiYsEeQCyQgFggAbFAAmKBBMQCCYgFEhALJCAWSEAskIBYIAGxQAJigQTE\nAgmIBRIQCySEFutysMNl9CbLjlyPdaBGzLXhf8vl4S3k7Z/Zv7vXmI/hQMtjXisJfAYOLFZx\nE+/j+kCXPNDh0Y65NvzjYN0y14X88V/m/VDEvPuKeatuJjMcbWHgsGL92uGW3Q42cQPvaW72\n7/H8U/vXirk6/LlYxv5CHv56P87PByx4i/kvf1zDxdus/3UqJBiOtjRwWLEu9vP389u+1gY6\nF3U/F0kj5trw3+WdoLyF/M4leDwfCeMtpvmd9audyojD0ZYGDivW2Z4b8JudPcV7LpJGzJXh\n79Uy9hby3+uJQ95ilnvrp6w+Yv6pX4o1HG1p4LBiNf7afPB4PoKzEXNl+JPdi67eQh4t+zrk\ne21/Mb/KXeGXn5i3bt9OtKWBkxbr+txMe1xj35lnscyKx8L4jJldn0fvh6u/mIjV4X44Z/7W\nWL619y7W8+D9n6+tS8FXfpr25a9OxGrzOJw6MVeFPz4HBbyL9TzGuj9P1r3FvD53hX+yXhHr\nxcGnWKdjN+aa8P/ys5+iq6eQrbXiLebRnodsj6esnmKWnYajLQ0cVqziFOPu46zwfjzduzHX\nhG/eetpTyNaoiLeY5j1m66ywG21p4LBifeWbhZ/xJ7I68WOnfsw14ZtieQpZ9b4/a/UWs9iE\n5GNj3ma9UWw32tLAYcXyNvJ+r73yOvJeLmNvIe/Ph8/+HQ99e4x5secHdxePo/mlWCmPvP8d\nIDw5TTec4F/jkQmNmKvDl8vYW8ivwUDrYp58x6yOn4ajLQwcWKzik/n1cZrP4mjEXB2+XMb+\nQv6cBgKtjDkYaEXMSqzhaAsDBxYL9gJigQTEAgmIBRIQCyQgFkhALJCAWCABsUACYoEExAIJ\niAUSEAskIBZIQCyQgFggAbFAAmKBBMQCCYgFEhALJCAWSEAskIBYIAGxQAJigQTEAgmIBRIQ\nCyQgFkhALJCAWCABsUACYo1Q3OruZ7rhz6s1FLAwRihvcTrZ7th87ADksDBGyFVx8AWl+rBI\nRkCs5bBIRiieJVM+HfP5EOpr8erj+HxQw8/ZitsJl03qZsey2f1sh/z5kT8ns5PDodongVgj\nNMU617c7z58Vdylv4f68j39TrFOj2aF8Ste1aHnddmYCg1gjNHaFP3Z6ZI9T8YDEU/Fky+/8\nYb91k+fP7/I5Dt9ls2vxLKXb8x0Pj1hPCMQaoWHNuXzmVv5Mw9/BJsVjk4onz5yqZsWmbGe7\nwSeINULLmvpRGPWh+v3n69QRq/egv+fPy9/O83brh/9oEGuEcbFO9UNXJsR6Phf674jrHn4G\nNgSxRuhY03z1+Zyo4/Xn7iTW387xcuQYCypax1g/rVfLX12xqmOs7uOlW1O7YF9zO5PSmuc+\nLD/d+zvJOzfE+s1up0aT3llhFeJYnD+yxYKS6nGIh6w6onoeKJViXcqDrt+qSW8cqwrxXTfc\nEYg1Qu7G7zEX6zmkbv+qLdOTf38C/eZ7vaJJOfJ+qEfe65/5yPu+vEIs0IBYIAGxQAJigQTE\nAgmIBRIQCyQgFkhALJCAWCABsUACYoEExAIJiAUSEAskIBZIQCyQgFggAbFAAmKBBMQCCYgF\nEhALJCAWSEAskIBYIAGxQAJigQTEAgmIBRIQCyT8B2I0GfBobR3vAAAAAElFTkSuQmCC",
      "text/plain": [
       "plot without title"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "post = metrapolis_hastings(n=n, ybar=ybar, n_iter=1e3, mu_init=0.0, cand_sd=0.05)\n",
    "\n",
    "traceplot(as.mcmc(post$mu))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAJYCAMAAACJuGjuAAAAMFBMVEUAAABNTU1oaGh8fHyM\njIyampqnp6eysrK9vb3Hx8fQ0NDZ2dnh4eHp6enw8PD////QFLu4AAAACXBIWXMAABJ0AAAS\ndAHeZh94AAAYzklEQVR4nO3diXrivLKG0WJoQhJA93+3HYMHyYOQbJWQ7fc7z/k3DVKVbK0G\nh06CGEIUIp9eANlmgEVUAiyiEmARlQCLqARYRCXAIioBFlEJsIhKgEVUAiyiEmARlQCLqARY\nRCXAIioBFlEJsIhKgEVUAiyiEmARlQCLqARYRCXAIioBFlEJsIhKgEVUAiyiEmARlQCLqARY\nRCXAIioBFlEJsIhKgEVUAiyiEmARlQCLqARYRCXAIioBFlEJsIhKgEVUAiyiEmARlQCLqARY\nRCXAIioBFlEJsIhKgEVUAiyiEmARlQCLqARYRCXAIioBFlEJsIhKgEVUAiyiEmARlQCLqARY\nRCXAIioBFlEJsIhKgEVUAiyiEmARlQCLqARYRCXAIioBFlEJsIhKgEVUAiyiEmARlQCLqARY\nRCXAIioBFlEJsIhKgEVUAiyiEmARlQCLqARYRCXAIioBFlEJsIhKgEVUAiyiEmARlQCLqARY\nRCXAIioBFlEJsIhKMsASsvrM2PX0kD7QgugGWEQlwCIqARZRCbCISoBFVAIsohJgEZUAi6gE\nWEQlwCIqARZRCbCISoC1nsz5joGPBVjrCbAWZ01nMF+AtThrOoP5AqzFWdMZzBdgLc6azmC+\nAGtxCjyDBexqAUsID7ACU8CuFrCE8AArMAXsagFLCA+wAlPArhawhPAAKzAF7GoBSwgPsAJT\nwK4WsITwACswBexqAUsID7ACU8CuFrCE8AArMAXsagFLCA+wAlPArhawhPAAKzAF7GoBSwgP\nsAJTwK4WsITwACswBexqAUsID7ACU8CuFrCE8AArMAXsagFLCA+wAlPArhawhPAAKzAF7GoB\nSwgPsAJTwK4WsITwACswBexqAUsID7ACU8CuFrCE8AArMAXsagFLCA+wAlPArhawhPAAKzAF\n7GoBSwgPsAJTwK4WsITwACswBexqAUsID7ACo7+rbxsAa3HynMGoj6gCVlx2DSumD7DiAqzQ\nwcCKCrBCBwMrKsAKHQysqAArdDCwogKs0MHAigqwQgcDKyrACh0MrKgAK3QwsKICrNDBwIoK\nsEIHAysqwAodDKyoACt0MLCiAqzQwcCKCrBCBwMrKsAKHQysqAArdDCwogKs0MHAigqwQgcD\nKyrACh0MrKgAK/3YmQHWso7Xo8j5R7VFYIClmIywXuflJM9cVFpELghYeskN6yKXhzH3i1w1\nWkQuCFh6yQ3rII/q9kOOGi0iF7QTWB/xmBtWc3L8JwlYowOANT6vmvivgXXQaBEXYCkmK6zz\n1/VHvv9uPi7+q3dgjQ5Qh5XwKi4rrFeeNw8PjRaRC9otrKlC64Rlbrfr9Xx+XsJfvK6ANT4A\nWEsDrNEBwFqaXLDCGwFrcVONKQW2MDuCNZi3YVhlvI8V2EgEWNFNs0wZKTI8WDspWrxfA7AS\ndRgrlWVKgS0MsNJ1GCuVZUqBLQyw0nUYK5VlSoEtDLDSdRgrlWVKgS0MsNJ1GCuVZUqBLczK\nYFVrANbSAGt0DUthtfO3BUsk+B2F3cESA6zZZ/gKLN9SzGD7B2sA1nhuh5N2i6iUA0uANW9K\nndubH85J0CImwJocuTx5L96vctNuERFgTY5cHr4qDBm4c1hzGgMrZCCw4ptmmVJgC5Mc1oJF\nA2vmlAJbGGB5Rgbd7W+aZUqBLQywPCOD7vY3zTKlwBYGWJ6RQXf7m2aZUmALAyzPyKC7/U2z\nTCmwhQGWZ2TQ3f6mWaYU2MLsD5YA6xOwPD0/Cqv+9/pyYIUtA1gBPT8LywArWYBllwVWsgDL\nLgusZNkHrG4CsGZOKbCFAdb0SO/dwAroYp0kYAErYRdgjY303g2skC7AsocBK1kXYFnDpg4S\nWPFdEsOK2ntgLZ5SYAsTBWvynPfLJYbVkX4La/pRYM1psWA1MbDeNmp+z3hk/2ayARawJopp\nwpIhrH4nYKVtsXVYr87lwXp3pMAKKJQY1lAQsGZOUWwBLLvGeD9gzWgBrKY9sJK2yADrtWnA\ncv4ELP9UDVjvlhQAq70NrKQJaiHhQ6cLlA7LhMASYKVsASzrJrDStWhO+aIuu4HlzAfWm0HA\nam4CK12LcmG9OenASjZFpcXU6Qxe3+5h9cEEwuqO1BNgvZ8GLGAFT24GAssaBqxmUHZY3kEG\nWKPLCw2wJgYZDVgGWAmTGdZ7MzuD1asFrPDJxvRNAQtY7aDR92YSwbJfJENhCbBGlhealcB6\nWwNYwJoaBKxmdjZYzShg+ZusGJbYt2bA6n9oZCCs5v+B5W0CLGuYBct5DFj2jfXCEvMGlrRz\nrAlJYTmrA5ZzY7uwBFhzA6zeGoGVJjNgDb/keXfswOp/mQmsZlCpsPwH0DwMrMVTVFrsB5a8\nEABrcZLB8h0+sIA1NSgG1qAmsIA1NUi6swOsFLCsle0dVnsuyocl9lC3A7DmTVFp4YclwHKr\nAyu0xTtYsitYnjfPgBXXYgSWtUGvW1qwhgvcICyngT11R7BkC7AsNz5Y0t5rzQ6D1bZ43Q2s\nqUHAaifUSsYaAyuuBbCs2cBK12L9sNydA1b8FI0WsjFYUiwsAdZSWNJ/uL6VBVbzeDAsRyKw\nZrUYOQxNWC6yQmEZZ4J0J2DQF1iTLaJgtbziYVknPh8s0YUlI7DsJzn3CKzj3zOsHqE5sOwX\nkd3AMsCq7wiCJaYHS0Jg7fEZC1jNHbNgOadk8Na0C8vZl3FYFpglsMQPy9pfYKUOsJw12rBE\ngDU/W4Rlv9QtgdWtHFjLW+jCMm2pDLAkMyxZBKs7odaKioP1+3WujkbOl9/IFstgDZB1TWxY\nkgeWSQmr2WMPLLN1WI+jdDnFtSgQ1qAasNxq0ZkL6yKH79vz1v3nIJeoFn5YzVmfhDU4PW2T\nImF1S/08rPaorBUVBusgt/b2TQ5RLd7AeveMNXIafbC6AvXAz8Iyn4UlpcNyluJf1zxYMg2r\nf3raJvqwmr7Aepdin7Ha8/xJWE0p94He/64NlnMyioP1d431c3/eSn6NVQgsaXZsu7C6NZQD\ny5ysrwqPj6gW/SdjEw9rhIJ5A6u/8WGwRANWB6Rb+QxYTsFZsNqHC4Jlfi/P97EO568l72Nl\nhNVtdRZYr+clDyyxjw1YCVpsAVZXbRKWaQpMwWoOvb6YG4ElI14WwRJg5YDV3hULSxLCso7A\nglVj6Z8+YPlaAAtYgUW8VYDlHCuwIoqM2LHiGd5daHwGln1141T7MCwzAkva6mlhdedmMp+C\nFdliDFZ7yjLDMruHZZ2byawcVncjPSwZwpJ4WM0ftGE1RwGsqBY+WK2JAmA1+y92l3XD6v5u\n2AvaDSwxZcKylxIBy0g3HVjpAixgqcQDS0ZhGT1YIgtgSfefxbBksJnAes1zE9XCgtWdc2AB\nq8p1x7CsQwZWwimv3A7+H6HwtMgPq/1jGljtWnLCkp3AMjf/t/d5WgBLA5Z7DGOwZB2w/l4N\nb+8HjbUYwhLrZjQsa0QkLJMIVgOivkcblmwc1uwW5cEabjKwemuLzSphidhV38Fq/s8Ly11k\nD5YNoxBY9fK7l3NgvYNlAmA5VdVhOTDKgNWchW7qJCxrHrDmwmq3fghL2nIWrKbPpmG1q+rD\nEnu2N2uFJd3NHixZCyxLz3phTfraICwzG1ZXUBVW93LaTtWA1R41sKZbqMHqtn7FsNrWLqz2\nUWBNtigS1mBbVgfLOdHAakaIdfaLgdWOcqvHw+q2eBJWZ7kHqzk+TVjN0Zrx7BWWWPu8OVim\n/WNTXQ/WFIfdwnJaAavtYq0aWCYbLNOe/zSw6qfOlcAS64x3sMZfDbcOq7GRGpZYr6bWnLbZ\nEFanpINVXxHNh9WuKAyWc4aBZQqF5axzDqz22IA1c8riFoth1eMaYsAaHEB3YNZxACsUlrHv\naFvlgeX0KhaWcxzvYEmvrF0tOsAqCFazx5+CZYC1U1gdFmCVD0uAtW9YprNhksIyY7BEDZZ0\nw/qnBViLWhQEq3vO6sFqlwos41YOzI5gNSeuB0t6fJwlA8s9H+EBFrCABSxg+VqkhSVOPeOF\nJWPXWAlgSTtjBJZzmMBKlzhYzd06sOyNt3cNWG7L/tLG7vSnSFjNeUkAS4DVHoddzz7GTcFq\ndMyB1VSdA6u7qQ2raZ4ZVrvWBbCGO9a/431WAcs+WbGw2o3YECwB1miLrLBMb38Xw3K2B1hT\n2SQsUYTljgHWVEqG1QCJh2WANRuWdTDD8xGesmE1g507w2FJOCx3nXGwJARW9xIKrETZGSyZ\ngGUfvRndPWAtbAGsSVgCrAUtVgSr139wWENYZgmsZu1tZWuJLSxJAqubPoA1eJEGlg3LqMOq\nt3MIq8OhAMu0z1rWzJSwhld/wMoMy2wCVq8lsEZgtcOWwXI3bRms5gkBWMnzcVjd9VYorP5i\nxw5kc7CkGzdYWv+O9wGW0224YmANjjwo+4Vlrc4Ly1j9B20GsMwELPtxowXLWgiwmhHA8sFy\nbE3BMtZx7RhWt6PbhSV5YNlDdw/L2tHNwuoWCiyFTJ1B67ANsIC1vIW0Z0obVn0qgQWsdk6i\nZ6ym4TphWRWA9b5FCCz7frFPVgAs6c5re48yLAeOPiyjBstsDpaR/pAgWCN/wz4Ay6wblgCr\nGyfdjI3AaovmhuX0ANa6YPUO7bOwrOtWZ5XAMplhychexMCyV6gCy0TBMoOD6e4ZqjKDwSEB\nVr9wbynSLXg2rMHygJUiuWGZTLDEAGs6n4H1OsRAWPJRWE0DHyxTDiznXAFrMMyB5fwt3DGs\n9kIPWOMtgDUTllVuCayRgwHWR2CZ8mD1DzoSVt8lsIA1ftBxsMYaA6utOguWmQtr7MXnk7AG\nK6ybd4WApQXLqhYCq8vIIQzuBJY3m4A1bDgHlv8QPgHLPh5gvW+xBJb1h0Ww7C+Q5sOqj2YT\nsNqrtNFHo1MCLPeM1UOau3RgGZmaPHUI24flHtZY9ZgUAWs4LATW9NYBa6+wzHxYzpkbbwAs\nYFnn7z0sEwireRRYwOoPm4Y1enOkJ7D2Dqs7Fe6Q5q6EsKSrvTNY9TnYF6yxYWuC1Xv0w7Dq\nW9GwpH9jpHpMgBV4CLGw7LnDeflhGWDVQ5q7tgprohGwnrkeRc4/kS2ao3e3ZG2w3IVMwxqM\nB5Z33nPi6fVvtpe4FpGwbCHA2gWsi1wextwvco1q0R09sEaKacGS6WMoDdZBHtXthxyjWgDr\nE7A8Xe2HioAlYv0hvAWwyoPlH5Mb1r/mcA5RLdpzBqyxYiOwrDnbh3X+uv7I99/Nx8V/9V4k\nrKkNmjyEtLCmGzXX2FaR5mudWbDsduuA1X4Tr8jhEdUiPyyzJlhdHffOWbDM2mCZ2+16PZ+r\n6YeL1xWweuPzwlrbS+GSFsDyNBqDJXNgNTN2D8sdkhqWZSMxLKMOyyyCFdraN2arsHp310WA\n1VuD2Rws/4KB5Y4GVniRoR07Iz2BNVVKRgbuFlZkC2B5Son9B2cKsN61yAFLgAWssTkLYbn3\nu7CMb4MmKxUOqxu5eli/X+fnFdT58hvZYuuwDLDmw3ocravzU1yLXLCcFQBrHbAucvi+PW/d\nfw6R/wjd3vchWM4kYAWMyQjrILf29i3y22ba+6Jh2fu3fVi9KfuA5W6St8pCWBO1lsIKnmw8\nk8dKAGvulGfyPWNN1NofLLfZVmH9XWP93J+3tK+xJmptE9bEGncEq/nJr2eOcd/oZ4AV92g4\nrPaB9cIyv5fn+1iH81fs+1gmMay3SQEr+FFlWG+XZH/ds0pYi1oAK+ZRYAW3AFbMo8AKbhEG\ny1cLWN4KwALW+0eBFdyiHFjBBUIfBdbcKSlavGsLrHd3AgtYpv+4rxawFrV42zZgALA8FYAV\nMak3AFieCsCKmNQbACxPBWBFTOoN2CSs4CnAKgBWXRBYwFoVrN5YYM2ckqIFsKLax8AywIqd\n1BuwUliJdglYwFpQbM6UaVhBs3cPK/KsZTjQuhGwZk5J0WI5rBSrUAmw5k5J0QJYcTWjBwMr\nV4AVOHvdsPIHWIGzgRUXYAXOBlZcgBU4G1hxAVbgbGDFBViBs4EVl5yrABawlHp9FFb7S1CB\nlSH7gdUtAFgZAqzAqcCKC7AWBFjTAdaCAKuMAAtYKgEWsFQCLGCpBFjAUgmwgKUSYAFLJcAC\nlkpifzQtoOCsBQBrYwEWsFRSCqyP9Z83pcAWhQVYwFIJsIClko/DSnd1NbM/sFQCLGCpJDms\n+BUAa4sBFrBUAixgqQRYwFIJsIClEmABSyXAApZKgAUslQALWCoBFrBUAixgqQRYwFIJsICl\nEmABSyXAApZKgAUslQALWCoBFrBUAixgqQRYwFIJsIClEmABSyXAApZKgAUslQALWCoBFrBU\n8nlYiQOsMgKsJbv++3WuPsBMzpdfrRZrTeIXos8nI6zHUbqcVFqsN8Cav+sXOXzfnrfuPwe5\naLRYb4A1f9cPcmtv3+Sg0YIUk4ywnL+U/r+hwFp9eMYiKsl7jfVzf97iGmv7yfl2w8n6qvD4\nUGlBSkne97Euz/exDucv3sfaenjnnagEWEQlwCIq+RQs3sfaeMqBJXZStCCfDC+FRCXAIiop\nFBZZfWbs+mwvwd/ol65l1porWWaph57hG/1StcxccyXLLPXQM3yjX6qWmWuuZJmlHnqGb5tJ\n1TJzzZUss9RDz/CNfqlaZq65kmWWeug8Y2UsuZaaH77GCv1Gv1QtM9dcyTJLPfQM3+iXrGXe\nmitZZqmHnuEb/dK1zFpzJcss9dA/8O8tOzq7K60JLM2aK1lmqYcOrIwl11JzpbDIHgIsohJg\nEZUAi6gEWEQlwCIqARZRCbCISoBFVAIsohJgEZUAi6gEWEQlwCIqARZRCbCISnLDuhzkcIn6\n2YuJXI9tIavm0vK/9flIVvL2T+TfPWnNx3ih+TWvDYKUhTPDev1sz3F5ocuz0OHh1lxa/nGQ\n/jKXlfxJv8z74VXznqrmrfl54/FqMwvnhfUrh5u5HSTy53qGucm/R/VX7Z9Tc3H58+scpyt5\n+Jv9OFc/d5ms5r/nT3Fekh3636QXgvFqcwvnhXWRn7//fsvX0kLn17qrU2LVXFr+u/5NUMlK\nfj8RPKqfFE9WU9Ie+lVOdcXxanML54V1luoJ/CbnRPWqU2LVXFj+3pzjZCX/db+IIFnN+tW6\nwpqi5h/9GtZ4tbmF88Ky/ralyKP6zVxWzYXlT3J/TU1W8ijm6/B81U5X86t+KfxKU/PWn9ur\nNrfwqmFdq6fphDv2bRLDEnn9tHjKmuZaXb0frulqAquX++Fs0u3Y89k+Oazq4v1fqmeXV76e\nX6Z9pVsnsNw8DqdezUXlj9WbAslhVddY9+qL9WQ1r9VL4R/WK7C6HFLCOh37NZeU//f86uc1\nNVFJZ1eS1TxKdcn2qLAmqllPGq82t3BeWK8vMe4pviq8H0/3fs0l5e1fPZ2opPOuSLKakrym\n81Vhv9rcwnlhfT2fFn7iflHbaH7aX9Vs1VxS3oaVqGQz+16tNVnN11PI872xZIduLbZfbW7h\nvLCSvfN+734FeMp33utznKzkvfqddH/XQ98Ja16k+oe7S8J382tYa37n3Rxn/F74sfyzPjLB\nqrm4fH2Ok5X8Gi20rOYpdc3m+mm82szCmWG9/mV+eR3rdcuuubh8fY7Tlfw5jRRaWHO00IKa\nDazxajMLZ4ZF9hJgEZUAi6gEWEQlwCIqARZRCbCISoBFVAIsohJgEZUAi6gEWEQlwCIqARZR\nCbCISoBFVAIsohJgEZUAi6gEWEQlwCIqARZRCbCISoBFVAIsohJgEZUAi6gEWEQlwCIqARZR\nCbCISoBFVAIsohJgefL6VXc/7wf+dKPJK5wMT+pfcfp23NH+2AHyDCfDkyeVAC+QGoZT4gmw\n5odT4snrs2TqT8esPoT6+rr3caw+qOHnLK9fJ1wPaYcd62H3sxyenx/5cxI5BVyqbSnA8sSG\ndW5/3fnzs+Iu9a9wr36Pvw3rZA071J/SdX2NvH72YDIHWJ5YL4U/cnqYx+n1AYmn1ydbfj8/\n7LcdUv33u/4ch+962PX1WUq36pEEH7G+ogDLE0vNuf7MrednGv6ODnl9bNLrk2dOzbDXU9nO\nXgarAMsTR037URjtpfr95+vUgzX4oL/qv5e/F8/bbVh+0wGWJ35Yp/ZDV97Aqj4X+u+K657/\nAD4YYHnSU2PfW31O1PH6cw+C9ffieDlyjUWaONdYP8699f/0YTXXWP2Pl3Zu7SL7OtrI1Gqq\n17Dnl3t/X+SdLVi/5nayhgy+KmxKHF9fP/KMReo0H4d4MM0VVXWhVMO61Bddv82QwftYTYnv\nduCOAixPnjZ+j09Y1Vvq8q95Zqry7w/Q7/NV7zWkfuf90L7z3v73+c77vlwBi+gEWEQlwCIq\nARZRCbCISoBFVAIsohJgEZUAi6gEWEQlwCIqARZRCbCISoBFVAIsohJgEZUAi6gEWEQlwCIq\nARZRCbCISoBFVAIsohJgEZUAi6gEWEQlwCIqARZRCbCISoBFVAIsopL/shj4b6eQ5bYAAAAA\nSUVORK5CYII=",
      "text/plain": [
       "plot without title"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "post = metrapolis_hastings(n=n, ybar=ybar, n_iter=1e3, mu_init=0.0, cand_sd=0.9)\n",
    "\n",
    "traceplot(as.mcmc(post$mu))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAJYCAMAAACJuGjuAAAAMFBMVEUAAABNTU1oaGh8fHyM\njIyampqnp6eysrK9vb3Hx8fQ0NDZ2dnh4eHp6enw8PD////QFLu4AAAACXBIWXMAABJ0AAAS\ndAHeZh94AAATCklEQVR4nO3djXaivAJG4fAjUitw/3c78h8UUSFvmjj7WevM57SaoO4CpXMa\n0wAC5q83AN+JsCBBWJAgLEgQFiQICxKEBQnCggRhQYKwIEFYkCAsSBAWJAgLEoQFCcKCBGFB\ngrAgQViQICxIEBYkCAsShAUJwoIEYUGCsCBBWJAgLEgQFiQICxKEBQnCggRhQYKwIEFYkCAs\nSBAWJAgLEoQFCcKCBGFBgrAgQViQICxIEBYkCAsShAUJwoIEYUGCsCBBWJAgLEgQFiQICxKE\nBQnCggRhQYKwIEFYkCAsSBAWJAgLEoQFCcKCBGFBgrAgQViQICxIEBYkCAsShAUJwoIEYUGC\nsCBBWJAgLEgQFiQICxKEBQnCggRhQYKwIEFYkCAsSBAWJAgLEoQFCcKCBGFBgrAgQViQICxI\nEBYkCAsShAUJwoIEYUFif1i/59y08uL3xRSI3ud57A2rTq1pM8kUCIbHsAqT/Fy7W9UlMYVi\nCgTDY1iJuU63ryZRTIFgeAxrcdjdPgYTVvTYY0HC7znWpepucY71/TyG1WTWd4VpLZkCofAZ\nVvNbdNexkvz86jrW7ikQCK9hhTQFtMIJ6+BlW4TlD8IqE5OW2inw53yGdc1NUjZnfqTzP/AY\n1rUrqjCnuqlys7nPIqzoeQzr1F67Kvoro7VJP5qC0mLj/Uc6Jrf+8vYUhBUb72H99MfAD3+k\nQ1ix8XooPI2X2+vThz/SIazY+PyHfsl0/DPbOyzCip/X61jFmFOyub8irC8QzpX3zSkIKzaE\nBQnCggRhQYKwIEFYkIgkrIZ/oxWZWMJipxUZwoIEYUGCsCBBWJAgLEgQFiQICxKEBQnCgkQ8\nYVFWVAgLEoQFCcKCBGFBgrAgQViQICxIEBYkCAsShAUJwoIEYUGCsCBBWJAgLEgQFiQICxKE\nBQnCggRhQcJ/WGVqTH7ZMQVhxcRjWP2vThsWht7+Re+EFT3fYRWmqJumKnYsK0dYMfEdVmK6\n9XQ+XVbuyLz4C77DGn+V6IfLyh2ZF3/Bd1inMazPlpU7Mi/+gtew8nN5MT+3m3Xx4bJyR+bF\nX/AaVq+7mdSbd3U5L/6Cz+tY12tZ5nl3Cl9sdkVY8ePKOyTCCcvYPM4LCZ9h1SdjsuGHOTsu\nN1BWTDyG1a4JfZP3gxDWd/MYVvdjnLpMsm4QwvpuHsNK+gdWSVrtDIuy4uH9XzfcdlpZti8s\ndlkR8RhWasaLV2m2LyzKiofHsEpzGm5VJiOsL+fzckMx1XR5sWAqYUXP6wXSaz7eqk6E9d3C\nufL+zhSEFY3IwqKsWMQVFrusaBAWJAgLEoQFCcKCBGFBgrAgQViQICxIEBYkCAsShAUJwoIE\nYUGCsCBBWJCILSz+qV8kIguLfVYsCAsShAUJwoIEYUGCsCBBWJAgLEgQFiQICxKEBQnCggRh\nQYKwIEFYkCAsSBAWJAgLEoQFCcKChNewfs95v7Jc8bt7CsKKg8ew6tRaQjXbOwVhxcFjWIVJ\nfq7dreqSmGLnFIQVB49hJeY63b6aZOcUhBUHj2Et/r+m+xZpOjA7/GKPBQm/51iXqrvFOdb3\n83m5IbO+K0zrrXsSVvT8XscquutYSX7mOta348o7JMIJy9i8zw7HfIZVn4zJLsMgXG74bj5/\npJP0PyjsByGs7+b1ckN5q6tMuh8TEtaX83qBtPtPlaQVYX29P/iRTp1lhPX1PIaVmvGiaJoR\n1rfzGFZpTsOtymT7w6KsKPi83FBMUVxe/FrtrU8SVhS8XiC95uOt6kRY3y2cK+9vT0FZMSAs\nSBAWJAgLEoQFCcKCBGFBgrAgQViQICxIEBYkCAsSMYZFWRGIMCx2WTEgLEgQFiQICxKEBQnC\nggRhQYKwIEFYkCAsSBAWJAgLEoQFiRjDoqwIEBYkCAsShAUJwoIEYUGCsCBBWJAgLEgQFiQI\nCxKEBQnCgoT/sMrUmPxyaArCCp/HsPrfuTCsOL65iD1hxc93WIUp6qapurULd09BWOHzHVbS\nLwFWm/TAFIQVPt9hjb+DaP8iTQe2AP74DmtcnMkkB6bgN2SFz2tY+bm8mJ/bzbrYPnt/NQVh\nBc9rWL3uZlJv3lW0CfDG53Ws67Us87w7hS82uyKs+EV55Z2wwhdOWMb2R9sAZ8IJ66MpCCt0\nhAUJwoLEH1xueOM0irCi5zGskrD+I16vYyWZqykIK3Rez7GuL/4Z1vtTEFbo/J68l+bqZgrC\nCh3fFUKCsCBBWJAgLEgQFiQICxKEBQnCggRhQYKwIEFYkCAsSBAWJAgLEoQFCcKCBGFBgrAg\nQViQICxIEBYkCAsShAWJSMPi99uGLtaw2GcFjrAgQViQiDYsygobYUGCsCBBWJAgLEgQFiQI\nCxKEBYmIw6KskMUbFrusoBEWJAgLEl7D+j3n3WonefHrYArCCpnHsOrUWklne/ETwoqex7AK\nk/z0C1NUl+TYKvaHtgM+eAwrsdY7uZrk8BSEFTKPYS0uPB1cVu7IdsAH9liQ8HuOdam6W5xj\nfT+flxsy67vCtD4+BWUFzO91rKK7jpXkZxfXsQgrZBFfeSeskIUT1tsrkas3BC6EE9bnUxBW\nwAgLEoQFCa9X3t8+jSKs6HkMqySs/4jPQ+E12f7HMp9OQVgB83qOdd3+Qc6nUxBWwPyevJfW\nz6GPT0FYAeO7QkgQFiQICxKEBYmYw6KsgBEWJAgLEoQFCcKCBGFBgrAgQViQICxIEBYkog6L\nssJFWJAgLEgQFiQICxKEBQnCggRhQYKwIBF5WCwBFqq4w2KfFazYw6KsQBEWJAgLEoQFCcKC\nBGFBgrAgQViQICxIEBYkCAsShAUJwoIEYUGCsCBBWJAgLEgQFiQICxKEBQnCgoT/sMrUmPzi\nbArCCpPHsPr/D2DWL9y7vSLmJ2FRVpB8h1WYom6aqjClqykIK0i+w0pM3d6uTepqCsIKku+w\nxiPX4xHs7SXuHW0OpHyHdRrDSlxNQVhB8hpWfi4v5ud2sy62z94JK3pew5oOc8YktaspCCtI\nPq9jXa9lmefdKXyx2RVhxS/6K++EFSbCggRhQYKwIPEFYfHTwhDFHxb7rCARFiQICxKEBQnC\nggRhQYKwIEFYkCAsSBAWJAgLEoQFCcKCBGFBgrAgQViQICxIEBYkCAsShAWJ7wiLtILzDWHx\nW/0C9CVhUVZoviIsTrPCQ1iQICxIEBYkCAsSXxLWf1BWZM+QsETWroDcf+ijjQ7mGb63IV8f\n1l9d4nIfVihlEVb/CTeb+/E12Nv9Hx7iOiyzevPlQJ/M+ubjV9aDcDOwY3u3am0v8c5Yr+8z\n/jzSLD/4fJYhLDvIl6Ftb8EHYfW/9nyt7Xd+sPriDu/siwnLnm/rrlthzYu0GPtmH1ZjrPut\nTvvqQ3dbMN/FPAnLDK/FeljbL4jZ3IZp+PtR/oOwHt689bHM+Kn+cYt1ftbvvxaWmYYw81xm\nXCqhGcKyWlvZ5pcfutuCp2FNGT2GtdhtPst8GtxM8axmRlj9x+/fjvmvVljz359MP45s7CmM\nmaKydph9UG+Gdb8i1ZMvA6ubcbC7XFbu8CSslS8js7hlpqdpVjZ6CGsx9PZh+l2Bh3Vf0BzW\n8qXs39Jx1/IY1vJkZBHW+MdmWIspGmuHMj+/sb1xM8btenxy1ptsdTMNMz4J+9HjaZb9Isyb\naO8B52c6PZ+7sMZRxhfMHtp+/PIF+1SgYZn5f/cfH4a8C6sZ353tsKYDwvByDm/e8EYuDhhW\nWIstmB5nPb85zLnvscfH5zb2N2+lMdNGzw9cPPY+rHEnOo7zLCwzHsvHx1gNrYW1tsVfGdbd\nsWV6HxanDOPbMf457TamV386slj7CzMd2qZHz91+GFbzPKw11r5x2v8Ze1sew5pfi3HP+CKs\nKSM7rPlJ34dlfWjl/Xj5jrl4iIcpxj2AFU9jnwrYe/bm7t2fD0nDp+yXzg5reuNXw7KnWm7Z\neBCZtnM8dE1fDtNnrQc/7mHnOOYvl2Z+oGnu3+IpN3O3idM+eNoDm3k77WGm/fLitbH3283D\nrOMUHwo0LOsr1t6PW+/jtF+4ewOnx0+3rK/o+WzD+rSx3qB5i9eatUY007DjUdBMb850jDV3\nY85/mXefxn6fm82w5kfY+9TFE5rSbsYvlZWwrL3xfGiYnuu3h2W/IuPX5JyQ/ZmHA+YiLPvE\nZLy13L89bueznaE9w7Tns3qfwpp2m9thjYU+hmWMeR3WvClPwpqyf5h5mnDc+d2/GPfzfizU\nsLo/7d3O8qtvGdbK4x/3M/PNxVFpNSz7ndp4se33czn3XNn6Y5/UOif37HHzwfCZ7bCW0/13\nYQ0PfB7WfJfH3VXzIqzmSVhPtuGdsB7mtvZB66M+n2/r8+NxcnOrrVfkyY7P2tLFs/hfwrJG\neBZW99l3Pvb0JXu1dZtvy/oD5rBejv7hpB8Pt36qNs0SZFi/57zbXeTFr2qKaYTxBH79qa+F\n9fZGOA/rrVFfcRjW9t4xuLDq1MwyyRTLQeYzAbcT7A9r7+7oHQ6H3gjrfq4QwipM8nPtblWX\nxBSKKQKx+yCqmHTPUO/vrEMIKzHX6fbVJIop/mcuw/rgrgGEtXJxZPlpy84p4ERkYbHH+kIh\nhHU7x7pU3a1vP8f6j+z+1tnNQ3qZdaxLa8kU8CyIsJrforuOleRn+XUs/LFvvfKOP0ZYkCAs\nSBAWJAgLEoQFCcKCBGFBgrAgQViQICxIEBYkAg0L0dvxrrsP6S+mFIwZyWaG+tQJy+OQsYxJ\nWMoxI9nMUJ86YXkcMpYxCUs5ZiSbGepTJyyPQ8YyJmEpx4xkM0N96oTlcchYxiQs5ZiRbGao\nT52wPA4Zy5iEpRwzks0M9akTlschYxkz0rDwPyAsSBAWJAgLEoQFCcKCBGFBgrAgQViQICxI\nEBYkCAsShAUJwoIEYUGCsCDhO6wiMUmxuULKm8p0Gsga8+jwv8Pr4WzI68mYU+V0zHp9oP1j\nlmMELgf2HFa/Ak96fKCiGyipl2MeHb5OzP1mHhvy4n4zq6Qfs3I15nX8ZTLro+0c2G9Yvya5\nNtfEvFh957WrOdXtl9ppMebh4fP+NXY3ZHJ7dJ23q6M5G/PUrbVWOHvqtwf1EayPtndgv2EV\n5nL788ecjw6U99vdviTWmEeH/xl+E5SzIX+6COp2PUdnYxq3T7002TDi+mh7B/YbVm7aHfjV\n5I7Ga18Sa8yDw1fja+xsyNO8XKizMYejdRurizFv6Q9hrY+2d2C/YVlfbS7UJluMeXD4zFT9\nQ50NmZrmnHRHbXdjnodD4dnNmNf7x96NtnfgqMMq2920w3fsp3EcljH9mo4ux2zK9uw9Kd2N\nSVh3qiRv3L1j3d7eeVjtyfvJ1d6ld+6+TTu7207CWqqT7G7MQ8On7UUB52G151hV+826szHL\n9lB4i7UkrFniMqwsvR/zyPCn7ruf/qGOhly8K87GTE17yla3sToac3jQ+mh7B/YbVv8tRuXi\nu8Iqzar7MY8Mb//qaUdDLq6KOBvTOB9z8V3h/Wh7B/Yb1rnbLVy672uOuZjsccwjw9thORpy\nfHTVbquzMftdSHdtzNlTtzb2frS9A/sNy9mV92rqyumV9+E1djbk7eyqbs+HfhyOWZj2B3eF\nw6v5Q1gxX3m/nSC0std3fOFkLZlgjXl4+OE1djbkeXWgY2Nmrsccz5/WR9s5sOew+p/MHx/H\nXovDGvPw8MNr7G7IS7Yy0MExVwc6MOYY1vpoOwf2HBb+F4QFCcKCBGFBgrAgQViQICxIEBYk\nCAsShAUJwoIEYUGCsCBBWJAgLEgQFiQICxKEBQnCggRhQYKwIEFYkCAsSBAWJAgLEoQFCcKC\nBGFBgrAgQViQICxIEBYkCAsShLWh/1V3l9d3vMz3Ro8XY8PwK05f3i+1lx1AhxdjQ5fKG72Q\n1CNekg2EtR8vyYZ+LZlhdcx2Eeqy/2idtgs1XHLT/zrh4S7T3dLhblVukm79yEtmTPbGqdo3\nIawNdlj59OvOu7XiiuFXuLe/x98OK7PulgyrdJX9Pcu/fTKeEdYG61B4MVnd1Fm/QGLWr2z5\n0y32O92l/fNnWMfhZ7hb2a+ldG0/42CJ9YgQ1garmnxYc6tb0/B39S79skn9yjPZeLd+V/af\nHQZbhLVhUc20FMZ0ql5dztldWA8L/bV/FreD5/X6OPxXI6wN22Fl06IrL8Jq14W+nXFV/p/A\nHyKsDXfV2B9t14lKy0v1Vli3g2ORco6F0eIc67L46PCf+7DGc6z75aUXt/4L/9ez/dBQTXsM\n677du32Tl1th/TbXzLrLw3eF4xBp//0jeywMxuUQk2Y8o2pPlIawiuGk63e8y8N1rHGIn+mO\n/xHC2tC18Zt2YbWX1M1p3DO1TreAfrujXn+X4cp7Ml15n/7srrz/X10RFjQICxKEBQnCggRh\nQYKwIEFYkCAsSBAWJAgLEoQFCcKCBGFBgrAgQViQICxIEBYkCAsShAUJwoIEYUGCsCBBWJAg\nLEgQFiQICxKEBQnCggRhQYKwIEFYkPgHzYwYHZ0QTnEAAAAASUVORK5CYII=",
      "text/plain": [
       "plot without title"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "post = metrapolis_hastings(n=n, ybar=ybar, n_iter=1e3, mu_init=30.0, cand_sd=0.9)\n",
    "\n",
    "traceplot(as.mcmc(post$mu))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAJYCAMAAACJuGjuAAAANlBMVEUAAAAA//9NTU1oaGh8\nfHyMjIyampqnp6eysrK9vb2+vr7Hx8fQ0NDZ2dnh4eHp6enw8PD///+UfXg/AAAACXBIWXMA\nABJ0AAASdAHeZh94AAAWy0lEQVR4nO2d62LqKhBGadSqp2617/+ypyZao83NwDdhcK0fre2W\nTIauDYQgCd8AAsLSJwBlglggAbFAAmKBBMQCCYgFEhALJCAWSEAskIBYIAGxQAJigQTEAgmI\nBRIQCyQgFkhALJCAWCABsUACYoEExAIJiAUSEAskIBZIQCyQgFggAbFAAmKBBMQCCYgFEhAL\nJCAWSEAskIBYIAGxQAJigQTEAgmIBRIQCyQgFkhALJCAWCABsUACYoEExAIJiAUSEAskIBZI\nQCyQgFggAbFAAmKBBMQCCYgFEhALJCAWSEAskIBYIAGxQAJigQTEAgmIBRIQCyQgFkhALJCA\nWCABsUACYoEExAIJiAUSEAskIBZIQCyQgFggAbFAAmKBBMQCCYgFEhALJCAWSEAskIBYIAGx\nQAJigQTEAgmIBRIQCyQgFkhALJCAWCABsUACYoEExAIJiAUSEAskIBZIQCyQYCBWAPfM+Kun\nF2mBEKAFsUACYoEExAIJiAUSEAskIBZIQCyQgFggAbFAAmKBBMQCCYgFEhDLiHk3/P2CWEaE\n7yLT6gWxTLi2VsXl1Y+lWOfPENaH60EGj1LcHyA8fHsHDMU6V/XKwk1zkLcSK/x5UTyGYm3D\n/seufbWuD/KeYhWXWS+GYlVNwVO1OiFW8RiKdXPpvF6/mVih82XZGIq1Cufbq/VbiRV6fygY\nQ7H24fP66hTWiFU4ltMN21+bDiPT0GXVPmLJilw5bm6vTp/vI1YY/LFUmHmXg1i6IhmGsOM5\nmaKS6wWx5CCWrkjHQd5mjPUnl5KS6ycfsSK3KskWxBIWyTCEFR2pFJRdP4glBrGURTIMYQVi\nKYtkGMIKxFIWyTCEFV2pFJReL4ilpTOTctLrx3Q91uQZhXJqvjuTcvLrxXTZDGIN/rYoLLvC\nY7PcXRkiOxBLWuTKMWzVIXIDsaRFbuzDUR0iL3oSKSa/frgqlNKXSDEJ9oJYUhBLWyTDECb0\n5lFKgv0glhLEEhfJMIQJiCUukmEIE/rzKCXDXhBLCWKJi2QYwgTEEhfJMIQFQ2kUkmIviCUE\nsdRFMgxhAWKpi2QYwgLEUhfJMIQFiKUukmEIC95kMWMniKVjOIsycuwFsXQglrxIhiEMQCx5\nkQxD6BlJoogc+0EsGYilL5JhCD1jSRSRZC+IJQOx9EUyDKEHsfRFMgyhZzSJIrLsA7FUjOdQ\nQpa9IJYKxDIokmEIOYhlUCTDEHIQy6BIhiHkIJZBkQxDyJmQQwlp9oFYKhDLoMiVf7tN8yD7\n7T9ViGyYkkIBafZiKNZ51doocnhvvwJqHLEsitRsQ/XV7Lt2OlTDe/sVUOOIZVGkpmpt53cM\nlSJEPkzKwH+a/RiK9bBRcum7Jk/LwH+evdBiaUAskyI1P2Osw6l+9QZjLMQyKdKwbl0Vrs6S\nENmAWCZFrvzb1vNY1WZX/DzWxAz8J9oHM+8aEMukSIYhtExNwH2ivSCWBMRaSqzC57G6Evj4\nmPa+MshHrMnPnPNARwIfF7We3XKfaC90hRL+JvDx8G3gjYWAWBL+JPDx9L33jaWAWBKeE/jo\neNX5xmJALAlPCXz0vC4g0z4QS0F/g4VYSYtkGELKQIP19JP3THsxXY/1Nk+xRyxLsfZvI9ag\nV29ilmVXeKyGP0KRIEQeIJbxGOs4vLwvRYgseDz9v3dy2r9xnmovtoP3fWt1sihEDrwilvdc\n++CqUMDD6Xfcen6HJguxBCAWYklon36XVw+/dZ5rH4glALEQS8G4V4iVqEiGIYRMEKv1e9+5\n9oJY6XlNLOfJ9oFY6WmdfZ9XiJWmSIYhhEwRq/i+ELHScz/7fq8QK0mRDEMIQaxvxBIwqSds\n/ZvrZHtBrORMarCKb7IQKzkTxSq8yUKs5CDWBcRKzcQhFmKlKJJhCBlTG6zCzUKs1CBWDWKl\nBrFqECs1iFWDWKn5Pfcxr+7v8JxuH4iVGsSqQazETO8Jy+4LESsxiNWAWIl5oScsui9ErMQg\nVgNiJQaxGhArMYjVgFiJeUms25sc59sHYiXmduqTvCq4yUKstLzWYCFWZJGG02eodt/f+1Wo\nRjZg81vPiHXFUKxzddl6dL+rdyAd3jTSbz2/KFa5ZhmKtb1sFLmtwuf5+7wt9ZnQrw2xECuy\nSE1VFwyhfhp0qU+xR6wrhmKFcP9a7PMKEevKAi3W5euZFuvxjX4z7mGBMdb2fH2dPsTyINYV\nrgqT8rJXiBVVpOEN5rFmi+U35R6YeU8KYt1ArJS8Oj3aeqvXlPtArJS83mAhVlSRjoMUOY+F\nWL/kI9bkhxlmDGL9QleYkjliFWoWYqUEsX5BrJRcz/slrxArosiVf7tNPYLabP+pQiwLYv1i\neUtn1RqdF3lLZ1ZPWOgUqelN6OqreXLv6VAVeRN6nlhlNlmmy2buD4Q+FrlsBrHumC/06/oh\nWYiFmTXEQqyIIjVv02Ih1rf1GOtwql8VPsZCrG/b6YZ166pwdZaEWJaZYhVplu081raex6o2\nuzLnsRDrDjPv6ZjrFWLNLpJhCAGI1SJSrNXulOxUekL4AbFaRIp1GYYL3HJZxTOnR+9FXGbd\nR6RY569PhVsuqxix2iQYY/3brVK75bKKZ/eERfaFaQbvx8tnUffxZzMQIn8Qq00SsQ7rCUth\n4kI4ALHaxIt13v00V6vD+ceuTZpzclrDiNUmVqx/l8H7trm7nO7DNS5rOEKsAkfvsfNYP43V\n/nbbb3jFwtwQbmhOepZXiPVcJGwOyU6lJ4QXYhqsAvvC2HmsZCfSG8ILiPVA9Mz79UWVrBt8\nDuEFxHogkVintJ+K91jBUWKVZ1aEWIeH3RZWC5/V4iDWAzEtVvtzgquRpXvys1qcmItCxHou\nItoXxmP9ItYDLPRLRFxPiFitIpfWSrSllcP6jRSruClSxEpEGrE8Zt4NXWEi4oZYiDWvSIYh\nEhPbYCHWU5H96vv7tEo82+CwehHriUixDpexVf0okzefx4oWq7TRe6RY6/D1fQyr76+Ey0e/\nPVYvYj2RYIL0eNng492vChHriQRibcIBsWIvCosbZEV3hcfDZeEoXWH9FbF+iR+8h7C7NFhJ\nV5K6q934nhCxnqYbmi3UVl8TCp63l9WAl4+3rkfe7q52E4hVmFmGE6Sn6mcg1jxmtbjtuBHr\nGUOxPsPm/PPl81Q/a7WsrSLjh1iINTvzEM7XL+U9xR6xnokVa7eavLqhfksVWj+kPKtlQaxn\nIsXavbBs5vOyHfeu2ZP7PDzIcle5KcQqa4o0UqzqhT1mjqHaHr831Y9Zh9Xw9IS3yk0xdi+s\nybJc836o7u3bLvVZLQpi/SFSrE146bPQX5/1kGwztkubt7pN0hMiVuv1qVqnXYn1N4QHEOsP\n0V0ha95T9YRljd4RKwGI9Zel1rwXNY+VVixv2XeTj1gPO0GkCGFHmiEWYj38dNjUi/3SbvTu\nrGpTiVVSXxgr1rppX0KV1CxnVZtYLG/pdxIp1j6szxex9uEz2Sl9u6tZxPpL9C2dczNcmjQs\n+rern1cYNtuinleYauyOWK3X4XuyWOf2flolLfRLJ1ZBg6xIsVbXFus4YUe/bai+mg3hC3sm\ndLKeELF+uY6xDlNWORT7FHvE6iD2qnAzqWtryk1fFuGrZpOL5Sz/TpLMY4XNlA/p0GJNoJgm\ny3Dm/WeMdWhmu8oaYyUcuyPWLNbtXZYH13G5qtiUDRZiNRzqhXvrsWmpG/+29ZCs2uxKmsdC\nrC5ixDrdm6D1G98rTCpWcxxX+XcTIda5ujz/8ufF6WuV7olyDyE8kHSIVU6TFSHWtjXHsB75\ndMTMEB5I22Ah1mXW/d7/nd53GyPE6iRCrBcmPOeG8IBELFc10AlixZJYrFKaLMSKpT7XdF4h\nFmLVpG6wEOvx0w+IlQrEQqwLiNWN5b3CrEKkIrlYhZiFWJGkHrsjlhJHtYpY3SBWJIjVDWLF\nkX6IhVhC/NRq+gYLsYT4qVWBWGWYhVhRCHpCxBLiplIRqw/EikIolp9K6ASxolAMsRBLh5s6\nlYhVRF+IWFEgVh+IFYNkiIVYOrzUqVQsN7XQCWLFoOkJEUuGlyoViVVCX4hYMSBWL4gVA2L1\nglgRiMbuiDU76NghnFSpWCwv1dAJYkWg6glLaLIMxXrh42I+alTWYCHWS/yrEGsq/vtCy67w\nvLnu/FdIV4hYA9iOsb5CuGzcjVhjINaLnNZhcy5KLIlXBQyyzK8Kd6E6lCGWsMFCrBkcV+M7\niLioUAOxfFREJ0vMY30i1iiIJcFFfSLWEIg1G+XY3f8gaymxCpggRawh8hFLtj2gCGlPiFga\nPNSniVguaqITxJrJ9RQRqwfEmom4wXJvlqlY/3bNI6Q3Yw84dFCbiDWMoVjnVWt0PvxMJwe1\niVjDGIq1DdVX87jxEp4JjVjDGIpV1FPs5V55n3AwXZrc90OyEHYg1gi0WPMwE8tBXXRiO8Y6\nNM9kLWCMhVgjWE43rFtXhauzJIQZiDWC7TzWtp7HqjY77/NYBl45H2Qx8z4LxBoDseagvlH4\ncPjcK6MbxJoDYo2CWHOwEcu1WYg1B5MhFmKlJ/eqRKxREGsGRj0hYiUn86pErHEQawZGPaFr\nsxBrBog1DmK1w14/eTYS3qonRKzk2NVkI9LzRxlvgnWfB2JN4K3FmrDtTccbbr8yEMvx7cK3\nFmsiofNHC68cN1nvKtZLn+J/6iebb4g1yHuK9fLmEO0CiDWFtxQr7vCWYvk1693EitnIpilr\n6hVipUUUInp/pMsBEGsSbyVWEhBrEm8jVrJDhqbhMxLLrVlvIla6PQKbAwUrrxArKalDJDye\nbU+IWGnJuBoRaxrli5V2p9z2jUKbP7bT24Wli5V6A+aHBstkd2enTVbhYqkGa7ee0EAtxEpI\ntrVouWTmMVS2VdJDuWJJGpNOscQV4nOQVapYoj6q+5pQ2yH6bLIsxTp/hrA+XA8yeJToSlTd\nbLx+/9MTKtVCrBHOzVPsN81BtGKJWGCI9e3ULEOxtmH/Y9e+qrd4F4ol/AMMi6UKjFjDVE3B\nU7U6CcVSdkpjDZYoNGKNlLsWPK/XOrGktT/aE2qsRqxhVuG2oe1qrRLLxCvjIZbPCQdDsfbh\n8/rqFNYeB+/TxBKcu8cmy3K6Yftr02FkkfCcEPq7KxOvCdOfCGKNcNzcXp0+E4tlcNNuek+Y\n/GQc9oWFzLxb1Pkyk1itkIgVTZZV+KJYSXPwZ1YBYhlV98vXhCnPy98oaymxkg3eTRbb1YGu\n36f3hAlPDbEmH+TPUUKb+YdRMWsSK93puesLC+gKbVhsdrQdNcN66cOvWHaNVRPu+v11r5Kc\np7u+0KtYxlpFNVhJztVbk2Uq1r/dplmStY19XqF5BcdNYiVQC7F6Oa9ao/O1JISOBWdH25Hz\nq5g+TBf6VV/N48ajnglt3Qk2QW8vIsSKPG9nTZbpQr8ET7FfRKtEl4Rx5+5s+L7AQr+/P0wP\nsVC1pppriFLLV5PlrcVahoUnsdqxc6uaPmzHWIdT/WrWGGuZTrAJfXuRxKu5eSBWH+vWVeHq\nPPTOvyEW1Cp5gzU3F1dm2c5jbet5rGqze3Uea9naTN4TzlMLsaLJq/JEI6zX93D2NJWFWKPc\nTyb50P1FtTw1WYg1ilCsF/E0lYVYY8i9eiFZR00WYo2hb7Cmd4iIFUlOVWfRE05Vy1FfiFgj\n2I2wpiTtp8lCrGFaZyIfuk9otxArjnxqzvaSsHnE+dA73JiFWINYNlj3oAPpuxllIdYgC81h\nDajlpclCrCEWabA6grdArBhyqbdlJ927P7nrxCzEGmDRBqvrLC4gVgSZVFsmdwmfWi4fZiFW\nP3k0WDd+5fJxYYhY/WTSYP1ya7hcNFmI1UfIq8H6JQTEmk0OlZapVz+Ej/rccqikfhBrwilk\nJtZ1iXIGlTQAYo2fQm5e3Z8cnEE99YFY42eQn1jt8ftLOyDagVijJ5ChV10XhpnphVid8XMX\nq3vKoVl1k4dgiDUWP0uvhiezxpd16UGskfCZejVh/v069lqoLhGrI7oHsabPvy8jGGINB8/W\nq9fv7NwEs6lcxBqMnbFXPyc3cy8HE78Q6zm0i37wStQuIdoGzFSsdNtxy/DTXtWk2H/mJlja\neQpDsTxsx+3Mq+Yc0wrxOxSLO6qhWIm24xYS3Hkl3TProSV7tUUzFCv7zW0fh1c+vDLfjW2q\naIZipdiOW4jH5qpm0X3+ev2ixbpGfNTKkVfXs83jDuEd2zFWzHbcSjxrdeFyxpmZZTndELMd\nt5DgXasL2allO481eztuGc9DBJ9aXfg580xWzNSUO/M+Wst/x50fd606S0f+4X6Ldw94W7+8\nv3wl5K3Vml5GaGKpYg3fruh4ovlHy6ru0pF3QH6Ldz9QvXX0+8tXQ37cmq1JZcq5pWMYInQd\nJ3T/TR+l6i8dd2q/xXvEun+9v5wR8pLLx7TJzMiEphxcXaTjINp5rFZP0ubxTR83ekuP/jLq\njHqOfn85J2T93o8W42f0wtFfPA95kY6D/DlKW4D/nvlQ8CfKjfD0vf+X0/kt3p1l6+j3l3NC\n/ikjqboJlbmQWOIQcf8b3bdYk8qU2WKpQ8SNHxyPsV4oU+QYSx0i7orH61Xha6dZzlWh6UK/\nuDkaj/NYHQdK9s7XT8KkSI2HhX6QCkOx8l/oB+kwFCvrZTOQGEOxMl/oB0mhxQIJtmOsXBf6\nQXIspxsyXegHCmznsfJb6AciSp15h4XJVCxwz4y/enqRsolZWhxXCSGWnziuEkIsP3FcJYRY\nfuK4Sgix/MRxlRBi+YnjKiHE8hPHVUKI5SeOq4QQy08cVwkhlp84rhJCLD9xXCWEWH7iuEqI\npQcgAbFAAmKBBMQCCYgFEhALJCAWSEAskIBYIAGxQAJigQTEAgmIBRIQCyQgFkhALJCwkFh7\nedxtFart4L5dqdDn0oRZmSR0/gzh8zj+vjGWEesof8hjsy/cShzlgj6Xmm2dUCU3q6rjxJu1\niFjHSv3H+Beq4yXMyJ5wCdDn0oQJn+dL4/gpjrO9RNiGTfSBlhBrH9bqP8Y2HH6+foWdNoxJ\nLjWbJog8VhXOacIs8mGKrbyCNuGy8+4xwf+8EQxyeQhnE2t4G+xph0hwGq9y1FdQMPoPbpFL\ni/PI02USsQ376GMsdFVYjFhWQRr2dQ8v5isMb68+DcRKFcuAUyXv2n/Yb6oEY1PEShVLz7ky\n6Qh/+IzvC00Hnvf9d9V/jKpEsdYW03I15/jRe6FiNVeFJ/1V4beZWKfV+mQS6EJ8ToV2hbt6\nlHtIMQodxUasg80FYTOPdYq/Z1GoWHYz70ZinWy8ambezxtXY6yHsOq4zZOFTf4YJmJ9zn9G\nxGtUiSquVLHO9eoGcZAGoysEI7Euy0JW8fOjrMcCDYgFEhALJCAWSEAskIBYIAGxQAJigQTE\nAgmIBRIQCyQgFkhALJCAWCABsUACYoEExAIJiAUSEAskIBZIQCyQgFggAbFAAmKBBMQCCYgF\nEhALJCAWSEAskIBYIAGxQAJigQTEAgmIBRIQK45m70bL5zQ5gRqJA7F6oEbiQKweqJE4EKsH\naiQOxOqBGokDsXqgRuJArB6okTgQqwdqJI4Q9t/HXRVOJ7OHCfoAseIIm8uTZw5VWNk9TNAF\niBXFkfrrgYqJYkf99UDFRFFRfz1QMVFsrJ4r7w7EAgmIBRIQCyQgFkhALJCAWCABsUACYoEE\nxAIJiAUSEAskIBZIQCyQgFggAbFAAmKBBMQCCYgFEhALJCAWSEAskIBYIAGxQAJigQTEAgmI\nBRIQCyQgFkhALJCAWCABsUDC/yp9/hNhKzmyAAAAAElFTkSuQmCC",
      "text/plain": [
       "Plot with title \"\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "post$mu_keep = post$mu[-c(1:100)] # discard the first 200 samples\n",
    "plot(density(post$mu_keep, adjust=2.0), main=\"\", xlim=c(-1.0, 3.0), xlab=expression(mu)) # plot density estimate of the posterior\n",
    "\n",
    "curve(dt(x=x, df=1), lty=2, add=TRUE) # prior distribution for mu\n",
    "points(y, rep(0,n), pch=1) # individual data points\n",
    "points(ybar, 0, pch=19) # sample mean\n",
    "\n",
    "curve(0.0141*exp(log_gdist(mu=x, n=n, ybar=ybar)), from=-1.0, to=3.0, add=TRUE, col=\"cyan\") # approximation to the true posterior in blue"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message:\n",
      "\"package 'pracma' was built under R version 3.4.4\""
     ]
    }
   ],
   "source": [
    "library(\"ggplot2\")\n",
    "library(\"pracma\") # for linspae"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = linspace(-1.0,3.0,n)\n",
    "x1 = linspace(-1.0,3.0,10000)\n",
    "y1 = dt(x1, df = 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message:\n",
      "\"Removed 1 rows containing missing values (geom_bar).\""
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA8AAAAJYCAMAAACaSn8zAAAASFBMVEUAAAAA//8zMzM/Pz9N\nTU1ZWVloaGh8fHyMjIyampqnp6eysrK9vb2+vr7Hx8fQ0NDZ2dnh4eHp6enr6+vw8PDy8vL/\nAAD///8ImqqkAAAACXBIWXMAABJ0AAASdAHeZh94AAAgAElEQVR4nO3dh3biyBaFYXmGcXY7\ntC+8/5teEOlIKFQ4p3RK/Hut6TGhtwtRX5cQqdkRQqpNs/QACCHpATAhFQfAhFQcABNScQBM\nSMUBMCEVB8CEVBwAE1JxtAD/nczMxRFRa3I4JJo0ipQmdC0BME21NQFYBMA01dYEYBEA01Rb\nE4BFAExTbU0AFgEwTbU1AVgEwDTV1gRgEQDTVFsTgEUATFNtTQAWATBNtTUBWATANNXWBGAR\nANNUWxOARQBMU21NABYBME21NQFYBMA01dYEYBEA01RbE4BFAExTbU0AFgEwTbU1AVgEwDTV\n1gRgkVDAm+tP+8j/H1PovgMwTQCWCQR8tbo5/bG5nGhT6L4DME0AlgkDvNkBmCYvTQAWid6F\nBjBNCzcBWEQD8D/7aI+LEBIQVuDFmraHqDTdZtVNrMAiAF6o6Yi3K3jpMVXSBGARAC/TdJEr\nCTvcTB6bACwC4GWatgM/LT2mWpoALALgRZqEWrEEO9xMHpsALBIH+PAnr8TKb+o89L2ecLiZ\nPDYBWITXQi/Q1Dv4DGDNIqUJXUsAvEDTduS0w83ksQnAIgAu37S9AbxNbBrNqpsALALg8k19\nv5dzHG4mj00AFgFw8aabBfhylsPN5LEJwCIALt506/d8nsPN5LEJwCIALt4EYNMipQldSwBc\numnI7+lch5vJYxOARQBcugnAtkVKE7qWALhw08AhrMvZDjeTxyYAiwC4cNOw3+P5DjeTxyYA\niwC4cBOAjYuUJnQtAXDZpjG/7T60w83ksQnAIgAu2zQK+C+AdYqUJnQtAXDRppFDWO1FLjeT\nxyYAiwC4aNO438NlDjeTxyYAiwC4aBOAzYuUJnQtAXDJpok9aAArFSlN6FoC4JJNU373lzrc\nTB6bACwC4JJNALYvUprQtQTAJZsAbF+kNKFrCYALNk0+BJ69OCYONziATQLggk1zQAGsUKQ0\noWsJgAs2AbhAkdKEriUALtc0u4ustw/tcIMD2CQALtc0y3MH4PwipQldSwBcrgnAJYqUJnQt\nAXC5JgCXKFKa0LUEwMWa5h/h7tQeBDvc4AA2CYCLNc3j3Kkdh3a4wQFsEgAXawJwkSKlCV1L\nAFysKQiwkmCHGxzAJgFwqaYAmju113I43OAANgmASzUBuEyR0oSuJQAu1QTgMkVKE7qWALhU\nUyBgHcEONziATQLgQk0hMg9NAM4sUprQtQTAhZpCYAJYoUhpQtcSABdqAjCALQLgQk0ABrBF\nAFymKejgVNukItjhBgewSQBcpimIJYAVipQmdC0BcJkmAAPYJAAu0wRgAJsEwEWawl6fsQu/\nakiTShw2AVgEwEWawlDuIq4b0qQRh00AFgFwkSYAA9gmAC7SBGAA2wTARZqiACsIdrjBAWwS\nAJdoCjR5agJwTpHShK4lAC7RFEgSwApFShO6lgC4RBOANZsALALgEk0A1mwCsAiASzRFAs4X\n7HCDA9gkAC7QFAry3ATgjCKlCV1LtACTiWyNr0/uN6zABZpYgVWbWIFFAFygCcCqTQAWAbB9\nU/BBqUtTtmCHGxzAJgGwfVMwRwArFClN6FoCYPsmAOs2AVgEwPZNANZtArAIgO2bAKzbBGAR\nAJs3hb+w6tqUK9jhBgewSQBs3hSOEcAKRUoTupYA2LwJwMpNABYBsHkTgJWbACwCYPMmACs3\nAVgEwOZNKYBzBTvc4AA2CYCtmyIoAlihSGlC1xIAWzcBWLsJwCIAtm4CsHYTgEUAbN0EYO0m\nAIsA2LopDXCmYIcbHMAmAbBxUwxEACsUKU3oWgJg4yYAqzcBWATAxk0AVm8CsAiAjZsArN4E\nYBEAGzclA84S7HCDA9gkALZtinLYaQJwWpHShK4lALZtimIIYIUipQldSwBs2wRg/SYAiwDY\ntgnA+k0AFgGwbVMG4BzBDjc4gE0CYNOmOIXdJgAnFSlN6FoCYNOmOIQAVihSmtC1BMCmTQA2\naAKwCIBNmwBs0ARgEQCbNgHYoAnAIgA2bcoCnCHY4QYHsEkAbNkUabDXBOCUIqUJXUsAbNkU\nSRDACkVKE7qWANiyCcAWTQAWAbBlE4AtmgAsAmDLJgBbNAFYBMCWTXmAMwQ73OAANgmADZti\nAQJYoUhpQtcSABs2AdikCcAiADZsArBJE4BFAGzYBGCTJgCLANiwKRtwsmCHGxzAJgGwYVMu\n4PQl2OEGB7BJAGzXFM0PwApFShO6lgDYrgnANk0AFgGwXROAbZoALAJguyYA2zQBWATAdk35\ngJMFO9zgADYJgO2aAGzTBGARAJs1xeMDsEKR0oSuJQA2awKwUROARQBs1gRgoyYAiwDYrAnA\nRk0AFgkDvNlH/Nie2IgzATwQDcCpgh1ucACbJAjw5vKHOGPTuUqh+w7AyU2JcdgEYJEkwDeg\nAXybBHoAVihSo1FHkgF3/QL4JgC2agKwSArg48nrQ+B/9lEfWPXZOmohq0064O55hf7xZQVO\nbkqMwyZWYJFkwL2fCt13dwc4UbDDDQ5gkyQA7h/NalPovgNwelNaHDYBWCQZMLvQ000p8ACs\nUKSpo4JkARZrcaH7DsDpTWlx2ARgkbhXYknJnRdiAbgfAJs1AViE10IbNQHYrAnAIgA2atIC\nnCbY4QYHsEkAbNQEYLMmAIsA2KYp6UsVAKxQpDShawmAbZr02AE4rkhpQtcSANs0AdiuCcAi\nALZpArBdE4BFAGzTpMhO7dF0Uhw2AVgEwDZNALZrArAIgE2a0r7ZF8AKRUoTupYA2KRJ8/VT\nAI4qUprQtQTAJk0ANmwCsAiATZoAbNgEYBEAmzSpAlb6aIC0OGwCsAiATZpU30ME4JgipQld\nSwBs0ZR2EBrAGkVKE7qWANiiSfezJAEcU6Q0oWsJgC2aAGzZBGARAFs0AdiyCcAiALZoUgas\n8jVpiXHYBGARAFs0KX+fAoAjipQmdC0BsEUTgC2bACwCYIOmRL8A1ihSmtC1BMAGTQA2bQKw\nCIANmgBs2gRgEQAbNKkDji50uMEBbBIAGzRpA44vdLjBAWwSABs0Adi0CcAiANZvSnwrA4BV\nipQmdC0BsH5Tql8AaxQpTehaAmD9JgDbNgFYBMD6TQC2bQKwCID1m/QBR1c63OAANgmA9ZsA\nbNsEYBEAqzclH4QGsEaR0oSuJQBWb0r2C2CNIqUJXUsArN4EYOMmAIsAWL0JwMZNABYBsHqT\nBeDYUocbHMAmAbB6E4CNmwAsAmD1JgAbNwFYBMDaTel+AaxRpDShawmAtZsAbN0EYBEAazcB\n2LoJwCIA1m4yAhxX63CDA9gkANZusgEcWetwgwPYJADWbgKwdROARQCs3JT+VgYAqxQpTeha\nAmDlpgy/ANYoUprQtQTAyk0ANm8CsAiAlZsAbN4EYBEAKzdZAY4rdrjBAWwSACs3Adi8CcAi\nAFZuArB5E4BFtACTY7YVNpOKwwqs25SzALMCaxQpTehaAmDdJgDbNwFYBMC6TYaAY6odbnAA\nmwTAuk12gKOqHW5wAJsEwLpNALZvArAIgFWbct7KAGCVIqUJXUsArNqU5RfAGkVKE7qWAFi1\nCcAFmgAsAmDVJlPAEeUONziATQJg1SZLwDHlDjc4gE0CYNUmABdoArAIgFWblAF/dLL9iE/W\neAbHtHwTgEUArNmU9yzSAOD/ZLb/RQfAaw+ANZvy/AJYo0hpQtcSAGs2GQOOFwzgtQfAmk22\ngBOWYACvPQDWbAJwiSYAiwBYswnAJZoALAJgxabMg9AA1ihSmtC1pAf48e0rrafQfeccsFrT\nKQCOL8o3UVV6gJum2bx8JvQUuu8ADGAAy/QA//553htunv78RPYUuu/uG3C8YACvPQOPgT9f\nN3vDj3HrcKH7DsAABrDM0EGsn9emXYZjegrddwAGcA7gpvlz/iF8bh+uO3X9y2Utm2bz+nPz\nG943s23tVWKGdaq7Oef7uV1+v56a54ieQvcdgAGcB3jzc/ohwkgc4H0++79BnBhra8/LB/z5\ndNl7jiordN+5Bpz7LBKANYqmp/t5xzIWcNDlxx9+XprN73TFKOD49J9Gaprn7/NFm5trj6fQ\nfecbsFrTOQCOL5qe7s1b8378YddSa15+2lPfm6fDec/73c6fx+b5APDr+bA7vLusmefldbf7\nPfy91ujPYT+1B3i3e2neTifeNs3j+3FpFr/l+Juefi5/5dR+O6yf5+MQzkWzgJvX76FrzafQ\nfXfvgGMFA7g3v5v9GnWG87tpH7L+Hk49NS/7Pw9PwfzZL2L7E7vPI9fXAcDt33vcnRqebwB/\nH5b5w4njwaT3i87jb9mdftPxN98AlsPanIZwLpoF3F/6ARzRZA84egkG8A3gn3Yn+shr/9PT\nkehxpX3Z/Tn8+Odw6ePheNf3WdXZ5tP+77wdSb0fG36fbgDvrn/rZ/d12JM9PcC9rudPv+ff\nLK9/M6z91d6Pf/9UdHuLbm5gm03M7jOAjwFwmaY8wHsR78cf2rX457CUNqdV+efwx2Vd/Pl8\ne+oBPvjd/7226vnSMAp407yI40nNZe1vxG+W1x8Zlii6vUXi500jAuDoJgCXacoEvFf4I1B2\nfxJ/PJ0hXK7x+3RaveUlNwexROfnntTjT6d19DdfrjJ04aXo9haJn9+F38EHzBMpdN8BGMDZ\ngA+r2zzgl+bx/fNHAv7ZHB4bBwD+OqzOxxPfj83mKx/wpej2Fg3cwJQUuu88A872C2CNoslp\nepzf781bf1/1cmHX0a8AvPfbHhA+7UIffxrchX4+76Qff1kzALjzm386Vxkc1rno9hbFOR1N\nofvu7gFHCgZwb7qfH8r2jxZdLhSAv04HqE7/nf3u/97+hz+Hv/x2PBjVBXx4Hnh3fuj6tfsW\nB7F2l+724Nfb4Rp/xG+5OYi16xfd3qLureMxcHpTCcCxSzCABwHvF9Om+3zNLeDXzo6ysHH6\ne833wNNIp3x1OtonhTddwOenkdprvB3PO11laFiXIgCrFwG4dFM+4OMTRfIVE7eA9w+Cm6ev\nIcDt33tqH5D+PHdfyNE+Qfz6e2l73TSbA7v3PuDDS0bag1L7a7y1512vMjSsc9HtLYqEOpZC\n9x2AAcy7kWQArNYE4EJNABbp7kLzPHB6U/ZbGQCsUqSpo4IAWKsp3y+ANYo0dVQQdqG1msoA\njhQM4LUHwFpNAC7VBGARnkbSagJwqSYAiwBYqwnApZoALMIutFYTgEs1ZQD+30SUJBQOgJWa\nFPwCWKNocpquH/Dv62PTPL5FfzJHofsOwHGCAdzN6gH/nN7Uv+GbGSKbAFysCcAiPcBP7Wfl\n/cR9JvQhhe47AAMYwDLDb+j/5Sh0bBOAizUBWKQH9bk5PvplBY5tAnCxJgCL9Ffa5+MudKzf\nuwes8FYGAKsUTU7TdQPmzQzpTRp+AwHHCAZwNwAeS6H7DsBxSzCAu1k34JwUuu8ADGBTwL3l\nS4uHXQCs0wTgck2mgHVZ2Kc/0tfBXejN5vpdK5vTCXkegAFcrgnAIr2Rvg4+Bt5c/rj+v3Pe\n3QNWOQgNYI2iyekeAbhp3523O+5V706n/aU3pk3z/dT8/D41nW9xAPBck4rfUMARggHczf/+\n9zCcAcDN6f/ytL/cvhLrrfnc/Z6+x/wUiXUzcN4OwGpNnQwCjlmCAdxNxEGsC1z5f3+5Bfwp\nv9jlmA7g80Pg63n/7GM9TufZ2tTmA7YZV7UJ3YXe1Qr4uflz+FqlrwnApz9YgWUTK3DBpgKP\ngXe1Aj7Ibb8X9UWe28MK4JsmABdsArBIf1ifj+2Xwrx2zgTwXBOACzYBWCRoWOxCzzTp+A0G\nHC4YwN0kAN7Jg1r+kgS4dxDrkEL3HYDjlmAAdzMLuLqEfSbW+VVX8hVYvBLr2gTgkk0AFuEz\nsTSaAFyyCcAifCaWRhOASzYBWITPxNJoAnDJJgCL8JlYCk06b2UAsErR5DRdPWA+EyulSclv\nMOBwwQDuZt2A+UidxCYAF23KALy+AFihCcBFm1iBRfhIHYUmABdtArAIgBWaAFy0CcAiw6/E\neuXbCSOatA5CA1ijaHKarh4wr8RKaNLyGw44WDCAu1k94JfzK7Fehq8+mkL3HYABDGCZ4Vdi\nxX+CV6H7DsAABrAMgPObAFy2CcAi7ELnNwG4bBOARTiIlR21g9AA1iianKazgCNexuTjIzp4\nGik7an4jAIcKBnA384CHWQyf5SG8kCM7AC7cBGCR/hv6Yx/7nlPovgMwgEsBPu9KH/8v//T0\nvUk3342U2FPovgMwgPMAb4dzC3joY2U735Pk43uTer/5++k19vDVMYXuOwADuNBBrP7nQY9+\nX9Kynxl98zwwbyeMjN5BaABrFE1O09Bd6N0V53l5Pe8xNwC2u/OWKNoa3rhxwIGCAdxNCuDr\nHnNnNfYJODmF7jsAA7g84O5j4B2Aze68JYoAXLqpLODBg1hD/18mnV/9/dQ0LxzEiguASzeV\nAdx7GmnXfRpp13mqyQng7+Oj3++UnkL3HYABzGuhZSTgl8O3ir5Ev4+hTaH7DsAXwEGCAdzN\nugG3ewO/zWbsulMpdN/5A7y1vHETgMOWYAB3cweAE19VUui+AzCAASwD4MwAuHgTgEUAnBkA\nF28CsAiAMwPg4k0AFukC5qtVorMc4BDBAO4GwGMpdN+5A7w1vXFTgIOWYACvPbyUMi8ALt8E\nYBEA5wXA5ZsALALgvAC4fBOARQCcFwCXbwKwCIDzsiDgEMEAXnsAnJWtWtPfoSYAxxcpTeha\nAuCsAHiBJgCLADgrk4A/8gPg6CKlCV1LAJyVacDzwOb8ATi6SGlC1xIAZ2VZwPOCAbz2ADgn\n20UBByzBAF57AJyT7WQTgG2aACwC4JwAeIkmAIsAOCcAXqIJwCIAzgmAl2gCsAiAc7I04FnB\nAF57AJyR7cKA55dgAK89AM7IdroJwDZNABYBcEYAvEgTgEUAnBEAL9IEYJEFv1et/mynLy4A\neE7wR5kNQRYLK3BGFl+BZ5dgVuC1B8Dp2QJ4kSYAiwA4PduZJgDbNAFYBMDpAfAyTQAWAXB6\nALxME4BFAJweD4BnBAN47QFwehwAnluCAbz2ADg527kmANs0AVgEwMkB8EJNABYBcHIAvFAT\ngEUAnBwfgKcFA3jtAXByXACeWYIBvPYAODVbAC/UBGARAKfm7BfApZsALALg1AB4qSYAiwA4\nNQBeqgnAIgBOjRPA04IBvPYAODUAXqoJwCIATszFL4BLNwFYBMCJAfBiTQAWAXBiALxYE4BF\nAJwYN4AnBQN47QFwYgC8WBOARQCcli2AF2sCsAiA03L1C+DSTQAWAXBaALxcE4BFAJwWT4An\nBAN47bkvwA/7nP7X/ykyjgBPLcEAXnvuBvBDq3W0KFLxFsDLNQFY5C4AD62yN0VRK7HwC+DS\nTQAWuQfAgzBHigIRA3jBJgCLrBzw+Ko6VhS2EPsCPC4YwGvPmgFPWpwoqg3wxBIM4LVnzYAz\niuaWYQAv2ARgkZUCnt8Pni2aqpB+AVy6CcAiqwQc8jA2a0gAXrIJwCIrBBx2HDnrXwJngMcF\nA3jtWR3g0Gdzs/bFAbxkE4BFVgY4/MUYwUMaKNwCeMkmAIusDLBNUf9fhY5fAJduArDIigDH\nvZo5aki9hR3AizYBWGQ1gGPfUxT9ki5xwh3gUcEAXnvWAjj6LYE5QwLwok0AFlkH4IQ39MYP\n6bLGbwG8aBOARdYAOOkd+Wlva2p/UdcvgEs3AVhkBYBT+KYOyS3gEcEAXnvqB5zmN+ejPRwC\nHluCAbz2VA44De9AUXgeHgC8bBOAReoGnLj63hbFZLvt/lYAF24CsEjNgJMOXg0VxWX7F8CL\nNgFYpGbACxW1e9Djn295TVHAw4IBvPZUCzhj8e0WRecI+Lr6uwA8sgQDeO2pFXDO3nOnKD7b\n3ggAXLgJwCJ1As7mmzGk3uuwAFy8CcAidQJeskj6fZhsKgt4UDCA154KAeevvqeitHQW4MOu\ngA/Aw0swgNeeMMCbffo/b+SZ5QBnPXcki1L/YncPej8aABduArBIEODN5Q/x86ZzlUL3ncaj\n32OUAO9HBODCTQAWqQ3w4o+B+37ZhS7eBGCRaMDnM7p+ywCeesAZG0XAY3v1ZQEPCgbw2pMO\n+PoQ+J991Ac2kIeHEr9lJtuhM4dH5gGw7cYgiycJ8Ob2PPt/fI/rnMcV+O/IQ3MPgBNvZffW\n6YQV2CLJgLs/3M1j4JuXcbh5DAzg06X3lRTAQz8Z33fzLzyOTSrg0abbR8KlAQ8IBvDakwB4\nM3CeMeCAdw7ERh3w7X50YcBDSzCA1554wOLpJLFXbXnfybXNM+C+YADbNAFYJO6VWOejzxtx\n3jGF7jvfgHsBsE0TgEXcvxa6t1u6MOABv90mubdQHPCtYACvPd4B9x9XegcsR1wa8MASDOC1\nxzfg2yO7/gFfBQPYpgnAIr4BGzbZAb4EwDZNABYBcEwGXsYx2HRchMsDvhEM4LXHL+Dhdwgs\nDDi46TD64oBvl2AArz1uAY+8w6cWwIfxA9imCcAiXgGPvXG/GsB/WYGtmgAs4hPw+MduLAp4\n8CHweNPDQ3HAfcEAXns8Ap762KtlAcc1fTzkEo5ewwEM4LQo3neTH3tVFeD//ssUDOD4IqUJ\nXUv8AZ7+2LrKAGcGwPFFShO6lngDPPehk0sCHn4IPAc4ZxGOB9wTDOC1xxfg+Q99XhRwZNOJ\nX8Yj4fg1HMBKE7qWuAIc8KHPFQLOWIQBHF+kNKFriSfAIR/aXiXg5AA4vkhpQtcST4Dz7zzT\nopGHwGGA0xbhBMBdwQBee7wADv3KlCUBxzZ1+aUQTljDAXxf8QH45uDVh8Nso/9Gl1aC4HzA\nDhI6CYIDYBEXgG+X3yUePs4VDH+FtuUQFAAvMYReQeAkCA+ARRwAdvClBiEFw9+gHTeE2GeU\nUgB3hulhO4a5iwiARZYH7OFbSUIKov0ODiFOcMqNAPBdZWnAIy/d8DDxeqd1AMcJBnDcbGov\nva8sDHjs4LOHidc7rQT4kGDDSYDlQD1sxzB3EQGwyKKAx5868jDxuifjHwJPAA4VnHQjAHxP\nWRDw1AufPUy87sl4vxNDCD2YBeDg2XS99L6yGODp9y14mHjdk6qAQ/eiARw4m+Sl95WlD2KN\nxMPE655UBvxfkOE0wGKoHrZj8iRIm01KE7qWLAN49nWTHiZe51TCQ+BZwPP70Wk3AsB3lCUA\nz7/r1yFgiyHMEgbw/Gy6ufS+sgDgkPcteJh4nVMmgI0axN6Ch+0Y5i4iABbhMXBQQcoedNgQ\nJhfhxBsB4PtJYcCh7xr0MPHkiRS/gUOY2o8G8ORsGr70vlIUcMCD31M8TDx5whDw1CIM4InZ\nNHbpfaUk4GC+dwZYveG6w+9hO4a5iwiARcoBDtf71x3gpIfAUUMYXoVTbwSA7yalAIfvPbfx\nMPHEz0l+44Yw+FAYwMOzafLS+0oZwJF87xHw4CKcDPg8YA/bMcxdRAAswtNIAQVpe9ApQ+gZ\nTr4RAL6XADigIM1vEuDuMgzg+NmkNKFrCYADCsoB7hFOB7zNLMgfwqUgeRKkzSalCV1LABxQ\nUBKwUgOA7yQAni9I9Js+hMsqDOD42aQ0oWsJgOcLigO+EM4AvM0dwikA9h0AzxcsAPh0PDqj\nAcD3EQDPFiQ+iZQ/hIcHAMfPJqUJXUsAPFuQ6jd/CBlfDX76Z8fDdkyeBGmzSWlC1xIAzxYs\nB/i/jxzBSkPILkieBGmzSWlC1xIAzxUk70Fr6UlEDOC7CIDnCpL9qulJ25Nu/+HxsB2TJ0Ha\nbFKa0LUEwHMFywNOXIQBfA8B8ExB+h60sp5owwC+hwB4piDdrzbg2D3pwz89HrZj8iRIm01K\nE7qWAHimwA3g/6L3pAF8B9ECrBwPE6/9M2MP2kZPBOL92D1sx6Xn0rrDCjxdkOHXSE/EnrQP\nwMmTIG02KU3oWgLgyYKcBdhMT7BgAK8/AJ4syPFrqicI8XbrYTsmT4K02aQ0oWsJgCcL3AIO\nOygN4NUHwFMFWXvQ5vuvAB689L4C4KmCLL8lHoDOLcN5NyBoCLMFyZMgbTYpTehaAuCpAveA\n53alAbz2AHiiIG8PuuQh4DHDmbcgZgijBcmTIG02KU3oWgLgiYLM2V8S8Mg6/AHglQfA4wW5\ny1fZGzEo+CN7CQaw7wB4vGDxuR/dcLMOO7gRADYNgMcLFp/7CQ29fWkHNwLApgHwaMHye5/p\nDWfEDh4HANg0AB4tWP74Twbg6yfDA3jVAfBYHDwDk9XQCnbwXBiATQPgsTh4DUR2Q7sQA3jN\nAfBI8hdgB4D3BQ8PC7+gG8CmAfBIPLwPQO0tVdGfp6U5hORJkDablCZ0LQHwcFy8lVbzUwlS\nCQPYdwA8HBcfZqH+uUAJKzGAfQfAg/HxeXAWn8wXSxjAvgPgwfj4RFazz8aNWIkB7DsAHoqT\nz0TXA3x7IDqUMIB9B8BDcfKtJGqAx54LfphfigHsOwAeiJcv9tMDPPFqjlbwuGIA+w6AB+Ll\nq3U1AU+/muNhdC0GsO8A+DbbFQIOfEHlrWIA+w6Ab7PNLcgfglaDABz4gsr+Sgxg3wHwTbar\nBBz3nobrLjWAfQfAN9nmFuQPQa1BAo5+T0N7dCv5RdSXISRPgrTZpDShawmA+9muFHDi+6s+\nMt4IcSxIngRps0lpQtcSAPdymeYO+OkCTntj8KnhYfw49VxB8iRIm01KE7qWALiXyyx3wE8b\ncIrg/hCiHQPYNADuZrtiwEmCh4bwEPPoGMCmAXAnYoo74KcNOEXwxBBOa/HcazGTJ0HabFKa\n0LUEwJ2ICe6AnzrghIfB80OYeXQMYNMAWGa7esDRgsOHMHKcC8CmAbBIZ3o74KcPOF5wwhC6\nkgFsGgBf053cDvgZAI4WnDOEk+TUSZA2m5QmdC0B8DXdqe2AnwXgWMEqL6XsSs4NgEUAfMn2\nPgBHHsjSfS30VXIGZwCLAPic/srkgODgBIgAAAaNSURBVJ8N4DjBZm9muEqOtQxgEQCfcrNn\n6YCfFeAYwSXejXSVHIQZwCIAPuZ2UjvgZwR4f2PDCS/wdsK5/WwAiwC4zcCMdsDPCnDMIrzs\n+4HlfvZldQawCIAPGZrPDvjZAQ4X7OwN/QfGABYB8H8js9kBP0PAwYKdAZ4vUprQtQTAY48I\nHfCzBBwqGMC+A+CxieyAnyngwENZAPaduwc8Oosd8LMFHLYIA9h37hzwxCrkgJ814JBFGMC+\nc9+Ap+avA37mgAMIA9h37hnw9OR1wK8A4Nn9aAD7zv0Cnlt7HPArAviwISx3RABsmnsFPP/o\nrzS/fw/RHkNggeWhAACb5i4BT685AQX5Q+jn338/DukRLgV4YpMA2HfuD3CQ3qmC/CHc5si3\nJaw6hrjXo5m8ohTApgkDvNmn/7M8rxbA21C9YwX5QxjM1W9XcFHA/522j/LnGgDYNEGAN5c/\nrj/L83bFADdNEzpveqcH5uZ078fkpSlDGP+tIYC7g+gNaWyE8fyaprulIhqGBwFg09QFuGmC\nLZ0atpfMTL2b3o/JSyOGMP9bpd+O4I+BKw8NaXSE0YDPTduwDTf0V/tDSJ4EabNJGYj3VAW4\nabqTZBuQ4KnXn3wfk5eGZEbPtTcAcHcQvSGNjzAW8G3TR8hGPmznsUEA2DQagP/ZZ7ogcBJE\n50Mnp7mXdKnGb+0Bnh9Eb0h6I0xusrqHh5I41dcZpyvwcE7TS6VLDmm6d+63pt64a28P8Pwg\nekO6HWH+mKKbZu8cVmCLVAX4r57f7pCme60m5qU3AHBvEL0h3Ywwf0zxTXN3DoAtUhfgv2p+\ne0Oa7rWamJfef2f99gfRG1J/hApjim+auXMAbJHKAJeaBYWbQgBHxdWt024CsAiAXTT9q+vX\n2a1TbgKwSNwrsTbyZ8tXYtk3ORvS5bXQGmXubp1uE4BFnL4W2r7J3ZDadyPpVDm8dZpNABYB\nME21NQFYBMA01dYEYBEA01RbE4BFAExTbU0AFgEwTbU1AVgEwDTV1gRgEQDTVFsTgEUATFNt\nTQAWATBNtTUBWATANNXWBGARANNUWxOARQBMU21NABYBME21NQFYBMA01dYEYBEA01RbE4BF\nAExTbU0AFtECPJ2Zj41eIg6HxJjC4nBIywXAjsKYguJwSMsFwI7CmILicEjLBcCOwpiC4nBI\nywXAjsKYguJwSMulDGBCiEkATEjFATAhFQfAhFQcABNScUoB3sxfpWS6X+zkJf6G5HA7ORzS\nkikE2Nk27321opM420g7l9vJ4ZAWTRnAG2eb3OUs8LaRdi63k8MhLZr73IV2Ogv8jegQh6Ny\nOKSlAmBH8TeiQxyOyuGQlgqAHcXfiHYOB8VBLBkAO4q/Ee18DsrjmBaKMeDLMX9fmxzAwfE4\nJp+DWiaswI7ib0T+huT0rlssAHYURjQfp3fdYrlPwE5fzuNuSJuNvw3lb0SLhtdCE1JxAExI\nxQEwIRUHwIRUHAATUnEATEjFATAhFQfAhFQcABNScQBMSMUBsGGaNpvXn+7Z7zcvBfxtHk8/\nPTbiyg33DpkJU8QwzTmfvbNvrvl0cvvTPE1fkZBOmCKGOQL8eWk2v7dnd/KneWv//9b8mb4i\nIZ0wRQxzBvjS8vx6PuxOH9dlcbLNeR/6sfkVlxyud+xo//x9aZqX9p+Ct03z+F72thCfAbBh\nzoC/DzvGn8e96dcT4MvJY57bfeifzhX7gDeHsw/QX9srIJgA2DSXXeDDD4+HvePvi8nryTaf\n7SLd7kF3r3gF/HYw/Xpw2+y1fzW8K5YA2DQdwPv19fPt6WrycvJ0lcPS2u5Bd694Bfx4/On5\nsBS/dA+LkbsNgA3TBfx03DU+n305eczLflX9OejsXfEK+HxIe79c73emH3/6v43cYwBsmDPP\nrwPMl+bx/fPnYvJ68pjDPvRb+3xT94pDgPd72I/N5qv8DSLuAmDDnHk+Hx+47g5Hm88mryfP\nV3487yR3rtie+rnuQp/zznNMZAdg01yfB25PfO1+nwTg88lTXprP5mV3c8VN8+d06vVwEOvP\n4Tj1Zn+Nbw5ikR2ATXN5JdZhb/e1uT603ciTp3yeX7Elr3g69Xb46bd9Gqn5Pl/jbambRRwF\nwIY5Qnx8Pb4O66Vpnr4OEt/bxfNy8pzN9XUfp0vaM143e6rHHen2gt3xvA1+yQ7AhFQdABNS\ncQBMSMUBMCEVB8CEVBwAE1JxAExIxQEwIRUHwIRUnP8DXf3SCrQ4yb0AAAAASUVORK5CYII=",
      "text/plain": [
       "plot without title"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "options(repr.plot.width = 8, repr.plot.height = 5)\n",
    "\n",
    "ggplot(data=data.frame(x=x, y=y)) + \n",
    " geom_histogram(aes(x=y, y=2*..count../sum(..count..)), bins = 10, col=\"gray\") + \n",
    "  geom_line(data=data.frame(u=x1, v=y1), aes(x=u,y=v, color=\"cyan\"), linetype=\"longdash\") + \n",
    "   geom_point(x=y, y=0, colour=\"black\", fill=\"black\") +\n",
    "    geom_point(x=ybar, y=0, shape=21, size=4, color=\"black\", fill=\"gray\", alpha=0.5) + \n",
    "     geom_line(data=data.frame(u=x1, v=y1), aes(x=u, y=0.0140*exp(log_gdist(mu=u, n=n, ybar=ybar)), color=\"red\")) +\n",
    "      scale_colour_manual(name=\"Normalized Distributions\",values =c('cyan'='cyan','red'='red'), labels = c('Prior','Posterior'))+\n",
    "       xlim(-1.0,3.0)+\n",
    "        ylab(\"Probability\")+\n",
    "         xlab(\"Data Values\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}