inhomPPLik.py 3.57 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
"""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
    """
44
    is_infected = t_inf_j < np.inf
Chris Jewell's avatar
Chris Jewell committed
45
    waifw = (events[0][:, None] < t_inf_j[None, :]) & (t_inf_j[None, :] < events[1][:, None])
46 47
    lambdaj = np.sum(np.multiply(B[:, is_infected], waifw[:, is_infected]), axis=0)
    I0 = np.argmin(t_inf_j[is_infected])
Chris Jewell's avatar
Chris Jewell committed
48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71
    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
72 73
    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, :])
Chris Jewell's avatar
Chris Jewell committed
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
    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)