"""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 """ is_infected = t_inf_j < np.inf waifw = (events[0][:, None] < t_inf_j[None, :]) & (t_inf_j[None, :] < events[1][:, None]) lambdaj = np.sum(np.multiply(B[:, is_infected], waifw[:, is_infected]), axis=0) I0 = np.argmin(t_inf_j[is_infected]) 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) 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)