Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- {
- "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
- }
Add Comment
Please, Sign In to add comment