diff options
Diffstat (limited to 'hawkes_experiments')
| -rw-r--r-- | hawkes_experiments/main.py | 75 |
1 files changed, 73 insertions, 2 deletions
diff --git a/hawkes_experiments/main.py b/hawkes_experiments/main.py index 9ed4205..0ec7c5d 100644 --- a/hawkes_experiments/main.py +++ b/hawkes_experiments/main.py @@ -4,6 +4,8 @@ from random import random, gauss import sys from itertools import product from multiprocessing import Pool +from scipy.optimize import minimize +import numpy as np GR = (sqrt(5) - 1) / 2 @@ -41,9 +43,34 @@ def approx(x): return x +def ll_gradient(lamb, alpha, mu): + l1 = [(1 + 0.43 * sin(0.0172 * t1 + 4.36), + sum(1. / d ** 2 * mu * exp(-mu * (t1 - t2)) for (n2, t2, d) in s)) + for ((n1, t1), s) in event_edges.iteritems()] + r1 = sum(log(lamb * a + alpha * b) for a, b in l1) + g1lambda = sum(a / (lamb * a + alpha * b) for a, b in l1) + g1alpha = sum(b / (lamb * a + alpha * b) for a, b in l1) + + t2 = sum(sum(1. / d ** 2 * approx(mu * (nodes[n2][0] - t1)) + for n2, d in edges[n1].iteritems() + if nodes[n2][0] > t1) + for (n1, t1) in iter_events(events)) + r2 = alpha * t2 + g2lambda = 0 + g2alpha = t2 + + t3 = sum(node[1] for node in nodes.itervalues()) + r3 = lamb * t3 + g3lambda = t3 + g3alpha = 0 + + return [-(r1 - r2 - r3), np.array([-(g1lambda - g2lambda - g3lambda), + -(g1alpha - g2alpha - g3alpha)])] + + def ll(lamb, alpha, mu): r1 = sum(log(lamb * (1 + 0.43 * sin(0.0172 * t1 + 4.36)) - + sum(alpha / d ** 2* mu * exp(-mu * (t1 - t2)) + + sum(alpha / d ** 2 * mu * exp(-mu * (t1 - t2)) for (n2, t2, d) in s)) for ((n1, t1), s) in event_edges.iteritems()) r2 = sum(sum(alpha / d ** 2 * approx(mu * (nodes[n2][0] - t1)) @@ -51,6 +78,7 @@ def ll(lamb, alpha, mu): if nodes[n2][0] > t1) for (n1, t1) in iter_events(events)) r3 = lamb * sum(node[1] for node in nodes.itervalues()) + return -(r1 - r2 - r3) @@ -81,6 +109,31 @@ def sa(x, y, z, sigma=0.5, niter=70, fc=None): return y, z, fo +def sa_bis(x, y, z, sigma=0.5, niter=30, fc=None): + T = 0.1 + e = 1.1 + if fc: + fo = fc + else: + fo = ll(x, y, z) + for _ in xrange(niter): + sys.stderr.write("sa: " + " ".join(map(str, [T, sigma, x, y, z, fo])) + + "\n") + sys.stderr.flush() + zn = max(z + gauss(0, sigma * z), 0) + fn = ll(x, y, zn) + if fn < fo or exp((fo - fn) / T) > random(): + z = zn + sigma *= 1.5 + fo = fn + else: + sigma /= e + if sigma < 1e-5: + break + T *= 0.99 + return z, fo + + def optimize_with_sa(x, y, z, niter=4): def f(x): @@ -95,6 +148,24 @@ def optimize_with_sa(x, y, z, niter=4): sys.stdout.flush() +def optimize_with_bfgs(x, y, z, niter=4): + + def f(x0, zum): + return ll_gradient(x0[0], x0[1], zum) + + z, fc = sa_bis(x, y, z) + for _ in xrange(niter): + result = minimize(f, np.array([x, y]), method='L-BFGS-B', jac=True, + bounds=[(0, None), (0, None)], + options={"maxiter": 50, "disp": True}, + args=[z]) + x, y = result.x + fc = result.fun + z, fc = sa_bis(x, y, z, fc=fc) + print x, y, z, fc + sys.stdout.flush() + + def optimize_with_gss(x, y, z, niter=5): def f(x): @@ -143,7 +214,7 @@ if __name__ == "__main__": with open(sys.argv[1]) as fh: l = [map(float, line.strip().split()[:3]) for line in fh] for e in l: - optimize_with_gss(*e) + optimize_with_bfgs(*e) # print ll(x, y, z) # coarse_search() |