diff options
| -rw-r--r-- | simulation/main.py | 28 | ||||
| -rw-r--r-- | simulation/vi.py | 26 | ||||
| -rw-r--r-- | simulation/vi_theano.py | 66 |
3 files changed, 92 insertions, 28 deletions
diff --git a/simulation/main.py b/simulation/main.py index 4fa8f6c..c2446d7 100644 --- a/simulation/main.py +++ b/simulation/main.py @@ -56,25 +56,17 @@ def build_cascade_list(cascades, collapse=False): return np.vstack(x), np.vstack(s) -def cascadeLkl(graph, infect, sus): - # There is a problem with the current implementation - # Note that you need to take into account the time diff between the label - # and the values being conditioned. Note also that the matrix if stacked as - # such will require to keep track of the state 0 of each cascade. - a = np.dot(infect, graph) - return np.log(1. - np.exp(-a[(infect[1:])*sus[1:]])).sum() \ - - a[(~infect[1:])*sus].sum() - - if __name__ == "__main__": g = np.array([[0, 0, 1], [0, 0, 0.5], [0, 0, 0]]) p = 0.5 g = np.log(1. / (1 - p * g)) - cascades = simulate_cascades(100, g) - cascade, y_obs = mn.build_matrix(cascades, 0) - conf = mn.bootstrap(x, y, n_iter=100) - - estimand = np.linalg.norm(np.delete(conf - g[0], 0, axis=1), axis=1) - error.append(mn.confidence_interval(*np.histogram(estimand, bins=50))) - plt.semilogx(sizes, error) - plt.show() + error = [] + sizes = [10, 10**2, 10**3] + for s in sizes: + cascades = simulate_cascades(s, g) + cascade, y_obs = mn.build_matrix(cascades, 0) + conf = mn.bootstrap(cascade, y_obs, n_iter=100) + estimand = np.linalg.norm(np.delete(conf - g[0], 0, axis=1), axis=1) + error.append(mn.confidence_interval(*np.histogram(estimand, bins=50))) + plt.semilogx(sizes, error) + plt.show() diff --git a/simulation/vi.py b/simulation/vi.py index 9604c7d..aeccb69 100644 --- a/simulation/vi.py +++ b/simulation/vi.py @@ -2,7 +2,10 @@ import time import main as mn import autograd.numpy as np from autograd import grad +import logging +logging.basicConfig(format='%(asctime)s - %(levelname)s - %(message)s', + level=logging.INFO) def g(m): assert (m > 0).all() @@ -47,29 +50,32 @@ def kl(params1, params0): grad_kl = grad(kl) -def sgd(mu1, sig1, mu0, sig0, cascades, n_e=100, lr=lambda t: 1e-2): +def sgd(mu1, sig1, mu0, sig0, cascades, n_e=100, lr=lambda t: 1e-2, n_print=10): g_mu1, g_sig1 = grad_kl((mu1, sig1), (mu0, sig0)) for t in xrange(n_e): lrt = lr(t) # learning rate mu1, sig1 = mu1 + lrt * g_mu1, sig1 + lrt * g_sig1 - for x, s in zip(*cascades): + for step, (x, s) in enumerate(zip(*cascades)): g_mu1, g_sig1 = grad_ll_full((mu1, sig1), x, s) mu1 = np.maximum(mu1 + lrt * g_mu1, 0) sig1 = np.maximum(sig1 + lrt * g_sig1, 1e-3) - res = np.sum(ll_full((mu1, sig1), x, s) for x, s in zip(*cascades)) + \ - kl((mu1, sig1), (mu0, sig0)) - print("Epoch: {}\t LB: {}\t Time: {}".format(t, res, time.time())) - print mu1 - print sig1 + res = np.sum(ll_full((mu1, sig1), x, s) for x, s in zip(*cascades))\ + + kl((mu1, sig1), (mu0, sig0)) + if step % n_print == 0: + logging.info("Epoch:{}\tStep:{}\tLB:{}\t".format(t, step, res)) + print mu1[0:2, 0:2] + print sig1[0:2, 0:2] if __name__ == '__main__': - graph = np.array([[0, 0, 1], [0, 0, 0.5], [0, 0, 0]]) + #graph = np.array([[0, 0, 1], [0, 0, 0.5], [0, 0, 0]]) + graph = np.random.binomial(2, p=.2, size=(10, 10)) p = 0.5 graph = np.log(1. / (1 - p * graph)) - cascades = mn.build_cascade_list(mn.simulate_cascades(1000, graph)) + print(graph[0:2, 0:2]) + cascades = mn.build_cascade_list(mn.simulate_cascades(500, graph)) mu0, sig0 = (1. + .2 * np.random.normal(size=graph.shape), 1 + .2 * np.random.normal(size=graph.shape)) mu1, sig1 = (1. + .2 * np.random.normal(size=graph.shape), 1 + .2 * np.random.normal(size=graph.shape)) - sgd(mu1, sig1, mu0, sig0, cascades, n_e=30) + sgd(mu1, sig1, mu0, sig0, cascades, n_e=30, n_print=1) diff --git a/simulation/vi_theano.py b/simulation/vi_theano.py new file mode 100644 index 0000000..562fa67 --- /dev/null +++ b/simulation/vi_theano.py @@ -0,0 +1,66 @@ +import main as mn +import theano +from theano import tensor as tsr +import theano.tensor.shared_randomstreams +import numpy as np + +n_cascades = 1000 +n_nodes = 4 +n_samples = 100 +srng = tsr.shared_randomstreams.RandomStreams(seed=123) +lr = 1e-2 +n_epochs = 10 + + +# Declare Theano variables +mu = theano.shared(.5 * np.random.rand(1, n_nodes, n_nodes), name="mu", + broadcastable=(True, False, False)) +sig = theano.shared(.3 * np.random.rand(1, n_nodes, n_nodes), name="sig", + broadcastable=(True, False, False)) +mu0 = theano.shared(.5 * np.random.rand(1, n_nodes, n_nodes), name="mu", + broadcastable=(True, False, False)) +sig0 = theano.shared(.3 * np.random.rand(1, n_nodes, n_nodes), name="sig", + broadcastable=(True, False, False)) +x = tsr.matrix(name='x', dtype='int8') +s = tsr.matrix(name='s', dtype='int8') + +# Construct Theano graph +theta = srng.normal((n_samples, n_nodes, n_nodes)) * sig + mu +y = tsr.clip(tsr.dot(x, theta), 1e-3, 1) +infect = tsr.log(1. - tsr.exp(-y[0:-1])).dimshuffle(1, 0, 2) +lkl_pos = tsr.sum(infect * (x[1:] & s[1:])) / n_samples +lkl_neg = tsr.sum(-y[0:-1].dimshuffle(1, 0, 2) * (~x[1:] & s[1:])) / n_samples + +lkl = lkl_pos + lkl_neg +kl = tsr.sum(tsr.log(sig / sig0) + (sig0**2 + (mu0 - mu)**2)/(2*sig)**2) +res = lkl + kl + +gmu, gsig = theano.gradient.grad(lkl, [mu, sig]) +gmukl, gsigkl = theano.grad(kl, [mu, sig]) + +# Compile into functions +loglkl_full = theano.function([x, s], lkl) +train = theano.function(inputs=[x, s], outputs=res, + updates=((mu, tsr.clip(mu + lr * gmu, 0, 1)), + (sig, tsr.clip(sig + lr * gsig, 1e-3, 1)))) +train_kl = theano.function(inputs=[], outputs=[], + updates=((mu, tsr.clip(mu + lr * gmukl, 0, 1)), + (sig, tsr.clip(sig + lr * gsigkl, 1e-3, 1)))) + + +if __name__ == "__main__": + graph = np.random.binomial(2, p=.2, size=(n_nodes, n_nodes)) + for k in range(len(graph)): + graph[k, k] = 0 + p = 0.5 + graph = np.log(1. / (1 - p * graph)) + cascades = mn.build_cascade_list(mn.simulate_cascades(n_cascades, graph), + collapse=True) + x_obs, s_obs = cascades[0], cascades[1] + for i in range(n_epochs): + train_kl() + for k in xrange(len(x_obs)/100): + cost = train(x_obs[k*100:(k+1)*100], s_obs[k*100:(k+1)*100]) + print(cost) + print(mu.get_value()) + print(graph) |