<jupyter_start><jupyter_text>Week 8: Probability & Statistics, Statistical Encoding & Decoding Tutorial 1 [insert your name]**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".To complete this tutorial, you should have watched Video 8.1 and 8.2.<jupyter_code># @markdown Imports

# Imports
import numpy as np
import matplotlib.pyplot as plt
import ipywidgets as widgets  # interactive display
import math
# @markdown Plotting functions
import numpy
from numpy.linalg import inv, eig
from math import ceil
from matplotlib import pyplot, ticker, get_backend, rc
from mpl_toolkits.mplot3d import Axes3D
from itertools import cycle


%config InlineBackend.figure_format = 'retina'
plt.style.use("https://raw.githubusercontent.com/NeuromatchAcademy/course-content/master/nma.mplstyle")

# @markdown Helper functions
def twoD_Gaussian(xdata_tuple, amplitude, xo, yo, sigma_x, sigma_y, theta, offset):
  """Create 2D Gaussian based on parameters

  Args:
    xdata_tuple (ndarray): grid of x and y values to compute Gaussian for
    amplitude (scalar): amplitude of Gaussian
    xo (scalar): center of Gaussian in x coordinates
    yo (scalar): center of Gaussian in y coordinates
    sigma_x (scalar): standard deviation of Gaussian in x direction
    sigma_y (scalar): standard deviation of Gaussian in y direction
    theta (scalar): rotation angle of Gaussian
    offset (scalar): offset of all Gaussian values

  Returns:
    ndarray: Gaussian values at every x/y point

  """
    (x, y) = xdata_tuple
    xo = float(xo)
    yo = float(yo)    
    a = (np.cos(theta)**2)/(2*sigma_x**2) + (np.sin(theta)**2)/(2*sigma_y**2)
    b = -(np.sin(2*theta))/(4*sigma_x**2) + (np.sin(2*theta))/(4*sigma_y**2)
    c = (np.sin(theta)**2)/(2*sigma_x**2) + (np.cos(theta)**2)/(2*sigma_y**2)
    g = offset + amplitude*np.exp( - (a*((x-xo)**2) + 2*b*(x-xo)*(y-yo)+c*((y-yo)**2)))
    return g.ravel()<jupyter_output><empty_output><jupyter_text>The dataIn 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. 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.<jupyter_code># @markdown Execute this cell to get and visualize example image of each type (may take a few min to download)

np.random.seed(123)
n_images = 10000

# Get WN images
WN_images = np.random.randn(n_images, 10*10)

# Get NS images
from sklearn.datasets import fetch_openml
mnist = fetch_openml(name='mnist_784')
mnist_images = mnist.data
mnist_images = mnist_images/255 
mnist_images = mnist_images - np.mean(mnist_images)

mnist_images = mnist_images.reshape((-1, 28, 28))[:, 4:24, 4:24]
mnist_images = mnist_images[:, fc00:db20:35b:7399::5, ::2]
NS_images = mnist_images[:n_images].reshape((-1, 10*10))

fig, axes = plt.subplots(1, 2, figsize=(10, 5))
axes[0].imshow(WN_images[0].reshape((10, 10)), vmin=-1, vmax=1, cmap='gray')
axes[1].imshow(NS_images[0].reshape((10, 10)), vmin=-1, vmax=1, cmap='gray')
axes[0].axis('Off')
axes[1].axis('Off')
axes[0].set(title='WN example image')
axes[1].set(title='NS example image');<jupyter_output><empty_output><jupyter_text>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 simulated 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...)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.<jupyter_code># @markdown Execute to simulate neural responses
np.random.seed(0)

x = np.arange(-5, 5, 1)
y = np.arange(-5, 5, 1)
x, y = np.meshgrid(x, y)

sd_x = 1.4
sd_y = .5
gauss = twoD_Gaussian((x,y), 1, 0, 0, sd_x, sd_y, 35, 0)
filter = gauss.reshape((10, 10))

WN_lambda = np.exp(np.dot(WN_images, filter.reshape((-1,))))
WN_spike_counts = np.random.poisson(WN_lambda)

NS_lambda = np.exp(np.dot(NS_images, filter.reshape((-1,))))
NS_spike_counts = np.random.poisson(NS_lambda)

fig, axes = plt.subplots(1, 2, figsize=(10, 5), sharex=True)
axes[0].plot(WN_spike_counts[0:100], 'ok')
axes[1].plot(NS_spike_counts[0:100], 'ok')
axes[0].set(ylabel='Spike counts', xlabel='Image number', title='WN')
axes[1].set(xlabel='Image number', title='NS');<jupyter_output><empty_output><jupyter_text>Exercise 1: Computing an STAWe 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.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.  AnswerFill out code below<jupyter_code>def compute_STA(images, spike_counts):

  STA = ...

  return STA


WN_STA = ...
NS_STA = ...

# @markdown Execute to visualize your computed STAs and the original filter

fig, axes = plt.subplots(1, 3, figsize=(15, 5))

axes[0].imshow(filter.reshape((10, 10)), vmin=-1, vmax=1, cmap='gray')
axes[1].imshow(WN_STA.reshape((10, 10)), vmin=-1, vmax=1, cmap='gray')
axes[2].imshow(NS_STA.reshape((10, 10)), vmin=-1, vmax=1, cmap='gray')

for i in range(3):
  axes[i].axis('Off')

axes[0].set(title='True filter')
axes[1].set(title='White noise STA')
axes[2].set(title='Natural scenes STA');<jupyter_output><empty_output><jupyter_text>Note that the white noise STA is a pretty good estimate for the true filter, but the natural scenes STA is not!  (Optional) Exercise: Estimate the nonlinearityEstimate the nonlinearity of the LNP model (so no longer predefine it as exponential) using the method discussed in Video 8.2.<jupyter_code><jupyter_output><empty_output><jupyter_text>Exercise 2: Numerically finding the filter with natural scenes dataThe 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.  A) Negative log likelihood equationTo implement gradient descent ourselves, we will need to compute the derivative of the negative log likelihood. 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.Make the final equation clear either with the green text below or by putting it in a box or otherwise highlighting it.  **Answer**Put NLL = ... equation here (show work above or below)   B) Negative log likelihood computationUse 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).  AnswerFill out code below<jupyter_code>def compute_NLL(k, x, y):

    NLL = ...
    
    return NLL<jupyter_output><empty_output><jupyter_text>C) Compute dNLL/dkTake 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}}$.When in confusion about dot products, my recommendation is to write out the first few elements of the dot product computation for clarifiation.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!  **Answer**Put dNLL/dk = ... equation here (show work above or below)   D) Implementing gradient descentWe now have all the tools we need to implement gradient descent to find an estimate of our filter k using the natural scenes data.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.  AnswerFill out code below<jupyter_code>def gradient_descent(x, y, init_guess, n_steps = 500, alpha=10**-6):

    k = init_guess
    NLL = np.zeros((n_steps,))
    for i_step in range(n_steps):

      # Update estimate of k (assign as k)
      # your code here

      # Compute NLL at each step
      NLL[i_step] = compute_NLL(k, x, y)

    return k, NLL

k, NLL = ...
# @markdown Execute to visualize negative log likelihood over gradient descent
fig, axes = plt.subplots()
axes.plot(NLL,'-ok')
axes.set(ylabel='NLL', xlabel='Gradient descent step', title='LNP fitting');
# @markdown Execute to visualize your estimated filter

fig, axes = plt.subplots(1, 2, figsize=(10, 5))

axes[0].imshow(filter.reshape((10, 10)), vmin=-1, vmax=1, cmap='gray')
axes[1].imshow(k.reshape((10, 10)), vmin=-1, vmax=1, cmap='gray')

for i in range(2):
  axes[i].axis('Off')

axes[0].set(title='True filter')
axes[1].set(title='Estimate of k using NS data');<jupyter_output><empty_output><jupyter_text>E) Larger stepsIn 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?<jupyter_code>k, NLL = ...

fig, axes = plt.subplots()
axes.plot(NLL,'-ok')
axes.set(ylabel='NLL', xlabel='Gradient descent step', title='LNP fitting');<jupyter_output><empty_output><jupyter_text>**Answer**Text answer here   Extra infoWe 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.<jupyter_code>from scipy.optimize import minimize
init_guess = np.zeros((10*10,))
outs = minimize(compute_NLL, init_guess, (NS_images, NS_spike_counts))


plt.imshow(outs.x.reshape((10, 10)), vmin=-1, vmax=1,cmap='gray')
plt.axis('Off');<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10569836419198274, 'fraction_numerical': 0.021211576487506743, 'mean_word_length': 4.959828602035351, 'pattern_counts': {'":': 0, '<': 37, '<?xml version=': 0, '>': 39, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '618783', 'license': 'MIT', 'n_tokens_mistral': 3698, 'n_tokens_neox': 3426, 'n_words': 1377, 'path': 'Week8/Week8Tutorial1.ipynb', 'repo_name': 'felipe-parodi/QuantTools4Neuro', 'size': '11108'}