{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# eda03_00a clear all variables and import vatious modules\n",
    "%reset -f\n",
    "import os\n",
    "from datetime import date\n",
    "from math import exp, pi, sin, sqrt, floor, ceil\n",
    "import numpy as np\n",
    "import scipy.linalg as la\n",
    "from matplotlib import pyplot as plt\n",
    "import matplotlib.cm as cm\n",
    "from matplotlib.colors import ListedColormap\n",
    "\n",
    "# eda_draw function makes a \"pictorial matrix equation\"\n",
    "# arguments are vectors, matrices and strings\n",
    "# which are plotted in the order that the appear\n",
    "# except that strings starting with 'title ' are plotted\n",
    "# under the subseqeunt matrix or vector\n",
    "# always returns a status of 1\n",
    "def eda_draw(*argv):\n",
    "    bw = np.zeros((256,4));\n",
    "    v = 0.9*(256 - np.linspace( 0, 255, 256 ))/255;\n",
    "    bw[:,0] = v;\n",
    "    bw[:,1] = v;\n",
    "    bw[:,2] = v;\n",
    "    bw[:,3] = np.ones(256);\n",
    "    bwcmap = ListedColormap(bw);\n",
    "    # size of plot\n",
    "    W = 16;\n",
    "    H = 4;\n",
    "    fig1 = plt.figure(1);\n",
    "    # figsize width and height in inches\n",
    "    fig1.set_size_inches(W,H);\n",
    "    ax1 = plt.subplot(1,1,1);\n",
    "    plt.axis([0, W, -H/2, H/2]);\n",
    "    plt.axis('off');\n",
    "    LM = W/6;    # matrix width and heoght\n",
    "    LV = W/40;   # vector width\n",
    "    FS = 0.12;    # character width\n",
    "    TO = 0.4;    # title vertical offset\n",
    "    SP = 0.2;    # space between objects\n",
    "    LS = 0.2;    # leading space\n",
    "    p = LS; # starting x-position\n",
    "    istitle=0; # flags presence of a title\n",
    "    for a in argv:\n",
    "        if isinstance(a,np.ndarray):\n",
    "            sh = np.shape(a);\n",
    "            if len(sh) == 1:  # conversion to nx1 array\n",
    "                n = sh[0];\n",
    "                m = 1;\n",
    "                ap = a;\n",
    "                a = np.zeros((n,1));\n",
    "                a[:,0] = ap;\n",
    "            else:\n",
    "                n = sh[0];\n",
    "                m = sh[1];\n",
    "            if m==1:\n",
    "                pold=p;\n",
    "                left=p;\n",
    "                right=p+LV;\n",
    "                bottom=-LM/2;\n",
    "                top=LM/2;\n",
    "                plt.imshow( a, cmap=bwcmap, vmin=np.min(a), vmax=np.max(a), extent=(left,right,bottom,top) );\n",
    "                p = p+LV;\n",
    "                pm = (p+pold)/2;\n",
    "                if istitle:\n",
    "                    plt.text(pm,-(LM/2)-TO,titlestr,horizontalalignment='center');\n",
    "                    istitle=0;\n",
    "                p = p+SP;\n",
    "            else:\n",
    "                pold=p;\n",
    "                left=p;\n",
    "                right=p+LM;\n",
    "                bottom=-LM/2;\n",
    "                top=LM/2;\n",
    "                plt.imshow( a, cmap=bwcmap, vmin=np.min(a), vmax=np.max(a), extent=(left,right,bottom,top) );\n",
    "                p = p+LM;\n",
    "                pm = (p+pold)/2;\n",
    "                if istitle:\n",
    "                    plt.text(pm,-(LM/2)-TO,titlestr,horizontalalignment='center');\n",
    "                    istitle=0;\n",
    "                p = p+SP;\n",
    "        elif isinstance(a,str):\n",
    "            ns = len(a);\n",
    "            istitle=0;\n",
    "            if( ns>=6 ):\n",
    "                if 'title ' in a[0:6]:\n",
    "                    istitle=1;\n",
    "                    titlestr=a[6:];\n",
    "            if( istitle != 1):\n",
    "                plt.text(p,0,a);\n",
    "                p = p + ns*FS + SP;\n",
    "    plt.show();\n",
    "    return 1;\n",
    "\n",
    "# function to make a numpy N-by-1 column vector\n",
    "# c=eda_cvec1(v) from entity v that is array-like,\n",
    "# including a number, a list of numbers, a tuple\n",
    "# of numbers, an N-by-0 np array and a N-by1 np array.\n",
    "# It promotes integers to floats, and integers\n",
    "# and floats to complex, by context. The function\n",
    "# also insures that, if v is an np array, that\n",
    "# c is a copy, as contrasted to a view, of it.\n",
    "# This version takes just one input argment.\n",
    "# whereas c=eda_cvec(v1,v2,...) concatenates\n",
    "# many argiments.\n",
    "def eda_cvec1(v):\n",
    "    if isinstance(v, int):\n",
    "        w = np.zeros((1,1),dtype=int);\n",
    "        w[0,0] = v;\n",
    "        return w;\n",
    "    elif isinstance(v, float):\n",
    "        w = np.zeros((1,1),dtype=float);\n",
    "        w[0,0] = v;\n",
    "        return w;\n",
    "    elif isinstance(v, complex):\n",
    "        w = np.zeros((1,1),dtype=complex);\n",
    "        w[0,0] = v;\n",
    "        return w;\n",
    "    elif isinstance(v, np.ndarray):\n",
    "        s = np.shape(v);\n",
    "        if len(s) == 1:\n",
    "            return np.copy(np.reshape(v,(s[0],1)));\n",
    "        else:\n",
    "            [r,c]=s;\n",
    "            if( c==1 ):\n",
    "                return(np.copy(v));\n",
    "            elif(r==1):\n",
    "                return(np.copy(v.T));\n",
    "            else:\n",
    "                raise TypeError(\"eda_cvec: %d by %d ndarray not allowed\" % (r, c));\n",
    "    elif isinstance(v, list):\n",
    "        r = len(v);\n",
    "        t = int;\n",
    "        for vi in v:\n",
    "            if isinstance(vi,int):\n",
    "                pass;\n",
    "            elif isinstance(vi,float):\n",
    "                t=float;\n",
    "            elif isinstance(vi,complex):\n",
    "                t=complex;\n",
    "                break;\n",
    "            else:\n",
    "                raise TypeError(\"eda_cvec: list contains unsupported type %s\" % type(vi));\n",
    "        w = np.zeros((r,1),dtype=t);\n",
    "        w[:,0] = v;\n",
    "        return w;\n",
    "    elif isinstance(v, tuple):\n",
    "        r = len(v);\n",
    "        t = int;\n",
    "        for vi in v:\n",
    "            if isinstance(vi,int):\n",
    "                pass;\n",
    "            elif isinstance(vi,float):\n",
    "                t=float;\n",
    "            elif isinstance(vi,complex):\n",
    "                t=complex;\n",
    "                break;\n",
    "            else:\n",
    "                raise TypeError(\"eda_cvec: tuple contains unsupported type %s\" % type(vi));\n",
    "        w = np.zeros((r,1),dtype=t);\n",
    "        w[:,0] = v;\n",
    "        return w;\n",
    "    else:\n",
    "        raise TypeError(\"eda_cvec: %s not supported\" % type(v));\n",
    "        \n",
    "# function to make a numpy N-by-1 column vector\n",
    "# c=eda_cvec(v1, v2, ...) from a list of several\n",
    "# array-like entities v1, v2, including a number\n",
    "# a list of numbers, a tuple of numbers, an N-by-0 np array\n",
    "# and a N-by-1 np array. The function\n",
    "# also insures that, if v is an np array, that\n",
    "# c is a copy, as contrasted to a view, of it\n",
    "# It promotes integers to floats, and integers\n",
    "# and floats to complex, by context.\n",
    "# This version concatenates many argments,\n",
    "# whereas c=eda_cvec1(v1) takes just one argument.\n",
    "# I recommend always using eda_cvec(v1, v2, ...)\n",
    "def eda_cvec(*argv):\n",
    "    t = int;\n",
    "    Nt = 0;\n",
    "    for a in argv:\n",
    "        v = eda_cvec1(a);\n",
    "        N,M = np.shape(v);\n",
    "        Nt = Nt + N;\n",
    "        if (t==int) and isinstance(v[0,0],float):\n",
    "            t=float;\n",
    "        elif isinstance(v[0,0],complex):\n",
    "            t=complex;\n",
    "    w = np.zeros((Nt,1),dtype=t);\n",
    "    Nt = 0;\n",
    "    for a in argv:\n",
    "        v = eda_cvec1(a);\n",
    "        N,M = np.shape(v);\n",
    "        w[Nt:Nt+N,0] = v[0:N,0];\n",
    "        Nt = Nt + N;\n",
    "    return w;\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "#eda03_00b: definition of the eda_draw() function\n",
    "\n",
    "# eda_draw function makes a \"pictorial matrix equation\"\n",
    "# arguments are vectors, matrices and strings\n",
    "# which are plotted in the order that the appear\n",
    "# except that strings starting with 'title ' are plotted\n",
    "# under the subseqeunt matrix or vector\n",
    "# always returns a status of 1\n",
    "\n",
    "def eda_draw(*argv):\n",
    "    bw = np.zeros((256,4));\n",
    "    v = 0.9*(256 - np.linspace( 0, 255, 256 ))/255;\n",
    "    bw[:,0] = v;\n",
    "    bw[:,1] = v;\n",
    "    bw[:,2] = v;\n",
    "    bw[:,3] = np.ones(256);\n",
    "    bwcmap = ListedColormap(bw);\n",
    "    # size of plot\n",
    "    W = 16;\n",
    "    H = 4;\n",
    "    fig1 = plt.figure(1);\n",
    "    # figsize width and height in inches\n",
    "    fig1.set_size_inches(W,H);\n",
    "    ax1 = plt.subplot(1,1,1);\n",
    "    plt.axis([0, W, -H/2, H/2]);\n",
    "    plt.axis('off');\n",
    "    LM = W/6;    # matrix width and heoght\n",
    "    LV = W/40;   # vector width\n",
    "    FS = 0.12;    # character width\n",
    "    TO = 0.4;    # title vertical offset\n",
    "    SP = 0.2;    # space between objects\n",
    "    LS = 0.2;    # leading space\n",
    "    p = LS; # starting x-position\n",
    "    istitle=0; # flags presence of a title\n",
    "    for a in argv:\n",
    "        if isinstance(a,np.ndarray):\n",
    "            sh = np.shape(a);\n",
    "            if len(sh) == 1:  # conversion to nx1 array\n",
    "                n = sh[0];\n",
    "                m = 1;\n",
    "                ap = a;\n",
    "                a = np.zeros((n,1));\n",
    "                a[:,0] = ap;\n",
    "            else:\n",
    "                n = sh[0];\n",
    "                m = sh[1];\n",
    "            if m==1:\n",
    "                pold=p;\n",
    "                left=p;\n",
    "                right=p+LV;\n",
    "                bottom=-LM/2;\n",
    "                top=LM/2;\n",
    "                plt.imshow( a, cmap=bwcmap, vmin=np.min(a), vmax=np.max(a), extent=(left,right,bottom,top) );\n",
    "                p = p+LV;\n",
    "                pm = (p+pold)/2;\n",
    "                if istitle:\n",
    "                    plt.text(pm,-(LM/2)-TO,titlestr,horizontalalignment='center');\n",
    "                    istitle=0;\n",
    "                p = p+SP;\n",
    "            else:\n",
    "                pold=p;\n",
    "                left=p;\n",
    "                right=p+LM;\n",
    "                bottom=-LM/2;\n",
    "                top=LM/2;\n",
    "                plt.imshow( a, cmap=bwcmap, vmin=np.min(a), vmax=np.max(a), extent=(left,right,bottom,top) );\n",
    "                p = p+LM;\n",
    "                pm = (p+pold)/2;\n",
    "                if istitle:\n",
    "                    plt.text(pm,-(LM/2)-TO,titlestr,horizontalalignment='center');\n",
    "                    istitle=0;\n",
    "                p = p+SP;\n",
    "        elif isinstance(a,str):\n",
    "            ns = len(a);\n",
    "            istitle=0;\n",
    "            if( ns>=6 ):\n",
    "                if 'title ' in a[0:6]:\n",
    "                    istitle=1;\n",
    "                    titlestr=a[6:];\n",
    "            if( istitle != 1):\n",
    "                plt.text(p,0,a);\n",
    "                p = p + ns*FS + SP;\n",
    "    plt.show();\n",
    "    return 1;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3QAAADnCAYAAABWihTLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAJeElEQVR4nO3dPYtkaRnH4ftZW1lxTTTQBgMNTGQRDAyFRTDSTIQ1UT+AaLKfwJdETDZXE1/jRTTxJRAEWVBYNdLAQNmt7e3Z2p7qLqurqx+DWWQFxTmnzpyaf9d1RdMM93OfbvrA/HhmplvvvQAAAMjzxKEfAAAAgHEEHQAAQChBBwAAEErQAQAAhBJ0AAAAoQQdAABAKEEHAAAQStABAACEEnQAAAChBB0AAEAoQQcAABBK0AEAAIQSdAAAAKEEHQAAQChBBwAAEErQAQAAhBJ0AAAAoQQdAABAKEEHAAAQStABAACEEnQAAAChBB0AAEAoQQcAABBK0AEAAIQSdAAAAKEEHQAAQChBBwAAEErQAQAAhBJ0AAAAoQQdAABAKEEHAAAQStABAACEEnQAAAChBB0AAEAoQQcAABBK0AEAAIQSdAAAAKEEHQAAQChBBwAAEErQAQAAhBJ0AAAAoQQdAABAqJOpDnr55Zf7VGe91enpaXsU5/4vj+rz4HgM+Z69K+8NAACH4YYOAAAglKADAAAIJegAAABCCToAAIBQgg4AACCUoAMAAAgl6AAAAEIJOgAAgFCCDgAAIJSgAwAACCXoAAAAQgk6AACAUCdTHbRYLKY66j+cnp4+knOndHt7O3hmt9uN2tV7n2XXXJ/TXJ/P2LkxX4ch37PH/N4AALA/N3QAAAChBB0AAEAoQQcAABBK0AEAAIQSdAAAAKEEHQAAQChBBwAAEErQAQAAhBJ0AAAAoQQdAABAKEEHAAAQStABAACEOpnqoMViMdVRcXa73eCZm5ubUbu22+0su+7anrG7xswMcczvDQAA+3NDBwAAEErQAQAAhBJ0AAAAoQQdAABAKEEHAAAQStABAACEEnQAAAChBB0AAEAoQQcAABBK0AEAAIQSdAAAAKEEHQAAQChBBwAAEOpkqoNeffXVqY6K03sfPLPdbkft2mw2g2eur69n2fM4z4ydG/O1G+KY3xsAAPbnhg4AACCUoAMAAAgl6AAAAEIJOgAAgFCCDgAAIJSgAwAACCXoAAAAQgk6AACAUIIOAAAglKADAAAIJegAAABCCToAAIBQJ1Md9Morr0x1VJzdbjd45ubmZtSu6+vrwTPr9XrwzNXV1Sx75pqpmu9zGuKY3xsAAPbnhg4AACCUoAMAAAgl6AAAAEIJOgAAgFCCDgAAIJSgAwAACCXoAAAAQgk6AACAUIIOAAAglKADAAAIJegAAABCCToAAIBQJ1MdtFgspjoqzu3t7eCZ7XY7atdmsxk8c3V1NXhmtVo9tnvGzIydu7y8HLXrYR3zewMAwP7c0AEAAISa7IbucdFa+01Vvfu//NZzvfdfzP08AAAAj8qdC7re+ycO/QwAAABzuHNB54YOAAA4Fncu6NzQAQAAx8J/igIAABBK0AEAAIQSdAAAAKEEHQAAQChBBwAAEErQAQAAhBJ0AAAAoSb7OXSLxWKqo+LsdrvBMzc3N6N2bTabwTPr9XrwzNXV1eCZN954Y/DM/fv3B89cXFwMnhk7N3bXwzrm9wYAgP25oQMAAAgl6AAAAEIJOgAAgFCCDgAAIJSgAwAACCXoAAAAQgk6AACAUIIOAAAglKADAAAIJegAAABCCToAAIBQgg4AACDUyVQHLRaLqY6K03sfPLPdbkft2mw2g2fW6/XgmdVqNXjm/v37g2def/31wTPL5XLwzNi5sbse1jG/NwAA7M8NHQAAQChBBwAAEErQAQAAhBJ0AAAAoQQdAABAKEEHAAAQStABAACEEnQAAAChBB0AAEAoQQcAABBK0AEAAIQSdAAAAKFOpjposVhMdVSc3W43eObm5mbUrs1mM3hmvV4PnlmtVoNnLi4uBs8sl8vBM/fu3Rs8U1V1fn4+266HdczvDQAA+3NDBwAAEErQAQAAhBJ0AAAAoQQdAABAKEEHAAAQStABAACEEnQAAAChBB0AAEAoQQcAABBK0AEAAIQSdAAAAKEEHQAAQChBBwAAEOpkqoMWi8VUR8XZ7XaDZ7bb7ahdm81m8MzV1dXgmdVqNXjm4uJi8MxyuRw8c35+Pnimquq1114bPHN2djZq18M65vcGAID9uaEDAAAIJegAAABCCToAAIBQgg4AACCUoAMAAAgl6AAAAEIJOgAAgFCCDgAAIJSgAwAACCXoAAAAQgk6AACAUIIOAAAg1MlUBy0Wi6mOinN7ezt4Zrvdjtp1fX09eGa9Xg+euby8HDxzcXExeGa5XA6euXfv3uCZqqqzs7PBM4/6+/qY3xsAAPbnhg4AACCUoAMAAAgl6AAAAEIJOgAAgFCCDgAAIJSgAwAACCXoAAAAQgk6AACAUIIOAAAglKADAAAIJegAAABCCToAAIBQrfd+6GcAAABgBDd0AAAAoQQdAABAKEEHAAAQStABAACEEnQAAAChBB0AAEAoQQcAABBK0AEAAIQSdAAAAKEEHQAAQChBBwAAEErQAQAAhBJ0AAAAoQQd/1drrbfWvv+Wj09aa2ettZ8e8rkAAODYCToexmVVPd1ae+ebH3+qqv5xwOcBAABqxqBrrX29tfbVt3z8zdbaV+baz95+XlWffvPXn6+qHx/wWY5Ga+3jrbWXWmtPttbe1Vr7c2vt6UM/FwAAj4c5b+i+W1VfrKpqrT1RVc9W1Q9n3M9+flJVz7bWnqyqj1bV7w78PEeh9/5iVb1QVd+oqm9V1Q9673867FMBAPC4OJlrUe/9b62189bax6rqfVX1h977+Vz72U/v/aXW2gfrwe3czw77NEfna1X1YlX9s6rcagMA8G+zBd2bvlNVX6qq91fV92bezf5eqKpvV9UzVfXewz7KUXlPVT1VVW+vqifrwb9pBACAar33+Za19o6q+mM9+IPph3vvu9mWM1prbdV7f6q19oGq+mzv/fnW2jNV9Vzv/TMHfrw7r7X2Qj34K68fqqrT3vuXD/xIAAA8Jma9oeu9X7fWfl1VSzGXp/f+96p6/tDPcUxaa1+oqpve+49aa2+rqt+21j7Ze//VoZ8NAIDDm/uG7omq+n1Vfa73/pfZFgMAANxBc/7Ygo9U1V+r6pdiDgAAYH+z3tABAAAwnTl/Dh0AAAATEnQAAAChBB0AAEAoQQcAABBK0AEAAIQSdAAAAKEEHQAAQChBBwAAEErQAQAAhBJ0AAAAoQQdAABAKEEHAAAQStABAACEEnQAAAChBB0AAEAoQQcAABBK0AEAAIQSdAAAAKH+BQmw7q9H1TVOAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda03_00c: example of calling eda_draw()\n",
    "\n",
    "# make the matrix equation y = M x\n",
    "N=11;\n",
    "x = eda_cvec( np.linspace(0,N-1,N) );\n",
    "M = np.matmul(x,x.T);\n",
    "y = np.matmul(M,x);\n",
    "\n",
    "# pictorially plot it\n",
    "eda_draw('title y',y, '=', 'title M', M, 'title x', x);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mode [1.5]\n"
     ]
    }
   ],
   "source": [
    "# eda03_01  mode of p(d)\n",
    "\n",
    "# define an exemplary p(d)\n",
    "N=15;\n",
    "Dd = 1.0;\n",
    "d = eda_cvec( Dd*np.linspace(0, N-1, N) + 0.5*Dd );\n",
    "p = eda_cvec( [5.00, 15.00, 10.00, 10.00, 10.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00] );\n",
    "\n",
    "# normalise so integral is unity\n",
    "p = p / (Dd*np.sum(p));\n",
    "\n",
    "# mode is d for which p is maximum\n",
    "# argmax returns the index k at which p is maximum\n",
    "k = np.argmax(p);\n",
    "# the mode is the d with index k\n",
    "dmode = d[k];\n",
    "print('mode', dmode );\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "median [4.5]\n"
     ]
    }
   ],
   "source": [
    "# eda03_02 median of p(d)\n",
    "\n",
    "# define an exemplary p(d)\n",
    "N=15;\n",
    "Dd = 1.0;\n",
    "d = eda_cvec( Dd*np.linspace(0, N-1, N) + 0.5*Dd );\n",
    "p = eda_cvec( [5.00, 15.00, 10.00, 10.00, 10.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00] );\n",
    "\n",
    "# normalise so integral is unity\n",
    "p = p / (Dd*np.sum(p));\n",
    "\n",
    "# P is thecumative probability\n",
    "# Approximate integral as summation\n",
    "P = np.cumsum(p);\n",
    "# median is 50% area point of P\n",
    "# search for first occurence of P>=0.5\n",
    "for i in range(N):\n",
    "    if P[i] >= 0.5:\n",
    "        dmedian = d[i];\n",
    "        break;\n",
    "\n",
    "print('median', dmedian );"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mean 6.30\n"
     ]
    }
   ],
   "source": [
    "# eda03_03 mean of p(d)\n",
    "\n",
    "# define an exemplary p(d)\n",
    "N=15;\n",
    "Dd = 1.0;\n",
    "d = eda_cvec( Dd*np.linspace(0, N-1, N) + 0.5*Dd );\n",
    "p = eda_cvec( [5.00, 15.00, 10.00, 10.00, 10.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00] );\n",
    "\n",
    "# normalise so integral is unity\n",
    "p = p / (Dd*np.sum(p));\n",
    "\n",
    "# mean = (integral) d p(d) dd\n",
    "# approximate integral as a summation\n",
    "dmean = Dd*np.sum(np.multiply(p,d));\n",
    "print('mean %.2f' % dmean );\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "sigma 4.3197\n"
     ]
    }
   ],
   "source": [
    "# eda03_04 variance of p(d)\n",
    "\n",
    "# define an exemplary p(d)\n",
    "N=15;\n",
    "Dd = 1.0;\n",
    "d = eda_cvec( Dd*np.linspace(0, N-1, N) + 0.5*Dd );\n",
    "p = eda_cvec( [5.00, 15.00, 10.00, 10.00, 10.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00] );\n",
    "\n",
    "# normalise so integral is unity\n",
    "p = p / (Dd*np.sum(p));\n",
    "\n",
    "# mean = (integral) (d-dbar)^2 p(d) dd\n",
    "# approximate integral as sum\n",
    "dbar = Dd*np.sum(np.multiply(p,d));\n",
    "\n",
    "# variance = (integral) (d-dbar)^2 p(d) dd = (integral) q p(d) dd with q=(d-dbar)^2\n",
    "q = np.power(d-dbar,2);\n",
    "# approximate integral as sum\n",
    "sigma2 = Dd*np.sum(np.multiply(q,p));\n",
    "sigma = sqrt(sigma2);\n",
    "print('sigma %.4f' % sigma );\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 4780 characters]\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda03_05 visualation of process of computing variance\n",
    "\n",
    "# set up d axis\n",
    "Nd = 101;\n",
    "d = eda_cvec(np.linspace(-30.0,30.0, Nd));\n",
    "Dd = d[2,0]-d[1,0];\n",
    "\n",
    "q = np.power(d,2);  # q is d-squared\n",
    "\n",
    "# probability at fixed sigma\n",
    "db = 0;      # mean\n",
    "sd = 5;      # sqrt of variance\n",
    "sd2 = sd**2; # variance\n",
    "p = (1/(sqrt(2*pi)*sd)) * np.exp( -0.5*np.power(d-db,2)/sd2  ); # Normal pdf\n",
    "\n",
    "# variance calcuation (just to check)\n",
    "# variance = (integral) (d-dbar)^2 p(d) dd = (integral) q p(d) dd with q=(d-dbar)^2\n",
    "# approximate integral with sum\n",
    "qp = np.multiply(q,p);\n",
    "sigma2 =Dd * np.sum(qp);\n",
    "sigma = sqrt(sigma2);\n",
    "\n",
    "# find nearest element of array to d=dbar-sigma, d=dbar and d=dbar+sigma\n",
    "v = np.abs(d-(db-sd));\n",
    "ilo = np.argmin(v);\n",
    "v = np.abs(d-db);\n",
    "imi = np.argmin(v);\n",
    "v = np.abs(d-(db+sd));\n",
    "ihi = np.argmin(v);\n",
    "\n",
    "# make a marker vector with non-zero elements corresponding to mean, mean+/-sigma\n",
    "mk = np.zeros((Nd,1));\n",
    "mk[ilo,0] = 1;\n",
    "mk[imi,0] = 1;\n",
    "mk[ihi,0] = 1;\n",
    "\n",
    "eda_draw(mk,'title p', p, 'title d**2', q, 'title q*p', qp, mk);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 5020 characters]\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda03_06 visualation of effect of changing mean on p.d.f.\n",
    "\n",
    "# list of means\n",
    "Ndbar=5;\n",
    "dbar = eda_cvec( np.linspace(10.0,30.0, Ndbar) );\n",
    "\n",
    "# set up d-axis\n",
    "Nd = 101;\n",
    "d = eda_cvec( np.linspace(0,40, Nd) );\n",
    "\n",
    "# probability at fixed sigma\n",
    "sd = 5;       # sqrt of variacne\n",
    "sd2 = sd**2;  # variance\n",
    "titles = ();  # titles for plot is built incrementally; starts as an empty tuple\n",
    "p = np.zeros((Nd,Ndbar));  # array, each column of which is a pdf\n",
    "for i in range(Ndbar): # loop over exemplary means\n",
    "    db = dbar[i,0]; # mean\n",
    "    g = (1/(sqrt(2*pi)*sd)) * np.exp( -0.5*np.power(d-db,2)/sd2  );  # Normal pdf\n",
    "    p[:,i] = g.ravel();\n",
    "    t = \"title %.0f\" % db;  # title of plot\n",
    "    titles = titles + (t,); # append title\n",
    "\n",
    "eda_draw(titles[0], p[:,0], titles[1], p[:,1], titles[2], p[:,2], titles[3], p[:,3], titles[4], p[:,4]);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 4992 characters]\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda03_07 visualation of effect of changing variance on p.d.f.\n",
    "\n",
    "# list of sqrt-variances\n",
    "Nsigma=5;\n",
    "sigma = eda_cvec( [2.5, 5, 20, 20, 40] );\n",
    "\n",
    "Nd = 101;\n",
    "d = eda_cvec( np.linspace(0.0,40.0, Nd) );\n",
    "\n",
    "# probability at fixed sigma\n",
    "p = np.zeros((Nd,Nsigma)); # array, each column of which is a pdf\n",
    "db = 20;     # mean\n",
    "titles = (); # build list of titles incrementally, starting with an empty tuple\n",
    "for i in range(Nsigma):\n",
    "    sd = sigma[i]; # sqrt of variance\n",
    "    sd2 = sd**2;   # variance\n",
    "    g = (1/(sqrt(2*pi)*sd)) * np.exp( -0.5*np.power(d-db,2)/sd2  );\n",
    "    p[:,i]=g.ravel();\n",
    "    t = \"title %.0f\" % sd\n",
    "    titles = titles + (t,);\n",
    "\n",
    "eda_draw(titles[0], p[:,0], titles[1], p[:,1], titles[2], p[:,2], titles[3], p[:,3], titles[4], p[:,4]);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3QAAADnCAYAAABWihTLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAIOElEQVR4nO3dPciddxnH8d+VF0wiZFCXaAhBjIJRQYSKpEIdXAUHQRHRriIOTi7iCy7dBcFBFFwEF1HBoWADigjWoRUqDkVwCNhNxWKr/Tv0VMNDSZt6Xvo7z+czndzc59zXk+3Ldd/nzForAAAA9Dlz6AEAAAB4bQQdAABAKUEHAABQStABAACUEnQAAAClBB0AAEApQQcAAFBK0AEAAJQSdAAAAKUEHQAAQClBBwAAUErQAQAAlBJ0AAAApQQdAABAKUEHAABQStABAACUEnQAAAClBB0AAEApQQcAAFBK0AEAAJQSdAAAAKUEHQAAQClBBwAAUErQAQAAlBJ0AAAApQQdAABAKUEHAABQStABAACUEnQAAAClBB0AAEApQQcAAFBK0AEAAJQSdAAAAKUEHQAAQClBBwAAUErQAQAAlBJ0AAAApQQdAABAKUEHAABQStABAACUEnQAAAClBB0AAECpc9v6oDt37qxtfdbdrly5Mq/23JnZyQxrrVc9AwAAwL5sLeheD27evHnoEQAAAPbmqILuxo0bhx4BAABgb44q6K5evXroEQAAAPbmqILu8uXLhx4BAABgb44q6C5evHjoEQAAAPbmqILu0qVLhx4BAABgb44q6C5cuHDoEQAAAPbmqILu/Pnzhx4BAABgbwQdAABAKUEHAABQ6qiC7ty5o/pzAAAA7umoCujs2bOHHgEAAGBvBB0AAEApQQcAAFBK0AEAAJQ6qqCbmUOPAAAAsDdHFXQ2dAAAwGlyVEF35syZQ48AAACwN4IOAACg1FEFnWfoAACA00TQAQAAlDqqoHPLJQAAcJoIOgAAgFJHFXRuuQQAAE4TQQcAAFBK0AEAAJQSdAAAAKV8iwgAAECpo9rQ+ZZLAADgNDmqoHPLJQAAcJoIOgAAgFLuUQQAACh1VBs6z9ABAACniQICAAAodVQbOs/QAQAAp4mgAwAAKCXoAAAASnmGDgAAoJSgAwAAKOWWSwAAgFI2dAAAAKVs6AAAAEoJOgAAgFJuuQQAACgl6AAAAEq55RIAAKCUoAMAACjllksAAIBSNnQAAAClbOgAAABK2dABAACUsqEDAAAoJegAAABKueUSAACglA0dAABAKRs6AACAUjZ0AAAApQQdAABAKUEHAABQStABAACUEnQAAAClfMslAABAKRs6AACAUjZ0AAAApWzoAAAASgk6AACAUoIOAACglGfoAAAAStnQAQAAlBJ0AAAApdxyCQAAUMqGDgAAoJSgAwAAKCXoAAAASnmGDgAAoJQNHQAAQClBBwAAUErQAQAAlPIMHQAAQCkbOgAAgFKCDgAAoJSgAwAAKOUZOgAAgFI2dAAAAKUEHQAAQClBBwAAUErQAQAAlPKlKAAAAKVs6AAAAEoJOgAAgFKCDgAAoJSgAwAAKOVLUQAAAErZ0AEAAJQSdAAAAKUEHQAAQClBBwAAUErQAQAAlBJ0AAAApQQdAABAKUEHAABQStABAACUEnQAAAClzh16gG2amUOPAAAAsDc2dAAAAKUEHQAAQClBBwAAUErQAQAAlBJ0AAAApQQdAABAKT9bAAAAUMqGDgAAoJSgAwAAKCXoAAAASgk6AACAUoIOAACglKADAAAo5WcLAAAAStnQAQAAlBJ0AAAApQQdAABAKUEHAABQStABAACUEnQAAAClBB0AAEApv0MHAABQyoYOAACglKADAAAoJegAAABKCToAAIBSgg4AAKCUoAMAACgl6AAAAEoJOgAAgFJ+WBwAAKCUDR0AAEApQQcAAFBK0AEAAJQSdAAAAKUEHQAAQClBBwAAUErQAQAAlBJ0AAAApQQdAABAqVlrHXoGAAAAXgMbOgAAgFKCDgAAoJSgAwAAKCXoAAAASgk6AACAUoIOAACglKADAAAoJegAAABKCToAAIBSgg4AAKCUoAMAACgl6AAAAEoJOgAAgFJ7DbqZ+dHMvP1ljn9uZr61ef2FmXl4X9e+x/nvnZnvbXsOAACAbdlb0M3MzSRn11pPv8Kp303yxQNd+7/WWk8muToz17Y5CwAAwLZsPehm5vrM/GFmvj8zT2w2Y5eSfDrJj+867+GZ+ePM3E5y66Xja61/JPnTzDyww2v/fWYemZnHZ+bRmXlgZh6bmadn5mN3feRPknzy/v8XAAAAdm9XG7p3JfnOWut9Sf6a5PN5MdoeT5KZuZLk65tjH03y7hPv/22SD+/i2htvTPLYWusDSf6W5JubOT6e5BtbmgMAAGCndhV0f15r/Wrz+gdJHkxyJckzm2MfzItB9cxa67kkPzzx/r8keeuOrp0kzyX5+eb1k0lur7We37y+vqU5AAAAdmpXQbde5t/PJrlwj3PudmFz/q6u/fxa66XzXkjyzyRZa72Q5NyW5gAAANipXQXdtZn50Ob1p5L8MslTSd6xOfabJA/NzJtn5nyST5x4/zuT/H5H174f/88cAAAAO7WroHsqyWdn5okkb0ry7SQ/S/JQkqy17iT5WpJfJ3k0ye9OvP/W5vjWr32fPrJ5LwAAwOvO/O/Owy194Mz1JD9da73nxPGLSX6R5NZa69/3eP/7k3xprfWZfV/7xHvekOR2kgfXWv+631kAAAB2bW+/Q7fWejbJV5O87RVOfUuSrxzo2ne7luTLYg4AAHi92vqGDgAAgP3Y24YOAACA7RJ0AAAApQQdAABAKUEHAABQStABAACUEnQAAAClBB0AAEApQQcAAFBK0AEAAJQSdAAAAKUEHQAAQClBBwAAUErQAQAAlBJ0AAAApQQdAABAKUEHAABQStABAACUEnQAAACl/gN4Bs0ocnwabgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda03_08 transformation of variables m = d^2\n",
    "\n",
    "Nd = 101;\n",
    "d = eda_cvec( np.linspace(0.01,1.0, Nd) );\n",
    "Dd = d[2,0] - d[1,0];\n",
    "\n",
    "# uniform pdf p(d)\n",
    "pd = np.ones((Nd,1))/Dd;\n",
    "\n",
    "# set up m-axis\n",
    "Nm = 101;\n",
    "m = eda_cvec( np.linspace(0.01,1, Nm) );\n",
    "Dm = m[2,0] - m[1,0];\n",
    "\n",
    "# transformed, non-uniform pdf p(m)\n",
    "pm = 0.5*np.reciprocal(np.sqrt(m));\n",
    "\n",
    "eda_draw( 'title p(d)', pd, 'title p(m)', pm);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 6616 characters]\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda03_09 joint p.f.d.\n",
    "\n",
    "# set up d1-axis\n",
    "Nd1 = 101;\n",
    "d1 = eda_cvec( np.linspace(0.0,5.0, Nd1) );\n",
    "Dd1 = d1[1,0] - d1[0,0];\n",
    "\n",
    "# set up d2-axis\n",
    "Nd2 = 101;\n",
    "d2 = eda_cvec( np.linspace(0,5,Nd2) );\n",
    "Dd2 = d2[1,0] - d2[0,0];\n",
    "\n",
    "# 2D normal p.f.d. setup\n",
    "dbar1 = 1.5;\n",
    "dbar2 = 2.5;\n",
    "sigma1 = 0.25;\n",
    "sigma2 = 0.25;\n",
    "\n",
    "# 2-d Normal p.d.f.\n",
    "norm1 = (1/(sqrt(2*pi)*sigma1));\n",
    "norm2 = (1/(sqrt(2*pi)*sigma2));\n",
    "p1 = norm1 * np.exp( - np.power(d1-dbar1,2) / (2*sigma1*sigma1) );\n",
    "p2 = norm2 * np.exp( - np.power(d2-dbar2,2) / (2*sigma2*sigma2) );\n",
    "# note use of tensor product to construct p(d1,d2)=p(d1)p(d2)\n",
    "p = np.matmul( p1, p2.T );  \n",
    "\n",
    "# integrate along d2-axis to get p2, appoximating integral as sum\n",
    "p_1 = np.zeros((Nd1,1));\n",
    "p_1[:,0] = Dd2*np.sum(p,axis=1);      \n",
    "        \n",
    "# integrate along d1-axis to get p2, appoximating integral as sum\n",
    "p_2 = np.zeros((Nd2,1));\n",
    "p_2[:,0] = Dd1*np.sum(p,axis=0);\n",
    "\n",
    "eda_draw( 'title p(d1,d2)', p, 'title p(d1)', p_1, 'title p(d2)', p_2);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "area 1.0000\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3QAAADnCAYAAABWihTLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAH6ElEQVR4nO3dT6jlZR3H8c+3jDQQotrckkEkFMqCCIlQw4ighQqJgtFCa9k6qHCRrtuFULkyqCAwpDKCFNT+UIINMhWai4gKBnLnQiEmnxb3N3Ic9Oq9MzB9Zl6v1e889/k9v+cs3zznnDtrrQAAANDnbed7AwAAAByNoAMAACgl6AAAAEoJOgAAgFKCDgAAoJSgAwAAKCXoAAAASgk6AACAUoIOAACglKADAAAoJegAAABKCToAAIBSgg4AAKCUoAMAACgl6AAAAEoJOgAAgFKCDgAAoJSgAwAAKCXoAAAASgk6AACAUoIOAACglKADAAAoJegAAABKCToAAIBSgg4AAKCUoAMAACgl6AAAAEoJOgAAgFKCDgAAoJSgAwAAKCXoAAAASgk6AACAUoIOAACglKADAAAoJegAAABKCToAAIBSgg4AAKCUoAMAACgl6AAAAEoJOgAAgFKCDgAAoJSgAwAAKHXJuVro5MmT61ytxcVpb29vzvceAACgiRM6AACAUoIOAACglKADAAAoJegAAABKCToAAIBSgg4AAKCUoAMAACgl6AAAAEoJOgAAgFKCDgAAoJSgAwAAKCXoAAAASgk6AACAUoIOAACglKADAAAoJegAAABKCToAAIBSgg4AAKCUoAMAACgl6AAAAEoJOgAAgFKCDgAAoJSgAwAAKCXoAAAASgk6AACAUoIOAACglKADAAAoJegAAABKCToAAIBSgg4AAKCUoAMAACgl6AAAAEoJOgAAgFKCDgAAoJSgAwAAKCXoAAAASgk6AACAUoIOAACglKADAAAoJegAAABKCToAAIBSgg4AAKCUoAMAACgl6AAAAEoJOgAAgFKCDgAAoJSgAwAAKCXoAAAASgk6AACAUoIOAACglKADAAAoJegAAABKCToAAIBSgg4AAKCUoAMAACgl6AAAAEoJOgAAgFKCDgAAoJSgAwAAKCXoAAAASgk6AACAUoIOAACglKADAAAoJegAAABKCToAAIBSgg4AAKCUoAMAACgl6AAAAEoJOgAAgFKCDgAAoJSgAwAAKCXoAAAASgk6AACAUoIOAACglKADAAAoJegAAABKCToAAIBSgg4AAKCUoAMAACgl6AAAAEoJOgAAgFKCDgAAoJSgAwAAKCXoAAAASgk6AACAUoIOAACglKADAAAoJegAAABKCToAAIBSgg4AAKCUoAMAACgl6AAAAEoJOgAAgFKCDgAAoJSgAwAAKCXoAAAASgk6AACAUoIOAACglKADAAAoJegAAABKCToAAIBSgg4AAKCUoAMAACgl6AAAAEoJOgAAgFKCDgAAoJSgAwAAKCXoAAAASgk6AACAUoIOAACglKADAAAoJegAAABKCToAAIBSgg4AAKCUoAMAACgl6AAAAEoJOgAAgFKCDgAAoJSgAwAAKCXoAAAASgk6AACAUoIOAACglKADAAAoJegAAABKCToAAIBSgg4AAKCUoAMAACgl6AAAAEoJOgAAgFKCDgAAoJSgAwAAKCXoAAAASgk6AACAUoIOAACglKADAAAoJegAAABKCToAAIBSgg4AAKCUoAMAACgl6AAAAEoJOgAAgFKCDgAAoJSgAwAAKCXoAAAASs1a63zvAQAAgCNwQgcAAFBK0AEAAJQSdAAAAKUEHQAAQClBBwAAUErQAQAAlBJ0AAAApQQdAABAKUEHAABQStABAACUEnQAAAClBB0AAEApQQcAAFBK0B3BzDw0M1e9zvjdM3P/dv2pmTk+M6dm5va3uO69M/PV7fpbM/PczJyYmYdn5t3b+Edm5sFz+HYAAIBSgu6QZubDSd6+1vrbm0z9R5K7k/zoiI96NMm1a62PJnk+yTeSZK31pyRXzMyxI64LAABcIATdG5iZK7cTsu9vp2QPzcy7knwxyU935n1pZp6fmSeTXH96fK3197XWiSSvvMlz7pmZv87MY0mu2bn/V2utU9vLPyS5Yue2nye586zfJAAAUE3QHeyaJA9sp2QvJvlK9qPtj0kyM3tJ7tvGPpvkQ4dZfGY+nv0w+1iS25Jc9wZTv5zklzuvn05y42GeBQAAXHgE3cH+udb63Xb9gyQ3JNlL8sI29okkT6y1Xlhr/SfJjw+5/o1JHl5rvbTWejHJz86cMDP3JDmV5Ic7w/9O8v5DPgsAALjAXHK+N/B/br3O65eTXHrAnLN9xqtm5q4kNyf5zFprd96l2z4AAICLmBO6gx2bmU9u119I8tskzyb54Db2VJKbZua9M/OOJHe8lUVn5rnt8tdJPj8zl83M5Ulu2ZnzuSRfS3LrWuulM5a4Osmfj/KGAACAC4egO9izSe6amRNJ3pPkO0l+keSmJFlrnUxyb5LfJ3ksyfHTN87MdTPzr+xH3vdm5i/b+PuSzHb/8ex/TPOZJD9J8pudZ9+f5PIkj87MMzPz3Z2/fXrbBwAAcBGb136Sj9Nm5sokj6y1rj1j/LIkjye5fq313yOse3OSq9Za3z7ivt6Z5MkkN+z8CiYAAHAR8h26Q1prvTwz30zygez/r7nD3v/IWW7hWJKvizkAAMAJHQAAQCnfoQMAACgl6AAAAEoJOgAAgFKCDgAAoJSgAwAAKCXoAAAASgk6AACAUoIOAACglKADAAAoJegAAABKCToAAIBSgg4AAKCUoAMAACgl6AAAAEoJOgAAgFKCDgAAoJSgAwAAKCXoAAAASv0Pt3m7jl/XJtMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda03_10  normalising white p.d.f.\n",
    "\n",
    "# set up d1-axis\n",
    "Nd1 = 101;\n",
    "d1 = eda_cvec( np.linspace(0.0,5.0, Nd1) );\n",
    "Dd1 = d1[1,0] - d1[0,0];\n",
    "\n",
    "# set up d2-axis\n",
    "Nd2 = 101;\n",
    "d2 = eda_cvec( np.linspace(0,5, Nd2) );\n",
    "Dd2 = d2[1,0] - d2[0,0];\n",
    "\n",
    "# calulae normalisation factor so that integral is unity\n",
    "# approximating integgral as sum\n",
    "p = np.ones((Nd1, Nd2));\n",
    "norm = Dd1*Dd2*np.sum(p);\n",
    "# normalize the pdf\n",
    "p = np.copy(p)/norm;\n",
    "\n",
    "# area of pdf (should be 1)   \n",
    "A = Dd1*Dd2*np.sum(p);\n",
    "print('area %.4f' % A);\n",
    "\n",
    "eda_draw( 'title p(d1,d2)', p);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "area   1.0000 1.0000 1.0000\n",
      "dbar1  1.5000 1.5000 1.5000\n",
      "dbar2  2.5000 2.5000 2.5000\n",
      "sigma1 0.2500 0.2500 0.2500\n",
      "sigma2 0.2500 0.2500 0.2500\n"
     ]
    }
   ],
   "source": [
    "# eda03_11, 12 and 13 constructing a joint p.f.d. and computing its mean and variance\n",
    "\n",
    "# set up d1 axis\n",
    "Nd1 = 101;\n",
    "d1 = eda_cvec( np.linspace(0.0,5.0, Nd1) );\n",
    "Dd1 = d1[1,0] - d1[0,0];\n",
    "\n",
    "# set up d2 axis\n",
    "Nd2 = 101;\n",
    "d2 = eda_cvec( np.linspace(0,5,Nd2) );\n",
    "Dd2 = d2[1,0] - d2[0,0];\n",
    "\n",
    "# 2D normal p.f.d. setup\n",
    "dbar1 = 1.5;\n",
    "dbar2 = 2.5;\n",
    "sigma1 = 0.25;\n",
    "sigma2 = 0.25;\n",
    "\n",
    "# since d1 and d2 are uncorrelated, the 2-d Normal\n",
    "# p.d.f. p(d1,d2) is product of p(d1) and p(d2)\n",
    "norm1 = (1/(sqrt(2*pi)*sigma1));   # normalization of p(d1) \n",
    "norm2 = (1/(sqrt(2*pi)*sigma2));   # normalization of p(d2) \n",
    "p1 = norm1 * np.exp( - np.power(d1-dbar1,2) / (2*sigma1*sigma1) );  # p(d1)\n",
    "p2 = norm2 * np.exp( - np.power(d2-dbar2,2) / (2*sigma2*sigma2) );  # p(d2)\n",
    "p = np.matmul( p1, p2.T );  # p(d1,d2) i= p(d1) p(d2)\n",
    "\n",
    "# integrate along d2-axis to get p2\n",
    "p_1 = eda_cvec( Dd2*np.sum(p,axis=1) );      \n",
    "        \n",
    "# integrate along d1-axis to get p2\n",
    "p_2 = eda_cvec( Dd1*np.sum(p,axis=0) );\n",
    "\n",
    "# areas (should all be 1) (approximate integrals as sums)\n",
    "A = Dd1*Dd2*np.sum(p);\n",
    "A1 = Dd1*np.sum(p_1);\n",
    "A2 = Dd2*np.sum(p_2);\n",
    "\n",
    "# vectors of ones\n",
    "w1 = np.ones((Nd1,1));\n",
    "w2 = np.ones((Nd2,1));\n",
    "\n",
    "# eda03_12: means\n",
    "\n",
    "# mean along 1 axis, computed from p\n",
    "D1 = np.matmul( d1, w2.T );\n",
    "dbar_1 = Dd1*Dd2*np.sum( np.multiply(D1, p) );\n",
    "\n",
    "# mean along 1 axis, computed from p_1\n",
    "dbar_1a = Dd1*np.sum( np.multiply(d1, p_1) );\n",
    "\n",
    "# mean along 2 axis, computed from p\n",
    "D2 = np.matmul( w1, d2.T  );\n",
    "dbar_2 = Dd1*Dd2*np.sum( np.multiply(D2, p) );\n",
    "\n",
    "# mean along 2 axis, computed from p_2\n",
    "dbar_2a = Dd2*np.sum( np.multiply(d2, p_2) );\n",
    "\n",
    "# eda03_13: variances\n",
    "\n",
    "# variance along 1 axis, computed from p(d1,d2)\n",
    "# use a tensor product to create a 2D version of the quantity (d1-d1bar)**2\n",
    "D1s = np.matmul( np.power(d1-dbar1,2), w2.T ); \n",
    "# compute the variance with the usual integral, approximated as a sum\n",
    "sigma_1_squared = Dd1*Dd2*np.sum( np.multiply(D1s, p) );\n",
    "sigma_1 = sqrt(sigma_1_squared);\n",
    "\n",
    "# variance along 1 axis, computed from p(d1) via the usual integral, approximated as a sum\n",
    "sigma_1a_squared = Dd1*np.sum( np.multiply( np.power(d1-dbar_1,2), p_1) );\n",
    "sigma_1a = sqrt(sigma_1_squared);\n",
    "\n",
    "# variance along 2 axis, computed from p(d1,d2)\n",
    "# use a tensor product to create a 2D version of the quantity (d2-d2bar)**2\n",
    "D2s = np.matmul( w1, np.power(d2-dbar2,2).T  );\n",
    "# compute the variance with the usual integral, approximated as a sum\n",
    "sigma2_squared = Dd1*Dd2*np.sum( np.multiply(D2s, p) );\n",
    "sigma_2 = sqrt(sigma2_squared);\n",
    "\n",
    "# variance along 2 axis, computed from p(d2) via the usual integral, approximated as a sum\n",
    "sigma_2a_squared = Dd2*np.sum( np.multiply( np.power(d2-dbar_2,2), p_2) );\n",
    "sigma_2a = sqrt(sigma_2a_squared);\n",
    "\n",
    "print('area   %.4f %.4f %.4f' % (A, A1, A2));\n",
    "print('dbar1  %.4f %.4f %.4f' % (dbar1, dbar_1, dbar_1a));\n",
    "print('dbar2  %.4f %.4f %.4f' % (dbar2, dbar_2, dbar_2a));\n",
    "print('sigma1 %.4f %.4f %.4f' % (sigma1, sigma_1, sigma_1a));\n",
    "print('sigma2 %.4f %.4f %.4f' % (sigma2, sigma_2, sigma_2a));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 12740 characters]\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda03_14 example of a conditional p.f.d.\n",
    "\n",
    "# set up d1-axis\n",
    "Nd1 = 101;\n",
    "d1 = eda_cvec( np.linspace(0.0,5.0, Nd1) );\n",
    "Dd1 = d1[1,0] - d1[0,0];\n",
    "\n",
    "# set up d2-axis\n",
    "Nd2 = 101;\n",
    "d2 = eda_cvec( np.linspace(0.0,5.0,Nd2) );\n",
    "Dd2 = d2[1,0] - d2[0,0];\n",
    "\n",
    "# 2D normal p.f.d. setup\n",
    "dbar1 = 1.5;\n",
    "dbar2 = 3.5;\n",
    "sigma1 = 0.75;\n",
    "sigma2 = 0.75;\n",
    "\n",
    "# since d1 and d2 are uncorrelated, the 2-d Normal\n",
    "# p.d.f. p(d1,d2) is product of p(d1) and p(d2)\n",
    "norm1 = (1/(sqrt(2*pi)*sigma1));  # p(d1) normalization\n",
    "norm2 = (1/(sqrt(2*pi)*sigma2));  # p(d2) normalization\n",
    "p1 = norm1 * np.exp( - np.power(d1-dbar1,2) / (2*sigma1*sigma1) ); # p(d1)\n",
    "p2 = norm2 * np.exp( - np.power(d2-dbar2,2) / (2*sigma2*sigma2) ); # p(d2)\n",
    "p = np.matmul( p1, p2.T );  # use tensor product to make p(d1,d2) = p(d1)p(d2)\n",
    "\n",
    "# integrate p(d1,d2) along d2-axis to get p2\n",
    "p_1 = eda_cvec( Dd2*np.sum(p,axis=1) );      \n",
    "        \n",
    "# integrate p(d1,d2) along d1-axis to get p2\n",
    "p_2 = eda_cvec( Dd1*np.sum(p,axis=0) );\n",
    "\n",
    "# vectors of ones\n",
    "w1 = np.ones((Nd1,1));\n",
    "w2 = np.ones((Nd2,1));\n",
    "\n",
    "# conditional p.d.f.: p(d1|d2) = p(d1,d2)/p2\n",
    "p1g2 = np.divide( p, np.matmul(w1,p_2.T));\n",
    "\n",
    "# conditional p.d.f.: p(d2|d1) = p(d1,d2)/p2\n",
    "p2g1 = np.divide( p, np.matmul(p_1,w2.T));\n",
    "\n",
    "eda_draw( 'title p(d1,d2)', p, 'title p(d1|d2)', p1g2, 'title p(d2|d1)', p2g1);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 5124 characters]\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda03_15  Constructing Normal p.d.f. without for loop\n",
    "\n",
    "# set up d1-axis\n",
    "Nd1 = 101;\n",
    "d1 = eda_cvec( np.linspace(0.0,5.0, Nd1) );\n",
    "Dd1 = d1[1,0] - d1[0,0];\n",
    "\n",
    "# set up d2-axis\n",
    "Nd2 = 101;\n",
    "d2 = eda_cvec( np.linspace(0,5,Nd2) );\n",
    "Dd2 = d2[1,0] - d2[0,0];\n",
    "\n",
    "# vectors of ones, later used in tensor products\n",
    "w1 = np.ones((Nd1,1));\n",
    "w2 = np.ones((Nd2,1));\n",
    "\n",
    "# 2D normal p.f.d. setup\n",
    "dbar1 = 1.5;\n",
    "dbar2 = 2.5;\n",
    "sigma1 = 0.25;\n",
    "sigma2 = 0.25;\n",
    "covar = 0.05;  # positive covariance\n",
    "covpos = covar;\n",
    "C = np.array( [ [sigma1**2, covar], [covar, sigma2**2]] ); # covariance matrix\n",
    "Ci = la.inv(C); # its inverse\n",
    "sqrtdet = sqrt(la.det(C)); # square root of its determinant\n",
    "norm=(1/(2*pi*sqrtdet));  # normalization\n",
    "\n",
    "# Normal p.d.f.\n",
    "dd1 = d1-dbar1;  # d1 minus its mean\n",
    "dd2 = d2-dbar2;  # d2 minus its mean\n",
    "dd11 = np.power(dd1,2); # (d1 minus its mean) squared\n",
    "dd22 = np.power(dd2,2); # (d2 minus its mean) squared\n",
    "# 2D Normal pdf p(d1,d2) using tensor products\n",
    "# my sense its that this expression is way to complicate to easy check\n",
    "# and that an impementation using for loops, though computationally inefficient,\n",
    "# would be clearer (see 03_16, below)\n",
    "p=norm*np.exp(-0.5*Ci[0,0]*np.matmul(dd11,w2.T)-0.5*Ci[1,1]*np.matmul(w1,dd22.T)-Ci[0,1]*np.matmul(dd1,dd2.T));\n",
    "\n",
    "eda_draw( 'title cov +', p); "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cov +: 0.0500 0.0500\n",
      "cov -: -0.0500 -0.0500\n",
      "cov 0: 0.0000 -0.0000\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 9480 characters]\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda03_16, 17  Constructing Normal p.d.f. with for loop and computing its covariance\n",
    "\n",
    "# set up d1-axis\n",
    "Nd1 = 101;\n",
    "d1 = eda_cvec( np.linspace(0,5, Nd1) );\n",
    "Dd1 = d1[1,0] - d1[0,0];\n",
    "\n",
    "# set up d2-axis\n",
    "Nd2 = 101;\n",
    "d2 = eda_cvec( np.linspace(0,5,Nd2) );\n",
    "Dd2 = d2[1,0] - d2[0,0];\n",
    "\n",
    "# Part 1: p(d1,d2) with positive covariance\n",
    "\n",
    "# 2D normal p.f.d. setup\n",
    "dbar1 = 1.5;\n",
    "dbar2 = 2.5;\n",
    "sigma1 = 0.25;\n",
    "sigma2 = 0.25;\n",
    "covar = 0.05;\n",
    "covpos = covar;\n",
    "Cd = [ [sigma1**2, covar], [covar, sigma2**2]];\n",
    "Cdinv = la.inv(Cd);\n",
    "sqrtdet = sqrt(la.det(Cd));\n",
    "\n",
    "# eda03_16: constructing Normal p.d.f. with covariance\n",
    "# 2-d Normal p.d.f.\n",
    "p = np.zeros((Nd1,Nd2));\n",
    "dd = np.zeros((2,1));\n",
    "for i in range(Nd1):\n",
    "    for j in range(Nd2):\n",
    "        dd[:,0] = [ d1[i,0] - dbar1, d2[j,0] - dbar2 ]; # vector of data minus their means\n",
    "        E = np.matmul( np.matmul(dd.T,Cdinv), dd ); # ddT inv(C) dd\n",
    "        p[i,j] = (1/(2*pi*sqrtdet)) * exp( -0.5*E ); # p(d1,d2)\n",
    "ppos = np.copy(p);\n",
    "\n",
    "# part 2: p(d1,d2) with negative covariance\n",
    "\n",
    "# 2D normal p.f.d. setup\n",
    "dbar1 = 1.5;\n",
    "dbar2 = 2.5;\n",
    "sigma1 = 0.25;\n",
    "sigma2 = 0.25;\n",
    "covar = -0.05;\n",
    "covneg = covar;\n",
    "Cd = [ [sigma1**2, covar], [covar, sigma2**2]]; # covariance matrix\n",
    "Cdinv = la.inv(Cd); # its inverse\n",
    "sqrtdet = sqrt(la.det(Cd)); # square root of it determinant\n",
    "\n",
    "p = np.zeros((Nd1,Nd2));\n",
    "dd = np.zeros((2,1));\n",
    "for i in range(Nd1):\n",
    "    for j in range(Nd2):\n",
    "        dd[:,0] = [ d1[i,0] - dbar1, d2[j,0] - dbar2 ]; # vector of data minus their means\n",
    "        E = np.matmul( np.matmul(dd.T,Cdinv), dd ); # ddT inv(C) dd\n",
    "        p[i,j] = (1/(2*pi*sqrtdet)) * exp( -0.5*E ); # p(d1,d2)\n",
    "pneg = np.copy(p);\n",
    "        \n",
    "# Part 3: p(d1,d2) with zero covariance\n",
    "\n",
    "# 2D normal p.f.d. setup\n",
    "dbar1 = 1.5;\n",
    "dbar2 = 2.5;\n",
    "sigma1 = 0.25;\n",
    "sigma2 = 0.25;\n",
    "covar = 0.0;\n",
    "covzero = covar;\n",
    "Cd = [ [sigma1**2, covar], [covar, sigma2**2]]; # covariance matrix\n",
    "Cdinv = la.inv(Cd); # its inverse\n",
    "sqrtdet = sqrt(la.det(Cd)); # square root of it determinant\n",
    "\n",
    "# 2-d Normal p.d.f.\n",
    "p = np.zeros((Nd1,Nd2));\n",
    "dd = np.zeros((2,1));\n",
    "for i in range(Nd1):\n",
    "    for j in range(Nd2):\n",
    "        dd[:,0] = [ d1[i,0] - dbar1, d2[j,0] - dbar2 ]; # vector of data minus their means\n",
    "        E = np.matmul( np.matmul(dd.T,Cdinv), dd ); # ddT inv(C) dd\n",
    "        p[i,j] = (1/(2*pi*sqrtdet)) * exp( -0.5*E ); # p(d1,d2)\n",
    "pzero=np.copy(p);\n",
    "\n",
    "# eda03_17\n",
    "# covariance, computed from p\n",
    "D12 = np.matmul( (d1-dbar1), (d2-dbar2).T ); # use tensor product to make (d1-dbar1)(d2-dbar2)\n",
    "covpos_a = Dd1*Dd2*np.sum( np.multiply(D12, ppos) ); # covariance intergal, + covariance case\n",
    "covneg_a = Dd1*Dd2*np.sum( np.multiply(D12, pneg) ); # covariance intergal, - covariance case\n",
    "covzero_a = Dd1*Dd2*np.sum( np.multiply(D12, pzero) ); # covariance intergal, 0 covariance case\n",
    "\n",
    "print( 'cov +: %.4f %.4f' % (covpos, covpos_a) );\n",
    "print( 'cov -: %.4f %.4f' % (covneg, covneg_a) );\n",
    "print( 'cov 0: %.4f %.4f' % (covzero, covzero_a) );\n",
    "\n",
    "eda_draw( 'title cov +', ppos, 'title cov -', pneg, 'title cov 0', pzero);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "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.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
