diff options
| -rw-r--r-- | simulation/main.py | 46 |
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() |