Initial commit

R/inhomPPLik.R

+# 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
+}
+
+

figPressure.tex

+% TikZ SIR model picture
+%
+% Requries postitioning,arrows,fit,calc,backgrounds TikZ libraries
+
+%\begin{hyphenrules}{nohyphenation}
+
+\begin{tikzpicture}[
+auto,
+node distance=24pt,
+circ/.style={
+circle,
+draw=black,
+align=center,
+minimum height=24pt,
+minimum width=24pt}]
+\node[circ, fill=ForestGreen] (S) {S};
+\node[circ, fill=Red, above right=of S] (I1) {I};
+\node[circ,fill=Red, above left=of S] (I2) {I};
+\node[circ,fill=Red, below left=of S] (I3) {I};
+\node[circ,fill=Red, below right=of S] (I4) {I};
+
+\draw[->,very thick] (I1) to node (I1S) {$\beta$} (S);
+\draw[->,very thick] (I2) to node (I2S) {$\beta$} (S);
+\draw[->,very thick] (I3) to node (I3S) {$\beta$} (S);
+\draw[->,very thick] (I4) to node (I4S) {$\beta$} (S);
+
+\end{tikzpicture}
+%\end{hyphenrules}
\ No newline at end of file

figPressureExt.tex

+\documentclass[xcolor=dvipsnames]{standalone}
+
+\usepackage[usenames, dvipsnames]{xcolor}
+\usepackage{graphicx}
+\usepackage{tikz}
+\usetikzlibrary{positioning,calc,fit,arrows,backgrounds}
+
+\begin{document}
+\input{figPressure}
+\end{document}
\ No newline at end of file

figSIRModel0.tex

+% TikZ SIR model picture
+%
+% Requries postitioning,arrows,fit,calc,backgrounds TikZ libraries
+
+\begin{hyphenrules}{nohyphenation}
+
+\begin{tikzpicture}[
+auto,
+box/.style={
+draw=black,
+align=center,
+minimum height=36pt,
+minimum width=36pt}]
+\node[box, fill=ForestGreen] (S) {S};
+\node[box, fill=Red, right=of S] (I) {I};
+\node[box, fill=Gray, right=of I] (R) {R};
+
+\draw[->,very thick] (S) to node (SI) {} (I);
+\draw[->,very thick] (I) to node (IR) {} (R);
+%\draw[-*,dashed] (I) to [out=270,in=270] (SI);
+
+\end{tikzpicture}
+\end{hyphenrules}
\ No newline at end of file

figSIRModel1.tex

+% TikZ SIR model picture
+%
+% Requries postitioning,arrows,fit,calc,backgrounds TikZ libraries
+
+\begin{hyphenrules}{nohyphenation}
+
+\begin{tikzpicture}[
+auto,
+node distance=48pt,
+box/.style={
+draw=black,
+anchor=center,
+align=center,
+minimum height=36pt,
+minimum width=36pt}]
+\node[box, fill=ForestGreen] (S) {$S(t)$};
+\node[box, fill=Red, right=of S] (I) {$I(t)$};
+\node[box, fill=Gray, right=of I] (R) {$R(t)$};
+
+\draw[->,very thick] (S) to node (SI) {$\beta I(t)S(t)$} (I);
+\draw[->,very thick] (I) to node (IR) {$\gamma I(t)$} (R);
+\draw[-*,dashed] (I) to [out=240,in=270,looseness=2] (SI);
+
+\end{tikzpicture}
+\end{hyphenrules}
\ No newline at end of file

python/inhomPPLik.py

+"""Name: inhomPPLik.py
+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 ;-)
+"""
+
+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)