example/18-hypothesis-eeibf-validation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Validate implementation of *Expected Encompassing Intrinsic Bayes Factor* (EEIBF) normal mean hypothesis testing with unknown variance and an equality hypothesis.\n",
    "\n",
    "The reference used to validate is the paper *Default Bayes Factors for Non-Nested Hypothesis Testing* by J. Berger and J. Mortera.\n",
    "\n",
    "1: Berger, J. and J. Mortera (1999). Default bayes factors for nonnested hypothesis testing. Journal of the American Statistical Association 94 (446), 542–554.\n",
    "Postscript for paper: https://www2.stat.duke.edu/~berger/papers/mortera.ps"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Notation\n",
    "From [1] section 4, the three hypotheses are\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "H_1&: \\mu=0 \\\\\n",
    "H_2&: \\mu<0 \\\\\n",
    "H_3&: \\mu>0\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "and $H_0$ denotes the encompassing hypothesis (i.e. $\\mu$ is unrestricted)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy\n",
    "from IPython.display import display, Math\n",
    "from bbai.stat import NormalMeanHypothesis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sample Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mu_ml = -0.6493372888888889\n",
      "sigma_ml = 0.5302624236732443\n"
     ]
    }
   ],
   "source": [
    "# define some sample data to validate the equations against\n",
    "data = [\n",
    "        -0.619487, -0.685728, -0.798718, \n",
    "        -0.353721, 0.36195, -0.0545616, \n",
    "        -1.09919, -1.28801, -1.30657\n",
    "]\n",
    "n = len(data)\n",
    "mu_ml = np.mean(data)\n",
    "sigma_ml = np.std(data)\n",
    "print('mu_ml =', mu_ml)\n",
    "print('sigma_ml =', sigma_ml)\n",
    "hypothesis = NormalMeanHypothesis().test(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Uncorrected Factors\n",
    "\n",
    "We'll first compute the uncorrected Bayes factors, $B_{i0}^N$. See equation (2) of [1].\n",
    "\n",
    "The reference prior for the normal distribution with the mean, $\\mu$, as the parameter of interest and unknown variance $\\sigma^2$ is given by\n",
    "$$\n",
    "\\pi(\\mu, \\sigma^2) = \\frac{1}{\\sigma^2}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### $H_1$ uncorrected factor"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "B_{10}^N = \\frac{\\int_0^\\infty P(y\\mid \\mu=0, \\sigma^2) \\frac{1}{\\sigma^2} d\\sigma^2}{\n",
    "\\int_{-\\infty}^\\infty \\int_0^\\infty P(y\\mid \\mu, \\sigma^2) \\frac{1}{\\sigma^2} d\\sigma^2 d\\mu\n",
    "}\n",
    "$$\n",
    "Now,\n",
    "$$\n",
    "\\int_0^\\infty P(y \\mid \\mu, \\sigma^2) \\times \\frac{1}{\\sigma^2} \\propto\n",
    "   \\mathcal{t}_{n-1}\\left[ \\sqrt{n-1} \\left(\\frac{\\mu - \\mu_{ml}}{\\sigma_{ml}}\\right)\\right]\n",
    "$$\n",
    "where $\\mathcal{t}_{n-1}$ is the PDF of the student t distribution with $n-1$ degrees of freedom.\n",
    "\n",
    "Thus,\n",
    "$$\n",
    "B_{10}^N = \\frac{\\sqrt{n-1}}{\\sigma_{ml}} \\times \\mathcal{t}_{n-1}\\left[ \\sqrt{n-1} \\left(\\frac{- \\mu_{ml}}{\\sigma_{ml}}\\right)\\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_B_10_N(y):\n",
    "    n = len(y)\n",
    "    t_dist = scipy.stats.t(n-1)\n",
    "    mu_ml = np.mean(y)\n",
    "    sigma_ml = np.std(y)\n",
    "    return np.sqrt(n-1) / sigma_ml * t_dist.pdf(-np.sqrt(n - 1) * mu_ml / sigma_ml)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle B_{10}^N = 0.03342368282076605\\;\\textrm{vs}\\;0.03342368282076605$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "B_10_N = compute_B_10_N(data)\n",
    "display(Math(r'B_{10}^N = ' + str(B_10_N) + r'\\;\\textrm{vs}\\;' + str(hypothesis.equal_to_star_N)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### $H_2$ uncorrected factor"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "B_{20}^N &= \\frac{\\int_{-\\infty}^0\\int_0^\\infty P(y\\mid \\mu, \\sigma^2) \\frac{1}{\\sigma^2} d\\sigma^2 d\\mu}{\n",
    "\\int_{-\\infty}^\\infty \\int_0^\\infty P(y\\mid \\mu, \\sigma^2) \\frac{1}{\\sigma^2} d\\sigma^2 d\\mu\n",
    "} \\\\\n",
    "         &= \\frac{\n",
    "                \\int_{-\\infty}^0 \\mathcal{t}_{n-1}\\left[\\sqrt{n-1}\\left(\\frac{\\mu - \\mu_{ml}}{\\sigma_{ml}}\\right)\\right]d\\mu\n",
    "             }{\n",
    "                 \\int_{-\\infty}^\\infty \\mathcal{t}_{n-1}\\left[\\sqrt{n-1}\\left(\\frac{\\mu - \\mu_{ml}}{\\sigma_{ml}}\\right)\\right]d\\mu\n",
    "             } \\\\\n",
    "          &= \\mathcal{T}_{n-1}\\left[\\sqrt{n-1}\\left(\\frac{- \\mu_{ml}}{\\sigma_{ml}}\\right)\\right]\n",
    "\\end{align*}\n",
    "$$\n",
    "where $\\mathcal{T}_{n-1}$ represents the CDF of the Student t distribution with $n-1$ degrees of freedom"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_B_20_N(y):\n",
    "    n = len(y)\n",
    "    t_dist = scipy.stats.t(n-1)\n",
    "    mu_ml = np.mean(y)\n",
    "    sigma_ml = np.std(y)\n",
    "    return t_dist.cdf(-np.sqrt(n - 1) * mu_ml / sigma_ml)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle B_{20}^N = 0.9957385642382396\\;\\textrm{vs}\\;0.9957385642382396$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "B_20_N = compute_B_20_N(data)\n",
    "display(Math(r'B_{20}^N = ' + str(B_20_N) + r'\\;\\textrm{vs}\\;' + str(hypothesis.left_to_star_N)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### $H_3$ uncorrected factor"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similarly, for $H_3$, we have\n",
    "$$\n",
    "B_{30}^N = 1 - \\mathcal{T}_{n-1}\\left[\\sqrt{n-1}\\left(\\frac{- \\mu_{ml}}{\\sigma_{ml}}\\right)\\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_B_30_N(y):\n",
    "    return 1 - compute_B_20_N(y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle B_{30}^N = 0.004261435761760413\\;\\textrm{vs}\\;0.004261435761760413$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "B_30_N = compute_B_30_N(data)\n",
    "display(Math(r'B_{30}^N = ' + str(B_30_N) + r'\\;\\textrm{vs}\\;' + str(hypothesis.right_to_star_N)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Correction Factors\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From section 2.4.2 of [1], for testing the mean of normally distributed data $\\boldsymbol{x}$ with unknown variance, the full EEIBF Bayes factor is given by\n",
    "\n",
    "$$\n",
    "B_{ji}^{EEI} = B_{ji}^N(\\boldsymbol{x}) \\times \n",
    "\\frac{\n",
    " \\mathbb{E}_\\hat{\\boldsymbol{\\theta}}^{H_0} \\left[B_{i0}^N(X_1, X_2)\\right]\n",
    " }{\n",
    " \\mathbb{E}_\\hat{\\boldsymbol{\\theta}}^{H_0} \\left[B_{j0}^N(X_1, X_2)\\right]\n",
    " }\n",
    "$$\n",
    "where $X_1$ and $X_2$ denote independent random variables distributed under the\n",
    "model for $H_0$ with parameters the $\\hat{\\boldsymbol{\\theta}}=(\\mu_{ml}, \\sigma_{ml})$ that maximize the likelhood of the observed data $\\boldsymbol{x}$.\n",
    "\n",
    "bbai uses an efficient deterministic algorithm to compute the values\n",
    "$$\n",
    "\\mathbb{E}_\\hat{\\boldsymbol{\\theta}}^{H_0} \\left[B_{i0}^N(X_1, X_2)\\right]\n",
    "$$\n",
    "But to validate, we'll use the slower but more straightforward method of generating a bunch of random values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate some random values\n",
    "np.random.seed(0)\n",
    "N = 10000 \n",
    "    # if we select N larger (e.g. 1000000), we'll get a closer match; \n",
    "    # but this should be reasonable for a ball park test\n",
    "rx = np.random.normal(mu_ml, sigma_ml, size=(N, 2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### $H_1$ correction factor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_E_B_10_Xl(rx):\n",
    "    N = len(rx)\n",
    "    res = 0\n",
    "    for xli in rx:\n",
    "        res += compute_B_10_N(xli)\n",
    "    return res / N"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\mathbb{E}_{\\boldsymbol{\\theta}}^{H_0}\\left[B_{10}^N(X_1, X_2)\\right] = 0.2687721494618709\\;\\textrm{vs}\\;0.27755975516051407$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "E_B_10_Xl = compute_E_B_10_Xl(rx)\n",
    "display(Math(r'\\mathbb{E}_{\\boldsymbol{\\theta}}^{H_0}\\left[B_{10}^N(X_1, X_2)\\right] = ' + str(E_B_10_Xl) + r'\\;\\textrm{vs}\\;' + str(hypothesis.correction_equal)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### $H_2$ correction factor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_E_B_20_Xl(rx):\n",
    "    N = len(rx)\n",
    "    res = 0\n",
    "    for xli in rx:\n",
    "        res += compute_B_20_N(xli)\n",
    "    return res / N"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\mathbb{E}_{\\boldsymbol{\\theta}}^{H_0}\\left[B_{20}^N(X_1, X_2)\\right] = 0.8351464188652707\\;\\textrm{vs}\\;0.8332747776394973$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "E_B_20_Xl = compute_E_B_20_Xl(rx)\n",
    "display(Math(r'\\mathbb{E}_{\\boldsymbol{\\theta}}^{H_0}\\left[B_{20}^N(X_1, X_2)\\right] = ' + str(E_B_20_Xl) + r'\\;\\textrm{vs}\\;' + str(hypothesis.correction_left)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### $H_3$ correction factor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_E_B_30_Xl(rx):\n",
    "    N = len(rx)\n",
    "    res = 0\n",
    "    for xli in rx:\n",
    "        res += compute_B_30_N(xli)\n",
    "    return res / N"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\mathbb{E}_{\\boldsymbol{\\theta}}^{H_0}\\left[B_{30}^N(X_1, X_2)\\right] = 0.16485358113473542\\;\\textrm{vs}\\;0.1667252223605027$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "E_B_30_Xl = compute_E_B_30_Xl(rx)\n",
    "display(Math(r'\\mathbb{E}_{\\boldsymbol{\\theta}}^{H_0}\\left[B_{30}^N(X_1, X_2)\\right] = ' + str(E_B_30_Xl) + r'\\;\\textrm{vs}\\;' + str(hypothesis.correction_right)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Full Bayes Factors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "B_{i0}^{EEI} = B_{i0}^N(\\boldsymbol{x}) \\times \n",
    "\\frac{\n",
    " 1\n",
    " }{\n",
    " \\mathbb{E}_\\hat{\\boldsymbol{\\theta}}^{H_0} \\left[B_{i0}^N(X_1, X_2)\\right]\n",
    " }\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle B_{10} = 0.12435694281452205\\;\\textrm{vs}\\;0.1204197733977571$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "B_10 = B_10_N / E_B_10_Xl\n",
    "display(Math(r'B_{10} = ' + str(B_10) + r'\\;\\textrm{vs}\\;' + str(hypothesis.factor_equal)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle B_{20} = 1.1922922037924422\\;\\textrm{vs}\\;1.1949702438599787$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "B_20 = B_20_N / E_B_20_Xl\n",
    "display(Math(r'B_{20} = ' + str(B_20) + r'\\;\\textrm{vs}\\;' + str(hypothesis.factor_left)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle B_{30} = 0.025849822202391383\\;\\textrm{vs}\\;0.025559634597734082$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "B_30 = B_30_N / E_B_30_Xl\n",
    "display(Math(r'B_{30} = ' + str(B_30) + r'\\;\\textrm{vs}\\;' + str(hypothesis.factor_right)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Probabilities from factors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.09263094 0.88811406 0.019255  ] vs [0.08980186 0.89113729 0.01906085]\n"
     ]
    }
   ],
   "source": [
    "T = B_10 + B_20 + B_30\n",
    "expected = np.array([B_10, B_20, B_30] / T)\n",
    "actual = np.array([hypothesis.equal, hypothesis.left, hypothesis.right])\n",
    "print(expected, 'vs', actual)"
   ]
  },
  {
   "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.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}