{
    "nbformat": 4,
    "nbformat_minor": 5,
    "metadata": {
        "kernelspec": {
            "display_name": "Python",
            "language": "python",
            "name": "python"
        },
        "language_info": {
            "name": "python"
        }
    },
    "cells": [
        {
            "id": "6b11a84a",
            "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": "eb6f4ea1",
            "cell_type": "markdown",
            "source": "A well-known test of the null hypothesis that data were drawn from a\ngiven distribution is the Kolmogorov-Smirnov (KS) test, available in SciPy\nas `scipy.stats.ks_1samp`. Suppose we wish to test whether the following\ndata:\n",
            "metadata": {}
        },
        {
            "id": "e8d55d71",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "import numpy as np\nfrom scipy import stats\nrng = np.random.default_rng()\nx = stats.uniform.rvs(size=75, random_state=rng)",
            "outputs": []
        },
        {
            "id": "71475537",
            "cell_type": "markdown",
            "source": "were sampled from a normal distribution. To perform a KS test, the\nempirical distribution function of the observed data will be compared\nagainst the (theoretical) cumulative distribution function of a normal\ndistribution. Of course, to do this, the normal distribution under the null\nhypothesis must be fully specified. This is commonly done by first fitting\nthe ``loc`` and ``scale`` parameters of the distribution to the observed\ndata, then performing the test.\n",
            "metadata": {}
        },
        {
            "id": "4600d20a",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "loc, scale = np.mean(x), np.std(x, ddof=1)\ncdf = stats.norm(loc, scale).cdf\nstats.ks_1samp(x, cdf)",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "KstestResult(statistic=0.1119257570456813,\n             pvalue=0.2827756409939257,\n             statistic_location=0.7751845155861765,\n             statistic_sign=-1)"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "3d648ed8",
            "cell_type": "markdown",
            "source": "An advantage of the KS-test is that the p-value - the probability of\nobtaining a value of the test statistic under the null hypothesis as\nextreme as the value obtained from the observed data - can be calculated\nexactly and efficiently. `goodness_of_fit` can only approximate these\nresults.\n",
            "metadata": {}
        },
        {
            "id": "812ace75",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "known_params = {'loc': loc, 'scale': scale}\nres = stats.goodness_of_fit(stats.norm, x, known_params=known_params,\n                            statistic='ks', rng=rng)\nres.statistic, res.pvalue",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "(0.1119257570456813, 0.2788)"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "bde2dcf7",
            "cell_type": "markdown",
            "source": "The statistic matches exactly, but the p-value is estimated by forming\na \"Monte Carlo null distribution\", that is, by explicitly drawing random\nsamples from `scipy.stats.norm` with the provided parameters and\ncalculating the stastic for each. The fraction of these statistic values\nat least as extreme as ``res.statistic`` approximates the exact p-value\ncalculated by `scipy.stats.ks_1samp`.\n\nHowever, in many cases, we would prefer to test only that the data were\nsampled from one of *any* member of the normal distribution family, not\nspecifically from the normal distribution with the location and scale\nfitted to the observed sample. In this case, Lilliefors [6] argued that\nthe KS test is far too conservative (that is, the p-value overstates\nthe actual probability of rejecting a true null hypothesis) and thus lacks\npower - the ability to reject the null hypothesis when the null hypothesis\nis actually false.\nIndeed, our p-value above is approximately 0.28, which is far too large\nto reject the null hypothesis at any common significance level.\n\nConsider why this might be. Note that in the KS test above, the statistic\nalways compares data against the CDF of a normal distribution fitted to the\n*observed data*. This tends to reduce the value of the statistic for the\nobserved data, but it is \"unfair\" when computing the statistic for other\nsamples, such as those we randomly draw to form the Monte Carlo null\ndistribution. It is easy to correct for this: whenever we compute the KS\nstatistic of a sample, we use the CDF of a normal distribution fitted\nto *that sample*. The null distribution in this case has not been\ncalculated exactly and is tyically approximated using Monte Carlo methods\nas described above. This is where `goodness_of_fit` excels.\n",
            "metadata": {}
        },
        {
            "id": "597bf9e0",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "res = stats.goodness_of_fit(stats.norm, x, statistic='ks',\n                            rng=rng)\nres.statistic, res.pvalue",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "(0.1119257570456813, 0.0196)"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "326a83fc",
            "cell_type": "markdown",
            "source": "Indeed, this p-value is much smaller, and small enough to (correctly)\nreject the null hypothesis at common significance levels, including 5% and\n2.5%.\n\nHowever, the KS statistic is not very sensitive to all deviations from\nnormality. The original advantage of the KS statistic was the ability\nto compute the null distribution theoretically, but a more sensitive\nstatistic - resulting in a higher test power - can be used now that we can\napproximate the null distribution\ncomputationally. The Anderson-Darling statistic [1] tends to be more\nsensitive, and critical values of the this statistic have been tabulated\nfor various significance levels and sample sizes using Monte Carlo methods.\n",
            "metadata": {}
        },
        {
            "id": "0953b0a7",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "res = stats.anderson(x, 'norm')\nprint(res.statistic)",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "1.2139573337497467"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "42864a4c",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "print(res.critical_values)",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "[0.549 0.625 0.75  0.875 1.041]"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "5ccfdba3",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "print(res.significance_level)",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "[15.  10.   5.   2.5  1. ]"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "1b520605",
            "cell_type": "markdown",
            "source": "Here, the observed value of the statistic exceeds the critical value\ncorresponding with a 1% significance level. This tells us that the p-value\nof the observed data is less than 1%, but what is it? We could interpolate\nfrom these (already-interpolated) values, but `goodness_of_fit` can\nestimate it directly.\n",
            "metadata": {}
        },
        {
            "id": "945ed209",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "res = stats.goodness_of_fit(stats.norm, x, statistic='ad',\n                            rng=rng)\nres.statistic, res.pvalue",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "(1.2139573337497467, 0.0034)"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "98e57f5b",
            "cell_type": "markdown",
            "source": "A further advantage is that use of `goodness_of_fit` is not limited to\na particular set of distributions or conditions on which parameters\nare known versus which must be estimated from data. Instead,\n`goodness_of_fit` can estimate p-values relatively quickly for any\ndistribution with a sufficiently fast and reliable ``fit`` method. For\ninstance, here we perform a goodness of fit test using the Cramer-von Mises\nstatistic against the Rayleigh distribution with known location and unknown\nscale.\n",
            "metadata": {}
        },
        {
            "id": "fa398aac",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "rng = np.random.default_rng()\nx = stats.chi(df=2.2, loc=0, scale=2).rvs(size=1000, random_state=rng)\nres = stats.goodness_of_fit(stats.rayleigh, x, statistic='cvm',\n                            known_params={'loc': 0}, rng=rng)",
            "outputs": []
        },
        {
            "id": "fe6b96ad",
            "cell_type": "markdown",
            "source": "This executes fairly quickly, but to check the reliability of the ``fit``\nmethod, we should inspect the fit result.\n",
            "metadata": {}
        },
        {
            "id": "824dfc0a",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "res.fit_result  # location is as specified, and scale is reasonable",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "  params: FitParams(loc=0.0, scale=2.1026719844231243)\n success: True\n message: 'The fit was performed successfully.'"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "d3a03862",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "import matplotlib.pyplot as plt  # matplotlib must be installed to plot\nres.fit_result.plot()\nplt.show()",
            "outputs": []
        },
        {
            "id": "7f818580",
            "cell_type": "markdown",
            "source": "If the distribution is not fit to the observed data as well as possible,\nthe test may not control the type I error rate, that is, the chance of\nrejecting the null hypothesis even when it is true.\n\nWe should also look for extreme outliers in the null distribution that\nmay be caused by unreliable fitting. These do not necessarily invalidate\nthe result, but they tend to reduce the test's power.\n",
            "metadata": {}
        },
        {
            "id": "56d42252",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "_, ax = plt.subplots()\nax.hist(np.log10(res.null_distribution))\nax.set_xlabel(\"log10 of CVM statistic under the null hypothesis\")\nax.set_ylabel(\"Frequency\")\nax.set_title(\"Histogram of the Monte Carlo null distribution\")\nplt.show()",
            "outputs": []
        },
        {
            "id": "f963a018",
            "cell_type": "markdown",
            "source": "This plot seems reassuring.\n\nIf ``fit`` method is working reliably, and if the distribution of the test\nstatistic is not particularly sensitive to the values of the fitted\nparameters, then the p-value provided by `goodness_of_fit` is expected to\nbe a good approximation.\n",
            "metadata": {}
        },
        {
            "id": "8df46801",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "res.statistic, res.pvalue",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "(0.2231991510248692, 0.0525)"
                    },
                    "execution_count": null
                }
            ]
        }
    ]
}