diff options
Diffstat (limited to 'experiments/ml.pyx')
| -rw-r--r-- | experiments/ml.pyx | 63 |
1 files changed, 30 insertions, 33 deletions
diff --git a/experiments/ml.pyx b/experiments/ml.pyx index c8c99cc..0d72c44 100644 --- a/experiments/ml.pyx +++ b/experiments/ml.pyx @@ -10,14 +10,9 @@ cdef DTYPE_t weight_success(int dist, int dt, DTYPE_t alpha, DTYPE_t delta, DTYPE_t w1, DTYPE_t w2, DTYPE_t w3): """weight for successful infection, exponential time model""" cdef DTYPE_t structural, temporal, result - structural = dist * log(delta) - # structural = plogis(w1,delta) * plogis(w2,delta) * plogis(w3,delta) - temporal = log(exp(alpha)-1.) - alpha*dt - # temporal = 1 - exp(-alpha*dt) - # if exp(-alpha*dt)==0.: print 'UNDERFLOW ERROR' - # temporal = 1. / (1. + (dt - 1.)/alpha)**0.01 - 1. / (1. + dt/alpha)**0.01 - result = structural + temporal - # print 'st', structural, temporal + structural = delta ** dist + temporal = log(exp(alpha)-1.) - alpha*dt/1. + result = log(structural) + temporal return result cdef DTYPE_t weight_failure(int dist, int dt, DTYPE_t alpha, DTYPE_t delta, @@ -25,24 +20,23 @@ cdef DTYPE_t weight_failure(int dist, int dt, DTYPE_t alpha, DTYPE_t delta, """weight for failed infection, exponential time model""" cdef DTYPE_t structural, temporal, result structural = delta ** dist - # structural = plogis(w1,delta) * plogis(w2,delta) * plogis(w3,delta) - temporal = exp(-alpha * dt) - # temporal = 1. - 1. / (1. + dt/alpha)**0.01 + temporal = exp(-alpha * dt/1.) result = log(1. - structural + structural * temporal) - # print 'stnv', structural, temporal return result def ml(dict root_victims, dict victims, dict non_victims, DTYPE_t age, DTYPE_t alpha, DTYPE_t delta): cdef: - int n_roots, n_victims, n_nodes, roots, i, dist, dt, t, l + int n_roots, n_victims, roots, i, dist, dt, t, l DTYPE_t beta, ll, beta_add, max_beta, max_beta_add list parents, failures, successes n_roots, n_victims = len(root_victims), len(victims) - n_nodes = 148152 + n_nodes = 5000 cdef: np.ndarray[DTYPE_t] probs = np.zeros(n_victims, dtype=DTYPE) np.ndarray[DTYPE_t] probs_fail = np.zeros(n_victims, dtype=DTYPE) + np.ndarray[DTYPE_t] parent_dists = np.zeros(n_victims, dtype=DTYPE) + np.ndarray[DTYPE_t] parent_dts = np.zeros(n_victims, dtype=DTYPE) np.ndarray[DTYPE_t] probs_nv = np.zeros(len(non_victims), dtype=DTYPE) # loop through victims @@ -51,45 +45,48 @@ def ml(dict root_victims, dict victims, dict non_victims, DTYPE_t age, # fail to infect it, also computes the probability that its most # likely parent infects it failures = [weight_failure(dist, dt, alpha, delta, w1, w2, w3) - for (dist, dt, w1, w2, w3) in parents] + for (prnt, dist, dt, w1, w2, w3) in parents] probs_fail[i] = sum(failures) successes = [weight_success(dist, dt, alpha, delta, w1, w2, w3) - for (dist, dt, w1, w2, w3) in parents] + for (prnt, dist, dt, w1, w2, w3) in parents] + dists = [dist for (prnt, dist, dt, w1, w2, w3) in parents] + dts = [dt for (prnt, dist, dt, w1, w2, w3) in parents] # find parent that maximizes log(p) - log(\tilde{p}) - probs[i] = max(s - failures[l] for l, s in enumerate(successes)) + # probs[i] = max(s - failures[l] for l, s in enumerate(successes)) + probs[i] = float("-inf") + for l, s in enumerate(successes): + prob = s - failures[l] + if prob > probs[i]: + probs[i] = prob + parent_dists[i] = dists[l] + parent_dts[i] = dts[l] + # probs_fail[i] = failures[l] # loop through non-victims for i, parents in enumerate(non_victims.itervalues()): # for each non victim node, compute the probability that all its # parents fail to infect it failures = [weight_failure(dist, dt, alpha, delta, w1, w2, w3) - for (dist, dt, w1, w2, w3) in parents] + for (prnt, dist, dt, w1, w2, w3) in parents] probs_nv[i] = sum(failures) - # print successes - # print failures - # print probs - # calculate log likelihood - # probs.sort(); probs = probs[::-1] # sort probs in descending order - # cdef: - # np.ndarray[DTYPE_t] cums = probs.cumsum() - ll = probs_fail.sum() # add probability that all edges to victims fail + ll = probs_fail.sum() # add probability that all edges to all victims fail ll += probs_nv.sum() # add probability that all edges to non_victims fail - # print 'probs', probs max_beta_add = float('-inf') # iterate over all victim nodes to find the optimal threshold - for beta in np.arange(0.001, .2, .002): - thresh = log(beta/(3012.*(1.-beta))) - # print 'beta:', beta, 'thresh:', thresh, 'infected:', len(probs[probs>=thresh]) - roots = n_roots + len(probs[probs<thresh]) + for beta in np.arange(0.01, 1., .01): + thresh = log(beta/(1000.*(1.-beta))) + seeds = probs<thresh + non_seeds = probs>=thresh + roots = n_roots + sum(seeds) beta_add = 0. # add probability for realized edges and subtract probability these edges fail - beta_add += (probs[probs>=thresh]).sum() + beta_add += (probs[non_seeds]).sum() # add probability for the seeds and non-seeds - beta_add += roots * log(beta/3012.) + (n_nodes-roots) * log(1. - beta) + beta_add += roots * log(beta/1000.) + (n_nodes-roots) * log(1. - beta) if beta_add > max_beta_add: max_beta = beta |