inhomPPLik.py 3.55 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
106
107
108
109
"""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)