{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Numerical integration\n",
    "## Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from scipy import stats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Fixed parameters\n",
    "I = 1e9  # mm^4\n",
    "L = 10_000  # mm\n",
    "Q = 100  # kN\n",
    "\n",
    "# Measurements\n",
    "d_meas = 50  # mm\n",
    "sigma_meas =  10  # mm\n",
    "\n",
    "# Prior\n",
    "E_mean = 60  # GPa\n",
    "E_std = 20  # GPa"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Forward model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def midspan_deflection(E):\n",
    "    return Q * L ** 3 / (48 * E * I)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bayesian functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def likelihood(E):\n",
    "    return ...\n",
    "\n",
    "\n",
    "def prior(E):\n",
    "    return ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Perform numerical integration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "E_values = np.linspace(-20, 140, 1000)\n",
    "\n",
    "# Calculate prior and likelihood values\n",
    "prior_values = ...\n",
    "likelihood_values = ...\n",
    "\n",
    "evidence = ... # use the trapezoidal rule to integrate\n",
    "posterior_values = ... "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Posterior summary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Mean\n",
    "mean = ... # use the trapezoidal rule to integrate\n",
    "print(f'Mean: {mean:.2f} GPa')\n",
    "\n",
    "# Median\n",
    "cdf = np.cumsum(posterior_values) * np.diff(E_values, prepend=0)\n",
    "median = E_values[np.argmin(np.abs(cdf - 0.5))]\n",
    "print(f'Median: {median:.2f} GPa')\n",
    "\n",
    "# Standard deviation\n",
    "std = ...  # use the trapezoidal rule to integrate\n",
    "print(f'Standard deviation: {std:.2f} GPa')\n",
    "\n",
    "# 5th and 95th percentiles\n",
    "percentiles = np.interp([0.05, 0.95], cdf, E_values)\n",
    "print(f'5th percentile: {percentiles[0]:.2f} GPa')\n",
    "print(f'95th percentile: {percentiles[1]:.2f} GPa')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(E_values, prior_values, label='Prior')\n",
    "plt.plot(E_values, likelihood_values, label='Likelihood')\n",
    "plt.plot(E_values, posterior_values, label='Posterior')\n",
    "plt.xlabel('E (GPa)')\n",
    "plt.ylabel('Density')\n",
    "plt.ylim([0, 0.045])\n",
    "plt.legend()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "venv",
   "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.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}