diff --git a/R/inhomPPLik.R b/R/inhomPPLik.R index 35f85fc68160179c2ea2017de24b982d25214be9..bbb4396a7975b4232913f762d5b1925c1f18ac25 100644 --- a/R/inhomPPLik.R +++ b/R/inhomPPLik.R @@ -52,9 +52,10 @@ Kernel = function(c_matrix) { #' @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) + is_infected = t_inf_j < Inf + waifw = sapply(t_inf_j, function(t) events[,1] < t & t < events[,2]) + lambdaj = colSums(B[,is_infected] * waifw[,is_infected]) + I0 = which.min(t_inf_j[is_infected]) print(waifw) sum(log(lambdaj[-I0])) } @@ -79,7 +80,7 @@ interval_intersect = function(interval_i, interval_j) { #' @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 + i_infected = events[,1] < Inf E = interval_intersect(events[i_infected,], cbind(0, t_inf_j)) integral = E * B[i_infected,] sum(integral) diff --git a/epiLikelihoods.pdf b/epiLikelihoods.pdf new file mode 100644 index 0000000000000000000000000000000000000000..fa574cc5d3c6f23e7ade1da8ca8d4a37d5d6216d Binary files /dev/null and b/epiLikelihoods.pdf differ diff --git a/epiLikelihoods.tex b/epiLikelihoods.tex index ce9877ca275b3d8a505af4a3eccd936224db32e4..64eb078a17f554b512582414fe5ac4478ecfbfcd 100644 --- a/epiLikelihoods.tex +++ b/epiLikelihoods.tex @@ -404,7 +404,7 @@ $$\begin{itemize} \item We can now calculate \lambda_i(t) for each member i of \mathcal{S}: -$$\lambda_j(t) = \beta \sum{i \in \mathcal{I}(t)} c_{ij}$$+$$\lambda_j(t) = \beta \sum_{i \in \mathcal{I}(t)} c_{ij}$$\item Likelihood takes a product over \emph{all} individuals in the population \mathcal{P}: \begin{eqnarray*} L(\bm{t}^{SI} | \beta) & = & \prod_{j \in \mathcal{P}} \left[ \lambda_j(s^-) e^{-\int_{I_0}^{s} \lambda_j(t) \mathrm{d}t} \right]^{t^{SI}_j < T_{max}} \left[ e^{-\int_{I_0}^{s} \lambda_j(t) \mathrm{d}t} \right]^{t^{SI}_j \geq T_{max}} \\ @@ -516,25 +516,23 @@$$W:\;w_{ij} = t^{inf}_i < t^{inf}_j < t^{rem}_i \begin{lstlisting}[language=R] # R 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) + is_infected = t_inf_j < Inf + waifw = sapply(t_inf_j, function(t) events[,1] < t & t < events[,2]) + lambdaj = colSums(B[,is_infected] * waifw[,is_infected]) + I0 = which.min(t_inf_j[is_infected]) sum(log(lambdaj[-I0])) -} - \end{lstlisting} +} \end{lstlisting} \begin{lstlisting}[language=Python] # Python import numpy as np - def prod_part(t_inf_j, events, B): - infec, remove = events[:,0], events[:,1] - waifw = (infec[:,None] < t_inf_j[None,:]) & (t_inf_j[None,:] < remove[:,None]) - lambdaj = np.sum(waifw * B, axis=0) - I0 = np.argmin(t_inf_j) - return np.sum(np.log(np.delete(lambdaj,I0))) + is_infected = t_inf_j < np.inf + waifw = (events[0][:, None] < t_inf_j[None, :]) & (t_inf_j[None, :] < events[1][:, None]) + lambdaj = np.sum(np.multiply(B[:, is_infected], waifw[:, is_infected]), axis=0) + I0 = np.argmin(t_inf_j[is_infected]) + return np.sum(np.log(np.delete(lambdaj, I0))) \end{lstlisting} - \end{frame} @@ -594,18 +592,18 @@ All we have to do to calculate the integral is multiply relevant rows of $B$ and \begin{lstlisting}[language=R] # R 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) + i_infected = events[,1] < Inf + E = interval_intersect(events[i_infected,], cbind(0, t_inf_j)) + integral = E * B[i_infected,] + sum(integral) } \end{lstlisting} \begin{lstlisting}[language=Python] # Python def integral_part(t_inf_j, events, B): - i_infected = events[:,1] < np.inf - E = interval_intersect((events[:,0],events[:,1]), (0., t_inf_j)) + 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 = E * B[i_infected,:] return np.sum(integral) \end{lstlisting} @@ -653,8 +651,8 @@ def log_likelihood(par, t_inf, t_rem, kernel): \vspace{12pt} \item Measure your code in terms of \textcolor{blue}{lines not written} \begin{itemize} - \item 28 lines of R - \item x lines of Python + \item 29 lines of R + \item 25 lines of Python \end{itemize} \end{itemize} diff --git a/python/inhomPPLik.py b/python/inhomPPLik.py index 2e7caf5cb81d51fd4d8a87510ed3ba759a62ad7d..11463d228ff2d2555b2e6b763fc6e7a04a62c435 100644 --- a/python/inhomPPLik.py +++ b/python/inhomPPLik.py @@ -41,10 +41,10 @@ def prod_part(t_inf_j, events, B): :param B: the full transmission matrix :return: scalar product part of the epidemic likelihoood """ + is_infected = t_inf_j < np.inf 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) + lambdaj = np.sum(np.multiply(B[:, is_infected], waifw[:, is_infected]), axis=0) + I0 = np.argmin(t_inf_j[is_infected]) return np.sum(np.log(np.delete(lambdaj, I0))) @@ -69,8 +69,8 @@ def integral_part(t_inf_j, events, B): :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,:]) + 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) @@ -85,9 +85,7 @@ def log_likelihood(par, t_inf, t_rem, kernel): :""" 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