{ "cells": [ { "cell_type": "markdown", "metadata": { "toc": "true" }, "source": [ "# Table of Contents\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Requirements and helper functions" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Requirements\n", "\n", "This notebook requires to have [`numpy`](https://www.numpy.org/) and [`matplotlib`](https://matplotlib.org/) installed.\n", "One function needs a function from [`scipy.special`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.comb.html#scipy.special.comb).\n", "[`joblib`](https://joblib.readthedocs.io/en/latest/) is used to have parallel computations (at the end).\n", "\n", "The bottleneck performance of the main functions are very simple functions, for which we can write efficient versions using either [`numba.jit`](https://numba.pydata.org/numba-doc/latest/reference/jit-compilation.html#numba.jit) or [`cython`](https://cython.readthedocs.io/en/latest/src/quickstart/build.html#jupyter-notebook)." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Requirement already satisfied: watermark in /usr/local/lib/python3.6/dist-packages (1.5.0)\n", "Requirement already satisfied: numpy in /usr/local/lib/python3.6/dist-packages (1.15.4)\n", "Requirement already satisfied: scipy in /usr/local/lib/python3.6/dist-packages (1.1.0)\n", "Requirement already satisfied: matplotlib in /usr/local/lib/python3.6/dist-packages (3.0.2)\n", "Requirement already satisfied: joblib in /usr/local/lib/python3.6/dist-packages (0.12.1)\n", "Requirement already satisfied: numba in /usr/local/lib/python3.6/dist-packages (0.41.0)\n", "Requirement already satisfied: cython in /usr/local/lib/python3.6/dist-packages (0.29.1)\n", "Requirement already satisfied: ipython in /usr/local/lib/python3.6/dist-packages (from watermark) (7.2.0)\n", "Requirement already satisfied: cycler>=0.10 in /usr/local/lib/python3.6/dist-packages (from matplotlib) (0.10.0)\n", "Requirement already satisfied: kiwisolver>=1.0.1 in /usr/local/lib/python3.6/dist-packages (from matplotlib) (1.0.1)\n", "Requirement already satisfied: pyparsing!=2.0.4,!=2.1.2,!=2.1.6,>=2.0.1 in /usr/local/lib/python3.6/dist-packages (from matplotlib) (2.3.0)\n", "Requirement already satisfied: python-dateutil>=2.1 in /usr/local/lib/python3.6/dist-packages (from matplotlib) (2.7.3)\n", "Requirement already satisfied: llvmlite>=0.26.0dev0 in /usr/local/lib/python3.6/dist-packages (from numba) (0.26.0)\n", "Requirement already satisfied: prompt-toolkit<2.1.0,>=2.0.0 in /usr/local/lib/python3.6/dist-packages (from ipython->watermark) (2.0.7)\n", "Requirement already satisfied: jedi>=0.10 in /usr/local/lib/python3.6/dist-packages (from ipython->watermark) (0.12.1)\n", "Requirement already satisfied: pygments in /usr/local/lib/python3.6/dist-packages (from ipython->watermark) (2.2.0)\n", "Requirement already satisfied: decorator in /usr/local/lib/python3.6/dist-packages (from ipython->watermark) (4.3.0)\n", "Requirement already satisfied: setuptools>=18.5 in /usr/local/lib/python3.6/dist-packages (from ipython->watermark) (40.5.0)\n", "Requirement already satisfied: traitlets>=4.2 in /usr/local/lib/python3.6/dist-packages (from ipython->watermark) (4.3.2)\n", "Requirement already satisfied: pickleshare in /usr/local/lib/python3.6/dist-packages (from ipython->watermark) (0.7.5)\n", "Requirement already satisfied: pexpect; sys_platform != \"win32\" in /usr/local/lib/python3.6/dist-packages (from ipython->watermark) (4.6.0)\n", "Requirement already satisfied: backcall in /usr/local/lib/python3.6/dist-packages (from ipython->watermark) (0.1.0)\n", "Requirement already satisfied: six in /home/lilian/.local/lib/python3.6/site-packages (from cycler>=0.10->matplotlib) (1.11.0)\n", "Requirement already satisfied: wcwidth in /usr/local/lib/python3.6/dist-packages (from prompt-toolkit<2.1.0,>=2.0.0->ipython->watermark) (0.1.7)\n", "Requirement already satisfied: parso>=0.3.0 in /usr/local/lib/python3.6/dist-packages (from jedi>=0.10->ipython->watermark) (0.3.1)\n", "Requirement already satisfied: ipython-genutils in /usr/local/lib/python3.6/dist-packages (from traitlets>=4.2->ipython->watermark) (0.2.0)\n", "Requirement already satisfied: ptyprocess>=0.5 in /usr/local/lib/python3.6/dist-packages (from pexpect; sys_platform != \"win32\"->ipython->watermark) (0.6.0)\n", "The watermark extension is already loaded. To reload it, use:\n", " %reload_ext watermark\n", "Lilian Besson and Emilie Kaufmann \n", "\n", "CPython 3.6.7\n", "IPython 7.2.0\n", "\n", "numpy 1.15.4\n", "scipy 1.1.0\n", "matplotlib 3.0.2\n", "joblib 0.12.1\n", "numba 0.41.0\n", "cython 0.29.1\n", "\n", "compiler : GCC 8.2.0\n", "system : Linux\n", "release : 4.15.0-42-generic\n", "machine : x86_64\n", "processor : x86_64\n", "CPU cores : 4\n", "interpreter: 64bit\n" ] } ], "source": [ "!pip3 install watermark numpy scipy matplotlib joblib numba cython\n", "%load_ext watermark\n", "%watermark -v -m -p numpy,scipy,matplotlib,joblib,numba,cython -a \"Lilian Besson and Emilie Kaufmann\"" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "import numpy as np" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "try:\n", " from tqdm import tqdm_notebook as tqdm\n", "except:\n", " def tqdm(iterator, *args, **kwargs):\n", " return iterator" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Mathematical notations for stationary problems\n", "\n", "We consider $K \\geq 1$ arms, which are distributions $\\nu_k$.\n", "We focus on Bernoulli distributions, which are characterized by their means, $\\nu_k = \\mathcal{B}(\\mu_k)$ for $\\mu_k\\in[0,1]$.\n", "A stationary bandit problem is defined here by the vector $[\\mu_1,\\dots,\\mu_K]$.\n", "\n", "For a fixed problem and a *horizon* $T\\in\\mathbb{N}$, $T\\geq1$, we draw samples from the $K$ distributions to get *data*: $\\forall t, r_k(t) \\sim \\nu_k$, ie, $\\mathbb{P}(r_k(t) = 1) = \\mu_k$ and $r_k(t) \\in \\{0,1\\}$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Generating stationary data\n", "\n", "Here we give some examples of stationary problems and examples of data we can draw from them." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "def bernoulli_samples(means, horizon=1000):\n", " if np.size(means) == 1:\n", " return np.random.binomial(1, means, size=horizon)\n", " else:\n", " results = np.zeros((np.size(means), horizon))\n", " for i, mean in enumerate(means):\n", " results[i] = np.random.binomial(1, mean, size=horizon)\n", " return results" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 0, 0])" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "problem1 = [0.5]\n", "\n", "bernoulli_samples(problem1, horizon=20)" ] }, { "cell_type": "code", "execution_count": 27, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "61.5 µs ± 2.59 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n" ] } ], "source": [ "%timeit bernoulli_samples(problem1, horizon=1000)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now for Gaussian data:" ] }, { "cell_type": "code", "execution_count": 49, "metadata": {}, "outputs": [], "source": [ "sigma = 0.25 # Bernoulli are 1/4-sub Gaussian too!" ] }, { "cell_type": "code", "execution_count": 50, "metadata": {}, "outputs": [], "source": [ "def gaussian_samples(means, horizon=1000, sigma=sigma):\n", " if np.size(means) == 1:\n", " return np.random.normal(loc=means, scale=sigma, size=horizon)\n", " else:\n", " results = np.zeros((np.size(means), horizon))\n", " for i, mean in enumerate(means):\n", " results[i] = np.random.normal(loc=mean, scale=sigma, size=horizon)\n", " return results" ] }, { "cell_type": "code", "execution_count": 51, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([ 0.62387502, 0.25629644, 0.67019174, 0.17045626, 0.64599757,\n", " 0.37851622, 0.7365035 , 0.59688851, 0.63485378, 0.59576263,\n", " 0.27670489, 0.69860233, 0.49817269, 0.26672791, -0.05480486,\n", " 0.43916024, 0.4352945 , 0.630749 , -0.11451167, 0.11158323])" ] }, "execution_count": 51, "metadata": {}, "output_type": "execute_result" } ], "source": [ "gaussian_samples(problem1, horizon=20)" ] }, { "cell_type": "code", "execution_count": 52, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "75.8 µs ± 4.3 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n" ] } ], "source": [ "%timeit gaussian_samples(problem1, horizon=1000)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For bandit problem with $K \\geq 2$ arms, the *goal* is to design an online learning algorithm that roughly do the following:\n", "\n", "- For time $t=1$ to $t=T$ (unknown horizon)\n", " 1. Algorithm $A$ decide to draw arm $A(t) \\in\\{1,\\dots,K\\}$,\n", " 2. Get the reward $r(t) = r_{A(t)}(t) \\sim \\nu_{A(t)}$ from the (Bernoulli) distribution of that arm,\n", " 3. Give this observation of reward $r(t)$ coming from arm $A(t)$ to the algorithm,\n", " 4. Update internal state of the algorithm\n", "\n", "An algorithm is efficient if it obtains a high (expected) sum reward, ie, $\\sum_{t=1}^T r(t)$.\n", "\n", "Note that I don't focus on bandit algorithm here." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0.,\n", " 0., 0., 0., 0.],\n", " [1., 1., 0., 1., 0., 0., 1., 0., 1., 0., 1., 1., 1., 0., 0., 1.,\n", " 0., 0., 1., 0.],\n", " [1., 0., 1., 1., 1., 0., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,\n", " 0., 1., 1., 1.]])" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "problem2 = [0.1, 0.5, 0.9]\n", "\n", "bernoulli_samples(problem2, horizon=20)" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[ 0.46905853, 0.08312363, -0.23166553, 0.59750825, 0.23356612,\n", " -0.02132874, -0.26851831, 0.18202564, 0.37409142, -0.14397223,\n", " -0.22241626, -0.18888125, -0.13127815, 0.2314255 , 0.49433066,\n", " 0.39645143, -0.11529154, -0.13321403, 0.34457873, -0.16463127],\n", " [ 0.09801182, 0.8898005 , 0.76216447, 0.72837837, 0.25624864,\n", " 0.71225959, 0.56964125, 0.6309913 , 1.06315774, 0.28064027,\n", " 0.2251329 , 0.61036472, 0.43995415, 0.57338936, 0.74570767,\n", " 0.89390964, 0.50845914, 0.62240528, 0.27255829, 0.03454648],\n", " [ 1.04280525, 0.38307288, 0.81820526, 0.95947436, 1.21587846,\n", " 0.82568302, 0.60475185, 1.09188758, 1.40998237, 0.7291101 ,\n", " 0.94030203, 0.95715236, 0.6893489 , 1.26446965, 0.80080259,\n", " 0.96966515, 0.79117913, 0.48041062, 0.92280092, 1.27883349]])" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "problem2 = [0.1, 0.5, 0.9]\n", "\n", "gaussian_samples(problem2, horizon=20)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For instance on these data, the best arm is clearly the third one, with expected reward of $\\mu^* = \\max_k \\mu_k = 0.9$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Mathematical notations for piecewise stationary problems\n", "\n", "Now we fix the horizon $T\\in\\mathbb{N}$, $T\\geq1$ and we also consider a set of $\\Upsilon_T$ *break points*, $\\tau_1,\\dots,\\tau_{\\Upsilon_T} \\in\\{1,\\dots,T\\}$. We denote $\\tau_0 = 0$ and $\\tau_{\\Upsilon_T+1} = T$ for convenience of notations.\n", "We can assume that breakpoints are far \"enough\" from each other, for instance that there exists an integer $N\\in\\mathbb{N},N\\geq1$ such that $\\min_{i=0}^{\\Upsilon_T} \\tau_{i+1} - \\tau_i \\geq N K$. That is, on each *stationary interval*, a uniform sampling of the $K$ arms gives at least $N$ samples by arm.\n", "\n", "Now, in any stationary interval $[\\tau_i + 1, \\tau_{i+1}]$, the $K \\geq 1$ arms are distributions $\\nu_k^{(i)}$.\n", "We focus on Bernoulli distributions, which are characterized by their means, $\\nu_k^{(i)} := \\mathcal{B}(\\mu_k^{(i)})$ for $\\mu_k^{(i)}\\in[0,1]$.\n", "A piecewise stationary bandit problem is defined here by the vector $[\\mu_k^{(i)}]_{1\\leq k \\leq K, 1 \\leq i \\leq \\Upsilon_T}$.\n", "\n", "For a fixed problem and a *horizon* $T\\in\\mathbb{N}$, $T\\geq1$, we draw samples from the $K$ distributions to get *data*: $\\forall t, r_k(t) \\sim \\nu_k^{(i)}$ for $i$ the unique index of stationary interval such that $t\\in[\\tau_i + 1, \\tau_{i+1}]$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Generating fake piecewise stationary data\n", "\n", "The format to define piecewise stationary problem will be the following. It is compact but generic!\n", "\n", "The first example considers a unique arm, with 2 breakpoints uniformly spaced.\n", "- On the first interval, for instance from $t=1$ to $t=500$, that is $\\tau_1 = 500$, $\\mu_1^{(1)} = 0.1$,\n", "- On the second interval, for instance from $t=501$ to $t=1000$, that is $\\tau_2 = 100$, $\\mu_1^{(2)} = 0.5$,\n", "- On the third interval, for instance from $t=1001$ to $t=1500$, that $\\mu_1^{(3)} = 0.9$." ] }, { "cell_type": "code", "execution_count": 53, "metadata": {}, "outputs": [], "source": [ "# With 1 arm only!\n", "problem_piecewise_0 = lambda horizon: {\n", " \"listOfMeans\": [\n", " [0.1], # 0 to 499\n", " [0.5], # 500 to 999\n", " [0.8], # 1000 to 1499\n", " ],\n", " \"changePoints\": [\n", " int(0 * horizon / 1500.0),\n", " int(500 * horizon / 1500.0),\n", " int(1000 * horizon / 1500.0),\n", " ],\n", "}" ] }, { "cell_type": "code", "execution_count": 54, "metadata": { "code_folding": [ 1 ] }, "outputs": [], "source": [ "# With 2 arms\n", "problem_piecewise_1 = lambda horizon: {\n", " \"listOfMeans\": [\n", " [0.1, 0.2], # 0 to 399\n", " [0.1, 0.3], # 400 to 799\n", " [0.5, 0.3], # 800 to 1199\n", " [0.4, 0.3], # 1200 to 1599\n", " [0.3, 0.9], # 1600 to end\n", " ],\n", " \"changePoints\": [\n", " int(0 * horizon / 2000.0),\n", " int(400 * horizon / 2000.0),\n", " int(800 * horizon / 2000.0),\n", " int(1200 * horizon / 2000.0),\n", " int(1600 * horizon / 2000.0),\n", " ],\n", "}" ] }, { "cell_type": "code", "execution_count": 55, "metadata": { "code_folding": [ 1 ] }, "outputs": [], "source": [ "# With 3 arms\n", "problem_piecewise_2 = lambda horizon: {\n", " \"listOfMeans\": [\n", " [0.2, 0.5, 0.9], # 0 to 399\n", " [0.2, 0.2, 0.9], # 400 to 799\n", " [0.2, 0.2, 0.1], # 800 to 1199\n", " [0.7, 0.2, 0.1], # 1200 to 1599\n", " [0.7, 0.5, 0.1], # 1600 to end\n", " ],\n", " \"changePoints\": [\n", " int(0 * horizon / 2000.0),\n", " int(400 * horizon / 2000.0),\n", " int(800 * horizon / 2000.0),\n", " int(1200 * horizon / 2000.0),\n", " int(1600 * horizon / 2000.0),\n", " ],\n", "}" ] }, { "cell_type": "code", "execution_count": 56, "metadata": { "code_folding": [ 1 ] }, "outputs": [], "source": [ "# With 3 arms\n", "problem_piecewise_3 = lambda horizon: {\n", " \"listOfMeans\": [\n", " [0.4, 0.5, 0.9], # 0 to 399\n", " [0.5, 0.4, 0.7], # 400 to 799\n", " [0.6, 0.3, 0.5], # 800 to 1199\n", " [0.7, 0.2, 0.3], # 1200 to 1599\n", " [0.8, 0.1, 0.1], # 1600 to end\n", " ],\n", " \"changePoints\": [\n", " int(0 * horizon / 2000.0),\n", " int(400 * horizon / 2000.0),\n", " int(800 * horizon / 2000.0),\n", " int(1200 * horizon / 2000.0),\n", " int(1600 * horizon / 2000.0),\n", " ],\n", "}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we can write a utility function that transform this compact representation into a full list of means." ] }, { "cell_type": "code", "execution_count": 57, "metadata": { "code_folding": [] }, "outputs": [], "source": [ "def getFullHistoryOfMeans(problem, horizon=2000):\n", " \"\"\"Return the vector of mean of the arms, for a piece-wise stationary MAB.\n", "\n", " - It is a numpy array of shape (nbArms, horizon).\n", " \"\"\"\n", " pb = problem(horizon)\n", " listOfMeans, changePoints = pb['listOfMeans'], pb['changePoints']\n", " nbArms = len(listOfMeans[0])\n", " if horizon is None:\n", " horizon = np.max(changePoints)\n", " meansOfArms = np.ones((nbArms, horizon))\n", " for armId in range(nbArms):\n", " nbChangePoint = 0\n", " for t in range(horizon):\n", " if nbChangePoint < len(changePoints) - 1 and t >= changePoints[nbChangePoint + 1]:\n", " nbChangePoint += 1\n", " meansOfArms[armId][t] = listOfMeans[nbChangePoint][armId]\n", " return meansOfArms" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For examples :" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1,\n", " 0.1, 0.1, 0.1, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5,\n", " 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8,\n", " 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8]])" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "getFullHistoryOfMeans(problem_piecewise_0, horizon=50)" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1,\n", " 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5,\n", " 0.5, 0.5, 0.5, 0.5, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4,\n", " 0.4, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3],\n", " [0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.3, 0.3, 0.3,\n", " 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3,\n", " 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3,\n", " 0.3, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9]])" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "getFullHistoryOfMeans(problem_piecewise_1, horizon=50)" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2,\n", " 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2,\n", " 0.2, 0.2, 0.2, 0.2, 0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0.7,\n", " 0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0.7],\n", " [0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.2, 0.2, 0.2,\n", " 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2,\n", " 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2,\n", " 0.2, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5],\n", " [0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9,\n", " 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1,\n", " 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1,\n", " 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1]])" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "getFullHistoryOfMeans(problem_piecewise_2, horizon=50)" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.5, 0.5, 0.5,\n", " 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6,\n", " 0.6, 0.6, 0.6, 0.6, 0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0.7,\n", " 0.7, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8],\n", " [0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.4, 0.4, 0.4,\n", " 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3,\n", " 0.3, 0.3, 0.3, 0.3, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2,\n", " 0.2, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1],\n", " [0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.7, 0.7, 0.7,\n", " 0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5,\n", " 0.5, 0.5, 0.5, 0.5, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3,\n", " 0.3, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1]])" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "getFullHistoryOfMeans(problem_piecewise_3, horizon=50)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And now we need to be able to generate samples from such distributions." ] }, { "cell_type": "code", "execution_count": 58, "metadata": {}, "outputs": [], "source": [ "def piecewise_bernoulli_samples(problem, horizon=1000):\n", " fullMeans = getFullHistoryOfMeans(problem, horizon=horizon)\n", " nbArms, horizon = np.shape(fullMeans)\n", " results = np.zeros((nbArms, horizon))\n", " for i in range(nbArms):\n", " mean_i = fullMeans[i, :]\n", " for t in range(horizon):\n", " mean_i_t = max(0, min(1, mean_i[t])) # crop to [0, 1] !\n", " results[i, t] = np.random.binomial(1, mean_i_t)\n", " return results" ] }, { "cell_type": "code", "execution_count": 59, "metadata": {}, "outputs": [], "source": [ "def piecewise_gaussian_samples(problem, horizon=1000, sigma=sigma):\n", " fullMeans = getFullHistoryOfMeans(problem, horizon=horizon)\n", " nbArms, horizon = np.shape(fullMeans)\n", " results = np.zeros((nbArms, horizon))\n", " for i in range(nbArms):\n", " mean_i = fullMeans[i, :]\n", " for t in range(horizon):\n", " mean_i_t = mean_i[t]\n", " results[i, t] = np.random.normal(loc=mean_i_t, scale=sigma, size=1)\n", " return results" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Examples:" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1,\n", " 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1,\n", " 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5,\n", " 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5,\n", " 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5,\n", " 0.5, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8,\n", " 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8,\n", " 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8]])" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" }, { "data": { "text/plain": [ "array([[0., 1., 0., 1., 0., 1., 1., 0., 0., 0., 1., 0., 0., 0., 0., 0.,\n", " 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n", " 1., 1., 1., 1., 0., 0., 1., 0., 0., 0., 1., 0., 0., 1., 0., 1.,\n", " 0., 1., 1., 1., 1., 1., 1., 0., 1., 1., 1., 0., 0., 0., 1., 0.,\n", " 0., 1., 1., 1., 1., 1., 0., 1., 1., 1., 1., 1., 1., 1., 1., 1.,\n", " 1., 1., 1., 1., 1., 1., 1., 0., 1., 1., 1., 1., 1., 1., 1., 0.,\n", " 0., 1., 0., 1.]])" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "getFullHistoryOfMeans(problem_piecewise_0, horizon=100)\n", "piecewise_bernoulli_samples(problem_piecewise_0, horizon=100)" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[ 0.00315559, 0.2159303 , 0.01763935, -0.09514176, 0.07699997,\n", " 0.12997326, 0.36806044, 0.26968036, -0.29470607, -0.12074206,\n", " 0.38758995, 0.09229431, -0.1390067 , 0.01674725, 0.29248761,\n", " -0.00743648, -0.32053746, -0.05338394, -0.22795571, -0.07397954,\n", " 0.12089186, -0.05089131, 0.07538026, 0.37811067, 0.28527238,\n", " -0.44775343, 0.05399375, 0.56363279, 0.12052245, 0.04308359,\n", " 0.05259644, -0.09940602, 0.59624864, 0.47309089, 0.73371138,\n", " 0.71548802, 0.58391054, 0.57838614, 0.84475862, 0.17150277,\n", " 0.08094113, 0.12821665, 0.9622493 , 0.46279281, 0.69436121,\n", " 0.67516175, 0.63060937, 0.3242004 , 0.46987869, 0.87444974,\n", " 0.72351635, 0.13288009, 0.60558045, 0.42477301, 0.08474411,\n", " 0.28108802, 0.38570242, 0.16206467, 0.45767589, 0.62731769,\n", " 0.50540122, 0.09523314, 0.84005644, 0.29788404, 0.22493549,\n", " 0.43068172, 1.01212527, 0.66926627, 1.12914167, 0.56228928,\n", " 0.29419603, 0.83918851, 0.55743231, 1.2773794 , 0.47145187,\n", " 0.5581184 , 0.69071478, 0.78587582, 0.66516129, 0.13434453,\n", " 0.70156 , 1.00476689, 0.70888586, 0.61272042, 0.47476143,\n", " 0.52100125, 0.54323956, 0.84789493, 0.54069643, 1.35941397,\n", " 0.93167352, 0.85845656, 1.00137879, 0.97582926, 0.83642737,\n", " 0.63083642, 0.78315512, 0.92622805, 0.63902527, 0.48292036]])" ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ "piecewise_gaussian_samples(problem_piecewise_0, horizon=100)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We easily spot the (approximate) location of the breakpoint!\n", "\n", "Another example:" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1.,\n", " 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0.,\n", " 0., 0., 0., 0., 0., 0., 1., 1., 0., 0., 0., 0., 0., 1., 0., 1.,\n", " 0., 0., 0., 0., 1., 1., 1., 1., 1., 0., 0., 1., 1., 0., 0., 1.,\n", " 0., 1., 0., 1., 1., 0., 1., 0., 0., 0., 0., 1., 0., 0., 0., 0.,\n", " 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n", " 0., 1., 1., 1.],\n", " [0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 1., 1., 0., 0., 0.,\n", " 0., 0., 0., 0., 0., 0., 1., 1., 1., 0., 0., 0., 1., 0., 0., 1.,\n", " 1., 0., 1., 0., 1., 0., 0., 1., 0., 0., 1., 0., 1., 0., 0., 0.,\n", " 0., 0., 0., 0., 1., 0., 0., 1., 0., 0., 1., 0., 1., 0., 0., 0.,\n", " 1., 0., 1., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n", " 1., 1., 1., 1., 1., 1., 1., 0., 1., 1., 1., 1., 1., 1., 1., 1.,\n", " 1., 1., 1., 1.]])" ] }, "execution_count": 24, "metadata": {}, "output_type": "execute_result" } ], "source": [ "piecewise_bernoulli_samples(problem_piecewise_1, horizon=100)" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[-0.03440385, 0.35339867, 0.39487031, 0.03900667, -0.04061467,\n", " 0.42847871, 0.12854273, 0.11698151, 0.53333388, 0.55173192,\n", " 0.24042019, 0.6838909 , 0.64939741, 0.9333824 , 0.9194079 ,\n", " 0.49387316, 0.32736459, 0.2525435 , 0.37381781, 0.55302675],\n", " [ 0.28821344, 0.18003026, 0.25372301, 0.40054737, 0.55935287,\n", " -0.38328144, 0.06330063, 0.12570344, 0.46626937, 0.12882954,\n", " 0.35432014, 0.31273292, 0.90022341, 0.23341473, 0.0748465 ,\n", " 0.4811396 , 1.02582958, 0.99077802, 1.07355786, 0.48020466]])" ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ "piecewise_gaussian_samples(problem_piecewise_1, horizon=20)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "----\n", "# Python implementations of some statistical tests\n", "\n", "I will implement here the following statistical tests.\n", "I give a link to the implementation of the correspond bandit policy in my framework [`SMPyBandits`](https://smpybandits.github.io/)\n", "\n", "- Monitored (based on a McDiarmid inequality), for Monitored-UCB or [`M-UCB`](),\n", "- CUSUM, for [`CUSUM-UCB`](https://smpybandits.github.io/docs/Policies.CD_UCB.html?highlight=cusum#Policies.CD_UCB.CUSUM_IndexPolicy),\n", "- PHT, for [`PHT-UCB`](https://smpybandits.github.io/docs/Policies.CD_UCB.html?highlight=cusum#Policies.CD_UCB.PHT_IndexPolicy),\n", "- Gaussian GLR, for [`GaussianGLR-UCB`](https://smpybandits.github.io/docs/Policies.CD_UCB.html?highlight=glr#Policies.CD_UCB.GaussianGLR_IndexPolicy),\n", "- Bernoulli GLR, for [`BernoulliGLR-UCB`](https://smpybandits.github.io/docs/Policies.CD_UCB.html?highlight=glr#Policies.CD_UCB.BernoulliGLR_IndexPolicy)." ] }, { "cell_type": "code", "execution_count": 60, "metadata": {}, "outputs": [], "source": [ "class ChangePointDetector(object):\n", " def __init__(self, **kwargs):\n", " self._kwargs = kwargs\n", " for key, value in kwargs.items():\n", " setattr(self, key, value)\n", "\n", " def __str__(self):\n", " return f\"{self.__class__.__name__}{f'({repr(self._kwargs)})' if self._kwargs else ''}\"\n", "\n", " def detect(self, all_data, t):\n", " raise NotImplementedError" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Having classes is simply to be able to pretty print the algorithms when they have parameters:" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "ChangePointDetector\n" ] } ], "source": [ "print(ChangePointDetector())" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "ChangePointDetector({'w': 10, 'b': 1})\n" ] } ], "source": [ "print(ChangePointDetector(w=10, b=1))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## `Monitored`\n", "\n", "It uses a McDiarmid inequality. For a (pair) window size $w\\in\\mathbb{N}$ and a threshold $b\\in\\mathbb{R}^+$.\n", "At time $t$, if there is at least $w$ data in the data vector $(X_i)_i$, then let $Y$ denote the last $w$ data.\n", "A change is detected if\n", "$$ |\\sum_{i=w/2+1}^{w} Y_i - \\sum_{i=1}^{w/2} Y_i | > b ? $$" ] }, { "cell_type": "code", "execution_count": 64, "metadata": {}, "outputs": [], "source": [ "NB_ARMS = 1\n", "WINDOW_SIZE = 80" ] }, { "cell_type": "code", "execution_count": 65, "metadata": {}, "outputs": [], "source": [ "import numba" ] }, { "cell_type": "code", "execution_count": 88, "metadata": {}, "outputs": [], "source": [ "class Monitored(ChangePointDetector):\n", " def __init__(self, window_size=WINDOW_SIZE, threshold_b=None):\n", " super().__init__(window_size=window_size, threshold_b=threshold_b)\n", "\n", " def __str__(self):\n", " if self.threshold_b:\n", " return f\"Monitored($w={self.window_size:.3g}$, $b={self.threshold_b:.3g}$)\"\n", " else:\n", " latexname = r\"\\sqrt{\\frac{w}{2} \\log(2 T^2)}\"\n", " return f\"Monitored($w={self.window_size:.3g}$, $b={latexname}$)\"\n", "\n", " def detect(self, all_data, t):\n", " r\"\"\" A change is detected for the current arm if the following test is true:\n", "\n", " .. math:: |\\sum_{i=w/2+1}^{w} Y_i - \\sum_{i=1}^{w/2} Y_i | > b ?\n", "\n", " - where :math:`Y_i` is the i-th data in the latest w data from this arm (ie, :math:`X_k(t)` for :math:`t = n_k - w + 1` to :math:`t = n_k` current number of samples from arm k).\n", " - where :attr:`threshold_b` is the threshold b of the test, and :attr:`window_size` is the window-size w.\n", " \"\"\"\n", " data = all_data[:t]\n", " # don't try to detect change if there is not enough data!\n", " if len(data) < self.window_size:\n", " return False\n", "\n", " # compute parameters\n", " horizon = len(all_data)\n", " threshold_b = self.threshold_b\n", " if threshold_b is None:\n", " threshold_b = np.sqrt(self.window_size/2 * np.log(2 * NB_ARMS * horizon**2))\n", "\n", " last_w_data = data[-self.window_size:]\n", " sum_first_half = np.sum(last_w_data[:self.window_size//2])\n", " sum_second_half = np.sum(last_w_data[self.window_size//2:])\n", " return abs(sum_first_half - sum_second_half) > threshold_b" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## `CUSUM`\n", "\n", "The two-sided CUSUM algorithm, from [Page, 1954], works like this:\n", "\n", "- For each *data* k, compute:\n", "\n", "$$\n", "s_k^- = (y_k - \\hat{u}_0 - \\varepsilon) 1(k > M),\\\\\n", "s_k^+ = (\\hat{u}_0 - y_k - \\varepsilon) 1(k > M),\\\\\n", "g_k^+ = \\max(0, g_{k-1}^+ + s_k^+),\\\\\n", "g_k^- = \\max(0, g_{k-1}^- + s_k^-).\n", "$$\n", "\n", "- The change is detected if $\\max(g_k^+, g_k^-) > h$, where $h=$`threshold_h` is the threshold of the test,\n", "- And $\\hat{u}_0 = \\frac{1}{M} \\sum_{k=1}^{M} y_k$ is the mean of the first M samples, where M is `M` the min number of observation between change points." ] }, { "cell_type": "code", "execution_count": 67, "metadata": { "code_folding": [] }, "outputs": [], "source": [ "#: Precision of the test.\n", "EPSILON = 0.5\n", "\n", "#: Default value of :math:`\\lambda`.\n", "LAMBDA = 1\n", "\n", "#: Hypothesis on the speed of changes: between two change points, there is at least :math:`M * K` time steps, where K is the number of arms, and M is this constant.\n", "MIN_NUMBER_OF_OBSERVATION_BETWEEN_CHANGE_POINT = 50\n", "\n", "MAX_NB_RANDOM_EVENTS = 1" ] }, { "cell_type": "code", "execution_count": 68, "metadata": {}, "outputs": [], "source": [ "from scipy.special import comb" ] }, { "cell_type": "code", "execution_count": 89, "metadata": { "code_folding": [] }, "outputs": [], "source": [ "def compute_h__CUSUM(horizon, \n", " verbose=False,\n", " M=MIN_NUMBER_OF_OBSERVATION_BETWEEN_CHANGE_POINT,\n", " max_nb_random_events=MAX_NB_RANDOM_EVENTS,\n", " nbArms=1,\n", " epsilon=EPSILON,\n", " lmbda=LAMBDA,\n", " ):\n", " r\"\"\" Compute the values :math:`C_1^+, C_1^-, C_1, C_2, h` from the formulas in Theorem 2 and Corollary 2 in the paper.\"\"\"\n", " T = int(max(1, horizon))\n", " UpsilonT = int(max(1, max_nb_random_events))\n", " K = int(max(1, nbArms))\n", " C1_minus = np.log(((4 * epsilon) / (1-epsilon)**2) * comb(M, int(np.floor(2 * epsilon * M))) * (2 * epsilon)**M + 1)\n", " C1_plus = np.log(((4 * epsilon) / (1+epsilon)**2) * comb(M, int(np.ceil(2 * epsilon * M))) * (2 * epsilon)**M + 1)\n", " C1 = min(C1_minus, C1_plus)\n", " if C1 == 0: C1 = 1 # FIXME\n", " h = 1/C1 * np.log(T / UpsilonT)\n", " return h" ] }, { "cell_type": "code", "execution_count": 151, "metadata": { "code_folding": [] }, "outputs": [], "source": [ "class CUSUM(ChangePointDetector):\n", " def __init__(self,\n", " epsilon=EPSILON,\n", " M=MIN_NUMBER_OF_OBSERVATION_BETWEEN_CHANGE_POINT,\n", " threshold_h=None,\n", " ):\n", " assert 0 < epsilon < 1, f\"Error: epsilon for CUSUM must be in (0, 1) but is {epsilon}.\"\n", " super().__init__(epsilon=epsilon, M=M, threshold_h=threshold_h)\n", " \n", " def __str__(self):\n", " if self.threshold_h:\n", " return fr\"CUSUM($\\varepsilon={self.epsilon:.3g}$, $M={self.M}$, $h={self.threshold_h:.3g}$)\"\n", " else:\n", " return fr\"CUSUM($\\varepsilon={self.epsilon:.3g}$, $M={self.M}$, $h=$'auto')\"\n", "\n", " def detect(self, all_data, t):\n", " r\"\"\" Detect a change in the current arm, using the two-sided CUSUM algorithm [Page, 1954].\n", "\n", " - For each *data* k, compute:\n", "\n", " .. math::\n", "\n", " s_k^- &= (y_k - \\hat{u}_0 - \\varepsilon) 1(k > M),\\\\\n", " s_k^+ &= (\\hat{u}_0 - y_k - \\varepsilon) 1(k > M),\\\\\n", " g_k^+ &= \\max(0, g_{k-1}^+ + s_k^+),\\\\\n", " g_k^- &= \\max(0, g_{k-1}^- + s_k^-).\n", "\n", " - The change is detected if :math:`\\max(g_k^+, g_k^-) > h`, where :attr:`threshold_h` is the threshold of the test,\n", " - And :math:`\\hat{u}_0 = \\frac{1}{M} \\sum_{k=1}^{M} y_k` is the mean of the first M samples, where M is :attr:`M` the min number of observation between change points.\n", " \"\"\"\n", " data = all_data[:t]\n", "\n", " # compute parameters\n", " horizon = len(all_data)\n", " threshold_h = self.threshold_h\n", " if self.threshold_h is None:\n", " threshold_h = compute_h__CUSUM(horizon, self.M, 1, epsilon=self.epsilon)\n", "\n", " gp, gm = 0, 0\n", " # First we use the first M samples to calculate the average :math:`\\hat{u_0}`.\n", " u0hat = np.mean(data[:self.M])\n", " for k in range(self.M + 1, len(data)):\n", " y_k = data[k]\n", " sp = u0hat - y_k - self.epsilon # no need to multiply by (k > self.M)\n", " sm = y_k - u0hat - self.epsilon # no need to multiply by (k > self.M)\n", " gp = max(0, gp + sp)\n", " gm = max(0, gm + sm)\n", " if max(gp, gm) >= threshold_h:\n", " return True\n", " return False" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## `PHT`\n", "\n", "The two-sided CUSUM algorithm, from [Hinkley, 1971], works like this:\n", "\n", "- For each *data* k, compute:\n", "\n", "$$\n", "s_k^- = y_k - \\hat{y}_k - \\varepsilon,\\\\\n", "s_k^+ = \\hat{y}_k - y_k - \\varepsilon,\\\\\n", "g_k^+ = \\max(0, g_{k-1}^+ + s_k^+),\\\\\n", "g_k^- = \\max(0, g_{k-1}^- + s_k^-).\n", "$$\n", "\n", "- The change is detected if $\\max(g_k^+, g_k^-) > h$, where $h=$`threshold_h` is the threshold of the test,\n", "- And $\\hat{y}_k = \\frac{1}{k} \\sum_{s=1}^{k} y_s$ is the mean of the first k samples." ] }, { "cell_type": "code", "execution_count": 152, "metadata": {}, "outputs": [], "source": [ "class PHT(ChangePointDetector):\n", " def __init__(self,\n", " epsilon=EPSILON,\n", " M=MIN_NUMBER_OF_OBSERVATION_BETWEEN_CHANGE_POINT,\n", " threshold_h=None,\n", " ):\n", " assert 0 < epsilon < 1, f\"Error: epsilon for CUSUM must be in (0, 1) but is {epsilon}.\"\n", " super().__init__(epsilon=epsilon, M=M, threshold_h=threshold_h)\n", " \n", " def __str__(self):\n", " if self.threshold_h:\n", " return fr\"PHT($\\varepsilon={self.epsilon:.3g}$, $M={self.M}$, $h={self.threshold_h:.3g}$)\"\n", " else:\n", " return fr\"PHT($\\varepsilon={self.epsilon:.3g}$, $M={self.M}$, $h=$'auto')\"\n", "\n", " def detect(self, all_data, t):\n", " r\"\"\" Detect a change in the current arm, using the two-sided PHT algorithm [Hinkley, 1971].\n", "\n", " - For each *data* k, compute:\n", "\n", " .. math::\n", "\n", " s_k^- &= y_k - \\hat{y}_k - \\varepsilon,\\\\\n", " s_k^+ &= \\hat{y}_k - y_k - \\varepsilon,\\\\\n", " g_k^+ &= \\max(0, g_{k-1}^+ + s_k^+),\\\\\n", " g_k^- &= \\max(0, g_{k-1}^- + s_k^-).\n", "\n", " - The change is detected if :math:`\\max(g_k^+, g_k^-) > h`, where :attr:`threshold_h` is the threshold of the test,\n", " - And :math:`\\hat{y}_k = \\frac{1}{k} \\sum_{s=1}^{k} y_s` is the mean of the first k samples.\n", " \"\"\"\n", " data = all_data[:t]\n", "\n", " # compute parameters\n", " horizon = len(all_data)\n", " threshold_h = self.threshold_h\n", " if threshold_h is None:\n", " threshold_h = compute_h__CUSUM(horizon, self.M, 1, epsilon=self.epsilon)\n", "\n", " gp, gm = 0, 0\n", " y_k_hat = 0\n", " # First we use the first M samples to calculate the average :math:`\\hat{u_0}`.\n", " for k, y_k in enumerate(data):\n", " y_k_hat = (k * y_k_hat + y_k) / (k + 1) # XXX smart formula to update the mean!\n", " sp = y_k_hat - y_k - self.epsilon\n", " sm = y_k - y_k_hat - self.epsilon\n", " gp = max(0, gp + sp)\n", " gm = max(0, gm + sm)\n", " if max(gp, gm) >= threshold_h:\n", " return True\n", " return False" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "## `Gaussian GLR`\n", "\n", "The Generalized Likelihood Ratio test (GLR) works with a one-dimensional exponential family, for which we have a function `kl` such that if $\\mu_1,\\mu_2$ are the means of two distributions $\\nu_1,\\nu_2$, then $\\mathrm{KL}(\\mathcal{D}(\\nu_1), \\mathcal{D}(\\nu_1))=$ `kl` $(\\mu_1,\\mu_2)$.\n", "\n", "- For each *time step* $s$ between $t_0=0$ and $t$, compute:\n", "$$G^{\\mathcal{N}_1}_{t_0:s:t} = (s-t_0+1) \\mathrm{kl}(\\mu_{t_0,s}, \\mu_{t_0,t}) + (t-s) \\mathrm{kl}(\\mu_{s+1,t}, \\mu_{t_0,t}).$$\n", "\n", "- The change is detected if there is a time $s$ such that $G^{\\mathcal{N}_1}_{t_0:s:t} > b(t_0, s, t, \\delta)$, where $b(t_0, s, t, \\delta)=$ `threshold_h` is the threshold of the test,\n", "- And $\\mu_{a,b} = \\frac{1}{b-a+1} \\sum_{s=a}^{b} y_s$ is the mean of the samples between $a$ and $b$.\n", "\n", "The threshold is computed as:\n", "$$ b(t_0, s, t, \\delta):= \\left(1 + \\frac{1}{t - t_0 + 1}\\right) 2 \\log\\left(\\frac{2 (t - t_0) \\sqrt{(t - t_0) + 2}}{\\delta}\\right).$$\n", "\n", "Another threshold we want to check is the following:\n", "$$ b(t_0, s, t, \\delta):= \\log\\left(\\frac{(s - t_0 + 1) (t - s)}{\\delta}\\right).$$" ] }, { "cell_type": "code", "execution_count": 124, "metadata": { "code_folding": [] }, "outputs": [], "source": [ "from math import log, isinf\n", "\n", "def compute_c__GLR_0(t0, s, t, horizon=None, delta=None):\n", " r\"\"\" Compute the values :math:`c` from the corollary of of Theorem 2 from [\"Sequential change-point detection: Laplace concentration of scan statistics and non-asymptotic delay bounds\", O.-A. Maillard, 2018].\n", "\n", " - The threshold is computed as:\n", "\n", " .. math:: h := \\left(1 + \\frac{1}{t - t_0 + 1}\\right) 2 \\log\\left(\\frac{2 (t - t_0) \\sqrt{(t - t_0) + 2}}{\\delta}\\right).\n", " \"\"\"\n", " if delta is None:\n", " T = int(max(1, horizon))\n", " delta = 1.0 / T\n", " t_m_t0 = abs(t - t0)\n", " c = (1 + (1 / (t_m_t0 + 1.0))) * 2 * log((2 * t_m_t0 * np.sqrt(t_m_t0 + 2)) / delta)\n", " if c < 0 and isinf(c): c = float('+inf')\n", " return c" ] }, { "cell_type": "code", "execution_count": 125, "metadata": { "code_folding": [] }, "outputs": [], "source": [ "from math import log, isinf\n", "\n", "def compute_c__GLR(t0, s, t, horizon=None, delta=None):\n", " r\"\"\" Compute the values :math:`c` from the corollary of of Theorem 2 from [\"Sequential change-point detection: Laplace concentration of scan statistics and non-asymptotic delay bounds\", O.-A. Maillard, 2018].\n", "\n", " - The threshold is computed as:\n", "\n", " .. math:: h := \\log\\left(\\frac{(s - t_0 + 1) (t - s)}{\\delta}\\right).\n", " \"\"\"\n", " if delta is None: \n", " T = int(max(1, horizon))\n", " delta = 1.0 / T\n", " arg = (s - t0 + 1) * (t - s) / delta\n", " if arg <= 0: c = float('+inf')\n", " else: c = log(arg)\n", " return c" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For Gaussian distributions of known variance, the Kullback-Leibler divergence is easy to compute:\n", "\n", "Kullback-Leibler divergence for Gaussian distributions of means $x$ and $y$ and variances $\\sigma^2_x = \\sigma^2_y$, $\\nu_1 = \\mathcal{N}(x, \\sigma_x^2)$ and $\\nu_2 = \\mathcal{N}(y, \\sigma_x^2)$ is:\n", "\n", "$$\\mathrm{KL}(\\nu_1, \\nu_2) = \\frac{(x - y)^2}{2 \\sigma_y^2} + \\frac{1}{2}\\left( \\frac{\\sigma_x^2}{\\sigma_y^2} - 1 \\log\\left(\\frac{\\sigma_x^2}{\\sigma_y^2}\\right) \\right).$$" ] }, { "cell_type": "code", "execution_count": 126, "metadata": { "code_folding": [ 0 ] }, "outputs": [], "source": [ "def klGauss(x, y, sig2x=0.25):\n", " r\"\"\" Kullback-Leibler divergence for Gaussian distributions of means ``x`` and ``y`` and variances ``sig2x`` and ``sig2y``, :math:`\\nu_1 = \\mathcal{N}(x, \\sigma_x^2)` and :math:`\\nu_2 = \\mathcal{N}(y, \\sigma_x^2)`:\n", "\n", " .. math:: \\mathrm{KL}(\\nu_1, \\nu_2) = \\frac{(x - y)^2}{2 \\sigma_y^2} + \\frac{1}{2}\\left( \\frac{\\sigma_x^2}{\\sigma_y^2} - 1 \\log\\left(\\frac{\\sigma_x^2}{\\sigma_y^2}\\right) \\right).\n", "\n", " See https://en.wikipedia.org/wiki/Normal_distribution#Other_properties\n", "\n", " - sig2y = sig2x (same variance).\n", " \"\"\"\n", " return (x - y) ** 2 / (2. * sig2x)" ] }, { "cell_type": "code", "execution_count": 127, "metadata": { "code_folding": [ 3 ] }, "outputs": [], "source": [ "class GaussianGLR(ChangePointDetector):\n", " def __init__(self, mult_threshold_h=1, delta=None):\n", " super().__init__(mult_threshold_h=mult_threshold_h, delta=delta)\n", " \n", " def __str__(self):\n", " return r\"Gaussian-GLR($h_0={}$, $\\delta={}$)\".format(\n", " f\"{self.mult_threshold_h:.3g}\" if self.mult_threshold_h is not None else 'auto',\n", " f\"{self.delta:.3g}\" if self.delta is not None else 'auto',\n", " )\n", "\n", " def detect(self, all_data, t):\n", " r\"\"\" Detect a change in the current arm, using the Generalized Likelihood Ratio test (GLR) and the :attr:`kl` function.\n", "\n", " - For each *time step* :math:`s` between :math:`t_0=0` and :math:`t`, compute:\n", "\n", " .. math::\n", "\n", " G^{\\mathcal{N}_1}_{t_0:s:t} = (s-t_0+1)(t-s) \\mathrm{kl}(\\mu_{s+1,t}, \\mu_{t_0,s}) / (t-t_0+1).\n", "\n", " - The change is detected if there is a time :math:`s` such that :math:`G^{\\mathcal{N}_1}_{t_0:s:t} > h`, where :attr:`threshold_h` is the threshold of the test,\n", " - And :math:`\\mu_{a,b} = \\frac{1}{b-a+1} \\sum_{s=a}^{b} y_s` is the mean of the samples between :math:`a` and :math:`b`.\n", " \"\"\"\n", " data = all_data[:t]\n", " t0 = 0\n", " horizon = len(all_data)\n", "\n", " # compute parameters\n", " mean_all = np.mean(data[t0 : t+1])\n", " mean_before = 0\n", " mean_after = mean_all\n", " for s in range(t0, t):\n", " # DONE okay this is efficient we don't compute the same means too many times!\n", " y = data[s]\n", " mean_before = (s * mean_before + y) / (s + 1)\n", " mean_after = ((t + 1 - s + t0) * mean_after - y) / (t - s + t0)\n", " kl_before = klGauss(mean_before, mean_all)\n", " kl_after = klGauss(mean_after, mean_all)\n", " threshold_h = self.mult_threshold_h * compute_c__GLR(t0, s, t, horizon=horizon, delta=self.delta)\n", " glr = (s - t0 + 1) * kl_before + (t - s) * kl_after\n", " if glr >= threshold_h:\n", " return True\n", " return False" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## `Bernoulli GLR`\n", "\n", "The same GLR algorithm but using the Bernoulli KL, given by:\n", "\n", "$$\\mathrm{KL}(\\mathcal{B}(x), \\mathcal{B}(y)) = x \\log(\\frac{x}{y}) + (1-x) \\log(\\frac{1-x}{1-y}).$$" ] }, { "cell_type": "code", "execution_count": 257, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The cython extension is already loaded. To reload it, use:\n", " %reload_ext cython\n" ] } ], "source": [ "import cython\n", "%load_ext cython" ] }, { "cell_type": "code", "execution_count": 265, "metadata": {}, "outputs": [], "source": [ "def klBern(x: float, y: float) -> float:\n", " r\"\"\" Kullback-Leibler divergence for Bernoulli distributions. https://en.wikipedia.org/wiki/Bernoulli_distribution#Kullback.E2.80.93Leibler_divergence\n", "\n", " .. math:: \\mathrm{KL}(\\mathcal{B}(x), \\mathcal{B}(y)) = x \\log(\\frac{x}{y}) + (1-x) \\log(\\frac{1-x}{1-y}).\"\"\"\n", " x = min(max(x, 1e-6), 1 - 1e-6)\n", " y = min(max(y, 1e-6), 1 - 1e-6)\n", " return x * log(x / y) + (1 - x) * log((1 - x) / (1 - y))" ] }, { "cell_type": "code", "execution_count": 259, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1.91 µs ± 107 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)\n" ] } ], "source": [ "%timeit klBern(np.random.random(), np.random.random())" ] }, { "cell_type": "code", "execution_count": 260, "metadata": { "code_folding": [ 2 ] }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\n", " \n", "Generated by Cython 0.29.1
\n", "\n",
" Yellow lines hint at Python interaction.
\n",
" Click on a line that starts with a \"+
\" to see the C code that Cython generated for it.\n",
"
01: from libc.math cimport log\n", "
+02: eps = 1e-6 #: Threshold value: everything in [0, 1] is truncated to [eps, 1 - eps]\n", "
if (PyDict_SetItem(__pyx_d, __pyx_n_s_eps, __pyx_float_1eneg_6) < 0) __PYX_ERR(0, 2, __pyx_L1_error)\n", "
03:
\n",
"+04: def klBern(float x, float y) -> float:\n", "
/* Python wrapper */\n", "static PyObject *__pyx_pw_46_cython_magic_e1f92aa869fc57ff4f0ff63e582d134a_1klBern(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds); /*proto*/\n", "static char __pyx_doc_46_cython_magic_e1f92aa869fc57ff4f0ff63e582d134a_klBern[] = \" Kullback-Leibler divergence for Bernoulli distributions. https://en.wikipedia.org/wiki/Bernoulli_distribution#Kullback.E2.80.93Leibler_divergence\\n\\n .. math:: \\\\mathrm{KL}(\\\\mathcal{B}(x), \\\\mathcal{B}(y)) = x \\\\log(\\\\frac{x}{y}) + (1-x) \\\\log(\\\\frac{1-x}{1-y}).\";\n", "static PyMethodDef __pyx_mdef_46_cython_magic_e1f92aa869fc57ff4f0ff63e582d134a_1klBern = {\"klBern\", (PyCFunction)(void*)(PyCFunctionWithKeywords)__pyx_pw_46_cython_magic_e1f92aa869fc57ff4f0ff63e582d134a_1klBern, METH_VARARGS|METH_KEYWORDS, __pyx_doc_46_cython_magic_e1f92aa869fc57ff4f0ff63e582d134a_klBern};\n", "static PyObject *__pyx_pw_46_cython_magic_e1f92aa869fc57ff4f0ff63e582d134a_1klBern(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds) {\n", " float __pyx_v_x;\n", " float __pyx_v_y;\n", " PyObject *__pyx_r = 0;\n", " __Pyx_RefNannyDeclarations\n", " __Pyx_RefNannySetupContext(\"klBern (wrapper)\", 0);\n", " {\n", " static PyObject **__pyx_pyargnames[] = {&__pyx_n_s_x,&__pyx_n_s_y,0};\n", " PyObject* values[2] = {0,0};\n", " if (unlikely(__pyx_kwds)) {\n", " Py_ssize_t kw_args;\n", " const Py_ssize_t pos_args = PyTuple_GET_SIZE(__pyx_args);\n", " switch (pos_args) {\n", " case 2: values[1] = PyTuple_GET_ITEM(__pyx_args, 1);\n", " CYTHON_FALLTHROUGH;\n", " case 1: values[0] = PyTuple_GET_ITEM(__pyx_args, 0);\n", " CYTHON_FALLTHROUGH;\n", " case 0: break;\n", " default: goto __pyx_L5_argtuple_error;\n", " }\n", " kw_args = PyDict_Size(__pyx_kwds);\n", " switch (pos_args) {\n", " case 0:\n", " if (likely((values[0] = __Pyx_PyDict_GetItemStr(__pyx_kwds, __pyx_n_s_x)) != 0)) kw_args--;\n", " else goto __pyx_L5_argtuple_error;\n", " CYTHON_FALLTHROUGH;\n", " case 1:\n", " if (likely((values[1] = __Pyx_PyDict_GetItemStr(__pyx_kwds, __pyx_n_s_y)) != 0)) kw_args--;\n", " else {\n", " __Pyx_RaiseArgtupleInvalid(\"klBern\", 1, 2, 2, 1); __PYX_ERR(0, 4, __pyx_L3_error)\n", " }\n", " }\n", " if (unlikely(kw_args > 0)) {\n", " if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, \"klBern\") < 0)) __PYX_ERR(0, 4, __pyx_L3_error)\n", " }\n", " } else if (PyTuple_GET_SIZE(__pyx_args) != 2) {\n", " goto __pyx_L5_argtuple_error;\n", " } else {\n", " values[0] = PyTuple_GET_ITEM(__pyx_args, 0);\n", " values[1] = PyTuple_GET_ITEM(__pyx_args, 1);\n", " }\n", " __pyx_v_x = __pyx_PyFloat_AsFloat(values[0]); if (unlikely((__pyx_v_x == (float)-1) && PyErr_Occurred())) __PYX_ERR(0, 4, __pyx_L3_error)\n", " __pyx_v_y = __pyx_PyFloat_AsFloat(values[1]); if (unlikely((__pyx_v_y == (float)-1) && PyErr_Occurred())) __PYX_ERR(0, 4, __pyx_L3_error)\n", " }\n", " goto __pyx_L4_argument_unpacking_done;\n", " __pyx_L5_argtuple_error:;\n", " __Pyx_RaiseArgtupleInvalid(\"klBern\", 1, 2, 2, PyTuple_GET_SIZE(__pyx_args)); __PYX_ERR(0, 4, __pyx_L3_error)\n", " __pyx_L3_error:;\n", " __Pyx_AddTraceback(\"_cython_magic_e1f92aa869fc57ff4f0ff63e582d134a.klBern\", __pyx_clineno, __pyx_lineno, __pyx_filename);\n", " __Pyx_RefNannyFinishContext();\n", " return NULL;\n", " __pyx_L4_argument_unpacking_done:;\n", " __pyx_r = __pyx_pf_46_cython_magic_e1f92aa869fc57ff4f0ff63e582d134a_klBern(__pyx_self, __pyx_v_x, __pyx_v_y);\n", "\n", " /* function exit code */\n", " __Pyx_RefNannyFinishContext();\n", " return __pyx_r;\n", "}\n", "\n", "static PyObject *__pyx_pf_46_cython_magic_e1f92aa869fc57ff4f0ff63e582d134a_klBern(CYTHON_UNUSED PyObject *__pyx_self, float __pyx_v_x, float __pyx_v_y) {\n", " PyObject *__pyx_r = NULL;\n", " __Pyx_RefNannyDeclarations\n", " __Pyx_RefNannySetupContext(\"klBern\", 0);\n", "/* … */\n", " /* function exit code */\n", " __pyx_L1_error:;\n", " __Pyx_XDECREF(__pyx_t_6);\n", " __Pyx_AddTraceback(\"_cython_magic_e1f92aa869fc57ff4f0ff63e582d134a.klBern\", __pyx_clineno, __pyx_lineno, __pyx_filename);\n", " __pyx_r = NULL;\n", " __pyx_L0:;\n", " __Pyx_XGIVEREF(__pyx_r);\n", " __Pyx_RefNannyFinishContext();\n", " return __pyx_r;\n", "}\n", "/* … */\n", " __pyx_tuple_ = PyTuple_Pack(2, __pyx_n_s_x, __pyx_n_s_y); if (unlikely(!__pyx_tuple_)) __PYX_ERR(0, 4, __pyx_L1_error)\n", " __Pyx_GOTREF(__pyx_tuple_);\n", " __Pyx_GIVEREF(__pyx_tuple_);\n", "/* … */\n", " __pyx_t_1 = PyCFunction_NewEx(&__pyx_mdef_46_cython_magic_e1f92aa869fc57ff4f0ff63e582d134a_1klBern, NULL, __pyx_n_s_cython_magic_e1f92aa869fc57ff4f); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 4, __pyx_L1_error)\n", " __Pyx_GOTREF(__pyx_t_1);\n", " if (PyDict_SetItem(__pyx_d, __pyx_n_s_klBern, __pyx_t_1) < 0) __PYX_ERR(0, 4, __pyx_L1_error)\n", " __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;\n", "
05: r""" Kullback-Leibler divergence for Bernoulli distributions. https://en.wikipedia.org/wiki/Bernoulli_distribution#Kullback.E2.80.93Leibler_divergence\n", "
06:
\n",
"07: .. math:: \\mathrm{KL}(\\mathcal{B}(x), \\mathcal{B}(y)) = x \\log(\\frac{x}{y}) + (1-x) \\log(\\frac{1-x}{1-y})."""\n", "
+08: x = min(max(x, 1e-6), 1 - 1e-6)\n", "
__pyx_t_1 = (1.0 - 1e-6);\n", " __pyx_t_2 = 1e-6;\n", " __pyx_t_3 = __pyx_v_x;\n", " if (((__pyx_t_2 > __pyx_t_3) != 0)) {\n", " __pyx_t_4 = __pyx_t_2;\n", " } else {\n", " __pyx_t_4 = __pyx_t_3;\n", " }\n", " __pyx_t_2 = __pyx_t_4;\n", " if (((__pyx_t_1 < __pyx_t_2) != 0)) {\n", " __pyx_t_4 = __pyx_t_1;\n", " } else {\n", " __pyx_t_4 = __pyx_t_2;\n", " }\n", " __pyx_v_x = __pyx_t_4;\n", "
+09: y = min(max(y, 1e-6), 1 - 1e-6)\n", "
__pyx_t_4 = (1.0 - 1e-6);\n", " __pyx_t_1 = 1e-6;\n", " __pyx_t_3 = __pyx_v_y;\n", " if (((__pyx_t_1 > __pyx_t_3) != 0)) {\n", " __pyx_t_2 = __pyx_t_1;\n", " } else {\n", " __pyx_t_2 = __pyx_t_3;\n", " }\n", " __pyx_t_1 = __pyx_t_2;\n", " if (((__pyx_t_4 < __pyx_t_1) != 0)) {\n", " __pyx_t_2 = __pyx_t_4;\n", " } else {\n", " __pyx_t_2 = __pyx_t_1;\n", " }\n", " __pyx_v_y = __pyx_t_2;\n", "
+10: return x * log(x / y) + (1 - x) * log((1 - x) / (1 - y))\n", "
__Pyx_XDECREF(__pyx_r);\n", " if (unlikely(__pyx_v_y == 0)) {\n", " PyErr_SetString(PyExc_ZeroDivisionError, \"float division\");\n", " __PYX_ERR(0, 10, __pyx_L1_error)\n", " }\n", " __pyx_t_3 = (1.0 - __pyx_v_x);\n", " __pyx_t_5 = (1.0 - __pyx_v_y);\n", " if (unlikely(__pyx_t_5 == 0)) {\n", " PyErr_SetString(PyExc_ZeroDivisionError, \"float division\");\n", " __PYX_ERR(0, 10, __pyx_L1_error)\n", " }\n", " __pyx_t_6 = PyFloat_FromDouble(((__pyx_v_x * log((__pyx_v_x / __pyx_v_y))) + ((1.0 - __pyx_v_x) * log((__pyx_t_3 / __pyx_t_5))))); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 10, __pyx_L1_error)\n", " __Pyx_GOTREF(__pyx_t_6);\n", " __pyx_r = __pyx_t_6;\n", " __pyx_t_6 = 0;\n", " goto __pyx_L0;\n", "