MambaLine.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "using Mamba, Gadfly"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Model and user-defined sampler specification"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "An object of type \"Mamba.Sampler{Dict{Any,Any}}\"\n",
       "Sampling Block Nodes:\n",
       "[:s2]\n",
       "\n",
       "AST(:($(Expr(:lambda, Any[:(model::Model),:(block::Integer)], Any[Any[Any[:model,:Any,18],Any[:block,:Any,18],Any[:f,:Any,18]],Any[],0,Any[]], :(begin \n",
       "        model = (top(typeassert))(model,Mamba.Model)\n",
       "        block = (top(typeassert))(block,Mamba.Integer)\n",
       "        f = (anonymous function)\n",
       "        return f((Mamba.getindex)(model,:mu),(Mamba.getindex)(model,:s2),(Mamba.getindex)(model,:y))\n",
       "    end)))))\n"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model = Model(\n",
    "\n",
    "  y = Stochastic(1,\n",
    "    (mu, s2) ->  MvNormal(mu, sqrt(s2)),\n",
    "    false\n",
    "  ),\n",
    "\n",
    "  mu = Logical(1,\n",
    "    (xmat, beta) -> xmat * beta,\n",
    "    false\n",
    "  ),\n",
    "\n",
    "  beta = Stochastic(1,\n",
    "    () -> MvNormal(2, sqrt(1000))\n",
    "  ),\n",
    "\n",
    "  s2 = Stochastic(\n",
    "    () -> InverseGamma(0.001, 0.001)\n",
    "  )\n",
    "\n",
    ")\n",
    "\n",
    "Gibbs_beta = Sampler([:beta],\n",
    "  (beta, s2, xmat, y) ->\n",
    "    begin\n",
    "      beta_mean = mean(beta.distr)\n",
    "      beta_invcov = invcov(beta.distr)\n",
    "      Sigma = inv(xmat' * xmat / s2 + beta_invcov)\n",
    "      mu = Sigma * (xmat' * y / s2 + beta_invcov * beta_mean)\n",
    "      rand(MvNormal(mu, Sigma))\n",
    "    end\n",
    ")\n",
    "\n",
    "Gibbs_s2 = Sampler([:s2],\n",
    "  (mu, s2, y) ->\n",
    "    begin\n",
    "      a = length(y) / 2.0 + shape(s2.distr)\n",
    "      b = sumabs2(y - mu) / 2.0 + scale(s2.distr)\n",
    "      rand(InverseGamma(a, b))\n",
    "    end\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Object of type \"Mamba.Model\"\n",
       "-------------------------------------------------------------------------------\n",
       "y:\n",
       "An unmonitored node of type \"0-element Mamba.ArrayStochastic{1}\"\n",
       "Float64[]\n",
       "-------------------------------------------------------------------------------\n",
       "s2:\n",
       "A monitored node of type \"Mamba.ScalarStochastic\"\n",
       "NaN\n",
       "-------------------------------------------------------------------------------\n",
       "beta:\n",
       "A monitored node of type \"0-element Mamba.ArrayStochastic{1}\"\n",
       "Float64[]\n",
       "-------------------------------------------------------------------------------\n",
       "mu:\n",
       "An unmonitored node of type \"0-element Mamba.ArrayLogical{1}\"\n",
       "Float64[]\n"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "scheme1 = [NUTS(:beta),\n",
    "           Slice(:s2, 3.0)]\n",
    "\n",
    "## No-U-Turn Sampling Scheme\n",
    "scheme2 = [NUTS([:beta, :s2])]\n",
    "\n",
    "## User-Defined Sampling Scheme\n",
    "scheme3 = [Gibbs_beta, Gibbs_s2]\n",
    "\n",
    "\n",
    "## Sampling Scheme Assignment\n",
    "setsamplers!(model, scheme1)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "digraph MambaModel {\n",
      "\t\"y\" [shape=\"ellipse\", fillcolor=\"gray85\", style=\"filled\"];\n",
      "\t\"mu\" [shape=\"diamond\", fillcolor=\"gray85\", style=\"filled\"];\n",
      "\t\t\"mu\" -> \"y\";\n",
      "\t\"s2\" [shape=\"ellipse\"];\n",
      "\t\t\"s2\" -> \"y\";\n",
      "\t\"beta\" [shape=\"ellipse\"];\n",
      "\t\t\"beta\" -> \"mu\";\n",
      "\t\"xmat\" [shape=\"box\", fillcolor=\"gray85\", style=\"filled\"];\n",
      "\t\t\"xmat\" -> \"mu\";\n",
      "}\n"
     ]
    }
   ],
   "source": [
    "draw(model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5x2 Array{Float64,2}:\n",
       " 1.0  1.0\n",
       " 1.0  2.0\n",
       " 1.0  3.0\n",
       " 1.0  4.0\n",
       " 1.0  5.0"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "line = Dict{Symbol, Any}(\n",
    "  :x => [1, 2, 3, 4, 5],\n",
    "  :y => [1, 3, 3, 3, 5]\n",
    ")\n",
    "line[:xmat] = [ones(5) line[:x]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{Dict{Symbol,Any},1}:\n",
       " Dict{Symbol,Any}(:y=>[1,3,3,3,5],:s2=>1.1830434911443408,:beta=>[1.1902678809862768,2.04817970778924])   \n",
       " Dict{Symbol,Any}(:y=>[1,3,3,3,5],:s2=>0.3057490088250573,:beta=>[1.142650902867199,0.45941562040708034]) \n",
       " Dict{Symbol,Any}(:y=>[1,3,3,3,5],:s2=>0.6108061237629374,:beta=>[-0.396679079295223,-0.6647125451916877])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "srand(123)\n",
    "\n",
    "\n",
    "## Initial Values\n",
    "inits = [\n",
    "  Dict{Symbol, Any}(\n",
    "    :y => line[:y],\n",
    "    :beta => rand(Normal(0, 1), 2),\n",
    "    :s2 => rand(Gamma(1, 1))\n",
    "    )\n",
    "  for i in 1:3\n",
    "]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MCMC Simulation of 10000 Iterations x 3 Chains...\n",
      "\n",
      "Chain 1:   0% [0:31:50 of 0:31:52 remaining]\n",
      "Chain 1:  10% [0:00:24 of 0:00:27 remaining]\n",
      "Chain 1:  20% [0:00:13 of 0:00:16 remaining]\n",
      "Chain 1:  30% [0:00:09 of 0:00:12 remaining]\n",
      "Chain 1:  40% [0:00:06 of 0:00:10 remaining]\n",
      "Chain 1:  50% [0:00:05 of 0:00:09 remaining]\n",
      "Chain 1:  60% [0:00:03 of 0:00:08 remaining]\n",
      "Chain 1:  70% [0:00:02 of 0:00:08 remaining]\n",
      "Chain 1:  80% [0:00:01 of 0:00:07 remaining]\n",
      "Chain 1:  90% [0:00:01 of 0:00:07 remaining]\n",
      "Chain 1: 100% [0:00:00 of 0:00:07 remaining]\n",
      "\n",
      "Chain 2:   0% [0:00:05 of 0:00:05 remaining]\n",
      "Chain 2:  10% [0:00:04 of 0:00:04 remaining]\n",
      "Chain 2:  20% [0:00:03 of 0:00:04 remaining]\n",
      "Chain 2:  30% [0:00:03 of 0:00:04 remaining]\n",
      "Chain 2:  40% [0:00:02 of 0:00:04 remaining]\n",
      "Chain 2:  50% [0:00:02 of 0:00:04 remaining]\n",
      "Chain 2:  60% [0:00:02 of 0:00:04 remaining]\n",
      "Chain 2:  70% [0:00:01 of 0:00:04 remaining]\n",
      "Chain 2:  80% [0:00:01 of 0:00:04 remaining]\n",
      "Chain 2:  90% [0:00:00 of 0:00:05 remaining]\n",
      "Chain 2: 100% [0:00:00 of 0:00:05 remaining]\n",
      "\n",
      "Chain 3:   0% [0:00:04 of 0:00:04 remaining]\n",
      "Chain 3:  10% [0:00:03 of 0:00:04 remaining]\n",
      "Chain 3:  20% [0:00:03 of 0:00:04 remaining]\n",
      "Chain 3:  30% [0:00:02 of 0:00:03 remaining]\n",
      "Chain 3:  40% [0:00:02 of 0:00:03 remaining]\n",
      "Chain 3:  50% [0:00:02 of 0:00:03 remaining]\n",
      "Chain 3:  60% [0:00:01 of 0:00:03 remaining]\n",
      "Chain 3:  70% [0:00:01 of 0:00:03 remaining]\n",
      "Chain 3:  80% [0:00:01 of 0:00:03 remaining]\n",
      "Chain 3:  90% [0:00:00 of 0:00:03 remaining]\n",
      "Chain 3: 100% [0:00:00 of 0:00:03 remaining]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "sim1 = mcmc(model, line, inits, 10000, burnin=250, thin=2, chains=3);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iterations = 252:10000\n",
      "Thinning interval = 2\n",
      "Chains = 1,2,3\n",
      "Samples per chain = 4875\n",
      "\n",
      "Geweke Diagnostic:\n",
      "First Window Fraction = 0.1\n",
      "Second Window Fraction = 0.5\n",
      "\n",
      "        Z-score p-value\n",
      "beta[1]   1.237  0.2162\n",
      "beta[2]  -1.568  0.1168\n",
      "     s2   1.710  0.0872\n",
      "\n",
      "        Z-score p-value\n",
      "beta[1]  -1.457  0.1452\n",
      "beta[2]   1.752  0.0797\n",
      "     s2  -1.428  0.1534\n",
      "\n",
      "        Z-score p-value\n",
      "beta[1]   0.550  0.5824\n",
      "beta[2]  -0.440  0.6597\n",
      "     s2   0.583  0.5596\n",
      "\n"
     ]
    }
   ],
   "source": [
    "gewekediag(sim1) |> showall"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iterations = 252:10000\n",
      "Thinning interval = 2\n",
      "Chains = 1,2,3\n",
      "Samples per chain = 4875\n",
      "\n",
      "Empirical Posterior Estimates:\n",
      "           Mean       SD       Naive SE       MCSE       ESS   \n",
      "beta[1] 0.5971183 1.14894446 0.0095006014 0.016925598 4607.9743\n",
      "beta[2] 0.8017036 0.34632566 0.0028637608 0.004793345 4875.0000\n",
      "     s2 1.2203777 2.00876760 0.0166104638 0.101798287  389.3843\n",
      "\n",
      "Quantiles:\n",
      "            2.5%       25.0%       50.0%     75.0%     97.5%  \n",
      "beta[1] -1.74343373 0.026573102 0.59122696 1.1878720 2.8308472\n",
      "beta[2]  0.12168742 0.628297573 0.80357822 0.9719441 1.5051573\n",
      "     s2  0.17091385 0.383671702 0.65371989 1.2206381 6.0313970\n",
      "\n"
     ]
    }
   ],
   "source": [
    "describe(sim1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[Gadfly.Plot([Gadfly.Layer(nothing,Dict{Any,Any}(),Gadfly.StatisticElement[],Gadfly.Geom.LineGeometry(Gadfly.Stat.Identity(),false,2,symbol(\"\")),nothing,0)],nothing,Data(\n",
      "  x=["
     ]
    }
   ],
   "source": [
    "show(plot(sim1))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.4.2",
   "language": "julia",
   "name": "julia-0.4"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.4.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}