{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Inference with pyLFI\nExample usage of the ``pyLFI`` for parameter identification.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\nimport numpy as np\nimport pylfi\nimport scipy.stats as stats"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "In the following, we demonstrate ``pyLFI`` on a toy example. We will infer\nthe model parameters of a univariate Gaussian distribution: the mean\n$\\mu$ and standard deviation $\\sigma$. In this toy example, the\nlikelihood is\n$p (y_\\mathrm{obs} \\mid \\mu, \\sigma) = \\mathrm{N(\\mu=163, \\sigma=15)}$,\nand the observed data are sampled from the likelihood:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "mu_true = 163\nsigma_true = 15\nlikelihood = stats.norm(loc=mu_true, scale=sigma_true)\nobs_data = likelihood.rvs(size=1000)\n\nx = np.linspace(103, 223, 1000)\nlikelihood_pdf = likelihood.pdf(x)\nfig, ax = plt.subplots(figsize=(8, 4), tight_layout=True)\npylfi.utils.densityplot(x, likelihood_pdf, ax=ax, label='Likelihood')\npylfi.utils.rugplot(obs_data, pos=-0.0005, ax=ax, label='Observed data')\nax.set(xlabel='x', ylabel='Density')\nax.legend()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We assume that the likelihood is unknown, and formulate a model to describe\nthe observed data. The model needs to be implemented as a Python `callable`,\ni.e., a function or a ``__call__`` method in a class, that is parametrized by\nthe unknown model parameters we aim to infer, here $\\mu$ and\n$\\sigma$:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def simulator(mu, sigma, size=1000):\n    y_sim = stats.norm(loc=mu, scale=sigma).rvs(size=size)\n    return y_sim"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Next, we need to reduce the raw data into low-dimensional summary statistics.\nThe summary statistics calculator also needs to be implemented as a Python\n`callable`. The function must return the summary statistics as a Python\n`list` or `numpy.ndarray`. Here, we take the mean and standard deviation to\nbe summary statistics of the data (these are actually sufficient summary\nstatistics):\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def stat_calc(y):\n    sum_stats = [numpy.mean(y), numpy.std(y)]\n    return sum_stats"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We then place priors over the unknown model parameters using the `.Prior`\nclass. In the present example, we define the priors:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "mu_prior = pylfi.Prior('norm',\n                       loc=165,\n                       scale=2,\n                       name='mu',\n                       tex=r'$\\mu$'\n                       )\n\nsigma_prior = pylfi.Prior('uniform',\n                          loc=12,\n                          scale=7,\n                          name='sigma',\n                          tex=r'$\\sigma$'\n                          )\n\npriors = [mu_prior, sigma_prior]\n\nfig, axes = plt.subplots(nrows=2, figsize=(8, 4), tight_layout=True)\nx = np.linspace(159, 171, 1000)\nmu_prior.plot_prior(x, ax=axes[0])\nx = np.linspace(11, 20, 1000)\nsigma_prior.plot_prior(x, color='C1', facecolor='wheat', ax=axes[1])"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "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.6"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}