inhomPPLik.R 3.15 KB
Newer Older
Chris Jewell's avatar
Chris Jewell committed
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 <c.jewell@lancaster.ac.uk>
# 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
}