March2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Topics for March 2, 2016\n",
    "\n",
    "My goal is to provide fast methods to use Markov-Chain Monte Carlo (MCMC) methods on linear mixed-effects models, with possible extensions to generalized linear mixed models.\n",
    "\n",
    "There are several MCMC frameworks for [Julia](http://julialang.org).  An interface to [Stan](http://mc-stan.org) is available and there are two native Julia implementations; [Mamba](https://github.com/brian-j-smith/Mamba.jl) and [Lora](https://github.com/JuliaStats/Lora.jl).  I prefer `Mamba` because of its flexibility.  The problem with Stan, BUGS and JAGS is that each of them reinvents all the data structures, input/output, data manipulation, distribution definitions, etc. in its own environment.  They also define a Domain Specific Language (DSL) for which they must provide parsers, interpreters, run-time environments, etc.\n",
    "\n",
    "A native implementation like Mamba can use all of the facilities of Julia and its packages.\n",
    "\n",
    "## Linear predictor\n",
    "\n",
    "Whenever you have a linear predictor (i.e. an $\\bf X\\beta$ type of expression) in a model there is a good chance that you can write out the full conditional distribution of $\\beta$ or obtain a good approximation to it.  If you can write out the conditional distribution you can use a multivariate Gibbs sampler to obtain a vector-valued sample from the condtional.  This helps to avoid one of the underlying problems of MCMC methods which is successive sampling from conditionals of correlated parameters.  Consider a case where you might have hundreds or thousands of random effects and dozens of fixed effects.  You don't want to sample sequentially in those cases when you can sample from the distribution of the entire vector in one step.\n",
    "\n",
    "## Multivariate normal conditionals\n",
    "\n",
    "In most cases the conditional distribution of the coefficients (i.e. both random and fixed effects) is a multivariate normal with known mean and covariance matrix.  It is worthwhile examining the representation in Julia of this distribution.  Not surprisingly the representation involved the mean and covariance but the form of the covariance is encoded in the type.  For example, a common prior distribution for coefficients is a zero-mean multivariate normal with a covariance that is a large multiple of the identity - a diffuse prior."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: readbytes is deprecated, use read instead.\n",
      " [inlined code] from ./error.jl:26\n",
      " in depwarn(::ASCIIString, ::Symbol) at ./deprecated.jl:64\n",
      " in readbytes(::Base.PipeEndpoint, ::Vararg{Any}) at ./deprecated.jl:30\n",
      " in send_stream(::Base.PipeEndpoint, ::ASCIIString) at /home/bates/.julia/v0.5/IJulia/src/stdio.jl:25\n",
      " in "
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Julia Version 0.5.0-dev+3149\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "watch_stream(::Base.PipeEndpoint, ::ASCIIString) at /home/bates/.julia/v0.5/IJulia/src/stdio.jl:41\n",
      " in (::IJulia.##4#8)() at ./task.jl:431\n",
      "while loading In[1], in expression starting on line 1\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Commit 3c72273* (2016-03-14 16:15 UTC)\n",
      "Platform Info:\n",
      "  System: Linux (x86_64-linux-gnu)\n",
      "  CPU: Intel(R) Core(TM) i5-3570 CPU @ 3.40GHz\n",
      "  WORD_SIZE: 64\n",
      "           Ubuntu 15.10\n",
      "  uname: Linux 4.2.0-34-generic #39-Ubuntu SMP Thu Mar 10 22:13:01 UTC 2016 x86_64 x86_64\n",
      "Memory: 15.561424255371094 GB (10432.33203125 MB free)\n",
      "Uptime: 3851.0 sec\n",
      "Load Avg:  0.65771484375  0.4658203125  0.5283203125\n",
      "Intel(R) Core(TM) i5-3570 CPU @ 3.40GHz: \n",
      "       speed         user         nice          sys         idle          irq\n",
      "#1  3798 MHz      21737 s        540 s       3672 s     341064 s          0 s\n",
      "#2  3698 MHz      33490 s        990 s       3509 s     329140 s          0 s\n",
      "#3  3689 MHz      19396 s        465 s       3304 s     350182 s          0 s\n",
      "#4  3692 MHz      20761 s        454 s       3654 s     348260 s          0 s\n",
      "\n",
      "  BLAS: libopenblas (NO_LAPACKE DYNAMIC_ARCH NO_AFFINITY Sandybridge)\n",
      "  LAPACK: liblapack\n",
      "  LIBM: libopenlibm\n",
      "  LLVM: libLLVM-3.7.1 (ORCJIT, ivybridge)\n",
      "Environment:\n",
      "  XDG_SESSION_PATH = /org/freedesktop/DisplayManager/Session0\n",
      "  MANDATORY_PATH = /usr/share/gconf/ubuntu.mandatory.path\n",
      "  DEFAULTS_PATH = /usr/share/gconf/ubuntu.default.path\n",
      "  COMPIZ_BIN_PATH = /usr/bin/\n",
      "  NODE_PATH = /usr/lib/nodejs:/usr/lib/node_modules:/usr/share/javascript\n",
      "  XDG_SEAT_PATH = /org/freedesktop/DisplayManager/Seat0\n",
      "  PATH = /home/bates/bin:/usr/local/sbin:/usr/local/bin:/usr/sbin:/usr/bin:/sbin:/bin:/usr/games:/usr/local/games\n",
      "  HOME = /home/bates\n",
      "  TERM = screen-256color-bce\n",
      "\n",
      "Package Directory: /home/bates/.julia/v0.5\n",
      "12 required packages:\n",
      " - FileIO                        0.0.4\n",
      " - GLM                           0.5.0\n",
      " - Gadfly                        0.4.2\n",
      " - HDF5                          0.5.8\n",
      " - IJulia                        1.1.8\n",
      " - Mamba                         0.9.0\n",
      " - MixedModels                   0.4.5+             master\n",
      " - PDMats                        0.4.0\n",
      " - ProfileView                   0.1.1\n",
      " - RCall                         0.3.2+             master\n",
      " - VideoIO                       0.0.14\n",
      " - ZipFile                       0.2.6\n",
      "59 additional packages:\n",
      " - ArrayViews                    0.6.4\n",
      " - BinDeps                       0.3.21\n",
      " - Blosc                         0.1.4\n",
      " - Cairo                         0.2.31\n",
      " - Calculus                      0.1.14\n",
      " - Codecs                        0.1.5\n",
      " - ColorTypes                    0.2.1\n",
      " - ColorVectorSpace              0.1.1\n",
      " - Colors                        0.6.3\n",
      " - Compat                        0.7.12\n",
      " - Compose                       0.4.2\n",
      " - Conda                         0.1.9\n",
      " - Contour                       0.1.0\n",
      " - Cxx                           0.0.0-             master (unregistered)\n",
      " - DataArrays                    0.2.20\n",
      " - DataFrames                    0.6.10\n",
      " - DataStructures                0.4.3\n",
      " - Dates                         0.4.4\n",
      " - Distances                     0.3.0\n",
      " - Distributions                 0.8.10\n",
      " - Docile                        0.5.23\n",
      " - DualNumbers                   0.2.2\n",
      " - FixedPointNumbers             0.1.2\n",
      " - FixedSizeArrays               0.0.10\n",
      " - GZip                          0.2.18\n",
      " - Glob                          1.0.2\n",
      " - Graphics                      0.1.3\n",
      " - Graphs                        0.6.0\n",
      " - Grid                          0.4.0\n",
      " - Gtk                           0.9.3\n",
      " - GtkUtilities                  0.0.8\n",
      " - Hexagons                      0.0.4\n",
      " - ImageView                     0.1.19\n",
      " - Images                        0.5.3\n",
      " - IniFile                       0.2.5\n",
      " - Iterators                     0.1.9\n",
      " - JSON                          0.5.0\n",
      " - KernelDensity                 0.1.2\n",
      " - Loess                         0.0.6\n",
      " - MathProgBase                  0.4.2\n",
      " - Measures                      0.0.2\n",
      " - NLopt                         0.3.1\n",
      " - NaNMath                       0.2.1\n",
      " - Nettle                        0.2.2\n",
      " - Optim                         0.4.4\n",
      " - Reexport                      0.0.3\n",
      " - SHA                           0.1.2\n",
      " - SIUnits                       0.0.6\n",
      " - Showoff                       0.0.6\n",
      " - SortingAlgorithms             0.0.6\n",
      " - StatsBase                     0.8.0\n",
      " - StatsFuns                     0.2.0\n",
      " - TexExtensions                 0.0.3\n",
      " - Tk                            0.3.7\n",
      " - URIParser                     0.1.3\n",
      " - Winston                       0.11.13\n",
      " - WoodburyMatrices              0.1.5\n",
      " - ZMQ                           0.3.1\n",
      " - Zlib                          0.1.12\n"
     ]
    }
   ],
   "source": [
    "versioninfo(true)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: New definition \n",
      "    write(GZip.GZipStream, Array{#T<:Any, N<:Any}) at /home/bates/.julia/v0.5/GZip/src/GZip.jl:456\n",
      "is ambiguous with: \n",
      "    write(Base.IO, Array{UInt8, N<:Any}) at io.jl:154.\n",
      "To fix, define \n",
      "    write(GZip.GZipStream, Array{UInt8, N<:Any})\n",
      "before the new definition.\n",
      "WARNING: Method definition (::Type{Graphs.KeyVertex})(Int64, #K<:Any) in module Graphs at /home/bates/.julia/v0.5/Graphs/src/common.jl:12 overwritten at /home/bates/.julia/v0.5/Graphs/src/common.jl:16.\n",
      "WARNING: Method definition (::Type{Graphs.Edge})(Int64, #V<:Any, #V<:Any) in module Graphs at /home/bates/.julia/v0.5/Graphs/src/common.jl:54 overwritten at /home/bates/.julia/v0.5/Graphs/src/common.jl:60.\n",
      "WARNING: Method definition (::Type{Graphs.ExEdge})(Int64, #V<:Any, #V<:Any, Base.Dict{UTF8String, Any}) in module Graphs at /home/bates/.julia/v0.5/Graphs/src/common.jl:72 overwritten at /home/bates/.julia/v0.5/Graphs/src/common.jl:83.\n",
      "WARNING: Method definition (::Type{Graphs.TargetIterator})(#G<:Graphs.AbstractGraph, #EList<:Any) in module Graphs at /home/bates/.julia/v0.5/Graphs/src/common.jl:123 overwritten at /home/bates/.julia/v0.5/Graphs/src/common.jl:127.\n",
      "WARNING: Method definition (::Type{Graphs.SourceIterator})(#G<:Graphs.AbstractGraph, #EList<:Any) in module Graphs at /home/bates/.julia/v0.5/Graphs/src/common.jl:141 overwritten at /home/bates/.julia/v0.5/Graphs/src/common.jl:145.\n",
      "WARNING: Method definition edge_property_requirement(Graphs.AbstractEdgePropertyInspector{#T<:Any}, Graphs.AbstractGraph{#V<:Any, E<:Any}) in module Graphs at /home/bates/.julia/v0.5/Graphs/src/common.jl:164 overwritten at /home/bates/.julia/v0.5/Graphs/src/common.jl:179.\n",
      "WARNING: Method definition vertex_index(#V<:Union{Graphs.ExVertex, Graphs.KeyVertex}, Graphs.GenericGraph{#V<:Union{Graphs.ExVertex, Graphs.KeyVertex}, E<:Any, VList<:Any, EList<:Any, IncList<:Any}) in module Graphs at /home/bates/.julia/v0.5/Graphs/src/graph.jl:65 overwritten at /home/bates/.julia/v0.5/Graphs/src/graph.jl:67.\n",
      "WARNING: Method definition (::Type{Graphs.GDistanceVisitor})(#G<:Graphs.AbstractGraph, #DMap<:Any) in module Graphs at /home/bates/.julia/v0.5/Graphs/src/breadth_first_visit.jl:107 overwritten at /home/bates/.julia/v0.5/Graphs/src/breadth_first_visit.jl:111.\n",
      "WARNING: Method definition reset!(Mamba.ChainProgress) in module Mamba at /home/bates/.julia/v0.5/Mamba/src/progress.jl:34 overwritten at /home/bates/.julia/v0.5/Mamba/src/progress.jl:41.\n",
      "WARNING: Method definition Slice(Union{Symbol, Array{Symbol, 1}}, Union{#T<:Real, Array{#T<:Real, 1}}) in module Mamba at /home/bates/.julia/v0.5/Mamba/src/samplers/slice.jl:51 overwritten at /home/bates/.julia/v0.5/Mamba/src/samplers/slice.jl:125.\n",
      "WARNING: Method definition #Slice(Array, Mamba.#Slice, Union{Symbol, Array{Symbol, 1}}, Union{Array{#T<:Real, 1}, #T<:Real}) in module Mamba overwritten.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "ZeroMeanIsoNormal(\n",
       "dim: 2\n",
       "μ: [0.0,0.0]\n",
       "Σ: 2x2 Array{Float64,2}:\n",
       " 1.0e6  0.0  \n",
       " 0.0    1.0e6\n",
       ")\n"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using Distributions, GLM, Mamba, PDMats, StatsBase\n",
    "d = MvNormal(2, 1000.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you take apart the representation itself you discover that the only values that are stored are the scalar $\\sigma^2$ and its inverse."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Array{Symbol,1}:\n",
       " :μ\n",
       " :Σ"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fieldnames(d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Distributions.ZeroVector{Float64}"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "typeof(d.μ)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "PDMats.ScalMat{Float64}"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "typeof(d.Σ)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{Symbol,1}:\n",
       " :dim      \n",
       " :value    \n",
       " :inv_value"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fieldnames(d.Σ)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2,1.0e6,1.0e-6)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(d.Σ.dim, d.Σ.value, d.Σ.inv_value)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sampling without a prior\n",
    "\n",
    "Let's defer the issue of a prior for a moment and consider that we have the matrix $\\bf X'X$, the vector $\\bf X'y$ and the scalar $\\sigma$ defining the log-likelihood.  Without the prior, the mean is"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "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": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X = hcat(ones(5), [1.:5;])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "PDMats.PDMat{Float64,Array{Float64,2}}(2,2x2 Array{Float64,2}:\n",
       "  5.0  15.0\n",
       " 15.0  55.0,Base.LinAlg.Cholesky{Float64,Array{Float64,2}} with factor:\n",
       "2x2 UpperTriangular{Float64,Array{Float64,2}}:\n",
       " 2.23607  6.7082 \n",
       "  ⋅       3.16228)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "XtX = PDMat(X'X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Array{Float64,1}:\n",
       " 15.0\n",
       " 53.0"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y = [1.,3,3,3,5];\n",
    "Xty = X'y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Array{Float64,1}:\n",
       " 0.6\n",
       " 0.8"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "β = XtX\\Xty"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We could obtain the residual sum of squares by forming $\\mu = \\bf X'\\beta$ but there is a short cut that we can use later with mixed-effects models.  We append the column of $\\bf y$ to the $\\bf X$ and then form the Cholesky decomposition."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5x3 Array{Float64,2}:\n",
       " 1.0  1.0  1.0\n",
       " 1.0  2.0  3.0\n",
       " 1.0  3.0  3.0\n",
       " 1.0  4.0  3.0\n",
       " 1.0  5.0  5.0"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "const Xy = hcat(X,y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "PDMats.PDMat{Float64,Array{Float64,2}}(3,3x3 Array{Float64,2}:\n",
       "  5.0  15.0  15.0\n",
       " 15.0  55.0  53.0\n",
       " 15.0  53.0  53.0,Base.LinAlg.Cholesky{Float64,Array{Float64,2}} with factor:\n",
       "3x3 UpperTriangular{Float64,Array{Float64,2}}:\n",
       " 2.23607  6.7082   6.7082 \n",
       "  ⋅       3.16228  2.52982\n",
       "  ⋅        ⋅       1.26491)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pp = PDMat(Xy'Xy)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This `PDMat` object contains the positive-definite matrix itself and its (upper) Cholesky factor."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{Symbol,1}:\n",
       " :dim \n",
       " :mat \n",
       " :chol"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fieldnames(pp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3x3 Array{Float64,2}:\n",
       "  5.0  15.0  15.0\n",
       " 15.0  55.0  53.0\n",
       " 15.0  53.0  53.0"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pp.mat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Base.LinAlg.Cholesky{Float64,Array{Float64,2}} with factor:\n",
       "3x3 UpperTriangular{Float64,Array{Float64,2}}:\n",
       " 2.23607  6.7082   6.7082 \n",
       "  ⋅       3.16228  2.52982\n",
       "  ⋅        ⋅       1.26491"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pp.chol"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As with many of the [factorizations](http://docs.julialang.org/en/latest/manual/linear-algebra/) available in Julia ([see also](http://docs.julialang.org/en/latest/stdlib/linalg/#stdlib-linalg)), the `Cholesky` factorization provides the factor itself by indexing with a symbol."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3x3 UpperTriangular{Float64,Array{Float64,2}}:\n",
       " 2.23607  6.7082   6.7082 \n",
       "  ⋅       3.16228  2.52982\n",
       "  ⋅        ⋅       1.26491"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "R = pp.chol[:U]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The [3, 3] element is the square root of the residual sum of squares."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.6000000000000012"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "abs2(R[3, 3]) # abs2 is x^2 as a function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.5999999999999992"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using StatsBase\n",
    "sqL2dist(y, X * β)  # squared L₂ distance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "true"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sqL2dist(y, X * β) ≈ abs2(R[3, 3])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The least squares estimates of the coefficients are obtained as"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "UpperTriangular{Float64,Array{Float64,2}}"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "R12 = UpperTriangular(R[1:2, 1:2]);\n",
    "typeof(R12)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Array{Float64,1}:\n",
       " 0.6\n",
       " 0.8"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bb = R12 \\ R[1:2, 3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Furthermore, the covariance matrix of the least squares parameter estimator is "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2x2 UpperTriangular{Float64,Array{Float64,2}}:\n",
       " 0.447214  -0.948683\n",
       "  ⋅         0.316228"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "R12inv = inv(R12)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "PDMats.PDMat{Float64,Array{Float64,2}}(2,2x2 Array{Float64,2}:\n",
       "  1.1  -0.3\n",
       " -0.3   0.1,Base.LinAlg.Cholesky{Float64,Array{Float64,2}} with factor:\n",
       "2x2 UpperTriangular{Float64,Array{Float64,2}}:\n",
       " 1.04881  -0.286039\n",
       "  ⋅        0.13484 )"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Σ₁ = unscaledΣ = PDMat(R12inv * R12inv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "PDMats.PDMat{Float64,Array{Float64,2}}(2,2x2 Array{Float64,2}:\n",
       "  1.1  -0.3\n",
       " -0.3   0.1,Base.LinAlg.Cholesky{Float64,Array{Float64,2}} with factor:\n",
       "2x2 UpperTriangular{Float64,Array{Float64,2}}:\n",
       " 1.04881  -0.286039\n",
       "  ⋅        0.13484 )"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "inv(PDMat(X'X))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We want to sample from a $\\mathcal{N}(\\beta, \\sigma^2 \\Sigma_1)$ distribution for the Gibbs update of $\\beta$.\n",
    "\n",
    "## Sampling from a multivariate normal distribution\n",
    "\n",
    "At this point we want to sample from\n",
    "If you trace back through methods for functions like `rand`, its mutating version `rand!` and its hidden, mutating version without dimension checks, `_rand!`, you will find that a call to\n",
    "```julia\n",
    "rand(MvNormal(β, abs2(σ) * Σ₁))\n",
    "```\n",
    "eventually calls\n",
    "```julia\n",
    "_rand!(d::MvNormal, x::VecOrMat{Float64}) = add!(unwhiten!(d.Σ, randn!(x)), d.μ)\n",
    "```\n",
    "where `randn!(x)` overwrites `x` with $\\mathcal{N}(0, I)$ random variates.\n",
    "\n",
    "This brings us to the question of what is \"unwhitening\"?  Well, it is the inverse of \"whitening\" which is producing \"white noise\", i.e. uncorrelated, constant variance multivariate normal disturbances, from correlated, non-constant variance disturbances.  \n",
    "\n",
    "This is the multivariate equivalent of the inverse of the univariate \"z transformation\"\n",
    "$$\n",
    "x = \\mu + \\sigma z\n",
    "$$\n",
    "\n",
    "For the multivariate case we need to multiply a vector $\\mathbf{z}$ by $\\sigma$ and a \"square root\" of $\\Sigma_1$, which is given by its Cholesky factor, $R_1$.  That is $\\mathrm{Var}(R_1\\mathbf{z})$ is $R_1'\\mathrm{Var}(\\mathbf{zz}')R_1=R_1'R_1=\\Sigma_1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.7302967433402218"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "σ = R[3, 3]/√3"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.5.0-dev",
   "language": "julia",
   "name": "julia-0.5"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.5.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}