MH.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Distributions, Random, StatPlots, KernelDensity\n",
    "pyplot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function propose(prev::Number)\n",
    "    rand(Normal(prev,1))\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function MH(lpdf::Function, iter=1000, init=0, propose=propose)\n",
    "    trace = zeros(iter)\n",
    "    point = init\n",
    "    for i=1:iter\n",
    "        trace[i] = point\n",
    "        propose_point = propose(point)\n",
    "        ratio = lpdf(propose_point) - lpdf(point)\n",
    "        if log(rand()) < ratio\n",
    "            point = propose_point\n",
    "        end\n",
    "    end\n",
    "    trace\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function likelihood(data)\n",
    "    n_head = count(x->x==1, data)\n",
    "    n = length(data)\n",
    "    function f(p)\n",
    "        if p>1 || p<0\n",
    "            -Inf\n",
    "        else\n",
    "            n_head*log(p) + (n-n_head)*log(1-p)\n",
    "        end\n",
    "    end\n",
    "    f\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = rand(Bernoulli(0.4),500);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "anim = @animate for k=1:length(data)\n",
    "    observed = data[1:k]\n",
    "    t = MH(likelihood(observed),10000,0.5)\n",
    "    d = kde(t)\n",
    "    plot(d,xlim=[0,0.8],ylim=[0,17])\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gif(anim, \"./anims.gif\", fps = 15);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.0.3 1.0.3",
   "language": "julia",
   "name": "julia-1.0.3-1.0"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.0.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}