Notation-and-derivations.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Notation and derivations for Pints\n",
    "\n",
    "This document is:\n",
    "\n",
    "1. An attempt to write down some standard notation for Pints\n",
    "2. A place to stick derivations of stuff used in Pints (e.g. loglikelihoods, not optimiser or sampling algorithms).\n",
    "3. A place to work out how to get loglikelihoods into NumPy, and document what you've done.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## List of symbols\n",
    "\n",
    "Please stick to the symbols below when editing this document.\n",
    "\n",
    "| Description                 | Mathematical version     | Code version          |\n",
    "|-----------------------------|--------------------------|-----------------------|\n",
    "| Time                        | $t$, $t_i$               | `times`               |\n",
    "| Observations (values)       | $v$, $v_i$               | `values`              |\n",
    "| Time series (data)          | $D$                      |                       |\n",
    "| Number of times/values      | $n_t$                    | `n_times`             |\n",
    "| Parameters                  | $x$, $x_i$               | `x`, `parameters`     |\n",
    "| Number of parameters        | $n_p$                    | `n_parameters`, `n_p` |\n",
    "| Forward model               | $m(t|x)$                 | `model`               |\n",
    "| Forward model values        | $m_i(x) = m(t_i|x)$      | `y`                   |\n",
    "| Number of model outputs     | $n_o$                    | `n_outputs`, `n_o`    |\n",
    "| PDF                         | $f(D|x)$                 |                       |\n",
    "| Likelihood                  | $l(x|D)$ or $l(x)$       |                       |\n",
    "| LogPDF                      |                          | `logpdf`              |\n",
    "| Loglikelihood               | $L(x)$                   | L                     |\n",
    "\n",
    "Prior?\n",
    "Logposterior?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Problem statement\n",
    "\n",
    "We have some noisy time-series data and a forward model (simulation) that can be used to replicate it.\n",
    "We'd like to find out which parameter values are compatible with the experimental evidence.\n",
    "\n",
    "- Observations $D = \\{(t_1, v_1),...,(t_{n_t}, v_{n_t})\\}$ where $v_i$ is the experimental measurement at time $t_i$, and the times are ordered $t_{i + 1} > t_i$, with $i = 1, 2, ..., n_t$.\n",
    "\n",
    "- A forward model $m(t|x)$ that takes a time $t$ as input, as well as a parameter vector $x$ of length $n_p$.\n",
    "\n",
    "- Observations can be scalars, or vectors of some fixed length $n_o$. If vector outputs are used, the forward model must also produce $n_o$ _outputs_. In general, this means $v \\in {\\rm I\\!R}^{n_o}$ and $m(t|x) \\to {\\rm I\\!R}^{n_o}$.\n",
    "\n",
    "Often, but not always\n",
    "\n",
    "- The parameters live in some bounded space $x \\in X \\subset{\\rm I\\!R}^{n_p}$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Noise/error models\n",
    "\n",
    "If we assume the model can perfectly describe the data with the correct parameters, we can interpret the remaining error as noise.\n",
    "If we have a probabilistic model for this noise, we can then write a probability density function (PDF) for the probability of a model with parameters $x$ generating the observations $D$:\n",
    "\n",
    "$$ f(D|x) $$\n",
    "\n",
    "We can rewrite this as a likelihood:\n",
    "\n",
    "$$ l(x|D) \\equiv f(D|x) $$\n",
    "\n",
    "As it turns out, it's usually easier to work with the natural logarithm of this function instead.\n",
    "\n",
    "$$ L(x|D) = \\log l(x|D)$$\n",
    "\n",
    "which we'll often shorten to\n",
    "\n",
    "$$ L(x) $$.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Gaussian independent noise\n",
    "\n",
    "In this section, we derive a loglikelihood for a parameter vector $x$ assuming Gaussian noise, that is independent from observation to observation.\n",
    "\n",
    "We start by looking at single-output models, and assume our noise is from a Gaussian distribution with mean 0 and standard distribution $\\sigma$.\n",
    "For now we'll assume we have some way of knowing $\\sigma$, e.g. by doing an independent measurement.\n",
    "\n",
    "We can then treat our observations as random variables of the form (model prediction + Gausian noise):\n",
    "\n",
    "$$ V_i \\sim m_i(x) + \\mathcal{N}(0, \\sigma^2) = \\mathcal{N}(m_i(x), \\sigma^2)$$\n",
    "\n",
    "Filling in the equation for the Gaussian distribution, we find\n",
    "\n",
    "$$ f_i(v_i | x) = \\frac{1}{\\sqrt{2\\pi\\sigma^2}} \\exp\n",
    "                  \\left( -\\frac{ \\left( m_i(x) - v_i \\right)^2}{2\\sigma^2} \\right) $$\n",
    "\n",
    "The independent noise assumption then gives\n",
    "\n",
    "$$ f(D | x) = \\prod_{i=1}^{n_t} \\frac{1}{\\sqrt{2\\pi\\sigma^2}} \\exp\n",
    "              \\left( -\\frac{ \\left( m_i(x) - v_i \\right)^2}{2\\sigma^2} \\right) $$\n",
    "\n",
    "To find $L(x)$ we use $L(x) = L(x|D) = \\log l(x|D) = \\log f(D|x)$:\n",
    "\n",
    "$$ L(x) = \n",
    "    - \\frac{n_t}{2} \\log(2\\pi)\n",
    "    - n_t \\log(\\sigma)\n",
    "    - \\frac{1}{2\\sigma^2} \\sum_{i=1}^{n_t} \\left(m_i(x) - v_i \\right)^2\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1 Multiple outputs\n",
    "\n",
    "To find a multi-output version, we define the observation in output $j$ at time $i$ as $v_{ij}$. The equivalent model variable is $m_{ij}(x)$.\n",
    "\n",
    "If we assume that the noise in any one output is independent from the noise in the others, we can write the total log-likelihood of the observations as the sum of the loglikelihoods in each output:\n",
    "\n",
    "$$ L(x) = \\sum_{j=1}^{n_o} L_j(x) $$\n",
    "$$ L(x) = \\sum_{j=1}^{n_o} - \\frac{n_t}{2} \\log(2\\pi)\n",
    "                           - n_t \\log(\\sigma)\n",
    "                           - \\frac{1}{2\\sigma^2} \\sum_{i=1}^{n_t} \\left(m_i(x) - v_i \\right)^2\n",
    "$$\n",
    "$$ L(x) = - \\frac{n_t n_o}{2} \\log(2\\pi) \n",
    "          - n_t \\sum_{j=1}^{n_o} \\log(\\sigma_j)\n",
    "          - \\sum_{j=1}^{n_o} \\left[ \\frac{1}{2\\sigma_j^2} \\sum_{i=1}^{n_t} \\left(m_{ij}(x) - v_{ij} \\right)^2 \\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2 Derivatives with respect to the parameters\n",
    "\n",
    "The partial derivatives with respect to parameter $x_k$, evaluated at $x$ are written as\n",
    "\n",
    "$$ \\left. \\frac{\\partial m_{ij}}{\\partial x_k} \\right|_x $$\n",
    "\n",
    "For the independent Gaussian loglikelihood, we then find\n",
    "\n",
    "$$ \\left. \\frac{\\partial L}{\\partial x_k} \\right|_x =\n",
    "    \\sum_{j=1}^{n_o} \\left[ \\frac{1}{\\sigma_j^2} \\sum_{i=1}^{n_t}\n",
    "    \\left(v_{ij} - m_{ij}(x) \\right) \\left. \\frac{\\partial m_{ij}}{\\partial x_k} \\right|_x \\right]\n",
    "$$\n",
    "\n",
    "NumPy implementation follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(n_times, n_outputs, n_parameters)\n",
      "(4, 3, 2)\n",
      "times\n",
      "(4,)\n",
      "[3 4 5 6]\n",
      "residuals\n",
      "(4, 3)\n",
      "[[ 3  1  2]\n",
      " [-1  3 -2]\n",
      " [-2  1  7]\n",
      " [ 3  4  5]]\n",
      "derivatives\n",
      "(4, 3, 2)\n",
      "[[[  0.   3.]\n",
      "  [  6.   9.]\n",
      "  [ 12.  15.]]\n",
      "\n",
      " [[ 18.  21.]\n",
      "  [ 24.  27.]\n",
      "  [ 30.  33.]]\n",
      "\n",
      " [[ 36.  39.]\n",
      "  [ 42.  45.]\n",
      "  [ 48.  51.]]\n",
      "\n",
      " [[ 54.  57.]\n",
      "  [ 60.  63.]\n",
      "  [ 66.  69.]]]\n",
      "\n",
      "L:\n",
      "-1933.61883471\n",
      "\n",
      "With sums, dL = \n",
      "23200.0\n",
      "25175.0\n",
      "\n",
      "With NumPy, dL = \n",
      "23200.0\n",
      "25175.0\n"
     ]
    }
   ],
   "source": [
    "# Multi-output case\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "nt = 4\n",
    "no = 3\n",
    "nx = 2 # np is numpy...\n",
    "\n",
    "print('(n_times, n_outputs, n_parameters)')\n",
    "print((nt, no, nx))\n",
    "\n",
    "times = 3 + np.arange(nt)\n",
    "print('times')\n",
    "print(times.shape)\n",
    "print(times)\n",
    "assert(times.shape == (nt, ))\n",
    "\n",
    "r = np.array([[3, 1, 2], [-1, 3, -2], [-2, 1, 7], [3, 4, 5]])\n",
    "print('residuals')\n",
    "print(r.shape)\n",
    "print(r)\n",
    "assert(r.shape == (nt, no))\n",
    "\n",
    "dm = np.zeros((nt, no, nx))\n",
    "z = 0\n",
    "for i in range(nt):\n",
    "    for j in range(no):\n",
    "        for k in range(nx):\n",
    "            dm[i, j, k] = z\n",
    "            z += 3\n",
    "del(i, j, k, z)\n",
    "print('derivatives')\n",
    "print(dm.shape)\n",
    "print(dm)\n",
    "assert(dm.shape == (nt, no, nx))\n",
    "\n",
    "# Pre-calculate fixed offset and multipliers\n",
    "sigma = np.array([0.1, 0.2, 0.3])\n",
    "offset = -0.5 * nt * np.log(2 * np.pi) - nt * np.log(sigma)\n",
    "multip = -1 / (2.0 * sigma**2)\n",
    "isigma2 = sigma**-2\n",
    "\n",
    "# Calculate log-likelihood\n",
    "L = np.sum(offset + multip * np.sum(r**2, axis=0))\n",
    "\n",
    "print()\n",
    "print('L:')\n",
    "print(L)\n",
    "\n",
    "# Derivative should map\n",
    "#\n",
    "# (nt, no, nx) --> (nt)\n",
    "#\n",
    "#  r: (nt, no)\n",
    "# dm: (nt, no, nx)\n",
    "#\n",
    "# Multication not allowed, need to do (r.T * dm.T).T\n",
    "#\n",
    "# r.T * dm.T --> (no, nt) * (nx, no, nt) --> (nx, no, nt)\n",
    "# (r.T * dm.T).T --> (nt, no, nx)\n",
    "#\n",
    "# Then sum over axis=0 to drop out t, leaving (no, nx)\n",
    "# Then multiply by sigma of shape (no, )\n",
    "# Then second sum to drop out no\n",
    "#\n",
    "\n",
    "# Derivative, with sums:\n",
    "dL = np.zeros(nx)\n",
    "for k in range(nx):\n",
    "    isum = 0\n",
    "    ijsum = np.zeros(no)\n",
    "    for i in range(nt):\n",
    "        ijsum += r[i] * dm[i, :, k]\n",
    "    for j in range(no):\n",
    "        isum += ijsum[j] / sigma[j]**2\n",
    "    dL[k] = isum\n",
    "\n",
    "print()\n",
    "print('With sums, dL = ')\n",
    "print(dL[0])\n",
    "print(dL[1])\n",
    "\n",
    "# Reshape, in case of single-output\n",
    "r = r.reshape(nt, no)\n",
    "dm = dm.reshape(nt, no, nx)\n",
    "\n",
    "# Calculate derivative\n",
    "dL = np.sum((isigma2 * np.sum((r.T * dm.T).T, axis=0).T).T, axis=0)\n",
    "\n",
    "print()\n",
    "print('With NumPy, dL = ')\n",
    "print(dL[0])\n",
    "print(dL[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(n_times, n_parameters)\n",
      "(3, 2)\n",
      "times\n",
      "(3,)\n",
      "[3 4 5]\n",
      "residuals\n",
      "(3,)\n",
      "[ 5.   5.5  6. ]\n",
      "derivatives\n",
      "(3, 2)\n",
      "[[  3.   5.]\n",
      " [  7.   9.]\n",
      " [ 11.  13.]]\n",
      "\n",
      "L:\n",
      "-4558.34906032\n",
      "\n",
      "With sums, dL = \n",
      "11950.0\n",
      "15250.0\n",
      "\n",
      "With NumPy, dL = \n",
      "11950.0\n",
      "15250.0\n"
     ]
    }
   ],
   "source": [
    "# Single-output case\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "nt = 3\n",
    "no = 1\n",
    "nx = 2 # np is numpy...\n",
    "\n",
    "print('(n_times, n_parameters)')\n",
    "print((nt, nx))\n",
    "\n",
    "times = 3 + np.arange(nt)\n",
    "print('times')\n",
    "print(times.shape)\n",
    "print(times)\n",
    "\n",
    "r = 5 + 0.5 * np.arange(nt)\n",
    "print('residuals')\n",
    "print(r.shape)\n",
    "print(r)\n",
    "\n",
    "dm = np.zeros((nt, nx))\n",
    "z = 3\n",
    "for i in range(nt):\n",
    "    for k in range(nx):\n",
    "        dm[i, k] = z\n",
    "        z += 2\n",
    "del(i, k, z)\n",
    "print('derivatives')\n",
    "print(dm.shape)\n",
    "print(dm)\n",
    "\n",
    "# Pre-calculate fixed offset and multipliers\n",
    "sigma = 0.1\n",
    "offset = -0.5 * nt * np.log(2 * np.pi) - nt * np.log(sigma)\n",
    "multip = -1 / (2.0 * sigma**2)\n",
    "isigma2 = sigma**-2\n",
    "\n",
    "# Calculate log-likelihood\n",
    "L = np.sum(offset + multip * np.sum(r**2, axis=0))\n",
    "print()\n",
    "print('L:')\n",
    "print(L)\n",
    "\n",
    "# Derivative, with sums:\n",
    "dL = np.zeros(nx)\n",
    "for k in range(nx):\n",
    "    isum = 0\n",
    "    for i in range(nt):\n",
    "        isum += r[i] * dm[i, k]\n",
    "    dL[k] = isum / sigma**2\n",
    "        \n",
    "print()\n",
    "print('With sums, dL = ')\n",
    "print(dL[0])\n",
    "print(dL[1])\n",
    "\n",
    "# Calculate with NumPy\n",
    "\n",
    "# Reshape, in case of single-output\n",
    "r = r.reshape(nt, no)\n",
    "dm = dm.reshape(nt, no, nx)\n",
    "\n",
    "dL = np.sum((isigma2 * np.sum((r.T * dm.T).T, axis=0).T).T, axis=0)\n",
    "\n",
    "print()\n",
    "print('With NumPy, dL = ')\n",
    "print(dL[0])\n",
    "print(dL[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.3 Unknown noise: treating $\\sigma$ as a parameter\n",
    "\n",
    "If we don't have accurate information about $\\sigma$, we can try treating it as a parameter an infering it along with the others.\n",
    "In code, we'll do this by simply extending $x$ to include $\\sigma$.\n",
    "In the descriptions below, we'll leave $\\sigma$ separate for clarity.\n",
    "\n",
    "The Gaussian independent log-likelihood becomes:\n",
    "\n",
    "$$ L(x, \\sigma) = - \\frac{n_t n_o}{2} \\log(2\\pi) \n",
    "          - n_t \\sum_{j=1}^{n_o} \\log(\\sigma_j)\n",
    "          - \\sum_{j=1}^{n_o} \\left[ \\frac{1}{2\\sigma_j^2} \\sum_{i=1}^{n_t} \\left(m_{ij}(x) - v_{ij} \\right)^2 \\right]\n",
    "$$\n",
    "\n",
    "where $x$ is the original parameter vector $x_1, x_2, ..., x_{n_p}$ and $\\sigma$ is a vector of the standard deviations in each output $\\sigma_1, \\sigma_2, ..., \\sigma_{n_o}$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.4.1 Derivatives\n",
    "\n",
    "For the derivatives, we now find $n_p$ equations of the previous form:\n",
    "\n",
    "$$ \\left. \\frac{\\partial L}{\\partial x_k} \\right|_{x, \\sigma} =\n",
    "    \\sum_{j=1}^{n_o} \\left[ \\frac{1}{\\sigma_j^2} \\sum_{i=1}^{n_t}\n",
    "    \\left(v_{ij} - m_{ij}(x) \\right) \\left. \\frac{\\partial m_{ij}}{\\partial x_k} \\right|_{x, \\sigma} \\right]\n",
    "$$\n",
    "\n",
    "where $k \\in 1, 2, ..., n_p$.\n",
    "In addition, we get $n_o$ equations:\n",
    "\n",
    "$$ \\left. \\frac{\\partial L}{\\partial x_m} \\right|_{x, \\sigma} =\n",
    "    - n_t / \\sigma_m\n",
    "    + \\sigma_m^{-3} \\sum_{i=1}^{n_t} \\left(m_{im}(x) - v_{im} \\right)^2\n",
    "$$\n",
    "\n",
    "where $m \\in 1, 2, ..., n_o$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Priors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.1 Gaussian prior, 1-dimensional\n",
    "\n",
    "We define a 1-dimensional Gaussian log-prior as \n",
    "\n",
    "$$\n",
    "\\frac{1}{\\sqrt{2 \\pi \\sigma^2}} - \\frac{(x - \\mu)^2}{2 \\sigma^2}\n",
    "$$\n",
    "\n",
    "where $\\mu$ is the prior's mean and $\\sigma$ is its standard deviation.\n",
    "\n",
    "The single parameter is denoted $x$.\n",
    "\n",
    "\n",
    "#### Derivatives\n",
    "\n",
    "The derivative of this prior with respect to $x$ is given by\n",
    "\n",
    "$$\n",
    "\\frac{\\mu - x}{\\sigma^2}\n",
    "$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. LogPosterior\n",
    "\n",
    "A logposterior is the sum of a logprior and a loglikelihood\n",
    "\n",
    "### Derivative\n",
    "\n",
    "The derivative of a logposterior is simply the sum of the derivatives of its logprior and of its loglikelihood."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6. Error measures\n",
    "\n",
    "All error measures implement some measure $E(x)$ that is minimised at $E(x_\\text{true})$, but is otherwise completely unrestrained in its values, smoothness, etc.\n",
    "\n",
    "Any LogLikelihood can be made into an error measure by reversing it's sign."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6.1 Sum of squares\n",
    "\n",
    "One of the simplest error measures is\n",
    "\n",
    "$$ E(x) = \\sum_{i=1}^{n_t} \\left( m_i(x) - v_i \\right)^2 $$\n",
    "\n",
    "For multiple outputs, all weighted equally, this becomes\n",
    "\n",
    "$$ E(x) = \\sum_{j=1}^{n_o} \\sum_{i=1}^{n_t} \\left( m_{ij}(x) - v_{ij} \\right)^2 $$\n",
    "\n",
    "For its derivative, we find\n",
    "\n",
    "$$ \\left. \\frac{\\partial E}{\\partial x_k} \\right|_x = \n",
    "    2 \\sum_{j=1}^{n_o} \\sum_{i=1}^{n_t} \n",
    "    \\left( m_{ij}(x) - v_{ij} \\right)\n",
    "    \\left. \\frac{\\partial m_{ij}}{\\partial x_k} \\right|_x\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(n_times, n_outputs, n_parameters)\n",
      "(3, 2, 2)\n",
      "times\n",
      "(3,)\n",
      "[0 1 2]\n",
      "residuals\n",
      "(3, 2)\n",
      "[[ 0.  0.]\n",
      " [ 1.  1.]\n",
      " [ 2.  2.]]\n",
      "derivatives\n",
      "(3, 2, 2)\n",
      "[[[  0.   1.]\n",
      "  [  2.   3.]]\n",
      "\n",
      " [[  4.   5.]\n",
      "  [  6.   7.]]\n",
      "\n",
      " [[  8.   9.]\n",
      "  [ 10.  11.]]]\n",
      "Expected output\n",
      "92.0\n",
      "104.0\n",
      "NumPy output\n",
      "[  92.  104.]\n"
     ]
    }
   ],
   "source": [
    "# Multi-output case\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "nt = 3\n",
    "no = 2\n",
    "nx = 2 # np is numpy...\n",
    "\n",
    "print('(n_times, n_outputs, n_parameters)')\n",
    "print((nt, no, nx))\n",
    "\n",
    "times = np.arange(nt)\n",
    "print('times')\n",
    "print(times.shape)\n",
    "print(times)\n",
    "\n",
    "r = (np.arange(nt).reshape((nt, 1))).dot(np.ones((1, no)) )\n",
    "print('residuals')\n",
    "print(r.shape)\n",
    "print(r)\n",
    "\n",
    "dm = np.zeros((nt, no, nx))\n",
    "z = 0\n",
    "for i in range(nt):\n",
    "    for j in range(no):\n",
    "        for k in range(nx):\n",
    "            dm[i, j, k] = z\n",
    "            z += 1\n",
    "del(i, j, k, z)\n",
    "print('derivatives')\n",
    "print(dm.shape)\n",
    "print(dm)\n",
    "\n",
    "# With sums\n",
    "print('Expected output')\n",
    "dE1 = 2 * np.sum([np.sum([np.sum(r[i, j] * dm[i, j, 0]) for i in range(nt)]) for j in range(no)])\n",
    "dE2 = 2 * np.sum([np.sum([np.sum(r[i, j] * dm[i, j, 1]) for i in range(nt)]) for j in range(no)])\n",
    "print(dE1)\n",
    "print(dE2)\n",
    "\n",
    "# NumPy way\n",
    "print('NumPy output')\n",
    "dE = 2 * np.sum((r.T * dm.T), axis=(1,2))\n",
    "print(dE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(n_times, n_parameters)\n",
      "(3, 2)\n",
      "times\n",
      "(3,)\n",
      "[0 1 2]\n",
      "residuals\n",
      "(3,)\n",
      "[0 1 2]\n",
      "derivatives\n",
      "(3, 2)\n",
      "[[ 0.  1.]\n",
      " [ 2.  3.]\n",
      " [ 4.  5.]]\n",
      "Expected output\n",
      "20\n",
      "26\n",
      "With sums\n",
      "20.0\n",
      "26.0\n",
      "NumPy output\n",
      "[ 20.  26.]\n"
     ]
    }
   ],
   "source": [
    "# Single-output case\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "nt = 3\n",
    "no = 1\n",
    "nx = 2 # np is numpy...\n",
    "\n",
    "print('(n_times, n_parameters)')\n",
    "print((nt, nx))\n",
    "\n",
    "times = np.arange(nt)\n",
    "print('times')\n",
    "print(times.shape)\n",
    "print(times)\n",
    "\n",
    "r = np.arange(nt)\n",
    "print('residuals')\n",
    "print(r.shape)\n",
    "print(r)\n",
    "\n",
    "dm = np.zeros((nt, nx))\n",
    "z = 0\n",
    "for i in range(nt):\n",
    "    for k in range(nx):\n",
    "        dm[i, k] = z\n",
    "        z += 1\n",
    "del(i, k, z)\n",
    "print('derivatives')\n",
    "print(dm.shape)\n",
    "print(dm)\n",
    "\n",
    "# Manual\n",
    "print('Expected output')\n",
    "dE1 = 2 * (0*0 + 1*2 + 2*4)\n",
    "dE2 = 2 * (0*1 + 1*3 + 2*5)\n",
    "print(dE1)\n",
    "print(dE2)\n",
    "\n",
    "# With sums\n",
    "print('With sums')\n",
    "dE1 = 2 * np.sum([np.sum(r[i] * dm[i, 0]) for i in range(nt)])\n",
    "dE2 = 2 * np.sum([np.sum(r[i] * dm[i, 1]) for i in range(nt)])\n",
    "print(dE1)\n",
    "print(dE2)\n",
    "\n",
    "# NumPy form\n",
    "print('NumPy output')\n",
    "dm = dm.reshape((nt, no, nx))\n",
    "dE = 2 * np.sum((r.T * dm.T), axis=(1,2))\n",
    "print(dE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6.2 Weighted sum of squares\n",
    "\n",
    "We can scale the error in each output by introducing weighing factors $w_j$ for $j \\in 1, 2, ..., n_o$. Then the error becomes\n",
    "\n",
    "$$ E(x) = \\sum_{j=1}^{n_o} w_j \\left( \\sum_{i=1}^{n_t} \\left( m_{ij}(x) - v_{ij} \\right)^2 \\right)$$\n",
    "\n",
    "For its derivative, we have\n",
    "\n",
    "$$ \\left. \\frac{\\partial E}{\\partial x_k} \\right|_x = \n",
    "    2 \\sum_{j=1}^{n_o} w_j \\left( \\sum_{i=1}^{n_t} \n",
    "    \\left( m_{ij}(x) - v_{ij} \\right)\n",
    "    \\left. \\frac{\\partial m_{ij}}{\\partial x_k} \\right|_x \\right)\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(n_times, n_outputs, n_parameters)\n",
      "(3, 2, 2)\n",
      "times\n",
      "(3,)\n",
      "[0 1 2]\n",
      "residuals\n",
      "(3, 2)\n",
      "[[ 0.  0.]\n",
      " [ 1.  1.]\n",
      " [ 2.  2.]]\n",
      "weights\n",
      "(2,)\n",
      "[1 2]\n",
      "derivatives\n",
      "(3, 2, 2)\n",
      "[[[  0.   1.]\n",
      "  [  2.   3.]]\n",
      "\n",
      " [[  4.   5.]\n",
      "  [  6.   7.]]\n",
      "\n",
      " [[  8.   9.]\n",
      "  [ 10.  11.]]]\n",
      "Expected output\n",
      "144.0\n",
      "162.0\n",
      "NumPy output\n",
      "[ 144.  162.]\n"
     ]
    }
   ],
   "source": [
    "# Multi-output case\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "nt = 3\n",
    "no = 2\n",
    "nx = 2 # np is numpy...\n",
    "\n",
    "print('(n_times, n_outputs, n_parameters)')\n",
    "print((nt, no, nx))\n",
    "\n",
    "times = np.arange(nt)\n",
    "print('times')\n",
    "print(times.shape)\n",
    "print(times)\n",
    "\n",
    "r = (np.arange(nt).reshape((nt, 1))).dot(np.ones((1, no)) )\n",
    "print('residuals')\n",
    "print(r.shape)\n",
    "print(r)\n",
    "\n",
    "w = np.arange(no) + 1\n",
    "print('weights')\n",
    "print(w.shape)\n",
    "print(w)\n",
    "\n",
    "dm = np.zeros((nt, no, nx))\n",
    "z = 0\n",
    "for i in range(nt):\n",
    "    for j in range(no):\n",
    "        for k in range(nx):\n",
    "            dm[i, j, k] = z\n",
    "            z += 1\n",
    "del(i, j, k, z)\n",
    "print('derivatives')\n",
    "print(dm.shape)\n",
    "print(dm)\n",
    "\n",
    "# With sums\n",
    "print('Expected output')\n",
    "dE1 = 2 * np.sum([np.sum([np.sum(r[i, j] * dm[i, j, 0]) for i in range(nt)]) * w[j] for j in range(no)])\n",
    "dE2 = 2 * np.sum([np.sum([np.sum(r[i, j] * dm[i, j, 1]) for i in range(nt)]) * w[j] for j in range(no)])\n",
    "print(dE1)\n",
    "print(dE2)\n",
    "\n",
    "# NumPy way\n",
    "print('NumPy output')\n",
    "dE = 2 * np.sum(np.sum((r.T * dm.T), axis=2) * w, axis=1)\n",
    "print(dE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6.3 Mean squared error\n",
    "\n",
    "In some cases, it can be desirable to make the error invariant to the length of the sample $n_t$.\n",
    "\n",
    "$$ E(x) = \\frac{1}{n_t} \\sum_{i=1}^{n_t} \\left( m_i(x) - v_i \\right)^2 $$\n",
    "\n",
    "For multiple outputs, with weighting $\\{w_j\\}$ for each output $j$, this becomes\n",
    "\n",
    "$$ E(x) = \\frac{1}{n_t n_o} \\sum_{j=1}^{n_o} w_j \\left( \\sum_{i=1}^{n_t} \\left( m_{ij}(x) - v_{ij} \\right)^2 \\right)$$\n",
    "\n",
    "For its derivative, we find\n",
    "\n",
    "$$ \\left. \\frac{\\partial E}{\\partial x_k} \\right|_x = \n",
    "\\frac{2}{n_t n_o} \\sum_{j=1}^{n_o} w_j \\left( \\sum_{i=1}^{n_t} \n",
    "    \\left( m_{ij}(x) - v_{ij} \\right)\n",
    "    \\left. \\frac{\\partial m_{ij}}{\\partial x_k} \\right|_x \\right)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(n_times, n_outputs, n_parameters)\n",
      "(3, 2, 2)\n",
      "times\n",
      "(3,)\n",
      "[0 1 2]\n",
      "residuals\n",
      "(3, 2)\n",
      "[[ 0.  0.]\n",
      " [ 1.  1.]\n",
      " [ 2.  2.]]\n",
      "weights\n",
      "(2,)\n",
      "[1 2]\n",
      "derivatives\n",
      "(3, 2, 2)\n",
      "[[[  0.   1.]\n",
      "  [  2.   3.]]\n",
      "\n",
      " [[  4.   5.]\n",
      "  [  6.   7.]]\n",
      "\n",
      " [[  8.   9.]\n",
      "  [ 10.  11.]]]\n",
      "Expected output\n",
      "24.0\n",
      "27.0\n",
      "NumPy output\n",
      "[ 24.  27.]\n"
     ]
    }
   ],
   "source": [
    "# Multi-output case\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "nt = 3\n",
    "no = 2\n",
    "nx = 2 # np is numpy...\n",
    "\n",
    "print('(n_times, n_outputs, n_parameters)')\n",
    "print((nt, no, nx))\n",
    "\n",
    "times = np.arange(nt)\n",
    "print('times')\n",
    "print(times.shape)\n",
    "print(times)\n",
    "\n",
    "r = (np.arange(nt).reshape((nt, 1))).dot(np.ones((1, no)) )\n",
    "print('residuals')\n",
    "print(r.shape)\n",
    "print(r)\n",
    "\n",
    "w = np.arange(no) + 1\n",
    "print('weights')\n",
    "print(w.shape)\n",
    "print(w)\n",
    "\n",
    "dm = np.zeros((nt, no, nx))\n",
    "z = 0\n",
    "for i in range(nt):\n",
    "    for j in range(no):\n",
    "        for k in range(nx):\n",
    "            dm[i, j, k] = z\n",
    "            z += 1\n",
    "del(i, j, k, z)\n",
    "print('derivatives')\n",
    "print(dm.shape)\n",
    "print(dm)\n",
    "\n",
    "# With sums\n",
    "print('Expected output')\n",
    "dE1 = 2. / nt / no * np.sum([np.sum([np.sum(r[i, j] * dm[i, j, 0]) for i in range(nt)]) * w[j] for j in range(no)])\n",
    "dE2 = 2. / nt / no * np.sum([np.sum([np.sum(r[i, j] * dm[i, j, 1]) for i in range(nt)]) * w[j] for j in range(no)])\n",
    "print(dE1)\n",
    "print(dE2)\n",
    "\n",
    "# NumPy way\n",
    "print('NumPy output')\n",
    "dE = 2. / nt / no * np.sum(np.sum((r.T * dm.T), axis=2) * w, axis=1)\n",
    "print(dE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(n_times, n_parameters)\n",
      "(3, 2)\n",
      "times\n",
      "(3,)\n",
      "[0 1 2]\n",
      "residuals\n",
      "(3,)\n",
      "[0 1 2]\n",
      "derivatives\n",
      "(3, 2)\n",
      "[[ 0.  1.]\n",
      " [ 2.  3.]\n",
      " [ 4.  5.]]\n",
      "Expected output\n",
      "20\n",
      "26\n",
      "With sums\n",
      "6.66666666667\n",
      "8.66666666667\n",
      "NumPy output\n",
      "[ 6.66666667  8.66666667]\n"
     ]
    }
   ],
   "source": [
    "# Single-output case\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "nt = 3\n",
    "no = 1\n",
    "nx = 2 # np is numpy...\n",
    "\n",
    "print('(n_times, n_parameters)')\n",
    "print((nt, nx))\n",
    "\n",
    "times = np.arange(nt)\n",
    "print('times')\n",
    "print(times.shape)\n",
    "print(times)\n",
    "\n",
    "r = np.arange(nt)\n",
    "print('residuals')\n",
    "print(r.shape)\n",
    "print(r)\n",
    "\n",
    "dm = np.zeros((nt, nx))\n",
    "z = 0\n",
    "for i in range(nt):\n",
    "    for k in range(nx):\n",
    "        dm[i, k] = z\n",
    "        z += 1\n",
    "del(i, k, z)\n",
    "print('derivatives')\n",
    "print(dm.shape)\n",
    "print(dm)\n",
    "\n",
    "# Manual\n",
    "print('Expected output')\n",
    "dE1 = 2 * (0*0 + 1*2 + 2*4)\n",
    "dE2 = 2 * (0*1 + 1*3 + 2*5)\n",
    "print(dE1)\n",
    "print(dE2)\n",
    "\n",
    "# With sums\n",
    "print('With sums')\n",
    "dE1 = 2 / nt / no * np.sum([np.sum(r[i] * dm[i, 0]) for i in range(nt)])\n",
    "dE2 = 2 / nt / no * np.sum([np.sum(r[i] * dm[i, 1]) for i in range(nt)])\n",
    "print(dE1)\n",
    "print(dE2)\n",
    "\n",
    "# NumPy form\n",
    "print('NumPy output')\n",
    "dm = dm.reshape((nt, no, nx))\n",
    "dE = 2 / nt / no * np.sum((r.T * dm.T), axis=(1,2))\n",
    "print(dE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6.4 Root-mean squared error\n",
    "\n",
    "A common error measure is\n",
    "\n",
    "$$ E(x) = \\sqrt{ \\frac{1}{n_t} \\sum_{i=1}^{n_t} \\left( m_i(x) - v_i \\right)^2 } $$\n",
    "\n",
    "For multiple outputs, all weighted equally, this becomes\n",
    "\n",
    "$$ E(x) = \\sqrt{  \\frac{1}{n_t n_o} \\sum_{j=1}^{n_o} \\sum_{i=1}^{n_t} \\left( m_{ij}(x) - v_{ij} \\right)^2 } $$\n",
    "\n",
    "For its derivative, we find\n",
    "\n",
    "$$ \\left. \\frac{\\partial E}{\\partial x_k} \\right|_x = \n",
    "    \\frac{\n",
    "        \\frac{2}{n_t n_o} \\sum_{j=1}^{n_o} \\sum_{i=1}^{n_t} \n",
    "        \\left( m_{ij}(x) - v_{ij} \\right)\n",
    "        \\left. \\frac{\\partial m_{ij}}{\\partial x_k} \\right|_x\n",
    "    }{\n",
    "        \\sqrt{ \\frac{1}{n_t n_o} \\sum_{j=1}^{n_o} \\sum_{i=1}^{n_t} \\left( m_{ij}(x) - v_{ij} \\right)^2 }\n",
    "    }\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  }
 ],
 "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.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}