{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "172d07be",
   "metadata": {},
   "source": [
    "# Kurtosis test\n",
    "\n",
    "The kurtosis test {func}`scipy.stats.kurtosistest` function tests the null\n",
    "hypothesis that the kurtosis of the population from which the sample was drawn\n",
    "is that of the normal distribution.\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": "54673265",
   "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": "1b300934",
   "metadata": {},
   "source": [
    "The kurtosis test from [^2] begins by computing a statistic based on the sample\n",
    "(excess/Fisher) kurtosis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cafe1f2a",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import stats\n",
    "res = stats.kurtosistest(x)\n",
    "res.statistic"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48cc2f07",
   "metadata": {},
   "source": [
    "(The test warns that our sample has too few observations to perform the test.\n",
    "We'll return to this at the end of the example.) Because normal distributions\n",
    "have zero excess kurtosis (by definition), the magnitude of this statistic tends\n",
    "to be low for samples 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 this test, the null distribution of the statistic for very large samples is\n",
    "the standard normal distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f928dc08",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "dist = stats.norm()\n",
    "kt_val = np.linspace(-5, 5, 100)\n",
    "pdf = dist.pdf(kt_val)\n",
    "fig, ax = plt.subplots(figsize=(8, 5))\n",
    "\n",
    "def kt_plot(ax):  # we'll reuse this\n",
    "    ax.plot(kt_val, pdf)\n",
    "    ax.set_title(\"Kurtosis Test Null Distribution\")\n",
    "    ax.set_xlabel(\"statistic\")\n",
    "    ax.set_ylabel(\"probability density\")\n",
    "\n",
    "kt_plot(ax)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ea013548",
   "metadata": {},
   "source": [
    "The comparison is quantified by the p-value: the proportion of values in the\n",
    "null distribution as extreme or more extreme than the observed value of the\n",
    "statistic. In a two-sided test in which the statistic is positive, elements of\n",
    "the null distribution greater than the observed statistic and elements of the\n",
    "null distribution less than the negative of the observed statistic are both\n",
    "considered \"more extreme\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cfdd0b2a",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(8, 5))\n",
    "kt_plot(ax)\n",
    "pvalue = dist.cdf(-res.statistic) + dist.sf(res.statistic)\n",
    "annotation = (f'p-value={pvalue:.3f}\\n(shaded area)')\n",
    "props = dict(facecolor='black', width=1, headwidth=5, headlength=8)\n",
    "_ = ax.annotate(annotation, (3, 0.005), (3.25, 0.02), arrowprops=props)\n",
    "i = kt_val >= res.statistic\n",
    "ax.fill_between(kt_val[i], y1=0, y2=pdf[i], color='C0')\n",
    "i = kt_val <= -res.statistic\n",
    "ax.fill_between(kt_val[i], y1=0, y2=pdf[i], color='C0')\n",
    "ax.set_xlim(-5, 5)\n",
    "ax.set_ylim(0, 0.1)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d633895d",
   "metadata": {},
   "outputs": [],
   "source": [
    "res.pvalue"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d4e624a8",
   "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\n",
    "  evidence for 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 [^3] with consideration of the\n",
    "  risks of both false positives (incorrectly rejecting the null hypothesis)\n",
    "  and false negatives (failure to reject a false null hypothesis).\n",
    "\n",
    "Note that the standard normal distribution provides an asymptotic approximation\n",
    "of the null distribution; it is only accurate for samples with many\n",
    "observations. This is the reason we received a warning at the beginning of the\n",
    "example; our sample is quite small. In this case,\n",
    "{class}`scipy.stats.monte_carlo_test` may provide a more accurate, albeit\n",
    "stochastic, approximation of the exact p-value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "843b1ff6",
   "metadata": {},
   "outputs": [],
   "source": [
    "def statistic(x, axis):\n",
    "    # get just the skewtest statistic; ignore the p-value\n",
    "    return stats.kurtosistest(x, axis=axis).statistic\n",
    "res = stats.monte_carlo_test(x, stats.norm.rvs, statistic)\n",
    "fig, ax = plt.subplots(figsize=(8, 5))\n",
    "kt_plot(ax)\n",
    "ax.hist(res.null_distribution, np.linspace(-5, 5, 50),\n",
    "        density=True)\n",
    "ax.legend(['asymptotic approximation\\n(many observations)',\n",
    "           'Monte Carlo approximation\\n(11 observations)'])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "55f233d9",
   "metadata": {},
   "outputs": [],
   "source": [
    "res.pvalue"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3ffe18bf",
   "metadata": {},
   "source": [
    "Furthermore, despite their stochastic nature, p-values computed in this way can\n",
    "be used to exactly control the rate of false rejections of the null hypothesis [^4].\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]: Anscombe, F. J. and Glynn, W. J. (1983). Distribution of the kurtosis\n",
    "statistic b2 for normal samples, Biometrika, vol. 70, pp. 227-234.\n",
    "[^3]: 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",
    "[^4]: 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
}
