derivation for mode of Normal Inverse-Chisq distribution

nix2-mode.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Mode of normal-inverse Chi-squared distribution\n",
    "\n",
    "The formula given in murphy seems to be wrong."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "using Distributions, ConjugatePriors\n",
    "using ConjugatePriors: NormalInverseChisq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "ConjugatePriors.NormalInverseChisq{Float64}(μ=0.0, σ2=10.0, κ=2.0, ν=3.0)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nix = NormalInverseChisq(0., 10., 2., 3.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.0, 15.0)"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "μ, σ2 = mode(nix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0040313337318566775"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pdf(nix, μ, σ2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(false, true)"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pdf.(nix, μ, σ2 .+ (-0.001, 0.001)) .< pdf(nix, μ, σ2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Math\n",
    "\n",
    "The mode for the mean is clearly $\\mu_0$.\n",
    "\n",
    "Let $f(\\sigma^2)$ be the pdf as a function of the variance.  Then based on Murphy (125),\n",
    "\n",
    "\\begin{align}\n",
    "f(x) \\propto& x^{-1/2} x^{-(\\nu/2 + 1)} \\exp\\left( \\frac{-1}{2} x^{-1} (\\nu\\sigma^2_0 + \\kappa_0(\\mu_0-\\mu)^2) \\right) \\\\\n",
    "     \\propto& x^{-\\frac{\\nu+3}{2}} \\exp\\left(\\frac{-1}{2} x^{-1} \\nu\\sigma^2_0 \\right)\n",
    "\\end{align}\n",
    "\n",
    "(dropping constants since we're solving for $f^\\prime = 0$).\n",
    "\n",
    "Let\n",
    "\\begin{align}\n",
    "A &= -\\frac{\\nu+3}{2} \\\\\n",
    "B &= -\\frac{\\nu\\sigma^2_0}{2}\n",
    "\\end{align}\n",
    "\n",
    "Then $f(x) \\propto x^A \\exp(Bx^{-1})$ and $f^\\prime(x) \\propto A x^{A-1} \\exp(Bx^{-1}) - x^A Bx^{-2}\\exp(Bx^{-1})$.  To find the mode,\n",
    "\n",
    "\\begin{align}\n",
    "f^\\prime(x) &= 0 \\\\\n",
    "Ax^{A-1} \\exp(Bx^{-1}) &= x^A Bx^{-2} \\exp(Bx^{-1}) \\\\\n",
    "A x^{-1} &= Bx^{-2} \\\\\n",
    "x &= \\frac{B}{A} \\\\\n",
    "  &= \\frac{\\nu\\sigma^2_0}{\\nu+3}\n",
    "\\end{align}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.0, 5.0)"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mode2(d::NormalInverseChisq) = (d.μ, d.ν * d.σ2 / (d.ν + 3))\n",
    "μ, σ2 = mode2(nix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.014730705695397627"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pdf(nix, μ, σ2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(true, true)"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pdf.(nix, μ, σ2 .+ (-0.001, 0.001)) .< pdf(nix, μ, σ2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Normal-Inverse gamma prior\n",
    "\n",
    "This is equivalent with the following (relevant) re-parametrizations:\n",
    "\n",
    "\\begin{align}\n",
    "a_0 &= \\nu_0/2 \\\\ \n",
    "b_0 &= \\nu_0\\sigma^2_0 / 2\n",
    "\\end{align}\n",
    "\n",
    "where $a_0$ is the shape parameter and $b_0$ is the scale.  With this parametrization, the mode is $\\frac{b_0}{a_0 + 3/2}$"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.6.2",
   "language": "julia",
   "name": "julia-0.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}