Skip to content

Instantly share code, notes, and snippets.

@iantruslove
Last active October 13, 2015 14:38
Show Gist options
  • Save iantruslove/193eb4258c5fa993d887 to your computer and use it in GitHub Desktop.
Save iantruslove/193eb4258c5fa993d887 to your computer and use it in GitHub Desktop.
Numerically stable math
(ns jupiter.test.utils.math
(:require [clojure.test :refer :all]
[jupiter.utils.math :refer :all]))
(defn float=
([x y] (float= x y 0.000001))
([x y epsilon]
(let [scale (if (or (zero? x) (zero? y)) 1 (Math/abs x))]
(<= (Math/abs (- x y)) (* scale epsilon)))))
(deftest test-k+
(testing "basics"
(is (= [12 0] (k+ 6 6)))
(is (= [12.1 0.0] (k+ [6 0.1] 6)))
(is (= [12.1 0.0] (k+ 6 [6 0.1])))
(is (= [12 0] (k+ [6 0] 6))))
(testing "numerical stability at double precision"
(let [iterations 300000 ; crank this up and look at the
; regular-addition-error-term grow...
regular-addition-error-term
(+ -1 (reduce + (repeat iterations (double (/ 1 iterations)))))
kahan-result (reduce k+ (repeat iterations (double (/ 1 iterations))))
kahan-addition-error-term (+ -1
(first kahan-result)
(second kahan-result))]
(is (> regular-addition-error-term 1E-12))
(is (< kahan-addition-error-term 1E-15)))))
(defn square [x]
(* x x))
(deftest test-incremental-arithmetic
(let [iterations 1000
denominator 10
integer-series (repeatedly iterations #(rand-int denominator))
number-series (map #(double (/ % denominator)) integer-series)
correctable-series (map #(vector % 1) number-series)
accuracy 0.0000000000000001]
(testing "mean"
(let [calculated-mean (double (/ (reduce + integer-series) iterations denominator))
naive-incremental-mean (double
(first
(reduce naive-incremental-mean
correctable-series)))
stable-incremental-mean (double
(correct
(first
(reduce stable-incremental-mean
correctable-series))))]
(is (float= calculated-mean naive-incremental-mean))
(is (not (float= calculated-mean naive-incremental-mean accuracy)))
(is (float= calculated-mean stable-incremental-mean accuracy)))
;; TODO: test that incremental averaging on partitions of the
;; series yields good results
)
(testing "standard deviation"
(let [data-series [2 4 4 4 5 5 7 9]
calculated-mean (/ (reduce + data-series) (count data-series))
differences (map (partial + (- calculated-mean)) data-series)
square-differences (map square differences)
variance (/ (reduce + square-differences) (count square-differences))
calculated-std-dev (Math/sqrt variance)]
))))
(ns jupiter.utils.math)
(defprotocol StableArithmetic
(stable+ [this] [this y]))
(defrecord Correctable [value correction]
StableArithmetic
(stable+ [this] this)
(stable+ [this y] (+ this y)))
(extend-type java.lang.Number
StableArithmetic
(stable+ [this] this)
(stable+ [this y] (clojure.core/+ this y)))
(defn correctable [n]
(cond
(number? n) [n 0]
(and (sequential? n)
(number? (first n))
(number? (second n)))
n
:else (throw (Exception. "Cannot convert that to a correctable number."))))
(defn correct [n]
(let [n (correctable n)]
(+ (first n) (second n))))
(defn k+
"Kahan Increment: numerically stable summation of two numbers, by
keeping track of correction terms. Accepts either numbers or
2-tuples of [value correction]. Returns a [value correction] tuple.
See http://researcher.watson.ibm.com/researcher/files/us-ytian/stability.pdf"
[x y]
(let [[s1 c1] (correctable x)
[s2 c2] (correctable y)
corrected-s2 (+ s2 c1 c2)
sum (stable+ s1 corrected-s2)
correction (- corrected-s2 (- sum s1))]
[sum correction]))
(defn naive-incremental-mean
[[mu-a n-a] [mu-b n-b]]
(let [n (+ n-a n-b)
sigma (- mu-b mu-a)]
[(+ mu-a (/ (* n-b sigma) n)) n]))
(defn stable-incremental-mean [[mu-a n-a] [mu-b n-b]]
(let [mu-a (correctable mu-a)
mu-b (correctable mu-b)
n (stable+ n-a n-b)
sigma (- (correct mu-b) (correct mu-a))]
[(k+ mu-a (/ (* n-b sigma) n)) n]))
(defn naive-incremental-std-dev
"Takes tuples of second central moment, mean, n of the subsets of
data to combine."
[[m2-a mu-a n-a] [m2-b mu-b n-b]]
(let [n (+ n-a n-b)
sigma (- mu-b mu-a)
mu (+ mu-a (/ (* n-b sigma) n))
M2-a (* m2-a n-a)
M2-b (* m2-b n-b)])
)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment