kalman-folding-in-clojure.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "# COMPOSABLE STATISTICS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some sample scalar data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 195,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "#'user/zs"
      ]
     },
     "execution_count": 195,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(def zs [-0.178654, 0.828305, 0.0592247, -0.0121089, -1.48014, \n",
    "         -0.315044, -0.324796, -0.676357, 0.16301, -0.858164])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "TODO: generate new random data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "TODO: incanter integration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## RUNNING COUNT"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The traditional and obvious way:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 196,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10"
      ]
     },
     "execution_count": 196,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(reduce (fn [c z] (inc c)) \n",
    "        0 zs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... with all intermediate results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 197,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0 1 2 3 4 5 6 7 8 9 10)"
      ]
     },
     "execution_count": 197,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(reductions (fn [c z] (inc c)) 0 zs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### THREAD-SAFE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A thread-safe way, overkill for now, but can be called on multiple threads reading from asynchronous streams. It also shows the _let-over-lambda_ (LOL) style: closing over stateful variables and using thread-safe atomic updates. This is sematically equivalent to an OOP style with mutual exclusion and transactional updates, but shorter and easier to verify."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 198,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10"
      ]
     },
     "execution_count": 198,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(let [initial-count 0]\n",
    "    (reduce ; Let-over-lambda follows\n",
    "        (let [running-count (atom initial-count)]\n",
    "            (fn [c z] ; closes over the atom\n",
    "                (swap! running-count inc)\n",
    "                @running-count)) ; ~~> new value for c\n",
    "        initial-count\n",
    "        zs))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Showing all intermediate results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 200,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0 1 2 3 4 5 6 7 8 9 10)"
      ]
     },
     "execution_count": 200,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(let [initial-count 0]\n",
    "    (reductions ; <-- this is the only difference to above\n",
    "        (let [running-count (atom initial-count)]\n",
    "            (fn [c z]\n",
    "                (swap! running-count inc)\n",
    "                @running-count)) ; ~~> new value for c\n",
    "        initial-count\n",
    "        zs))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### AVOIDING _REDUCE_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A thread-safe way, without `reduce`. This reduces input coupling to the environment by absorbing `initial-count` into the LOL. The environment only provides data. In this case, the environment provides data through mapping over a sequence. Below, the environment provides data through  asynchronous streams, but the LOL stays the same.\n",
    "\n",
    "The output is still coupled to the environment through `println`. We get rid of that, too, below. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 240,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 2 3 4 5 6 7 8 9 10 "
     ]
    }
   ],
   "source": [
    "(dorun \n",
    "    (map \n",
    "        (let [running-count (atom 0)]\n",
    "            (fn [z]\n",
    "                (swap! running-count inc)\n",
    "                (print (str @running-count \" \"))))\n",
    "        zs))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## RUNNING MEAN"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is a beautiful expression for the running mean, $x$. The new running mean is the old running mean plus a correction. The correction depends on new inputs $z$ but only on old statistics. That's very nice because everything except $z$ on the right-hand side of equations is prior information. The correction is a _gain_ times the _residual_. The residual is the difference betwween the new measurement $z$ and the old mean $x$. The gain is $1/(n+1)$. $n$ is a statistic, so it is an _old_ value, too: the count of observations / measurements / $z$ values seen heretofore, that is, not including the current measurement $z$. \n",
    "\n",
    "We can write the recurrence elegantly without subscripts as follows: $x\\leftarrow{x+K\\,(z-x)}$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 241,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{:count 1, :mean -0.178654}\n",
      "{:count 2, :mean 0.3248255}\n",
      "{:count 3, :mean 0.2362919}\n",
      "{:count 4, :mean 0.1741917}\n",
      "{:count 5, :mean -0.15667464000000003}\n",
      "{:count 6, :mean -0.18306953333333337}\n",
      "{:count 7, :mean -0.20331617142857145}\n",
      "{:count 8, :mean -0.262446275}\n",
      "{:count 9, :mean -0.21517335555555556}\n",
      "{:count 10, :mean -0.27947242}\n"
     ]
    }
   ],
   "source": [
    "(dorun \n",
    "    (map\n",
    "        (let [running-stats (atom {:count 0, :mean 0})]\n",
    "            (fn [z]\n",
    "                (let [{x :mean, n :count} @running-stats\n",
    "                      n+1 (inc n)\n",
    "                      K   (/ 1.0 n+1)]\n",
    "                    (swap! running-stats conj\n",
    "                           [:count n+1]\n",
    "                           [:mean (+ x (* K (- z x)))]))\n",
    "                (println @running-stats)))\n",
    "        zs))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### REMOVING OUTPUT COUPLING"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can remove `println` from inside the function and lift it into a parameter. Now the function is completely decoupled from its environment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 242,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "({:count 1, :mean -0.178654}\n",
      " {:count 2, :mean 0.3248255}\n",
      " {:count 3, :mean 0.2362919}\n",
      " {:count 4, :mean 0.1741917}\n",
      " {:count 5, :mean -0.15667464000000003}\n",
      " {:count 6, :mean -0.18306953333333337}\n",
      " {:count 7, :mean -0.20331617142857145}\n",
      " {:count 8, :mean -0.262446275}\n",
      " {:count 9, :mean -0.21517335555555556}\n",
      " {:count 10, :mean -0.27947242})\n"
     ]
    }
   ],
   "source": [
    "(defn make-running-stats-mapper []\n",
    "    (let [running-stats (atom {:count 0 :mean 0})]\n",
    "        (fn [z]\n",
    "            (let [{x :mean, n :count} @running-stats\n",
    "                  n+1 (inc n)\n",
    "                  K   (/ 1.0 n+1)]\n",
    "                (swap! running-stats conj\n",
    "                       [:count n+1]\n",
    "                       [:mean (+ x (* K (- z x)))]))\n",
    "            @running-stats)))\n",
    "(clojure.pprint/pprint (map (make-running-stats-mapper) zs))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### NUMERICAL CHECK"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 222,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.27947242"
      ]
     },
     "execution_count": 222,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(defn mean [zs] (/ (reduce + zs) (count zs)))\n",
    "(mean zs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## CORE.ASYNC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 223,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "(require \n",
    "    '[clojure.core.async \n",
    "      :refer \n",
    "      [sliding-buffer dropping-buffer buffer\n",
    "       <!! <! >! >!! go chan onto-chan close! \n",
    "       thread alts! alts!! timeout]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can write on the asynch thread and read on the UI thread:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 151,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "42\n"
     ]
    }
   ],
   "source": [
    "(let [c (chan)]\n",
    "    (go (>! c 42))\n",
    "    (println (<!! c))\n",
    "    (close! c))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can read on the async thread and write on the UI thread:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 152,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "42\n"
     ]
    }
   ],
   "source": [
    "(let [c (chan)]\n",
    "    (go (println (<! c)))\n",
    "    (>!! c 42)\n",
    "    (close! c))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In all cases, we must do any blocking call without timeout on the UI thread last. The following will hang:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 153,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#_(let [c (chan)]\n",
    "    (>!! c 42)\n",
    "    (go (println (<! c)))\n",
    "    (close! c))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We're won't block if we add a timeout, but we don't know how to write to a timeout channel:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 153,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "IllegalArgumentException No implementation of method: :put! of protocol: #'clojure.core.async.impl.protocols/WritePort found for class: clojure.lang.PersistentVector  clojure.core/-cache-protocol-fn (core_deftype.clj:568)\n"
     ]
    },
    {
     "ename": "class java.lang.IllegalArgumentException",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "       core_deftype.clj:   568 clojure.core$_cache_protocol_fn/invokeStatic                      ",
      "       core_deftype.clj:   560 clojure.core$_cache_protocol_fn/invoke                            ",
      "          protocols.clj:    18 clojure.core.async.impl.protocols$fn__10594$G__10589__10603/invoke",
      "              async.clj:   142 clojure.core.async$_GT__BANG__BANG_/invokeStatic                  ",
      "              async.clj:   137 clojure.core.async$_GT__BANG__BANG_/invoke                        ",
      "                     []:     2 user$eval6468/invokeStatic                                        ",
      "                     []:     1 user$eval6468/invoke                                              ",
      "               core.clj:  3105 clojure.core$eval/invokeStatic                                    ",
      "               core.clj:  3101 clojure.core$eval/invoke                                          ",
      "               main.clj:   240 clojure.main$repl$read_eval_print__7408$fn__7411/invoke           ",
      "               main.clj:   258 clojure.main$repl$fn__7417/invoke                                 ",
      "               main.clj:   174 clojure.main$repl/doInvoke                                        ",
      "            RestFn.java:  1523 clojure.lang.RestFn/invoke                                        ",
      "               AFn.java:   152 clojure.lang.AFn/applyToHelper                                    ",
      "               AFn.java:   144 clojure.lang.AFn/applyTo                                          ",
      "               core.clj:   646 clojure.core$apply/invokeStatic                                   ",
      "               core.clj:  1881 clojure.core$with_bindings_STAR_/invokeStatic                     ",
      "            RestFn.java:   425 clojure.lang.RestFn/invoke                                        ",
      "               AFn.java:    22 clojure.lang.AFn/run                                              ",
      "ThreadPoolExecutor.java:  1142 java.util.concurrent.ThreadPoolExecutor/runWorker                 ",
      "ThreadPoolExecutor.java:   617 java.util.concurrent.ThreadPoolExecutor$Worker/run                ",
      "            Thread.java:   748 java.lang.Thread/run                                              "
     ]
    }
   ],
   "source": [
    "(let [c (chan)]\n",
    "    (>!! (alts!! [c (timeout 500)]) 42)\n",
    "    (go (println (<! c)))\n",
    "    (close! c))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "We can't read from a `timeout` channel either, but at least we won't hang. Here, we block-read too early on the UI thread:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 206,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "IllegalArgumentException No implementation of method: :take! of protocol: #'clojure.core.async.impl.protocols/ReadPort found for class: clojure.lang.PersistentVector  clojure.core/-cache-protocol-fn (core_deftype.clj:568)\n"
     ]
    },
    {
     "ename": "class java.lang.IllegalArgumentException",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "       core_deftype.clj:   568 clojure.core$_cache_protocol_fn/invokeStatic                      ",
      "       core_deftype.clj:   560 clojure.core$_cache_protocol_fn/invoke                            ",
      "          protocols.clj:    15 clojure.core.async.impl.protocols$fn__10576$G__10571__10583/invoke",
      "              async.clj:   111 clojure.core.async$_LT__BANG__BANG_/invokeStatic                  ",
      "              async.clj:   106 clojure.core.async$_LT__BANG__BANG_/invoke                        ",
      "                     []:     2 user$eval7053/invokeStatic                                        ",
      "                     []:     1 user$eval7053/invoke                                              ",
      "               core.clj:  3105 clojure.core$eval/invokeStatic                                    ",
      "               core.clj:  3101 clojure.core$eval/invoke                                          ",
      "               main.clj:   240 clojure.main$repl$read_eval_print__7408$fn__7411/invoke           ",
      "               main.clj:   258 clojure.main$repl$fn__7417/invoke                                 ",
      "               main.clj:   174 clojure.main$repl/doInvoke                                        ",
      "            RestFn.java:  1523 clojure.lang.RestFn/invoke                                        ",
      "               AFn.java:   152 clojure.lang.AFn/applyToHelper                                    ",
      "               AFn.java:   144 clojure.lang.AFn/applyTo                                          ",
      "               core.clj:   646 clojure.core$apply/invokeStatic                                   ",
      "               core.clj:  1881 clojure.core$with_bindings_STAR_/invokeStatic                     ",
      "            RestFn.java:   425 clojure.lang.RestFn/invoke                                        ",
      "               AFn.java:    22 clojure.lang.AFn/run                                              ",
      "ThreadPoolExecutor.java:  1142 java.util.concurrent.ThreadPoolExecutor/runWorker                 ",
      "ThreadPoolExecutor.java:   617 java.util.concurrent.ThreadPoolExecutor$Worker/run                ",
      "            Thread.java:   748 java.lang.Thread/run                                              "
     ]
    }
   ],
   "source": [
    "(let [c (chan)]\n",
    "    (println (<!! (alts!! [c (timeout 500)])))\n",
    "    (go (>! c 42))\n",
    "    (close! c))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following illustrates writing to a blocking channel at random times and reads some data on the UI thread. It will leave values in the channel and thus leak the channel according to the documentation for `close!` here https://clojure.github.io/core.async/api-index.html#C. (we have"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 207,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-0.178654\n",
      "-0.676357\n",
      "0.0592247\n"
     ]
    }
   ],
   "source": [
    "(def echo-chan (chan))\n",
    "\n",
    "(doseq   [z zs] (go (Thread/sleep (rand 100)) (>! echo-chan z)))\n",
    "(dotimes [_ 3] (println (<!! echo-chan)))\n",
    "\n",
    "(close! echo-chan)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can chain channels. Reading from `echo-chan` may hang the UI thread because the UI thread races the internal thread that reads `echo-chan`, but the timeout trick works here as above. This will also leak channels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 207,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "IllegalArgumentException No implementation of method: :take! of protocol: #'clojure.core.async.impl.protocols/ReadPort found for class: clojure.lang.PersistentVector  clojure.core/-cache-protocol-fn (core_deftype.clj:568)\n"
     ]
    },
    {
     "ename": "class java.lang.IllegalArgumentException",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "       core_deftype.clj:   568 clojure.core$_cache_protocol_fn/invokeStatic                      ",
      "       core_deftype.clj:   560 clojure.core$_cache_protocol_fn/invoke                            ",
      "          protocols.clj:    15 clojure.core.async.impl.protocols$fn__10576$G__10571__10583/invoke",
      "              async.clj:   111 clojure.core.async$_LT__BANG__BANG_/invokeStatic                  ",
      "              async.clj:   106 clojure.core.async$_LT__BANG__BANG_/invoke                        ",
      "                     []:     6 user$eval7222/invokeStatic                                        ",
      "               core.clj:  3105 clojure.core$eval/invokeStatic                                    ",
      "               core.clj:  3101 clojure.core$eval/invoke                                          ",
      "               main.clj:   240 clojure.main$repl$read_eval_print__7408$fn__7411/invoke           ",
      "               main.clj:   258 clojure.main$repl$fn__7417/invoke                                 ",
      "               main.clj:   174 clojure.main$repl/doInvoke                                        ",
      "            RestFn.java:  1523 clojure.lang.RestFn/invoke                                        ",
      "               AFn.java:   152 clojure.lang.AFn/applyToHelper                                    ",
      "               AFn.java:   144 clojure.lang.AFn/applyTo                                          ",
      "               core.clj:   646 clojure.core$apply/invokeStatic                                   ",
      "               core.clj:  1881 clojure.core$with_bindings_STAR_/invokeStatic                     ",
      "            RestFn.java:   425 clojure.lang.RestFn/invoke                                        ",
      "               AFn.java:    22 clojure.lang.AFn/run                                              ",
      "ThreadPoolExecutor.java:  1142 java.util.concurrent.ThreadPoolExecutor/runWorker                 ",
      "ThreadPoolExecutor.java:   617 java.util.concurrent.ThreadPoolExecutor$Worker/run                ",
      "            Thread.java:   748 java.lang.Thread/run                                              "
     ]
    }
   ],
   "source": [
    "(def echo-chan (chan))\n",
    "(def repl-chan (chan))\n",
    "\n",
    "(dotimes [_ 10] (go (>! repl-chan (<! echo-chan))))\n",
    "(doseq   [z zs] (go (Thread/sleep (rand 100)) (>! echo-chan z)))\n",
    "(dotimes [_ 3] (println (<!! (alts!! [echo-chan (timeout 500)]))))\n",
    "\n",
    "(close! echo-chan)\n",
    "(close! repl-chan)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`println` on a `go` process works if we wait long enough. This, of course, is bad practice or \"code smell.\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 208,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-0.324796\n",
      "-1.48014\n",
      "-0.676357\n"
     ]
    }
   ],
   "source": [
    "(def echo-chan (chan))\n",
    "\n",
    "(doseq   [z zs] (go (Thread/sleep (rand 100)) (>! echo-chan z)))\n",
    "(dotimes [_ 3]  (go (println (<! echo-chan))))\n",
    "\n",
    "(Thread/sleep 500) ; no visible output if you remove this line.\n",
    "(close! echo-chan)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ASYNC RUNNING MEAN"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We want our `running-stats` function called with the data delivered at random times and in random order. A _transducer_, `(map mapper)`, lets us collect items off the buffer. The size of the buffer does not matter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 243,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{:count 1, :mean 0.0592247}\n",
      "{:count 2, :mean -0.30856615}\n",
      "{:count 3, :mean -0.20974706666666665}\n",
      "{:count 4, :mean -0.23607129999999998}\n",
      "{:count 5, :mean -0.023196039999999973}\n",
      "{:count 6, :mean 0.00783830000000002}\n",
      "{:count 7, :mean -0.20473002857142852}\n",
      "{:count 8, :mean -0.20147052499999996}\n",
      "{:count 9, :mean -0.27443646666666666}\n",
      "{:count 10, :mean -0.27947242}\n"
     ]
    }
   ],
   "source": [
    "(defn async-scan [zs mapper effector]\n",
    "    (let [transducer (map mapper)\n",
    "          echo-chan (chan (buffer 1) transducer)]\n",
    "        (doseq [z zs] (go (Thread/sleep (rand 100)) (>! echo-chan z)))\n",
    "        (dotimes [_ (count zs)] (effector (<!! echo-chan)))\n",
    "        (close! echo-chan)))\n",
    "(async-scan zs (make-running-stats-mapper) println)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now have complete parity between space (vector `zs`) and time (values on `echo-chan` in random order)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## RUNNING STDDEV"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### BRUTE-FORCE (SCALAR VERSION)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The definition of variance is the following, for $N>1$:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\frac{1}{N-1}\\sum\\limits_{i=1}^{N}\\left({z_i-\\bar{z}_N}\\right)^2\\tag{1}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### SSR: SUM OF SQUARED RESIDUALS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 210,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.5566483654807355"
      ]
     },
     "execution_count": 210,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(defn ssr [sequ]\n",
    "    (let [m (mean sequ)]\n",
    "        (reduce #(+ %1 (* (- %2 m) (- %2 m)))\n",
    "                0 sequ)))\n",
    "(ssr zs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### VARIANCE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 211,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.3951831517200817"
      ]
     },
     "execution_count": 211,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(defn variance [sequ]\n",
    "    (let [n (count sequ)]\n",
    "        (case n\n",
    "            0 0\n",
    "            1 (first sequ)\n",
    "            #_default (/ (ssr sequ) (- n 1.0)))))\n",
    "(variance zs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's do a smaller example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 212,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2017.0"
      ]
     },
     "execution_count": 212,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(def z2s [55. 89. 144.])\n",
    "(variance z2s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### REALLY DUMB VARIANCE FORMULA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is a really dumb recurrent form. It's dumb because\n",
    "1. it requires the whole sequence up front, so it cannot function in constant memory\n",
    "3. the intermediate values are meaningless because they refer to the final mean and count, not to the intermediate ones\n",
    "\n",
    "But, the final value is correct."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 214,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0 840.5 865.0 2017.0)"
      ]
     },
     "execution_count": 214,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(reductions \n",
    "    (let [m (mean z2s) ; uh-oh, we refer to _all_ the data ??\n",
    "          c (count z2s)]\n",
    "        (fn [var z] (+ var (let [r (- z m)] ; residual\n",
    "                               (/ (* r r) (- c 1.0))))))\n",
    "    0 z2s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is sufficiently dumb that we won't bother with a thread-safe, stateful, or asynchronous form."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### SCHOOL VARIANCE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\frac{1}{N-1}\\sum\\limits_{i=1}^{N}\\left({z_i-\\bar{z}_N}\\right)^2 =\n",
    "\\frac{1}{N-1}\\left(\\sum\\limits_{i=1}^{N}\\left(z_i^2\\right)-N\\,{\\bar{z}_N^2}\\right)\\tag{2}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Instead of the sum of squared residuals, $ssr$, accumulate the sum of squares, $ssq$, which grows quickly. The final result is exposed to _catastrophic cancellation_, but works for our small example. We get a clue that something is not optimal with this form by the fact that we don't use the old variance to compute the new variance. We do better below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 245,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "({:count 1, :mean 55.0, :variance 0.0, :ssq 3025.0}\n",
      " {:count 2, :mean 72.0, :variance 578.0, :ssq 10946.0}\n",
      " {:count 3, :mean 96.0, :variance 2017.0, :ssq 31682.0})\n"
     ]
    }
   ],
   "source": [
    "(defn make-school-stats-mapper []\n",
    "    (let [running-stats (atom {:count 0, :mean 0, \n",
    "                               :variance 0, :ssq 0})]\n",
    "        (fn [z]\n",
    "            (let [{x :mean, n :count, s :ssq} @running-stats\n",
    "                  n+1 (inc n)\n",
    "                  K   (/ 1.0 n+1)\n",
    "                  r   (- z x)\n",
    "                  x2  (+ x (* K r))\n",
    "                  s2  (+ s (* z z))]\n",
    "                (swap! running-stats conj\n",
    "                       [:count    n+1]\n",
    "                       [:mean     x2 ]\n",
    "                       [:ssq      s2]\n",
    "                       [:variance (/ (- s2 (* n+1 x2 x2)) (max 1 n))]))\n",
    "            @running-stats)))\n",
    "(clojure.pprint/pprint (map (make-school-stats-mapper) z2s))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### RECURRENT VARIANCE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We already know the recurrence for the mean:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$x\\leftarrow{x+K\\cdot(z-x)=x+\\frac{1}{n+1}(z-x)}\\tag{3}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The recurrence for the variance takes a little work to prove:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$v\\leftarrow\\frac{\\left(n-1\\right)v+K\\,n\\,\\left(z-x\\right)^2}{\\max(1,n)}\\tag{4}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 246,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "({:count 1, :mean 55.0, :variance 0.0, :ssq 0}\n",
      " {:count 2, :mean 72.0, :variance 578.0, :ssq 0}\n",
      " {:count 3, :mean 96.0, :variance 2017.0, :ssq 0})\n"
     ]
    }
   ],
   "source": [
    "(defn make-recurrent-stats-mapper []\n",
    "    (let [running-stats (atom {:count 0, :mean 0, \n",
    "                               :variance 0, :ssq 0})]\n",
    "        (fn [z]\n",
    "            (let [{x :mean, n :count, v :variance} @running-stats\n",
    "                  n+1 (inc n)\n",
    "                  K   (/ 1.0 (inc n))\n",
    "                  r   (- z x)\n",
    "                  x2  (+ x (* K r))\n",
    "                  ssr (+ (* (- n 1) v) ; old ssr is (* (- n 1) v)\n",
    "                         (* K n r r))]\n",
    "                (swap! running-stats conj\n",
    "                       [:count    n+1]\n",
    "                       [:mean     x2 ]\n",
    "                       [:variance (/ ssr  (max 1 n))]))\n",
    "            @running-stats)))\n",
    "(clojure.pprint/pprint (map (make-recurrent-stats-mapper) z2s))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Of course, this works asynchronously, with potentially different intermediate values because the order is random."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 247,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{:count 1, :mean 144.0, :variance 0.0, :ssq 0}\n",
      "{:count 2, :mean 116.5, :variance 1512.5, :ssq 0}\n",
      "{:count 3, :mean 96.0, :variance 2017.0, :ssq 0}\n"
     ]
    }
   ],
   "source": [
    "(async-scan z2s (make-recurrent-stats-mapper) println)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### WELFORD'S VARIANCE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above is equivalent, algebraically and numerically, to Welford's famous recurrence for the sum of squared residuals $S$. Remember that, in recurrences, we want everything on the right-hand sides of equations or left arrows to be be old, _prior_ statistics except for the new observation / measurement / input $z$. Welford's requires the new, _posterior_ mean on the right-hand side, so it's not as elegant as our recurrence above. However, it is easier to remember!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$S\\leftarrow{S} + \\left(z-x_N\\right)\\left(z-x_{N+1}\\right)=S+\\left(z-x\\right)\\left(z-\\left(x+K\\,\\left(z-x\\right)\\right)\\right)\\tag{5}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 248,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "({:count 1, :mean 55.0, :variance 0.0}\n",
      " {:count 2, :mean 72.0, :variance 578.0}\n",
      " {:count 3, :mean 96.0, :variance 2017.0})\n"
     ]
    }
   ],
   "source": [
    "(defn make-welfords-stats-mapper []\n",
    "    (let [running-stats (atom {:count 0, :mean 0, :variance 0})]\n",
    "        (fn [z]\n",
    "            (let [{x :mean, n :count, v :variance} @running-stats\n",
    "                  n+1 (inc n)\n",
    "                  K   (/ 1.0 n+1)\n",
    "                  r   (- z x)\n",
    "                  x2  (+ x (* K r))\n",
    "                  ssr (+ (* (- n 1) v) \n",
    "                         (* (- z x) (- z x2)))] ; <-- only difference to recurrent variance\n",
    "                (swap! running-stats conj\n",
    "                       [:count    n+1]\n",
    "                       [:mean     x2 ]\n",
    "                       [:variance (/ ssr  (max 1 n))]))\n",
    "            @running-stats)))\n",
    "(clojure.pprint/pprint (map (make-welfords-stats-mapper) z2s))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## WINDOWED STATISTICS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppose we want running statistics looking part way back. For example, suppose we have $N=10$ ten data and we want the statitics in a window of length $w=3$ three behind the current value, inclusively. When the first datum arrives, the window and the total include one datum. The window overhangs the left until the third datum. When the fourth datum arrives, the window contains three data and the total contain four data. After the tenth datum, we may consider three more steps marching the window \"off the cliff\" to the right. The following figure illustrates (the first row corresponds to $n=0$, not to $n=1$):"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"sliding-window.png\" style=\"width: 200px;\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We won't go through the derivation of the following formulas, but rather say that they have been vetted independently (in a C program and in a Mathematica program). The following table shows a unit test that we must reproduce. The notation is explained after the table."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"sliding-window-unit-test-001.png\" style=\"width: 800px;\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use the language of recurrences and avoid indices other than for datum $z_j$. The best algorithm we have found for tracking that datum (first datum to the left of the window) is to keep a FIFO queue of length $w$. This is still constant memory because it depends only on the length $w$ of the window, not on the length of the data stream."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "variable | description \n",
    "---------|-------------\n",
    "$n$      | prior count of data points; equals $0$ when considering the first point\n",
    "$z$      | current data point\n",
    "$w$      | fixed, constant, maximum width of window; $w\\geq{1}$\n",
    "$j$      | posterior number of points left of the window; $j\\geq{0}$\n",
    "$u$      | posterior number of points including $z$ in the running window; $1\\leq{u}\\leq{w}$\n",
    "$m$      | prior mean of all points, not including $z$\n",
    "$m'$     | posterior mean of all points including $z$\n",
    "$m_j$    | prior mean of points left of the window, lagging $w$ behind $m$\n",
    "$m'_j$   | posterior mean of points left of the window\n",
    "$m'_w$   | posterior mean of points in the window\n",
    "$v$      | prior variance\n",
    "$v'$     | posterior variance of points including $z$\n",
    "$v_j$    | prior variance of points left of the window, lagging $w$ behind $u_n$\n",
    "$v'_j$   | posterior variance of points left of the window\n",
    "$v'_w$   | posterior variance of points within the window"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\begin{align}\n",
    "j     &= \\max(0,n+1-w)               \\tag{6}\\\\\n",
    "u     &= n-j+1                       \\tag{7}\\\\\n",
    "m'    &= m+\\frac{z-m}{n+1}           \\tag{8}\\\\\n",
    "m'_j  &= \\begin{cases} \n",
    "  m_j+\\frac{z_j-m_j}{j} & j>0 \\\\\n",
    "  0 & \\mathrm{otherwise}\n",
    "\\end{cases}                            \\tag{9}\\\\\n",
    "m'_w  &= \\frac{(n+1)\\,m'-j\\,m'_j}{u} \\tag{10}\\\\\n",
    "v'    &= \\frac{(n-1)\\,v+\\frac{n}{n+1}\\left(z-m\\right)^2}{\\max(1,n)}               \\tag{11}\\\\\n",
    "v'_j  &= \\begin{cases}\n",
    "  \\frac{j-2}{j-1}\\,v_j+\\frac{1}{j}\\,\\left(z_j-m_j\\right)^2 & j>1 \\\\\n",
    "  0 & \\mathrm{otherwise}\n",
    "\\end{cases}                                                                         \\tag{12}\\\\\n",
    "v'_w  &= \\frac{n\\,v'+(n-w)\\,v'_j+(n+1)\\,{m'}^2-j\\,{m'_j}^2-u\\,{m'_w}^2}{\\max(1,u-1)}    \\tag{13}\n",
    "\\end{align}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 229,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "#'user/z3s"
      ]
     },
     "execution_count": 229,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(def z3s [0.857454, 0.312454, 0.705325, 0.839363, 1.63781, 0.699257, -0.340016, -0.213596, -0.0418609, 0.054705])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 274,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "#'user/push-to-back"
      ]
     },
     "execution_count": 274,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(defn push-to-back [item vek]\n",
    "    (conj (vec (drop 1 vek)) item))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 280,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "({:n 1.0, :m 0.857454, :v 0.0, :win [0 0 0.857454], :mj 0, :vj 0}\n",
      " {:n 2.0,\n",
      "  :m 0.584954,\n",
      "  :v 0.14851250000000002,\n",
      "  :win [0 0.857454 0.312454],\n",
      "  :mj 0,\n",
      "  :vj 0}\n",
      " {:n 3.0,\n",
      "  :m 0.6250776666666666,\n",
      "  :v 0.07908597588033334,\n",
      "  :win [0.857454 0.312454 0.705325],\n",
      "  :mj 0,\n",
      "  :vj 0}\n",
      " {:n 4.0,\n",
      "  :m 0.678649,\n",
      "  :v 0.06420353494066668,\n",
      "  :win [0.312454 0.705325 0.839363],\n",
      "  :mj 0.857454,\n",
      "  :vj 0}\n",
      " {:n 5.0,\n",
      "  :m 0.8704812,\n",
      "  :v 0.23215061598970005,\n",
      "  :win [0.705325 0.839363 1.63781],\n",
      "  :mj 0.584954,\n",
      "  :vj 0.14851250000000002}\n",
      " {:n 6.0,\n",
      "  :m 0.8419438333333333,\n",
      "  :v 0.1906067805693667,\n",
      "  :win [0.839363 1.63781 0.699257],\n",
      "  :mj 0.6250776666666666,\n",
      "  :vj 0.07908597588033334}\n",
      " {:n 7.0,\n",
      "  :m 0.6730924285714286,\n",
      "  :v 0.35841456203828576,\n",
      "  :win [1.63781 0.699257 -0.340016],\n",
      "  :mj 0.678649,\n",
      "  :vj 0.06420353494066666}\n",
      " {:n 8.0,\n",
      "  :m 0.562256375,\n",
      "  :v 0.4054895279174108,\n",
      "  :win [0.699257 -0.340016 -0.213596],\n",
      "  :mj 0.8704812,\n",
      "  :vj 0.2321506159897}\n",
      " {:n 9.0,\n",
      "  :m 0.49513223333333334,\n",
      "  :v 0.39535419047811504,\n",
      "  :win [-0.340016 -0.213596 -0.0418609],\n",
      "  :mj 0.8419438333333333,\n",
      "  :vj 0.19060678056936667}\n",
      " {:n 10.0,\n",
      "  :m 0.45108950999999997,\n",
      "  :v 0.3708235618778233,\n",
      "  :win [-0.213596 -0.0418609 0.054705],\n",
      "  :mj 0.6730924285714286,\n",
      "  :vj 0.35841456203828576})\n"
     ]
    }
   ],
   "source": [
    "(defn make-sliding-stats-mapper [w]\n",
    "    (let [running-stats (atom {:n 0, :m 0, :v 0,\n",
    "                               :win (vec (repeat w 0))\n",
    "                               :mj 0, :vj 0})]\n",
    "        (fn [z]\n",
    "            (let [{:keys [m n v win mj vj]} @running-stats\n",
    "                  zj   (first win)\n",
    "                  win' (push-to-back z win)\n",
    "                  n+1  (double (inc n))\n",
    "                  n-1  (double (dec n))\n",
    "                  K    (/ 1.0 n+1)\n",
    "                  Kv   (* n K)\n",
    "                  r    (- z m)\n",
    "                  j    (max 0, (- n+1 w))\n",
    "                  u    (- n+1 j)\n",
    "                  m'   (+ m (* K r))\n",
    "                  rj   (- zj mj)\n",
    "                  mj'  (if (> j 0), (+ mj (/ rj j)), 0)\n",
    "                  mw'  (/ (- (* n+1 m') (* j mj')) u)\n",
    "                  v'   (/  (+ (* n-1 v) (* Kv r r))\n",
    "                           (max 1 n))\n",
    "                  vj'  (if (> j 1)\n",
    "                           (let [j21 (/ (- j 2.0) \n",
    "                                        (- j 1.0))]\n",
    "                               (+ (* j21 vj) \n",
    "                                  (/ (* rj rj) j)))\n",
    "                           0)\n",
    "                  vw'  (let [t1 (- (* n v')\n",
    "                                   (* (- n w) vj'))\n",
    "                             t2 (- (* n+1 m' m')\n",
    "                                   (* j mj' mj'))\n",
    "                             t3 (- (* u mw' mw'))]\n",
    "                           (/  (+ t1 t2 t3)\n",
    "                               (max 1 (- u 1))))\n",
    "                  ]\n",
    "                (swap! running-stats conj\n",
    "                       [:n    n+1 ]\n",
    "                       [:m    m'  ]\n",
    "                       [:v    v'  ]\n",
    "                       [:mj   mj' ]\n",
    "                       [:vj   vj' ]\n",
    "                       [:win  win']))\n",
    "            @running-stats)))\n",
    "(clojure.pprint/pprint (map (make-sliding-stats-mapper 3) z3s))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This passes the unit test."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Clojure",
   "language": "clojure",
   "name": "clojure"
  },
  "language_info": {
   "file_extension": ".clj",
   "mimetype": "text/x-clojure",
   "name": "clojure",
   "version": "1.8.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}