{
    "nbformat": 4,
    "nbformat_minor": 5,
    "metadata": {
        "kernelspec": {
            "display_name": "Python",
            "language": "python",
            "name": "python"
        },
        "language_info": {
            "name": "python"
        }
    },
    "cells": [
        {
            "id": "abbfdba4",
            "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": "b97395ee",
            "cell_type": "markdown",
            "source": "Suppose we wish to simulate the power of the independent sample t-test\nunder the following conditions:\n\n- The first sample has 10 observations drawn from a normal distribution\n  with mean 0.\n- The second sample has 12 observations drawn from a normal distribution\n  with mean 1.0.\n- The threshold on p-values for significance is 0.05.\n",
            "metadata": {}
        },
        {
            "id": "ebf03bce",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "import numpy as np\nfrom scipy import stats\nrng = np.random.default_rng()\n\ntest = stats.ttest_ind\nn_observations = (10, 12)\nrvs1 = rng.normal\nrvs2 = lambda size: rng.normal(loc=1, size=size)\nrvs = (rvs1, rvs2)\nres = stats.power(test, rvs, n_observations, significance=0.05)\nres.power",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "0.6116"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "b331d146",
            "cell_type": "markdown",
            "source": "With samples of size 10 and 12, respectively, the power of the t-test\nwith a significance threshold of 0.05 is approximately 60% under the chosen\nalternative. We can investigate the effect of sample size on the power\nby passing sample size arrays.\n",
            "metadata": {}
        },
        {
            "id": "646b3fb8",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "import matplotlib.pyplot as plt\nnobs_x = np.arange(5, 21)\nnobs_y = nobs_x\nn_observations = (nobs_x, nobs_y)\nres = stats.power(test, rvs, n_observations, significance=0.05)\nax = plt.subplot()\nax.plot(nobs_x, res.power)\nax.set_xlabel('Sample Size')\nax.set_ylabel('Simulated Power')\nax.set_title('Simulated Power of `ttest_ind` with Equal Sample Sizes')\nplt.show()",
            "outputs": []
        },
        {
            "id": "0025a9fd",
            "cell_type": "markdown",
            "source": "Alternatively, we can investigate the impact that effect size has on the power.\nIn this case, the effect size is the location of the distribution underlying\nthe second sample.\n",
            "metadata": {}
        },
        {
            "id": "9e25736f",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "n_observations = (10, 12)\nloc = np.linspace(0, 1, 20)\nrvs2 = lambda size, loc: rng.normal(loc=loc, size=size)\nrvs = (rvs1, rvs2)\nres = stats.power(test, rvs, n_observations, significance=0.05,\n                  kwargs={'loc': loc})\nax = plt.subplot()\nax.plot(loc, res.power)\nax.set_xlabel('Effect Size')\nax.set_ylabel('Simulated Power')\nax.set_title('Simulated Power of `ttest_ind`, Varying Effect Size')\nplt.show()",
            "outputs": []
        },
        {
            "id": "28d94b34",
            "cell_type": "markdown",
            "source": "We can also use `power` to estimate the Type I error rate (also referred to by the\nambiguous term \"size\") of a test and assess whether it matches the nominal level.\nFor example, the null hypothesis of `jarque_bera` is that the sample was drawn from\na distribution with the same skewness and kurtosis as the normal distribution. To\nestimate the Type I error rate, we can consider the null hypothesis to be a true\n*alternative* hypothesis and calculate the power.\n",
            "metadata": {}
        },
        {
            "id": "caa853f9",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "test = stats.jarque_bera\nn_observations = 10\nrvs = rng.normal\nsignificance = np.linspace(0.0001, 0.1, 1000)\nres = stats.power(test, rvs, n_observations, significance=significance)\nsize = res.power",
            "outputs": []
        },
        {
            "id": "0aff8307",
            "cell_type": "markdown",
            "source": "As shown below, the Type I error rate of the test is far below the nominal level\nfor such a small sample, as mentioned in its documentation.\n",
            "metadata": {}
        },
        {
            "id": "fcf717cf",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "ax = plt.subplot()\nax.plot(significance, size)\nax.plot([0, 0.1], [0, 0.1], '--')\nax.set_xlabel('nominal significance level')\nax.set_ylabel('estimated test size (Type I error rate)')\nax.set_title('Estimated test size vs nominal significance level')\nax.set_aspect('equal', 'box')\nax.legend(('`ttest_1samp`', 'ideal test'))\nplt.show()",
            "outputs": []
        },
        {
            "id": "7322ca43",
            "cell_type": "markdown",
            "source": "As one might expect from such a conservative test, the power is quite low with\nrespect to some alternatives. For example, the power of the test under the\nalternative that the sample was drawn from the Laplace distribution may not\nbe much greater than the Type I error rate.\n",
            "metadata": {}
        },
        {
            "id": "0deb549f",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "rvs = rng.laplace\nsignificance = np.linspace(0.0001, 0.1, 1000)\nres = stats.power(test, rvs, n_observations, significance=0.05)\nprint(res.power)",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "0.0587"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "43230e94",
            "cell_type": "markdown",
            "source": "This is not a mistake in SciPy's implementation; it is simply due to the fact\nthat the null distribution of the test statistic is derived under the assumption\nthat the sample size is large (i.e. approaches infinity), and this asymptotic\napproximation is not accurate for small samples. In such cases, resampling\nand Monte Carlo methods (e.g. `permutation_test`, `goodness_of_fit`,\n`monte_carlo_test`) may be more appropriate.\n",
            "metadata": {}
        }
    ]
}