{ "nbformat": 4, "nbformat_minor": 0, "metadata": { "colab": { "provenance": [], "collapsed_sections": [], "toc_visible": true, "authorship_tag": "ABX9TyOq1WKMwitCV31H9X5WPjPT", "include_colab_link": true }, "kernelspec": { "name": "python3", "display_name": "Python 3" } }, "cells": [ { "cell_type": "markdown", "metadata": { "id": "view-in-github", "colab_type": "text" }, "source": [ "\"Open" ] }, { "cell_type": "markdown", "metadata": { "id": "S9l0lyqLvkW7" }, "source": [ "# Tutorial 1\n", "\n", "**Probability & Statistics III: Statistical Encoding & Decoding**\n", "\n", "\n", "**[insert your name]**\n", "\n", "**Important reminders**: Before starting, click \"File -> Save a copy in Drive\". Produce a pdf for submission by \"File -> Print\" and then choose \"Save to PDF\".\n", "\n", "To complete this tutorial, you should have watched Video 8.1 and 8.2.\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " Imports\n" ] }, { "cell_type": "code", "metadata": { "id": "cv9HSBNPyLV9", "cellView": "form", "tags": [ "hide-input" ] }, "source": [ "# @markdown Imports\n", "\n", "# Imports\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import ipywidgets as widgets # interactive display\n", "import math" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ " Plotting functions\n" ] }, { "cell_type": "code", "metadata": { "id": "ZIdPVYl9TzmK", "cellView": "form", "tags": [ "hide-input" ] }, "source": [ "# @markdown Plotting functions\n", "import numpy\n", "from numpy.linalg import inv, eig\n", "from math import ceil\n", "from matplotlib import pyplot, ticker, get_backend, rc\n", "from mpl_toolkits.mplot3d import Axes3D\n", "from itertools import cycle\n", "\n", "\n", "%config InlineBackend.figure_format = 'retina'\n", "plt.style.use(\"https://raw.githubusercontent.com/NeuromatchAcademy/course-content/master/nma.mplstyle\")\n" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ " Helper functions\n" ] }, { "cell_type": "code", "metadata": { "id": "0TCUlgD2L2y7", "cellView": "form", "tags": [ "hide-input" ] }, "source": [ "# @markdown Helper functions\n", "def twoD_Gaussian(xdata_tuple, amplitude, xo, yo, sigma_x, sigma_y, theta, offset):\n", " \"\"\"Create 2D Gaussian based on parameters\n", "\n", " Args:\n", " xdata_tuple (ndarray): grid of x and y values to compute Gaussian for\n", " amplitude (scalar): amplitude of Gaussian\n", " xo (scalar): center of Gaussian in x coordinates\n", " yo (scalar): center of Gaussian in y coordinates\n", " sigma_x (scalar): standard deviation of Gaussian in x direction\n", " sigma_y (scalar): standard deviation of Gaussian in y direction\n", " theta (scalar): rotation angle of Gaussian\n", " offset (scalar): offset of all Gaussian values\n", "\n", " Returns:\n", " ndarray: Gaussian values at every x/y point\n", "\n", " \"\"\"\n", " (x, y) = xdata_tuple\n", " xo = float(xo)\n", " yo = float(yo) \n", " a = (np.cos(theta)**2)/(2*sigma_x**2) + (np.sin(theta)**2)/(2*sigma_y**2)\n", " b = -(np.sin(2*theta))/(4*sigma_x**2) + (np.sin(2*theta))/(4*sigma_y**2)\n", " c = (np.sin(theta)**2)/(2*sigma_x**2) + (np.cos(theta)**2)/(2*sigma_y**2)\n", " g = offset + amplitude*np.exp( - (a*((x-xo)**2) + 2*b*(x-xo)*(y-yo)+c*((y-yo)**2)))\n", " return g.ravel()\n" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "PcrR0diWJj76" }, "source": [ "# The data\n", "\n", "In this tutorial, we will be working with simulated neural data from a visual neuron in response to Gaussian white noise and MNIST images (we'll call these natural scenes for ease even though they aren't very natural). We will be fitting LNP models for both types of stimuli separately. We have 10000 images for each type of stimuli and each image is 10 x 10 pixels. \n", "\n", "The next cell gives you `WN_images` and `NS_images`, the white noise and MNIST images respectively. Each is 10000 x 100 so the images have already been vectorized." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " Execute this cell to get and visualize example image of each type (may take a few min to download)\n" ] }, { "cell_type": "code", "metadata": { "id": "UqWBI2gK0Jsu", "cellView": "form", "tags": [ "hide-input" ] }, "source": [ "# @markdown Execute this cell to get and visualize example image of each type (may take a few min to download)\n", "\n", "np.random.seed(123)\n", "n_images = 10000\n", "\n", "# Get WN images\n", "WN_images = np.random.randn(n_images, 10*10)\n", "\n", "# Get NS images\n", "from sklearn.datasets import fetch_openml\n", "mnist = fetch_openml(name='mnist_784')\n", "mnist_images = np.array(mnist.data)\n", "mnist_images = mnist_images/255 \n", "mnist_images = mnist_images - np.mean(mnist_images)\n", "\n", "mnist_images = mnist_images.reshape((-1, 28, 28))[:, 4:24, 4:24]\n", "mnist_images = mnist_images[:, ::2, ::2]\n", "NS_images = mnist_images[:n_images].reshape((-1, 10*10))\n", "\n", "fig, axes = plt.subplots(1, 2, figsize=(10, 5))\n", "axes[0].imshow(WN_images[0].reshape((10, 10)), vmin=-1, vmax=1, cmap='gray')\n", "axes[1].imshow(NS_images[0].reshape((10, 10)), vmin=-1, vmax=1, cmap='gray')\n", "axes[0].axis('Off')\n", "axes[1].axis('Off')\n", "axes[0].set(title='WN example image')\n", "axes[1].set(title='NS example image');\n" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "qi7I3zdCKo9B" }, "source": [ "The response to each image is a summed spike count response so we do not have to worry about accounting for the time lags of the stimuli etc. We are simulating the neuron as an LNP model with an exponential nonlinearity so a linear filter, then exponential, then Poisson draws. Note that this means our LNP fits will be really good because we are using the correct model (this will literally never happen in real life...)\n", "\n", "Execute the next cell to simulate our neuron and get `WN_spike_counts` and `NS_spike_counts`. `filter` is the true linear filter of this neuron." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " Execute to simulate neural responses\n" ] }, { "cell_type": "code", "metadata": { "id": "0MM_kuaEyoAR", "cellView": "form", "tags": [ "hide-input" ] }, "source": [ "# @markdown Execute to simulate neural responses\n", "np.random.seed(0)\n", "\n", "x = np.arange(-5, 5, 1)\n", "y = np.arange(-5, 5, 1)\n", "x, y = np.meshgrid(x, y)\n", "\n", "sd_x = 1.4\n", "sd_y = .5\n", "gauss = twoD_Gaussian((x,y), 1, 0, 0, sd_x, sd_y, 35, 0)\n", "filter = gauss.reshape((10, 10))\n", "\n", "WN_lambda = np.exp(np.dot(WN_images, filter.reshape((-1,))))\n", "WN_spike_counts = np.random.poisson(WN_lambda)\n", "\n", "NS_lambda = np.exp(np.dot(NS_images, filter.reshape((-1,))))\n", "NS_spike_counts = np.random.poisson(NS_lambda)\n", "\n", "fig, axes = plt.subplots(1, 2, figsize=(10, 5), sharex=True)\n", "axes[0].plot(WN_spike_counts[0:100], 'ok')\n", "axes[1].plot(NS_spike_counts[0:100], 'ok')\n", "axes[0].set(ylabel='Spike counts', xlabel='Image number', title='WN')\n", "axes[1].set(xlabel='Image number', title='NS');\n" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "HAqXLN6LMhVa" }, "source": [ "# Exercise 1: Computing an STA\n", "\n", "We want to fit an LNP model for each type of stimulus. Since our white noise is stochastic and spherically distributed, we know we can compute a spike triggered average and it will be an unbiased estimator for our linear filter. In fact, we will assume an exponential nonlinearity so it will be the maximum likelihood estimator for our linear filter.\n", "\n", "Fill out the code below to create a function that computes the STA from a set of images and associated spike counts. Compute this STA for both white noise and natural scenes. Run the next cell to visualize your computed STAs next to the original (true) linear filter." ] }, { "cell_type": "markdown", "metadata": { "id": "ZiQHKd-lPWZj" }, "source": [ "## Answer\n", "Fill out code below" ] }, { "cell_type": "code", "source": [ "WN_spike_counts.shape" ], "metadata": { "id": "5nAQk_-S9MyW", "outputId": "c2a9fb64-c12a-4b1c-c7a1-faf22c7c9722", "colab": { "base_uri": "https://localhost:8080/", "height": 166 } }, "execution_count": null, "outputs": [ { "output_type": "error", "ename": "NameError", "evalue": "ignored", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)", "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mWN_spike_counts\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m", "\u001b[0;31mNameError\u001b[0m: name 'WN_spike_counts' is not defined" ] } ] }, { "cell_type": "code", "metadata": { "id": "2gaUc9YwM8qC" }, "source": [ "def compute_STA(images, spike_counts):\n", "\n", " STA = ...\n", "\n", " return STA\n", "\n", "\n", "WN_STA = ...\n", "NS_STA = ...\n" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ " Execute to visualize your computed STAs and the original filter\n" ] }, { "cell_type": "code", "metadata": { "id": "DM5driqtNT7Q", "cellView": "form", "tags": [ "hide-input" ] }, "source": [ "# @markdown Execute to visualize your computed STAs and the original filter\n", "\n", "fig, axes = plt.subplots(1, 3, figsize=(15, 5))\n", "\n", "axes[0].imshow(filter.reshape((10, 10)), vmin=-1, vmax=1, cmap='gray')\n", "axes[1].imshow(WN_STA.reshape((10, 10)), vmin=-1, vmax=1, cmap='gray')\n", "axes[2].imshow(NS_STA.reshape((10, 10)), vmin=-1, vmax=1, cmap='gray')\n", "\n", "for i in range(3):\n", " axes[i].axis('Off')\n", "\n", "axes[0].set(title='True filter')\n", "axes[1].set(title='White noise STA')\n", "axes[2].set(title='Natural scenes STA');\n" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "ANcFszyGNyUu" }, "source": [ "Note that the white noise STA is a pretty good estimate for the true filter, but the natural scenes STA is not!" ] }, { "cell_type": "markdown", "metadata": { "id": "LLKi7xRMN27Y" }, "source": [ "# (Optional) Exercise: Estimate the nonlinearity\n", "\n", "Estimate the nonlinearity of the LNP model (so no longer predefine it as exponential) using the method discussed in Video 8.2." ] }, { "cell_type": "code", "metadata": { "id": "Akj6UNGgN91f" }, "source": [], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "_smPEDZ9N_gx" }, "source": [ "# Exercise 2: Numerically finding the filter with natural scenes data\n", "\n", "The STA was a very convenient estimate of our linear filter of an LNP model for the white noise stimuli. Unfortunately, it is a bad estimator for the natural scenes stimuli so we will have to use a numerical approach to estimate the filter using this data. In this exercise, we will implement gradient descent ourselves." ] }, { "cell_type": "markdown", "metadata": { "id": "HB9YYg8aOYqJ" }, "source": [ "## A) Negative log likelihood equation\n", "\n", "To implement gradient descent ourselves, we will need to compute the derivative of the negative log likelihood. \n", "\n", "Write out the negative log likelihood equation for our LNP model with an exponential nonlinearity. Simplify as much as possible. Drop constants that don't depend on the filter (so we won't compute the true NLL but the relative NLL for different filters). Show the math! Use y for the spike counts and x for the images.\n", "\n", "Make the final equation clear either with the green text below or by putting it in a box or otherwise highlighting it.\n", "\n" ] }, { "cell_type": "markdown", "metadata": { "id": "zbdNpTacO26N" }, "source": [ "### **Answer**\n", "\n", "Put NLL = ... equation here (show work above or below)\n", " " ] }, { "cell_type": "markdown", "metadata": { "id": "SoOmR2s7O-7c" }, "source": [ "## B) Negative log likelihood computation\n", "\n", "Use your equation in part A to fill out the code below to compute the negative log likelihood for a given filter (k) and set of images (x) and spike counts (y). x and y would be the same shape as WN_images and WN_spike_counts for example" ] }, { "cell_type": "markdown", "metadata": { "id": "2BOfzCG1Pa86" }, "source": [ "### Answer\n", "Fill out code below" ] }, { "cell_type": "code", "metadata": { "id": "qYTZnYAuPaVV" }, "source": [ "def compute_NLL(k, x, y):\n", "\n", " NLL = ...\n", " \n", " return NLL\n", "\n", "fake_k = np.zeros((100,))\n", "NLL = compute_NLL(fake_k, WN_images, WN_spike_counts)\n", "print(NLL)" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "ElK1Y8S9Piuc" }, "source": [ "## C) Compute dNLL/dk\n", "\n", "Take your answer in part A and now take the derivative with respect to $\\bar{k}$. Note that $\\bar{k}$ is a vector so this can get tricky! I would take the derivative with respect to $\\bar{k}_o$ first (the first element of $k$). Since each entry of $\\bar{k}$ is present in the negative log likelihood equation in a similar manner, you should be able to extend your calculation for $\\frac{dNLL}{d\\bar{k}_0}$ to figure out the whole vector $\\frac{dNLL}{d\\bar{k}}$.\n", "\n", "When in confusion about dot products, my recommendation is to write out the first few elements of the dot product computation for clarifiation.\n", "\n", "\n", "Make the final equation clear either with the green text below or by putting it in a box or otherwise highlighting it. Show your work!" ] }, { "cell_type": "markdown", "metadata": { "id": "BrkEaiZoQXsz" }, "source": [ "### **Answer**\n", "\n", "Put dNLL/dk = ... equation here (show work above or below)\n", " " ] }, { "cell_type": "markdown", "metadata": { "id": "czEDnIkZQqwA" }, "source": [ "## D) Implementing gradient descent\n", "\n", "We now have all the tools we need to implement gradient descent to find an estimate of our filter k using the natural scenes data.\n", "\n", "Fill out the following code to perform gradient descent and then call it for the natural scenes data. The following cells plot the loss function (negative log likelihood) over step of the gradient descent algorithm and the fitted filter." ] }, { "cell_type": "markdown", "metadata": { "id": "-W0_xn--TY0p" }, "source": [ "### Answer\n", "Fill out code below" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " Execute to visualize negative log likelihood over gradient descent\n" ] }, { "cell_type": "code", "metadata": { "id": "GtMUWZyXO9_o", "tags": [ "hide-input" ] }, "source": [ "def gradient_descent(x, y, init_guess, n_steps = 500, alpha=10**-6):\n", "\n", " k = init_guess\n", " NLL = np.zeros((n_steps,))\n", " for i_step in range(n_steps):\n", "\n", " # Update estimate of k (assign as k)\n", " # your code here\n", "\n", " # Compute NLL at each step\n", " NLL[i_step] = compute_NLL(k, x, y)\n", "\n", " return k, NLL\n", "\n", "k, NLL = ..." ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ " Execute to visualize your estimated filter\n" ] }, { "cell_type": "code", "metadata": { "id": "Z66V8ZXQRgYR", "cellView": "form", "tags": [ "hide-input" ] }, "source": [ "# @markdown Execute to visualize negative log likelihood over gradient descent\n", "fig, axes = plt.subplots()\n", "axes.plot(NLL,'-ok')\n", "axes.set(ylabel='NLL', xlabel='Gradient descent step', title='LNP fitting');" ], "execution_count": null, "outputs": [] }, { "cell_type": "code", "metadata": { "id": "by6E5eLOSOAP", "cellView": "form" }, "source": [ "# @markdown Execute to visualize your estimated filter\n", "\n", "fig, axes = plt.subplots(1, 2, figsize=(10, 5))\n", "\n", "axes[0].imshow(filter.reshape((10, 10)), vmin=-1, vmax=1, cmap='gray')\n", "axes[1].imshow(k.reshape((10, 10)), vmin=-1, vmax=1, cmap='gray')\n", "\n", "for i in range(2):\n", " axes[i].axis('Off')\n", "\n", "axes[0].set(title='True filter')\n", "axes[1].set(title='Estimate of k using NS data');\n" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "XJqFvwGkTpN9" }, "source": [ "## E) Larger steps\n", "\n", "In the next cell, try out performing gradient descent using your function above and step size alpha = $10^{-5}$ instead of $10^{-6}$. What happens with the negative log likelihood over time? Why is this happening?" ] }, { "cell_type": "code", "metadata": { "id": "E3DcEyJET4FR" }, "source": [ "k, NLL = ...\n", "\n", "fig, axes = plt.subplots()\n", "axes.plot(NLL,'-ok')\n", "axes.set(ylabel='NLL', xlabel='Gradient descent step', title='LNP fitting');" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "8rlT0lkYUFDJ" }, "source": [ "### **Answer**\n", "\n", "Text answer here\n", " " ] }, { "cell_type": "markdown", "metadata": { "id": "31vWDOMtSYrf" }, "source": [ "## Extra info\n", "\n", "We didn't need to compute gradient descent ourselves. We could have used an optimizer from scipy as shown in the following code. We computed our gradient by hand for practice and to really look \"under the hood\" of gradient descent." ] }, { "cell_type": "code", "metadata": { "id": "qouEIifr8r1Y" }, "source": [ "from scipy.optimize import minimize\n", "init_guess = np.zeros((10*10,))\n", "outs = minimize(compute_NLL, init_guess, (NS_images, NS_spike_counts))\n", "\n", "\n", "plt.imshow(outs.x.reshape((10, 10)), vmin=-1, vmax=1,cmap='gray')\n", "plt.axis('Off');" ], "execution_count": null, "outputs": [] } ] }