{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# edapy_03_00 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( N==0 ):\n",
    "            continue; # skip vector of zero length\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": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mode 1.5\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 13724 characters]\n",
      "text/plain": [
       "<Figure size 720x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edapy_03_01: mode calculated from 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,0];\n",
    "print('mode', dmode );\n",
    "\n",
    "# plot the p.d.f with the mode as a dotted line\n",
    "fig1 = plt.figure(1,figsize=(10,4));\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "plt.axis([np.min(d), np.max(d), 0.0, np.max(p)]);\n",
    "plt.plot(d, p, 'k-');\n",
    "plt.plot([dmode,dmode],[0,np.max(p)/10],'k:');\n",
    "plt.xlabel('d');\n",
    "plt.ylabel('p(d)');\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "median 4.5\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 13724 characters]\n",
      "text/plain": [
       "<Figure size 720x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edapy_03_02: median calculated from 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 the cumulative probability\n",
    "# approximate integral as summation\n",
    "P = eda_cvec( Dd*np.cumsum(p) );\n",
    "\n",
    "# median is d at the 50% area point of P\n",
    "# search for first occurrence of P>=0.5\n",
    "for i in range(N):\n",
    "    if P[i,0] >= 0.5:\n",
    "        dmedian = d[i,0];\n",
    "        break;\n",
    "\n",
    "print('median', dmedian );\n",
    "\n",
    "# plot the p.d.f with the median as a dotted line\n",
    "fig1 = plt.figure(1,figsize=(10,4));\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "plt.axis([np.min(d), np.max(d), 0.0, np.max(p)]);\n",
    "plt.plot(d, p, 'k-');\n",
    "plt.plot([dmedian,dmedian],[0,np.max(p)/10],'k:');\n",
    "plt.xlabel('d');\n",
    "plt.ylabel('p(d)');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mean 6.30\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 13720 characters]\n",
      "text/plain": [
       "<Figure size 720x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edaps_03_03: mean calculated from 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",
    "# perform the integal for the mean\n",
    "dmean = Dd*np.sum(np.multiply(p,d));\n",
    "print('mean %.2f' % dmean );\n",
    "\n",
    "# plot the p.d.f with the mean as a dotted line\n",
    "fig1 = plt.figure(1,figsize=(10,4));\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "plt.axis([np.min(d), np.max(d), 0.0, np.max(p)]);\n",
    "plt.plot(d, p, 'k-'); # plot data d1\n",
    "plt.plot([dmean,dmean],[0,np.max(p)/10],'k:');\n",
    "plt.xlabel('d');\n",
    "plt.ylabel('p(d)');\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "sigma 4.3197\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 13848 characters]\n",
      "text/plain": [
       "<Figure size 720x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edapy_03_04: variance calculated from 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",
    "dmean = Dd*np.sum(np.multiply(p,d));\n",
    "\n",
    "# perform the integral for the variance\n",
    "q = np.power(d-dmean,2);\n",
    "sigma2 = Dd*np.sum(np.multiply(q,p));\n",
    "sigma = sqrt(sigma2);\n",
    "print('sigma %.4f' % sigma );\n",
    "\n",
    "# plot the p.d.f with the mean as a solid line and +/- 1 sigma as dotted line\n",
    "fig1 = plt.figure(1,figsize=(10,4));\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "plt.axis([np.min(d), np.max(d), 0.0, np.max(p)]);\n",
    "plt.plot(d, p, 'k-'); # plot data d1\n",
    "plt.plot([dmean,dmean],[0,np.max(p)/10],'k-'); # plot zero line\n",
    "plt.plot([dmean-sigma,dmean-sigma],[0,np.max(p)/10],'k:'); # plot zero line\n",
    "plt.plot([dmean+sigma,dmean+sigma],[0,np.max(p)/10],'k:'); # plot zero line\n",
    "plt.xlabel('d');\n",
    "plt.ylabel('p(d)');\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 4444 characters]\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edapy_03_05: visualation of process of computing variance\n",
    "\n",
    "# The plot had three vectors, p(d), q(d) and p(d)q(d)\n",
    "# depicted as grey-shaded coun vectors, side-by-side\n",
    "# with mean and +/- one sigma levels marked\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 q', q, 'title qp', qp, mk);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "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": [
    "# edapy_03_06: Normal p.d.f.'s with different means\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": 80,
   "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": [
    "# edapy_03_07: Normal p.d.f.'s with different variances\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": 81,
   "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": [
    "# edapy_03_08: transformation of variables m = d^2\n",
    "\n",
    "# compare original p.d.f. is uniform; tat is, p(d)=1 on 0<d<1\n",
    "# the transformation of variables is m = d^2\n",
    "# with transformed p.d.f. is p(m)=0.5*m^(-0.5) on 0<m<1\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": 88,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 6040 characters]\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edapy_03_09: The joint p.d.f. p(d1,d2) and the two univariate\n",
    "# p.d.f.'s p(d1) and p(d2) formed from it\n",
    "\n",
    "# set up the d1 and d2 axes\n",
    "L = 40;\n",
    "Dd1 = 1.0;\n",
    "Ddq = 1.0;\n",
    "d1 = eda_cvec( Dd1*np.linspace(0,L-1,L) );\n",
    "d2 = eda_cvec( Dd2*np.linspace(0,L-1,L) );\n",
    "\n",
    "# 2D normal p.f.d. setup\n",
    "dbar1 = 10.0;\n",
    "dbar2 = 20.0;\n",
    "sigma1 = 5.0;\n",
    "sigma2 = 5.0;\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 = eda_cvec( Dd*np.sum(p,axis=1) );\n",
    "        \n",
    "# integrate along d1-axis to get p2, appoximating integral as sum\n",
    "p_2 = eda_cvec( Dd*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": 15,
   "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": [
    "# edapy_03_10  The uniform p.d.f. p(d1,d2)=constant\n",
    "\n",
    "# set up the d1 and d2 axes\n",
    "Nd1 = 40;\n",
    "Nd2 = 40;\n",
    "Dd1 = 1.0;\n",
    "Dd2 = 1.0;\n",
    "d1 = eda_cvec( Dd1*np.linspace(0,Nd1-1,Nd1) );\n",
    "d2 = eda_cvec( Dd2*np.linspace(0,Nd2-1,Nd2) );\n",
    "\n",
    "# make a uniform p.d.f. and normalize it\n",
    "P = np.ones((Nd1, Nd2));\n",
    "# calculate the area under the p.d.f.\n",
    "norm = Dd1*Dd2*np.sum(P);\n",
    "# normalize the p.d.f. by dividing by the area\n",
    "P = P/norm;\n",
    "\n",
    "# area of pdf (should be 1)   \n",
    "A = Dd*Dd*np.sum(P);\n",
    "print('area %.4f' % A);\n",
    "\n",
    "eda_draw( 'title p(d1,d2)', P);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 4600 characters]\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "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": [
    "# edapy_03_11: construct a joint Normal p.f.d.\n",
    "# and compute 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",
    "# p(d1)\n",
    "p1=norm1*np.exp(-np.power(d1-dbar1,2)/(2*sigma1*sigma1));\n",
    "# p(d2)\n",
    "p2=norm2*np.exp(-np.power(d2-dbar2,2)/(2*sigma2*sigma2));\n",
    "P=np.matmul(p1,p2.T); # p(d1,d2) = p(d1) p(d2)\n",
    "\n",
    "eda_draw(p);\n",
    "\n",
    "# integrate along d2-axis to get p1\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",
    "# 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": 20,
   "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": [
    "# edapy_03_12: example of a conditional 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.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",
    "# conditional p.d.f.: p(d1|d2) = p(d1,d2)/p2\n",
    "p1g2 = np.divide( P, np.matmul( np.ones((Nd1,1)) ,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,np.ones((Nd2,1)).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": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "d1bar 0.002 d2bar 0.002 covariance 0.910\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 30744 characters]\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edama_03_13: covariace calculated from p.d.f. p(d1,d2)\n",
    "\n",
    "# set up d1-axis\n",
    "Nd1 = 101;\n",
    "d1 = eda_cvec( np.linspace(-2.5,2.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(-2.5,2.5,Nd2) );\n",
    "Dd2 = d2[1,0] - d2[0,0];\n",
    "\n",
    "# create a highly and positively correlated p.d.f.\n",
    "P = np.zeros((Nd1,Nd2));\n",
    "for i in range(Nd1):\n",
    "    for j in range(Nd2):\n",
    "        d_lo = -1.0+d1[i,0];\n",
    "        d_hi =  1.0+d1[i,0];\n",
    "        r2 = 2.0 + d1[i,0]**2 + d2[j,0]**2;\n",
    "        P[i,j] = (d2[j,0]>=d_lo)*(d2[j,0]<=d_hi)/r2;\n",
    "norm = Dd1*Dd2*np.sum(P);\n",
    "P = P/norm;\n",
    "\n",
    "# compute means\n",
    "d1w = np.matmul(d1,np.ones((1,Nd2)));\n",
    "wd2 = np.matmul(np.ones((Nd1,1)),d2.T);\n",
    "d1bar = Dd1*Dd2*np.sum(np.multiply(d1w,P));\n",
    "d2bar = Dd1*Dd2*np.sum(np.multiply(wd2,P));\n",
    "\n",
    "# make the alternating sign function\n",
    "S = np.matmul((d1-d1bar),(d2-d2bar).T);\n",
    "# perform covariance integral\n",
    "sigma12 = (Dd1*Dd2)*np.sum(np.multiply(S,P));\n",
    "\n",
    "print(\"d1bar %.3f d2bar %.3f covariance %.3f\" % (d1bar, d2bar, sigma12) );\n",
    "\n",
    "eda_draw('title s(d1,d2)', S, 'title p(d1,d2)', P, 'title s(d1,d2) p(d1,d2)', np.multiply(S,P));\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 9688 characters]\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edapy_03_14: 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.0,5.0,Nd2) );\n",
    "Dd2 = d2[1,0] - d2[0,0];\n",
    "\n",
    "# 2D normal p.f.d. setup, positively correlated\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",
    "# vectors of ones\n",
    "w1 = np.ones((Nd1,1));\n",
    "w2 = np.ones((Nd2,1));\n",
    "\n",
    "# Normal p.d.f., positive covariance\n",
    "dd1 = d1-dbar1;  # d1 minus its mean\n",
    "dd2 = d2-dbar2;  # d2 minus its mean\n",
    "dd12 = 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(dd12,w2.T)\n",
    "              -0.5*Ci[1,1]*np.matmul(w1,dd22.T)\n",
    "              -Ci[0,1]*np.matmul(dd1,dd2.T));\n",
    "pp=P;\n",
    "\n",
    "# 2D normal p.f.d. setup, negatively correlated\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., negative covariance\n",
    "dd1 = d1-dbar1;  # d1 minus its mean\n",
    "dd2 = d2-dbar2;  # d2 minus its mean\n",
    "dd12 = 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",
    "pn=norm*np.exp(-0.5*Ci[0,0]*np.matmul(dd12,w2.T)\n",
    "               -0.5*Ci[1,1]*np.matmul(w1,dd22.T)\n",
    "               -Ci[0,1]*np.matmul(dd1,dd2.T));\n",
    "\n",
    "# 2D normal p.f.d. setup, uncorrelated\n",
    "dbar1 = 1.5;\n",
    "dbar2 = 2.5;\n",
    "sigma1 = 0.25;\n",
    "sigma2 = 0.25;\n",
    "covar = 0.0;  # zero 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., uncorrelated\n",
    "dd1 = d1-dbar1;  # d1 minus its mean\n",
    "dd2 = d2-dbar2;  # d2 minus its mean\n",
    "dd12 = 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",
    "pu=norm*np.exp(-0.5*Ci[0,0]*np.matmul(dd12,w2.T)\n",
    "               -0.5*Ci[1,1]*np.matmul(w1,dd22.T)\n",
    "               -Ci[0,1]*np.matmul(dd1,dd2.T));\n",
    "\n",
    "eda_draw( 'title cov +', pp, 'title cov -', pn, 'title uncor', pu); "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "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": [
    "# edapy_03_15: Constructing Normal p.d.f. with for loops\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",
    "# 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",
    "# 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",
    "        # vector of data minus their means\n",
    "        dd[:,0] = [ d1[i,0] - dbar1, d2[j,0] - dbar2 ];\n",
    "        # ddT inv(C) dd\n",
    "        E = np.matmul( np.matmul(dd.T,Cdinv), dd ); \n",
    "        # p(d1,d2)\n",
    "        P[i,j] = (1/(2*pi*sqrtdet)) * exp( -0.5*E ); \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",
    "eda_draw( 'title cov +', ppos, 'title cov -', pneg, 'title cov 0', pzero);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.5 2.5\n",
      "cov +: true 0.0500 est 0.0500\n",
      "cov -: true -0.0500 est -0.0500\n",
      "cov 0: true 0.0000 est -0.0000\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 17304 characters]\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edapy_03_16: compute covariance of p(d1,d2)\n",
    "# warning: you must have run the previous cell before this one will work\n",
    "\n",
    "S = np.matmul( (d1-dbar1), (d2-dbar2).T ); # use tensor product to make (d1-dbar1)(d2-dbar2)\n",
    "covpos_est = Dd1*Dd2*np.sum( np.multiply(S, ppos) ); # covariance intergal, + covariance case\n",
    "covneg_est = Dd1*Dd2*np.sum( np.multiply(S, pneg) ); # covariance intergal, - covariance case\n",
    "covzero_est = Dd1*Dd2*np.sum( np.multiply(S, pzero) ); # covariance intergal, 0 covariance case\n",
    "\n",
    "print( 'cov +: true %.4f est %.4f' % (covpos, covpos_est) );\n",
    "print( 'cov -: true %.4f est %.4f' % (covneg, covneg_est) );\n",
    "print( 'cov 0: true %.4f est %.4f' % (covzero, covzero_est) );\n",
    "\n",
    "eda_draw( 'title p(d1,d2)', ppos, 'title s(d1,d2)', S, 'title sp', np.multiply(S,P));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mbar1=3.00 mbar2=2.00\n",
      "varm1=0.12 varm2=0.31 covm=0.19\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 9784 characters]\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edapy_03_17: error propagation from data to model parameters\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, uncorrelated\n",
    "dbar1 = 1.0;\n",
    "dbar2 = 4.0;\n",
    "sigma1 = 0.25;\n",
    "sigma2 = 0.25;\n",
    "covar = 0.0;  # uncorrelated\n",
    "Cd = np.array( [ [sigma1**2, covar], [covar, sigma2**2]] ); # covariance matrix\n",
    "Cdi = la.inv(Cd); # its inverse\n",
    "sqrtdet = sqrt(la.det(Cd)); # square root of its determinant\n",
    "norm=(1/(2*pi*sqrtdet));  # normalization\n",
    "\n",
    "# vectors of ones\n",
    "w1 = np.ones((Nd1,1));\n",
    "w2 = np.ones((Nd2,1));\n",
    "\n",
    "# Normal p.d.f., positive covariance\n",
    "dd1 = d1-dbar1;  # d1 minus its mean\n",
    "dd2 = d2-dbar2;  # d2 minus its mean\n",
    "dd12 = np.power(dd1,2); # (d1 minus its mean) squared\n",
    "dd22 = np.power(dd2,2); # (d2 minus its mean) squared\n",
    "pd=norm*np.exp(-0.5*Cdi[0,0]*np.matmul(dd12,w2.T)-0.5*Cdi[1,1]*np.matmul(w1,dd22.T)-Cdi[0,1]*np.matmul(dd1,dd2.T));\n",
    "\n",
    "# model is related to data by m=Md\n",
    "M = np.array( [ [-1, 1], [-2, 1] ] );\n",
    "\n",
    "# error propagation\n",
    "mbar = np.matmul(M, eda_cvec( dbar1, dbar2) );\n",
    "mbar1 = mbar[0,0];\n",
    "mbar2 = mbar[1,0];\n",
    "Cm = np.matmul(M,np.matmul(Cd,M.transpose()));\n",
    "Cmi = la.inv(Cm);\n",
    "sqrtdet = sqrt(la.det(Cm)); # square root of its determinant\n",
    "norm=(1/(2*pi*sqrtdet));  # normalization\n",
    "\n",
    "# print results\n",
    "print('mbar1=%.2f mbar2=%.2f' % (mbar1,mbar2) );\n",
    "print('varm1=%.2f varm2=%.2f covm=%.2f' % (Cm[0,0], Cm[1,1], Cm[0,1]) );\n",
    "\n",
    "# set up m1-axis\n",
    "Nm1 = 101;\n",
    "m1 = eda_cvec( np.linspace(0.0,5.0, Nm1) );\n",
    "Dm1 = m1[1,0] - m1[0,0];\n",
    "\n",
    "# set up m2-axis\n",
    "Nm2 = 101;\n",
    "m2 = eda_cvec( np.linspace(0.0,5.0,Nm2) );\n",
    "Dm2 = m2[1,0] - m2[0,0];\n",
    "\n",
    "# vectors of ones\n",
    "w1 = np.ones((Nm1,1));\n",
    "w2 = np.ones((Nm2,1));\n",
    "\n",
    "# Normal p.d.f.\n",
    "dm1 = m1-mbar1;  # d1 minus its mean\n",
    "dm2 = m2-mbar2;  # d2 minus its mean\n",
    "dm12 = np.power(dm1,2); # (d1 minus its mean) squared\n",
    "dm22 = np.power(dm2,2); # (d2 minus its mean) squared\n",
    "pm=norm*np.exp(-0.5*Cmi[0,0]*np.matmul(dm12,w2.T)-0.5*Cmi[1,1]*np.matmul(w1,dm22.T)-Cmi[0,1]*np.matmul(dm1,dm2.T));\n",
    "\n",
    "eda_draw( 'title p(d)', pd, 'title p(m)', pm);"
   ]
  },
  {
   "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
}
