{
    "nbformat": 4,
    "nbformat_minor": 5,
    "metadata": {
        "kernelspec": {
            "display_name": "Python",
            "language": "python",
            "name": "python"
        },
        "language_info": {
            "name": "python"
        }
    },
    "cells": [
        {
            "id": "7c1fea13",
            "cell_type": "markdown",
            "source": "<div class='alert alert-warning'>\n\nSciPy's interactive examples with Jupyterlite are experimental and may not always work as expected. Execution of cells containing imports may result in large downloads (up to 60MB of content for the first import from SciPy). Load times when importing from SciPy may take roughly 10-20 seconds. If you notice any problems, feel free to open an [issue](https://github.com/scipy/scipy/issues/new/choose).\n\n</div>",
            "metadata": {}
        },
        {
            "id": "61e9040d",
            "cell_type": "markdown",
            "source": "Suppose we wish to test whether two samples are drawn from the same\ndistribution. Assume that the underlying distributions are unknown to us,\nand that before observing the data, we hypothesized that the mean of the\nfirst sample would be less than that of the second sample. We decide that\nwe will use the difference between the sample means as a test statistic,\nand we will consider a p-value of 0.05 to be statistically significant.\n\nFor efficiency, we write the function defining the test statistic in a\nvectorized fashion: the samples ``x`` and ``y`` can be ND arrays, and the\nstatistic will be calculated for each axis-slice along `axis`.\n",
            "metadata": {}
        },
        {
            "id": "50198588",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "import numpy as np\ndef statistic(x, y, axis):\n    return np.mean(x, axis=axis) - np.mean(y, axis=axis)",
            "outputs": []
        },
        {
            "id": "91740ec4",
            "cell_type": "markdown",
            "source": "After collecting our data, we calculate the observed value of the test\nstatistic.\n",
            "metadata": {}
        },
        {
            "id": "4c42ec0c",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "from scipy.stats import norm\nrng = np.random.default_rng()\nx = norm.rvs(size=5, random_state=rng)\ny = norm.rvs(size=6, loc = 3, random_state=rng)\nstatistic(x, y, 0)",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "-3.5411688580987266"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "2c51180b",
            "cell_type": "markdown",
            "source": "Indeed, the test statistic is negative, suggesting that the true mean of\nthe distribution underlying ``x`` is less than that of the distribution\nunderlying ``y``. To determine the probability of this occurring by chance\nif the two samples were drawn from the same distribution, we perform\na permutation test.\n",
            "metadata": {}
        },
        {
            "id": "f877b087",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "from scipy.stats import permutation_test\n# because our statistic is vectorized, we pass `vectorized=True`\n# `n_resamples=np.inf` indicates that an exact test is to be performed\nres = permutation_test((x, y), statistic, vectorized=True,\n                       n_resamples=np.inf, alternative='less')\nprint(res.statistic)",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "-3.5411688580987266"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "55119095",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "print(res.pvalue)",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "0.004329004329004329"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "2ddc6a0d",
            "cell_type": "markdown",
            "source": "The probability of obtaining a test statistic less than or equal to the\nobserved value under the null hypothesis is 0.4329%. This is less than our\nchosen threshold of 5%, so we consider this to be significant evidence\nagainst the null hypothesis in favor of the alternative.\n\nBecause the size of the samples above was small, `permutation_test` could\nperform an exact test. For larger samples, we resort to a randomized\npermutation test.\n",
            "metadata": {}
        },
        {
            "id": "b741d364",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "x = norm.rvs(size=100, random_state=rng)\ny = norm.rvs(size=120, loc=0.2, random_state=rng)\nres = permutation_test((x, y), statistic, n_resamples=9999,\n                       vectorized=True, alternative='less',\n                       rng=rng)\nprint(res.statistic)",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "-0.4230459671240913"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "b558999e",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "print(res.pvalue)",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "0.0015"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "94572e7b",
            "cell_type": "markdown",
            "source": "The approximate probability of obtaining a test statistic less than or\nequal to the observed value under the null hypothesis is 0.0225%. This is\nagain less than our chosen threshold of 5%, so again we have significant\nevidence to reject the null hypothesis in favor of the alternative.\n\nFor large samples and number of permutations, the result is comparable to\nthat of the corresponding asymptotic test, the independent sample t-test.\n",
            "metadata": {}
        },
        {
            "id": "abb86489",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "from scipy.stats import ttest_ind\nres_asymptotic = ttest_ind(x, y, alternative='less')\nprint(res_asymptotic.pvalue)",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "0.0014669545224902675"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "ca8e1a37",
            "cell_type": "markdown",
            "source": "The permutation distribution of the test statistic is provided for\nfurther investigation.\n",
            "metadata": {}
        },
        {
            "id": "9ba865a2",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "import matplotlib.pyplot as plt\nplt.hist(res.null_distribution, bins=50)\nplt.title(\"Permutation distribution of test statistic\")\nplt.xlabel(\"Value of Statistic\")\nplt.ylabel(\"Frequency\")\nplt.show()",
            "outputs": []
        },
        {
            "id": "e57c1f0d",
            "cell_type": "markdown",
            "source": "Inspection of the null distribution is essential if the statistic suffers\nfrom inaccuracy due to limited machine precision. Consider the following\ncase:\n",
            "metadata": {}
        },
        {
            "id": "c6c8b0b8",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "from scipy.stats import pearsonr\nx = [1, 2, 4, 3]\ny = [2, 4, 6, 8]\ndef statistic(x, y, axis=-1):\n    return pearsonr(x, y, axis=axis).statistic\nres = permutation_test((x, y), statistic, vectorized=True,\n                       permutation_type='pairings',\n                       alternative='greater')\nr, pvalue, null = res.statistic, res.pvalue, res.null_distribution",
            "outputs": []
        },
        {
            "id": "22ca3f53",
            "cell_type": "markdown",
            "source": "In this case, some elements of the null distribution differ from the\nobserved value of the correlation coefficient ``r`` due to numerical noise.\nWe manually inspect the elements of the null distribution that are nearly\nthe same as the observed value of the test statistic.\n",
            "metadata": {}
        },
        {
            "id": "beb05c5c",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "r",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "0.7999999999999999"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "9c5bbf38",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "unique = np.unique(null)\nunique",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([-1. , -1. , -0.8, -0.8, -0.8, -0.6, -0.4, -0.4, -0.2, -0.2, -0.2,\n    0. ,  0.2,  0.2,  0.2,  0.4,  0.4,  0.6,  0.8,  0.8,  0.8,  1. ,\n    1. ])  # may vary"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "e0552008",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "unique[np.isclose(r, unique)].tolist()",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "[0.7999999999999998, 0.7999999999999999, 0.8]  # may vary"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "650a0276",
            "cell_type": "markdown",
            "source": "If `permutation_test` were to perform the comparison naively, the\nelements of the null distribution with value ``0.7999999999999998`` would\nnot be considered as extreme or more extreme as the observed value of the\nstatistic, so the calculated p-value would be too small.\n",
            "metadata": {}
        },
        {
            "id": "06daad83",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "incorrect_pvalue = np.count_nonzero(null >= r) / len(null)\nincorrect_pvalue",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "0.14583333333333334  # may vary"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "16a67153",
            "cell_type": "markdown",
            "source": "Instead, `permutation_test` treats elements of the null distribution that\nare within ``max(1e-14, abs(r)*1e-14)`` of the observed value of the\nstatistic ``r`` to be equal to ``r``.\n",
            "metadata": {}
        },
        {
            "id": "77982495",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "correct_pvalue = np.count_nonzero(null >= r - 1e-14) / len(null)\ncorrect_pvalue",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "0.16666666666666666"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "65eb7089",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "res.pvalue == correct_pvalue",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "True"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "3769b79d",
            "cell_type": "markdown",
            "source": "This method of comparison is expected to be accurate in most practical\nsituations, but the user is advised to assess this by inspecting the\nelements of the null distribution that are close to the observed value\nof the statistic. Also, consider the use of statistics that can be\ncalculated using exact arithmetic (e.g. integer statistics).\n",
            "metadata": {}
        }
    ]
}