inhomPPLik.py 3.57 KB
 Chris Jewell committed Jan 16, 2019 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 ``````"""Name: inhomPPLik.py Author: Chris Jewell Created: 16th January 2019 Copyright: Chris Jewell 2019 License: MIT Warranty: None! Use at your own risk ;-) """ import numpy as np class Kernel: def __init__(self, c_matrix): """Represents a transmission kernel matrix #' #' Encapsulated data using a 'closure' technique #' and calling the returned function returns the full #' transmission matrix #' #' @example #' > kernel = Kernel(C) #' > kernel(0.5) #' [,1] [,2] [,3] [,4] [,5] #'[1,] 0.0 0.5 0.0 0.5 0.5 #'[2,] 0.5 0.0 0.5 0.0 0.5 #'[3,] 0.0 0.5 0.0 0.5 0.0 #'[4,] 0.5 0.0 0.5 0.0 0.0 #'[5,] 0.5 0.5 0.0 0.0 0.0 """ self.c_matrix = c_matrix def __call__(self,beta): return beta * self.c_matrix def prod_part(t_inf_j, events, B): """Returns product part of infection log likelihood for SIR epidemic :param t_inf_j: infection times of everyone beginning in S :param events: tuple of arrival and departure times of individuals in I state :param B: the full transmission matrix :return: scalar product part of the epidemic likelihoood """ `````` Chris Jewell committed Jan 16, 2019 44 `````` is_infected = t_inf_j < np.inf `````` Chris Jewell committed Jan 16, 2019 45 `````` waifw = (events[0][:, None] < t_inf_j[None, :]) & (t_inf_j[None, :] < events[1][:, None]) `````` Chris Jewell committed Jan 16, 2019 46 47 `````` lambdaj = np.sum(np.multiply(B[:, is_infected], waifw[:, is_infected]), axis=0) I0 = np.argmin(t_inf_j[is_infected]) `````` Chris Jewell committed Jan 16, 2019 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 `````` return np.sum(np.log(np.delete(lambdaj, I0))) def interval_intersect(interval_i, interval_j): """Computes intersection between two intervals :param interval_i: a tuple of lower and upper bounds of interval i :param interval_j: a tuple of lower and upper bounds of interval j :return: a matrix of interval intersections """ int_start = np.maximum(interval_i[0][:,None], interval_j[0][None,:]) int_end = np.minimum(interval_i[1][:,None], interval_j[1][None,:]) return np.maximum(0., int_end - int_start) def integral_part(t_inf_j, events, B): """Computes integral part of the infection log likelihood for SIR epidemic :param t_inf_j: vector of infection times for transitting out of S :param events: tuple of arrival and departure times of individuals from I state :param B: the full transmission matrix :return: the integrated infection pressure for the SI epidemic""" i_infected = events[0] < np.inf `````` Chris Jewell committed Jan 16, 2019 72 73 `````` E = interval_intersect((events[0][i_infected], events[1][i_infected]), (np.zeros_like(t_inf_j), t_inf_j)) integral = np.multiply(E, B[i_infected, :]) `````` Chris Jewell committed Jan 16, 2019 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 `````` return np.sum(integral) def log_likelihood(par, t_inf, t_rem, kernel): """Computes the infection log likelihood of a homogeneous SIR epidemic :param par: a (vector of) parameter(s) to be passed to the kernel :param t_inf: a vector of infection times (np.inf if individual is never infected) :param t_rem: a vector of removal times (np.inf if individuals is never removed) :param kernel: an instance of a Kernel class :return: the infection process log likelihood for the SIR model :""" B = kernel(par) prod = prod_part(t_inf, (t_inf, t_rem), B) integral = integral_part(t_inf, (t_inf, t_rem), B) return prod - integral if __name__=="__main__": C = np.matrix([[0, 1, 0, 1, 1], [1, 0, 1, 0, 1], [0, 1, 0, 1, 0], [1, 0, 1, 0, 0], [1, 1, 0, 0, 0]], dtype=np.float64) t_inf = np.array([0, 3, 6, np.inf, np.inf]) t_rem = np.array([5, 8, 10, np.inf, np.inf]) beta = 0.5 kernel = Kernel(C) ll = log_likelihood(beta, t_inf, t_rem, kernel) print("log likelihood =",ll)``````