LA4DS_ch04.ipynb

{
 "nbformat": 4,
 "nbformat_minor": 0,
 "metadata": {
  "colab": {
   "name": "LA4DS_ch04.ipynb",
   "provenance": [],
   "collapsed_sections": [],
   "mount_file_id": "1CNEC5tpXXVD9_kPZZ72mm-uy6dwQQfhF",
   "authorship_tag": "ABX9TyOCDwyXDHmxkg4K+oZIczNv"
  },
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3"
  },
  "language_info": {
   "name": "python"
  }
 },
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# Practical Linear Algebra for Data Science\n",
    "## Mike X Cohen (sincxpress.com)\n",
    "### https://www.oreilly.com/library/view/practical-linear-algebra/9781098120603/\n",
    "\n",
    "#### Code for chapter 4"
   ],
   "metadata": {
    "id": "SbGFWGzkd44U"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# NOTE: these lines define global figure properties used for publication.\n",
    "import matplotlib_inline.backend_inline\n",
    "matplotlib_inline.backend_inline.set_matplotlib_formats('svg') # print figures in svg format\n",
    "plt.rcParams.update({'font.size':14}) # set global font size"
   ],
   "metadata": {
    "id": "IvinJbZP_CLv"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "\n",
    "N = 30\n",
    "\n",
    "# correlated random variables\n",
    "x = np.linspace(0,10,N) + np.random.randn(N)\n",
    "y = x + np.random.randn(N)\n",
    "\n",
    "\n",
    "# set up figure\n",
    "_,axs = plt.subplots(2,2,figsize=(6,6))\n",
    "\n",
    "# positive correlation\n",
    "axs[0,0].plot(x,y,'ko')\n",
    "axs[0,0].set_title('Positive correlation',fontweight='bold')\n",
    "axs[0,0].set_xlabel('Variable x')\n",
    "axs[0,0].set_ylabel('Variable y')\n",
    "axs[0,0].set_xticks([])\n",
    "axs[0,0].set_yticks([])\n",
    "axs[0,0].axis('square')\n",
    "\n",
    "\n",
    "# negative correlation\n",
    "axs[0,1].plot(x,-y,'ko')\n",
    "axs[0,1].set_title('Negative correlation',fontweight='bold')\n",
    "axs[0,1].set_xlabel('Variable x')\n",
    "axs[0,1].set_ylabel('Variable y')\n",
    "axs[0,1].set_xticks([])\n",
    "axs[0,1].set_yticks([])\n",
    "axs[0,1].axis('square')\n",
    "\n",
    "\n",
    "# zero correlation, part 1\n",
    "axs[1,0].plot(np.random.randn(N),np.random.randn(N),'ko')\n",
    "axs[1,0].set_title('Zero correlation',fontweight='bold')\n",
    "axs[1,0].set_xlabel('Variable x')\n",
    "axs[1,0].set_ylabel('Variable y')\n",
    "axs[1,0].set_xticks([])\n",
    "axs[1,0].set_yticks([])\n",
    "axs[1,0].axis('square')\n",
    "\n",
    "\n",
    "# zero correlation, part 2\n",
    "x = np.cos(np.linspace(0,2*np.pi,N)) + np.random.randn(N)/20\n",
    "y = np.sin(np.linspace(0,2*np.pi,N)) + np.random.randn(N)/20\n",
    "axs[1,1].plot(x,y,'ko')\n",
    "axs[1,1].set_title('Zero correlation',fontweight='bold')\n",
    "axs[1,1].set_xlabel('Variable x')\n",
    "axs[1,1].set_ylabel('Variable y')\n",
    "axs[1,1].set_xticks([])\n",
    "axs[1,1].set_yticks([])\n",
    "axs[1,1].axis('square')\n",
    "\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('Figure_04_01.png',dpi=300) # write out the fig to a file\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "OU4BBQmyUdOE"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    ""
   ],
   "metadata": {
    "id": "Vg5BHVwt0gQG"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "### Note: The code for k-means is in Exercise 7 below."
   ],
   "metadata": {
    "id": "-puIEx1q0gtR"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    ""
   ],
   "metadata": {
    "id": "klzkPQ8D0gwc"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 1"
   ],
   "metadata": {
    "id": "nuvWclOc0g1J"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# the function\n",
    "def corrAndCosine(x,y):\n",
    "\n",
    "  # compute cosine similarity\n",
    "  num = np.dot(x,y) # numerator\n",
    "  den = np.linalg.norm(x) * np.linalg.norm(y) # denominator\n",
    "  cos = num / den\n",
    "\n",
    "  # compute correlation (similar to above but mean-centered!)\n",
    "  xm  = x-np.mean(x)\n",
    "  ym  = y-np.mean(y)\n",
    "  num = np.dot(xm,ym) # numerator\n",
    "  den = np.linalg.norm(xm) * np.linalg.norm(ym) # denominator\n",
    "  cor = num / den\n",
    "\n",
    "  return cor,cos\n",
    "\n",
    "\n",
    "# test it\n",
    "a = np.random.randn(15)\n",
    "b = np.random.randn(15)\n",
    "\n",
    "# compute the correlation and cosine\n",
    "r,c = corrAndCosine(a,b)\n",
    "\n",
    "# confirm that the correlation matches with np.corrcoef\n",
    "print(r,np.corrcoef(a,b)[0,1])"
   ],
   "metadata": {
    "id": "l9cyAOCn_lt-"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# compare r and c without mean-centering\n",
    "a = np.random.randn(15) + 10 # note the offset!\n",
    "b = np.random.randn(15)\n",
    "\n",
    "# mean-center\n",
    "aNoMean = a - np.mean(a)\n",
    "bNoMean = b - np.mean(b)\n",
    "\n",
    "\n",
    "# show the results with and without mean-centering\n",
    "print('Without mean-centering (should differ):')\n",
    "print( np.round(corrAndCosine(a,b),4) )\n",
    "print(' ')\n",
    "\n",
    "print('With mean-centering (should be the same):')\n",
    "print( np.round(corrAndCosine(aNoMean,bNoMean),4) )\n",
    "\n",
    "# NOTE: In the printing code above, I rounded to 4 significant digits just for visual clarity."
   ],
   "metadata": {
    "id": "TCHVd-TBCOc-"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Exercise 2"
   ],
   "metadata": {
    "id": "XfDHd5pgIG6I"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# create the variables\n",
    "a = np.arange(4,dtype=float)\n",
    "offsets = np.arange(-50,51)\n",
    "\n",
    "# initialize the results\n",
    "results = np.zeros((len(offsets),2))\n",
    "\n",
    "# run the simulation!\n",
    "for i in range(len(offsets)):\n",
    "    results[i,:] = corrAndCosine(a,a+offsets[i])\n",
    "\n",
    "\n",
    "# plot the results!\n",
    "plt.figure(figsize=(8,4))\n",
    "h = plt.plot(offsets,results)\n",
    "h[0].set_color('k')\n",
    "h[0].set_marker('o')\n",
    "h[1].set_color([.7,.7,.7])\n",
    "h[1].set_marker('s')\n",
    "\n",
    "plt.xlabel('Mean offset')\n",
    "plt.ylabel('r or c')\n",
    "plt.legend(['Pearson','Cosine sim.'])\n",
    "plt.savefig('Figure_04_02.png',dpi=300) # write out the fig to a file\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "GL26OrcAp6Hm"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    ""
   ],
   "metadata": {
    "id": "L5YGtmh4q1G6"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Exercise 3"
   ],
   "metadata": {
    "id": "VrJ4GQmzI0RU"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# import the function\n",
    "from scipy.stats import pearsonr\n",
    "\n",
    "# inspect the source code\n",
    "??pearsonr"
   ],
   "metadata": {
    "id": "cHjKrhM1DOdD"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    ""
   ],
   "metadata": {
    "id": "ZXMHlxZu00Fs"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Exercise 4"
   ],
   "metadata": {
    "id": "-mp0Tn-6Mx7E"
   }
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "oo6WdWnsdcne"
   },
   "source": [
    "# a bare-bones correlation function\n",
    "def rho(x,y):\n",
    "  xm = x-np.mean(x)\n",
    "  ym = y-np.mean(y)\n",
    "  n  = np.dot(xm,ym)\n",
    "  d  = np.linalg.norm(xm) * np.linalg.norm(ym)\n",
    "  return n/d\n",
    "\n",
    "\n",
    "# import the time library\n",
    "import time\n",
    "\n",
    "# experiment parameters\n",
    "numIters  = 1000\n",
    "varLength =  500\n",
    "\n",
    "# clock my custom-written function\n",
    "tic = time.time()\n",
    "for i in range(numIters):\n",
    "  x = np.random.randn(varLength,2)\n",
    "  rho(x[:,0],x[:,1])\n",
    "t1 = time.time() - tic\n",
    "\n",
    "\n",
    "# now for numpy's corrcoef function\n",
    "tic = time.time()\n",
    "for i in range(numIters):\n",
    "  x = np.random.randn(varLength,2)\n",
    "  pearsonr(x[:,0],x[:,1])\n",
    "t2 = time.time() - tic\n",
    "\n",
    "\n",
    "# print the results!\n",
    "# Note: time() returns seconds, so I multiply by 1000 for ms\n",
    "print(f'My function took {t1*1000:.2f} ms')\n",
    "print(f'   pearsonr took {t2*1000:.2f} ms')"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    ""
   ],
   "metadata": {
    "id": "gy0-joJm03P4"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 5"
   ],
   "metadata": {
    "id": "bvHwvToaN1tv"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# create the kernel (in the book figure I used +1.5)\n",
    "kernel = np.array([-1,1])\n",
    "\n",
    "# and the \"signal\" (a plateau)\n",
    "signal = np.zeros(30)\n",
    "signal[10:20] = 1\n",
    "\n",
    "\n",
    "# plot them\n",
    "_,axs = plt.subplots(1,2,figsize=(12,4))\n",
    "axs[0].plot(kernel,'ks-')\n",
    "axs[0].set_title('Kernel')\n",
    "axs[0].set_xlim([-15,15])\n",
    "\n",
    "axs[1].plot(signal,'ks-')\n",
    "axs[1].set_title('Time series signal')\n",
    "\n",
    "plt.savefig('Figure_04_04ab.png',dpi=300)\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "hF-G_i8IN3uU"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# initialize the feature map as zeros\n",
    "featureMap = np.zeros(len(signal))\n",
    "\n",
    "# loop over the signal and do template-matching (via dot products!)\n",
    "for t in range(1,len(signal)-1):\n",
    "  featureMap[t] = np.dot(kernel,signal[t-1:t+1])\n",
    "\n",
    "\n",
    "# plot the result\n",
    "_,axs = plt.subplots(1,2,figsize=(12,4))\n",
    "axs[0].plot(kernel,'ks-')\n",
    "axs[0].set_title('Kernel')\n",
    "axs[0].set_xlim([-15,15])\n",
    "\n",
    "\n",
    "axs[1].plot(signal,'ks-',label='Signal',linewidth=3)\n",
    "markers,stemlines,_ = axs[1].stem(range(len(featureMap)),featureMap,\n",
    "                                  basefmt=' ',linefmt='',markerfmt='o',\n",
    "                                  label='Edge detection')\n",
    "\n",
    "plt.setp(stemlines,'color',[.7,.7,.7])\n",
    "plt.setp(markers,'color',[.7,.7,.7])\n",
    "\n",
    "axs[1].legend()\n",
    "plt.savefig('Figure_04_04c.png',dpi=300)\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "VU4mT29IN2EE"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    ""
   ],
   "metadata": {
    "id": "FHS9GLbFN2Hd"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 6"
   ],
   "metadata": {
    "id": "MboACnUwN2J-"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# define the kernel (a sorta-kinda Gaussian)\n",
    "kernel = np.array([0,.1,.3,.8,1,.8,.3,.1,0])\n",
    "kernel = kernel / np.sum(kernel)\n",
    "\n",
    "# some handy length parameters\n",
    "Nkernel = len(kernel)\n",
    "halfKrn = Nkernel//2\n",
    "\n",
    "\n",
    "# and the signal\n",
    "Nsignal = 100\n",
    "timeseries = np.random.randn(Nsignal)\n",
    "\n",
    "\n",
    "# plot them\n",
    "_,axs = plt.subplots(1,2,figsize=(12,4))\n",
    "axs[0].plot(kernel,'ks-')\n",
    "axs[0].set_title('Kernel')\n",
    "axs[0].set_xlim([-1,Nsignal])\n",
    "\n",
    "axs[1].plot(timeseries,'ks-')\n",
    "axs[1].set_title('Time series signal')\n",
    "\n",
    "plt.savefig('Figure_04_06ab.png',dpi=300)\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "LyLAL7L0RBDt"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# make a copy of the signal for filtering\n",
    "filtsig = timeseries.copy()\n",
    "\n",
    "# loop over the signal time points\n",
    "for t in range(halfKrn+1,Nsignal-halfKrn):\n",
    "  filtsig[t] = np.dot(kernel,timeseries[t-halfKrn-1:t+halfKrn])\n",
    "\n",
    "\n",
    "# and plot\n",
    "_,axs = plt.subplots(1,2,figsize=(12,4))\n",
    "axs[0].plot(kernel,'ks-')\n",
    "axs[0].set_title('Kernel')\n",
    "axs[0].set_xlim([-1,Nsignal])\n",
    "\n",
    "axs[1].plot(timeseries,color='k',label='Original',linewidth=1)\n",
    "axs[1].plot(filtsig,'--',color=[.6,.6,.6],label='Smoothed',linewidth=2)\n",
    "axs[1].legend()\n",
    "\n",
    "plt.savefig('Figure_04_06c.png',dpi=300)\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "1HDFchK6RBGm"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 7"
   ],
   "metadata": {
    "id": "fpe9_N2_RBJZ"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# define the kernel (a sorta-kinda Gaussian)\n",
    "kernel = np.array([0,.1,.3,.8,-1,.8,.3,.1,0])\n",
    "kernel /= np.sum(kernel)\n",
    "kernel -= np.mean(kernel)\n",
    "\n",
    "# plot them\n",
    "_,axs = plt.subplots(1,2,figsize=(12,4))\n",
    "axs[0].plot(kernel,'s-')\n",
    "axs[0].set_title('Kernel')\n",
    "axs[0].set_xlim([-1,Nsignal])\n",
    "\n",
    "axs[1].plot(timeseries,'s-')\n",
    "axs[1].set_title('Time series signal')\n",
    "plt.show()\n",
    "\n",
    "\n",
    "\n",
    "# loop over the signal time points\n",
    "filtsig2 = timeseries.copy()\n",
    "for t in range(halfKrn+1,Nsignal-halfKrn):\n",
    "  filtsig2[t] = np.dot(kernel,timeseries[t-halfKrn-1:t+halfKrn])\n",
    "\n",
    "plt.plot(timeseries,color='k',label='Original',linewidth=1)\n",
    "plt.plot(filtsig2,color=[.9,.2,.7],label='Sharpened',linewidth=1)\n",
    "plt.legend()\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "poWxOslmcbAz"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    ""
   ],
   "metadata": {
    "id": "t9Ud0MCwcbGf"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 8"
   ],
   "metadata": {
    "id": "7w7jf9jXcbPL"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "## Create data\n",
    "nPerClust = 50\n",
    "\n",
    "# blur around centroid (std units)\n",
    "blur = 1\n",
    "\n",
    "# XY centroid locations\n",
    "A = [  1, 1 ]\n",
    "B = [ -3, 1 ]\n",
    "C = [  3, 3 ]\n",
    "\n",
    "# generate data\n",
    "a = [ A[0]+np.random.randn(nPerClust)*blur , A[1]+np.random.randn(nPerClust)*blur ]\n",
    "b = [ B[0]+np.random.randn(nPerClust)*blur , B[1]+np.random.randn(nPerClust)*blur ]\n",
    "c = [ C[0]+np.random.randn(nPerClust)*blur , C[1]+np.random.randn(nPerClust)*blur ]\n",
    "\n",
    "# concatanate into a matrix\n",
    "data = np.transpose( np.concatenate((a,b,c),axis=1) )\n",
    "\n",
    "\n",
    "# plot data\n",
    "plt.plot(data[:,0],data[:,1],'ko',markerfacecolor='w')\n",
    "plt.title('Raw (preclustered) data')\n",
    "plt.xticks([])\n",
    "plt.yticks([])\n",
    "\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "z5hEt8nigBL1"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "## initialize random cluster centroids\n",
    "k = 3 # extract three clusters\n",
    "\n",
    "# random cluster centers (randomly sampled data points)\n",
    "ridx = np.random.choice(range(len(data)),k,replace=False)\n",
    "centroids = data[ridx,:]\n",
    "\n",
    "\n",
    "# setup the figure\n",
    "fig,axs = plt.subplots(2,2,figsize=(6,6))\n",
    "axs = axs.flatten()\n",
    "lineColors = [ [0,0,0],[.4,.4,.4],[.8,.8,.8] ]#'rbm'\n",
    "\n",
    "\n",
    "# plot data with initial random cluster centroids\n",
    "axs[0].plot(data[:,0],data[:,1],'ko',markerfacecolor='w')\n",
    "axs[0].plot(centroids[:,0],centroids[:,1],'ko')\n",
    "axs[0].set_title('Iteration 0')\n",
    "axs[0].set_xticks([])\n",
    "axs[0].set_yticks([])\n",
    "\n",
    "\n",
    "\n",
    "# loop over iterations\n",
    "for iteri in range(3):\n",
    "    \n",
    "  # step 1: compute distances\n",
    "  dists = np.zeros((data.shape[0],k))\n",
    "  for ci in range(k):\n",
    "    dists[:,ci] = np.sum((data-centroids[ci,:])**2,axis=1)\n",
    "        \n",
    "  # step 2: assign to group based on minimum distance\n",
    "  groupidx = np.argmin(dists,axis=1)\n",
    "    \n",
    "  # step 3: recompute centers\n",
    "  for ki in range(k):\n",
    "    centroids[ki,:] = [ np.mean(data[groupidx==ki,0]), np.mean(data[groupidx==ki,1]) ]\n",
    "  \n",
    "\n",
    "  # plot data points\n",
    "  for i in range(len(data)):\n",
    "    axs[iteri+1].plot([ data[i,0],centroids[groupidx[i],0] ],[ data[i,1],centroids[groupidx[i],1] ],color=lineColors[groupidx[i]])\n",
    "  axs[iteri+1].plot(centroids[:,0],centroids[:,1],'ko')\n",
    "  axs[iteri+1].set_title(f'Iteration {iteri+1}')\n",
    "  axs[iteri+1].set_xticks([])\n",
    "  axs[iteri+1].set_yticks([])\n",
    "\n",
    "\n",
    "plt.savefig('Figure_04_03.png',dpi=300)\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "2nKemxTUgHMM"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    ""
   ],
   "metadata": {
    "id": "kUMVrvJ2gBPY"
   },
   "execution_count": null,
   "outputs": []
  }
 ]
}