lectures/lecture10/lectures10and11.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lectures 10 and 11: Data analysis with `pandas` and `scikit-learn`\n",
    "\n",
    "One of the strengths of the python language is the amazing collection of packages built on top of the core library. Here we'll learn about two libraries that are useful for \"data science\". \n",
    "\n",
    "`pandas` is a multi-purpose data analysis library built on top of numpy and matplotlib. It's great for analyzing tabular data.\n",
    "\n",
    "https://pandas.pydata.org/\n",
    "\n",
    "Some useful introductory `pandas` docs:\n",
    "\n",
    "https://pandas.pydata.org/pandas-docs/stable/getting_started/10min.html\n",
    "\n",
    "https://pandas.pydata.org/pandas-docs/stable/getting_started/basics.html\n",
    "\n",
    "https://pandas.pydata.org/pandas-docs/stable/user_guide/cookbook.html\n",
    "\n",
    "\n",
    "`scikit-learn` is a nice library of machine learning tools with a clean and consistent interface. \n",
    "\n",
    "http://scikit-learn.org/stable/\n",
    "\n",
    "Both pandas and scikit-learn are built on top of the core python numeric library `numpy`, so we'll start with a quick intro to `numpy`.\n",
    "\n",
    "https://numpy.org\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Warmup: indexing into lists"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# how do we get the \"third\" element of a_list, i.e. the 9?\n",
    "a_list = [2, 7, 9]\n",
    "\n",
    "a_list[2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# lists can contain other lists. Here's a list made up of two lists:\n",
    "b_list = [[0, 11, 1], [3, 6, 2]]\n",
    "\n",
    "# how do we get the \"first\" list, ie [0, 11, 1] ?\n",
    "b_list[0]\n",
    "\n",
    "# how do we get the \"second\" element of the \"first\" list, i.e. the 11?\n",
    "b_list[0][1]\n",
    "\n",
    "# python expressions can be composed of a complex sequence of operations:\n",
    "import math\n",
    "math.log((b_list[0][1])**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Quick intro to numpy\n",
    "`numpy` is the foundation on which most of python data analysis is built. It provides highly optimized routines for operating on multidimensional numeric data. The only things we really need to understand at this point are:\n",
    "* there are things called `numpy` arrays (class `ndarray`)\n",
    "* we index into arrays like we index into lists, but with multiple indices if there are multiple dimensions\n",
    "* `numpy` has all sorts of fast numeric operations, so if we are dealing with lots of numeric data we are usually better off getting it into a numpy array than writing `for` loops ourselves. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# this is how we import numpy\n",
    "import numpy as np\n",
    "\n",
    "# do this so that our plots show up in the notebook\n",
    "# may not be necessary in all environments\n",
    "#\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we can create a numpy array from a list\n",
    "# here is a 1-dimensional array of integers\n",
    "A1 = np.array([0,4,5,6])\n",
    "print(A1)\n",
    "print('A1 has data type:', A1.dtype)\n",
    "print('A1 has shape:', A1.shape)\n",
    "\n",
    "print(A1[2]) # indexing into a 1-D array is like indexing into a list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we can also create a 2D array, which we can think of as a list of lists.\n",
    "# Here we create a 2D array from a list of lists\n",
    "A2_list = [ [0,3,4], [2,7,1.5] ]\n",
    "print('A2_list:\\n', A2_list)\n",
    "\n",
    "A2 = np.array(A2_list)\n",
    "print('A2:\\n', A2) # All the numbers were turned into floats\n",
    "print('A2 has data type:', A2.dtype) # arrays have a single consistent datatype\n",
    "print('A2 has shape:', A2.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We can think of A2 as being a list of the individual rows\n",
    "# so this gives us the first row:\n",
    "print('first row:', A2[0])\n",
    "\n",
    "# and this gives us the second element of the first row\n",
    "print('second element of first row:', A2[0][1])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# To avoid typing all those [][], numpy gives us a shorter way\n",
    "# of indexing into multidimensional arrays:\n",
    "print('The element in row 0, column 1 is',A2[0][1]) # old way, as a \"list of lists\"\n",
    "print('The element in row 0, column 1 is',A2[0,1]) # new way, with a comma\n",
    "print('The element in row 1, column 2 is',A2[1,2])\n",
    "\n",
    "# we can take entire rows or columns by using : instead of a number\n",
    "print('Row 0 of A2 is', A2[0,:])\n",
    "print('Column 1 of A2 is', A2[:,1])\n",
    "\n",
    "# we don't have to take the entire row:\n",
    "print('The first two elements of row 1 of A1 are', A2[0,:2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# np.zeros is handy for creating and initializing an array:\n",
    "nrows = 4\n",
    "ncols = 5\n",
    "A = np.zeros((nrows, ncols))\n",
    "print(A)\n",
    "\n",
    "for i in range(nrows):\n",
    "    for j in range(ncols):\n",
    "        A[i,j] = i*j\n",
    "print(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# matplotlib is the core plotting library for python\n",
    "# pandas plotting and other fancy libraries like seaborn are built on top of matplotlib\n",
    "#\n",
    "import matplotlib.pyplot as plt # this is how we get easy access to all of matplotlib's plotting commands\n",
    "plt.imshow(A) # show the array A as a colored image\n",
    "plt.colorbar() # show the color scale"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# numpy has lots of builtin functions that can be applied to arrays\n",
    "\n",
    "print('the max value in A is',np.max(A))\n",
    "print('the min value in A is',np.min(A)) \n",
    "print('the mean and standard deviation of the values in A are',np.mean(A), np.std(A)) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# here we use a builtin numpy routine to take the sine of the entire array A\n",
    "B = np.sin(A)\n",
    "plt.imshow(B)\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Practice time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create a 3 (rows) by 5 (columns) numpy array filled with random numbers between 0 and 1\n",
    "# use the random() method below to get the random numbers\n",
    "#\n",
    "# find the mean and standard deviation of the array values\n",
    "# find the mean of the 1st row (ie row index 0)\n",
    "# \n",
    "import random\n",
    "num = random.random()\n",
    "print('A single random number between 0 and 1, uniformly distributed:', num)\n",
    "\n",
    "# now create and fill the whole array\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# this is the lazy answer \n",
    "#np.random.rand(3,5)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# figure out what each of these four statements is doing. You could insert print() statements to \n",
    "# see what A looks like after each.\n",
    "\n",
    "A = np.zeros((10,2))\n",
    "\n",
    "A[:,0] = range(10)\n",
    "\n",
    "A[:,1] = np.cos(A[:,0])\n",
    "\n",
    "plt.plot(A[:,0], A[:,1])\n",
    "         \n",
    "# plot is very configurable\n",
    "#plt.plot(A[:,0], A[:,1], c='r', marker='o', alpha=0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# On to `pandas`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# this is how we import pandas\n",
    "#\n",
    "import pandas as pd\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will typically get our data by reading from files. `pandas` has a variety of file-reading functions with names like `read_csv` and `read_excel` . Find them all by typing `pd.read_[TAB]`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read in Fisher's famous iris dataset, which contains anatomical data on a variety of flowers\n",
    "# read_csv returns a pandas dataframe\n",
    "iris = pd.read_csv('data/iris_dataset.csv')\n",
    "\n",
    "iris"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Note that this is a \"tidy\" dataset (as seen in lecture 5)\n",
    "Each row corresponds to a single sample (in this case, a single flower that was measured). The columns are different measurement types, and a species identifier. The machine learning work flows we will learn about expect data arranged in this form. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "type(iris)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pandas has two main datatypes. A `DataFrame` like the `iris` object represents multiple columns of data (variables) for multiple rows of observations. A `Series` represents a single column or row of data.\n",
    "\n",
    "A DataFrame object has an `index` which represents and provides named access to the rows, and a `columns` attribute which represents the columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This index was created automatically when the data was read in. It's just a counter.\n",
    "# Some datasets might have a natural index such as a date column, an observation identifier, etc.\n",
    "iris.index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# these are the columns of the iris dataset\n",
    "iris.columns\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# here we get a single column, using dictionary-style indexing\n",
    "# we could also just write iris.sepal_length\n",
    "#\n",
    "iris['sepal_length']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# a single column is of type Series\n",
    "type(iris['sepal_length'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# super-useful function for getting summary statistics on a dataframe\n",
    "iris.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `describe` function just returns another `pandas DataFrame`, so we can use it for downstream analysis (more details on plotting below)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# make a new dataframe object with the results of iris.describe()\n",
    "df = iris.describe()\n",
    "print(df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get rid of the \"count\" row because it screws up the y-scale on the bar plot\n",
    "# note that df.drop returns a copy, it does not modify df in-place\n",
    "df = df.drop('count')\n",
    "print(df)\n",
    "\n",
    "# now make a bar plot. Lots more on plotting later\n",
    "df.plot.bar()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Digression: taking subsets of a dataframe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# the == expression below creates a column (pd.Series) of True's and False's, with the value in a given\n",
    "# row being True if the species for that row is 'setosa' and False otherwise\n",
    "iris['species'] == 'setosa'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we can create a new dataframe by indexing into an existing dataframe with this column of True's and False's\n",
    "iris_subset = iris[ iris['species'] == 'setosa' ]\n",
    "\n",
    "# here are the stats for just the setosa species flowers\n",
    "iris_subset.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we can loop over all the unique elements in the 'species' column and \n",
    "#  get stats for each subset of the dataframe\n",
    "#\n",
    "# here we use a python \"set\" object to remove redundant items (reducing from 150 strings to just 3)\n",
    "# Background: python sets are special-purpose containers that hold unordered sets of elements. \n",
    "# It is much faster to check if an element is in a set than to check if it's in a list. Creating\n",
    "# a set from a list or other kind of sequence data removes all the duplicates\n",
    "\n",
    "species_set = set(iris['species'])\n",
    "print('species_set=', species_set)\n",
    "\n",
    "for species in species_set:\n",
    "    print('stats for species', species)\n",
    "    print(iris[ iris['species'] == species].describe())\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Digression `groupby`\n",
    "Another handy dataframe method is `groupby` which collects rows based on their value for a given column and returns an object that can be used to access summary stats about those subsets. Here we are using a groupby operation to look at the rows grouped by species."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get a table (dataframe) of the mean column values for the different species \n",
    "iris.groupby('species').mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get full stats on the species subsets using the describe() function\n",
    "iris.groupby('species').describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Saving tabular data and images"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "DataFrame has `to_csv` and other `to_*` functions (see them using `df.to_<TAB>`) for saving the contents to a file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We can save a dataframe to a file using the to_csv function:\n",
    "df = iris.groupby('species').mean()\n",
    "df.to_csv('species_mean_values.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`plt.savefig` saves the most recently plotted figure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# for non-numeric columns we can count the number of times each value appears using value_counts():\n",
    "df = iris['species'].value_counts()\n",
    "print(df)\n",
    "\n",
    "# this is kind of silly\n",
    "df.plot.bar()\n",
    "\n",
    "# this saves the silly bar plot to a pngfile named silly_bar.png\n",
    "plt.savefig('silly_bar.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# plotting with pandas\n",
    "Pandas has all sorts of built-in plotting routines. If `dataframe` is a pandas `DataFrame`, `dataframe.plot.TAB` will tab complete to a list of plotting functions.\n",
    "\n",
    "Pandas plotting is built on top of `matplotlib`, so plotting functions typically return `matplotlib` objects and `matplotlib.pyplot` commands can be used to tweak features of the plots just as when using plain-vanilla matplotlib. \n",
    "\n",
    "See the pandas visualization page:\n",
    "https://pandas.pydata.org/pandas-docs/stable/visualization.html\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The dataframe object has a couple plotting methods: `plot` and `hist` and maybe one or two others"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "iris.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "iris.hist()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Many more plotting methods are available in `dataframe.plot.*`, which you can browse using tab completion "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#iris.plot.<TAB>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "iris.plot.hist(figsize=(8,6),alpha=0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Maybe a histogram isn't the best visualization, how about kernel density estimation, which produces smoothed density plots?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "iris.plot.density(figsize=(8,6))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Scatter plots of one column against another are easy enough:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "iris.plot.scatter( x='sepal_length', y='sepal_width', figsize=(8,6) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Practice time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a hexbin plot of petal_length versus petal_width\n",
    "# Try out different gridsizes\n",
    "# Compare the output to a scatterplot of petal_length versus petal_width\n",
    "# hexbin is one of the builtin plotting functions: iris.plot.hexbin?\n",
    "#\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# What does this plot represent? \n",
    "# You could first try typing iris.mean() to see what that dataframe looks like...\n",
    "\n",
    "iris.mean().plot.pie()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### But what about the `species` structure in this dataset?\n",
    "\n",
    "Let's add some color. This is super-easy to do using the plotting routines in the `seaborn` package, but for right now we'll stick with `pandas` and `matplotlib`.\n",
    "\n",
    "Definitely check out the [seaborn gallery](https://seaborn.pydata.org/examples/index.html), though.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Here we are going to set up a dictionary that maps from iris species to plotting colors\n",
    "color_dict = {'virginica': 'blue', 'versicolor': 'orange', 'setosa': 'green'}\n",
    "\n",
    "# now let's use this dictionary to make a new column for the iris dataframe, containing the colors\n",
    "colors = []\n",
    "for species in iris['species']:\n",
    "    colors.append(color_dict[species])\n",
    "\n",
    "print(colors[:5])\n",
    "\n",
    "# or colors = [color_dict[s] for s in iris['species']] # list comprehension\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is how we add a column to a dataframe:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "iris['color'] = colors \n",
    "\n",
    "#or iris['color'] = iris['species'].map(color_dict)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# confirm the new column:\n",
    "iris.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now we can see the different flower species\n",
    "#\n",
    "iris.plot.scatter( x='sepal_length', y='sepal_width', c='color')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# this is a nice way of comparing the values in different columns of a DataFrame\n",
    "pd.plotting.scatter_matrix(iris, diagonal='kde', figsize=(12,12), c=iris['color']);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Digression: list comprehensions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# quick introduction to list comprehensions:\n",
    "evens = [ 2*x for x in range(10)]\n",
    "print('evens one way:', evens)\n",
    "\n",
    "evens = [ x for x in range(20) if x%2 == 0]\n",
    "print('evens another way:', evens)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Digression: seaborn\n",
    "As mentioned above, seaborn is a nice high-level plotting library that is built on top of matplotlib. \n",
    "Seaborn understands pandas dataframes and tidy data and can make generating complex\n",
    "visualizations really easy.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import seaborn as sns\n",
    "\n",
    "# it identifies the numeric columns, picks a nice color palette, adds nice density plots\n",
    "# along the diagonal and a nice legend\n",
    "sns.pairplot(iris, hue='species')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Practice time\n",
    "Read the data in the csv file `data/vdj.csv` into a pandas dataframe using `pd.read_csv` like we did for the iris dataset. How many columns are there? How many numeric columns? Make a scatter plot of 'reads' versus 'umis'. Count how many times each of the different gene names occurs in the `v_gene` column using the `.value_counts()` method. Can you make a bar plot or pie chart of those numbers (like we did for the species numbers above?)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# start by reading the data into a dataframe. \n",
    "#df = pd.read_csv(...)\n",
    "\n",
    "# .head() is a nice way to get a peak at a dataset\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Machine learning with the `scikit-learn` package\n",
    "\n",
    "borrowing heavily from \"The Python Data Science Handbook\" by Jake VanderPlas\n",
    "\n",
    "https://jakevdp.github.io/PythonDataScienceHandbook/\n",
    "\n",
    "## categories of machine learning\n",
    "Broadly speaking, machine learning for data science can be broken down into **supervised** and **unsupervised** methods. \n",
    "\n",
    "In **supervised** machine learning, we have a set of labeled training data, and we'd like to train a classifier or predictor in order to make predictions on some new set of testing data. We may be trying to predict something discrete, like different categories of data (for example, biological species from flower dimensions in the iris dataset; this is called *classification*), or predicting something continuous-valued like a response variable as a function of predictor variable(s) (this is called *regression*). \n",
    "\n",
    "In **unsupervised** machine learning, we are trying to discover hidden structure in a dataset. Two broad approaches are *clustering*, where we look for subgroups of similar data points, and *dimensionality reduction*, where we project a high-dimensional dataset into a lower-dimensional space while trying to preserve important features of the data points (e.g., keeping points that were nearby in the high-dimensional space also nearby in the low-dimensional projection)  \n",
    "\n",
    "`scikit-learn` has a nice flow-chart with even more information:\n",
    "\n",
    "https://scikit-learn.org/stable/tutorial/machine_learning_map/index.html\n",
    "\n",
    "![logo](https://scikit-learn.org/stable/_static/ml_map.png)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Unsupervised dimensionality reduction on the iris dataset with PCA\n",
    "\n",
    "`scikit-learn` has several methods for dimensionality reduction. One popular method is principal components analysis (PCA). \n",
    "\n",
    "Quoting from [the wikipedia entry for PCA](https://en.wikipedia.org/wiki/Principal_component_analysis), \"Principal component analysis (PCA) is a statistical procedure that uses an orthogonal transformation to convert a set of observations of possibly correlated variables (entities each of which takes on various numerical values) into a set of values of linearly uncorrelated variables called principal components. This transformation is defined in such a way that the first principal component has the largest possible variance (that is, accounts for as much of the variability in the data as possible), and each succeeding component in turn has the highest variance possible under the constraint that it is orthogonal to the preceding components\"\n",
    "\n",
    "This image shows a 2-D dataset with vectors indicating the directions of the first two principal components:\n",
    "\n",
    "![logo](https://upload.wikimedia.org/wikipedia/commons/f/f5/GaussianScatterPCA.svg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To apply PCA to the iris dataset, we need to create an array with just the numeric features we'll be using. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'iris' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m/workspaces/tfcb_2022/lectures/lecture10/lectures10and11.ipynb Cell 75\u001b[0m in \u001b[0;36m<cell line: 1>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> <a href='vscode-notebook-cell://dev-container%2B2f55736572732f616e616761696e692f4465736b746f702f746663625f323032322f746663625f32303232/workspaces/tfcb_2022/lectures/lecture10/lectures10and11.ipynb#Y134sdnNjb2RlLXJlbW90ZQ%3D%3D?line=0'>1</a>\u001b[0m iris\u001b[39m.\u001b[39mhead()\n",
      "\u001b[0;31mNameError\u001b[0m: name 'iris' is not defined"
     ]
    }
   ],
   "source": [
    "iris.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create a features array with just the measurement data\n",
    "\n",
    "X_iris = iris.drop(columns=['species','color'])\n",
    "X_iris.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_iris.head()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A standard scikit-learn workflow\n",
    "\n",
    "### Step 1. create and configure the model (here `PCA`)\n",
    "\n",
    "### Step 2. `fit` the model to the data\n",
    "\n",
    "### Step 3. `transform` the data using the model. \n",
    "For classification and clustering, where we have discrete labels, this third step is typically called `predict`, as we'll see below for K means clustering."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'X_iris' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m/workspaces/tfcb_2022/lectures/lecture10/lectures10and11.ipynb Cell 79\u001b[0m in \u001b[0;36m<cell line: 8>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      <a href='vscode-notebook-cell://dev-container%2B2f55736572732f616e616761696e692f4465736b746f702f746663625f323032322f746663625f32303232/workspaces/tfcb_2022/lectures/lecture10/lectures10and11.ipynb#Y141sdnNjb2RlLXJlbW90ZQ%3D%3D?line=4'>5</a>\u001b[0m pca \u001b[39m=\u001b[39m PCA(n_components\u001b[39m=\u001b[39m\u001b[39m4\u001b[39m)\n\u001b[1;32m      <a href='vscode-notebook-cell://dev-container%2B2f55736572732f616e616761696e692f4465736b746f702f746663625f323032322f746663625f32303232/workspaces/tfcb_2022/lectures/lecture10/lectures10and11.ipynb#Y141sdnNjb2RlLXJlbW90ZQ%3D%3D?line=6'>7</a>\u001b[0m \u001b[39m# Step 2. fit the model to the data\u001b[39;00m\n\u001b[0;32m----> <a href='vscode-notebook-cell://dev-container%2B2f55736572732f616e616761696e692f4465736b746f702f746663625f323032322f746663625f32303232/workspaces/tfcb_2022/lectures/lecture10/lectures10and11.ipynb#Y141sdnNjb2RlLXJlbW90ZQ%3D%3D?line=7'>8</a>\u001b[0m pca\u001b[39m.\u001b[39mfit(X_iris)\n",
      "\u001b[0;31mNameError\u001b[0m: name 'X_iris' is not defined"
     ]
    }
   ],
   "source": [
    "from sklearn.decomposition import PCA # sklearn == scikit-learn package name\n",
    "\n",
    "# This is the general scikit-learn workflow\n",
    "# Step 1. create and configure the model:\n",
    "pca = PCA(n_components=4)\n",
    "\n",
    "# Step 2. fit the model to the data\n",
    "pca.fit(X_iris)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here are the four principal component vectors. They can be thought as vectors of weights, one weight for each of the columns in the original feature matrix, with the magnitude of each weight reflecting the contribution of that feature (e.g., petal_width) to the component. The first component gives the direction in feature space along which the greatest variation exists in the dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(pca.components_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The principal components are all unit (length 1) vectors and they are \n",
    "#  orthogonal (perpendicular) to one another. Here we test that a bit.\n",
    "#\n",
    "\n",
    "print('The length of the first principal component is',np.linalg.norm(pca.components_[0]))\n",
    "print('The dot product of the first and second principal components is',\n",
    "    np.dot( pca.components_[0], pca.components_[1] ))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('The fraction of the total variance explained by each component is:', pca.explained_variance_ratio_)\n",
    "print('The total explained variance is :', np.sum(pca.explained_variance_ratio_))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 3. Now we can `transform` the data using the fitted model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 3. Apply the model to 'transform' the data, in this case by mapping\n",
    "#    the input feature matrix onto the principal components\n",
    "#\n",
    "\n",
    "X_pca = pca.transform(X_iris)\n",
    "X_pca.shape\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now the first two columns of X_pca represent the axes of \n",
    "#  greatest variation, so we can plot just those and retain\n",
    "#  the maximum information.\n",
    "\n",
    "plt.scatter( X_pca[:,0], X_pca[:,1], c=colors )\n",
    "plt.title('First two principal components')\n",
    "plt.xlabel('PCA1')\n",
    "plt.ylabel('PCA2'); # add ; so notebook doesn't print the Text object returned by plt.ylabel\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Digression on rescaling the data\n",
    "\n",
    "This seems to work pretty well, but in some situations the multidimensional variables we are measuring are not directly comparable (for example because they have different units), in which case PCA will be dominated by the features that have the largest raw standard deviation. We can correct for this by first normalizing all columns to have mean 0 and standard deviation 1, using one of `scikit-learn`'s many preprocessing routines: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# compare the standard deviations of the different features. \n",
    "#  They aren't too far apart...\n",
    "#\n",
    "X_iris.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import StandardScaler\n",
    "\n",
    "# Step 1: create the sklearn object\n",
    "scaler = StandardScaler()\n",
    "scaler"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 2: fit the model (in this case, that means estimating the mean and variance)\n",
    "scaler.fit( X_iris )\n",
    "print('mean:', scaler.mean_)\n",
    "print('variance:', scaler.var_)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 3: apply the fitted model to transform the data (in this case,\n",
    "#  that means subtract the mean and divide by the standard deviation)\n",
    "X_iris_scaled = scaler.transform( X_iris )\n",
    "\n",
    "# check the mean, std-dev of the first column\n",
    "print( np.mean( X_iris_scaled[:,0]), np.std( X_iris_scaled[:,0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 1: create PCA model\n",
    "pca2 = PCA(n_components=4)\n",
    "\n",
    "# Step 2: solve for PCs\n",
    "pca2.fit(X_iris_scaled)\n",
    "\n",
    "# Step 3: transform the data to the new PC coordinate system\n",
    "X_pca2 = pca2.transform( X_iris_scaled )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set a larger figure size\n",
    "plt.figure(figsize=(14,8))\n",
    "\n",
    "# create a multi-panel plot:\n",
    "nrows=1 # one row,\n",
    "ncols=2 #  two columns\n",
    "\n",
    "plt.subplot(nrows, ncols, 1)\n",
    "plt.scatter( X_pca[:,0], X_pca[:,1], c=colors)\n",
    "plt.title('raw')\n",
    "plt.xlabel('PCA1', fontsize=6)\n",
    "plt.ylabel('PCA2')\n",
    "\n",
    "plt.subplot(nrows, ncols, 2)\n",
    "plt.scatter( X_pca2[:,0], X_pca2[:,1], c=colors)\n",
    "plt.title('scaled')\n",
    "plt.xlabel('PCA1')\n",
    "\n",
    "\n",
    "print('old explained variance:',pca.explained_variance_ratio_)\n",
    "print('old components:',pca.components_)\n",
    "print('new explained variance:',pca2.explained_variance_ratio_)\n",
    "print('new components:',pca2.components_)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Unsupervised clustering on the iris dataset using K-means clustering\n",
    "\n",
    "K-means is an iterative clustering approach in which we repeatedly (1) choose K cluster centroids, (2) assign each data point to the closest centroid, and (3) update the centroids to reflect the means of the new clusters. The starting cluster centroids are chosen randomly (hence the use of `random_state=10` below to get reproducible results). Here's a snapshot from the [wikipedia page for K-means clustering]( https://en.wikipedia.org/wiki/K-means_clustering).\n",
    "\n",
    "![logo](kmeans.png)\n",
    "\n",
    "The naming convention I'm following here is to call the input (2-D) data `X` and the (1-D) classification/clusters `y`. This loosely follows the math-y convention where `x` is the independent variable and `y` is the dependent variable: `y=f(x)` \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.cluster import KMeans\n",
    "\n",
    "# STEP 1: create and configure the model\n",
    "# create the classifier object, tell it we are looking for 3 clusters\n",
    "#  and set the random_state for reproducibility\n",
    "kmeans = KMeans(n_clusters=3, random_state=10)\n",
    "\n",
    "# STEP 2: fit the model on the data, using \"model.fit\"\n",
    "# learn the cluster assignments\n",
    "kmeans.fit(X_iris)\n",
    "\n",
    "# STEP 3: apply the model to the data, using \"model.predict\"\n",
    "# create a new array with the cluster assignments, represented by the integers 0, 1, and 2.\n",
    "y_kmeans = kmeans.predict(X_iris)\n",
    "y_kmeans"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# compare the species/colors array to the kmeans cluster assignments\n",
    "np.array(colors) # converting to a numpy array gives a more compact view"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Here we convert the kmeans clusters into colors so we can compare them with\n",
    "#  the species colors we set up earlier. Looking at the two arrays printed \n",
    "#  above, we can see a rough correspondence between kmeans cluster # and\n",
    "#  species color\n",
    "#\n",
    "#  cluster   color\n",
    "#    0         green\n",
    "#    1         orange\n",
    "#    2         blue\n",
    "\n",
    "# use a dictionary to store this mapping\n",
    "cluster_to_color = {0:'green', 1:'orange', 2:'blue'}\n",
    "\n",
    "kmeans_colors = []\n",
    "\n",
    "for cluster in y_kmeans:\n",
    "    kmeans_colors.append(cluster_to_color[cluster])\n",
    "    \n",
    "    \n",
    "error_colors = []\n",
    "for color, kmeans_color in zip( colors, kmeans_colors ):\n",
    "    if color == kmeans_color:\n",
    "        # correct cluster assignment\n",
    "        error_colors.append('b')\n",
    "    else:\n",
    "        # incorrect cluster assignment\n",
    "        error_colors.append('r')\n",
    "print('error_colors:', error_colors)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# put three plots side-by-side\n",
    "plt.figure(figsize=(12,4))\n",
    "plt.subplot(1,3,1)\n",
    "plt.scatter( X_pca[:,0], X_pca[:,1], c=colors)\n",
    "plt.title('species colors')\n",
    "\n",
    "plt.subplot(1,3,2)\n",
    "# a more succinct way to set up kmeans_colors using a list comprehension\n",
    "#kmeans_colors = [ 'bgr'[x] for x in y_kmeans]\n",
    "plt.scatter( X_pca[:,0], X_pca[:,1], c=kmeans_colors)\n",
    "plt.title('kmeans cluster colors')\n",
    "\n",
    "plt.subplot(1,3,3)\n",
    "# a more succinct way to set up error_colors using a list comprehension, zip, and conditional assignment\n",
    "#error_colors = [ ('red' if x!=y else 'blue') for x,y in zip( colors, kmeans_colors) ]\n",
    "plt.scatter( X_pca[:,0], X_pca[:,1], c=error_colors, alpha=0.5)\n",
    "plt.title('species/cluster agreement\\nblue=match, red=mismatch'); # semicolon so it doesn't print out Text"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Supervised classification on the iris dataset using Gaussian naive bayes classifier. \n",
    "\n",
    "A very simple approach to classification is to assume that the feature values for each class are drawn independently of one another from one dimensional gaussian distributions (bell-shaped curves), and to choose the predicted class labels that maximize the likelihood of the observed feature values. The Gaussian naive bayes classifier implements this method. Fitting the model involves estimating, for each feature (like sepal width) and each class (like 'setosa'), the mean and standard deviation of the feature distribution over the class members.\n",
    "\n",
    "For further details check out the [wikipedia page](https://en.wikipedia.org/wiki/Naive_Bayes_classifier#Gaussian_naive_Bayes)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## we need an array with the categories we are trying to classify:\n",
    "y_iris = iris['species']\n",
    "print(y_iris.head())\n",
    "y_iris.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## choose the model\n",
    "from sklearn.naive_bayes import GaussianNB\n",
    "\n",
    "## create (and optionally configure) the model\n",
    "model = GaussianNB()\n",
    "\n",
    "## train the model\n",
    "model.fit(X_iris, y_iris)\n",
    "\n",
    "## see how well we do:\n",
    "y_model = model.predict(X_iris)\n",
    "y_model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can see how well we did:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import accuracy_score\n",
    "print('The fraction of correctly classified samples is',accuracy_score(y_iris, y_model))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# or we can do this \"by hand\"\n",
    "num_correct = 0\n",
    "\n",
    "# note that the python \"zip\" function works here even though y_model is a numpy array \n",
    "#  and y_iris is a pandas Series \n",
    "for prediction, correct_answer in zip(y_model, y_iris):\n",
    "    if prediction == correct_answer:\n",
    "        num_correct += 1\n",
    "        \n",
    "print(num_correct/len(y_iris))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### What's the problem with the above analysis?\n",
    "We trained and tested on the same set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## redo, correctly\n",
    "\n",
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "Xtrain, Xtest, ytrain, ytest = train_test_split(X_iris, y_iris,\n",
    "                                                random_state=30)\n",
    "\n",
    "# create the model\n",
    "model2 = GaussianNB()\n",
    "\n",
    "# train the model\n",
    "model2.fit(Xtrain, ytrain)\n",
    "\n",
    "# predict on NEW DATA that wasn't used for training\n",
    "y_model = model2.predict(Xtest)\n",
    "\n",
    "# assess accuracy\n",
    "accuracy_score(ytest, y_model)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear regression using scikit-learn\n",
    "\n",
    "In linear regression we model a continuous-valued response variable as a linear combination of predictor variables. The parameters we need to fit are the coefficients for each predictor variable (the *slope* in the 1-D case), and if desired a constant offset (the *intercept*)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's start by making some slightly noisy data using a linear function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "slope = 2.\n",
    "intercept = -1\n",
    "\n",
    "noise_factor = 0.1\n",
    "\n",
    "num_points = 50\n",
    "\n",
    "# create uniformly distributed x-values between 0 and 1\n",
    "x = np.random.rand( num_points)\n",
    "print(x)\n",
    "\n",
    "# now create y values that are a linear function of the x-values, with some random noise added in\n",
    "# see how easy it is to multiply and add entire arrays of data with numpy\n",
    "# \n",
    "y = slope * x - intercept + noise_factor * np.random.randn(50)\n",
    "\n",
    "plt.scatter(x, y);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LinearRegression\n",
    "\n",
    "\n",
    "model = LinearRegression()\n",
    "model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A slightly tricky point is that scikit-learn routines generally expect a two-dimensional data array with one row for each observation (datapoint) and columns for each of the data values. Our set of x-values is actually a one-dimensional array (basically a list). So we need to change it's shape to make it two-dimensional, using the \"reshape\" command: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"old 1D shape:\",x.shape)\n",
    "X = x.reshape( (num_points,1)) # note: upper-case X\n",
    "print(\"new 2D shape:\", X.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# now we can fit the model to the data:\n",
    "model.fit(X, y)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# this is the slope value we recovered from the model: not too far from the actual value!\n",
    "model.coef_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# and this is the intercept:\n",
    "model.intercept_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can plot the fitted line:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create a linear array of new test points\n",
    "xfit = np.linspace(0,1,10)\n",
    "\n",
    "# reshape into a 10-rows by 1-column 2D array to make scikit-learn happy\n",
    "Xfit = xfit.reshape((10,1))\n",
    "\n",
    "# predict the y-values using our fitted model:\n",
    "yfit = model.predict(Xfit)\n",
    "\n",
    "# scatter plot of original data \n",
    "plt.scatter(x, y)\n",
    "\n",
    "# line plot of fitted line\n",
    "plt.plot(xfit, yfit, c='r');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# Random junk below here, probably unused"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# another useful way of creating an array is to use the numpy random module\n",
    "# here we create an array of normally distributed values\n",
    "nrows = 10\n",
    "ncols = 10\n",
    "A = np.random.randn(nrows,ncols)\n",
    "print(A)\n",
    "plt.imshow(A)\n",
    "plt.colorbar()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math # import the math module\n",
    "nrows = 1000\n",
    "ncols = 1000\n",
    "A = np.random.randn(nrows,ncols)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%%timeit\n",
    "B = np.sin(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit\n",
    "B = np.zeros((nrows,ncols))\n",
    "for i in range(nrows):\n",
    "    for j in range(ncols):\n",
    "        B[i,j] = math.sin( A[i,j])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# do this so the plots appear inline in the notebook:\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "iris.sort_index(axis=1).head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "iris.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "iris.sort_values('petal_width',axis=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sorting: note that the default in most pandas operations is to return *copies* not to change the existing dataframe. This can usually be changed by passing inplace=True as one of the arguments. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "iris.groupby('species').mean()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#X_normed = X_iris.transform(lambda x: (x - x.mean()) / x.std())\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.9.7 ('tfcb2022_rna')",
   "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.9.7"
  },
  "vscode": {
   "interpreter": {
    "hash": "3011704b6779c0b541a689cf71206a950289813906dd73e2862ef7e407814574"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}