inhomPPLik.py 3.55 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 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 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 108 109 ``````"""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 """ waifw = (events[0][:, None] < t_inf_j[None, :]) & (t_inf_j[None, :] < events[1][:, None]) lambdaj = np.sum(waifw * B, axis=0) I0 = np.argmin(t_inf_j) print(waifw) 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 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,:]) 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) print("prod = ", prod) integral = integral_part(t_inf, (t_inf, t_rem), B) print('integral = ', integral) 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)``````