inhomPPLik.R 3.15 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 ``````# Name: inhomPPLik.R # Author: Chris Jewell # Created: 16th January 2019 # Copyright: Chris Jewell 2019 # License: MIT # Warranty: None! Use at your own risk ;-) ################################# # Start with defining some data # ################################# # A population is defined by a contact matrix: C = matrix(c(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), nrow= 5, byrow=T) # Data.frame of (made up) infection and removal times event_times = data.frame(t_inf = c(0, 3, 6, Inf, Inf), t_rem = c(5, 8, 10, Inf, Inf)) #' 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 Kernel = function(c_matrix) { function(beta) { beta * c_matrix } } #' Returns product part of epidemic likelihood #' #' @param t_inf_j -- infection times of everyone beginning in S #' @param events -- the arrival and departure times of individuals in I state #' @param B -- the full transmission matrix #' @return scalar product part of the epidemic likelihoood prod_part = function(t_inf_j, events, B) { waifw = sapply(t_inf_j, function(t) events[,1] < t & events[,2]) lambdaj = colSums(B * waifw) I0 = which.min(t_inf_j) print(waifw) sum(log(lambdaj[-I0])) } #' Computes intersection between two intervals #' #' @param interval_i -- 2D matrix of intervals (rows in result) #' @param interval_j -- 2D matrix of intervals (cols in result) #' @return an n_i X n_j matrix of interval intersections interval_intersect = function(interval_i, interval_j) { int_start = sapply(interval_j[,1], function(x) pmax(x, interval_i[,1])) int_end = sapply(interval_j[,2], function(x) pmin(x, interval_i[,2])) pmax(int_end - int_start, 0) } #' Computes the integral part of the log likelihood #' #' @param t_inf_j -- infection times of everyone beginning in S #' @param events -- the arrival and departure times of individuals in I state #' @param B -- the full transmission matrix #' @return -- scalar integrated infection pressure for the epidemic integral_part = function(t_inf_j, events, B) { i_infected = events[,1] < Inf # Infection time finite E = interval_intersect(events[i_infected,], cbind(0, t_inf_j)) integral = E * B[i_infected,] sum(integral) } #' Computes the log likelihood of an SI epidemic #' #' @param par -- a (vector of) parameter(s) to be passed to the kernel #' @param t_inf -- a vector of infection times (Inf if individual is never infected) #' @param t_rem -- a vector of removal times (Inf if individual is never removed) #' @param kernel -- an instance of a Kernel closure #' @return the log likelihood for the SI epidemic log_likelihood = function(par, t_inf, t_rem, kernel) { B = kernel(par) prod = prod_part(t_inf, cbind(t_inf, t_rem), B) print(prod) integral = integral_part(t_inf, cbind(t_inf, t_rem), B) print(integral) prod - integral } ``````