{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Bayesian Regression Models for Whole-Genome Analyses\n",
    "====================================================\n",
    "\n",
    "Meuwissen et al. (2001) introduced three regression\n",
    "models for whole-genome prediction of breeding value of the form\n",
    "\n",
    "\n",
    "$$y_{i}=\\mu+\\sum_{j=1}^{k}X_{ij}\\alpha_{j}+e_{i},$$ \n",
    "\n",
    "where $y_{i}$ is the\n",
    "phenotypic value, $\\mu$ is the intercept, $X_{ij}$ is $j^{th}$ marker\n",
    "covariate of animal $i$, $\\alpha_{j}$ is the partial regression\n",
    "coefficient of $X_{ij}$, and $e_{i}$ are identically and independently\n",
    "distributed residuals with mean zero and variance $\\sigma_{e}^{2}.$ In\n",
    "most current analyses, $X_{ij}$ are SNP genotype covariates that can be\n",
    "coded as 0, 1 and 2, depending on the number of B alleles at SNP locus\n",
    "$j$.\n",
    "\n",
    "In all three of their models, a flat prior was used for the intercept\n",
    "and a scaled inverted chi-square distribution for $\\sigma_{e}^{2}.$ The\n",
    "three models introduced by Meuwissen et al. (2001)\n",
    "differ only in the prior used for $\\alpha_{j}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "BLUP\n",
    "----\n",
    "\n",
    "In their first model, which they called BLUP, a normal distribution with\n",
    "mean zero and known variance, $\\sigma_{\\alpha}^{2}$, is used as the\n",
    "prior for $\\alpha_{j}$.\n",
    "\n",
    "### The meaning of $\\sigma_{\\alpha}^{2}$\n",
    "\n",
    "Assume the QTL are in the marker panel. Then, the genotypic value\n",
    "$g_{i}$ for a randomly sampled animal $i$ can be written as\n",
    "$$g_{i}=\\mu+\\mathbf{x}'_{i}\\boldsymbol{\\alpha},$$ where\n",
    "$\\mathbf{x}'_{i}$ is the vector of SNP genotype covariates and\n",
    "$\\boldsymbol{\\alpha}$ is the vector of regression coefficients. Note\n",
    "that randomly sampled animals differ only in $\\mathbf{x}'_{i}$ and have\n",
    "$\\boldsymbol{\\alpha}$ in common. Thus, genotypic variability is entirely\n",
    "due to variability in the genotypes of animals. So,\n",
    "$\\sigma_{\\alpha}^{2}$ is not the genetic variance at a locus\n",
    "(Fernando et al., 2007; Gianola et al., 2009)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Relationship of $\\sigma_{\\alpha}^{2}$ to genetic variance\n",
    "\n",
    "Assume loci with effect on trait are in linkage equilibrium. Then, the\n",
    "additive genetic variance is\n",
    "$$V_{A}=\\sum_{j}^{k}2p_{j}q_{j}\\alpha_{j}^{2},$$ where $p_{j}=1-q_{j}$\n",
    "is gene frequency at SNP locus $j$. Letting $U_{j}=2p_{j}q_{j}$ and\n",
    "$V_{j}=\\alpha_{j}^{2}$, $$V_{A}=\\sum_{j}^{k}U_{j}V_{j}.$$ For a randomly"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "sampled locus, covariance between $U_{j}$ and $V_{j}$ is\n",
    "$$C_{UV}=\\frac{\\sum_{j}U_{j}V_{j}}{k}-(\\frac{\\sum_{j}U_{j}}{k})(\\frac{\\sum_{j}V_{j}}{k})$$\n",
    "Rearranging this expression for $C_{UV}$ gives\n",
    "$$\\sum_{j}U_{j}V_{j}=kC_{UV}+(\\sum_{j}U_{j})(\\frac{\\sum_{j}V_{j}}{k})$$\n",
    "So,"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$V_{A}=kC_{UV}+(\\sum_{j}2p_{j}q_{j})(\\frac{\\sum_{j}\\alpha_{j}^{2}}{k}).$$\n",
    "Letting $\\sigma_{\\alpha}^{2}=\\frac{\\sum_{j}\\alpha_{j}^{2}}{k}$ gives\n",
    "$$V_{A}=kC_{UV}+(\\sum_{j}2p_{j}q_{j})\\sigma_{\\alpha}^{2}$$ and"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\sigma_{\\alpha}^{2}=\\frac{V_{A}-kC_{UV}}{\\sum_{j}2p_{j}q_{j}},$$ which\n",
    "gives $$\\sigma_{\\alpha}^{2}=\\frac{V_{A}}{\\sum_{j}2p_{j}q_{j}},$$ if gene\n",
    "frequency is independent of the effect of the gene."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Full-conditionals:\n",
    "\n",
    "The joint posterior for all the parameters is proportional to\n",
    "\n",
    "$$\\begin{aligned}\n",
    "f(\\boldsymbol{\\theta}|\\mathbf{y}) & \\propto  f(\\mathbf{y}|\\boldsymbol{\\theta})f(\\boldsymbol{\\theta})\\\\\n",
    " & \\propto  \\left(\\sigma_{e}^{2}\\right)^{-n/2}\\exp\\left\\{ \n",
    " -\\frac{(\\mathbf{y}-\\mathbf{1}\\mu-\\sum\\mathbf{X}_{j}\\alpha_{j})'(\\mathbf{y}-\\mathbf{1}\\mu-\\sum\\mathbf{X}_{j}\\alpha_{j})}{2\\sigma_{e}^{2}}\\right\\} \\\\\n",
    " & \\times \\prod_{j=1}^{k}\\left(\\sigma_{\\alpha}^{2}\\right)^{-1/2}\\exp\\left\\{ -\\frac{\\alpha_{j}^{2}}{2\\sigma_{\\alpha}^{2}}\\right\\} \\\\\n",
    " & \\times (\\sigma_{\\alpha}^{2})^{-(\\nu_{\\alpha}+2)/2}\\exp\\left\\{ -\\frac{\\nu_{\\alpha}S_{\\alpha}^{2}}{2\\sigma_{\\alpha}^{2}}\\right\\} \\\\\n",
    " & \\times (\\sigma_{e}^{2})^{-(2+\\nu_{e})/2}\\exp\\left\\{ -\\frac{\\nu_{e}S_{e}^{2}}{2\\sigma_{e}^{2}}\\right\\} ,\\end{aligned}$$\n",
    "\n",
    "where $\\boldsymbol{\\theta}$ denotes all the unknowns."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Full-conditional for $\\mu$\n",
    "\n",
    "The full-conditional for $\\mu$ is a normal distribution with mean\n",
    "$\\hat{\\mu}$ and variance $\\frac{\\sigma_{e}^{2}}{n}$, where $\\hat{\\mu}$\n",
    "is the least-squares estimate of $\\mu$ in the model\n",
    "$$\\mathbf{y-\\sum_{j=1}^{k}}\\mathbf{X}_{j}\\alpha_{j}=\\mathbf{1}\\mu+\\mathbf{e},$$\n",
    "and $\\frac{\\sigma_{e}^{2}}{n}$ is the variance of this estimator ($n$ is\n",
    "the number of observations)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Full-conditional for $\\alpha_{j}$\n",
    "\n",
    "$$\\begin{aligned}\n",
    "f(\\alpha_{j}|\\text{ELSE}) \n",
    "& \\propto \\exp\\left\\{ -\\frac{(\\mathbf{w}_{j}-\\mathbf{X}_{j}\\alpha_{j})'(\\mathbf{w}_{j}-\\mathbf{X}_{j}\\alpha_{j})}{2\\sigma_{e}^{2}}\\right\\} \\\\\n",
    "& \\times  \\exp\\left\\{ -\\frac{\\alpha_{j}^{2}}{2\\sigma_{\\alpha}^{2}}\\right\\} \\\\\n",
    "& \\propto \\exp\\left\\{ -\\frac{[\\mathbf{w}'_{j}\\mathbf{w}_{j}-2\\mathbf{w}'_{j}\\mathbf{X}_{j}\\alpha_{j}+\\alpha_{j}^{2}(\\mathbf{x}'_{j}\\mathbf{x}_{j}+\\sigma_{e}^{2}/\\sigma_{\\alpha}^{2})]}{2\\sigma_{e}^{2}}\\right\\} \\\\\n",
    "& \\propto \\exp\\left\\{ -\\frac{(\\alpha_{j}-\\hat{\\alpha_{j}})^{2}}{\\frac{2\\sigma_{e}^{2}}{(\\mathbf{x}'_{j}\\mathbf{x}_{j}+\\sigma_{e}^{2}/\\sigma_{\\alpha}^{2})}}\\right\\} ,\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where\n",
    "$$\\mathbf{w}_{j}=\\mathbf{y}-\\mathbf{1}\\mu-\\sum_{l\\neq j}\\mathbf{X}_{l}\\alpha_{l}.$$\n",
    "So, the full-conditional for $\\alpha_{j}$ is a normal distribution with\n",
    "mean\n",
    "$$\\hat{\\alpha}_{j}=\\frac{\\mathbf{X}'_{j}\\mathbf{w}_{j}}{(\\mathbf{x}'_{j}\\mathbf{x}_{j}+\\sigma_{e}^{2}/\\sigma_{\\alpha}^{2})}$$\n",
    "and variance\n",
    "$\\frac{\\sigma_{e}^{2}}{(\\mathbf{x}'_{j}\\mathbf{x}_{j}+\\sigma_{e}^{2}/\\sigma_{\\alpha}^{2})}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Full-conditional for $\\sigma_{\\alpha}^{2}$\n",
    "\n",
    "$$\\begin{aligned}\n",
    "f(\\sigma_{\\alpha}^{2}|\\text{ELSE}) \n",
    "& \\propto \\prod_{j=1}^{k}\\left(\\sigma_{\\alpha}^{2}\\right)^{-1/2}\\exp\\left\\{ -\\frac{\\alpha_{j}^{2}}{2\\sigma_{\\alpha}^{2}}\\right\\} \\\\\n",
    "& \\times  (\\sigma_{\\alpha}^{2})^{-(\\nu_{\\alpha}+2)/2}\\exp\\left\\{ -\\frac{\\nu_{\\alpha}S_{\\alpha}^{2}}{2\\sigma_{\\alpha}^{2}}\\right\\} \\\\\n",
    "& \\propto (\\sigma_{\\alpha}^{2})^{-(k+\\nu_{\\alpha}+2)/2}\\exp\\left\\{ -\\frac{\\sum_{j=1}^{k}\\alpha_{j}^{2}+\\nu_{\\alpha}S_{\\beta\\alpha}^{2}}{2\\sigma_{\\alpha}^{2}}\\right\\} ,\\end{aligned}$$\n",
    "\n",
    "and this is proportional to a scaled inverted chi-square distribution\n",
    "with $\\tilde{\\nu}_{\\alpha}=\\nu_{\\alpha}+k$ and scale parameter\n",
    "$\\tilde{S}_{\\alpha}^{2}=(\\sum_{k}\\alpha_{j}^{2}+\\nu_{\\alpha}S_{\\alpha}^{2})/\\tilde{\\nu}_{\\alpha}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Full-conditional for $\\sigma_{e}^{2}$\n",
    "\n",
    "$$\\begin{aligned}\n",
    "f(\\sigma_{e}^{2}|\\text{ELSE}) \n",
    "& \\propto \\left(\\sigma_{e}^{2}\\right)^{-n/2}\n",
    "\\exp\\left\\{ -\\frac{\n",
    "(\\mathbf{y}-\\mathbf{1}\\mu-\\sum\\mathbf{X}_{j}\\alpha_{j})'\n",
    "(\\mathbf{y}-\\mathbf{1}\\mu-\\sum\\mathbf{X}_{j}\\alpha_{j})\n",
    "}{2\\sigma_{e}^{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",
    "& \\propto (\\sigma_{e}^{2})^{-(n+2+\\nu_{e})/2}\\exp\\left\\{ -\\frac{(\\mathbf{y}-\\mathbf{1}\\mu-\\sum\\mathbf{X}_{j}\\alpha_{j})'(\\mathbf{y}-\\mathbf{1}\\mu-\\sum\\mathbf{X}_{j}\\alpha_{j})+\\nu_{e}S_{e}^{2}}{2\\sigma_{e}^{2}}\\right\\} ,\\end{aligned}$$\n",
    "\n",
    "which is proportional to a scaled inverted chi-square density with\n",
    "$\\tilde{\\nu}_{e}=n+\\nu_{e}$ degrees of freedom and\n",
    "$\\tilde{S_{e}^{2}}=\\frac{(\\mathbf{y}-\\mathbf{1}\\mu-\\sum\\mathbf{X}_{j}\\alpha_{j})'(\\mathbf{y}-\\mathbf{1}\\mu-\\sum\\mathbf{X}_{j}\\alpha_{j})+\\nu_{e}S_{e}^{2}}{\\tilde{\\nu_{e}}}$\n",
    "scale parameter."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "BayesB\n",
    "------\n",
    "\n",
    "### Model\n",
    "\n",
    "The usual model for BayesB is:\n",
    "\n",
    "$$y_{i}=\\mu+\\sum_{j=1}^{k}X_{ij}\\alpha_{j}+e_{i},$$ where the prior\n",
    "$\\mu$ is flat and the prior for $\\alpha_{j}$ is a mixture distribution:\n",
    "\n",
    "$$\\alpha_{j}=\\begin{cases}\n",
    "0 & \\text{probability}\\,\\pi\\\\\n",
    "\\sim N(0,\\sigma_{j}^{2}) & \\text{probability}\\,(1-\\pi)\n",
    "\\end{cases},\n",
    "$$ \n",
    "\n",
    "where $\\sigma_{j}^{2}$ has a scaled inverted chi-square prior\n",
    "with scale parameter $S_{\\alpha}^{2}$ and $\\nu_{\\alpha}$ degrees of\n",
    "freedom. The residual is normally distributed with mean zero and\n",
    "variance $\\sigma_{e}^{2}$, which has a scaled inverted chi-square prior\n",
    "with scale parameter $S_{e}^{2}$ and $\\nu_{e}$ degrees of freedom.\n",
    "Meuwissen et al. @Meuwissen.THE.ea.2001a gave a Metropolis-Hastings\n",
    "sampler to jointly sample $\\sigma_{j}^{2}$ and $\\alpha_{j}$. Here, we\n",
    "will show how the Gibbs sampler can be used in BayesB.\n",
    "\n",
    "In order to use the Gibbs sampler, the model is written as\n",
    "\n",
    "$$y_{i}=\\mu+\\sum_{j=1}^{k}X_{ij}\\beta_{j}\\delta_{j}+e_{i},$$ where\n",
    "$\\beta_{j}\\sim N(0,\\sigma_{j}^{2})$ and $\\delta_{j}$ is\n",
    "Bernoulli($1-\\pi$):\n",
    "\n",
    "$$\\begin{aligned}\n",
    "\\delta_{j} & = & \\begin{cases}\n",
    "0 & \\text{probability}\\,\\pi\\\\\n",
    "1 & \\text{probability}\\,(1-\\pi)\n",
    "\\end{cases}.\\\\\\end{aligned}$$\n",
    "\n",
    "Other priors are the same as in the usual model. Note that in this\n",
    "model, $\\alpha_{j}=\\beta_{j}\\delta_{j}$ has a mixture distribution as in\n",
    "the usual BayesB model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Full-conditionals:\n",
    "\n",
    "The joint posterior for all the parameters is proportional to\n",
    "\n",
    "$$\\begin{aligned}\n",
    "f(\\boldsymbol{\\theta}|\\mathbf{y}) \n",
    "& \\propto f(\\mathbf{y}|\\boldsymbol{\\theta})f(\\boldsymbol{\\theta})\\\\\n",
    "& \\propto \\left(\\sigma_{e}^{2}\\right)^{-n/2}\\exp\\left\\{ -\\frac{\n",
    " (\\mathbf{y}-\\mathbf{1}\\mu-\\sum\\mathbf{X}_{j}\\beta_{j}\\delta_{j})'\n",
    " (\\mathbf{y}-\\mathbf{1}\\mu-\\sum\\mathbf{X}_{j}\\beta_{j}\\delta_{j}) \n",
    " }{2\\sigma_{e}^{2}}\\right\\} \\\\\n",
    "& \\times \\prod_{j=1}^{k}\\left(\\sigma_{j}^{2}\\right)^{-1/2}\\exp\\left\\{ -\\frac{\\beta_{j}^{2}}{2\\sigma_{j}^{2}}\\right\\} \\\\\n",
    "& \\times \\prod_{j=1}^{k}\\pi^{(1-\\delta_{j})}(1-\\pi)^{\\delta_{j}}\\\\\n",
    "& \\times \\prod_{j=1}^{k}(\\sigma_{j}^{2})^{-(\\nu_{\\beta}+2)/2}\\exp\\left\\{ -\\frac{\\nu_{\\beta}S_{\\beta}^{2}}{2\\sigma_{j}^{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}\n",
    "$$\n",
    "\n",
    "where $\\boldsymbol{\\theta}$ denotes all the unknowns."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Full-conditional for $\\mu$\n",
    "\n",
    "The full-conditional for $\\mu$ is a normal distribution with mean\n",
    "$\\hat{\\mu}$ and variance $\\frac{\\sigma_{e}^{2}}{n}$, where $\\hat{\\mu}$\n",
    "is the least-squares estimate of $\\mu$ in the model\n",
    "$$\\mathbf{y-\\sum_{j=1}^{k}}\\mathbf{X}_{j}\\beta_{j}\\delta_{j}=\\mathbf{1}\\mu+\\mathbf{e},$$\n",
    "and $\\frac{\\sigma_{e}^{2}}{n}$ is the variance of this estimator ($n$ is\n",
    "the number of observations)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Full-conditional for $\\beta_{j}$\n",
    "\n",
    "$$\\begin{aligned}\n",
    "f(\\beta_{j}|\\text{ELSE}) \n",
    "& \\propto \\exp\\left\\{ -\\frac{(\\mathbf{w}_{j}-\\mathbf{X}_{j}\\beta_{j}\\delta_{j})'(\\mathbf{w}_{j}-\\mathbf{X}_{j}\\beta_{j}\\delta_{j})}{2\\sigma_{e}^{2}}\\right\\} \\\\\n",
    "& \\times \\exp\\left\\{ -\\frac{\\beta_{j}^{2}}{2\\sigma_{j}^{2}}\\right\\} \\\\\n",
    "& \\propto \\exp\\left\\{ -\\frac{[\\mathbf{w}'_{j}\\mathbf{w}_{j}-2\\mathbf{w}'_{j}\\mathbf{X}_{j}\\beta_{j}\\delta_{j}+\\beta_{j}^{2}(\\mathbf{x}'_{j}\\mathbf{x}_{j}\\delta_{j}+\\sigma_{e}^{2}/\\sigma_{j}^{2})]}{2\\sigma_{e}^{2}}\\right\\} \\\\\n",
    "& \\propto \\exp\\left\\{ -\\frac{(\\beta_{j}-\\hat{\\beta_{j}})^{2}}{\\frac{2\\sigma_{e}^{2}}{(\\mathbf{x}'_{j}\\mathbf{x}_{j}\\delta_{j}+\\sigma_{e}^{2}/\\sigma_{j}^{2})}}\\right\\} ,\\end{aligned}$$\n",
    "\n",
    "where\n",
    "$$\\mathbf{w}_{j}=\\mathbf{y}-\\mathbf{1}\\mu-\\sum_{l\\neq j}\\mathbf{X}_{l}\\beta_{l}\\delta_{l}.$$\n",
    "So, the full-conditional for $\\beta_{j}$ is a normal distribution with\n",
    "mean\n",
    "$$\\hat{\\beta}_{j}=\\frac{\\mathbf{X}'_{j}\\mathbf{w}_{j}\\delta_{j}}{(\\mathbf{x}'_{j}\\mathbf{x}_{j}\\delta_{j}+\\sigma_{e}^{2}/\\sigma_{j}^{2})}$$\n",
    "and variance\n",
    "$\\frac{\\sigma_{e}^{2}}{(\\mathbf{x}'_{j}\\mathbf{x}_{j}\\delta_{j}+\\sigma_{e}^{2}/\\sigma_{j}^{2})}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Full-conditional for $\\delta_{j}$\n",
    "$$\n",
    "\\Pr(\\delta_{j}=1|\\text{ELSE})\\propto\\frac{h(\\delta_{j}=1)}{h(\\delta_{j}=1)+h(\\delta_{j}=0)},\n",
    "$$\n",
    "\n",
    "where\n",
    "$$\n",
    "h(\\delta_{j})=\\pi^{(1-\\delta_{j})}(1-\\pi)^{\\delta_{j}}\\exp\\left\\{-\\frac{\n",
    "(\\mathbf{w}_{j}-\\mathbf{X}_{j}\\beta_{j}\\delta_{j})'\n",
    "(\\mathbf{w}_{j}-\\mathbf{X}_{j}\\beta_{j}\\delta_{j})\n",
    "}{2\\sigma_{e}^{2}}\\right\\}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Full-conditional for $\\sigma_{j}^{2}$\n",
    "\n",
    "$$\\begin{aligned}\n",
    "f(\\sigma_{j}^{2}|\\text{ELSE}) \n",
    "&\\propto \\left(\\sigma_{j}^{2}\\right)^{-1/2}\\exp\\left\\{ -\\frac{\\beta_{j}^{2}}{2\\sigma_{j}^{2}}\\right\\} \\\\\n",
    "&\\times (\\sigma_{j}^{2})^{-(\\nu_{\\beta}+2)/2}\\exp\\left\\{ -\\frac{\\nu_{\\beta}S_{\\beta}^{2}}{2\\sigma_{j}^{2}}\\right\\} \\\\\n",
    "&\\propto (\\sigma_{j}^{2})^{-(1+\\nu_{\\beta}+2)/2}\\exp\\left\\{ -\\frac{\\beta_{j}^{2}+\\nu_{\\beta}S_{\\beta}^{2}}{2\\sigma_{j}^{2}}\\right\\} ,\\end{aligned}$$\n",
    "\n",
    "and this is proportional to a scaled inverted chi-square distribution\n",
    "with $\\tilde{\\nu}_{j}=\\nu_{\\beta}+1$ and scale parameter\n",
    "$\\tilde{S}_{j}^{2}=(\\beta_{j}^{2}+\\nu_{\\beta}S_{\\beta}^{2})/\\tilde{\\nu}_{j}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Full-conditional for $\\sigma_{e}^{2}$\n",
    "\n",
    "$$\\begin{aligned}\n",
    "f(\\sigma_{e}^{2}|\\text{ELSE}) \n",
    "& \\propto \\left(\\sigma_{e}^{2}\\right)^{-n/2}\\exp\\left\\{ -\\frac{\n",
    "(\\mathbf{y}-\\mathbf{1}\\mu-\\sum\\mathbf{X}_{j}\\beta_{j}\\delta_{j})'\n",
    "(\\mathbf{y}-\\mathbf{1}\\mu-\\sum\\mathbf{X}_{j}\\beta_{j}\\delta_{j})\n",
    "}{2\\sigma_{e}^{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",
    "& \\propto (\\sigma_{e}^{2})^{-(n+2+\\nu_{e})/2}\\exp\\left\\{ -\\frac{(\\mathbf{y}-\\mathbf{1}\\mu-\\sum\\mathbf{X}_{j}\\beta_{j}\\delta_{j})'(\\mathbf{y}-\\mathbf{1}\\mu-\\sum\\mathbf{X}_{j}\\beta_{j}\\delta_{j})+\\nu_{e}S_{e}^{2}}{2\\sigma_{e}^{2}}\\right\\} ,\\end{aligned}$$\n",
    "\n",
    "which is proportional to a scaled inverted chi-square density with\n",
    "$\\tilde{\\nu}_{e}=n+\\nu_{e}$ degrees of freedom and\n",
    "$\\tilde{S_{e}^{2}}=\\frac{(\\mathbf{y}-\\mathbf{1}\\mu-\\sum\\mathbf{X}_{j}\\beta_{j}\\delta_{j})'(\\mathbf{y}-\\mathbf{1}\\mu-\\sum\\mathbf{X}_{j}\\beta_{j}\\delta_{j})+\\nu_{e}S_{e}^{2}}{\\tilde{\\nu_{e}}}$\n",
    "scale parameter."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### References\n",
    "\n",
    "1. R. L. Fernando, D. Habier, C. Stricker, J. C. M. Dekkers, and L. R. Totir. Genomic selection. Acta Agriculturae Scandinavica, Section A - Animal Science, 57(4):192–195, 2007. \n",
    "\n",
    "1. D Gianola, G de los Campos, W G Hill, E Manfredi, and R Fernando. Additive genetic variability and the bayesian alphabet. Genetics, 183(1):347– 363, Sep 2009. \n",
    "\n",
    "1. T. H. E. Meuwissen, B. J. Hayes, and M. E. Goddard. Prediction of total genetic value using genome-wide dense marker maps. Genetics, 157:1819– 1829, 2001."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "BayesB with Unknown $S^2_{\\beta}$\n",
    "--------------------------------"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, the scale parameter, $S^2_{\\beta}$, in scaled inverse chi-squarred prior for $\\sigma^2_j$ is treated as an unknown. \n",
    "\n",
    "The density function for this prior distribution is:\n",
    "$$\n",
    "f(\\sigma_{j}^{2}|S_{\\beta}^{2}, \\nu_{\\beta}) = \n",
    "\\frac{(S_{\\beta}^{2}\\nu_{\\beta}/2)^{\\nu_{\\beta}/2}}\n",
    "{\\Gamma(\\nu_{\\beta}/2)}(\\sigma_{j}^{2})^{-(2+\\nu_{\\beta})/2}\\exp\\left\\{ -\\frac{\\nu_{\\beta}S_{\\beta}^{2}}{2\\sigma_{j}^{2}}\\right\\},\n",
    "$$\n",
    "\n",
    "where $S_{\\beta}^{2}$ and $\\nu_{\\beta}$ are the scale and the degrees of freedom\n",
    "parameters for this distribution. The scale parameter $S_{\\beta}^{2}$ does not appear anywhere else in  $f(\\boldsymbol{\\theta}|\\mathbf{y})$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Thus, the Gamma$(S^2_{\\beta}| a,b)$ distribution with density function:\n",
    "$$ \n",
    "f(S^2_{\\beta}|a,b) =  \n",
    "\\frac\n",
    "{b^a(S^2_{\\beta})^{a-1}}\n",
    "{\\Gamma(a)}\n",
    "\\exp\\{-bS^2_{\\beta}\\}\n",
    "$$\n",
    "is a conjugate prior for $S^2_{\\beta}$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, using Gamma$(S^2_{\\beta}| a,b)$ as the prior for $S^2_{\\beta}$, its full conditional can be written as\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "f(S^2_{\\beta}|\\text{ELSE})\n",
    "& \\propto \n",
    "\\prod_{j=1}^{k}(S^2_{\\beta})^{\\nu_{\\beta}/2}\\exp\\left\\{ -\\frac{\\nu_{\\beta}S_{\\beta}^{2}}{2\\sigma_{j}^{2}}\\right\\} \\\\\n",
    "& \\times \n",
    "b^a(S^2_{\\beta})^{a-1}\\exp\\{-bS^2_{\\beta}\\}\\\\\n",
    "&\\times\n",
    "(S^2_{\\beta})^{(k\\nu_{\\beta})/2 + a - 1}\n",
    "\\exp\\left\\{ -S_{\\beta}^{2}(\\sum_{j=1}^k\\frac{\\nu_{\\beta}}{2\\sigma_{j}^{2}} + b)\\right\\},\n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "which can be recognized as a Gamma distribution with shape $a^*= (k\\nu_{\\beta})/2 + a$ and rate $b^*=(\\sum_{j=1}^k\\frac{\\nu_{\\beta}}{2\\sigma_{j}^{2}} + b)$ parameters.\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.3.7",
   "language": "julia",
   "name": "julia 0.3"
  },
  "language": "Julia",
  "language_info": {
   "name": "julia",
   "version": "0.3.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}