From 8cadf3d79c0a6a5ab8f7cf32c73b5471a9bc5f38 Mon Sep 17 00:00:00 2001
From: Vitalik Buterin <v@buterin.com>
Date: Sun, 9 Nov 2014 14:55:24 -0500
Subject: [PATCH] More estimators, and csvgen file

---
 stability/csvgen.py | 31 +++++++++++++++++++++++++++++++
 stability/fit.py    | 39 ++++++++++++++++++++++-----------------
 2 files changed, 53 insertions(+), 17 deletions(-)
 create mode 100644 stability/csvgen.py

diff --git a/stability/csvgen.py b/stability/csvgen.py
new file mode 100644
index 0000000..e0427d5
--- /dev/null
+++ b/stability/csvgen.py
@@ -0,0 +1,31 @@
+import fit
+import spread
+
+tests = [
+    [fit.simple_estimator, [1], [1]],
+    [fit.simple_estimator, [1], [1.2]],
+    [fit.minimax_estimator, [1], [1.2]],
+    [fit.diff_estimator, [1, 0, 0.001], [1.2, 10, 1]],
+    [fit.diff_estimator, [1, 0, 0.001, 0], [1.2, 10, 1, 5]],
+    [fit.ndiff_estimator, [1, 0, 0, 0.001], [1.2, 10, 1, 1]],
+]
+
+vals = [fit.optimize(t, mi, ma, rate=0.4, rounds=12000, tries=10)
+        for t, mi, ma in tests]
+
+estimators = [t[0](*v) for t, v in zip(tests, vals)]
+
+scores = [fit.evaluate_estimates(e) for e in estimators]
+
+for t, v, e, s in zip(tests, vals, estimators, scores):
+    print v, s
+
+adjestimators = [[e/est[0] for e in est] for est in estimators]
+adjprices = [[p/e*est[0] for e, p in 
+             zip(est, fit.prices)] for est in estimators]
+output = [['Price'] + fit.prices]
+for i in range(len(tests)):
+    output.append(['Estimator %d' % (i+1)] + adjestimators[i])
+    output.append(['Adjusted price %d' % (i+1)] + adjprices[i])
+
+spread.save('o.csv', zip(*output))
diff --git a/stability/fit.py b/stability/fit.py
index 7c17ede..4045fde 100644
--- a/stability/fit.py
+++ b/stability/fit.py
@@ -27,7 +27,7 @@ def minimax_estimator(fac):
     return o
 
 
-def diff_estimator(fac, dw, mf):
+def diff_estimator(fac, dw, mf, exp=1):
     o = [1]
     derivs = [0] * 14
     for i in range(14, len(diffs)):
@@ -37,11 +37,11 @@ def diff_estimator(fac, dw, mf):
     vals = [max(diffs[i] + derivs[i] * dw, diffs[i] * mf) for i in range(len(diffs))]
     for i in range(1, len(diffs)):
         if vals[i] * 1.0 / vals[i-1] > fac:
-            o.append(o[-1] * vals[i] * 1.0 / vals[i-1] / fac)
+            o.append(o[-1] * 1.0 / fac * (vals[i] / vals[i-1])**exp)
         elif vals[i] > vals[i-1]:
             o.append(o[-1])
         else:
-            o.append(o[-1] * vals[i] * 1.0 / vals[i-1])
+            o.append(o[-1] * 1.0 * (vals[i] / vals[i-1])**exp)
     return o
 
 
@@ -101,21 +101,26 @@ def evaluate_estimates(estimates, crossvalidate=False):
 
 
 # Simulated annealing optimizer
-def optimize(producer, floors, ceilings, rate=0.7):
-    vals = [f*0.5+c*0.5 for f, c in zip(floors, ceilings)]
-    y = evaluate_estimates(producer(*vals))
-    for i in range(1, 5000):
-        stepsizes = [(f*0.5-c*0.5) / i**rate for f, c in zip(floors, ceilings)]
-        steps = [(random.random() * 2 - 1) * s for s in stepsizes]
-        newvals = [max(mi, min(ma, v+s)) for v, s, mi, ma in zip(vals, steps, floors, ceilings)]
-        newy = evaluate_estimates(producer(*newvals))
-        if newy < y:
-            vals = newvals
-            y = newy
-        if not i % 1000:
-            print i, vals, y
+def optimize(producer, floors, ceilings, rate=0.7, rounds=5000, tries=1):
+    bestvals, besty = None, 999999999999999
+    for t in range(tries):
+        print 'Starting test %d of %d' % (t + 1, tries)
+        vals = [f*0.5+c*0.5 for f, c in zip(floors, ceilings)]
+        y = evaluate_estimates(producer(*vals))
+        for i in range(1, rounds):
+            stepsizes = [(f*0.5-c*0.5) / i**rate for f, c in zip(floors, ceilings)]
+            steps = [(random.random() * 2 - 1) * s for s in stepsizes]
+            newvals = [max(mi, min(ma, v+s)) for v, s, mi, ma in zip(vals, steps, floors, ceilings)]
+            newy = evaluate_estimates(producer(*newvals))
+            if newy < y:
+                vals = newvals
+                y = newy
+            if not i % 1000:
+                print i, vals, y
+        if y < besty:
+            bestvals, besty = vals, y
         
-    return vals
+    return bestvals
 
 
 def score(producer, *vals):