{
    "nbformat": 4,
    "nbformat_minor": 5,
    "metadata": {
        "kernelspec": {
            "display_name": "Python",
            "language": "python",
            "name": "python"
        },
        "language_info": {
            "name": "python"
        }
    },
    "cells": [
        {
            "id": "e0e4ecf9",
            "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": "5d37c6f8",
            "cell_type": "markdown",
            "source": "Suppose we wish to fit a distribution to the following data.\n",
            "metadata": {}
        },
        {
            "id": "c2736ec8",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "import numpy as np\nfrom scipy import stats\nrng = np.random.default_rng()\ndist = stats.nbinom\nshapes = (5, 0.5)\ndata = dist.rvs(*shapes, size=1000, random_state=rng)",
            "outputs": []
        },
        {
            "id": "8e6066c5",
            "cell_type": "markdown",
            "source": "Suppose we do not know how the data were generated, but we suspect that\nit follows a negative binomial distribution with parameters *n* and *p*\\.\n(See `scipy.stats.nbinom`.) We believe that the parameter *n* was fewer\nthan 30, and we know that the parameter *p* must lie on the interval\n[0, 1]. We record this information in a variable `bounds` and pass\nthis information to `fit`.\n",
            "metadata": {}
        },
        {
            "id": "8e5b762f",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "bounds = [(0, 30), (0, 1)]\nres = stats.fit(dist, data, bounds)",
            "outputs": []
        },
        {
            "id": "fa95551a",
            "cell_type": "markdown",
            "source": "`fit` searches within the user-specified `bounds` for the\nvalues that best match the data (in the sense of maximum likelihood\nestimation). In this case, it found shape values similar to those\nfrom which the data were actually generated.\n",
            "metadata": {}
        },
        {
            "id": "4d11f339",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "res.params",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "FitParams(n=5.0, p=0.5028157644634368, loc=0.0)  # may vary"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "fb40b327",
            "cell_type": "markdown",
            "source": "We can visualize the results by superposing the probability mass function\nof the distribution (with the shapes fit to the data) over a normalized\nhistogram of the data.\n",
            "metadata": {}
        },
        {
            "id": "afe0fee5",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "import matplotlib.pyplot as plt  # matplotlib must be installed to plot\nres.plot()\nplt.show()",
            "outputs": []
        },
        {
            "id": "e79d41b5",
            "cell_type": "markdown",
            "source": "Note that the estimate for *n* was exactly integral; this is because\nthe domain of the `nbinom` PMF includes only integral *n*, and the `nbinom`\nobject \"knows\" that. `nbinom` also knows that the shape *p* must be a\nvalue between 0 and 1. In such a case - when the domain of the distribution\nwith respect to a parameter is finite - we are not required to specify\nbounds for the parameter.\n",
            "metadata": {}
        },
        {
            "id": "21a142e5",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "bounds = {'n': (0, 30)}  # omit parameter p using a `dict`\nres2 = stats.fit(dist, data, bounds)\nres2.params",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "FitParams(n=5.0, p=0.5016492009232932, loc=0.0)  # may vary"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "249b5421",
            "cell_type": "markdown",
            "source": "If we wish to force the distribution to be fit with *n* fixed at 6, we can\nset both the lower and upper bounds on *n* to 6. Note, however, that the\nvalue of the objective function being optimized is typically worse (higher)\nin this case.\n",
            "metadata": {}
        },
        {
            "id": "13659f7f",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "bounds = {'n': (6, 6)}  # fix parameter `n`\nres3 = stats.fit(dist, data, bounds)\nres3.params",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "FitParams(n=6.0, p=0.5486556076755706, loc=0.0)  # may vary"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "7fec33bc",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "res3.nllf() > res.nllf()",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "True  # may vary"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "5f1e5f80",
            "cell_type": "markdown",
            "source": "Note that the numerical results of the previous examples are typical, but\nthey may vary because the default optimizer used by `fit`,\n`scipy.optimize.differential_evolution`, is stochastic. However, we can\ncustomize the settings used by the optimizer to ensure reproducibility -\nor even use a different optimizer entirely - using the `optimizer`\nparameter.\n",
            "metadata": {}
        },
        {
            "id": "22e49245",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "from scipy.optimize import differential_evolution\nrng = np.random.default_rng()\ndef optimizer(fun, bounds, *, integrality):\n    return differential_evolution(fun, bounds, strategy='best2bin',\n                                  rng=rng, integrality=integrality)\nbounds = [(0, 30), (0, 1)]\nres4 = stats.fit(dist, data, bounds, optimizer=optimizer)\nres4.params",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "FitParams(n=5.0, p=0.5015183149259951, loc=0.0)"
                    },
                    "execution_count": null
                }
            ]
        }
    ]
}