1 files changed, 33 insertions, 13 deletions
diff --git a/simulation/main.py b/simulation/main.py
index 2425cf6..bdc5c97 100644
--- a/simulation/main.py
+++ b/simulation/main.py
@@ -1,8 +1,10 @@
 import numpy as np
+from numpy.linalg import norm
 import numpy.random as nr
 from scipy.optimize import minimize
 import matplotlib.pyplot as plt
 import seaborn
+from random import random, randint
 
 seaborn.set_style("white")
 
@@ -39,7 +41,7 @@ def infer(x, y):
         l, g = likelihood_gradient(p, x, y)
         return -l, -g
     x0 = np.ones(x.shape[1])
-    bounds = [(0, None)] * x.shape[1]
+    bounds = [(1e-10, None)] * x.shape[1]
     return minimize(f, x0, jac=True, bounds=bounds, method="L-BFGS-B").x
 
 
@@ -106,24 +108,42 @@ def simulate_cascade(x, graph):
         yield x, susc
 
 
-def simulate_cascades(n, graph):
+def uniform_source(graph):
+    x0 = np.zeros(graph.shape[0], dtype=bool)
+    x0[nr.randint(0, graph.shape[0])] = True
+    return x0
+
+
+def simulate_cascades(n, graph, source=uniform_source):
     for _ in xrange(n):
-        x0 = np.zeros(graph.shape[0], dtype=bool)
-        x0[nr.randint(0, graph.shape[0])] = True
+        x0 = source(graph)
         yield simulate_cascade(x0, graph)
 
 
 if __name__ == "__main__":
-    g = np.array([[0, 1, 1, 0], [1, 0, 0, 1], [1, 0, 0, 1], [0, 1, 1, 0]])
+    # g = np.array([[0, 1, 1, 0], [1, 0, 0, 1], [1, 0, 0, 1], [0, 1, 1, 0]])
+    g = np.array([[0, 0, 1], [0, 0, 0.5], [0, 0, 0]])
     p = 0.5
     g = np.log(1. / (1 - p * g))
     sizes = [100, 500, 1000, 5000, 10000]
-    error = []
+    # error = []
+
+    def source(graph):
+        x0 = np.zeros(graph.shape[0], dtype=bool)
+        a = randint(0, 1)
+        x0[a] = True
+        if random() > 0.01:
+            x0[1-a] = True
+        return x0
+
     for i in sizes:
-        cascades = simulate_cascades(i, g)
-        x, y = build_matrix(cascades, 0)
-        conf = bootstrap(x, y, n_iter=100)
-        estimand = np.linalg.norm(np.delete(conf - g[0], 0, axis=1), axis=1)
-        error.append(confidence_interval(*np.histogram(estimand, bins=50)))
-    plt.semilogx(sizes, error)
-    plt.show()
+        cascades = simulate_cascades(i, g, source=source)
+        x, y = build_matrix(cascades, 2)
+        e = infer(x, y)
+        print e
+
+        # conf = bootstrap(x, y, n_iter=100)
+        # estimand = np.linalg.norm(np.delete(conf - g[0], 0, axis=1), axis=1)
+        # error.append(confidence_interval(*np.histogram(estimand, bins=50)))
+    # plt.semilogx(sizes, error)
+    # plt.show()