{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Bayesian Inference by Application to Simple Linear Regression\n",
"=============================================================\n",
"\n",
"Simple linear regression is used to illustrate Bayesian inference, using\n",
"the Gibbs sampler. The Gibbs sampler is used to draw samples from the\n",
"posterior distribution of the intercept, the slope and the residual\n",
"variance.\n",
"\n",
"The Model\n",
"---------\n",
"\n",
"Consider the linear model:\n",
"\n",
"$$y_{i}=\\beta_{0}+x_{i}\\beta_{1}+e_{i}. \\tag{35}$$\n",
"\n",
"where for observation $i$, $y_{i}$ is the value of the dependent\n",
"variable, $\\beta_{0}$ is the intercept, $x_{i}$ is the value of the\n",
"independent variable and $e_{i}$ is a residual. Flat priors are used for\n",
"the intercept and slope, and the residuals are assumed to be identically\n",
"and independently distributed normal random variables with mean zero and\n",
"variance $\\sigma_{e}^{2}.$ A scaled inverted chi-square prior is used\n",
"for $\\sigma_{e}^{2}.$\n",
"\n",
"Simulation of Data\n",
"------------------"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"using Distributions\n",
"using StatsBase"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"n = 200 #number of observations\n",
"k = 1 #number of covariates\n",
"\n",
"x = sample([0,1,2],(n,k))\n",
"X = [ones(Int64,n) x]\n",
"\n",
"betaTrue = [1,2]\n",
"y = X*betaTrue+ randn(n);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Least Squares Estimation\n",
"------------------------\n",
"\n",
"In matrix notation, the model (35) is\n",
"\n",
"$$\\mathbf{y}=\\mathbf{X}\\mathbf{\\beta}+\\mathbf{e},$$ where\n",
"$$\\mathbf{X}=\\left[\\begin{array}{cc}\n",
"1 & x_{1}\\\\\n",
"1 & x_{2}\\\\\n",
"\\vdots & \\vdots\\\\\n",
"1 & x_{n}\n",
"\\end{array}\\right].$$ Then, the least-squares estimator of\n",
"$\\mathbf{\\beta}$ is\n",
"$$\\mathbf{\\hat{\\mathbf{\\beta}}}=(\\mathbf{X}'\\mathbf{X})^{-1}\\mathbf{X}'\\mathbf{y},$$\n",
"and the variance of this estimator is\n",
"$$Var(\\hat{\\mathbf{\\mathbf{\\beta}}})=(\\mathbf{X}'\\mathbf{X})^{-1}\\sigma_{e}^{2}.$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Calculations in Julia:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[0.014328960645812312 -0.008678102926337033\n",
" -0.008678102926337033 0.008072653884964682]\n"
]
}
],
"source": [
"XPX = X'X\n",
"rhs = X'y\n",
"XPXi= inv(XPX)\n",
"println(XPXi)"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[0.8871642409759968,2.0840804836041267]\n"
]
}
],
"source": [
"betaHat = XPXi*rhs\n",
"println(betaHat)"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[1.02798325102907]\n"
]
}
],
"source": [
"eHat = y - X*betaHat\n",
"resVar = eHat'eHat/(n-2)\n",
"println(resVar)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Bayesian Inference\n",
"------------------\n",
"\n",
"Consider making inferences about $\\mathbf{\\beta}$ from\n",
"$f(\\mathbf{\\beta}|\\mathbf{y},\\sigma_{e}^{2}).$ By using the Bayes\n",
"theorem, this conditional density is written as\n",
"\n",
"\n",
"$$\\begin{eqnarray}\n",
"f(\\mathbf{\\beta}|\\mathbf{y},\\sigma_{e}^{2}) & = & \\frac{f(\\mathbf{y}|\\mathbf{\\beta},\\sigma_{e}^{2})f(\\mathbf{\\beta})f(\\sigma_{e}^{2})}{f(\\mathbf{y},\\sigma_{e}^{2})}\\nonumber \\\\\n",
" & \\propto & f(\\mathbf{y}|\\mathbf{\\beta},\\sigma_{e}^{2})f(\\mathbf{\\beta})f(\\sigma_{e}^{2})\\nonumber \\\\\n",
" & \\propto & f(\\mathbf{y}|\\mathbf{\\beta},\\sigma_{e}^{2})\\nonumber \\\\\n",
" & = & (2\\pi\\sigma_{e}^{2})^{-n/2}\\exp\\left\\{ -\\frac{1}{2}\\frac{(\\mathbf{y}-\\mathbf{X}\\mathbf{\\beta})'(\\mathbf{y}-\\mathbf{X}\\mathbf{\\beta})}{\\sigma_{e}^{2}}\\right\\}\\tag{36}\n",
"\\end{eqnarray}$$\n",
"\n",
"which looks like the $n$-dimensional normal density of $\\mathbf{y}$ with\n",
"mean $\\mathbf{X}\\mathbf{\\beta}$ and covariance matrix\n",
"$\\mathbf{I}\\sigma_{e}^{2}.$ But,\n",
"$f(\\mathbf{\\beta}|\\mathbf{y},\\sigma_{e}^{2})$ should be a\n",
"two-dimensional density. So, the quadratic\n",
"$Q=(\\mathbf{y}-\\mathbf{X}\\mathbf{\\beta})'(\\mathbf{y}-\\mathbf{X}\\mathbf{\\beta})$\n",
"in the exponent of (36) is rearranged as\n",
"\n",
"$$\\begin{eqnarray}\n",
"Q & = & (\\mathbf{y}-\\mathbf{X}\\mathbf{\\beta})'(\\mathbf{y}-\\mathbf{X}\\mathbf{\\beta})\\\\\n",
" & = & \\mathbf{y}'\\mathbf{y}-2\\mathbf{y}'\\mbox{X}\\mathbf{\\beta}+\\mathbf{\\beta}'(\\mathbf{X}'\\mathbf{X})\\mathbf{\\beta}\\\\\n",
" & = & \\mathbf{y}'\\mathbf{y}+(\\mathbf{\\beta}-\\hat{\\mathbf{\\beta}})'(\\mathbf{X}'\\mathbf{X})(\\mathbf{\\beta}-\\hat{\\mathbf{\\beta}})-\\hat{\\mathbf{\\beta}}'(\\mathbf{X}'\\mathbf{X})\\hat{\\mathbf{\\beta}},\n",
"\\end{eqnarray}$$\n",
"\n",
"where $\\hat{\\mathbf{\\beta}}$ is the solution to\n",
"$(\\mathbf{X}'\\mathbf{X})\\hat{\\mathbf{\\beta}}=\\mathbf{X}'\\mathbf{y}$,\n",
"which is the least-squares estimator of $\\mathbf{\\beta}$. In this\n",
"expression, only the second term depends on $\\mathbf{\\beta}$. Thus,\n",
"$f(\\mathbf{\\beta}|\\mathbf{y},\\sigma_{e}^{2})$ can be written as\n",
"\n",
"$$\n",
"f(\\mathbf{\\beta}|\\mathbf{y},\\sigma_{e}^{2})\\propto\\exp\\left\\{ -\\frac{1}{2}\\frac{(\\mathbf{\\beta}-\\hat{\\mathbf{\\beta}})'(\\mathbf{X}'\\mathbf{X})(\\mathbf{\\beta}-\\hat{\\mathbf{\\beta}})}{\\sigma_{e}^{2}}\\right\\},\n",
"$$\n",
"\n",
"which can be recognized as proportional to the density for a\n",
"two-dimensional normal distribution with mean $\\hat{\\mathbf{\\beta}}$\n",
"and variance $(\\mathbf{X}'\\mathbf{X})^{-1}\\sigma_{e}^{2}$. Thus, in this\n",
"simple setting, the posterior mean of $\\mathbf{\\beta}$ is given by\n",
"the least-squares estimate, and drawing samples from the posterior are\n",
"not needed. But, to illustrate the Gibbs sampler, we will apply it to\n",
"this simple example."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Gibbs Sampler for $\\mathbf{\\beta}$\n",
"\n",
"The simple regression model can be written as\n",
"$$\\mathbf{y=}\\mathbf{1}\\beta_{0}+\\mathbf{x}\\beta_{1}+\\mathbf{e}.$$ In\n",
"the Gibbs sampler, $\\beta_{0}$ is sampled from its full-conditional\n",
"posterior: $f(\\beta_{0}|\\mathbf{y},\\beta_{1},\\sigma_{e}^{2}).$ This\n",
"conditional distribution is computed for the current values of\n",
"$\\beta_{1}$ and $\\sigma_{e}^{2}.$ So, we can write the model as\n",
"$$\\mathbf{w}_{0}=\\mathbf{1}\\beta_{0}+\\mathbf{e},$$ where\n",
"$\\mathbf{w}_{0}=\\mathbf{y}-\\mathbf{x}\\beta_{1}.$ Then, the least-squares\n",
"estimator of $\\beta_{0}$ is\n",
"$$\\hat{\\beta}_{0}=\\frac{\\mathbf{1}'\\mathbf{w}_{0}}{\\mathbf{1}'\\mathbf{1}},$$\n",
"and the variance of this estimator is\n",
"$$Var(\\hat{\\beta_{0}})=\\frac{\\sigma_{e}^{2}}{\\mathbf{1}'\\mathbf{1}}.$$\n",
"By applying the strategy used to derive\n",
"$f(\\mathbf{\\beta}|\\mathbf{y},\\sigma_{e}^{2})$ above, the\n",
"full-conditional posterior for $\\beta_{0}$ can be shown to be a normal\n",
"distribution with mean $\\hat{\\beta_{0}}$ and variance\n",
"$\\frac{\\sigma_{e}^{2}}{\\mathbf{1}'\\mathbf{1}}.$ Similarly, the\n",
"full-conditional posterior for $\\beta_{1}$ is a normal distribution with\n",
"mean\n",
"$$\\hat{\\beta}_{1}=\\frac{\\mathbf{x}'\\mathbf{w}_{1}}{\\mathbf{x}'\\mathbf{x}}$$\n",
"and variance $\\frac{\\sigma_{e}^{2}}{\\mathbf{x}'\\mathbf{x}}$, where\n",
"$\\mathbf{w}_{1}=\\mathbf{y}-1\\beta_{0}.$ In the calculations below, we\n",
"will use the true value of $\\sigma_{e}^{2}.$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Calculations in Julia:"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Intercept = 0.904 \n",
"Slope = 2.071 \n",
"Intercept = 0.899 \n",
"Slope = 2.074 \n",
"Intercept = 0.896 \n",
"Slope = 2.076 \n",
"Intercept = 0.892 \n",
"Slope = 2.080 \n",
"Intercept = 0.892 \n",
"Slope = 2.080 \n",
"Intercept = 0.893 \n",
"Slope = 2.079 \n",
"Intercept = 0.892 \n",
"Slope = 2.081 \n",
"Intercept = 0.891 \n",
"Slope = 2.081 \n",
"Intercept = 0.891 \n",
"Slope = 2.081 \n",
"Intercept = 0.889 \n",
"Slope = 2.082 \n"
]
}
],
"source": [
"# loop for Gibbs sampler\n",
"niter = 10000 # number of samples\n",
"b = [0.0, 0.0] # initial value of b\n",
"meanB = [0.0, 0.0]\n",
"a=Float64[]\n",
"\n",
"\n",
"for iter = 1:niter\n",
" \n",
" # sampling intercept\n",
" w = y - X[:,2] * b[2]\n",
" x = X[:,1]\n",
" xpxi = 1/(x'x)[1]\n",
" bHat = (xpxi*x'w)[1]\n",
" b[1] = rand(Normal(bHat, sqrt(xpxi))) # using residual var = 1 \n",
" \n",
" # sampling slope\n",
" w = y - X[:,1]*b[1]\n",
" x = X[:,2]\n",
" xpxi = 1/(x'x)[1]\n",
" bHat = (xpxi*x'w)[1]\n",
" b[2] = rand(Normal(bHat, sqrt(xpxi))) # using residual var = 1 \n",
" meanB = meanB + b \n",
" push!(a,b[2])\n",
" \n",
" if ((iter%1000) == 0)\n",
" @printf(\"Intercept = %6.3f \\n\", meanB[1]/iter)\n",
" @printf(\"Slope = %6.3f \\n\", meanB[2]/iter)\n",
" end\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"using Gadfly"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
],
"text/html": [
"\n",
"\n"
],
"text/plain": [
"Plot(...)"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(x=a, Geom.histogram, \n",
"Guide.title(\"Posterior distribution of β1\"),\n",
"Guide.ylabel(\"Frequency\"),\n",
"Guide.xlabel(\"β1\"))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Full-conditional Posterior for $\\sigma_{e}^{2}$\n",
"\n",
"Recall that we assumed a scaled inverted chi-square prior for\n",
"$\\sigma_{e}^{2}$. The density function for this is:\n",
"\n",
"$$f(\\sigma_{e}^{2})=\\frac{(S_{e}^{2}\\nu_{e}/2)^{\\nu_{e}/2}}{\\Gamma(\\nu_{e}/2)}(\\sigma_{e}^{2})^{-(2+\\nu_{e})/2}\\exp\\left\\{ -\\frac{\\nu_{e}S_{e}^{2}}{2\\sigma_{e}^{2}}\\right\\},\\tag{37}$$\n",
"\n",
"where $S_{e}^{2}$ and $\\nu_{e}$ are the scale and the degrees of freedom\n",
"parameters for this distribution. Applying Bayes theorem to combine this\n",
"prior with the “likelihood” (given in (36)), the full-conditional\n",
"posterior for the residual variance can be written as\n",
"\n",
"$$\\begin{aligned}\n",
"f(\\sigma_{e}^{2}|\\mathbf{y},\\mathbf{\\beta}) \n",
"& = \\frac{f(\\mathbf{y}|\\mathbf{\\beta},\\sigma_{e}^{2})f(\\mathbf{\\beta})f(\\sigma_{e}^{2})}{f(\\mathbf{y},\\mathbf{\\beta})}\\nonumber \\\\\n",
"& \\propto f(\\mathbf{y}|\\mathbf{\\beta},\\sigma_{e}^{2})f(\\mathbf{\\beta})f(\\sigma_{e}^{2})\\nonumber \\\\\n",
"& \\propto (\\sigma_{e}^{2})^{-n/2}\\exp\\left\\{ -\\frac{1}{2}\\frac{(\\mathbf{y}-\\mathbf{X}\\mathbf{\\beta})'(\\mathbf{y}-\\mathbf{X}\\mathbf{\\beta})}{\\sigma_{e}^{2}}\\right\\} \\nonumber \\\\\n",
"& \\times (\\sigma_{e}^{2})^{-(2+\\nu_{e})/2}\\exp\\left\\{ -\\frac{\\nu_{e}S_{e}^{2}}{2\\sigma_{e}^{2}}\\right\\} \\nonumber \\\\\n",
"& = (\\sigma_{e}^{2})^{-(n+2+\\nu_{e})/2}\\exp\\left\\{ -\\frac{(\\mathbf{y}-\\mathbf{X}\\mathbf{\\beta})'(\\mathbf{y}-\\mathbf{X}\\mathbf{\\beta})+\\nu_{e}S_{e}^{2}}{2\\sigma_{e}^{2}}\\right\\}.\n",
"\\end{aligned}\\tag{38}$$\n",
"\n",
"Comparing (38) with (37), can see that it is proportional\n",
"to a scaled inverse chi-squared density with $\\tilde{\\nu}_{e}=n+\\nu_{e}$\n",
"degrees of freedom and\n",
"$\\tilde{S_{e}^{2}}=\\frac{(\\mathbf{y}-\\mathbf{X}\\mathbf{\\beta})'(\\mathbf{y}-\\mathbf{X}\\mathbf{\\beta})+\\nu_{e}S_{e}^{2}}{\\tilde{\\nu_{e}}}$\n",
"scale parameter. A sample from this density can be obtained as\n",
"$\\frac{(\\mathbf{y}-\\mathbf{X}\\mathbf{\\beta})'(\\mathbf{y}-\\mathbf{X}\\mathbf{\\beta})+\\nu_{e}S_{e}^{2}}{\\chi_{\\tilde{\\nu_{e}}}^{2}}$,\n",
"where $\\chi_{\\tilde{\\nu_{e}}}^{2}$ is a chi-squared random variable with\n",
"$\\tilde{\\nu_{e}}$ degrees of freedom. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Exercise\n",
"\n",
"In the Julia script given here, the simulated value of the residual variance\n",
"was used in the sampling of $\\mathbf{\\beta}$. Extend this script to\n",
"also sample $\\sigma_{e}^{2}$ from its full-conditional posterior given\n",
"above. In Julia, rand(Chisq($\\nu$),1) gives a chi-squared random variable with\n",
"$\\nu$ degrees of freedom."
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
],
"text/html": [
"\n",
"\n"
],
"text/plain": [
"Plot(...)"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# loop for Gibbs sampler\n",
"niter = 10000 # number of samples\n",
"b = [0.0, 0.0] # initial value of b\n",
"meanB = [0.0, 0.0]\n",
"a = Float64[]\n",
"vRes = 1.0\n",
"νRes = 5\n",
"S2Res = vRes*(νRes - 1)/νRes\n",
"n = size(X,1)\n",
"for iter = 1:niter\n",
" \n",
" # sampling intercept\n",
" w = y - X[:,2] * b[2]\n",
" x = X[:,1]\n",
" xpxi = 1/(x'x)[1]\n",
" bHat = (xpxi*x'w)[1]\n",
" b[1] = rand(Normal(bHat, sqrt(xpxi*vRes))) \n",
" \n",
" # sampling slope\n",
" w = y - X[:,1]*b[1]\n",
" x = X[:,2]\n",
" xpxi = 1/(x'x)[1]\n",
" bHat = (xpxi*x'w)[1]\n",
" b[2] = rand(Normal(bHat, sqrt(xpxi*vRes)))\n",
" meanB = meanB + b \n",
" \n",
" # sampling vRes\n",
" e = y - X*b\n",
" SSE = (e'e)[1]\n",
" vRes = (SSE + νRes*S2Res)/rand(Chisq(n+νRes),1)[1]\n",
" push!(a,vRes)\n",
"\n",
"end\n",
"plot(x=a, Geom.histogram, \n",
"Guide.title(\"Posterior distribution of vREs\"),\n",
"Guide.ylabel(\"Frequency\"),\n",
"Guide.xlabel(\"vRes\"))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Model with Normal Prior for Slope\n",
"---------------------------------\n",
"\n",
"Consider the simple regression model that can be written as\n",
"$$\\mathbf{y=}\\mathbf{1}\\beta_{0}+\\mathbf{x}\\beta_{1}+\\mathbf{e}.$$ \n",
"Here we consider a model with a flat prior for $\\beta_{0}$ and a normal\n",
"prior for the slope: \n",
"\n",
"$$\\beta_{1}\\sim N(0,\\sigma_{\\beta}^{2}),$$ \n",
"\n",
"where $\\sigma_{\\beta}^{2}$ is assumed to be known. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Then, the full-conditional\n",
"posterior for $\\theta'=[\\mathbf{\\beta},\\sigma_{e}^{2}]$ is\n",
"\n",
"$$\\begin{aligned}\n",
"f(\\mathbf{\\theta}|\\mathbf{y}) \n",
"& \\propto f(\\mathbf{y}|\\mathbf{\\theta})f(\\mathbf{\\theta})\\\\\n",
"& \\propto \\left(\\sigma_{e}^{2}\\right)^{-n/2}\\exp\\left\\{ -\\frac{(\\mathbf{y}-\\mathbf{1}\\beta_{0}-\\mathbf{x}\\beta_{1})'(\\mathbf{y}-\\mathbf{1}\\beta_{0}-\\mathbf{x}\\beta_{1})}{2\\sigma_{e}^{2}}\\right\\} \\\\\n",
"& \\times \\left(\\sigma_{\\beta}^{2}\\right)^{-1/2}\\exp\\left\\{ -\\frac{\\beta_{1}^{2}}{2\\sigma_{\\beta}^{2}}\\right\\} \\\\\n",
"& \\times (\\sigma_{e}^{2})^{-(2+\\nu_{e})/2}\\exp\\left\\{ -\\frac{\\nu_{e}S_{e}^{2}}{2\\sigma_{e}^{2}}\\right\\} .\\\\\n",
"\\end{aligned}$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Full-conditional for $\\beta_{1}$:\n",
"\n",
"The full-conditional for $\\beta_{1}$ is obtained by dropping all terms\n",
"and factors that do not involve $\\beta_{1}$:\n",
"\n",
"$$\\begin{aligned}\n",
"f(\\beta_{1}|\\text{ELSE}) \n",
"& \\propto \\exp\\left\\{ -\\frac{(\\mathbf{y}-\\mathbf{1}\\beta_{0}-\\mathbf{x}\\beta_{1})'(\\mathbf{y}-\\mathbf{1}\\beta_{0}-\\mathbf{x}\\beta_{1})}{2\\sigma_{e}^{2}}\\right\\} \\times \\exp\\left\\{ -\\frac{\\beta_{1}^{2}}{2\\sigma_{\\beta}^{2}}\\right\\} \\\\\n",
"& \\propto \\exp\\left\\{ -\\frac{\\mathbf{w}'\\mathbf{w}-2\\mathbf{w}'\\mathbf{x}\\beta_{1}+\\beta_{1}^{2}(\\mathbf{x}'\\mathbf{x}+\\sigma_{e}^{2}/\\sigma_{\\beta}^{2})}{2\\sigma_{e}^{2}}\\right\\} \\\\\n",
"& \\propto \\exp\\left\\{ -\\frac{\\mathbf{w}'\\mathbf{w}-(\\beta_{1}-\\hat{\\beta}_{1})^{2}(\\mathbf{x}'\\mathbf{x}+\\sigma_{e}^{2}/\\sigma_{\\beta}^{2})-\\hat{\\beta}_{1}^{2}(\\mathbf{x}'\\mathbf{x}+\\sigma_{e}^{2}/\\sigma_{\\beta}^{2})}{2\\sigma_{e}^{2}}\\right\\} \\\\\n",
"& \\propto \\exp\\left\\{ -\\frac{(\\beta_{1}-\\hat{\\beta}_{1})^{2}}{\\frac{2\\sigma_{e}^{2}}{(\\mathbf{x}'\\mathbf{x}+\\sigma_{e}^{2}/\\sigma_{\\beta}^{2})}}\\right\\} ,\n",
"\\end{aligned}$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"where\n",
"$$\\hat{\\beta}_{1}=\\frac{\\mathbf{x}'\\mathbf{w}}{(\\mathbf{x}'\\mathbf{x}+\\sigma_{e}^{2}/\\sigma_{\\beta}^{2})},$$\n",
"and $\\mathbf{w}=\\mathbf{y}-\\mathbf{1}\\beta_{0}.$ So, the\n",
"full-conditional posterior for $\\beta_{1}$ is a normal distribution with\n",
"$ $mean $\\hat{\\beta}_{1}$ and variance\n",
"$\\frac{\\sigma_{e}^{2}}{(\\mathbf{x}'\\mathbf{x}+\\sigma_{e}^{2}/\\sigma_{\\beta}^{2})}.$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Exercise\n",
"\n",
" 1. In the following we will use a normal distribution with mean zero and variance 3 as the prior \n",
" for $\\beta_{1}$. Use Julia to simulate a vector of 1000 values for $\\beta_{1}$ from a\n",
" this prior. Plot a histogram of these values. \n",
"\n",
"2. Use $\\beta_{0}=1$, $\\beta_{1}=2$ and $\\sigma_{e}^{2}=5$, to generate a vector of observations, $\\mathbf{y},$ \n",
" that follows a simple linear regression model. \n",
"\n",
"3. Use the Gibbs sampler to draw 10,000 samples for $\\beta_{1}$ from\n",
" its posterior distribution.\n",
"\n",
" 1. Compute the mean and variance of the sampled values.\n",
"\n",
" 2. Draw a histogram of the sampled values. Compare with prior. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Simulation of the data"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"n = 200 # number of observations\n",
"k = 1 # number of covariates\n",
"\n",
"x = sample([0,1,2],(n,k))\n",
"X = [ones(Int64,n) x]\n",
"\n",
"betaTrue = [1,2]\n",
"y = X*betaTrue+ randn(n)*sqrt(5.0);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Prior Distibution of $\\beta_1$"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {
"collapsed": false,
"scrolled": true
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
],
"text/html": [
"\n",
"\n"
],
"text/plain": [
"Plot(...)"
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"prior = randn(10000)*sqrt(3.0)\n",
"plot(x=prior, Geom.histogram, \n",
"Guide.title(\"Prior distribution of β1\"),\n",
"Guide.ylabel(\"Frequency\"),\n",
"Guide.xlabel(\"β1\"))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Posterior Distibution of $\\beta_1$"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {
"collapsed": false,
"scrolled": true
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
],
"text/html": [
"\n",
"\n"
],
"text/plain": [
"Plot(...)"
]
},
"execution_count": 33,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# loop for Gibbs sampler\n",
"niter = 10000 # number of samples\n",
"b = [0.0, 0.0] # initial value of b\n",
"meanB = [0.0, 0.0]\n",
"a = Float64[]\n",
"vB1 = 3.0\n",
"vRes = 1.0\n",
"νRes = 5\n",
"S2Res = vRes*(νRes - 1)/νRes\n",
"n = size(X,1)\n",
"for iter = 1:niter\n",
" \n",
" λ = vRes/vB1\n",
" \n",
" # sampling intercept\n",
" w = y - X[:,2] * b[2]\n",
" x = X[:,1]\n",
" xpxi = 1/(x'x)[1]\n",
" bHat = (xpxi*x'w)[1]\n",
" b[1] = rand(Normal(bHat, sqrt(xpxi*vRes))) \n",
" \n",
" # sampling slope\n",
" w = y - X[:,1]*b[1]\n",
" x = X[:,2]\n",
" xpxi = 1/(x'x + λ)[1]\n",
" bHat = (xpxi*x'w)[1]\n",
" b[2] = rand(Normal(bHat, sqrt(xpxi*vRes)))\n",
" meanB = meanB + b \n",
" \n",
" # sampling vRes\n",
" e = y - X*b\n",
" SSE = (e'e)[1]\n",
" vRes = (SSE + νRes*S2Res)/rand(Chisq(n+νRes),1)[1]\n",
" push!(a,b[2])\n",
"\n",
"end\n",
"plot(x=a, Geom.histogram, \n",
"Guide.title(\"Posterior distribution of β1\"),\n",
"Guide.ylabel(\"Frequency\"),\n",
"Guide.xlabel(\"β1\"))"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"using DataFrames\n",
"df = DataFrame()\n",
"df[:PP] = [fill(\"Prior\",10000); fill(\"Posterior\",10000)]\n",
"df[:b1] = [prior; a];"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"