docs/usage/binned-fit.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": true,
    "hideOutput": true,
    "hidePrompt": true,
    "jupyter": {
     "source_hidden": true
    },
    "slideshow": {
     "slide_type": "skip"
    },
    "tags": [
     "remove-cell",
     "skip-execution"
    ]
   },
   "outputs": [],
   "source": [
    "# WARNING: advised to install a specific version, e.g. tensorwaves==0.1.2\n",
    "%pip install -q tensorwaves[doc,jax,pwa,viz] IPython"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": true,
    "hideOutput": true,
    "hidePrompt": true,
    "jupyter": {
     "source_hidden": true
    },
    "slideshow": {
     "slide_type": "skip"
    },
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [],
   "source": [
    "import os\n",
    "\n",
    "STATIC_WEB_PAGE = {\"EXECUTE_NB\", \"READTHEDOCS\"}.intersection(os.environ)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```{autolink-concat}\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Binned fit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Imagine we have the following data distribution over $x$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "sample_size = 1_000\n",
    "rng = np.random.default_rng(seed=0)\n",
    "x_distribution = rng.power(3, size=sample_size)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ":::{margin}\n",
    "\n",
    "Alternatively, we could have sampled the distribution with {func}`numpy.histogram`.\n",
    "\n",
    ":::\n",
    "\n",
    "We sample the distribution into a histogram with {func}`matplotlib.pyplot.hist`. The bin edges will serve as $x$-values and the bin heights will be the observed $y$-values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "%config InlineBackend.figure_formats = ['svg']\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(8, 5))\n",
    "ax.set_xlabel(\"$x$\")\n",
    "n_bins = 50\n",
    "bin_values, bin_edges, _ = ax.hist(x_distribution, bins=n_bins)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "x_values = (bin_edges[1:] + bin_edges[:-1]) / 2\n",
    "y_values = bin_values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "The data distribution has been generated by {func}`numpy.random.power` and can therefore be described by the following expression:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import sympy as sp\n",
    "\n",
    "n, a, x = sp.symbols(\"n a x\")\n",
    "expression = n * x ** (a - 1)\n",
    "expression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We would like to find which values for $a, n$ describe this distribution best. For fast computations, we first formulate this expression as a {class}`.ParametrizedFunction` with {mod}`numpy` as computational backend ({mod}`numpy` is fast enough for the small data samples with which we are working here). We also provide a first guess for the values of $a$ and $n$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from tensorwaves.function.sympy import create_parametrized_function\n",
    "\n",
    "function = create_parametrized_function(\n",
    "    expression,\n",
    "    parameters={a: 2, n: 15},\n",
    "    backend=\"numpy\",\n",
    ")\n",
    "initial_parameters = function.parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "function.update_parameters(initial_parameters)\n",
    "fig, ax = plt.subplots(figsize=(8, 5))\n",
    "ax.set_xlabel(\"$x$\")\n",
    "ax.hist(x_distribution, bins=n_bins, label=\"Data distribution\")\n",
    "ax.plot(x_values, function({\"x\": x_values}), label=\"Initial fit model\", c=\"red\")\n",
    "ax.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This **binned distribution** lends itself well to be optimized with a {class}`.ChiSquared` estimator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from tensorwaves.estimator import ChiSquared\n",
    "from tensorwaves.optimizer import Minuit2\n",
    "\n",
    "estimator = ChiSquared(\n",
    "    function,\n",
    "    domain={\"x\": x_values},\n",
    "    observed_values=y_values,\n",
    "    backend=\"numpy\",\n",
    ")\n",
    "optimizer = Minuit2()\n",
    "fit_result = optimizer.optimize(estimator, initial_parameters)\n",
    "fit_result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [],
   "source": [
    "assert fit_result.minimum_valid"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The optimized value of $a$ indeed lies close to the value of $a=3$ with which we generated the distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "function.update_parameters(fit_result.parameter_values)\n",
    "fig, ax = plt.subplots(figsize=(8, 5))\n",
    "ax.set_xlabel(\"$x$\")\n",
    "ax.hist(x_distribution, bins=n_bins, label=\"Data distribution\")\n",
    "ax.plot(x_values, function({\"x\": x_values}), label=\"Optimized model\", c=\"red\")\n",
    "ax.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Also note that the integral over the expression with these optimized parameter values is close to the surface of the bins:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "bin_width = bin_edges[1] - bin_edges[0]\n",
    "surface = bin_values.sum() * bin_width\n",
    "surface"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "from IPython.display import Math\n",
    "\n",
    "substitutions = {\n",
    "    a: fit_result.parameter_values[\"a\"],\n",
    "    n: fit_result.parameter_values[\"n\"],\n",
    "}\n",
    "integral = sp.Integral(expression, (x, 0, 1))\n",
    "evaluated_integral = integral.subs(substitutions).doit()\n",
    "latex = sp.multiline_latex(integral, evaluated_integral, environment=\"eqnarray\")\n",
    "Math(latex)"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "toc_visible": true
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.17"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}