{
    "nbformat": 4,
    "nbformat_minor": 5,
    "metadata": {
        "kernelspec": {
            "display_name": "Python",
            "language": "python",
            "name": "python"
        },
        "language_info": {
            "name": "python"
        }
    },
    "cells": [
        {
            "id": "3eae0e04",
            "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": "edf35a4d",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "from scipy.stats.sampling import DiscreteGuideTable\nimport numpy as np",
            "outputs": []
        },
        {
            "id": "f28446a7",
            "cell_type": "markdown",
            "source": "To create a random number generator using a probability vector, use:\n",
            "metadata": {}
        },
        {
            "id": "ee360cf3",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "pv = [0.1, 0.3, 0.6]\nurng = np.random.default_rng()\nrng = DiscreteGuideTable(pv, random_state=urng)",
            "outputs": []
        },
        {
            "id": "b2da05f4",
            "cell_type": "markdown",
            "source": "The RNG has been setup. Now, we can now use the `rvs` method to\ngenerate samples from the distribution:\n",
            "metadata": {}
        },
        {
            "id": "ef7031e3",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "rvs = rng.rvs(size=1000)",
            "outputs": []
        },
        {
            "id": "3ecf4189",
            "cell_type": "markdown",
            "source": "To verify that the random variates follow the given distribution, we can\nuse the chi-squared test (as a measure of goodness-of-fit):\n",
            "metadata": {}
        },
        {
            "id": "064163ca",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "from scipy.stats import chisquare\n_, freqs = np.unique(rvs, return_counts=True)\nfreqs = freqs / np.sum(freqs)\nfreqs",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([0.092, 0.355, 0.553])"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "71b32315",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "chisquare(freqs, pv).pvalue",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "0.9987382966178464"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "75ab7114",
            "cell_type": "markdown",
            "source": "As the p-value is very high, we fail to reject the null hypothesis that\nthe observed frequencies are the same as the expected frequencies. Hence,\nwe can safely assume that the variates have been generated from the given\ndistribution. Note that this just gives the correctness of the algorithm\nand not the quality of the samples.\n\nIf a PV is not available, an instance of a class with a PMF method and a\nfinite domain can also be passed.\n",
            "metadata": {}
        },
        {
            "id": "3f9504b9",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "urng = np.random.default_rng()\nfrom scipy.stats import binom\nn, p = 10, 0.2\ndist = binom(n, p)\nrng = DiscreteGuideTable(dist, random_state=urng)",
            "outputs": []
        },
        {
            "id": "dfbc0b51",
            "cell_type": "markdown",
            "source": "Now, we can sample from the distribution using the `rvs` method\nand also measure the goodness-of-fit of the samples:\n",
            "metadata": {}
        },
        {
            "id": "7a3939e3",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "rvs = rng.rvs(1000)\n_, freqs = np.unique(rvs, return_counts=True)\nfreqs = freqs / np.sum(freqs)\nobs_freqs = np.zeros(11)  # some frequencies may be zero.\nobs_freqs[:freqs.size] = freqs\npv = [dist.pmf(i) for i in range(0, 11)]\npv = np.asarray(pv) / np.sum(pv)\nchisquare(obs_freqs, pv).pvalue",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "0.9999999999999989"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "c8840152",
            "cell_type": "markdown",
            "source": "To check that the samples have been drawn from the correct distribution,\nwe can visualize the histogram of the samples:\n",
            "metadata": {}
        },
        {
            "id": "27b9b659",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "import matplotlib.pyplot as plt\nrvs = rng.rvs(1000)\nfig = plt.figure()\nax = fig.add_subplot(111)\nx = np.arange(0, n+1)\nfx = dist.pmf(x)\nfx = fx / fx.sum()\nax.plot(x, fx, 'bo', label='true distribution')\nax.vlines(x, 0, fx, lw=2)\nax.hist(rvs, bins=np.r_[x, n+1]-0.5, density=True, alpha=0.5,\n        color='r', label='samples')\nax.set_xlabel('x')\nax.set_ylabel('PMF(x)')\nax.set_title('Discrete Guide Table Samples')\nplt.legend()\nplt.show()",
            "outputs": []
        },
        {
            "id": "d2dc8011",
            "cell_type": "markdown",
            "source": "To set the size of the guide table use the `guide_factor` keyword argument.\nThis sets the size of the guide table relative to the probability vector\n",
            "metadata": {}
        },
        {
            "id": "9dcdc9f9",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "rng = DiscreteGuideTable(pv, guide_factor=1, random_state=urng)",
            "outputs": []
        },
        {
            "id": "d91af9bc",
            "cell_type": "markdown",
            "source": "To calculate the PPF of a binomial distribution with $n=4$ and\n$p=0.1$: we can set up a guide table as follows:\n",
            "metadata": {}
        },
        {
            "id": "5dfc22d4",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "n, p = 4, 0.1\ndist = binom(n, p)\nrng = DiscreteGuideTable(dist, random_state=42)\nrng.ppf(0.5)",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "0.0"
                    },
                    "execution_count": null
                }
            ]
        }
    ]
}