{
    "nbformat": 4,
    "nbformat_minor": 5,
    "metadata": {
        "kernelspec": {
            "display_name": "Python",
            "language": "python",
            "name": "python"
        },
        "language_info": {
            "name": "python"
        }
    },
    "cells": [
        {
            "id": "20466fff",
            "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": "0c36e6b5",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "from scipy.stats.sampling import NumericalInverseHermite\nfrom scipy.stats import norm, genexpon\nfrom scipy.special import ndtr\nimport numpy as np",
            "outputs": []
        },
        {
            "id": "c25a5e43",
            "cell_type": "markdown",
            "source": "To create a generator to sample from the standard normal distribution, do:\n",
            "metadata": {}
        },
        {
            "id": "da8cddb0",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "class StandardNormal:\n    def pdf(self, x):\n       return 1/np.sqrt(2*np.pi) * np.exp(-x**2 / 2)\n    def cdf(self, x):\n       return ndtr(x)\n\ndist = StandardNormal()\nurng = np.random.default_rng()\nrng = NumericalInverseHermite(dist, random_state=urng)",
            "outputs": []
        },
        {
            "id": "d6a79ae1",
            "cell_type": "markdown",
            "source": "The `NumericalInverseHermite` has a method that approximates the PPF of the\ndistribution.\n",
            "metadata": {}
        },
        {
            "id": "553f7c87",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "rng = NumericalInverseHermite(dist)\np = np.linspace(0.01, 0.99, 99) # percentiles from 1% to 99%\nnp.allclose(rng.ppf(p), norm.ppf(p))",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "True"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "4dc48bb0",
            "cell_type": "markdown",
            "source": "Depending on the implementation of the distribution's random sampling\nmethod, the random variates generated may be nearly identical, given\nthe same random state.\n",
            "metadata": {}
        },
        {
            "id": "f439de3d",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "dist = genexpon(9, 16, 3)\nrng = NumericalInverseHermite(dist)\n# `seed` ensures identical random streams are used by each `rvs` method\nseed = 500072020\nrvs1 = dist.rvs(size=100, random_state=np.random.default_rng(seed))\nrvs2 = rng.rvs(size=100, random_state=np.random.default_rng(seed))\nnp.allclose(rvs1, rvs2)",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "True"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "1c2f05d4",
            "cell_type": "markdown",
            "source": "To check that the random variates closely follow the given distribution, we can\nlook at its histogram:\n",
            "metadata": {}
        },
        {
            "id": "3c234b28",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "import matplotlib.pyplot as plt\ndist = StandardNormal()\nrng = NumericalInverseHermite(dist)\nrvs = rng.rvs(10000)\nx = np.linspace(rvs.min()-0.1, rvs.max()+0.1, 1000)\nfx = norm.pdf(x)\nplt.plot(x, fx, 'r-', lw=2, label='true distribution')\nplt.hist(rvs, bins=20, density=True, alpha=0.8, label='random variates')\nplt.xlabel('x')\nplt.ylabel('PDF(x)')\nplt.title('Numerical Inverse Hermite Samples')\nplt.legend()\nplt.show()",
            "outputs": []
        },
        {
            "id": "7e08273b",
            "cell_type": "markdown",
            "source": "Given the derivative of the PDF w.r.t the variate (i.e. ``x``), we can use\nquintic Hermite interpolation to approximate the PPF by passing the `order`\nparameter:\n",
            "metadata": {}
        },
        {
            "id": "459e38fe",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "class StandardNormal:\n    def pdf(self, x):\n       return 1/np.sqrt(2*np.pi) * np.exp(-x**2 / 2)\n    def dpdf(self, x):\n       return -1/np.sqrt(2*np.pi) * x * np.exp(-x**2 / 2)\n    def cdf(self, x):\n       return ndtr(x)\n\ndist = StandardNormal()\nurng = np.random.default_rng()\nrng = NumericalInverseHermite(dist, order=5, random_state=urng)",
            "outputs": []
        },
        {
            "id": "1259144b",
            "cell_type": "markdown",
            "source": "Higher orders result in a fewer number of intervals:\n",
            "metadata": {}
        },
        {
            "id": "db6526bc",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "rng3 = NumericalInverseHermite(dist, order=3)\nrng5 = NumericalInverseHermite(dist, order=5)\nrng3.intervals, rng5.intervals",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "(3000, 522)"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "efa7554f",
            "cell_type": "markdown",
            "source": "The u-error can be estimated by calling the `u_error` method. It runs a small\nMonte-Carlo simulation to estimate the u-error. By default, 100,000 samples are\nused. This can be changed by passing the `sample_size` argument:\n",
            "metadata": {}
        },
        {
            "id": "19b1b9cb",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "rng1 = NumericalInverseHermite(dist, u_resolution=1e-10)\nrng1.u_error(sample_size=1000000)  # uses one million samples",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "UError(max_error=9.53167544892608e-11, mean_absolute_error=2.2450136432146864e-11)"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "8c9f76df",
            "cell_type": "markdown",
            "source": "This returns a namedtuple which contains the maximum u-error and the mean\nabsolute u-error.\n\nThe u-error can be reduced by decreasing the u-resolution (maximum allowed u-error):\n",
            "metadata": {}
        },
        {
            "id": "bcaaef8b",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "rng2 = NumericalInverseHermite(dist, u_resolution=1e-13)\nrng2.u_error(sample_size=1000000)",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "UError(max_error=9.32027892364129e-14, mean_absolute_error=1.5194172675685075e-14)"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "9661948d",
            "cell_type": "markdown",
            "source": "Note that this comes at the cost of increased setup time and number of intervals.\n",
            "metadata": {}
        },
        {
            "id": "3318b7a5",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "rng1.intervals",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "1022"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "9f1cc2ba",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "rng2.intervals",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "5687"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "d9b4f8c4",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "from timeit import timeit\nf = lambda: NumericalInverseHermite(dist, u_resolution=1e-10)\ntimeit(f, number=1)",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "0.017409582000254886  # may vary"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "3fa7bdbe",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "f = lambda: NumericalInverseHermite(dist, u_resolution=1e-13)\ntimeit(f, number=1)",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "0.08671202100003939  # may vary"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "c886dbd8",
            "cell_type": "markdown",
            "source": "Since the PPF of the normal distribution is available as a special function, we\ncan also check the x-error, i.e. the difference between the approximated PPF and\nexact PPF\n```\n\n```",
            "metadata": {}
        },
        {
            "id": "2c38f263",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "import matplotlib.pyplot as plt\nu = np.linspace(0.01, 0.99, 1000)\napproxppf = rng.ppf(u)\nexactppf = norm.ppf(u)\nerror = np.abs(exactppf - approxppf)\nplt.plot(u, error)\nplt.xlabel('u')\nplt.ylabel('error')\nplt.title('Error between exact and approximated PPF (x-error)')\nplt.show()",
            "outputs": []
        }
    ]
}