Chapter 7 - Revisit Linear Regression With MCMC.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy as sp\n",
    "import scipy.stats\n",
    "import matplotlib.pyplot as plt\n",
    "import math as mt\n",
    "import scipy.special\n",
    "import seaborn as sns\n",
    "plt.style.use('fivethirtyeight')\n",
    "from statsmodels.graphics.tsaplots import plot_acf\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define some color in this cell.\n",
    "$$\n",
    "\\require{color}\n",
    "\\definecolor{red}{RGB}{240,5,5}\n",
    "\\definecolor{blue}{RGB}{5,5,240}\n",
    "\\definecolor{green}{RGB}{4,240,5}\n",
    "\\definecolor{black}{RGB}{0,0,0}\n",
    "\\definecolor{dsb}{RGB}{72, 61, 139}\n",
    "\\definecolor{Maroon}{RGB}{128,0,0}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "# <font face=\"gotham\" color=\"RebeccaPurple\"> Linear Regression Model Setting </font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The common matrix form of linear regression is\n",
    "$$\n",
    "y =    X \\beta+u\n",
    "$$\n",
    "where $u \\sim N( {0}, \\sigma^2 I)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The covariance matrix of disturbance term $u$ is \n",
    "$$\n",
    "\\begin{aligned}\n",
    "&\\operatorname{var}(u) \\equiv\\left[\\begin{array}{cccc}\n",
    "\\operatorname{var}\\left(u_{1}\\right) & \\operatorname{cov}\\left(u_{1}, u_{2}\\right) & \\ldots & \\operatorname{cov}\\left(u_{1}, u_{N}\\right) \\\\\n",
    "\\operatorname{cov}\\left(u_{1}, u_{2}\\right) & \\operatorname{var}\\left(u_{2}\\right) & \\ldots & . \\\\\n",
    "\\cdot & \\operatorname{cov}\\left(u_{2}, u_{3}\\right) & \\ldots & . \\\\\n",
    "\\cdot & \\cdot & \\ldots & \\operatorname{cov}\\left(u_{N-1}, u_{N}\\right) \\\\\n",
    "\\operatorname{cov}\\left(u_{1}, u_{N}\\right) & . & \\ldots & \\operatorname{var}\\left(u_{N}\\right)\n",
    "\\end{array}\\right]=\\left[\\begin{array}{cccc}\n",
    "h^{-1} & 0 & \\ldots & 0 \\\\\n",
    "0 & h^{-1} & \\ldots & . \\\\\n",
    ". & . & \\ldots & . \\\\\n",
    ". & . & \\cdots & 0 \\\\\n",
    "0 & . & . & h^{-1}\n",
    "\\end{array}\\right]\n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some derivation usually make use of $h = 1/\\sigma^2$, which is the **precision**. The diagonal form of covariance matrix represents two assumptions: **no serial correlation** and **homoscedasticity**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Also with the assumption that $X$ is exogenous, we can construct the joint probability density \n",
    "$$\n",
    "P(y, X \\mid \\beta, h)=P(y \\mid X, \\beta, h) P(X)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However because $X$ does not depend $\\beta$ and $h$, we narrow our interest onto\n",
    "$$\n",
    "P(y \\mid X, \\beta, h)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recall that multivariate normal distribution takes the form"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "f(X)=(2 \\pi)^{-N / 2}|\\Sigma|^{-1 / 2} \\exp \\left(-\\frac{1}{2}(X-\\mu)^T \\Sigma^{-1}(X-\\mu)\\right) \\text { for } \\sigma>0\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where $|\\Sigma|$ is the determinant of covariance matrix, $\\Sigma^{-1}$ is the inverse of covariance matrix."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font face=\"gotham\" color=\"DodgerBlue\"> Normal-Gamma Conjugacy</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the linear regression context, the determinant of covariance matrix is\n",
    "$$\n",
    "|\\text{var}(u)|^{-1/2}=\\left(\\prod_{i=1}^Nh^{-1}\\right)^{-1/2} = (h^{-N})^{-1/2}=h^{N/2}=|\\Sigma|^{-1/2}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The inverse matrix of covariance matrix\n",
    "$$\n",
    "\\text{var}(u)^{-1} = (h^{-1}I)^{-1}=hI=\\Sigma^{-1}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## <font face=\"gotham\" color=\"DodgerBlue\"> Likelihood Function </font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div style=\"background-color:DodgerBlue; color:FloralWhite; padding:30px;\">\n",
    "With all previous setting, the the likelihood function $P(y \\mid X, \\beta, h)$ is simplified to\n",
    "$$\n",
    "P(y \\mid  X, \\beta, h)=(2 \\pi)^{-\\frac{N}{2}} h^{\\frac{N}{2}} \\exp \\left(-\\frac{h}{2}(y-\\mathrm{X} \\beta)^T(y-\\mathrm{X} \\beta)\\right)\n",
    "$$\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <font face=\"gotham\" color=\"DodgerBlue\"> Kernel Decomposition </font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{aligned}\n",
    "( {y}- {X}  {\\beta})^{T}( {y}- {X}  {\\beta})=&(\\overbrace{ {y}- {X} \\hat{ {\\beta}}}^{A}+\\overbrace{{X} \\hat{{\\beta}}- {X}  {\\beta}}^{B})^{T} (\\overbrace{{y}- {X} \\hat{ {\\beta}}}^{A}+\\overbrace{ {X} \\hat{ {\\beta}}- {X}  {\\beta}}^{B}) \\\\\n",
    "=& A^{T} A+B^{T} B+ A^{T} B+B^TA\\\\\n",
    "=& A^{T} A+B^{T} B+2 A^{T} B \\\\\n",
    "&\\text{(We use the fact }A^TB=B^TA \\text{, because both terms are scalar, transposition is itself)}\\\\\n",
    "=&A^{T} A+B^{T} B-2\\overbrace{( {y}- {X} \\hat{ {\\beta}})^{T}( {X} \\hat{ {\\beta}}- {X}  {\\beta})}^{0} \\\\\n",
    "=&({y}-{X} \\hat{{\\beta}})^{T}({y}-{X} \\hat{{\\beta}})+(\\hat{\\beta}^TX^T-\\beta^TX^T)(X\\hat{\\beta}-X\\beta) \\\\\n",
    "=&({y}-{X} \\hat{{\\beta}})^{T}({y}-{X} \\hat{{\\beta}})+(\\hat{{\\beta}}-{\\beta})^{T} {X}^{T} {X}(\\hat{{\\beta}}-{\\beta}) \n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div style=\"background-color:DodgerBlue; color:FloralWhite; padding:30px;\">\n",
    "The kernel decomposition of likelihood function is\n",
    "$$\n",
    "( {y}- {X}  {\\beta})^{T}( {y}- {X}  {\\beta})=({y}-{X} \\hat{{\\beta}})^{T}({y}-{X} \\hat{{\\beta}})+(\\hat{{\\beta}}-{\\beta})^{T} {X}^{T} {X}(\\hat{{\\beta}}-{\\beta}) \n",
    "$$\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plug in back to likelihood function\n",
    "$$\n",
    "P(y \\mid  X, \\beta, h)=(2 \\pi)^{-\\frac{N}{2}} h^{\\frac{N}{2}} \\exp \\left[-\\frac{h}{2}\\left(({y}-{X} \\hat{{\\beta}})^{T}({y}-{X} \\hat{{\\beta}})+(\\beta-\\hat{\\beta})^T {X}^T {X}(\\beta-\\hat{\\beta})\\right)\\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div style=\"background-color:DodgerBlue; color:FloralWhite; padding:30px;\">\n",
    "Separate the exponential part\n",
    "$$\n",
    "p(y \\mid \\mathrm{X}, \\beta, h)=(2 \\pi)^{-\\frac{N}{2}} h^{\\frac{N}{2}} \\exp \\left[-\\frac{h}{2}\\left(\\beta-\\widehat{\\beta}\\right)^T\\mathrm{X}^T \\mathrm{X}\\left(\\beta-\\widehat{\\beta}\\right)\\right] \\exp \\left[-\\frac{h}{2}({y}-{X} \\hat{{\\beta}})^{T}({y}-{X} \\hat{{\\beta}})\\right]\n",
    "$$\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This form of likelihood function suggest a **natural conjugate prior** which has the same function form as likelihood and also yields a posterior within the same class of distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font face=\"gotham\" color=\"DodgerBlue\"> Prior </font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our prior should be a joint distribution $P(\\beta, h)$, however in order to transform into Normal-Gamma distribution, it proves convenient to write \n",
    "$$\n",
    "P(\\beta, h)=P(\\beta \\mid h) P(h)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And it is the **Normal-Gamma distribution** we are talking about, the right-hand side has the follow distribution\n",
    "$$\n",
    "\\begin{gathered}\n",
    "\\beta \\mid h \\sim N\\left(\\mu, (h V)^{-1}\\right) \\\\\n",
    "h \\sim \\operatorname{Gamma}(m, v)\n",
    "\\end{gathered}\n",
    "$$\n",
    "The advantage of this distribution is that $\\beta$ is on real number's domain and $h$ is on positive number's domain.\n",
    "\n",
    "Recall the function form of NG distribution.\n",
    "$$\n",
    "f(X, h)=(2 \\pi)^{-\\frac{N}{2}}(h)^{-\\frac{N}{2}}|\\Sigma|^{-1 / 2} \\exp \\left(-\\frac{h}{2}(X-\\mu)^T \\Sigma^{-1}(X-\\mu)\\right) \\frac{1}{\\left(\\frac{2 m}{v}\\right)^{v / 2}} \\frac{1}{\\Gamma\\left(\\frac{v}{2}\\right)} h^{\\frac{v-2}{2}} \\exp \\left[-\\frac{h v}{2 m}\\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div style=\"background-color:DodgerBlue; color:FloralWhite; padding:30px;\">\n",
    "Prior for normal part $P(\\beta\\mid h)$ \n",
    "$$\n",
    "\\begin{gathered}\n",
    "P(\\beta \\mid h)=(2 \\pi)^{-\\frac{k}{2}} h^{\\frac{k}{2}}|V|^{\\frac{1}{2}} \\exp \\left(-\\frac{h}{2}(\\beta-\\mu)^T V(\\beta-\\mu)\\right)\n",
    "\\end{gathered}\n",
    "$$\n",
    "where $k$ is the number of parameters in $\\beta$. We use the rule of determinant $|V^{-1}|=|V|^{-1}$.\n",
    "<br><br>\n",
    "Prior for gamma part $P(h)$ \n",
    "$$\n",
    "P(h)=\\frac{1}{\\left(\\frac{2 m}{v}\\right)^{v / 2}} \\frac{1}{\\Gamma\\left(\\frac{v}{2}\\right)} h^{\\frac{v-2}{2}} \\exp \\left[-\\frac{h v}{2 m}\\right]\n",
    "$$\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font face=\"gotham\" color=\"DodgerBlue\"> Posterior </font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The posterior is formulated by Bayes' Theorem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "P(\\beta, h \\mid y, X)\\propto P(y \\mid  X, \\beta, h)P(\\beta\\mid h) P(h)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "P(\\beta, h \\mid Y, {X}) \\propto \\ h^{\\frac{k}{2}} \\exp \\left[-\\frac{h}{2}\\left[\\left(\\beta-\\hat{\\beta}\\right)^T{X}^{T} {X}\\left(\\beta-\\hat{\\beta}\\right)+(\\beta-\\mu)^T V(\\beta-\\mu)\\right]\\right]  h^{\\frac{N+v-2}{2}} \\exp \\left[-\\frac{h}{2}\\left(({y}-{X} \\hat{{\\beta}})^{T}({y}-{X} \\hat{{\\beta}})+\\frac{v}{m}\\right)\\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Expand both terms, again the $2 \\hat{ {\\beta}}^T  {X}^{T}  {X}  {\\beta}$ and $2  {\\mu}^{T}  {\\Lambda}  {\\beta}$ exist because cross terms are scalars.\n",
    "\\begin{align}\n",
    "({\\beta}-\\hat{\\beta})^{T} {X}^{T} {X}({\\beta}-\\hat{\\beta})+\\left({\\beta}-{\\mu}\\right)^{T} {V}\\left({\\beta}-{\\mu}\\right)&= {\\beta}^{T}  {X}^{T}  {X}  {\\beta}+\\hat{ {\\beta}}^T  {X}^{T}  {X} \\hat{ {\\beta}}-2 \\hat{ {\\beta}}^T  {X}^{T}  {X}  {\\beta}+ {\\beta}^{T}  V  {\\beta}+ {\\mu}^{T}  V {\\mu}-2  {\\mu}^{T} V  {\\beta}\\\\\n",
    "&=\\color{red}{\\beta}^{T}  {X}^{T}  {X}  {\\beta}\\color{black} + \\color{red}{\\beta}^{T}  {V}  {\\beta} \\color{black}-\\color{blue}2 \\hat{ {\\beta}}^T  {X}^{T}  {X}  {\\beta}-\\color{blue}2  {\\mu}^{T}  V  {\\beta}+\\color{black} {\\mu}^{T} V  {\\mu}+\\hat{ {\\beta}}^T  {X}^{T}  {X} \\hat{ {\\beta}}\\\\\n",
    "&=\\color{red}\\beta^T\\underbrace{(X^TX+V)}_{W}\\beta\\color{black}-\\color{blue}2\\underbrace{(\\hat{\\beta}^TX^TX+\\mu^TV)}_{w^T}\\beta+\\color{black} {\\mu}^{T}  {V}  {\\mu}+\\hat{ {\\beta}}^T  {X}^{T}  {X} \\hat{ {\\beta}}\\\\\n",
    "&=\\underbrace{\\color{red}\\beta^TW\\beta\\color{black}-\\color{blue}2w^T\\beta}_{\\text{needs completing the square}}+\\color{black} {\\mu}^{T}  {V}  {\\mu}+\\hat{ {\\beta}}^T  {X}^{T}  {X} \\hat{ {\\beta}}\\\\\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where the variable are replace for completing the square\n",
    "$$\n",
    "W = (X^TX+V)\\\\\n",
    "w^T=\\hat{\\beta}^TX^TX+\\mu^TV\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Details of completing the square\n",
    "$$\n",
    "\\begin{align}\n",
    "\\color{red}\\beta^TW\\beta\\color{black}-\\color{blue}2w^T\\beta\\color{black}+\\overbrace{\\color{green}w^TW^{-1}w-w^TW^{-1}w}^{0}\n",
    "&=\\overbrace{\\color{red}\\beta^TW\\beta\\color{black}-\\color{blue}w^T\\beta\\color{black}-\\color{blue}\\beta^Tw\\color{black}+\\color{green}w^TW^{-1}w}^{\\text{complete the square}}-\\color{green}w^TW^{-1}w\\\\\n",
    "&=\\color{Maroon}(\\beta^TW-w^T)(\\beta-W^{-1}w)\\color{black}-\\color{green}w^TW^{-1}w\\\\\n",
    "&=\\color{Maroon}(\\beta^T-w^TW^{-1})W(\\beta-W^{-1}w)\\color{black}-\\color{green}w^TW^{-1}w\\\\\n",
    "&=\\color{Maroon}(\\beta-W^{-1}w)^TW(\\beta-W^{-1}w)\\color{black}-\\color{green}w^TW^{-1}w\n",
    "\\end{align}\n",
    "$$\n",
    "We used the fact that $(M^{-1})^T=M^{-1}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div style=\"background-color:DodgerBlue; color:FloralWhite; padding:30px;\">\n",
    "So this is the new expression we are working on\n",
    "$$\n",
    "(\\hat{{\\beta}}-{\\beta})^{T} {X}^{T} {X}(\\hat{{\\beta}}-{\\beta})+\\left({\\beta}-{\\mu}\\right)^{T} {V}\\left({\\beta}-{\\mu}\\right)=\\color{Maroon}\\underbrace{\\color{Maroon}(\\beta-W^{-1}w)^TW(\\beta-W^{-1}w)}_{\\text{Posterior normal kernal}}\\color{black}-\\color{green}w^TW^{-1}w+\\color{black} {\\mu}^{T}  {V}  {\\mu}+\\hat{ {\\beta}}^T  {X}^{T}  {X} \\hat{ {\\beta}}\n",
    "$$\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the mean in the posterior normal kernel, $N$ subscript in $\\mu_N^{\\ }$ means it is from the normal part\n",
    "$$\n",
    "\\begin{align}\n",
    "\\mu_{N}^{\\ } = W^{-1}w &= W^{-1}(X^TX\\hat{\\beta}+V^T\\mu)\\\\\n",
    "&=W^{-1}(X^TX(X^TX)^{-1}X^Xy+V^T\\mu)\\\\\n",
    "&=W^{-1}\\underbrace{(X^Ty+V^T\\mu)}_{w}\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div style=\"background-color:DodgerBlue; color:FloralWhite; padding:30px;\">\n",
    "So this is the new expression we are working on\n",
    "$$\n",
    "(\\hat{{\\beta}}-{\\beta})^{T} {X}^{T} {X}(\\hat{{\\beta}}-{\\beta})+\\left({\\beta}-{\\mu}\\right)^{T} {V}\\left({\\beta}-{\\mu}\\right)=\\color{Maroon}\\underbrace{\\color{Maroon}(\\beta-\\mu_{N}^{\\ })^TW(\\beta-\\mu_{N}^{\\ })}_{\\text{Posterior normal kernel}}\\color{black}-\\color{green}w^TW^{-1}w+\\color{black} {\\mu}^{T}  {V}  {\\mu}+\\hat{ {\\beta}}^T  {X}^{T}  {X} \\hat{ {\\beta}}\n",
    "$$\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div style=\"background-color:DodgerBlue; color:FloralWhite; padding:30px;\">\n",
    "So the undecomposed the likelihood kernel with $P(\\beta\\mid h)$ kernel is\n",
    "$$\n",
    "\\underbrace{({y}-{X} {\\beta})^{T}({y}-{X} {\\beta})}_{({y}-{X} \\hat{{\\beta}})^{T}({y}-{X} \\hat{{\\beta}})+(\\hat{{\\beta}}-{\\beta})^{T} {X}^{T} {X}(\\hat{{\\beta}}-{\\beta}) }+\\left({\\beta}-{\\mu}\\right)^{T} {V}\\left({\\beta}-{\\mu}\\right)=({y}-{X} \\hat{\\beta})^{T}({y}-{X} \\hat{\\beta})+\\color{Maroon}\\underbrace{\\color{Maroon}(\\beta-\\mu_{N}^{\\ })^TW(\\beta-\\mu_{N}^{\\ })}_{\\text{Posterior normal kernel}}\\color{black}-\\color{green}w^TW^{-1}w+\\color{black} {\\mu}^{T}  {V}  {\\mu}+\\hat{ {\\beta}}^T  {X}^{T}  {X} \\hat{ {\\beta}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Simplify each on the right-hand side"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{aligned}\n",
    "( {y}- {X} \\hat{ {\\beta}})^{T}( {y}- {X} \\hat{ {\\beta}})&= {y}^{T}  {y}+\\hat{ {\\beta}}^{T}  {X}^{T}  {X} \\hat{ {\\beta}} -2  {y}^{T}  {X} \\hat{ {\\beta}} \\\\\n",
    "&= {y}^{T}  {y}+\\hat{ {\\beta}}^{T}  {X}^{T}  {X} (\\overbrace{X^TX)^{-1}X^Ty}^{\\hat{\\beta}} -2  {y}^{T}  {X} \\hat{ {\\beta}}  \\\\\n",
    "&= {y}^{T}  {y}+\\hat{ {\\beta}}^{T}  X^Ty -2  {y}^{T}  {X} \\hat{ {\\beta}}  \\\\\n",
    "&= {y}^{T}  {y}+(\\overbrace{(X^TX)^{-1}X^Ty}^{\\hat{\\beta}})^{T}  X^Ty -2  {y}^{T}  {X} \\hat{ {\\beta}}  \\\\\n",
    "&= y^Ty+y^TX(X^TX)^{-1}X^Ty - 2y^TX\\hat{\\beta}\\\\\n",
    "&=  {y}^{T}  {y}+{y}^{T}  {X} \\hat{ {\\beta}}-2{y}^{T}  {X} \\hat{ {\\beta}}\\\\\n",
    "&={y}^{T}  {y}-{y}^{T}  {X} \\hat{ {\\beta}}\\\\\n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align}\n",
    "w^TW^{-1}w &= (X^Ty+V^T\\mu)^T\\underbrace{W^{-1}(X^Ty+V^T\\mu)}_{\\mu_{N}^{\\ }}\\\\\n",
    "&= (X^Ty+\\Lambda^T\\mu)^T{\\mu_{N}^{\\ }}\\\\\n",
    "&= \\mu_N^TW \\mu_{N}^{\\ }\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align}\n",
    "\\hat{ {\\beta}}^T  {X}^{T}  {X} \\hat{ {\\beta}}&=\\hat{ {\\beta}}^T  {X}^{T}  {X}(X^TX)^{-1}X^Ty\\\\\n",
    "&=\\hat{ {\\beta}}^T  X^Ty\\\\\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div style=\"background-color:DodgerBlue; color:FloralWhite; padding:30px;\">\n",
    "Join all the terms of undecomposed likelihood kernel plus $P(\\beta\\mid h)$\n",
    "$$\n",
    "\\begin{align}\n",
    "\\underbrace{({y}-{X} {\\beta})^{T}({y}-{X} {\\beta})}_{({y}-{X} \\hat{{\\beta}})^{T}({y}-{X} \\hat{{\\beta}})+(\\hat{{\\beta}}-{\\beta})^{T} {X}^{T} {X}(\\hat{{\\beta}}-{\\beta}) }+\\left({\\beta}-{\\mu}\\right)^{T} {V}\\left({\\beta}-{\\mu}\\right)&=y^T{y}-{y}^{T}  {X} \\hat{ {\\beta}}+\\color{Maroon}\\underbrace{\\color{Maroon}(\\beta-\\mu_{N}^{\\ })^TW(\\beta-\\mu_{N}^{\\ })}_{\\text{Posterior normal kernal} } \\color{black}-\\mu_N^T W \\mu_{N}^{\\ }+\\color{black} {\\mu}^{T}  {V}  {\\mu}+\\hat{ {\\beta}}^T  {X}^{T}  {X} \\hat{ {\\beta}}\\\\\n",
    "&= \\color{Maroon}\\underbrace{\\color{Maroon}(\\beta-\\mu_{N}^{\\ })^TW(\\beta-\\mu_{N}^{\\ })}_{\\text{Posterior normal kernal}}\\color{black}+{y}^{T}  {y} -\\mu_N^T W \\mu_{N}^{\\ }+\\color{black} {\\mu}^{T}  {V}  {\\mu}\n",
    "\\end{align}\n",
    "$$\n",
    "where \n",
    "$$\n",
    "W=X^TX+V\\\\\n",
    "\\mu_{N}^{\\ }=W^{-1}\\underbrace{(X^Ty+V^T\\mu)}_{w}\n",
    "$$\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align}\n",
    "P(\\beta, h \\mid {X}, y) \\propto&(2 \\pi)^{-\\frac{N}{2}} h^{\\frac{N}{2}} \\exp{\\left(-\\frac{h}{2}(\\beta-\\mu_{N}^{\\ })^TW(\\beta-\\mu_{N}^{\\ })\\right)}\\\\\n",
    "&\\times(2 \\pi)^{-\\frac{k}{2}} h^{\\frac{k}{2}}|V|^{\\frac{1}{2}}\\exp \\left(-\\frac{h}{2}({y}^{T}  {y} -\\mu_N^T W \\mu_{N}^{\\ }+\\color{black} {\\mu}^{T}  {V}  {\\mu})\\right)\\\\\n",
    "&\\times\\frac{1}{\\beta^{\\alpha}}\\frac{1}{\\Gamma(\\alpha)} h^{\\alpha-1}\\exp{\\left(-\\frac{h}{\\beta}\\right)}\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Omit the terms that do not depend on $\\beta$ or $h$\n",
    "$$\n",
    "\\begin{align}\n",
    "P(\\beta, h \\mid {X}, y) \\propto&(2 \\pi)^{-\\frac{N}{2}} h^{\\frac{N}{2}} \\exp{\\left(-\\frac{h}{2}(\\beta-\\mu_{N}^{\\ })^TW(\\beta-\\mu_{N}^{\\ })\\right)}\\\\\n",
    "&\\times h^{\\frac{k}{2}}|V|^{\\frac{1}{2}}\\exp \\left({y}^{T}  {y} -\\mu_N^T W \\mu_{N}^{\\ }+\\color{black} {\\mu}^{T}  {V}  {\\mu}\\right)\\\\\n",
    "&\\times\\frac{1}{\\beta^{\\alpha}}\\frac{1}{\\Gamma(\\alpha)} h^{\\alpha-1}\\exp{\\left(-\\frac{h}{\\beta}\\right)}\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Omit the terms that do not depend on $\\beta$ or $h$\n",
    "$$\n",
    "\\begin{align}\n",
    "P(\\beta, h \\mid {X}, y) \\propto&(2 \\pi)^{-\\frac{N}{2}} h^{\\frac{N}{2}} \\exp{\\left(-\\frac{h}{2}(\\beta-\\mu_{N}^{\\ })^TW(\\beta-\\mu_{N}^{\\ })\\right)}\\\\\n",
    "&\\times h^{\\alpha+\\frac{k}{2}-1}|V|^{\\frac{1}{2}}\\exp \\left(-\\frac{h}{\\beta}+{y}^{T}  {y} -\\mu_N^T W \\mu_{N}^{\\ }+\\color{black} {\\mu}^{T}  {V}  {\\mu}\\right)\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First ignore the constant parts such as $2\\pi$ and gamma function, then combine $h$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align}\n",
    "P(\\beta, h \\mid Y, \\mathrm{X}) &\\propto h^{\\frac{k}{2}} \\exp \\left[-\\frac{h}{2}\\left(\\beta-\\widehat{\\beta}\\right)^T\\mathrm{X}^T \\mathrm{X}\\left(\\beta-\\widehat{\\beta}\\right)\\right] \\exp \\left[-\\frac{h}{2}\\hat{u}^T \\hat{u}\\right]\\\\\n",
    "&\\times \\exp \\left[-\\frac{h}{2}(\\beta-\\mu)^{\\prime} V^{-1}(\\beta-\\mu)\\right] h^{\\frac{N+v-2}{2}} \\exp \\left[-\\frac{h v}{2 m}\\right]\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Join exponential terms\n",
    "$$\n",
    "\\begin{align}\n",
    "P(\\beta, h \\mid Y, \\mathrm{X}) &\\propto h^{\\frac{k}{2}} \\exp \\left[-\\frac{h}{2}\\left(\\beta-\\widehat{\\beta}\\right)^T\\mathrm{X}^T \\mathrm{X}\\left(\\beta-\\widehat{\\beta}\\right)-\\frac{h}{2}(\\beta-\\mu)^{\\prime} V^{-1}(\\beta-\\mu)\\right] h^{\\frac{N+v-2}{2}} \\exp \\left[-\\frac{h}{2}\\hat{u}^T \\hat{u}-\\frac{h v}{2 m}\\right]\\\\\n",
    "&\\propto\\ h^{\\frac{k}{2}} \\exp \\left[-\\frac{h}{2}\\left[\\left(\\beta-\\widehat{\\beta}\\right)^T\\mathrm{X}^{\\prime} \\mathrm{X}\\left(\\beta-\\widehat{\\beta}\\right)+(\\beta-\\mu)^T V^{-1}(\\beta-\\mu)\\right]\\right]  h^{\\frac{N+v-2}{2}} \\exp \\left[-\\frac{h}{2}\\left(\\hat{u}^T \\hat{u}+\\frac{v}{m}\\right)\\right]\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font face=\"gotham\" color=\"orange\"> Normal-Inverse Gamma Conjugacy</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recall multivariate normal distribution density function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "f(X)=(2 \\pi)^{-\\frac{N}{2}}|\\Sigma|^{-\\frac{1}{2}} \\exp \\left(-\\frac{1}{2}(X-\\mu)^T \\Sigma^{-1}(X-\\mu)\\right) \\text { for } \\sigma>0\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We assume likelihood ${y} \\sim {\\mathcal{N}}\\left({ {X}\\beta}, \\sigma^{2} {I}\\right)$ follows exactly the same form as density function above"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "p\\left({y} \\mid {X}, {\\beta}, \\sigma^{2}\\right)=(2 \\pi)^{-\\frac{N}{2}}\\left|\\sigma^{2} {I}\\right|^{-\\frac{1}{2}} \\exp \\left(-\\frac{1}{2}\\left({y}-{ {X}\\beta}\\right)^{T}\\left(\\sigma^{2} {I}\\right)^{-1}\\left({y}-{ {X}\\beta}\\right)\\right)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where $N$ is the number of observation in the data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use determinants rules\n",
    "$$\n",
    "\\left|\\sigma^{2} {I}\\right|=\\sigma^{2N}\\\\\n",
    "\\left(\\sigma^{2} {I}\\right)^{-1}=\\sigma^{-2}I\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div style=\"background-color:LightCoral; color:DarkBlue; padding:30px;\">\n",
    "    Likelihood simplifed to\n",
    "$$\n",
    "p\\left({y} \\mid {X}, {\\beta}, \\sigma^{2}\\right)=\\left(2 \\pi \\sigma^{2}\\right)^{-N / 2} \\exp \\left(-\\frac{1}{2 \\sigma^{2}}\\left({y}- {X}{\\beta}\\right)^{T}\\left( {y}- {X}{\\beta}\\right)\\right)\n",
    "$$\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Prior is usually decomposed $p\\left({\\beta}, \\sigma^{2}\\right)=p\\left({\\beta} \\mid \\sigma^{2}\\right) p\\left(\\sigma^{2}\\right)$ where ${\\beta} \\mid \\sigma^{2} \\sim {N}\\left({0}, \\sigma^{2} {\\Lambda}^{-1}\\right)$ and $\\sigma^{2} \\sim \\operatorname{InvGamma}\\left(\\alpha, \\beta\\right)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "p\\left( {\\beta} \\mid \\sigma^{2}\\right)=(2 \\pi)^{-\\frac{k}{2}}\\left|\\sigma^{2} {\\Lambda}^{-1}\\right|^{-\\frac{1}{2}} \\exp \\left(-\\frac{1}{2}\\left({\\beta}-{\\mu}\\right)^T\\left(\\sigma^2 {\\Lambda}^{-1}\\right)^{-1}\\left({\\beta}-{\\mu}\\right)\\right)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where $k$ is the number of $\\beta$ parameters in the model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use determinants rules\n",
    "$$\n",
    "|\\sigma^2\\Lambda^{-1}|^{-\\frac{1}{2}}=(\\sigma^{2k}|\\Lambda^{-1}|)^{-\\frac{1}{2}}=(\\sigma^{2k}|\\Lambda|^{-1})^{-\\frac{1}{2}}=\\sigma^{-k}|\\Lambda|^{\\frac{1}{2}}\\\\\n",
    "$$\n",
    "then\n",
    "$$\n",
    "(2 \\pi)^{-\\frac{k}{2}}\\left|\\sigma^{2} {\\Lambda}^{-1}\\right|^{-\\frac{1}{2}}=(2\\pi\\sigma^2)^{-\\frac{k}{2}}|\\Lambda|^\\frac{1}{2}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div style=\"background-color:LightCoral; color:DarkBlue; padding:30px;\">\n",
    "First part of prior $P\\left( {\\beta} \\mid \\sigma^{2}\\right)$ simplified to\n",
    "$$\n",
    "p\\left( {\\beta} \\mid \\sigma^{2}\\right)=\\left(2 \\pi \\sigma^{2}\\right)^{-\\frac{k}{2}}\\left|{\\Lambda}\\right|^{\\frac{1}{2}} \\exp \\left(-\\frac{1}{2 \\sigma^{2}}\\left({\\beta}-{\\mu}\\right)^{T} {\\Lambda}\\left({\\beta}-{\\mu}\\right)\\right) \n",
    "$$\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recall the inverse gamma distribution is\n",
    "\n",
    "$$\n",
    "f(x ; \\alpha, \\beta)=\\frac{\\beta^{\\alpha}}{\\Gamma(\\alpha)}(1 / x)^{\\alpha+1} \\exp (-\\beta / x)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div style=\"background-color:LightCoral; color:DarkBlue; padding:30px;\">\n",
    "following the function form, the second part of prior $P(\\sigma^2)$ is\n",
    "$$\n",
    "P\\left(\\sigma^{2}\\right)=\\frac{\\beta^{\\alpha}}{\\Gamma\\left(\\alpha\\right)}\\left(\\sigma^{2}\\right)^{-\\left(\\alpha+1\\right)} \\exp \\left(-\\beta / \\sigma^{2}\\right)\n",
    "$$\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font face=\"gotham\" color=\"orange\"> Kernel Decomposition </font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To move forward, there will be horrible amount of linear algebraic manipulation, the first is kernel decomposition of likelihood function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{aligned}\n",
    "( {y}- {X}  {\\beta})^{T}( {y}- {X}  {\\beta})=&(\\overbrace{ {y}- {X} \\hat{ {\\beta}}}^{A}+\\overbrace{{X} \\hat{{\\beta}}- {X}  {\\beta}}^{B})^{T} (\\overbrace{{y}- {X} \\hat{ {\\beta}}}^{A}+\\overbrace{ {X} \\hat{ {\\beta}}- {X}  {\\beta}}^{B}) \\\\\n",
    "=& A^{T} A+B^{T} B+ A^{T} B+B^TA\\\\\n",
    "=& A^{T} A+B^{T} B+2 A^{T} B \\\\\n",
    "&\\text{(We use the fact }A^TB=B^TA \\text{, because both terms are scalar, transposition is itself)}\\\\\n",
    "=&A^{T} A+B^{T} B-2\\overbrace{( {y}- {X} \\hat{ {\\beta}})^{T}( {X} \\hat{ {\\beta}}- {X}  {\\beta})}^{0} \\\\\n",
    "=&({y}-{X} \\hat{{\\beta}})^{T}({y}-{X} \\hat{{\\beta}})+(\\hat{\\beta}^TX^T-\\beta^TX^T)(X\\hat{\\beta}-X\\beta) \\\\\n",
    "=&({y}-{X} \\hat{{\\beta}})^{T}({y}-{X} \\hat{{\\beta}})+(\\hat{{\\beta}}-{\\beta})^{T} {X}^{T} {X}(\\hat{{\\beta}}-{\\beta}) \n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div style=\"background-color:LightCoral; color:DarkBlue; padding:30px;\">\n",
    "The kernel decomposition of likelihood function is\n",
    "$$\n",
    "( {y}- {X}  {\\beta})^{T}( {y}- {X}  {\\beta})=({y}-{X} \\hat{{\\beta}})^{T}({y}-{X} \\hat{{\\beta}})+(\\hat{{\\beta}}-{\\beta})^{T} {X}^{T} {X}(\\hat{{\\beta}}-{\\beta}) \n",
    "$$\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The quadratic form $(\\hat{{\\beta}}-{\\beta})^{T} {X}^{T} {X}(\\hat{{\\beta}}-{\\beta})$ is created be combined with the kernal of prior $P\\left( {\\beta} \\mid \\sigma^{2}\\right)$ that is $\\left({\\beta}-{\\mu}\\right)^{T} {\\Lambda}\\left({\\beta}-{\\mu}\\right)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Joining second decomposed terms in likelihood and prior $P(\\beta\\mid\\sigma^2)$ will end up as an addition of them due the feature of exponential term\n",
    "$$\n",
    "(\\hat{{\\beta}}-{\\beta})^{T} {X}^{T} {X}(\\hat{{\\beta}}-{\\beta})+\\left({\\beta}-{\\mu}\\right)^{T} {\\Lambda}\\left({\\beta}-{\\mu}\\right)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Expand both terms, again the $2 \\hat{ {\\beta}}^T  {X}^{T}  {X}  {\\beta}$ and $2  {\\mu}^{T}  {\\Lambda}  {\\beta}$ exist because cross terms are scalars.\n",
    "\\begin{align}\n",
    "(\\hat{{\\beta}}-{\\beta})^{T} {X}^{T} {X}(\\hat{{\\beta}}-{\\beta})+\\left({\\beta}-{\\mu}\\right)^{T} {\\Lambda}\\left({\\beta}-{\\mu}\\right)&= {\\beta}^{T}  {X}^{T}  {X}  {\\beta}+\\hat{ {\\beta}}^T  {X}^{T}  {X} \\hat{ {\\beta}}-2 \\hat{ {\\beta}}^T  {X}^{T}  {X}  {\\beta}+ {\\beta}^{T}  {\\Lambda}  {\\beta}+ {\\mu}^{T}  {\\Lambda}  {\\mu}-2  {\\mu}^{T}  {\\Lambda}  {\\beta}\\\\\n",
    "&=\\color{red}{\\beta}^{T}  {X}^{T}  {X}  {\\beta}\\color{black} + \\color{red}{\\beta}^{T}  {\\Lambda}  {\\beta} \\color{black}-\\color{blue}2 \\hat{ {\\beta}}^T  {X}^{T}  {X}  {\\beta}-\\color{blue}2  {\\mu}^{T}  {\\Lambda}  {\\beta}+\\color{black} {\\mu}^{T}  {\\Lambda}  {\\mu}+\\hat{ {\\beta}}^T  {X}^{T}  {X} \\hat{ {\\beta}}\\\\\n",
    "&=\\color{red}\\beta^T\\underbrace{(X^TX+\\Lambda)}_{M}\\beta\\color{black}-\\color{blue}2\\underbrace{(\\hat{\\beta}^TX^TX+\\mu^T\\Lambda)}_{m^T}\\beta+\\color{black} {\\mu}^{T}  {\\Lambda}  {\\mu}+\\hat{ {\\beta}}^T  {X}^{T}  {X} \\hat{ {\\beta}}\\\\\n",
    "&=\\underbrace{\\color{red}\\beta^TM\\beta\\color{black}-\\color{blue}2m^T\\beta}_{\\text{needs completing the square}}+\\color{black} {\\mu}^{T}  {\\Lambda}  {\\mu}+\\hat{ {\\beta}}^T  {X}^{T}  {X} \\hat{ {\\beta}}\\\\\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have defined \n",
    "$$\n",
    "M=X^TX+\\Lambda\\\\\n",
    "m^T = \\hat{\\beta}^TX^TX+\\mu^T\\Lambda\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Continue to complete the square of colored term\n",
    "$$\n",
    "\\begin{align}\n",
    "\\color{red}\\beta^TM\\beta\\color{black}-\\color{blue}2m^T\\beta\\color{black}+\\overbrace{\\color{green}m^TM^{-1}m-m^TM^{-1}m}^{0}\n",
    "&=\\overbrace{\\color{red}\\beta^TM\\beta\\color{black}-\\color{blue}m^T\\beta\\color{black}-\\color{blue}\\beta^Tm\\color{black}+\\color{green}m^TM^{-1}m}^{\\text{complete the square}}-\\color{green}m^TM^{-1}m\\\\\n",
    "&=\\color{Maroon}(\\beta^TM-m^T)(\\beta-M^{-1}m)\\color{black}-\\color{green}m^TM^{-1}m\\\\\n",
    "&=\\color{Maroon}(\\beta^T-m^TM^{-1})M(\\beta-M^{-1}m)\\color{black}-\\color{green}m^TM^{-1}m\\\\\n",
    "&=\\color{Maroon}(\\beta-M^{-1}m)^TM(\\beta-M^{-1}m)\\color{black}-\\color{green}m^TM^{-1}m\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Last step above made use of the fact $(M^{-1})^T=M^{-1}$, invoking the inverse matrix rule $(M^{-1})^T = (M^{T})^{-1}$ ."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div style=\"background-color:LightCoral; color:DarkBlue; padding:30px;\">\n",
    "So this is the new expression we are working on\n",
    "$$\n",
    "(\\hat{{\\beta}}-{\\beta})^{T} {X}^{T} {X}(\\hat{{\\beta}}-{\\beta})+\\left({\\beta}-{\\mu}\\right)^{T} {\\Lambda}\\left({\\beta}-{\\mu}\\right)=\\color{Maroon}\\underbrace{\\color{Maroon}(\\beta-M^{-1}m)^TM(\\beta-M^{-1}m)}_{\\text{Posterior normal kernal}}\\color{black}-\\underbrace{\\color{green}m^TM^{-1}m+\\color{black} {\\mu}^{T}  {\\Lambda}  {\\mu}+\\hat{ {\\beta}}^T  {X}^{T}  {X} \\hat{ {\\beta}}}_{\\text{shuffle into inverse-gamma part}}\n",
    "$$\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So far we have achieved the kernel for the normal distribution part in posterior. And the rest will be shuffle into the inverse-gamma part."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the mean in the posterior normal kernel, $N$ subscript in $\\mu_N^{\\ }$ shows it is from the normal part\n",
    "$$\n",
    "\\begin{align}\n",
    "\\mu_{N}^{\\ } = M^{-1}m &= M^{-1}(X^TX\\hat{\\beta}+\\Lambda^T\\mu)\\\\\n",
    "&=M^{-1}(X^TX(X^TX)^{-1}X^Xy+\\Lambda^T\\mu)\\\\\n",
    "&=M^{-1}\\underbrace{(X^Ty+\\Lambda^T\\mu)}_{m}\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have defined $m^T$ above, here we derived its transpose form i.e. $m$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div style=\"background-color:LightCoral; color:DarkBlue; padding:30px;\">\n",
    "So this is the new expression we are working on\n",
    "$$\n",
    "(\\hat{{\\beta}}-{\\beta})^{T} {X}^{T} {X}(\\hat{{\\beta}}-{\\beta})+\\left({\\beta}-{\\mu}\\right)^{T} {\\Lambda}\\left({\\beta}-{\\mu}\\right)=\\color{Maroon}\\underbrace{\\color{Maroon}(\\beta-\\mu_{N}^{\\ })^TM(\\beta-\\mu_{N}^{\\ })}_{\\text{Posterior normal kernel}}\\color{black}-\\color{green}m^TM^{-1}m+\\color{black} {\\mu}^{T}  {\\Lambda}  {\\mu}+\\hat{ {\\beta}}^T  {X}^{T}  {X} \\hat{ {\\beta}}\n",
    "$$\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div style=\"background-color:LightCoral; color:DarkBlue; padding:30px;\">\n",
    "To summarize, the undecomposed the likelihood kernel with $P(\\beta\\mid\\sigma^2)$ kernel is\n",
    "$$\n",
    "\\underbrace{({y}-{X} {\\beta})^{T}({y}-{X} {\\beta})}_{({y}-{X} \\hat{{\\beta}})^{T}({y}-{X} \\hat{{\\beta}})+(\\hat{{\\beta}}-{\\beta})^{T} {X}^{T} {X}(\\hat{{\\beta}}-{\\beta}) }+\\left({\\beta}-{\\mu}\\right)^{T} {\\Lambda}\\left({\\beta}-{\\mu}\\right)=({y}-{X} \\hat{\\beta})^{T}({y}-{X} \\hat{\\beta})+\\color{Maroon}\\underbrace{\\color{Maroon}(\\beta-\\mu_{N}^{\\ })^TM(\\beta-\\mu_{N}^{\\ })}_{\\text{Posterior normal kernel}}\\color{black}-\\color{green}m^TM^{-1}m+\\color{black} {\\mu}^{T}  {\\Lambda}  {\\mu}+\\hat{ {\\beta}}^T  {X}^{T}  {X} \\hat{ {\\beta}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Simplify each term on the right hand side"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{aligned}\n",
    "( {y}- {X} \\hat{ {\\beta}})^{T}( {y}- {X} \\hat{ {\\beta}})&= {y}^{T}  {y}+\\hat{ {\\beta}}^{T}  {X}^{T}  {X} \\hat{ {\\beta}} -2  {y}^{T}  {X} \\hat{ {\\beta}} \\\\\n",
    "&= {y}^{T}  {y}+\\hat{ {\\beta}}^{T}  {X}^{T}  {X} (\\overbrace{X^TX)^{-1}X^Ty}^{\\hat{\\beta}} -2  {y}^{T}  {X} \\hat{ {\\beta}}  \\\\\n",
    "&= {y}^{T}  {y}+\\hat{ {\\beta}}^{T}  X^Ty -2  {y}^{T}  {X} \\hat{ {\\beta}}  \\\\\n",
    "&= {y}^{T}  {y}+(\\overbrace{(X^TX)^{-1}X^Ty}^{\\hat{\\beta}})^{T}  X^Ty -2  {y}^{T}  {X} \\hat{ {\\beta}}  \\\\\n",
    "&= y^Ty+y^TX(X^TX)^{-1}X^Ty - 2y^TX\\hat{\\beta}\\\\\n",
    "&=  {y}^{T}  {y}+{y}^{T}  {X} \\hat{ {\\beta}}-2{y}^{T}  {X} \\hat{ {\\beta}}\\\\\n",
    "&={y}^{T}  {y}-{y}^{T}  {X} \\hat{ {\\beta}}\\\\\n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align}\n",
    "m^TM^{-1}m &= (X^Ty+\\Lambda_\\mu^T\\mu)^T\\underbrace{M^{-1}(X^Ty+\\Lambda^T\\mu)}_{\\mu_{N}^{\\ }}\\\\\n",
    "&= (X^Ty+\\Lambda^T\\mu)^T{\\mu_{N}^{\\ }}\\\\\n",
    "&= \\mu_N^T M \\mu_{N}^{\\ }\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align}\n",
    "\\hat{ {\\beta}}^T  {X}^{T}  {X} \\hat{ {\\beta}}&=\\hat{ {\\beta}}^T  {X}^{T}  {X}(X^TX)^{-1}X^Ty\\\\\n",
    "&=\\hat{ {\\beta}}^T  X^Ty\\\\\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Join them"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\underbrace{({y}-{X} {\\beta})^{T}({y}-{X} {\\beta})}_{({y}-{X} \\hat{{\\beta}})^{T}({y}-{X} \\hat{{\\beta}})+(\\hat{{\\beta}}-{\\beta})^{T} {X}^{T} {X}(\\hat{{\\beta}}-{\\beta}) }+\\left({\\beta}-{\\mu}\\right)^{T} {\\Lambda}\\left({\\beta}-{\\mu}\\right)=y^T{y}-{y}^{T}  {X} \\hat{ {\\beta}}+\\color{Maroon}\\underbrace{\\color{Maroon}(\\beta-\\mu_{N}^{\\ })^TM(\\beta-\\mu_{N}^{\\ })}_{\\text{Posterior normal kernal} } \\color{black}-\\mu_N^T M \\mu_{N}^{\\ }+\\color{black} {\\mu}^{T}  {\\Lambda}  {\\mu}+\\hat{ {\\beta}}^T  {X}^{T}  {X} \\hat{ {\\beta}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div style=\"background-color:LightCoral; color:DarkBlue; padding:30px;\">\n",
    "Final results of undecomposed likelihood kernel plus $P(\\beta\\mid \\sigma^2)$ kernel is\n",
    "$$\n",
    "\\begin{align}\n",
    "&(y-X\\beta)^T(y-X\\beta)+(\\beta-\\mu)^T\\Lambda(\\beta-\\mu)\\\\\n",
    "&=(y-X\\hat{\\beta})^T(y-X\\hat{\\beta})+(\\hat{{\\beta}}-{\\beta})^{T} {X}^{T} {X}(\\hat{{\\beta}}-{\\beta})+\\left({\\beta}-{\\mu}\\right)^{T} {\\Lambda}\\left({\\beta}-{\\mu}\\right)\\\\\n",
    "&={y}^{T}  {y}-{y}^{T}  {X} \\hat{ {\\beta}} + \\color{Maroon}\\underbrace{\\color{Maroon}(\\beta-\\mu_N)^TM(\\beta-\\mu_N)}_{\\text{Posterior normal kernal}}\\color{black}-\\mu_N^T M \\mu_{N}^{\\ }+\\color{black} {\\mu}^{T}  {\\Lambda}  {\\mu}+\\underbrace{\\hat{ {\\beta}}^T  {X}^{T}  {X} \\hat{ {\\beta}}}_{\\hat{ {\\beta}}^T  {X}^{T} y}\\\\\n",
    "&= \\color{Maroon}\\underbrace{\\color{Maroon}(\\beta-\\mu_{N}^{\\ })^TM(\\beta-\\mu_{N}^{\\ })}_{\\text{Posterior normal kernal}}\\color{black}+{y}^{T}  {y} -\\mu_N^T M \\mu_{N}^{\\ }+\\color{black} {\\mu}^{T}  {\\Lambda}  {\\mu}\n",
    "\\end{align}\n",
    "$$\n",
    "    \n",
    "where \n",
    "$$\n",
    "M=(X^TX+\\Lambda)\\\\\n",
    "\\mu_{N}^{\\ }=M^{-1}\\underbrace{(X^Ty+\\Lambda^T\\mu)}_{m}\n",
    "$$\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We obtain the normal part kernel in posterior, so the rest of terms in results above can combine with the inverse gamma prior $P(\\sigma^2)$ to achieve the posterior gamma kernel."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align}\n",
    "P(\\beta,\\sigma^2 |X, y) \\propto&\\left(2 \\pi \\sigma^{2}\\right)^{-N / 2}\\left|\\ {\\Lambda} \\right|^{1 / 2} \\exp \\left(-\\frac{1}{2 \\sigma^{2}}\\left[\\left(\\ {\\beta}-{\\mu}_{N}^{\\ }\\right)^{T} M\\left({\\beta}-{\\mu}_{N}\\right)\\right]\\right) \\\\\n",
    "&\\left(2 \\pi \\sigma^{2}\\right)^{-k / 2} \\exp \\left(-\\frac{1}{2 \\sigma^{2}}\\left[ {y}^{T}  {y}-{\\mu}_{N}^{T} M{\\mu}_{N}^{\\ }+ {\\mu}^T  {\\Lambda} {\\mu} \\right]\\right) \\\\\n",
    "&\\frac{\\beta ^{\\alpha}}{\\Gamma\\left(a \\right)}\\left(\\sigma^{2}\\right)^{-(\\alpha+1)} \\exp \\left(-\\beta  / \\sigma^{2}\\right)\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because it is proportional, we can safely ignore the normalizer in last two lines\n",
    "$$\n",
    "\\frac{\\beta ^{\\alpha}}{\\Gamma\\left(a \\right)}\\ \\text{ and }\\ 2\\pi^{-k/2}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align}\n",
    "P(\\beta,\\sigma^2 |X, y)\\propto&\\left(2 \\pi \\sigma^{2}\\right)^{-N / 2}\\left|\\ {\\Lambda} \\right|^{1 / 2} \\exp \\left(-\\frac{1}{2 \\sigma^{2}}\\left[\\left(\\ {\\beta}-{\\mu}_{N}^{\\ }\\right)^{T} M\\left({\\beta}-{\\mu}_{N}\\right)\\right]\\right) \\\\\n",
    "&(\\sigma^{2})^{-k / 2} \\exp \\left(-\\frac{1}{2 \\sigma^{2}}\\left[ {y}^{T}  {y}-{\\mu}_{N}^{T} M{\\mu}_{N}^{\\ }+ {\\mu}^T  {\\Lambda} {\\mu} \\right]\\right) \\\\\n",
    "&\\left(\\sigma^{2}\\right)^{-(\\alpha+1)} \\exp \\left(-\\beta  / \\sigma^{2}\\right)\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "which gives us chance to combine terms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align}\n",
    "P(\\beta,\\sigma^2 |X, y) \\propto&\\left(2 \\pi \\sigma^{2}\\right)^{-N / 2}\\left|\\ {\\Lambda} \\right|^{1 / 2} \\exp \\left(-\\frac{1}{2 \\sigma^{2}}\\left[\\left(\\ {\\beta}-{\\mu}_{N}^{\\ }\\right)^{T} M\\left({\\beta}-{\\mu}_{N}\\right)\\right]\\right) \\\\\n",
    "&(\\sigma^{2})^{-(\\alpha+\\frac{k}{2}+1)} \\exp \\left\\{-\\frac{1}{\\sigma^{2}}\\left[\\beta+\\frac{1}{2} \\left({y}^{T}  {y}-{\\mu}_{N}^{T} M{\\mu}_{N}^{\\ }+ {\\mu}^T  {\\Lambda} {\\mu} \\right)\\right]\\right\\}\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define \n",
    "$$\n",
    "\\alpha_G^{\\ } = \\alpha+\\frac{k}{2}\\\\\n",
    "\\beta_G^{\\ } = \\beta+ \\frac{1}{2}(y^Ty + \\mu^T\\Lambda\\mu -\\mu^T_{N}M \\mu_{N}^{\\ })\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div style=\"background-color:LightCoral; color:DarkBlue; padding:30px;\">\n",
    "Final conjugate posterior result is\n",
    "$$\n",
    "\\begin{align}\n",
    "P(\\beta,\\sigma^2 |X, y) \\propto&\\left(2 \\pi \\sigma^{2}\\right)^{-N / 2}\\left|\\ {\\Lambda} \\right|^{1 / 2} \\exp \\left(-\\frac{1}{2 \\sigma^{2}}\\left[\\left(\\ {\\beta}-{\\mu}_{N}^{\\ }\\right)^{T} M\\left({\\beta}-{\\mu}_{N}\\right)\\right]\\right) \\\\\n",
    "&(\\sigma^{2})^{-(\\alpha_G^{\\ }+1)} \\exp \\left\\{-\\frac{\\beta_G^{\\ }}{\\sigma^{2}}\\right\\}\n",
    "\\end{align}\n",
    "$$\n",
    "</div>\n",
    "We have successfully shown that natural conjugate prior derivation for the posterior which is a Normal-InvGamma distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{gathered}\n",
    "P\\left({\\beta}, \\sigma^{2} \\mid {X}, {y}\\right) \\propto P\\left({\\beta} \\mid {X}, {y}, \\sigma^{2}\\right) P\\left(\\sigma^{2} \\mid {X}, \\mathbf{y}\\right) \\\\\n",
    "\\text { where } \\\\\n",
    "{\\beta} \\mid {X}, {y}, \\sigma^{2} \\sim {N}_{N}\\left({\\mu}_{N}, \\sigma^{2} M^{-1}\\right) \\\\\n",
    "\\sigma^{2} \\mid {y}, {X} \\sim \\operatorname{InvGamma}\\left(\\alpha_{G}^{\\ }, \\beta_{G}^{\\ }\\right)\n",
    "\\end{gathered}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}