{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# eda04_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 argiment.\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",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1.  0.]\n",
      " [ 1.  1.]\n",
      " [ 1.  2.]\n",
      " [ 1.  3.]\n",
      " [ 1.  4.]\n",
      " [ 1.  5.]\n",
      " [ 1.  6.]\n",
      " [ 1.  7.]\n",
      " [ 1.  8.]\n",
      " [ 1.  9.]\n",
      " [ 1. 10.]]\n"
     ]
    }
   ],
   "source": [
    "# eda04_01: data kernel for straight-line linear model\n",
    "\n",
    "# make the auxillary variable x\n",
    "N=11;\n",
    "x = eda_cvec( np.linspace(0,N-1,N) );\n",
    "\n",
    "G = np.zeros((N,2));\n",
    "# G[:,0] as has size 11, but np.ones((N,1)) and x has size 11x1,\n",
    "# so a ravel() is required to match size.  (Yuck!)\n",
    "G[:,0] = np.ones((N,1)).ravel();\n",
    "G[:,1] = x.ravel();\n",
    "\n",
    "print(G);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[  1.   0.   0.]\n",
      " [  1.   1.   1.]\n",
      " [  1.   2.   4.]\n",
      " [  1.   3.   9.]\n",
      " [  1.   4.  16.]\n",
      " [  1.   5.  25.]\n",
      " [  1.   6.  36.]\n",
      " [  1.   7.  49.]\n",
      " [  1.   8.  64.]\n",
      " [  1.   9.  81.]\n",
      " [  1.  10. 100.]]\n"
     ]
    }
   ],
   "source": [
    "# eda04_02: data kernel for quadratic-curve linear model\n",
    "\n",
    "# make the auxillary variable x\n",
    "N=11;\n",
    "x = eda_cvec( np.linspace(0,N-1,N) );\n",
    "\n",
    "G = np.zeros((N,3));\n",
    "# G[:,0] as has size 11, but np.ones((N,1)) and x has size 11x1,\n",
    "# so a ravel() is required to match size.  (Yuck!)\n",
    "G[:,0] = np.ones((N,1)).ravel();\n",
    "G[:,1] = x.ravel();\n",
    "G[:,2] = np.power(x,2).ravel();\n",
    "\n",
    "print(G);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3QAAADnCAYAAABWihTLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAJ90lEQVR4nO3dPY5cTRUG4Kr56ZmxnSJkicgpGUgmJPhWQMYWEHtgC0RIrACJjBUQkBAgkAhICZ2QeOzxTM9P9/0iJLpu2XNpdd/2236erEpn7pQd+VX5nKrDMBQAAADynBz6AAAAAGxHoAMAAAgl0AEAAIQS6AAAAEIJdAAAAKEEOgAAgFACHQAAQCiBDgAAIJRABwAAEEqgAwAACCXQAQAAhBLoAAAAQgl0AAAAoQQ6AACAUAIdAABAKIEOAAAglEAHAAAQSqADAAAIJdABAACEEugAAABCCXQAAAChBDoAAIBQAh0AAEAogQ4AACCUQAcAABBKoAMAAAgl0AEAAIQS6AAAAEIJdAAAAKEEOgAAgFACHQAAQCiBDgAAIJRABwAAEEqgAwAACCXQAQAAhBLoAAAAQgl0AAAAoQQ6AACAUAIdAABAKIEOAAAglEAHAAAQSqADAAAIdbarD9Vah3bv3bt3u/o834DXr1/XQ58BAACSuKEDAAAIJdABAACEEugAAABCCXQAAAChdjYUpef+/n5jfXFxsc9fBwAA8E1xQwcAABBKoAMAAAgl0AEAAISatYdusViMamr1ljQAAMA23NABAACEEugAAABCCXQAAAChBDoAAIBQex2KslwuN9avXr0a1RiKAgAAsB03dAAAAKEEOgAAgFACHQAAQKhZe+jW6/Wo5uREpgQAANiGNAUAABBKoAMAAAgl0AEAAIQS6AAAAELtdSjK/f39xnoYhlFNu+ehcQAAgGnc0AEAAIQS6AAAAEIJdAAAAKFm7aFbrVbjA5zt9QgAAABHyw0dAABAKIEOAAAglEAHAAAQSqADAAAINetQlKenp1HNYrHYWHtYHAAAYBo3dAAAAKEEOgAAgFACHQAAQKi99tAtl8uNde9h8fV6vbHu9dDpqwMAABhzQwcAABBKoAMAAAgl0AEAAIQS6AAAAELNOhSl97D4MAz7PAIAAMDRckMHAAAQSqADAAAIJdABAACE2msP3f39/ca610PXPize66nzsDgAAMCYGzoAAIBQAh0AAEAogQ4AACCUQAcAABBq1qEoDw8Po5rVarWxPj09HdW0Q1EMSQEAAHBDBwAAEEugAwAACCXQAQAAhNprD91yudxYPz4+jmraHrrew+Ltnh46AAAAN3QAAACxBDoAAIBQAh0AAEAogQ4AACDUrA+Lt+tSthuK0qsxKAUAAPjWuKEDAAAIJdABAACEEugAAABCzfqw+MPDw6imfWy87akrpZSTk83cqYcOAADADR0AAEAsgQ4AACCUQAcAABBKoAMAAAg168PiU4aiPD09jWrOzjaPOeXxcUNSAACAY+eGDgAAIJRABwAAEEqgAwAACDXrw+LtupRxX13vYfF2r9cf1+ura+mrAwAAjokbOgAAgFACHQAAQCiBDgAAIJRABwAAEGrWh8XbdSnjoSjtQ+OllHJ+fr6xPj09HdVMGYrSMiQFAABI5oYOAAAglEAHAAAQSqADAAAIJdABAACEmnUoyt3d3bM1vaEoq9Xqi+uek5NxVm2HoPQGqRiUAgAApHBDBwAAEEqgAwAACCXQAQAAhPrqeujah8Z7e1P646bofaftq9NTBwAAfK3c0AEAAIQS6AAAAEIJdAAAAKEEOgAAgFB7HYrSDkHpDUVp95bL5ajm4uJiY312Nj52O+CkN8ykHXiyXq9HNVOGoBiUAgAAfA3c0AEAAIQS6AAAAEIJdAAAAKFmfVj89vZ2VNPu9XroLi8vN9anp6ejmt4j4a22967XCzflO1PoswMAAPbNDR0AAEAogQ4AACCUQAcAABBKoAMAAAh18KEonz59eramfVi8NxRlmyEkve9M+a7BKQAAwNfADR0AAEAogQ4AACCUQAcAABBq1h66tl+ulFJubm6erdlVD90wDF9c977d65drf25XPXU9+uwAAIDPcUMHAAAQSqADAAAIJdABAACEEugAAABC7XUoynK53Fi3A1B6ex8/fhzVnJ+fb6ynDApZr9fP7i0Wi1FNOxTl7Gz8V9QOQekNV2nP2Dvz1L1tGKYCAADHzw0dAABAKIEOAAAglEAHAAAQ6uAPi19fX2+sr66uRjVtH9uU/rDVajXaa3voejXt72r790oZ99n1Hjpv++x6j49PeZB8Sp9dr6bt69u2p04vHgAAfL3c0AEAAIQS6AAAAEIJdAAAAKEEOgAAgFCzPizeezS8HYpyeXk5qmmHh/QeDX96etpYPzw8jGravYuLi1FN+/t7D4u3g1J6Ne2glN4AlCnDVHpDSabUtLZ9xHxXNQAAwO65oQMAAAgl0AEAAIQS6AAAAELN+rB4r4fu/fv3mwfq9KO1PXOPj4/P/q67u7tRzYsXLzbWvX69dq/XZ9f20PUeH2//HFP67Hp7u+qz6/XwTXmgfErNcz+zbQ0AAPBlbugAAABCCXQAAAChBDoAAIBQAh0AAECoWR8Wv7m5GdUsFouN9TAMo5r2QfDb29tRTTtw5eXLl6OadijK1dXVqGbKUJR2r/0zlLLd4JRe3baDU6Y8bD5lcMq+hqv0vgMAAPx//KsaAAAglEAHAAAQSqADAAAINevD4tfX16Oatmeu/ZlSxr132/bHtXvbPiy+TQ/dlJre3rZ9du3PbVvT9rrts1/vzZs3oz0AAODz3NABAACEEugAAABCCXQAAAChBDoAAIBQsw5FmVLz4cOHUc2UISTtXm+YyZTHvtvvTBlcsqvhJr263s9NqWkHlUz5zq4eMe/9rinDVd6+fTvaAwAAPs8NHQAAQCiBDgAAIJRABwAAEKq2D3tv/aFad/MhvlnDMNRDnwEAAJK4oQMAAAgl0AEAAIQS6AAAAEIJdAAAAKF2NhQFAACAebmhAwAACCXQAQAAhBLoAAAAQgl0AAAAoQQ6AACAUAIdAABAKIEOAAAglEAHAAAQSqADAAAIJdABAACEEugAAABCCXQAAAChBDoAAIBQAt2e1Vp/WGv9Q63137XWv9da/1pr/cWhzwUAAOQT6Pao1lpLKX8qpfxlGIY3wzD8tJTyy1LKjw57MgAA4BjUYRgOfYajVWv9rpTym2EYfn7oswAAAMfHDd1+/biU8o9DHwIAADhOAt2Maq2/q7X+s9b6t0OfBQAAyCfQ7de/Sik/+e9iGIZfl1K+K6X84GAnAgAAjoZAt19/LqVc1lp/9T97Lw51GAAA4LgYirJntdbXpZTfllJ+Vkr5TynlUynl98Mw/PGgBwMAAOIJdAAAAKH8l0sAAIBQAh0AAEAogQ4AACCUQAcAABBKoAMAAAgl0AEAAIQS6AAAAEIJdAAAAKEEOgAAgFACHQAAQCiBDgAAIJRABwAAEEqgAwAACCXQAQAAhBLoAAAAQgl0AAAAoQQ6AACAUAIdAABAKIEOAAAglEAHAAAQ6nvSIM38CIdJAgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda04_03: data kernel for quadratic-curve linear model\n",
    "\n",
    "# make the auxillary variable x on varange (0,1)\n",
    "N=32;\n",
    "x = eda_cvec( np.linspace(0,1,N) );\n",
    "Dx = x[1,0]-x[0,0];\n",
    "\n",
    "M=N;  # square G\n",
    "\n",
    "G = np.zeros((N,M));\n",
    "G[:,0] = np.ones((N,1)).ravel();\n",
    "for i in range(1,M):\n",
    "    G[:,i] = np.power(x,i).ravel();\n",
    "\n",
    "eda_draw('title G',G);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 7920 characters]\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda04_04: data kernel for Fourier sum\n",
    "\n",
    "# make the auxillary variable x on varange (0,1)\n",
    "N=32;\n",
    "x = eda_cvec( np.linspace(0,1,N) );\n",
    "Dx = x[1,0]-x[0,0];\n",
    "xmax = x[N-1,0];\n",
    "\n",
    "M=N-2;  # G no suqare, since lambda=infinity case is omitted\n",
    "Mo2 = floor(M/2);\n",
    "\n",
    "# wavelengths lambda\n",
    "lam= eda_cvec( 2*xmax * np.reciprocal( np.linspace(1,Mo2,Mo2) ) );\n",
    "\n",
    "G = np.zeros((N,M));\n",
    "for i in range(1,Mo2):\n",
    "    icos = 2*i;\n",
    "    isin = 2*i+1;\n",
    "    G[:,icos] = np.cos(2*pi*x/lam[i,0]).ravel();\n",
    "    G[:,isin] = np.sin(2*pi*x/lam[i,0]).ravel();\n",
    "\n",
    "eda_draw('title G',G);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 9804 characters]\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda04_05: data kernel viewed as concatenation of columns\n",
    "\n",
    "# make the auxillary variable x on varange (0,1)\n",
    "N=32;\n",
    "x = eda_cvec( np.linspace(0,1,N) );\n",
    "Dx = x[1,0]-x[0,0];\n",
    "xmax = x[N-1,0];\n",
    "\n",
    "M=N-2;  # G no suqare, since lambda=infinity case is omitted\n",
    "Mo2 = floor(M/2);\n",
    "\n",
    "# wavelengths lambda\n",
    "lam = eda_cvec( 2*xmax * np.reciprocal( np.linspace(1,Mo2,Mo2) ) );\n",
    "\n",
    "G = np.zeros((N,M));\n",
    "for i in range(1,Mo2):\n",
    "    icos = 2*i;\n",
    "    isin = 2*i+1;\n",
    "    G[:,icos] = np.cos(2*pi*x/lam[i,0]).ravel();\n",
    "    G[:,isin] = np.sin(2*pi*x/lam[i,0]).ravel();\n",
    "    \n",
    "g1 = eda_cvec( G[:,0] );\n",
    "g2 = eda_cvec( G[:,1] );\n",
    "g3 = eda_cvec( G[:,2] );\n",
    "g4 = eda_cvec( G[:,3] );\n",
    "gM = eda_cvec( G[:,M-1] );\n",
    "\n",
    "eda_draw('title G',G,'=','title g1',g1,'title g2',g2,'title g3',g3,'title g4',g4,'. . .','title gM',gM);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3QAAADnCAYAAABWihTLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAKKUlEQVR4nO3dwatcdxUH8HPUgrYURFF5oqU7ESSUIgVRIgUxKMmi7ly5C2rBFvIXuBSSTTdKwYXgwoXiIhWLXVShQm0Vk1it4CYgdFFEixYXUvm5yFTa6TR9M+/emfed+Xzg8WZe7p13Fz8CX8475/QYowAAAMjzrl0/AAAAAJsR6AAAAEIJdAAAAKEEOgAAgFACHQAAQCiBDgAAIJRABwAAEEqgAwAACCXQAQAAhBLoAAAAQgl0AAAAoQQ6AACAUAIdAABAKIEOAAAglEAHAAAQSqADAAAIJdABAACEEugAAABCCXQAAAChBDoAAIBQAh0AAEAogQ4AACCUQAcAABBKoAMAAAgl0AEAAIQS6AAAAEIJdAAAAKEEOgAAgFACHQAAQCiBDgAAIJRABwAAEEqgAwAACCXQAQAAhBLoAAAAQgl0AAAAoQQ6AACAUAIdAABAKIEOAAAglEAHAAAQSqADAAAIJdABAACEEugAAABCvWeqD+rusfyzS5cu3fY9vNHR0VEf99pV5+3ChQtven/+/Pm33Ld8DYfrpOftzJkzt32/6TXsp3XOGwCsQ4UOAAAglEAHAAAQSqADAAAIJdABAACEmmwoyqqBJ1euXNnoPngnq4abXL16dZLPgWWrBpfcuHFja78LAODtqNABAACEEugAAABCCXQAAAChZu2hW6anjqmsWhq+bJOeuip9dbzVcfra5uqpO+7vBwAOkwodAABAKIEOAAAglEAHAAAQSqADAAAINdlQlFUMSmEumw4usXycTWw6lMTycQBgbip0AAAAoQQ6AACAUAIdAABAqFl76JbpqWNOm/S6WT7OcUzVs2b5OAAwNRU6AACAUAIdAABAKIEOAAAglEAHAAAQaqtDUVYxKIW5WD7OVOYcOGL5OABwEip0AAAAoQQ6AACAUAIdAABAqJ330C3TU8ecLB9nKnP1qFk+DgCsQ4UOAAAglEAHAAAQSqADAAAIJdABAACEOnVDUVYxKIW5WD7OVCwfBwB2QYUOAAAglEAHAAAQSqADAAAIFdFDt0xPHXOyfJypWD4OAMxNhQ4AACCUQAcAABBKoAMAAAgl0AEAAISKHIqyikEpzMXycaZi+TgAMDUVOgAAgFACHQAAQCiBDgAAINTe9NAt01PHnCwfZyqWjwMAJ6FCBwAAEEqgAwAACCXQAQAAhBLoAAAAQu3tUJRVDEphLpaPMxXLxwGAdajQAQAAhBLoAAAAQgl0AAAAoQ6qh26ZnjrmZPk4U7F8HAB4Oyp0AAAAoQQ6AACAUAIdAABAKIEOAAAg1EEPRVnFoBTmYvk4U0lcPl5lUAoAzEGFDgAAIJRABwAAEEqgAwAACKWH7h3oqWNOlo8zFcvHAeAwqdABAACEEugAAABCCXQAAAChBDoAAIBQhqJswKAU5rLN5eMn+X2cfpaPA8BhUKEDAAAIJdABAACEEugAAABC6aGbwFQ9dcf9LA6L5eNMxfJxANg/KnQAAAChBDoAAIBQAh0AAEAogQ4AACCUoSgzsXycuVg+zlQsHweAfCp0AAAAoQQ6AACAUAIdAABAKIEOAAAglKEoW3Lc4SYGpbCJbQ5KMSRlvyUOSjEkBYBDpkIHAAAQSqADAAAIJdABAACE0kO3Q5aPMxfLx5lKYk9dlb46AA6HCh0AAEAogQ4AACCUQAcAABBKoAMAAAhlKMopY1AKc7F8nKkkDkoxJAWAfaVCBwAAEEqgAwAACCXQAQAAhNJDd8rpqWNOm/S6WT7OssSeuip9dQDsBxU6AACAUAIdAABAKIEOAAAglEAHAAAQylCUQAalMBfLx5lK4qAUQ1IASKRCBwAAEEqgAwAACCXQAQAAhNJDtwf01DEny8eZylw9apaPA3DIVOgAAABCCXQAAAChBDoAAIBQAh0AAEAoQ1H2lEEpzMXycaZi+TgAnJwKHQAAQCiBDgAAIJRABwAAEEoP3YHQU8ecLB9nKpaPA8B6VOgAAABCCXQAAAChBDoAAIBQAh0AAEAoQ1EOmEEpzMXycaZi+TgA3J4KHQAAQCiBDgAAIJRABwAAEEoPHf+np445WT7OVBKXjx8dHc322QAcNhU6AACAUAIdAABAKIEOAAAglEAHAAAQylAUbsugFOZi+ThTSVg+fu7cuUk+BwCWqdABAACEEugAAABCCXQAAACh9NCxFj11zGmby8cvXry40X1kSFw+DgCbUKEDAAAIJdABAACEEugAAABCCXQAAAChDEXhxKYalHL58uUpHoc9MufycUNRDkvC8nEA2IQKHQAAQCiBDgAAIJRABwAAEEoPHZPbtKdODx3Hsc3l4+w3y8cB2AcqdAAAAKEEOgAAgFACHQAAQCiBDgAAIFSPMXb9DAAAAGxAhQ4AACCUQAcAABBKoAMAAAgl0AEAAIQS6AAAAEIJdAAAAKEEOgAAgFACHQAAQCiBDgAAIJRABwAAEEqgAwAACCXQAQAAhBLoAAAAQgl0E+ju93b3c919vbv/2N3fXvP+R7v7zrmej/3T3Te7+w/dfa27f7vmvc4ba+nuR7r7hcX/b4+uea/zBgAz6jHGrp8hXnd3Vd01xni1u++oqmeq6pExxrPHvP9mVX16jPG3GR+TPXKSM+O8sY7u/lRV/aiqHqiq/1TVk1X1jTHGX455/81y3gBgNip0Exi3vLp4e8fi601Jubvv7e4/d/cPuvtGd/+4u+/s7m9V1Uer6unufnrLj86ect6Y0Cer6tkxxr/HGK9V1a+q6qE3XuC8AcDuCHQT6e53d/e1qnq5qp4aY/xmxWWfqKrHxxhnquqfVfXNMcZjVfVSVT04xnhwe09MuFFVv+ju33X3xbe5xnljCi9U1dnu/uDiTye/XFUfX3Gd8wYAOyDQTWSM8d8xxn1V9bGqemDxZ0rL/jrG+PXi9Q+r6nNbe0D2zWfHGPdX1Zeq6uHuPrviGueNExtjvFhV36mqp+rWn1ter6rXVlzqvAHADgh0ExtjvFJVv6yqhxYDK65199df/+fly7f6cOyNMcZLi+8vV9VPq+rzzhtzGWN8f4xx/xjjbFX9vapecd4A4HQQ6CbQ3R/q7vcvXr+vqr5QVb8fY9y3+Pre4tJ7uvszi9dfrVvDU6qq/lVVd2/1oYnV3Xd1992vv66qL1bV884bc+nuDy++31NVX6mqx5w3ADgdBLppHNWtpv8bVfV83eqhe2LFdS9W1dcW132gqr67+PnjVfVzQwM4po9U1TPdfb2qnquqn40xnlxxnfPGVH7S3X+qqqtV9fAY4x8rrnHeAGAHrC3Yku6+t6qeGGOs6q2DSTlvbJPzBgC7o0IHAAAQSoUOAAAglAodAABAKIEOAAAglEAHAAAQSqADAAAIJdABAACEEugAAABCCXQAAAChBDoAAIBQAh0AAEAogQ4AACCUQAcAABBKoAMAAAgl0AEAAIQS6AAAAEIJdAAAAKEEOgAAgFACHQAAQCiBDgAAINT/AAourln0eqXJAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda04_06: data kernel fro three point weighted avarage\n",
    "\n",
    "N=32;\n",
    "M=N; # G is square\n",
    "\n",
    "# case 1: 3-point average\n",
    "\n",
    "# symmetric averaging with weights w, where w[0,0] is the central weight\n",
    "Lw = 2;\n",
    "w  = eda_cvec( [2,1] );\n",
    "\n",
    "# normalization, to ensure row sums are unity\n",
    "n = 2*np.sum(w)-w[0,0];\n",
    "w = w/n;\n",
    "\n",
    "# construct data kernel\n",
    "r = np.zeros((N,1));\n",
    "c = np.zeros((M,1));\n",
    "r[0:Lw,0]=w.ravel();\n",
    "c[0:Lw,0]=w.ravel();\n",
    "G3 = la.toeplitz(r,c.T);\n",
    "\n",
    "# case 2: 5-point average\n",
    "\n",
    "# symmetric averaging with weights w, where w[0,0] is the central weight\n",
    "Lw = 3;\n",
    "w = eda_cvec( [3,2,1] );\n",
    "\n",
    "# normalization, to ensure row sums are unity\n",
    "n = 2*np.sum(w)-w[0,0];\n",
    "w = w/n;\n",
    "\n",
    "# construct data kernel\n",
    "r = np.zeros((N,1));\n",
    "c = np.zeros((M,1));\n",
    "r[0:Lw,0]=w.ravel();\n",
    "c[0:Lw,0]=w.ravel();\n",
    "G5 = la.toeplitz(r,c.T);\n",
    "\n",
    "# case 3: 5-point average\n",
    "\n",
    "# symmetric averaging with weights w, where w[0,0] is the central weight\n",
    "Lw = 5;\n",
    "w = eda_cvec( [5,4,3,2,1] );\n",
    "\n",
    "# normalization, to ensure row sums are unity\n",
    "n = 2*np.sum(w)-w[0,0];\n",
    "w = w/n;\n",
    "\n",
    "# construct data kernel\n",
    "r = np.zeros((N,1));\n",
    "c = np.zeros((M,1));\n",
    "r[0:Lw,0]=w.ravel();\n",
    "c[0:Lw,0]=w.ravel();\n",
    "G9 = la.toeplitz(r,c.T);\n",
    "\n",
    "eda_draw('title 3-pt',G3,'title 5-pt',G5,'title 9-pt',G9);\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 18424 characters]\n",
      "text/plain": [
       "<Figure size 576x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda04_07\n",
    "# reads in noisy, linear data from file\n",
    "#  and computes and plots prediction error for a particular\n",
    "# choice of intercept and slope, m1 and m2\n",
    "\n",
    "# model paramerters are intercept and slope\n",
    "M=2;\n",
    "m = eda_cvec( [1.1, 2.3] );\n",
    "\n",
    "# load data from file\n",
    "T = np.genfromtxt('linedata01.txt', delimiter='\\t')\n",
    "[N, K]=T.shape;\n",
    "x = eda_cvec( T[:,0] );\n",
    "dobs = eda_cvec( T[:,1] );\n",
    "\n",
    "# linear data kernel\n",
    "G = np.zeros((N,2));\n",
    "G[:,0] = np.ones((N,1)).ravel();\n",
    "G[:,1] = x.ravel();\n",
    "dpre = np.matmul(G,m);\n",
    "\n",
    "#  plot data\n",
    "fig1 = plt.figure(1,figsize=(8,4));\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "plt.axis([1.1*np.min(x), 1.1*np.max(x), 1.1*np.min(dobs), 1.1*np.max(dobs)])\n",
    "plt.plot(x,dobs,'ko');\n",
    "plt.plot(x,dpre,'k-');\n",
    "plt.xlabel('x)');\n",
    "plt.ylabel('d');\n",
    "titlestr = \"intercept %.2f slope %.2f\" % (m[0,0], m[1,0]);\n",
    "plt.title(titlestr);\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 10420 characters]\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 18428 characters]\n",
      "text/plain": [
       "<Figure size 576x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda04_08: grid search for intersept and slope\n",
    "# for problem of fittng straight line to noisy data\n",
    "\n",
    "# model paramerters\n",
    "M=2;\n",
    "m = np.zeros((M,1));\n",
    "\n",
    "# load data from file\n",
    "T = np.genfromtxt('linedata01.txt', delimiter='\\t')\n",
    "[N, K]=T.shape;\n",
    "x = eda_cvec( T[:,0] );\n",
    "dobs = eda_cvec( T[:,1] );\n",
    "\n",
    "# linear data kernel\n",
    "G = np.zeros((N,2));\n",
    "G[:,0] = np.ones((N,1)).ravel();\n",
    "G[:,1] = x.ravel();\n",
    "dpre = np.matmul(G,m);\n",
    "\n",
    "# define grid of intercepts and slopes\n",
    "L1 = 101;\n",
    "L2 = 101;\n",
    "m1min=0;\n",
    "m1max=4;\n",
    "m2min=0;\n",
    "m2max=4;\n",
    "m1 = eda_cvec( np.linspace(m1min,m1max,L1) );\n",
    "m2 = eda_cvec( np.linspace(m2min,m2max,L2) );\n",
    "e = np.zeros((N,1));\n",
    "E = np.zeros((L1,L2));\n",
    "\n",
    "for i in range(L1):\n",
    "    for j in range(L2):\n",
    "        m = eda_cvec( [m1[i,0], m2[j,0]] );\n",
    "        dpre = np.matmul(G,m);\n",
    "        e = dobs - dpre;\n",
    "        E[i,j] = np.matmul(e.T,e);\n",
    "\n",
    "# best estimate is the one with smallest error\n",
    "k = np.argmin( E, axis=0 );\n",
    "Etemp = np.min( E, axis=0 );\n",
    "j = np.argmin(Etemp);\n",
    "i=k[j];\n",
    "Ebest = E[i,j];\n",
    "mbest = eda_cvec( [ m1[i,0],  m2[j,0] ] );\n",
    "dpre = np.matmul(G,mbest);\n",
    "\n",
    "eda_draw('title E(m1,m2)',E);\n",
    "\n",
    "#  plot data\n",
    "fig2 = plt.figure(2,figsize=(8,4));\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "plt.axis([1.1*np.min(x), 1.1*np.max(x), 1.1*np.min(dobs), 1.1*np.max(dobs)])\n",
    "plt.plot(x,dobs,'ko');\n",
    "plt.plot(x,dpre,'k-');\n",
    "plt.xlabel('x)');\n",
    "plt.ylabel('d');\n",
    "titlestr = \"intercept %.2f slope %.2f\" % (mbest[0,0], mbest[1,0]);\n",
    "plt.title(titlestr);\n",
    "plt.show();\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 27844 characters]\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda04_09: grid search for intersept and slope\n",
    "# for problem of fittng straight line to noisy data\n",
    "\n",
    "# model paramerters\n",
    "M=2;\n",
    "m = np.zeros((M,1));\n",
    "\n",
    "# load data from file\n",
    "T = np.genfromtxt('linedata01.txt', delimiter='\\t')\n",
    "[N, K]=T.shape;\n",
    "x = eda_cvec( T[:,0] );\n",
    "dobs = eda_cvec( T[:,1] );\n",
    "\n",
    "# Part 1: shift left\n",
    "xshift=-2;\n",
    "xs = x+xshift;\n",
    "\n",
    "# linear data kernel\n",
    "G = np.zeros((N,2));\n",
    "G[:,0] = np.ones((N,1)).ravel();\n",
    "G[:,1] = xs.ravel();\n",
    "dpre = np.matmul(G,m);\n",
    "\n",
    "# define grid of intercepts and slopes\n",
    "L1 = 101;\n",
    "L2 = 101;\n",
    "m1min=-10;\n",
    "m1max=10;\n",
    "m2min=-10;\n",
    "m2max=10;\n",
    "m1 = eda_cvec( np.linspace(m1min,m1max,L1) );\n",
    "m2 = eda_cvec( np.linspace(m2min,m2max,L2) );\n",
    "e = np.zeros((N,1));\n",
    "E = np.zeros((L1,L2));\n",
    "\n",
    "for i in range(L1):\n",
    "    for j in range(L2):\n",
    "        m[:,0] = [m1[i], m2[j]];\n",
    "        dpre = np.matmul(G,m);\n",
    "        e = dobs - dpre;\n",
    "        E0 = np.matmul(e.T,e);\n",
    "        E[i,j] = E0;\n",
    "    \n",
    "# best estimate is the one with smallest error\n",
    "k = np.argmin( E, axis=0 );\n",
    "Etemp = np.min( E, axis=0 );\n",
    "j = np.argmin(Etemp);\n",
    "i=k[j];\n",
    "Ebest = E[i,j];\n",
    "mbest = eda_cvec( [ m1[i,0],  m2[j,0] ] );\n",
    "dpre = np.matmul(G,mbest);\n",
    "Esl = E;\n",
    "\n",
    "\n",
    "# Part 2: no shift \n",
    "xshift=0;\n",
    "xs = x+xshift;\n",
    "\n",
    "# linear data kernel\n",
    "G = np.zeros((N,2));\n",
    "G[:,0] = np.ones((N,1)).ravel();\n",
    "G[:,1] = xs.ravel();\n",
    "dpre = np.matmul(G,m);\n",
    "\n",
    "# define grid of intercepts and slopes\n",
    "L1 = 101;\n",
    "L2 = 101;\n",
    "m1min=-10;\n",
    "m1max=10;\n",
    "m2min=-10;\n",
    "m2max=10;\n",
    "m1 = eda_cvec( np.linspace(m1min,m1max,L1) );\n",
    "m2 = eda_cvec( np.linspace(m2min,m2max,L2) );\n",
    "e = np.zeros((N,1));\n",
    "E = np.zeros((L1,L2));\n",
    "\n",
    "for i in range(L1):\n",
    "    for j in range(L2):\n",
    "        m = eda_cvec( [m1[i,0], m2[j,0]] );\n",
    "        dpre = np.matmul(G,m);\n",
    "        e = dobs-dpre;\n",
    "        E0 = np.matmul(e.T,e);\n",
    "        E[i,j] = E0;\n",
    "    \n",
    "# best estimate is the one with smallest error\n",
    "k = np.argmin( E, axis=0 );\n",
    "Etemp = np.min( E, axis=0 );\n",
    "j = np.argmin(Etemp);\n",
    "i=k[j];\n",
    "Ebest = E[i,j];\n",
    "mbest = eda_cvec( [ m1[i,0],  m2[j,0] ] );\n",
    "dpre = np.matmul(G,mbest);\n",
    "Ens = E;\n",
    "\n",
    "# Part 3: shift right\n",
    "xshift=2;\n",
    "xs = x+xshift;\n",
    "\n",
    "# linear data kernel\n",
    "G = np.zeros((N,2));\n",
    "G[:,0] = np.ones((N,1)).ravel();\n",
    "G[:,1] = xs.ravel();\n",
    "dpre = np.matmul(G,m);\n",
    "\n",
    "# define grid of intercepts and slopes\n",
    "L1 = 101;\n",
    "L2 = 101;\n",
    "m1min=-10;\n",
    "m1max=10;\n",
    "m2min=-10;\n",
    "m2max=10;\n",
    "m1 = eda_cvec( np.linspace(m1min,m1max,L1) );\n",
    "m2 = eda_cvec( np.linspace(m2min,m2max,L2) );\n",
    "e = np.zeros((N,1));\n",
    "E = np.zeros((L1,L2));\n",
    "\n",
    "for i in range(L1):\n",
    "    for j in range(L2):\n",
    "        m[:,0] = [m1[i], m2[j]];\n",
    "        dpre = np.matmul(G,m);\n",
    "        e = dobs-dpre;\n",
    "        E0 = np.matmul(e.T,e);\n",
    "        E[i,j] = E0;\n",
    "    \n",
    "# best estimate is the one with smallest error\n",
    "k = np.argmin( E, axis=0 );\n",
    "Etemp = np.min( E, axis=0 );\n",
    "j = np.argmin(Etemp);\n",
    "i=k[j];\n",
    "Ebest = E[i,j];\n",
    "mbest = eda_cvec( [ m1[i,0],  m2[j,0] ] );\n",
    "dpre = np.matmul(G,mbest);\n",
    "Esr = E;\n",
    "\n",
    "# lake log of error to reduce range of variation for plotting purposes\n",
    "eda_draw('title E shift left',np.log(Esl),'title E no shift',np.log(Ens),'title E shift right',np.log(Esr));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "intercept 0.941 +/- 0.974 (%95)\n",
      "slope     2.019 +/- 0.308 (%95)\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 19200 characters]\n",
      "text/plain": [
       "<Figure size 576x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda04_10: least-squares for straight line fit\n",
    "\n",
    "# model paramerters\n",
    "M=2;\n",
    "mest = np.zeros((M,1));\n",
    "\n",
    "# load data from file\n",
    "T = np.genfromtxt('linedata01.txt', delimiter='\\t')\n",
    "[N, K]=T.shape;\n",
    "x = eda_cvec( T[:,0] );\n",
    "dobs = eda_cvec( T[:,1] );\n",
    "\n",
    "# linear data kernel\n",
    "G = np.zeros((N,2));\n",
    "G[:,0] = np.ones((N,1)).ravel();\n",
    "G[:,1] = x.ravel();\n",
    "\n",
    "# siple least squares solution\n",
    "GTG = np.matmul(G.T,G);\n",
    "mest = la.solve( GTG, np.matmul(G.T, dobs)  );\n",
    "\n",
    "# eda04_10: variance calculation\n",
    "dpre = np.matmul(G, mest); # predicted data\n",
    "e = dobs - dpre;           # error\n",
    "E = np.matmul(e.T,e);      # total error\n",
    "sigmad2 = E/(N-M);         # posterior varaince of data\n",
    "Cm = sigmad2*la.inv(GTG);  # variance of model parameters\n",
    "sigmam1 = sqrt( Cm[0,0] ); # sqrt variance of model parameter 1\n",
    "sigmam2 = sqrt( Cm[1,1] ); # sqrt variance of model parameter 2\n",
    "\n",
    "# print model parameters and theoir 95% confidence intervals\n",
    "str1 = \"intercept %.3f +/- %.3f (%%95)\" % (mest[0,0],2*sigmam1);\n",
    "str2 = \"slope     %.3f +/- %.3f (%%95)\" % (mest[1,0],2*sigmam2);\n",
    "print(str1);\n",
    "print(str2);\n",
    "\n",
    "#  plot data\n",
    "fig2 = plt.figure(2,figsize=(8,6));\n",
    "ax1 = plt.subplot(2,1,1);\n",
    "plt.axis([1.1*np.min(x), 1.1*np.max(x), 1.1*np.min(dobs), 1.1*np.max(dobs)])\n",
    "plt.plot(x,dobs,'ko');\n",
    "plt.plot(x,dpre,'k-');\n",
    "plt.xlabel('x');\n",
    "plt.ylabel('d');\n",
    "titlestr = \"intercept %.2f slope %.2f\" % (mest[0,0], mest[1,0]);\n",
    "plt.title(titlestr);\n",
    "#  plot error, with different vertical scale\n",
    "ax2 = plt.subplot(2,1,2);\n",
    "plt.axis([1.1*np.min(x), 1.1*np.max(x), -3*np.amax(e), 3*np.amax(e)])\n",
    "plt.plot(x,e,'ko');\n",
    "plt.plot( [1.1*np.min(x), 1.1*np.max(x)],[0,0],'k:');\n",
    "plt.xlabel('x');\n",
    "plt.ylabel('e');\n",
    "plt.show();\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "slope -8.590021e-05 +/- 1.251757e-05  C per day (95%)\n",
      "slope -3.137505e-02 +/- 4.572041e-03  C per year (95%)\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 73304 characters]\n",
      "text/plain": [
       "<Figure size 1152x432 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda04_11 least squares for secular trend and annual cycles in black rock forest data\n",
    "\n",
    "# load data from file\n",
    "D = np.genfromtxt('brf_temp.txt', delimiter='\\t')\n",
    "[Nraw, K]=np.shape(D);\n",
    "traw = eda_cvec( D[:,0] );\n",
    "draw = eda_cvec( D[:,1] );\n",
    "\n",
    "# select only good data\n",
    "# the where() returns a list of indices that satisfy the condition\n",
    "# draw not equal to zero\n",
    "# and\n",
    "# draw great than (-40)\n",
    "# and\n",
    "# draw < 38\n",
    "k = np.where( (draw[:,0]!=0)&(draw[:,0]>-40)&(draw[:,0]<38) );\n",
    "(K,N)=np.shape(k);\n",
    "# now select only those elements of traw and draw that passes the test\n",
    "t = eda_cvec( traw[k] );\n",
    "dobs = eda_cvec( draw[k] );\n",
    "\n",
    "Ty = 365.25;\n",
    "M= 4;\n",
    "G = np.zeros((N,M));\n",
    "G[:,0] = np.ones( (N,1) ).ravel();\n",
    "G[:,1] = t.ravel();\n",
    "G[:,2] = np.cos(2*pi*t/Ty).ravel();\n",
    "G[:,3] = np.sin(2*pi*t/Ty).ravel();\n",
    "\n",
    "# simple least squares solution\n",
    "GTG = np.matmul(G.T,G);\n",
    "mest = la.solve( GTG, np.matmul(G.T, dobs)  );\n",
    "\n",
    "# eda04_10: variance calculation\n",
    "dpre = np.matmul(G, mest); # predicted data\n",
    "e = dobs - dpre;           # error\n",
    "E = np.matmul(e.T,e);      # total error\n",
    "sigmad2 = E/(N-M);         # posterior varaince of data\n",
    "Cm = sigmad2*la.inv(GTG);  # variance of model parameters\n",
    "sigmam1 = sqrt( Cm[0,0] ); # sqrt variance of model parameter 1\n",
    "sigmam2 = sqrt( Cm[1,1] ); # sqrt variance of model parameter 2\n",
    "sigmam3 = sqrt( Cm[2,2] ); # sqrt variance of model parameter 3\n",
    "sigmam4 = sqrt( Cm[3,3] ); # sqrt variance of model parameter 4\n",
    "\n",
    "# write out slope and its confidence\n",
    "slopestr1 = \"slope %e +/- %e  C per day (95%%)\" % (mest[1,0], sigmam2);\n",
    "print(slopestr1);\n",
    "\n",
    "# convert to per year\n",
    "slopestr2 = \"slope %e +/- %e  C per year (95%%)\" % (Ty*mest[1,0], Ty*sigmam2);\n",
    "print(slopestr2);\n",
    "\n",
    "#  plot data\n",
    "fig1 = plt.figure(1,figsize=(16,6));\n",
    "ax1 = plt.subplot(3,1,1);\n",
    "plt.axis([t[0,0], t[N-1,0], -40, 40])\n",
    "plt.plot(t,dobs,'k.');\n",
    "plt.xlabel('time (days)');\n",
    "plt.ylabel('temperature (C)');\n",
    "ax2 = plt.subplot(3,1,2);\n",
    "plt.axis([t[0,0], t[N-1,0], -40, 40])\n",
    "plt.plot(t,dpre,'k-');\n",
    "plt.xlabel('time (days)');\n",
    "plt.ylabel('temperature (C)');\n",
    "ax3 = plt.subplot(3,1,3);\n",
    "plt.axis([t[0,0], t[N-1,0], -40, 40])\n",
    "plt.plot(t,e,'k-');\n",
    "plt.xlabel('time (days)');\n",
    "plt.ylabel('error (C)');\n"
   ]
  },
  {
   "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
}
