{
    "nbformat": 4,
    "nbformat_minor": 5,
    "metadata": {
        "kernelspec": {
            "display_name": "Python",
            "language": "python",
            "name": "python"
        },
        "language_info": {
            "name": "python"
        }
    },
    "cells": [
        {
            "id": "d6d17426",
            "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": "b23c7dbe",
            "cell_type": "markdown",
            "source": "Suppose we wish to test whether a small sample has been drawn from a normal\ndistribution. We decide that we will use the skew of the sample as a\ntest statistic, and we will consider a p-value of 0.05 to be statistically\nsignificant.\n",
            "metadata": {}
        },
        {
            "id": "6a4590a9",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "import numpy as np\nfrom scipy import stats\ndef statistic(x, axis):\n    return stats.skew(x, axis)",
            "outputs": []
        },
        {
            "id": "da6276fd",
            "cell_type": "markdown",
            "source": "After collecting our data, we calculate the observed value of the test\nstatistic.\n",
            "metadata": {}
        },
        {
            "id": "d1dfa263",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "rng = np.random.default_rng()\nx = stats.skewnorm.rvs(a=1, size=50, random_state=rng)\nstatistic(x, axis=0)",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "0.12457412450240658"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "3fe76db2",
            "cell_type": "markdown",
            "source": "To determine the probability of observing such an extreme value of the\nskewness by chance if the sample were drawn from the normal distribution,\nwe can perform a Monte Carlo hypothesis test. The test will draw many\nsamples at random from their normal distribution, calculate the skewness\nof each sample, and compare our original skewness against this\ndistribution to determine an approximate p-value.\n",
            "metadata": {}
        },
        {
            "id": "2b58dcd8",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "from scipy.stats import monte_carlo_test\n# because our statistic is vectorized, we pass `vectorized=True`\nrvs = lambda size: stats.norm.rvs(size=size, random_state=rng)\nres = monte_carlo_test(x, rvs, statistic, vectorized=True)\nprint(res.statistic)",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "0.12457412450240658"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "3d71b4f0",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "print(res.pvalue)",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "0.7012"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "83078f01",
            "cell_type": "markdown",
            "source": "The probability of obtaining a test statistic less than or equal to the\nobserved value under the null hypothesis is ~70%. This is greater than\nour chosen threshold of 5%, so we cannot consider this to be significant\nevidence against the null hypothesis.\n\nNote that this p-value essentially matches that of\n`scipy.stats.skewtest`, which relies on an asymptotic distribution of a\ntest statistic based on the sample skewness.\n",
            "metadata": {}
        },
        {
            "id": "a5abb7dd",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "stats.skewtest(x).pvalue",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "0.6892046027110614"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "0c9bca8a",
            "cell_type": "markdown",
            "source": "This asymptotic approximation is not valid for small sample sizes, but\n`monte_carlo_test` can be used with samples of any size.\n",
            "metadata": {}
        },
        {
            "id": "2035619f",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "x = stats.skewnorm.rvs(a=1, size=7, random_state=rng)\n# stats.skewtest(x) would produce an error due to small sample\nres = monte_carlo_test(x, rvs, statistic, vectorized=True)",
            "outputs": []
        },
        {
            "id": "83107613",
            "cell_type": "markdown",
            "source": "The Monte Carlo distribution of the test statistic is provided for\nfurther investigation.\n",
            "metadata": {}
        },
        {
            "id": "e8cfd080",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "import matplotlib.pyplot as plt\nfig, ax = plt.subplots()\nax.hist(res.null_distribution, bins=50)\nax.set_title(\"Monte Carlo distribution of test statistic\")\nax.set_xlabel(\"Value of Statistic\")\nax.set_ylabel(\"Frequency\")\nplt.show()",
            "outputs": []
        }
    ]
}