{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "95f9180b",
   "metadata": {},
   "source": [
    "# Jarque-Bera goodness of fit test\n",
    "\n",
    "Suppose we wish to infer from measurements whether the weights of adult human\n",
    "males in a medical study are not normally distributed [^1]. The weights (lbs)\n",
    "are recorded in the array `x` below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "df770b3b",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "x = np.array([148, 154, 158, 160, 161, 162, 166, 170, 182, 195, 236])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3a6447d8",
   "metadata": {},
   "source": [
    "The Jarque-Bera test {class}`scipy.stats.jarque_bera` begins by computing a\n",
    "statistic based on the sample skewness and kurtosis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "176cf240",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import stats\n",
    "res = stats.jarque_bera(x)\n",
    "res.statistic"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e20b2a23",
   "metadata": {},
   "source": [
    "Because the normal distribution has zero skewness and zero (\"excess\" or\n",
    "\"Fisher\") kurtosis, the value of this statistic tends to be low for samples\n",
    "drawn from a normal distribution.\n",
    "\n",
    "The test is performed by comparing the observed value of the statistic against\n",
    "the null distribution: the distribution of statistic values derived under the\n",
    "null hypothesis that the weights were drawn from a normal distribution.\n",
    "\n",
    "For the Jarque-Bera test, the null distribution for very large samples is the\n",
    "{class}`chi-squared distribution <scipy.stats.chi2>` with two degrees of freedom."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8a54f9fa",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "dist = stats.chi2(df=2)\n",
    "jb_val = np.linspace(0, 11, 100)\n",
    "pdf = dist.pdf(jb_val)\n",
    "fig, ax = plt.subplots(figsize=(8, 5))\n",
    "\n",
    "def jb_plot(ax):  # we'll reuse this\n",
    "    ax.plot(jb_val, pdf)\n",
    "    ax.set_title(\"Jarque-Bera Null Distribution\")\n",
    "    ax.set_xlabel(\"statistic\")\n",
    "    ax.set_ylabel(\"probability density\")\n",
    "\n",
    "jb_plot(ax)\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "83b44e1d",
   "metadata": {},
   "source": [
    "The comparison is quantified by the p-value: the proportion of values in the\n",
    "null distribution greater than or equal to the observed value of the statistic."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "18029d6b",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(8, 5))\n",
    "jb_plot(ax)\n",
    "pvalue = dist.sf(res.statistic)\n",
    "annotation = (f'p-value={pvalue:.6f}\\n(shaded area)')\n",
    "props = dict(facecolor='black', width=1, headwidth=5, headlength=8)\n",
    "_ = ax.annotate(annotation, (7.5, 0.01), (8, 0.05), arrowprops=props)\n",
    "i = jb_val >= res.statistic  # indices of more extreme statistic values\n",
    "ax.fill_between(jb_val[i], y1=0, y2=pdf[i])\n",
    "ax.set_xlim(0, 11)\n",
    "ax.set_ylim(0, 0.3)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "89bc45f5",
   "metadata": {},
   "outputs": [],
   "source": [
    "res.pvalue"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "72e8629c",
   "metadata": {},
   "source": [
    "If the p-value is \"small\" - that is, if there is a low probability of sampling\n",
    "data from a normally distributed population that produces such an extreme value\n",
    "of the statistic - this may be taken as evidence against the null hypothesis in\n",
    "favor of the alternative: the weights were not drawn from a normal distribution.\n",
    "Note that:\n",
    "\n",
    "- The inverse is not true; that is, the test is not used to provide evidence for\n",
    "  the null hypothesis.\n",
    "- The threshold for values that will be considered \"small\" is a choice that\n",
    "  should be made before the data is analyzed [^2] with consideration of the\n",
    "  risks of both false positives (incorrectly rejecting the null hypothesis) and\n",
    "  false negatives (failure to reject a false null hypothesis).\n",
    "\n",
    "Note that the chi-squared distribution provides an asymptotic approximation\n",
    "of the null distribution; it is only accurate for samples with many\n",
    "observations. For small samples like ours, {class}`scipy.stats.monte_carlo_test`\n",
    "may provide a more accurate, albeit stochastic, approximation of the\n",
    "exact p-value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0ca7a7a7",
   "metadata": {},
   "outputs": [],
   "source": [
    "def statistic(x, axis):\n",
    "    # underlying calculation of the Jarque Bera statistic\n",
    "    s = stats.skew(x, axis=axis)\n",
    "    k = stats.kurtosis(x, axis=axis)\n",
    "    return x.shape[axis]/6 * (s**2 + k**2/4)\n",
    "\n",
    "res = stats.monte_carlo_test(x, stats.norm.rvs, statistic,\n",
    "                             alternative='greater')\n",
    "fig, ax = plt.subplots(figsize=(8, 5))\n",
    "jb_plot(ax)\n",
    "ax.hist(res.null_distribution, np.linspace(0, 10, 50),\n",
    "        density=True)\n",
    "ax.legend(['asymptotic approximation (many observations)',\n",
    "           'Monte Carlo approximation (11 observations)'])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3281b678",
   "metadata": {},
   "outputs": [],
   "source": [
    "res.pvalue"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f7f62846",
   "metadata": {},
   "source": [
    "Furthermore, despite their stochastic nature, p-values computed in this way\n",
    "can be used to exactly control the rate of false rejections of the null\n",
    "hypothesis [^3].\n",
    "\n",
    "## References\n",
    "\n",
    "[^1]: Shapiro, S. S., & Wilk, M. B. (1965). An analysis of variance test for\n",
    "normality (complete samples). Biometrika, 52(3/4), 591-611\n",
    "[^2]: Phipson, B. and Smyth, G. K. (2010) Permutation P-values Should Never Be\n",
    "Zero: Calculating Exact P-values When Permutations Are Randomly Drawn.\n",
    "Statistical Applications in Genetics and Molecular Biology 9.1.\n",
    "[^3]: Panagiotakos, D. B. (2008). The value of p-value in biomedical research.\n",
    "The open cardiovascular medicine journal, 2, 97."
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "text_representation": {
    "extension": ".md",
    "format_name": "myst",
    "format_version": 0.13,
    "jupytext_version": "1.16.1"
   }
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
