{
    "nbformat": 4,
    "nbformat_minor": 5,
    "metadata": {
        "kernelspec": {
            "display_name": "Python",
            "language": "python",
            "name": "python"
        },
        "language_info": {
            "name": "python"
        }
    },
    "cells": [
        {
            "id": "779c5853",
            "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": "08ec86aa",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "from scipy.stats.sampling import TransformedDensityRejection\nimport numpy as np",
            "outputs": []
        },
        {
            "id": "1c3394e8",
            "cell_type": "markdown",
            "source": "Suppose we have a density:\n\n$$ f(x) = \\begin{cases} 1 - x^2,  &  -1 \\leq x \\leq 1 \\\\ 0,        &  \\text{otherwise} \\end{cases} $$\nThe derivative of this density function is:\n\n$$ \\frac{df(x)}{dx} = \\begin{cases} -2x,  &  -1 \\leq x \\leq 1 \\\\ 0,    &  \\text{otherwise} \\end{cases} $$\nNotice that the PDF doesn't integrate to 1. As this is a rejection based\nmethod, we need not have a normalized PDF. To initialize the generator,\nwe can use:\n",
            "metadata": {}
        },
        {
            "id": "e7d6bea7",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "urng = np.random.default_rng()\nclass MyDist:\n    def pdf(self, x):\n        return 1-x*x\n    def dpdf(self, x):\n        return -2*x\n\ndist = MyDist()\nrng = TransformedDensityRejection(dist, domain=(-1, 1),\n                                  random_state=urng)",
            "outputs": []
        },
        {
            "id": "b9766a54",
            "cell_type": "markdown",
            "source": "Domain can be very useful to truncate the distribution but to avoid passing\nit every time to the constructor, a default domain can be set by providing a\n`support` method in the distribution object (`dist`):\n",
            "metadata": {}
        },
        {
            "id": "e583b48d",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "class MyDist:\n    def pdf(self, x):\n        return 1-x*x\n    def dpdf(self, x):\n        return -2*x\n    def support(self):\n        return (-1, 1)\n\ndist = MyDist()\nrng = TransformedDensityRejection(dist, random_state=urng)",
            "outputs": []
        },
        {
            "id": "4bdd1fc9",
            "cell_type": "markdown",
            "source": "Now, we can use the `rvs` method to generate samples from the distribution:\n",
            "metadata": {}
        },
        {
            "id": "f5594c1c",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "rvs = rng.rvs(1000)",
            "outputs": []
        },
        {
            "id": "35a41492",
            "cell_type": "markdown",
            "source": "We can check that the samples are from the given distribution by visualizing\nits histogram:\n",
            "metadata": {}
        },
        {
            "id": "6e4a86ed",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "import matplotlib.pyplot as plt\nx = np.linspace(-1, 1, 1000)\nfx = 3/4 * dist.pdf(x)  # 3/4 is the normalizing constant\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('Transformed Density Rejection Samples')\nplt.legend()\nplt.show()",
            "outputs": []
        }
    ]
}