Supervised.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Chapter 3: Supervised Learning"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$ % Latex macros\n",
    "\\newcommand{\\mat}[1]{\\begin{pmatrix} #1 \\end{pmatrix}}\n",
    "\\newcommand{\\p}[2]{\\frac{\\partial #1}{\\partial #2}}\n",
    "\\newcommand{\\b}[1]{\\boldsymbol{#1}}\n",
    "\\newcommand{\\w}{\\boldsymbol{w}}\n",
    "\\newcommand{\\x}{\\boldsymbol{x}}\n",
    "\\newcommand{\\y}{\\boldsymbol{y}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Contents\n",
    "* 3.1 Function Approximation (Bishop 2006, Chater 3)\n",
    "    * Linear regression\n",
    "    * Least squares – Maximum likelihood\n",
    "    * Online learning – Stochastic gradient descent\n",
    "\n",
    "\n",
    "* 3.2 Pattern Recognition (Bishop 2006, Chapter 4)\n",
    "    * Perceptron\n",
    "    * Logistic regression  \n",
    "    * Support vector machine\n",
    "    \n",
    "\n",
    "* 3.3 Classical Conditining\n",
    "\n",
    "\n",
    "* 3.4 Cerebellum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "from mpl_toolkits.mplot3d import Axes3D"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The aim of *supervised learning* is to construct an input-output mapping\n",
    "\n",
    "$$ \\y = f(\\x) $$\n",
    "\n",
    "from pairs of samples $(\\x_1,\\y_1),...,(\\x_N,\\y_N)$.\n",
    "\n",
    "When $\\y$ is continuous, it is called *function approximation* or *regression*.\n",
    "\n",
    "When $\\y$ is discrete, it is called *pattern recognition* or *classification*.\n",
    "\n",
    "Here we focus on the case where the output is one dimension and computed from weighted sum of the inputs $x_1,...,x_D$ \n",
    "\n",
    "$$\n",
    "    y = f(\\x;\\w) = g( w_0 + w_1 x_1 + ... + w_D x_D).\n",
    "$$\n",
    "\n",
    "This is considered as an artificial neuron unit with input synaptic weights $w_1,...,w_D$, bias $w_0$ (or threshold $-w_0$), and the output function $g(\\ )$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Illustration of a neuron unit\n",
    "D = 5\n",
    "x = [r\"$x_1$\", r\"$x_2$\", r\"$\\vdots$ \",r\"$x_D$\"]\n",
    "w = [r\"$w_1$\", r\"$w_2$\", \"\",r\"$w_D$\"]\n",
    "for i in range(4):\n",
    "    plt.text(-0.2, 1.5-i, x[i], fontsize=18)  # inputs\n",
    "    plt.plot([0,1], [1.5-i,0], \"k\")  # inputs\n",
    "    plt.text(0.5, (1.5-i)*0.6, w[i], fontsize=18)  # weights\n",
    "plt.plot([1,2], [0,0], \"k\")   # output\n",
    "plt.plot(1, 0, \"ko\", markersize=50, markerfacecolor=\"w\")\n",
    "plt.text(1.6, 0.1, r\"$y$\", fontsize=18)  # output\n",
    "plt.axis(\"off\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.1 Function Approximation\n",
    "\n",
    "A typical case of function approximation in the brain is learning motor-sensory mapping; how does your arm respond to the activation of the muscles. In this case $\\x$ is the acitivation pattern of the arm muscles and $\\y$ is the changes in the arm joint angles, for example. Such an *internal model* of the body dynamics is very much helpful in motor control. In this case, the supervisor is your musculoskeletal system and sensory neurons."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear Regression\n",
    "\n",
    "In the simplest case, the output function is identity $g(u)=u$. This is the case of *linear regression*.\n",
    "\n",
    "For $D$-dimensional input \n",
    "\n",
    "$$\\x = (x_1,...,x_D)^T,$$\n",
    "\n",
    "we take a weight vector\n",
    "\n",
    "$$\\w = (w_0,w_1,...,w_D)^T$$\n",
    "\n",
    "and give a scalar output\n",
    "\n",
    "$$\n",
    "y = f(\\x;\\w) = w_0 + w_1 x_1 + ... + w_D x_D.\n",
    "$$\n",
    "\n",
    "By considering a constant input $x_0=1$ and redefining $\\x$ as \n",
    "\n",
    "$$\\x = (1,x_1,...,x_D)^T,$$\n",
    "\n",
    "we have a simple vector notation\n",
    "\n",
    "$$\n",
    "y = f(\\x;\\w) = \\w^T \\x = \\x^T \\w.\n",
    "$$\n",
    "\n",
    "We represent a set of input data as a matrix\n",
    "\n",
    "$$X = \\pmatrix{\\x_1^T \\\\ \\vdots \\\\ \\x_N^T}\n",
    "    = \\pmatrix{1 & x_{11} & \\cdots & x_{1D}\\\\\n",
    "    \\vdots &\\vdots & & \\vdots\\\\\n",
    "    1 & x_{N1} & \\cdots & x_{ND}} \n",
    "$$\n",
    "\n",
    "and the set of outputs is given in a vector form as\n",
    "\n",
    "$$\n",
    "\\y = \\pmatrix{y_1\\\\ \\vdots \\\\y_N} = X \\w.\n",
    "$$\n",
    "\n",
    "When a set of target outputs\n",
    "\n",
    "$$ \\y^* = \\pmatrix{y^*_1\\\\ \\vdots \\\\y^*_N}$$\n",
    "\n",
    "is given, how can we find the appropriate weight vector $\\w$ so that $\\y$ becomes close to $\\y^*$?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Least square solution\n",
    "\n",
    "A basic solution is to minimize the squared error between the target output and the model output\n",
    "\n",
    "$$ E(\\w) = \\frac{1}{2}\\sum_{n=1}^N (\\w^T \\x_n - y^*_n)^2 \n",
    " = \\frac{1}{2} ||X\\w - \\y^*||^2. $$\n",
    "\n",
    "To minimize this, we differentiate this error function with each weight parameter and set it equal to zero:\n",
    "\n",
    "$$\n",
    "    \\p{E(\\w)}{w_i} = \\sum_{n=1}^N (\\w^T \\x_n - y^*_n)x_{ni} = 0.\n",
    "$$\n",
    "\n",
    "These are represented in a vector form as\n",
    "\n",
    "$$\n",
    "    \\p{E(\\w)}{\\w} = X^T(X\\w - \\y^*) = (X^T X)\\w - X^T\\y^* = \\b{0}.\n",
    "$$\n",
    "\n",
    "Thus the solution to this minimization problem is given by\n",
    "\n",
    "$$\n",
    "    \\b{\\hat{w}} = (X^T X)^{-1} X^T \\y^*.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Maximum likelihood estimate\n",
    "\n",
    "A statstician would assume that the target output is generated by a linear model with additional noise\n",
    "\n",
    "$$\n",
    "    y^*_n = \\w^T\\x_n + \\epsilon_n\n",
    "$$\n",
    "\n",
    "where $\\epsilon_n$ is assumed to follow a Gaussian distribution $\\mathcal{N}(0,\\sigma^2)$.\n",
    "This is also re-written as\n",
    "\n",
    "$$\n",
    "    y^*_n \\sim \\mathcal{N}(\\w^T\\x_n,\\sigma^2)\n",
    "$$\n",
    "\n",
    "> The Gaussian distribution is defined as:\n",
    ">\n",
    ">$$\n",
    "    p(x;\\mu,\\sigma^2) = \\frac{1}{\\sqrt{2\\pi\\sigma^2}}e^{-\\frac{(x-\\mu)^2}{2\\sigma^2}}\n",
    "$$\n",
    "\n",
    "In this setup, the standard way of selecting the parameter is to find the one with the *maximum likelihood*, i.e. the probability of producing the data.\n",
    "\n",
    "The likelihood of weights $\\w$ for the set of observed data is the product of the likelihood for each data:\n",
    "\n",
    "$$\n",
    "    L(\\w) = p(\\y^*|X,\\w,\\sigma^2)\n",
    "    = \\prod_{n=1}^N (2\\pi\\sigma^2)^{-\\frac{1}{2}}\n",
    "    e^{-\\frac{(y^*_n-\\w^T\\x_n)^2}{2\\sigma^2}}.\n",
    "$$\n",
    "\n",
    "In maximizing the likelihood, it is mathematically and computationally more convenient to take its logarithm\n",
    "\n",
    "$$\n",
    "    l(\\w) = \\log p(\\y^*|X,\\w,\\sigma^2)\n",
    "    = \\sum_{n=1}^N -\\frac{1}{2}\\log(2\\pi\\sigma^2)\n",
    "    - \\frac{(y^*_n-\\w^T\\x_n)^2}{2\\sigma^2}\n",
    "$$\n",
    "\n",
    "$$\n",
    "    = - \\frac{N}{2}\\log(2\\pi) - N\\log\\sigma  - \\frac{1}{\\sigma^2} E(\\w).\n",
    "$$\n",
    "\n",
    "In the above, only the last term depends on the weights $\\w$.\n",
    "Thus the maximizing the likelihood is equivalent to minimizing the sum of squared errors $E(\\w)$.\n",
    "\n",
    "This link between minimal error and maximum likelihood is helpful when we consider regularization of parameters and Bayesian perspective in Chapter 6."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example: Arm dynamics\n",
    "\n",
    "Let us simulate simple second-order dynamics of an arm hanging down and take the data.\n",
    "Then we will make a linear regression model to predict the angular accleration from the angle and angular velocity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.integrate import odeint"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Illustration of an arm hanging down\n",
    "l = 0.3  # arm length: m\n",
    "th = 0.3 # angle: rad\n",
    "plt.plot([0,l*np.sin(th)], [0,-l*np.cos(th)], lw=10)\n",
    "plt.plot([0,0], [0,-l], \"k-\")  # downward line\n",
    "tex = plt.text(0.02, -0.5*l, r\"$\\theta$\", fontsize=16)\n",
    "plt.axis('equal');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Dynamics of the arm\n",
    "m = 5     # arm mass: kg\n",
    "mu = 0.1  # damping: Nm/(rad/s)\n",
    "g = 9.8   # gravity: N/s^2\n",
    "I = (m*l**2)/3 # inertia for a rod around an end: kg m^2\n",
    "def arm(x, t):\n",
    "    \"\"\"arm dynamics for odeint: x=[th,om]\"\"\"\n",
    "    th, om = np.array(x)  # for readability\n",
    "    # angular acceleration\n",
    "    aa = (-m*g*l/2*np.sin(th) - mu*om)/I\n",
    "    return np.array([om, aa])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Simulate for 10 sec.\n",
    "dt = 0.1  # time step\n",
    "t = np.arange(0, 10, dt)  # time points\n",
    "X = odeint(arm, [0.1, 0], t)\n",
    "plt.plot(t, X)\n",
    "plt.xlabel(\"time\")\n",
    "plt.legend((\"theta\",\"omega\"));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Acceleration by differentiation\n",
    "Y = (X[1:,1] - X[:-1,1])/dt  # temporal difference\n",
    "X = X[:-1,:]    # omit the last point\n",
    "N, D = X.shape  # data count and dimension\n",
    "# add observation noise\n",
    "sig = 1.0   # noise size\n",
    "Y = Y + sig*np.random.randn(N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# You may want to use notebook mode for rotating\n",
    "#%matplotlib notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = plt.figure().add_subplot(111, projection='3d')\n",
    "# Scatter plot in 3D\n",
    "ax.scatter(X[:,0], X[:,1], Y, c=Y);\n",
    "plt.xlabel(\"theta\"); plt.ylabel(\"omega\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Prepare data matrix\n",
    "X1 = np.c_[np.ones(N), X] #  add a column of 1's\n",
    "# Compute the weights: W = (X^T X)^(-1) X^T Y\n",
    "w = np.linalg.inv(X1.T@X1) @ X1.T@Y\n",
    "#w = np.linalg.solve(X1.T@X1, X1.T@Y)\n",
    "print(\"w =\", w)\n",
    "print([0, -m*g*l/2/I, -mu/I])  # analytic values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Squared error\n",
    "mse = np.sum((X1@w - Y)**2)/N  # mean squared error\n",
    "print(\"mse =\", mse)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Show regression surface\n",
    "ax = plt.figure().add_subplot(111, projection='3d')\n",
    "#  make a grid\n",
    "X2, Y2 = np.meshgrid(np.linspace(-0.1,0.1,5),np.linspace(-1,1,5))\n",
    "# show wireframe\n",
    "Z2 = w[0] + w[1]*X2 + w[2]*Y2\n",
    "ax.plot_wireframe(X2, Y2, Z2);\n",
    "# show the data\n",
    "ax.scatter(X[:,0], X[:,1], Y, c=Y);\n",
    "plt.xlabel(\"theta\"); plt.ylabel(\"omega\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Online Learning\n",
    "\n",
    "In the above, we assumed that all the data are available at once. However, data come in sequence and we would rather learn from the data as they come in.\n",
    "\n",
    "The basic way of online learning is to minimize the output error for each input \n",
    "\n",
    "$$\n",
    "    E_n(\\w) = \\frac{1}{2}(y^*_n - \\w^T \\x_n)^2\n",
    "$$\n",
    "\n",
    "and move $\\w$ down to its gradient\n",
    "\n",
    "$$\n",
    "    \\Delta \\w = - \\alpha \\p{E_n(\\w)}{\\w} = \\alpha(y^*_n - \\w^T \\x_n)\\x_n \n",
    "$$\n",
    "\n",
    "where $\\alpha>0$ is a learning rate paramter.\n",
    "\n",
    "This is called *stochastic gradient descent (SGD)*, by assuming that $x_n$ are stochastic samples from the input data space.\n",
    "\n",
    "Online learning has several advantages over *batch* algorithms like linear regrassion that processes all the data at once.\n",
    "Online learning does not require matrix inversion, which is computationally expensive when dealing with high dimensional data (large $D$).\n",
    "For a very large dataset (large $N$), simply storing a huge $N\\times D$ data matrix $X$ in the memory can be costly. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basis Functions\n",
    "\n",
    "For approximating a nonlinear function of $\\x$, a standard way is to prepare a set of nonlinear functions \n",
    "$$\n",
    "    \\b{\\phi}(\\x) = (\\phi_1(\\x),...,\\phi_M(\\x))^T,\n",
    "$$\n",
    "called basis functions, and represent the target function by\n",
    "$$\n",
    "    f(\\x;\\w) = \\sum_{j=1}^M w_j \\phi_j(\\x) = \\w^T\\b{\\phi}(\\x).\n",
    "$$\n",
    "\n",
    "Classically, polynomials $(1, x, x^2, x^3,...)$ or sinusoids $(1, \\cos x, \\cos 2x,..., \\sin x, \\sin 2x,...)$ were often used as basis functions, motivated by polynomial expansion or Fourier expansion.\n",
    "\n",
    "However, these functions can grow so lage or become so steep when we include higher order terms, which can make learned function chagning wildly."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Radial Basis Functions\n",
    "\n",
    "In sensory nervous systems, neurons tend to respond to signals in a limited range, called *receptive field*. For example, each visual neuron responds to light stimuli on a small area in the retina. Each auditory neurons repond to a limited range of frequency.\n",
    "\n",
    "Partly inspired by such local receptive field properties, a popular class of basis functions is called *radial basis function (RBF)*. \n",
    "An RBF is give by a decreasing function of the distance from a center point:\n",
    "\n",
    "$$\n",
    "    \\phi_j(\\x) = g(||\\x-\\x_j||)\n",
    "$$\n",
    "\n",
    "A typical example is Gaussian basis functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "# Gaussian basis functions in 1D\n",
    "def phi(x, c, s=1):\n",
    "    \"\"\"Gaussian centered at x=c with scale s\"\"\"\n",
    "    return np.exp(-((x - c)/s)**2)\n",
    "M = 10  # x range\n",
    "x = np.linspace(0, M, 100)\n",
    "for c in range(M):\n",
    "    plt.plot(x, phi(x, c, 1))\n",
    "plt.xlabel(\"x\"); plt.ylabel(\"$\\phi(x)$\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here you can generate random samples of functions using RBF."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "w = np.random.randn(M)  # random weights\n",
    "y = x*0  # output to cover the same range as input\n",
    "for c in range(M):\n",
    "    y = y + w[c]*phi(x, c, s=1)  # try chaning the scale\n",
    "plt.plot(x, y);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here are RBFs in 2D input space on a grid."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 2D Gaussian basis functions\n",
    "ax = plt.figure().add_subplot(111, projection='3d')\n",
    "# prepare a surface grid\n",
    "x = np.linspace(0, 10, 31)\n",
    "X, Y = np.meshgrid(x, x)\n",
    "# centers on another grid\n",
    "for cx in range(2, 10, 3):\n",
    "    for cy in range(2, 10, 3):\n",
    "        Z = phi(X, cx, 1)*phi(Y, cy, 3)\n",
    "        ax.plot_surface(X, Y, Z)\n",
    "plt.xlabel(\"x1\"); plt.ylabel(\"x2\"); plt.title(\"$\\phi(x)$\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Example: non-linear pendulum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.2 Pattern Recognition\n",
    "\n",
    "Here we consider a problem of classifying input vectors $\\x$ into $K$ discrete classes $C_k$.\n",
    "\n",
    "There are three major approaches in pattern classification.\n",
    "\n",
    "* Learn a discriminant function: $\\x \\rightarrow C_k$\n",
    "\n",
    "* Learn a conditional probability: $p(C_k|\\x)$\n",
    "\n",
    "* Learn a generative model $p(\\x|C_k)$ and then use Bayes' theorem:\n",
    "\n",
    "$$ p(C_k|\\x) = \\frac{p(\\x|C_k)p(C_k)}{p(\\x)} $$\n",
    "\n",
    "Supervised pattern recognition has been highly successful in applications of image and speech recognition. \n",
    "However, it is a question whether suprevised learning is relevant for our vision and speech, becuase we do not usually receive labeles for objects or words when we learn to see or hear as infants.\n",
    "*Unsupervised learning*, covered in Chapter 5, might be a more plausible way of human sensory learning."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Perceptron\n",
    "\n",
    "The first pattern classification learning machine was called *Perceptron* (Rosenblatt 1962). \n",
    "\n",
    "% ![Perceptron](figures/Perceptron.jpg)\n",
    "% The structure of a three-layer perceptron (from Rosenblatt, 1962).\n",
    "\n",
    "The original Perceptron consisted of three layers of binary units: S(sensory)-units, A(associative)-units connected randomly with S-units, and R(response)-units."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we formulate the Perceptron in a generalized form. \n",
    "The input vector $\\x$ is converted by a set of basis functions $\\phi_i(\\x)$ into a feature vector\n",
    "\n",
    "$$\n",
    "    \\b{\\phi}(\\x)=(\\phi_1(\\x),...,\\phi_M(\\x))^T.\n",
    "$$\n",
    "\n",
    "In the simplest case, the feature vector is the same as the input vector $\\b{\\phi}(\\x)=\\x$, or just augumented by $1$ to represent bias as $\\b{\\phi}(\\x)=\\pmatrix{1 \\\\ \\x}$.\n",
    "This is called *linear Perceptron*.\n",
    "\n",
    "The output is given by\n",
    "\n",
    "$$\n",
    "    y = f( \\sum_{i=1}^M w_{i} \\phi_i(\\x)) = f( \\w^T \\b{\\phi}(\\x))\n",
    "$$\n",
    "where $\\w=(w_1,...,w_M)^T$ is the output connection weights.\n",
    "\n",
    "The output function takes +1 or -1:\n",
    "\n",
    "$$\n",
    "    f(u) = \\begin{cases} 1 & \\mbox{if } u \\ge 0 \\\\ -1 & \\mbox{if } u<0. \\end{cases}\n",
    "$$\n",
    "\n",
    "For each input $\\x_n$, the target output $y^*_n \\in \\{+1,-1\\}$ is given."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Perceptron Learning Rule\n",
    "\n",
    "Learning of Perceptron is based on the error function\n",
    "\n",
    "$$\n",
    "    E(\\w) = \\sum_{y_n \\ne y^*_n} -y^*_n\\w^T\\b{\\phi}(\\x_n)\n",
    "$$\n",
    "which takes a positive value for each missclassified output $y_n \\ne y^*_n$.\n",
    "\n",
    "The perceptron learning rule is the stochastic gradient\n",
    "\n",
    "$$\n",
    "    \\Delta \\w = -\\alpha\\p{E(\\w)}{\\w}\n",
    "$$\n",
    "\n",
    "$$\n",
    "    = \\alpha y^*_n\\b{\\phi}(\\x_n) \\ \\mbox{ if } y_n\\ne y^*_n\n",
    "$$\n",
    "\n",
    "$$\n",
    "    = \\frac{\\alpha}{2}(y^*_n-y_n)\\b{\\phi}(\\x_n)\n",
    "$$\n",
    "which is a product of the output error and the input to each weight.\n",
    "\n",
    "When two classes are *linearly separable*, the *Preceptron converence theorem* assures that learning converges to find a proper hyperplane to separate two classes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example: Perceptron in 2D feature space "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Perceptron:\n",
    "    \"\"\"Linear perceptron: phi(x)=[1,x]\"\"\"\n",
    "    \n",
    "    def __init__(self, D=2):\n",
    "        \"\"\"Create a new perceptron\"\"\"\n",
    "        # self.w = np.random.randn(D+1)  # output weight\n",
    "        self.w = np.zeros(D+1)  \n",
    "        self.vis = False\n",
    "\n",
    "    def output(self, x):\n",
    "        \"\"\"predict an output from input x\"\"\"\n",
    "        u = self.w@np.r_[1,x]\n",
    "        y = 1 if u>0 else -1\n",
    "        if self.vis:\n",
    "            plt.plot(x[0], x[1], \"yo\" if y>0 else \"bo\")\n",
    "        return y\n",
    "\n",
    "    def learn(self, x, yt, alpha=0.1):\n",
    "        \"\"\"learn from (x, yt) pair\"\"\"\n",
    "        y = self.output(x)\n",
    "        error = y - yt\n",
    "        if error != 0:\n",
    "            self.w += alpha*yt*np.r_[1,x]\n",
    "        if self.vis:\n",
    "            plt.plot(x[0], x[1], \"y*\" if yt>0 else \"b*\")\n",
    "            self.plot_boundary()\n",
    "        return error\n",
    "    \n",
    "    def plot_boundary(self, u=0, col=\"k\", l=5):\n",
    "        \"\"\"plot decision boundary with shift u, length l\"\"\"\n",
    "        # weight vector\n",
    "        plt.plot([0,self.w[1]], [0,self.w[2]], \"r\", lw=2)\n",
    "        plt.plot(self.w[1], self.w[2], \"r*\")  # arrowhead\n",
    "        # decision boundary\n",
    "        s = self.w[1]**2 + self.w[2]**2  # weight size\n",
    "        x = self.w[1]*(u-self.w[0])/s+self.w[2]*np.array([l,-l])/np.sqrt(s)\n",
    "        y = self.w[2]*(u-self.w[0])/s+self.w[1]*np.array([-l,l])/np.sqrt(s)\n",
    "        plt.plot(x, y, col, lw=0.5)\n",
    "        plt.xlabel(\"x1\"); plt.ylabel(\"x2\"); plt.axis(\"equal\"); \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Test by 2D gaussian data\n",
    "N = 10  # sample size\n",
    "D = 2\n",
    "Xp = np.random.randn(N,D) + np.array([3,2])  # positive data\n",
    "Xn = np.random.randn(N,D)   # negative data\n",
    "# concatenate positive/negative data\n",
    "X = np.r_[Xp,Xn]\n",
    "Yt = np.r_[np.ones(N),-np.ones(N)]  # target\n",
    "N = 2*N\n",
    "# plot the data\n",
    "plt.scatter(X[:,0], X[:,1], c=Yt, marker=\"o\")\n",
    "plt.xlabel(\"x1\"); plt.ylabel(\"x2\"); plt.axis(\"equal\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a Perceptron\n",
    "perc = Perceptron(D=2)  # perceptron with 2D input\n",
    "perc.vis = True  # turn on visualization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Learning: run this cell several times\n",
    "err = 0\n",
    "# shuffle data order\n",
    "for i in np.random.permutation(N):\n",
    "    err += perc.learn(X[i], Yt[i])**2\n",
    "print(\"mse =\", err/N, \"; w = \", perc.w)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Support Vector Machine\n",
    "\n",
    "With Perceptron learning rule, the line separating two classess can end up in any configuration so that two classes are separated, which may not be good for *generalization*, i.e., classification of new data that were not used for training.\n",
    "\n",
    "The *support vector machine (SVM)* (Bishop 2006, section 7.1) tries to find a separation line that has the largest margin to the positive and negative data by minimizing the *hinge* objective function\n",
    "\n",
    "$$\n",
    "    E(\\w) = \\sum_n \\max(0, 1-y^*_n\\w^T\\b{\\phi}(\\x_n)) + \\lambda||\\w||.\n",
    "$$\n",
    "\n",
    "This tries to make the weighted input sum $\\w^T\\b{\\phi}(\\x_n)>1$ for $y^*_n=1$ and $\\w^T\\b{\\phi}(\\x_n)<1$ for $y^*_n=-1$,\n",
    "while keeping $\\w$ small to create a large margin between the lines for $\\w^T\\b{\\phi}(\\x_n))=1$ and $\\w^T\\b{\\phi}(\\x_n)) =-1$.\n",
    "\n",
    "A standard way for solving this optimization problem is to use a batch optimization method called *quadratic programming*.\n",
    "It is often the case that the basis functions are allocated around each data point $\\x_n$ using so-called *kernel* function\n",
    "\n",
    "$$\n",
    "    \\phi_i(\\x) = K(\\x,\\x_i).\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "source": [
    "### Online SVM\n",
    "\n",
    "Here we consider an online version of SVM called *Pagasos* (Shalev-Shwartz et al. 2010) that has a closer link with Perceptron learning.\n",
    "\n",
    "The weights are updated by\n",
    "\n",
    "$$\n",
    "    \\Delta\\w = \\alpha\\{1[1-y^*_n\\w^T\\b{\\phi}(\\x_n)] y^*_n\\b{\\phi}(\\x_n) - \\lambda\\w\\}\n",
    "$$\n",
    "where $1[\\ ]$ represents an indicator function; $1[u]=1$ if $u>0$ and $1[u]=0$ otherwise.\n",
    "\n",
    "The learning rate is gradually decreased as $\\alpha=\\frac{1}{\\lambda n}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Pegasos(Perceptron):\n",
    "    \"\"\"Pagasos, online SVM with linear kernel\"\"\"\n",
    "    \n",
    "    def __init__(self, D):\n",
    "        super().__init__(D)\n",
    "        self.n = 0  # counter for rate scheduling\n",
    "        \n",
    "    def learn(self, x, yt, lamb=1.0):\n",
    "        \"\"\"learn from (x, y) pair\"\"\"\n",
    "        y = self.output(x)\n",
    "        self.n += 1  # data count\n",
    "        alpha = 1/(lamb*self.n)  # adapt learning rate\n",
    "        u = self.w@np.r_[1,x]  # input sum\n",
    "        # hinge loss and regularization except bias\n",
    "        self.w += alpha*(((yt*u)<1)*yt*np.r_[1,x] - lamb*np.r_[0,self.w[1:]])\n",
    "        if self.vis:  # for first 2D\n",
    "            # target output\n",
    "            plt.plot(x[0], x[1], \"y*\" if yt>0 else \"b*\")\n",
    "            self.plot_boundary(u=1, col=\"g\")\n",
    "            self.plot_boundary(u=0, col=\"k\")\n",
    "            self.plot_boundary(u=-1, col=\"m\")\n",
    "        return max(0, 1 - yt*u)  # hinge loss"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a Pagasos\n",
    "pega = Pegasos(D=2)  # 2D input\n",
    "pega.vis = True  # turn on visualization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Learning: run this cell several times\n",
    "err = 0\n",
    "# shuffle data order\n",
    "for i in np.random.permutation(N):\n",
    "    err += pega.learn(X[i], Yt[i], lamb=1)\n",
    "print(\"n =\", pega.n, \"; err =\", err/N, \"; w =\", pega.w)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Logistic Regression\n",
    "\n",
    "In Perceptron and SVM, the output is binary. We sometimes want to express the certainty in the output by the probability of the data belonging to a class $p(C_k|\\x)$.\n",
    "\n",
    "*Logistic regression* is a probabilistic generalization of linear regression (Bishop 2006, Section 4.3). Its probabilistic output is given by\n",
    "\n",
    "$$\n",
    "    p(C_1|\\x) = y = g(\\w^T\\x))\n",
    "$$\n",
    "and $p(C_0|\\x) = 1 - p(C_1|\\x)$. \n",
    "\n",
    "The function $g(\\ )$ is called *logistic sigmoid function*\n",
    "\n",
    "$$\n",
    "    g(u) = \\frac{1}{1+e^{-u}}.\n",
    "$$\n",
    "\n",
    "> The derivative of the logistic sigmoid function is represented as\n",
    ">\n",
    ">$$\n",
    "    g'(u) = \\frac{-1}{(1+e^{-u})^2}(-e^{-u})\n",
    "    = \\frac{1}{1+e^{-u}}\\frac{e^{-u}}{1+e^{-u}}\n",
    "    = g(u)(1-g(u))\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Maximum likelihood by iterative least squares\n",
    "\n",
    "We can find the weights of logistic regression by the principle of maximum likelihood, the probability of reproducing the data.\n",
    "\n",
    "The likelihood of weights $\\w$ for a set of observed data $(X,\\y^*)$ is\n",
    "\n",
    "$$\n",
    "    L(\\w) = p(\\y^*|X,\\w)\n",
    "    = \\prod_{y^*_n=1} y_n \\prod_{y^*_n=0}(1-y_n)\n",
    "    = \\prod_{n=1}^N y_n^{y^*_n} (1-y_n)^{1-y^*_n}.\n",
    "$$\n",
    "\n",
    "The negative log likelihood is often called *cross entropy error*\n",
    "\n",
    "$$\n",
    "    E(\\w) = -l(\\w) = -\\log p(\\y^*|X,\\w)\n",
    "    = -\\sum_{n=1}^N\\{y^*_n\\log y_n + (1-y^*_n)\\log(1-y_n)\\}.\n",
    "$$\n",
    "\n",
    "From\n",
    "\n",
    "$$\n",
    "    \\p{y}{\\w} = g'(\\w^T\\x)\\p{\\w^T\\x}{\\w} = y(1-y)\\x,\n",
    "$$\n",
    "\n",
    "the gradient of the cross entropy error is given as\n",
    "\n",
    "$$\n",
    "    \\p{E(\\w)}{\\w} \n",
    "    = -\\sum_{n=1}^N\\{y^*_n(1-y_n)\\x_n - (1-y^*_n)y_n\\x_n\\}\n",
    "$$\n",
    "\n",
    "$$\n",
    "    = \\sum_{n=1}^N (y_n - y^*_n)\\x_n\n",
    "    = X^T(\\y - \\y^*) \n",
    "$$\n",
    "which takes the same form as in linear regression."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "source": [
    "Unlike in linear regression, $\\p{E(\\w)}{\\w}=0$ does not have a closed form solution as in linear regression. \n",
    "Instead, we can apply Newton method using the *Hessian matrix*\n",
    "\n",
    "$$\n",
    "    H = \\p{}{\\w}\\p{E(\\w)}{\\w} = X^T R X\n",
    "$$\n",
    "where $R$ is a diagonal matrix made of the derivative of the sigmoid functions\n",
    "\n",
    "$$\n",
    "    R = \\mbox{diag}(y_1(1-y_1),...,y_n(1-y_n)).\n",
    "$$\n",
    "\n",
    "The update is made by\n",
    "\n",
    "$$\n",
    "    \\w := \\w - (X^TRX)^{-1}X^T(\\y - \\y^*)\n",
    "    = (X^TRX)^{-1}X^T R \\b{z}\n",
    "$$\n",
    "where $\\b{z}$ is an effective target value\n",
    "\n",
    "$$\n",
    "    \\b{z} = X\\w - R^{-1}(\\y - \\y^*).\n",
    "$$\n",
    "\n",
    "This algorithm is called *iterative reweighted least squares*."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Online logistic regression\n",
    "\n",
    "From the gradient of negative log likelihood for each data point\n",
    "\n",
    "$$\n",
    "    \\p{E_n(\\w)}{\\w} = (y_n - y^*_n)\\x_n,\n",
    "$$\n",
    "\n",
    "we can also apply stochastic gradient descent\n",
    "\n",
    "$$\n",
    "    \\Delta \\w = - \\alpha \\p{E_n(\\w)}{\\w} = \\alpha (y^*_n - y_n)\\x_n.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class LogisticRegression(Perceptron):\n",
    "    \"\"\"Logistic regression\"\"\"\n",
    "    \n",
    "    def sigmoid(self, u):\n",
    "        return 1/(1+np.exp(-u))\n",
    "    \n",
    "    def output(self, x):\n",
    "        \"\"\"output for vector/matrix input x\"\"\"\n",
    "        if x.ndim == 1:  # x is a vector\n",
    "            self.X = np.r_[1,x]\n",
    "        else:   # x is a matrix\n",
    "            self.X = np.c_[np.ones(x.shape[0]),x]\n",
    "        u = self.X@self.w\n",
    "        y = self.sigmoid(u)   # sigmoid output\n",
    "        if self.vis:\n",
    "            if x.ndim == 1:  # x is a vector\n",
    "                plt.plot(x[0], x[1], \"yo\" if y>0.5 else \"bo\")\n",
    "            else:   # x is a matrix\n",
    "                plt.scatter(x[:,0], x[:,1], c=y)\n",
    "        return y\n",
    "    \n",
    "    def rls(self, x, yt, alpha=0.1):\n",
    "        \"\"\"reweighted least square with (x, y) pairs\"\"\"\n",
    "        y = self.output(x)  # also set self.X\n",
    "        error = y - yt\n",
    "        R = np.diag(y*(1-y)) # weighting matrix\n",
    "        self.w -= np.linalg.inv(self.X.T@R@self.X) @ self.X.T@error\n",
    "        if self.vis:\n",
    "            plt.scatter(x[:,0], x[:,1], c=yt, marker=\"*\")\n",
    "            self.plot_boundary(u=1, col=\"g\")\n",
    "            self.plot_boundary(u=0, col=\"k\")\n",
    "            self.plot_boundary(u=-1, col=\"m\")\n",
    "        return error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "logreg = LogisticRegression(D=2)\n",
    "logreg.vis = True\n",
    "Yb = np.clip(Yt, 0, 1)  # target in [0,1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Repeat for Iterative Reweighted least squares\n",
    "err = logreg.rls(X, Yb)\n",
    "print(\" err =\", np.mean(err**2), \"; w =\", logreg.w)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.3 Classical Conditioning\n",
    "Animals have innate behaviors to avoid dangers or to acquire rewards, such as blinking when a strong light or wind hits the eyes or drooling when a food is expected. \n",
    "Such response is called a *unconditioned response (UR)* and the sensory cue to elicit UR is called *unconditioned stimulus (US)*.\n",
    "\n",
    "When an US is repreatedly preceded by another stimulus, such as a sound, animals start to take the same response as UR even before US is presented.\n",
    "In this case the response is called a *conditioned response (CR)* and the sensory cue to elicit CR is called *conditioned stimulus (CS)*.\n",
    "This type of learning is called *classical conditioning* or *Pavlovian conditioning* and can be considered as an example of supervised learning.\n",
    "\n",
    "The learned mapping can be CS to CR, or CS to US which elicits UR=CR.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Eye Blink Conditioning\n",
    "\n",
    "When an airpuff (US) is given around eyes, an animal would blink its eyes (UR). If the airpuff is preceded with a consistent interval by a cue, such as tone (CS), the animal learns to blink (CR) just before the airpuff is given.\n",
    "\n",
    "% ![EyeBlink](figures/McCormick1984jns.jpg)\n",
    "% An example of eye blink conditioning. The curve on top is the response of the eyelid. The PSTH below shows the firing of an output neuron in the cerebellum (from McCormick & Thompson, 1984).\n",
    "\n",
    "Richard Thompson and colleagues investigated the neural circuit behind this eye-blink conditioning and identified the cerebellum as the major locus of learning (Thompson, 1986).\n",
    "\n",
    "While blockade of the cerebellum does not affect eye-blink response itself (US-UR), it blocks learning (CS-CR). The activity of the output neurons of the cerebellum increase their activities as learning progress, as shown above.\n",
    "\n",
    "% ![CerebellumConditioning](figures/Thompson1986s.jpg)\n",
    "% The neural circuit of eye-blink conditioning (from Thompson, 1986). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.4 The Cerebellum\n",
    "\n",
    "The neural circuit of the cerebellum has a distinctive orthogonal structure. There are a massive number of *granule cells*, each of which collects inputs from a major input to the cerebellum, *the mossy fibers*.\n",
    "Granule cell axons run parallelly in the lateral direction in the cerebellar cortex, called *parallel fibers*.\n",
    "The output neurons of the cerebellum, *Purkinje cells*, spread fan-shaped dendrites in the longitudinal direction and receive a large nubmer of parallel fibers.\n",
    "Furthermore, each Purkinje cell is twined by a single axon called *climbing fiber* from a brain stem nucleus called the *inferior olive*.\n",
    "\n",
    "% ![CerebellarCortex](figures/Eccles1967.jpg)\n",
    "% The circuit of the cerebellar cortex (from Eccles et al. 1967)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cerebellar Perceptron Hyopthesis\n",
    "\n",
    "This peculiar organization inspired theorists around 1970 to come up with an idea that the cerebellum may work similar to the Perceptron (Marr 1969; Albus 1971).\n",
    "The massive number of granule cells may create a high-dimensional feature representation of the mossy fiber input.\n",
    "The single climbing fiber input may sever as the supervising signal to induce plasticity in the parallel-fiber input synapses to the Purkinje cell.\n",
    "\n",
    "% ![CerebellarPerceptron](figures/CerebellarPerceptron.jpg)\n",
    "% The cerebellar perceptron model (from Albus 1971).\n",
    "\n",
    "This hypothesis further motivated neurobiologists to test that experimentally.\n",
    "Indeed, cerebellar synaptic plasticity guided by the climbing fiber input was experimentally confirmed by Masao Ito (1984, 2000)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cerebellar Internal Model Hyopthesis\n",
    "\n",
    "The cerebellum is known to be important for movement control.\n",
    "A possible way in which suprevised learning can be helpful is to provide an *internal model* of the motor apparatus (Wolpert et al. 1998).\n",
    "\n",
    ">{figure-md} Fig_FEL\n",
    "<img src=\"https://pub.mdpi-res.com/biomimetics/biomimetics-05-00019/article_deploy/html/images/biomimetics-05-00019-g005.png?1590743622\" alt=\"Feedback Error Learning\" width=\"500px\">\n",
    ">\n",
    ">The feedback error learning model (Kawato et al. 1987). The innate feedback controller in the brainstem or spinal cord provides corrective motor signal, which is used as the error signal for supervised learning of *inverse model* in the cerebellum (from Ellery 2020)\n",
    "\n",
    "The molecular mechanism for the association of the parallel fiber input and the climbing fiber input (output error) at the Purkinje cell synapses has also been studied in detail (Ito 2000; Ogasawara et al. 2008)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "source": [
    "## References\n",
    "\n",
    "* Bishop CM (2006) Pattern Recognition and Machine Learning. Springer.  https://www.microsoft.com/en-us/research/people/cmbishop/prml-book/  \n",
    "    * Chapter 3: Linear models for regression  \n",
    "    * Chapter 4: Linear models for classification\n",
    "\n",
    "* Rosenblatt F (1962) Principles of Neurodynamics: Perceptrons and the Theory of Brain Mechanisms. Spartan. https://babel.hathitrust.org/cgi/pt?id=mdp.39015039846566&view=1up&seq=5\n",
    "\n",
    "### Classical conditioning\n",
    "\n",
    "* McCormick DA, Thompson RF (1984) Neuronal responses of the rabbit cerebellum during acquisition and performance of a classically conditioned nictitating membrane-eyelid response. Journal of Neuroscience, 4 (11) 2811-2822. https://doi.org/10.1523/JNEUROSCI.04-11-02811.1984\n",
    "* Thompson RF (1986) The neurobiology of learning and memory. Sicence, 233, 941-947. https://doi.org/10.1126/science.3738519\n",
    "\n",
    "### Cerebellar perceptron hypothesis\n",
    "\n",
    "* Marr D (1969) A theory of cerebellar cortex. Journal of Physiology, 202:437–470.  https://doi.org/10.1113/jphysiol.1969.sp008820\n",
    "* Albus JS (1971) A theory of cerebellar function. Mathematical Bioscience, 10:25–61, 1971. https://doi.org/10.1016/0025-5564(71)90051-4\n",
    "* Ito M (1984) The Cerebellum and Neural Control, Raven Press. https://archive.org/details/cerebellumneural0000itom\n",
    "\n",
    "### Cerebellar internal model hypothesis\n",
    "* Kawato M, Furukawa K, Suzuki R (1987). A hierarchical neural-network model for control and learning of voluntary movement. Biol Cybern, 57, 169-85. https://doi.org/10.1007/BF00364149\n",
    "* Wolpert DM, Miall RC, Kawato M (1998) Internal models in the cerebellum. Trends in Cognitive Sciences, 2:338–347. https://doi.org/10.1016/S1364-6613(98)01221-2\n",
    "* Ito M (2000) Mechanisms of motor learning in the cerebellum. Brain Research 886, 237–245. https://doi.org/10.1016/S0006-8993(00)03142-5\n",
    "* Ogasawara H, Doi T, Kawato M (2008) Systems biology perspectives on cerebellar long-term depression. Neurosignals, 16, 300–317.\n",
    "https://doi.org/10.1159/000123040"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}