{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Efficiency of different stimulus ensembles for systems identification\n", "\n", "In reverse correlation, we approximate a nonlinear, stochastic function $y \\sim f(X)$ locally by a linear approximation weighted by a gaussian window. We compute this approximation by computing the response $y$ of the system for a variety of normally distributed inputs $X \\sim N(0, \\sigma^2)$. An estimate of the response is given by a sum of the stimuli weighted by the responses, $\\hat w = \\frac{1}{N} X^T y$. When these responses are action potentials, or spikes measures from biological neurons, the estimate $\\hat w$ is also called the spike-triggered average.\n", "\n", "How good is this estimator? Not very good. It can, however, be significantly improved by careful consideration of the properties of this Monte Carlo estimator, namely by changing the input ensemble to use antithetic sampling, or by shifting the distribution of the response.\n", "\n", "To demonstrate this, I use an example nonlinear function to be estimated via this black-box method. It consists of weighting the input with a windowed sinusoid, followed by an expansive nonlinearity that drives a Poisson process." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "pycharm": { "is_executing": false } }, "outputs": [], "source": [ "%config InlineBackend.figure_format = 'retina'\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import pandas as pd\n", "import plotnine\n", "import seaborn as sns\n", "\n", "sns.set(style=\"darkgrid\")\n", "\n", "\n", "class LNP:\n", " \"\"\"A simple LNP model neuron.\"\"\"\n", " def __init__(self):\n", " rg = np.arange(-31.5, 32.5)\n", " self.w = np.cos(rg / 20.0 * 2 * np.pi) * np.exp(-rg ** 2 / 2 / 10 ** 2)\n", " self.w /= 4\n", " self.input_size = len(self.w)\n", " self.rate_multiplier = 1\n", " \n", " def forward(self, X):\n", " \"\"\"The nonlinearity is $\\log(1 + \\exp(x))$\"\"\"\n", " return np.random.poisson(self.rate_multiplier * np.log(1 + np.exp(X.dot(self.w) - .5)))" ] }, { "cell_type": "code", "execution_count": 2, "outputs": [ { "data": { "text/plain": "Text(0.5, 1.0, 'Model weights')" }, "metadata": {}, "output_type": "execute_result", "execution_count": 2 }, { "data": { "text/plain": "